|
mFES - molecular Finite Element Solver
0.4
|
00001 00015 #include <solve.hpp> 00016 #include <iostream> 00017 #include <string> 00018 #include <vector> 00019 #include <sstream> 00020 #include <set> 00021 00022 #include "Atom.h" 00023 #include "Charge.h" 00024 #include "Residue.h" 00025 #include "Model.h" 00026 #include "ST.h" 00027 #include "tCycle.h" 00028 #include "PDE.h" 00029 00030 using namespace std; 00031 00033 namespace netgen 00034 { 00035 int h_argc; 00036 char ** h_argv; 00037 } 00038 00039 char **argv; 00040 ngsolve::MyMPI mympi(0, argv); 00041 00043 typedef boost::property_tree::ptree INI; 00044 00058 class PQR { 00059 public: 00067 PQR(string _fileName): 00068 fileName(_fileName) 00069 { 00070 if(parse()) 00071 cout << atomList.size() << " atom(s) found." << endl; 00072 else { 00073 cout << "An error occured while parsing file..." << endl; 00074 exit(0); 00075 } 00076 } 00077 00086 void addInfo(INI &ini){ 00087 string _stFilesFolder = ini.get<string>("pka.st_folder"); 00088 stFilesFolder = _stFilesFolder; 00089 00090 string _calcCTE = ini.get<string>("pka.calc_cte"); 00091 if (_calcCTE == "yes") 00092 calcCTE = true; 00093 00094 string _calcNTE = ini.get<string>("pka.calc_nte"); 00095 if (_calcNTE == "yes") 00096 calcNTE = true; 00097 00098 string _explicitModels = ini.get<string>("pka.explicit_models"); 00099 if (_explicitModels == "yes") 00100 explicitModels = true; 00101 00102 } 00103 00104 00112 void parseSTFolder(){ 00113 boost::filesystem::path full_path( stFilesFolder ); 00114 00115 for ( boost::filesystem::directory_iterator it( full_path ); 00116 it != boost::filesystem::directory_iterator(); ++it ) 00117 { 00118 string titGroupName = it->path().stem().string(); 00119 string ext = it->path().extension().string(); 00120 boost::to_lower(ext); 00121 00124 if ( boost::filesystem::is_regular_file( it->status() ) 00125 && ( ext == ".st" )) 00126 { 00128 cout << "Parsing... " << it->path().string() << endl; 00129 ifstream in(it->path().string().c_str()); 00130 if (!in) { 00131 in.close(); 00132 cout << "Cannot open st file:" << it->path().string() << endl; 00133 exit(0); 00134 } 00135 string currentLine, temp; 00136 00137 char state; float deltaG; float shiftCycle0 = NOT_IN_ST; 00138 string atomName; double charge; 00139 vector<Charge> rules; 00140 unsigned int stateNr = 0; 00141 while( !in.eof() ) { 00142 getline(in, currentLine); 00143 if ( currentLine.find("ATOM", 0)==string::npos ) { 00144 if (stateNr > 0){ 00145 ST newST(titGroupName, state, stateNr, deltaG, rules, shiftCycle0); 00146 stList.push_back(newST); 00147 rules.clear(); 00148 } 00149 istringstream ss(currentLine); 00150 ss >> deltaG >> temp>> state >> temp >> shiftCycle0; 00151 stateNr++; 00152 } else { 00153 00154 temp = currentLine.substr(12,4); 00155 istringstream ss(temp); 00156 ss >> atomName; 00157 00158 temp = currentLine.substr(54,6); 00159 if( sscanf(temp.c_str(), "%lg", &(charge))!=1 ){ 00160 cout << "could not load charge in: " << endl; 00161 cout << currentLine << endl; 00162 exit(0); 00163 } 00164 Charge newCharge(atomName, charge); 00165 rules.push_back(newCharge); 00166 } 00167 00168 } 00169 00170 } 00171 } 00172 00173 cout << stList.size() << " patch rules added." << endl; 00174 } 00175 00176 00185 void parseSTFiles(){ 00186 parseSTFolder(); // we have all ST file information 00187 STparsed = true; 00188 } 00189 00190 00197 int parse(){ 00198 vector<Atom> atomList; 00199 string temp; 00200 00201 ifstream in(fileName.c_str()); 00202 00203 if (!in) { 00204 in.close(); 00205 cout << "Cannot open molecule file:" << fileName << endl; 00206 exit(0); 00207 } 00208 string pqrline; 00209 00210 while( !in.eof() ) { 00211 getline(in, pqrline); 00212 if ( pqrline.find("ATOM", 0)==string::npos ) continue; 00213 00214 if ( pqrline.size() < 72 ){ 00215 cout << "Molecule file is not well structured. Example: " << endl; 00216 cout << "ATOM 1 N VAL A 3 16.875 37.901 43.478-0.470 1.850 A" << endl; 00217 exit(0); 00218 } 00219 00220 int atomNumber; string atomName; string residueName; char confID; char chainID; 00221 int residueNumber; double charge; double radius; string segName; 00222 00223 // atom number 00224 temp = pqrline.substr(6,5); 00225 if( sscanf(temp.c_str(), "%i", &(atomNumber))!=1 ) 00226 return 0; 00227 00228 // atom name 00229 temp = pqrline.substr(12,4); 00230 istringstream ss(temp); 00231 ss >> atomName; 00232 00233 // conformer id 00234 temp = pqrline.substr(16,1); 00235 if( sscanf(temp.c_str(), "%c", &(confID))!=1 ) 00236 confID = ' '; 00237 00238 // residue name 00239 residueName = pqrline.substr(17,3); 00240 00241 // chain id 00242 temp = pqrline.substr(21,1); 00243 if( sscanf(temp.c_str(), "%c", &(chainID))!=1 ) 00244 chainID = ' '; 00245 00246 // residue number 00247 temp = pqrline.substr(22,4); 00248 if( sscanf(temp.c_str(), "%i", &(residueNumber))!=1 ) 00249 return 0; 00250 00251 // coordinates 00252 double x, y, z; 00253 temp = pqrline.substr(30,8); 00254 if( sscanf(temp.c_str(), "%lg", &x)!=1 ) 00255 return 0; 00256 00257 temp = pqrline.substr(38,8); 00258 if( sscanf(temp.c_str(), "%lg", &y)!=1 ) 00259 return 0; 00260 temp = pqrline.substr(46,8); 00261 if( sscanf(temp.c_str(), "%lg", &z)!=1 ) 00262 return 0; 00263 00264 Point3<float> coord(x, y, z); 00265 00266 // charge 00267 temp = pqrline.substr(54,6); 00268 if( sscanf(temp.c_str(), "%lg", &(charge))!=1 ) 00269 return 0; 00270 00271 // radius 00272 temp = pqrline.substr(60,6); 00273 if( sscanf(temp.c_str(), "%lg", &(radius))!=1 ) 00274 return 0; 00275 00276 // segment name 00277 try 00278 { 00279 temp = pqrline.substr(72,4); 00280 istringstream ss(temp); 00281 ss >> segName; 00282 } 00283 catch(out_of_range& x) 00284 { 00285 segName = " "; 00286 } 00287 Atom newAtom(atomNumber, atomName, confID, residueName, chainID, residueNumber, coord, charge, radius, segName); 00288 this->atomList.push_back(newAtom); 00289 } 00290 00291 return 1; 00292 00293 } 00294 00303 void calcModel(INI& ini){ 00305 if (ini.get<string>("experiment.cavity") == "yes") 00306 molecule.calcModel(atomList, ini, "cavity.vol", true); 00307 00309 boost::optional<string> ionc = ini.get_optional<string>("experiment.ionc"); 00310 if(ionc){ 00311 if(ini.get<float>("experiment.ionc") != 0){ 00312 molecule.calcIonLayer(atomList, ini, "exclusion.vol"); 00313 } 00314 } 00315 00317 molecule.calcModel(atomList, ini); 00318 } 00319 00327 void calcDeltaG(){ 00328 cout << "Calculating deltaG ..." << endl; 00329 ngsolve::PDE pde; 00330 00331 string pdeFile = "energy.pde"; 00332 00333 try { 00334 pde.LoadPDE (pdeFile.c_str()); 00335 pde.Solve(); 00336 } catch(ngstd::Exception & e) { 00337 std::cout << "Caught exception: " << std::endl 00338 << e.What() << std::endl; 00339 } 00340 00341 } 00342 00352 void calcResidues(INI &ini){ 00353 Residue newResidue; 00354 vector<Atom> currentAtomList; 00355 Atom lastAtom; 00356 int oldResID = -1; 00357 int currentResID = oldResID; 00358 00359 for (unsigned int i = 0; i < atomList.size(); i++){ 00360 Atom currentAtom = atomList.at(i); 00361 currentResID = currentAtom.getResidueNumber(); 00362 if (oldResID == -1){ // First element 00363 oldResID = currentResID; 00364 lastAtom = currentAtom; 00365 } 00366 //if (currentAtom is last element){ // Last element 00367 if (i == (atomList.size()-1)){ 00368 oldResID = -2; 00369 currentAtomList.push_back(currentAtom); 00370 } 00371 00372 if (oldResID != currentResID){ 00373 newResidue.setResidueName(lastAtom.getResidueName()); 00374 newResidue.setResidueNumber(lastAtom.getResidueNumber()); 00375 newResidue.setChainID(lastAtom.getChainID()); 00376 newResidue.setAtomList(currentAtomList); 00377 residueList.push_back(newResidue); 00378 // if (newResidue.isTitGroup()){ 00379 // titGroupList.push_back(newResidue); 00380 // } 00381 currentAtomList.clear(); 00382 } 00383 currentAtomList.push_back(currentAtom); 00384 lastAtom = currentAtom; 00385 oldResID = currentResID; 00386 } 00387 00388 // Read in sites file 00389 string sitesFile = ini.get<string>("pka.sites_file"); 00390 set<string> sitesData; 00391 ifstream in(sitesFile.c_str()); 00392 if (!in) { 00393 in.close(); 00394 cout << "Cannot open sites file:" << sitesFile << endl; 00395 exit(0); 00396 } 00397 string currentLine, _chainID, _resNumber, _resName; 00398 00399 while( !in.eof() ) { 00400 getline(in, currentLine); 00401 stringstream ss(currentLine); 00402 ss >> _chainID >> _resNumber >> _resName ; 00403 string name = _resName+"-"+_resNumber+"_"+_chainID; 00404 sitesData.insert(name); 00405 00406 } 00407 // Insert residues into titGroupList 00408 for (unsigned int i = 0; i < residueList.size(); i++){ 00409 Residue currentResidue = residueList.at(i); 00410 if (currentResidue.isTitGroup() && sitesData.find(currentResidue.getIdentifier()) != sitesData.end()){ 00411 vector<Atom> currentAtomList = currentResidue.getAtomList(); 00412 if (i >= 1){ 00413 Residue before = residueList.at(i-1); 00414 vector<Atom> newAtomList = before.getAtomList(); 00415 unsigned int atomListSize = newAtomList.size(); 00416 Atom C = newAtomList.at(atomListSize-2); 00417 Atom O = newAtomList.at(atomListSize-1);; 00418 currentAtomList.insert(currentAtomList.begin(), O); 00419 currentAtomList.insert(currentAtomList.begin(), C); 00420 } 00421 00422 if (i+1 < residueList.size()){ 00423 Residue after = residueList.at(i+1); 00424 vector<Atom> newAtomList = after.getAtomList(); 00425 Atom N = newAtomList.at(0); 00426 Atom HN = newAtomList.at(1); 00427 Atom CA = newAtomList.at(2); 00428 currentAtomList.insert(currentAtomList.end(), N); 00429 currentAtomList.insert(currentAtomList.end(), HN); 00430 currentAtomList.insert(currentAtomList.end(), CA); 00431 } 00432 currentResidue.setAtomList(currentAtomList); 00433 titGroupList.push_back(currentResidue); 00434 } 00435 } 00436 00437 // Patching also N-terminus? 00438 if (calcNTE){ 00439 Residue firstResidue = residueList.at(0); 00440 firstResidue.setResidueName("NTE"); 00441 int residueNumber = firstResidue.getResidueNumber(); 00442 firstResidue.setResidueNumber(residueNumber--); 00443 vector<Atom> currentAtomList = firstResidue.getAtomList(); 00444 00445 Residue after = residueList.at(1); 00446 vector<Atom> newAtomList = after.getAtomList(); 00447 Atom N = newAtomList.at(0); 00448 Atom HN = newAtomList.at(1); 00449 Atom CA = newAtomList.at(2); 00450 for (unsigned int i = 0; i < currentAtomList.size(); i++) 00451 currentAtomList.at(i).setResidueName("NTE"); 00452 00453 currentAtomList.insert(currentAtomList.end(), N); 00454 currentAtomList.insert(currentAtomList.end(), HN); 00455 currentAtomList.insert(currentAtomList.end(), CA); 00456 firstResidue.setAtomList(currentAtomList); 00457 titGroupList.insert(titGroupList.begin(), firstResidue); 00458 } 00459 00460 // Patching also ending C-terminus? 00461 if (calcCTE){ 00462 Residue lastResidue = residueList.at(residueList.size()-1); 00463 lastResidue.setResidueName("CTE"); 00464 int residueNumber = lastResidue.getResidueNumber(); 00465 lastResidue.setResidueNumber(residueNumber++); 00466 vector<Atom> currentAtomList = lastResidue.getAtomList(); 00467 00468 Residue before = residueList.at(residueList.size()-2); 00469 vector<Atom> newAtomList = before.getAtomList(); 00470 unsigned int atomListSize = newAtomList.size(); 00471 Atom C = newAtomList.at(atomListSize-2); 00472 Atom O = newAtomList.at(atomListSize-1); 00473 for (unsigned int i = 0; i < currentAtomList.size(); i++) 00474 currentAtomList.at(i).setResidueName("CTE"); 00475 00476 currentAtomList.insert(currentAtomList.begin(), O); 00477 currentAtomList.insert(currentAtomList.begin(), C); 00478 lastResidue.setAtomList(currentAtomList); 00479 titGroupList.insert(titGroupList.end(), lastResidue); 00480 } 00481 00482 cout << residueList.size() << " residue(s) with " << titGroupList.size() << " titratable group(s) found." << endl; 00483 /* cout << "titration groups ordering" << endl; 00484 for (unsigned int i = 0; i < titGroupList.size(); i++){ 00485 Residue currentTitGroup = titGroupList.at(i); 00486 cout << currentTitGroup.getResidueName() << endl; 00487 }*/ 00488 00489 00490 } 00491 00502 void calcExplicitModels(INI &ini){ 00503 for (unsigned int i = 0; i < titGroupList.size(); i++){ 00504 Residue currentTitGroup = titGroupList.at(i); 00505 // Write out PQR files 00506 for (unsigned int j = 1; j <= currentTitGroup.getNrStates(); j++){ 00507 writePQR(currentTitGroup, j); 00508 } 00509 00510 if (getCycle0Shift(currentTitGroup.getResidueName(), 1) == NOT_IN_ST || explicitModels){ 00511 string fname = currentTitGroup.getIdentifier()+".vol"; 00512 vector<Atom> currentAtomList = currentTitGroup.getAtomList(); 00513 boost::optional<string> ionc = ini.get_optional<string>("experiment.ionc"); 00514 if(ionc){ 00515 if(ini.get<float>("experiment.ionc") != 0){ 00516 molecule.calcIonLayer(atomList, ini, currentTitGroup.getIdentifier()+".vol_exclusion"); 00517 } 00518 } 00519 00520 molecule.calcModel(currentAtomList, ini, fname); 00521 } 00522 00523 } 00524 00525 } 00526 00527 00537 void writePQR(vector<Atom> &aList, string fileName){ 00538 ofstream pqr; 00539 pqr.open (fileName); 00540 for (unsigned int i = 0; i < aList.size(); i++){ 00541 Atom currentAtom = aList.at(i); 00542 pqr << currentAtom.pqrLine() << endl; 00543 } 00544 pqr.close(); 00545 00546 } 00547 00557 void writePQR(Residue &titGroup, int stateNr){ 00558 for (unsigned int i = 0; i < stList.size(); i++){ 00559 ST currentST = stList.at(i); 00560 if (currentST.getTitGroupName() == titGroup.getResidueName() && 00561 currentST.getStateNr() == stateNr){ 00562 vector<Atom> currentAtomList = titGroup.getAtomList(); 00563 patch(currentAtomList, currentST); 00564 stringstream ss; 00565 ss << titGroup.getIdentifier() << ".state" << stateNr << ".pqr"; 00566 ofstream pqr; 00567 pqr.open (ss.str()); 00568 for (unsigned int j = 0; j < currentAtomList.size(); j++){ 00569 Atom currentAtom = currentAtomList.at(j); 00570 pqr << currentAtom.pqrLine() << endl; 00571 } 00572 pqr.close(); 00573 break; 00574 } 00575 } 00576 } 00577 00588 void switchState(Residue &titGroup, int stateNr){ 00589 for (unsigned int i = 0; i < stList.size(); i++){ 00590 ST currentST = stList.at(i); 00591 if (currentST.getTitGroupName() == titGroup.getResidueName() && 00592 currentST.getStateNr() == stateNr){ 00593 vector<Atom> currentAtomList = titGroup.getAtomList(); 00594 patch(currentAtomList, currentST); 00595 titGroup.setAtomList(currentAtomList); 00596 break; 00597 } 00598 } 00599 00600 } 00601 00612 void neutralState(Residue &titGroup, int stateNr){ 00613 for (unsigned int i = 0; i < stList.size(); i++){ 00614 ST currentST = stList.at(i); 00615 if (currentST.getTitGroupName() == titGroup.getResidueName() && 00616 currentST.getStateNr() == stateNr){ 00617 vector<Atom> currentAtomList = titGroup.getAtomList(); 00618 patchNeutral(currentAtomList, currentST); 00619 titGroup.setAtomList(currentAtomList); 00620 break; 00621 } 00622 } 00623 00624 } 00625 00626 00638 void patch(vector<Atom> &atomList, ST &st){ 00639 vector<Charge> rules = st.getRules(); 00640 for (unsigned int i = 0; i < atomList.size(); i++){ 00641 Atom currentAtom = atomList.at(i); 00642 bool found = false; 00643 for (unsigned int j = 0; j < rules.size(); j++){ 00644 Charge currentRule = rules.at(j); 00645 if (currentRule.getAtomName() == currentAtom.getAtomName() && 00646 st.getTitGroupName() == currentAtom.getResidueName()){ 00647 float currentCharge = currentRule.getCharge(); 00648 if ( currentCharge == 0){ 00649 currentCharge = 0.01; 00650 } 00651 00652 currentAtom.setCharge(currentCharge); 00653 found = true; 00654 break; 00655 } 00656 } 00657 if (!found){ 00658 currentAtom.setCharge(0); 00659 } 00660 atomList.at(i) = currentAtom; 00661 } 00662 } 00663 00674 void patchNeutral(vector<Atom> &atomList, ST &st){ 00675 vector<Charge> rules = st.getRules(); 00676 for (unsigned int i = 0; i < atomList.size(); i++){ 00677 Atom currentAtom = atomList.at(i); 00678 bool found = false; 00679 for (unsigned int j = 0; j < rules.size(); j++){ 00680 Charge currentRule = rules.at(j); 00681 if (currentRule.getAtomName() == currentAtom.getAtomName() && 00682 st.getTitGroupName() == currentAtom.getResidueName()){ 00683 float currentCharge = currentRule.getCharge(); 00684 // There should be a better way than this. 00685 if ( currentCharge == 0){ 00686 currentCharge = 0.01; 00687 } 00688 00689 currentAtom.setCharge(currentCharge); 00690 found = true; 00691 break; 00692 } 00693 } 00694 if (!found){ 00695 //currentAtom.setCharge(0); 00696 } 00697 //currentAtom.print(); 00698 atomList.at(i) = currentAtom; 00699 } 00700 } 00701 00702 00710 void writeOutW(string fileName){ 00711 cout << "Writing out W matrix .. " << endl; 00712 ofstream out( fileName.c_str() ); 00713 out.setf(ios_base::scientific); 00714 out.precision(6); 00715 00716 cout << "Writing W matrix into '" << fileName << "' ... "; 00717 cout.flush(); 00718 00719 if (!out) 00720 { 00721 out.close(); 00722 cout << "error writing: " << fileName << endl; 00723 exit(0); 00724 // throw FileIOException(filename, "PQR::writeOutW"); 00725 } 00726 00727 // run over all matrix elements 00728 for(unsigned int i=0; i<titGroupList.size(); i++) 00729 { 00730 for(unsigned int j=0; j<titGroupList.size(); j++) 00731 { 00732 for(unsigned int state1 = 1; state1<titGroupList.at(i).getNrStates(); state1++) 00733 { 00734 for(unsigned int state2 = 1; state2<titGroupList.at(j).getNrStates(); state2++) 00735 { 00736 out.width(4); 00737 out << i+1 << " " << state1+1 << " "; 00738 out.width(4); 00739 out << j+1 << " " << state2+1 << " "; 00740 out << wmatrix[i][j][state1-1][state2-1] << endl; 00741 } 00742 } 00743 } 00744 } 00745 00746 out.close(); 00747 00748 cout << "finished." << endl; 00749 00750 } 00751 00764 double calcWmv(int state, Residue &ref, Residue &mref, Residue &mstat, string cycle) 00765 { 00766 00767 // double wmv = 0; 00768 double wmv_c0 = 0; 00769 double wmv_c1 = 0; 00770 00771 Atom charge_s0, charge_s1; 00772 00773 vector<Atom> mstatAtomList = mstat.getAtomList(); 00774 vector<Atom> mrefAtomList = mref.getAtomList(); 00775 00776 tCycle pp_s1 = ref.getTCycle(); 00777 00778 // mstat size equals to mref size 00779 unsigned int i = 0; 00780 00781 00782 cycle = "cycle1"; 00783 while (mrefAtomList.size() > i){ 00784 Atom charge_s1 = mstatAtomList.at(i); 00785 Atom charge_s0 = mrefAtomList.at(i); 00786 i++; 00787 00788 if (charge_s0.getCharge() == 0 && charge_s1.getCharge() == 0 ) 00789 continue; 00790 00791 if (cycle == "cycle0" || cycle == "diff"){ 00792 wmv_c0 += (charge_s1.getCharge() * ( 00793 pp_s1.m[state][charge_s1.getCoord()] 00794 - pp_s1.m[1][charge_s1.getCoord()] 00795 ) 00796 - charge_s0.getCharge() * ( 00797 pp_s1.m[state][charge_s0.getCoord()] 00798 - pp_s1.m[1][charge_s0.getCoord()] 00799 ) 00800 ); 00801 } 00802 00803 if (cycle == "cycle1" || cycle == "diff"){ 00804 // cout << "Lookint at charge_s1 (" << charge_s1.getCoord().X() << ", " << charge_s1.getCoord().Y() << ", " << charge_s1.getCoord().Z() << ")" << endl; 00805 // cout << "Lookint at charge_s0 (" << charge_s0.getCoord().X() << ", " << charge_s0.getCoord().Y() << ", " << charge_s0.getCoord().Z() << ")" << endl; 00806 float temp = (charge_s1.getCharge() - charge_s0.getCharge()) * ( 00807 pp_s1.p[state][charge_s1.getCoord()] - pp_s1.p[0][charge_s1.getCoord()] 00808 ); 00809 wmv_c1 += (charge_s1.getCharge() - charge_s0.getCharge()) * ( 00810 pp_s1.p[state][charge_s1.getCoord()] - pp_s1.p[0][charge_s1.getCoord()] 00811 ); 00812 00813 // cout << "(" << charge_s1.getCharge() << " - " << charge_s0.getCharge() << ") * (" << pp_s1.p[state][charge_s1.getCoord()] << " - " << pp_s1.p[0][charge_s1.getCoord()] << ") = " << temp << endl; 00814 00815 } 00816 00817 } 00818 // cout << "wmv_c1 = " << wmv_c1 << ", wmv_c0 = " << wmv_c0 << endl; 00819 // cout << "returning wmv = " << wmv_c1 << endl; 00820 return (wmv_c1 - wmv_c0); 00821 } 00822 00834 void calcW(string cycle){ 00835 cout << "Calculation of W matrix started.. " << endl; 00836 double G; 00837 00838 // create wmatrix 00839 for(unsigned int i = 0; i<titGroupList.size(); i++) 00840 { 00841 wmatrix.push_back( vector< vector< vector<double> > >() ); 00842 for(unsigned int j = 0; j<titGroupList.size(); j++) 00843 { 00844 wmatrix.at(i).push_back( vector< vector<double> >() ); 00845 /* there are no matrix entries for reference states 00846 that's why we take only the number of states-1 */ 00847 for(unsigned int state1 = 0; state1<titGroupList.at(i).getNrStates()-1; state1++) 00848 { 00849 wmatrix.at(i).at(j).push_back( vector<double>() ); 00850 for(unsigned int state2 = 0; state2<titGroupList[j].getNrStates()-1; state2++) 00851 { 00852 wmatrix.at(i).at(j).at(state1).push_back(0.0); 00853 } 00854 } 00855 } 00856 } 00857 00858 // have to make diff of both cycles or just read one cycle 00859 for(unsigned int i = 0; i<titGroupList.size(); i++) 00860 { 00861 Residue currentTitGroup_i = titGroupList.at(i); 00862 // cout << "i has name: " << currentTitGroup_i.getResidueName() << endl; 00863 // read potentials of residue 1 00864 for(unsigned int j = 0; j<titGroupList.size(); j++) 00865 { 00866 // cout << "at i,j: " << i << ", " << j << endl; 00867 Residue currentTitGroup_j = titGroupList.at(j); 00868 Residue currentTitGroup_j_Ref = titGroupList.at(j); 00869 // cout << "j has name: " << currentTitGroup_j.getResidueName() << endl; 00870 // diagonal elements to be zero 00871 if ( currentTitGroup_i.getResidueName() == currentTitGroup_j.getResidueName() && i == j) 00872 { 00873 // cout << "both have same name: " << currentTitGroup_i.getResidueName() << endl; 00874 00875 for(unsigned int state1 = 1; state1<currentTitGroup_i.getNrStates(); state1++) 00876 { 00877 for(unsigned int state2 = 1; state2<currentTitGroup_j.getNrStates(); state2++) 00878 { 00879 wmatrix[i][j][state1-1][state2-1] = 0.0; 00880 } 00881 } 00882 continue; 00883 } 00884 00885 00886 // run over all non-reference states 00887 for(unsigned int state1 = 1; state1<currentTitGroup_i.getNrStates(); state1++) 00888 { 00889 for(unsigned int state2 = 1; state2<currentTitGroup_j.getNrStates(); state2++) 00890 { 00891 00892 // cout << "state1 is " << state1 << endl; 00893 // cout << "state2 is " << state2 << endl; 00894 00895 // get charge of current state in residue 2 00896 // Switch titGroup j to current state2 (!= (reference which is 0)) 00897 // Switch titGroup i to reference state 00898 switchState(currentTitGroup_i, 1); 00899 switchState(currentTitGroup_j, state2+1); 00900 switchState(currentTitGroup_j_Ref, 1); 00901 00902 // Wmv 00903 // G = calcWmv(state1, currentTitGroup_j, currentTitGroup_i, cycle); 00904 00905 // cout << "G cycle0 is " << calcWmv(state1, pstat, pref, "cycle0"); 00906 // cout << "G cycle1 is " << calcWmv(state1, pstat, pref, "cycle1"); 00907 00908 G = calcWmv(state1, currentTitGroup_i, currentTitGroup_j_Ref, currentTitGroup_j, cycle); // in kJ/mol 00909 00910 G *= E2A; 00911 00912 // cout << currentTitGroup_i.getResidueName() << " " << i << " with " << currentTitGroup_j.getResidueName() << " " << j ; 00913 00914 // cout << "G: cycle1 - cycl0 difference is " << G << endl; 00915 00916 // cout << "writing at " << i << ", " << j << ", " << state1-1 << ", " << state2-1 << endl; 00917 // update matrix 00918 wmatrix[i][j][state1-1][state2-1] = G; 00919 00920 } 00921 } 00922 } 00923 } 00924 cout << "finished." << endl; 00925 00926 // symmetrize matrix 00927 double dev=0.0, w1=0.0, w2=0.0, ave=0.0; 00928 00929 cout << "Symmetrizing W matrix ... "; 00930 00931 // setup output stream for symmetry check 00932 ostringstream log; 00933 log.setf(ios_base::scientific); 00934 log.precision(6); 00935 00936 for(unsigned int i=0; i<titGroupList.size(); i++) 00937 { 00938 for(unsigned int j=i+1; (j>i) && (j<titGroupList.size()); j++) 00939 { 00940 for(unsigned int state1 = 1; state1<titGroupList.at(i).getNrStates(); state1++) 00941 { 00942 for(unsigned int state2 = 1; state2<titGroupList.at(j).getNrStates(); state2++) 00943 { 00944 w1 = wmatrix[i][j][state1-1][state2-1]; 00945 w2 = wmatrix[j][i][state2-1][state1-1]; 00946 dev = (w1-w2); 00947 ave = 0.5*(w1+w2); 00948 wmatrix[i][j][state1-1][state2-1] = ave; 00949 wmatrix[j][i][state2-1][state1-1] = ave; 00950 00951 00952 // cout << titGroupList.at(i).getResidueName() << " " << i << " with " << titGroupList.at(j).getResidueName() << " " << j; 00953 // cout << "w1: " << w1 << ", w2: " << w2 << endl; 00954 00955 // write symmetry check to log file 00956 log << titGroupList.at(i).getResidueName() << "[" << state1+1 << "]/"; 00957 log << titGroupList.at(j).getResidueName() << "[" << state2+1 << "]: "; 00958 log << "(" << w1 << "," << w2 << "), diff = " << dev; 00959 log << ", ave = " << ave << endl; 00960 } 00961 } 00962 } 00963 } 00964 00965 cout << "finished." << endl; 00966 00967 // cout << log.str() << endl; 00968 } 00969 00970 00971 00981 void calcPkint(string cycleName, string refPQRFileName){ 00982 00983 for (unsigned int i = 0; i < titGroupList.size(); i++){ 00984 Residue currentTitGroup = titGroupList.at(i); 00985 if (cycleName == "cycle0" && (getCycle0Shift(titGroupList.at(i).getResidueName(), 1) != NOT_IN_ST && !explicitModels) ){ 00986 tCycle cycle; 00987 currentTitGroup.setTCycle(cycle); 00988 } else { 00989 00990 string fileName = cycleName+"."+currentTitGroup.getIdentifier()+".potat"; 00991 cout << "Reading " << fileName << " .... "; 00992 00993 tCycle cycle; 00994 unsigned int nrStates = 0; 00995 unsigned int nrAtoms = 0; 00996 float potential = 0; 00997 unsigned int potNr = 0; 00998 00999 ifstream in(fileName.c_str()); 01000 if (!in) { 01001 in.close(); 01002 cout << "Cannot open potat file:" << fileName << endl; 01003 exit(0); 01004 } 01005 string currentLine; 01006 while( !in.eof() ) { 01007 getline(in, currentLine); 01008 nrStates = atoi(currentLine.c_str()); 01009 getline(in, currentLine); 01010 nrAtoms = atoi(currentLine.c_str()); 01011 for(unsigned int j = 0; j < nrStates; j++){ 01012 for(unsigned int k = 0; k < nrAtoms; k++){ 01013 Point3<float> newPoint; 01014 getline(in, currentLine); 01015 newPoint.X() = atof(currentLine.c_str()); 01016 getline(in, currentLine); 01017 newPoint.Y() = atof(currentLine.c_str()); 01018 getline(in, currentLine); 01019 newPoint.Z() = atof(currentLine.c_str()); 01020 getline(in, currentLine); 01021 potential = atof(currentLine.c_str()); 01022 cycle.p[j].insert( {newPoint, potential} ); 01023 potNr++; 01024 } 01025 getline(in, currentLine); 01026 } 01027 nrStates = 0; 01028 nrAtoms = 0; 01029 nrStates = atoi(currentLine.c_str()); 01030 getline(in, currentLine); 01031 nrAtoms = atoi(currentLine.c_str()); 01032 for(unsigned int i = 0; i < nrStates; i++){ 01033 for(unsigned int j = 0; j < nrAtoms; j++){ 01034 Point3<float> newPoint; 01035 getline(in, currentLine); 01036 newPoint.X() = atof(currentLine.c_str()); 01037 getline(in, currentLine); 01038 newPoint.Y() = atof(currentLine.c_str()); 01039 getline(in, currentLine); 01040 newPoint.Z() = atof(currentLine.c_str()); 01041 getline(in, currentLine); 01042 potential = atof(currentLine.c_str()); 01043 cycle.m[i].insert( {newPoint, potential} ); 01044 potNr++; 01045 } 01046 getline(in, currentLine); 01047 01048 } 01049 cout << potNr << " potentials read in." << endl; 01050 01051 } 01052 currentTitGroup.setTCycle(cycle); 01053 titGroupList.at(i) = currentTitGroup; 01054 in.close(); 01055 01056 } 01057 } 01058 01059 vector<Atom> refAllAtomList; 01060 for (unsigned int i = 0; i < residueList.size(); i++){ 01061 Residue currentResidue = residueList.at(i); 01062 01063 if (i == 0 && calcNTE){ 01064 Residue nteResidue = residueList.at(0); 01065 vector<Atom> list = nteResidue.getAtomList(); 01066 for (unsigned int k = 0; k < list.size(); k++){ 01067 Atom currentAtom = list.at(k); 01068 currentAtom.setResidueName("NTE"); 01069 list.at(k) = currentAtom; 01070 } 01071 nteResidue.setResidueName("NTE"); 01072 nteResidue.setAtomList(list); 01073 neutralState(nteResidue, 1); 01074 list = nteResidue.getAtomList(); 01075 for (unsigned int k = 0; k < list.size(); k++){ 01076 Atom currentAtom = list.at(k); 01077 currentAtom.setResidueName(currentResidue.getResidueName()); 01078 list.at(k) = currentAtom; 01079 } 01080 currentResidue.setAtomList(list); 01081 } else if ( i == residueList.size()-1 && calcCTE ){ 01082 Residue cteResidue = residueList.at(residueList.size()-1); 01083 vector<Atom> list = cteResidue.getAtomList(); 01084 for (unsigned int k = 0; k < list.size(); k++){ 01085 Atom currentAtom = list.at(k); 01086 currentAtom.setResidueName("CTE"); 01087 list.at(k) = currentAtom; 01088 } 01089 cteResidue.setResidueName("CTE"); 01090 cteResidue.setAtomList(list); 01091 neutralState(cteResidue, 1); 01092 list = cteResidue.getAtomList(); 01093 for (unsigned int k = 0; k < list.size(); k++){ 01094 Atom currentAtom = list.at(k); 01095 currentAtom.setResidueName(currentResidue.getResidueName()); 01096 list.at(k) = currentAtom; 01097 } 01098 currentResidue.setAtomList(list); 01099 01100 } 01101 01102 bool found = false; 01103 for (unsigned int j = 0; j < titGroupList.size(); j++){ 01104 Residue currentTGroup = titGroupList.at(j); 01105 if (currentTGroup.getResidueName() == currentResidue.getResidueName() && currentTGroup.getResidueNumber() == currentResidue.getResidueNumber()){ 01106 found = true; 01107 cout << "reference: setting " << currentTGroup.getIdentifier() << " to neutral state" << endl; 01108 01109 break; 01110 } 01111 } 01112 if (found){ 01113 neutralState(currentResidue, 1); 01114 } 01115 vector<Atom> residueAtoms = currentResidue.getAtomList(); 01116 refAllAtomList.insert(refAllAtomList.end(), residueAtoms.begin(), residueAtoms.end()); 01117 } 01118 01119 writePQR(refAllAtomList, refPQRFileName); 01120 01121 01122 float bornEner = 0; 01123 float backEner = 0; 01124 01125 for (unsigned int i = 0; i < titGroupList.size(); i++){ 01126 Residue currentOrigTitGroup = titGroupList.at(i); 01127 Residue currentTitGroup = titGroupList.at(i); 01128 unsigned int nrStates = currentTitGroup.getNrStates(); 01129 01130 01131 for (unsigned int stateNr = 1; stateNr < nrStates; stateNr++){ 01132 tCycle cycle = currentTitGroup.getTCycle(); 01133 01134 vector<Atom> refAtomList = getAtoms(currentTitGroup, 0); 01135 vector<Atom> compareAtomList = getAtoms(currentTitGroup, stateNr); 01136 01137 float born0 = 0; 01138 float born1 = 0; 01139 01140 for (unsigned int j = 0; j < refAtomList.size(); j++){ 01141 Atom currentRefAtom = refAtomList.at(j); 01142 Atom currentCompareAtom = compareAtomList.at(j); 01143 if (currentRefAtom.getCharge() != 0){ 01144 born0 += currentRefAtom.getCharge() * (cycle.p[0][currentRefAtom.getCoord()] - cycle.m[0][currentRefAtom.getCoord()]); 01145 } 01146 01147 if (currentCompareAtom.getCharge() != 0){ 01148 born1 += currentCompareAtom.getCharge() * (cycle.p[stateNr][currentCompareAtom.getCoord()] - cycle.m[stateNr][currentCompareAtom.getCoord()]); 01149 } 01150 } 01151 bornEner = 0.5*(born1-born0); 01152 01153 float back0 = 0; 01154 float back1 = 0; 01155 01156 // cout << "Calculating BACK" << endl; 01157 // cout << "================" << endl; 01158 01159 //Residue titGroup = titGroup ! 01160 // vector<Atom> pAtomList = deleteAtoms(refAllAtomList, refAtomList); // from whole protein without patched atoms in currentTitGroup 01161 vector<Atom> pAtomList = deleteAtoms(refAllAtomList, refAtomList); 01162 stringstream ss; 01163 ss << "pAtomList_" << currentTitGroup.getIdentifier() << ".pqr"; 01164 writePQR(pAtomList, ss.str()); 01165 01166 01167 vector<Atom> mAtomList = deleteAtoms(currentOrigTitGroup.getAtomList(), refAtomList); // not patched 01168 stringstream ss2; 01169 ss2 << "mAtomList_" << currentTitGroup.getIdentifier() << ".pqr"; 01170 writePQR(mAtomList, ss2.str()); 01171 01172 01173 for (unsigned int j = 0; j < pAtomList.size(); j++){ 01174 Atom currentAtom = pAtomList.at(j); 01175 if (currentAtom.getCharge() != 0){ 01176 back1 += currentAtom.getCharge() * (cycle.p[0][currentAtom.getCoord()] - cycle.p[stateNr][currentAtom.getCoord()]); 01177 if (cycle.p[0][currentAtom.getCoord()] != 0){ 01178 Point3<float> c = currentAtom.getCoord(); 01179 // cout << "Looking at pcharge (" << c.X() << ", " << c.Y() << ", " << c.Z() << endl; 01180 // cout << "back1 = " << currentAtom.getCharge() << " * ( " << cycle.p[0][currentAtom.getCoord()] << " - " << cycle.p[stateNr][currentAtom.getCoord()] << " ) = " << back1 << endl; 01181 } 01182 } 01183 } 01184 01185 for (unsigned int j = 0; j < mAtomList.size(); j++){ 01186 Atom currentAtom = mAtomList.at(j); 01187 if (currentAtom.getCharge() != 0){ 01188 back0 += currentAtom.getCharge() * (cycle.m[0][currentAtom.getCoord()] - cycle.m[stateNr][currentAtom.getCoord()]); 01189 if (cycle.m[0][currentAtom.getCoord()] != 0){ 01190 Point3<float> c = currentAtom.getCoord(); 01191 01192 // cout << "Looking at mcharge (" << c.X() << ", " << c.Y() << ", " << c.Z() << endl; 01193 // cout << "back2 = " << currentAtom.getCharge() << " * ( " << cycle.m[0][currentAtom.getCoord()] << " - " << cycle.m[stateNr][currentAtom.getCoord()] << " ) = " << back0 << endl; 01194 } 01195 } 01196 } 01197 backEner = (back1-back0); 01198 // cout << "returning (" << back1 << "-" << back0 << ") =" << backEner << endl; 01199 01200 unsigned int cycleNr = 0; 01201 if (cycleName == "cycle1") 01202 cycleNr = 1; 01203 01204 currentTitGroup.setBorn(bornEner, cycleNr, stateNr); // real titGroupList object 01205 currentTitGroup.setBack(backEner, cycleNr, stateNr); 01206 titGroupList.at(i) = currentTitGroup; 01207 } 01208 } 01209 } 01210 01227 vector<Atom> deleteAtoms(vector<Atom> result, vector<Atom>& compareAtomList){ 01228 // compareAtomList is << currentAtomList 01229 for (unsigned int i = 0; i < compareAtomList.size(); i++){ 01230 Atom currentCompareAtom = compareAtomList.at(i); 01231 // currentCompareAtom.print(); 01232 if (currentCompareAtom.getCharge() != 0){ 01233 unsigned int pos = 0; 01234 for (unsigned int j = 0; j < result.size(); j++){ 01235 Atom currentAtom = result.at(j); 01236 if (currentAtom.getAtomNumber() == currentCompareAtom.getAtomNumber()){ 01237 result.erase(result.begin()+j-pos); 01238 pos++; 01239 } 01240 } 01241 } 01242 } 01243 return result; 01244 } 01257 vector<Atom> getAtoms(Residue currentTitGroup, int stateNr){ 01258 vector<Atom> currentAtomList = currentTitGroup.getAtomList(); 01259 01260 for (unsigned int i = 0; i < stList.size(); i++){ 01261 ST currentST = stList.at(i); 01262 if (currentST.getStateNr() == stateNr+1 && currentST.getTitGroupName() == currentTitGroup.getResidueName()){ 01263 patch(currentAtomList, currentST); 01264 } 01265 } 01266 01267 return currentAtomList; 01268 } 01269 01270 01280 float getCycle0Shift(string residueName, int state){ 01281 float shift = NOT_IN_ST; 01282 for (unsigned int i=0; stList.size(); i++){ 01283 ST currentST = stList.at(i); 01284 if (currentST.getTitGroupName() == residueName && currentST.getStateNr() == state){ 01285 shift = currentST.getShiftCycle0(); 01286 return currentST.getShiftCycle0(); 01287 } 01288 } 01289 return shift; 01290 } 01291 01301 void calcPotat(string cycleName){ 01302 string pdeFile; 01303 try { 01304 if (cycleName == "cycle1"){ 01305 ngsolve::PDE pde; 01306 bool skip = false; 01307 pdeFile = "pka_"+cycleName+".pde"; 01308 boost::filesystem::wpath file(pdeFile); 01309 01310 if(boost::filesystem::exists(file) && !skip){ 01311 pde.LoadPDE (pdeFile.c_str()); 01312 pde.Solve(); 01313 } else { 01314 cout << "cycle1 computation skipped manually." << endl; } 01315 01316 01317 } else { 01318 for (unsigned int i = 0; i < titGroupList.size(); i++){ 01319 ngsolve::PDE pde; 01320 if (getCycle0Shift(titGroupList.at(i).getResidueName(), 1) == NOT_IN_ST || explicitModels){ 01321 string prefix = titGroupList.at(i).getIdentifier(); 01322 pdeFile = "pka_"+cycleName+"_"+prefix+".pde"; 01323 boost::filesystem::wpath file(pdeFile); 01324 string potatFile = cycleName+"."+prefix+".potat"; 01325 boost::filesystem::wpath potfile(potatFile); 01326 01327 if(boost::filesystem::exists(file)){ 01328 //if(boost::filesystem::exists(file) 01329 // && !boost::filesystem::exists(potfile)){ 01330 01331 pde.LoadPDE (pdeFile.c_str()); 01332 pde.Solve(); 01333 } 01334 // else { 01335 // cout << potatFile << ".. already calculated." << endl; 01336 //} 01337 } 01338 } 01339 } 01340 01341 } catch(ngstd::Exception & e) { 01342 std::cout << "Caught exception: " << std::endl 01343 << e.What() << std::endl; 01344 } 01345 } 01346 01347 01357 void writePDE(INI &ini, string mode = "explicit"){ 01358 PDE currentPDE; 01359 if (mode == "explicit"){ 01360 currentPDE.writeCycle0(titGroupList, ini); 01361 currentPDE.writeCycle1(fileName, titGroupList, ini); 01362 } else { 01363 currentPDE.writeEnergy(fileName, ini); 01364 } 01365 } 01366 01379 void writeOutPkint(string cycleName, string fileName, INI &ini){ 01380 string eps_in = ini.get<string>("experiment.eps_in"); 01381 string eps_out = ini.get<string>("experiment.eps_out"); 01382 bool eps_same = false; 01383 if (eps_in == eps_out) 01384 eps_same = true; 01385 01386 ofstream out( fileName.c_str() ); 01387 out.setf(ios_base::scientific); 01388 out.precision(6); 01389 01390 if (!out) 01391 { 01392 out.close(); 01393 cout << "could not open: " << fileName << endl; 01394 } 01395 01396 01397 unsigned int cycleNr = 0; 01398 if (cycleName == "cycle1") 01399 cycleNr = 1; 01400 01401 for (unsigned int i = 0; i < titGroupList.size(); i++){ 01402 Residue currentTitGroup = titGroupList.at(i); 01403 // write reference energy [kJ/mol] to pkint file 01404 out << 0.0 << " " << getState(currentTitGroup.getResidueName(), 0) << " "; 01405 01406 unsigned int nrStates = currentTitGroup.getNrStates(); 01407 for (unsigned int stateNr = 1; stateNr < nrStates; stateNr++){ 01408 float born, back, shift, result; 01409 float cycle0Shift = getCycle0Shift(currentTitGroup.getResidueName(), stateNr+1); 01410 if ( cycle0Shift != NOT_IN_ST && !explicitModels){ 01411 if (cycleName == "cycle0"){ 01412 cout << "using precalculated values." << endl; 01413 cycle0Shift *= CONVERT; 01414 born = 0; //currentTitGroup.getBorn(cycleNr, stateNr); 01415 back = 0; //currentTitGroup.getBack(cycleNr, stateNr); 01416 shift = getShift(currentTitGroup.getResidueName(), stateNr)*CONVERT; 01417 if (eps_same) 01418 result = 0; 01419 else 01420 result = cycle0Shift + ( born - back ); 01421 01422 cout << currentTitGroup.getIdentifier() << ": R -> " << getState(currentTitGroup.getResidueName(), stateNr) << endl; 01423 cout << "Gborn: " << born << endl; 01424 cout << "Gback: " << back << endl; 01425 cout << "intrinsic pka = "<< shift << " + " << (cycle0Shift) << " + ( " << born << " - " << back << ") = " << result << " [kJ/mol]" << endl; 01426 cout << "-> " << result/CONVERT << " []" << endl; 01427 01428 } else if (cycleName == "cycle1") { 01429 born = currentTitGroup.getBorn(cycleNr, stateNr); 01430 back = currentTitGroup.getBack(cycleNr, stateNr); 01431 shift = 0; 01432 result = shift + ( born - back ); 01433 01434 cout << currentTitGroup.getIdentifier() << ": R -> " << getState(currentTitGroup.getResidueName(), stateNr) << endl; 01435 cout << "Gborn: " << born << endl; 01436 cout << "Gback: " << back << endl; 01437 cout << "intrinsic pka = "<< shift << " + ( " << born << " - " << back << ") = " << result << " [kJ/mol]" << endl; 01438 cout << "-> " << result/CONVERT << " []" << endl; 01439 01440 } else { 01441 cout << "using precalculated values." << endl; 01442 cycle0Shift *= CONVERT; 01443 born = currentTitGroup.getBorn(1, stateNr) - 0; //currentTitGroup.getBorn(0, stateNr); 01444 back = currentTitGroup.getBack(1, stateNr) - 0; //currentTitGroup.getBack(0, stateNr); 01445 shift = getShift(currentTitGroup.getResidueName(), stateNr)*CONVERT; 01446 if (eps_same) 01447 result = 0; 01448 else 01449 result = shift+(born - back)-(cycle0Shift); 01450 01451 cout << currentTitGroup.getIdentifier() << ": R -> " << getState(currentTitGroup.getResidueName(), stateNr) << endl; 01452 cout << "Gborn: " << born << endl; 01453 cout << "Gback: " << back << endl; 01454 cout << "intrinsic pka = "<< shift << " + ( " << born << " - " << back << ") - " << (cycle0Shift) << " = " << result << " [kJ/mol]" << endl; 01455 cout << "-> " << result/CONVERT << " []" << endl; 01456 } 01457 01458 // write pkint [kJ/mol] and transition identifier to pkint file 01459 out << result << " " << getState(currentTitGroup.getResidueName(), stateNr) << " "; 01460 01461 } else { 01462 01463 if (cycleName == "cycle0"){ 01464 born = currentTitGroup.getBorn(cycleNr, stateNr); 01465 back = currentTitGroup.getBack(cycleNr, stateNr); 01466 shift = 0; 01467 result = shift + ( born - back ); 01468 } else if (cycleName == "cycle1") { 01469 born = currentTitGroup.getBorn(cycleNr, stateNr); 01470 back = currentTitGroup.getBack(cycleNr, stateNr); 01471 shift = 0; 01472 result = shift + ( born - back ); 01473 } else { 01474 born = currentTitGroup.getBorn(1, stateNr) - currentTitGroup.getBorn(0, stateNr); 01475 back = currentTitGroup.getBack(1, stateNr) - currentTitGroup.getBack(0, stateNr); 01476 shift = getShift(currentTitGroup.getResidueName(), stateNr); 01477 shift *= CONVERT; 01478 result = shift+(born - back); 01479 } 01480 cout << currentTitGroup.getIdentifier() << ": R -> " << getState(currentTitGroup.getResidueName(), stateNr) << endl; 01481 cout << "Gborn: " << born << endl; 01482 cout << "Gback: " << back << endl; 01483 cout << "intrinsic pka = "<< shift << " + ( " << born << " - " << back << ") = " << result << " [kJ/mol]" << endl; 01484 cout << "-> " << result/CONVERT << " []" << endl; 01485 // write pkint [kJ/mol] and transition identifier to pkint file 01486 out << result << " " << getState(currentTitGroup.getResidueName(), stateNr) << " "; 01487 01488 } 01489 01490 } 01491 01492 out << currentTitGroup.getIdentifier() << endl; 01493 01494 } 01495 01496 } 01497 01507 char getState(string residueName, int stateNr){ 01508 for (unsigned int i = 0; i < stList.size(); i++){ 01509 ST currentST = stList.at(i); 01510 if (currentST.getTitGroupName() == residueName && currentST.getStateNr() == (stateNr+1)){ 01511 return currentST.getState(); 01512 } 01513 } 01514 return '?'; 01515 } 01516 01526 float getShift(string residueName, int stateNr){ 01527 for (unsigned int i = 0; i < stList.size(); i++){ 01528 ST currentST = stList.at(i); 01529 if (currentST.getTitGroupName() == residueName && currentST.getStateNr() == (stateNr+1)){ 01530 return currentST.getShift(); 01531 } 01532 } 01533 return 0; 01534 } 01535 01548 float getShiftCycle0(string residueName, int stateNr){ 01549 for (unsigned int i = 0; i < stList.size(); i++){ 01550 ST currentST = stList.at(i); 01551 if (currentST.getTitGroupName() == residueName && currentST.getStateNr() == (stateNr+1)){ 01552 return currentST.getShiftCycle0(); 01553 } 01554 } 01555 return 0; 01556 } 01557 01559 static bool STparsed; 01561 static bool calcNTE; 01563 static bool calcCTE; 01565 static bool explicitModels; 01566 01567 01568 private: 01570 string fileName; 01573 string stFilesFolder; 01575 Model molecule; 01577 vector<Atom> atomList; 01579 vector<Residue> residueList; 01582 vector<Residue> titGroupList; 01584 vector<ST> stList; 01586 tCycle cycle; 01588 vector< vector< vector< vector<double> > > > wmatrix; 01589 }; 01590 01591 bool PQR::STparsed = false; 01592 bool PQR::explicitModels = false; 01593 bool PQR::calcCTE = false; 01594 bool PQR::calcNTE = false;
1.7.6.1