|
mFES - molecular Finite Element Solver
0.4
|
00001 00011 #ifndef PDE_H 00012 #define PDE_H 00013 00014 #include "Residue.h" 00015 #include <vector> 00016 #include "Defs.h" 00017 00032 class PDE { 00033 public: 00034 00035 PDE(){ 00036 ; 00037 } 00038 void writeEnergy(string fileName, INI &ini){ 00039 string solOrder = ini.get<string>("solver.solution_order"); 00040 string maxsteps = ini.get<string>("solver.maxsteps"); 00041 string eps_in = ini.get<string>("experiment.eps_in"); 00042 string eps_out = ini.get<string>("experiment.eps_out"); 00043 00044 float ionc = 0; 00045 boost::optional<string> ionc_opt = ini.get_optional<string>("experiment.ionc"); 00046 00047 if(ionc_opt){ 00048 ionc = ini.get<float>("experiment.ionc"); 00049 } 00050 00051 00052 ofstream energyFile; 00053 energyFile.open("energy.pde"); 00054 00055 energyFile << "mesh = protein.vol" << endl; 00056 energyFile << endl; 00057 energyFile << "shared = pointcharges" << endl; 00058 energyFile << "shared = energydiff" << endl; 00059 energyFile << "shared = writePotatAscii" << endl; 00060 energyFile << endl; 00061 energyFile << "define constant eps0 = 8.8541878e-22" << endl; 00062 energyFile << "define constant q0 = 1.60217646e-19" << endl; 00063 energyFile << "define constant ionc = " << ionc << endl; 00064 energyFile << endl; 00065 energyFile << "define constant heapsize = " << HEAPSIZE << endl; 00066 energyFile << endl; 00067 energyFile << "define coefficient epsilon_solv" << endl; 00068 if (ionc == 0) 00069 energyFile << "(" << eps_out << "*eps0),(" << eps_in << "*eps0),(" << eps_out << "*eps0)" << endl; 00070 else 00071 energyFile << "(" << eps_out << "*eps0),(" << eps_out << "*eps0),(" << eps_in << "*eps0),("<< eps_out <<"*eps0)" << endl; 00072 00073 energyFile << endl; 00074 energyFile << "define coefficient epsilon_ref" << endl; 00075 energyFile << "(" << eps_in << "*eps0),(" << eps_in << "*eps0),("<< eps_in <<"*eps0),("<< eps_in <<"*eps0)" << endl; 00076 energyFile << endl; 00077 energyFile << "define coefficient kappa" << endl; 00078 energyFile << "(("<< DL <<"*eps0)*ionc),0,0,0" << endl; 00079 energyFile << endl; 00080 energyFile << "define fespace v -order="<<solOrder<<" -type=h1ho -dirichlet=[1]" << endl; 00081 energyFile << endl; 00082 energyFile << "define gridfunction u_solv -fespace=v" << endl; 00083 energyFile << "define gridfunction u_ref -fespace=v" << endl; 00084 energyFile << endl; 00085 energyFile << "define bilinearform a_solv -fespace=v -symmetric" << endl; 00086 energyFile << "laplace epsilon_solv" << endl; 00087 energyFile << "mass kappa" << endl; 00088 energyFile << endl; 00089 energyFile << "define bilinearform a_ref -fespace=v -symmetric" << endl; 00090 energyFile << "laplace epsilon_ref" << endl; 00091 energyFile << "mass kappa" << endl; 00092 energyFile << endl; 00093 energyFile << "define linearform f -fespace=v" << endl; 00094 energyFile << endl; 00095 energyFile << "numproc pointcharges ps1 -linearform=f -pqrfile="<<fileName<<" -interpolate" << endl; 00096 energyFile << endl; 00097 energyFile << "define preconditioner c -type=multigrid -bilinearform=a_solv -inverse=mumps" << endl; 00098 energyFile << "define preconditioner c0 -type=multigrid -bilinearform=a_ref -inverse=mumps" << endl; 00099 energyFile << endl; 00100 energyFile << "numproc bvp np1 -gridfunction=u_solv -bilinearform=a_solv -linearform=f -preconditioner=c -maxsteps=" << maxsteps << endl; 00101 energyFile << "numproc bvp np10 -gridfunction=u_ref -bilinearform=a_ref -linearform=f -preconditioner=c0 -maxsteps=" << maxsteps << endl; 00102 energyFile << endl; 00103 energyFile << "numproc energydiff npeval -gridfunction=u_solv -gridfunction0=u_ref -pqrfile=" << fileName << endl; 00104 energyFile.close(); 00105 } 00106 00107 void writeCycle0(vector<Residue> titGroupList, INI &ini){ 00108 string solOrder = ini.get<string>("solver.solution_order"); 00109 string maxsteps = ini.get<string>("solver.maxsteps"); 00110 string eps_in = ini.get<string>("experiment.eps_in"); 00111 string eps_out = ini.get<string>("experiment.eps_out"); 00112 00113 string extend_preconditioner = " -cylce=6 -smoother=block -coarsetype=direct -coarsesmoothingsteps=0 -notest"; 00114 00115 float ionc = 0; 00116 boost::optional<string> ionc_opt = ini.get_optional<string>("experiment.ionc"); 00117 00118 if(ionc_opt){ 00119 ionc = ini.get<float>("experiment.ionc"); 00120 } 00121 00122 00123 for (unsigned int i = 0; i < titGroupList.size(); i++){ 00124 Residue currentTitGroup = titGroupList.at(i); 00125 string prefix = currentTitGroup.getIdentifier(); 00126 00127 // if ( !boost::filesystem::exists( "cycle0."+prefix+".potat" ) ){ 00128 ofstream potfile; 00129 00130 // if (i > 0){ 00131 // potfile.open("pka_cycle0_"+prefix+".pde", ios::in | ios::out | ios::app); 00132 //} 00133 //else { 00134 if(boost::filesystem::exists("pka_cycle0_"+prefix+".pde")){ 00135 boost::filesystem::remove("pka_cycle0_"+prefix+".pde"); 00136 cout << "removing pde file: " << "pka_cycle0_"+prefix+".pde" << endl; 00137 } 00138 potfile.open("pka_cycle0_"+prefix+".pde", ios::out); 00139 // } 00140 potfile << "###################################################" << endl; 00141 potfile << "#" << endl; 00142 potfile << "# Initialization: " << prefix << endl; 00143 potfile << "#" << endl; 00144 potfile << "###################################################" << endl; 00145 potfile << endl; 00146 potfile << "mesh = " << prefix << ".vol" << endl; 00147 potfile << endl; 00148 potfile << "shared = pointcharges" << endl; 00149 potfile << "shared = energydiff" << endl; 00150 potfile << "shared = writePotatAscii" << endl; 00151 potfile << "shared = precalculation" << endl; 00152 potfile << endl; 00153 potfile << "define constant eps0 = 8.8541878e-22" << endl; 00154 potfile << "define constant q0 = 1.60217646e-19" << endl; 00155 potfile << "define constant ionc = " << ionc << endl; 00156 potfile << endl; 00157 potfile << "define constant heapsize = " << HEAPSIZE << endl; 00158 potfile << endl; 00159 potfile << "define coefficient epsilon_solv" << endl; 00160 if (ionc == 0) 00161 potfile << "(" << eps_out << "*eps0),(" << eps_in << "*eps0),(" << eps_out << "*eps0)" << endl; 00162 else 00163 potfile << "(" << eps_out << "*eps0),(" << eps_out << "*eps0),(" << eps_in << "*eps0),("<< eps_out <<"*eps0)" << endl; 00164 potfile << endl; 00165 potfile << "define coefficient epsilon_ref" << endl; 00166 potfile << "(" << eps_in << "*eps0),(" << eps_in << "*eps0),("<< eps_in <<"*eps0),("<< eps_in <<"*eps0)" << endl; 00167 potfile << endl; 00168 potfile << "define coefficient kappa" << endl; 00169 potfile << "((" << DL << "*eps0)*ionc),0,0,0" << endl; 00170 potfile << endl; 00171 00172 00173 unsigned int nrStates = currentTitGroup.getNrStates(); 00174 string create = ""; 00175 string stateIndex = ""; 00176 potfile << endl; 00177 potfile << "define fespace v -order="<< solOrder <<" -type=h1ho -dirichlet=[1]" << endl; 00178 potfile << endl; 00179 potfile << "define gridfunction u_solv_"<<prefix<<" -fespace=v" << endl; 00180 potfile << "define gridfunction u_ref_"<<prefix<<" -fespace=v" << endl; 00181 potfile << endl; 00182 potfile << "define bilinearform a_solv_"<<prefix<<" -fespace=v -symmetric" << endl; 00183 potfile << "laplace epsilon_solv" << endl; 00184 potfile << "mass kappa" << endl; 00185 potfile << endl; 00186 potfile << "define bilinearform a_ref_"<<prefix<<" -fespace=v -symmetric" << endl; 00187 potfile << "laplace epsilon_ref" << endl; 00188 // potfile << "mass kappa" << endl; 00189 00190 potfile << endl; 00191 potfile << "define preconditioner c_solv_"<<prefix<<" -type=multigrid -bilinearform=a_solv_"<<prefix<<" -inverse=mumps" << extend_preconditioner << endl; 00192 potfile << "define preconditioner c_ref_"<<prefix<<" -type=multigrid -bilinearform=a_ref_"<<prefix<<" -inverse=mumps" << extend_preconditioner << endl << endl; 00193 00194 potfile << endl; 00195 potfile << "###################################################" << endl; 00196 potfile << "#" << endl; 00197 potfile << "# Solvation calculations: " << prefix << endl; 00198 potfile << "#" << endl; 00199 potfile << "###################################################" << endl; 00200 potfile << "" << endl; 00201 potfile << "" << endl; 00202 for (unsigned int j = 1; j <= nrStates; j++){ 00203 if (j == 1){ 00204 create = "-file=create"; 00205 stringstream ss; 00206 ss << nrStates; 00207 stateIndex = ss.str(); 00208 } else { 00209 create = ""; 00210 stateIndex = "-1"; 00211 } 00212 potfile << "# state " << j << endl; 00213 potfile << "define linearform f_state_solv_" << j << "_"<<prefix<<" -fespace=v" << endl; 00214 potfile << "numproc pointcharges ps1_solv_"<<prefix<<" -linearform=f_state_solv_" << j << "_"<< prefix << " -pqrfile=" << prefix <<".state" << j <<".pqr -interpolate " << endl; 00215 potfile << "numproc bvp np_solv_"<< j <<"_"<<prefix<<" -gridfunction=u_solv_" << prefix <<" -bilinearform=a_solv_" << prefix <<" -linearform=f_state_solv_"<< j <<"_"<<prefix<<" -preconditioner=c_solv_"<< prefix <<" -maxsteps=" << maxsteps << endl; 00216 potfile << "numproc writepotat npeval_solv_"<<j<<"_"<<prefix<<" -gridfunction=u_solv_"<<prefix<<" -pqrfile=" << prefix << ".state" << j <<".pqr -potatfile=cycle0."<<prefix<<".potat -statenr="<< stateIndex << " " << create << endl; 00217 potfile << endl; 00218 // Debug 00219 // potfile << "numproc energydiff npeval_solv_diff_"<<j<<"_"<<prefix<<" -gridfunction=u_solv_"<<prefix<<" -gridfunction0=u_ref_"<<prefix<<" -pqrfile=" << prefix << ".state" << j <<".pqr" <<endl; 00220 00221 00222 } 00223 potfile << endl; 00224 potfile << "###################################################" << endl; 00225 potfile << "#" << endl; 00226 potfile << "# Reference calculation: " << prefix << endl; 00227 potfile << "#" << endl; 00228 potfile << "###################################################" << endl; 00229 potfile << "" << endl; 00230 for (unsigned int j = 1; j <= nrStates; j++){ 00231 if (j == 1){ 00232 stringstream ss; 00233 ss << nrStates; 00234 stateIndex = ss.str(); 00235 } else { 00236 create = ""; 00237 stateIndex = "-1"; 00238 } 00239 potfile << "# state " << j << endl; 00240 potfile << "define linearform f_state_ref_" << j <<"_"<<prefix<<" -fespace=v" << endl; 00241 potfile << "numproc pointcharges ps_ref_" << j <<"_"<<prefix<<" -linearform=f_state_ref_" << j << "_"<<prefix<<" -pqrfile=" << prefix <<".state" << j << ".pqr -interpolate" << endl; 00242 potfile << "numproc bvp np_ref_" << j << "_"<<prefix<<" -gridfunction=u_ref_"<<prefix<<" -bilinearform=a_ref_"<<prefix<<" -linearform=f_state_ref_" << j << "_"<<prefix<<" -preconditioner=c_ref_"<<prefix<<" -maxsteps=" << maxsteps << endl; 00243 potfile << "numproc writepotat npeval_ref_"<<j<<"_"<<prefix<<" -gridfunction=u_ref_"<<prefix<<" -pqrfile=" << prefix << ".state" << j <<".pqr -potatfile=cycle0." << prefix << ".potat -statenr=" << stateIndex << endl; 00244 00245 // Debug 00246 // potfile << "numproc energydiff npeval_ref_diff_"<<j<<"_"<<prefix<<" -gridfunction=u_solv_"<<prefix<<" -gridfunction0=u_ref_"<<prefix<<" -pqrfile=" << prefix << ".state" << j <<".pqr" <<endl; 00247 } 00248 potfile << endl; 00249 potfile << endl; 00250 potfile.close(); 00251 00252 // } 00253 00254 } 00255 00256 } 00257 00258 void writeCycle1(string molecule, vector<Residue> titGroupList, INI &ini){ 00259 string mol = "protein"; 00260 string solOrder = ini.get<string>("solver.solution_order"); 00261 string maxsteps = ini.get<string>("solver.maxsteps"); 00262 string eps_in = ini.get<string>("experiment.eps_in"); 00263 string eps_out = ini.get<string>("experiment.eps_out"); 00264 00265 float ionc = 0; 00266 boost::optional<string> ionc_opt = ini.get_optional<string>("experiment.ionc"); 00267 00268 if(ionc_opt){ 00269 ionc = ini.get<float>("experiment.ionc"); 00270 } 00271 00272 00273 string extend_preconditioner = " -cylce=6 -smoother=block -coarsetype=direct -coarsesmoothingsteps=0 -notest"; 00274 00275 bool init = true; 00276 00277 if(boost::filesystem::exists("pka_cycle1.pde")){ 00278 boost::filesystem::remove("pka_cycle1.pde"); 00279 cout << "removing pde file: " << "pka_cycle1.pde" << endl; 00280 } 00281 00282 for (unsigned int i = 0; i < titGroupList.size(); i++){ 00283 Residue currentTitGroup = titGroupList.at(i); 00284 string prefix = currentTitGroup.getIdentifier(); 00285 00286 // if ( !boost::filesystem::exists( "cycle1."+prefix+".potat" ) ){ 00287 ofstream potfile; 00288 potfile.open("pka_cycle1.pde", ios::in | ios::out | ios::app); 00289 00290 if (init){ 00291 00292 potfile << "###################################################" << endl; 00293 potfile << "#" << endl; 00294 potfile << "# Initialization: " << mol << endl; 00295 potfile << "#" << endl; 00296 potfile << "###################################################" << endl; 00297 potfile << endl; 00298 potfile << endl; 00299 potfile << "mesh = "<<mol<<".vol" << endl; 00300 potfile << endl; 00301 00302 potfile << "shared = pointcharges" << endl; 00303 potfile << "shared = energydiff" << endl; 00304 potfile << "shared = writePotatAscii" << endl; 00305 potfile << "shared = precalculation" << endl; 00306 potfile << endl; 00307 potfile << "define constant eps0 = 8.8541878e-22" << endl; 00308 potfile << "define constant q0 = 1.60217646e-19" << endl; 00309 potfile << "define constant ionc = " << ionc << endl; 00310 00311 potfile << endl; 00312 potfile << "define constant heapsize = " << HEAPSIZE << endl; 00313 potfile << endl; 00314 potfile << "define coefficient epsilon_solv" << endl; 00315 if (ionc == 0) 00316 potfile << "(" << eps_out << "*eps0),(" << eps_in << "*eps0),(" << eps_out << "*eps0)" << endl; 00317 else 00318 potfile << "(" << eps_out << "*eps0),(" << eps_out << "*eps0),(" << eps_in << "*eps0),("<< eps_out <<"*eps0)" << endl; 00319 potfile << endl; 00320 potfile << "define coefficient epsilon_ref" << endl; 00321 potfile << "(" << eps_in << "*eps0),(" << eps_in << "*eps0),("<< eps_in <<"*eps0),("<< eps_in <<"*eps0)" << endl; 00322 potfile << endl; 00323 potfile << "define coefficient kappa" << endl; 00324 potfile << "((" << DL << "*eps0)*ionc),0,0,0" << endl; 00325 potfile << endl; 00326 00327 potfile << endl; 00328 potfile << "define fespace v -order="<< solOrder<<" -type=h1ho -dirichlet=[1]" << endl; 00329 potfile << endl; 00330 potfile << "define gridfunction u_solv_"<<mol<<" -fespace=v" << endl; 00331 potfile << "define gridfunction u_ref_"<<mol<<" -fespace=v" << endl; 00332 potfile << endl; 00333 potfile << "define bilinearform a_solv_"<<mol<<" -fespace=v -symmetric" << endl; 00334 potfile << "laplace epsilon_solv" << endl; 00335 potfile << "mass kappa" << endl; 00336 00337 potfile << endl; 00338 potfile << "define bilinearform a_ref_"<<mol<<" -fespace=v -symmetric" << endl; 00339 potfile << "laplace epsilon_ref" << endl; 00340 // potfile << "mass kappa" << endl; 00341 00342 potfile << endl; 00343 potfile << "define preconditioner c_solv_"<<mol<<" -type=multigrid -bilinearform=a_solv_"<<mol<<" -inverse=mumps" << extend_preconditioner << endl; 00344 potfile << "define preconditioner c_ref_"<<mol<<" -type=multigrid -bilinearform=a_ref_"<<mol<<" -inverse=mumps" << extend_preconditioner << endl << endl; 00345 potfile << "numproc precalculation pre -pqrfile=" << molecule << endl; 00346 potfile << endl; 00347 init = false; 00348 } 00349 00350 unsigned int nrStates = currentTitGroup.getNrStates(); 00351 string create = ""; 00352 string stateIndex = ""; 00353 00354 potfile << "###################################################" << endl; 00355 potfile << "#" << endl; 00356 potfile << "# Solvation calculations: " << prefix << endl; 00357 potfile << "#" << endl; 00358 potfile << "###################################################" << endl; 00359 potfile << "" << endl; 00360 potfile << "" << endl; 00361 for (unsigned int j = 1; j <= nrStates; j++){ 00362 if (j == 1){ 00363 create = "-file=create"; 00364 stringstream ss; 00365 ss << nrStates; 00366 stateIndex = ss.str(); 00367 } else { 00368 create = ""; 00369 stateIndex = "-1"; 00370 } 00371 potfile << "# state " << j << endl; 00372 potfile << "define linearform f_state_solv_" << j << "_"<<prefix<<" -fespace=v" << endl; 00373 potfile << "numproc pointcharges ps1_solv_"<<prefix<<" -linearform=f_state_solv_" << j << "_"<< prefix << " -pqrfile=" << prefix <<".state" << j <<".pqr -interpolate " << endl; 00374 potfile << "numproc bvp np_solv_"<< j <<"_"<<prefix<<" -gridfunction=u_solv_" << mol <<" -bilinearform=a_solv_" << mol <<" -linearform=f_state_solv_"<< j <<"_"<<prefix<<" -preconditioner=c_solv_"<< mol <<" -maxsteps=" << maxsteps << endl; 00375 potfile << "numproc writepotat npeval_solv_"<<j<<"_"<<prefix<<" -gridfunction=u_solv_"<<mol<<" -pqrfile=" << molecule <<" -potatfile=cycle1."<<prefix<<".potat -statenr="<< stateIndex << " " << create << " -predef" << endl; 00376 00377 // Debug 00378 // potfile << "numproc energydiff npeval_diff_"<<j<<"_"<<prefix<<" -gridfunction=u_solv_"<<mol<<" -gridfunction0=u_ref_"<<mol<<" -pqrfile="<<molecule<<endl; 00379 potfile << endl; 00380 00381 } 00382 potfile << endl; 00383 potfile << "###################################################" << endl; 00384 potfile << "#" << endl; 00385 potfile << "# Reference calculation: " << prefix << endl; 00386 potfile << "#" << endl; 00387 potfile << "###################################################" << endl; 00388 potfile << "" << endl; 00389 for (unsigned int j = 1; j <= nrStates; j++){ 00390 if (j == 1){ 00391 stringstream ss; 00392 ss << nrStates; 00393 stateIndex = ss.str(); 00394 } else { 00395 create = ""; 00396 stateIndex = "-1"; 00397 } 00398 potfile << "# state " << j << endl; 00399 potfile << "define linearform f_state_ref_" << j <<"_"<<prefix<<" -fespace=v" << endl; 00400 potfile << "numproc pointcharges ps_ref_" << j <<"_"<<prefix<<" -linearform=f_state_ref_" << j << "_"<<prefix<<" -pqrfile=" << prefix <<".state" << j << ".pqr -interpolate" << endl; 00401 potfile << "numproc bvp np_ref_" << j << "_"<<prefix<<" -gridfunction=u_ref_"<<mol<<" -bilinearform=a_ref_"<<mol<<" -linearform=f_state_ref_" << j << "_"<<prefix<<" -preconditioner=c_ref_"<<mol<<" -maxsteps=" << maxsteps << endl; 00402 potfile << "numproc writepotat npeval_ref_"<<j<<"_"<<prefix<<" -gridfunction=u_ref_"<<mol<<" -pqrfile=" << molecule << " -potatfile=cycle1." << prefix << ".potat -statenr=" << stateIndex << " -predef" << endl; 00403 00404 // Debug 00405 // potfile << "numproc energydiff npeval_diff_"<<j<<"_"<<prefix<<" -gridfunction=u_solv_"<<mol<<" -gridfunction0=u_ref_"<<mol<<" -pqrfile="<<molecule<<endl; 00406 00407 } 00408 potfile << endl; 00409 potfile << endl; 00410 potfile.close(); 00411 // } 00412 00413 } 00414 00415 } 00416 00417 }; 00418 00419 00420 #endif
1.7.6.1