1 #include "buildTopology.h" 2 3 #include "ExclusionTable.h" 4 #include "ExclusionType.h" 5 #include "GenericTopology.h" 6 #include "LennardJonesParameterTable.h" 7 #include "PAR.h" 8 #include "PSF.h" 9 #include "Report.h" 10 #include "pmconstants.h" 11 #include "mathutilities.h" 12 #include "stringutilities.h" 13 #include "topologyutilities.h" 14 //#include "Molecule.h" 15 //#include "Vector3D.h" 16 17 #include <algorithm> 18 #include <list> 19 #include <map> 20 #include <set> 21 #include <string> 22 #include <vector> 23 24 using std::list; 25 using std::map; 26 using std::set; 27 using std::pair; 28 using std::string; 29 using std::vector; 30 using std::swap; 31 32 using namespace ProtoMol::Report; 33 34 //#define DEBUG_PRINT_MOLECULETABLE 35 36 namespace ProtoMol { 37 38 static void findNextNeighbor(int a, vector<int>& v, vector<PairInt>& p, vector<bool>& unused, const vector<vector<int> >& graph, set<PairInt>& pairs); 39 // bool cmpSize(const vector<int>& m1, const vector<int>& m2); 40 41 42 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 43 // 44 // buildTopology 45 // 46 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 47 buildTopology(GenericTopology * topo,const PSF & psf,const PAR & par,bool dihedralMultPSF)48 void buildTopology(GenericTopology* topo,const PSF& psf, const PAR& par, bool dihedralMultPSF){ 49 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 50 // First, generate the array of atomtypes 51 // Each time a new atom comes up, we need to check if it is 52 // already in the vector.... 53 // NOTE: this may take a while for large systems; however, it will cut 54 // down on the size of the atomTypes vector, and therefore, the amount 55 // access time in the back end. 56 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 57 58 topo->atoms.clear(); 59 topo->atomTypes.clear(); 60 topo->bonds.clear(); 61 topo->angles.clear(); 62 topo->dihedrals.clear(); 63 topo->impropers.clear(); 64 65 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 66 // Get the atoms 67 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 68 69 map<string,int> atomLookUpTable; 70 71 // loop over all atoms in the PSF object 72 for (vector<PSF::Atom>::const_iterator atom = psf.atoms.begin(); 73 atom != psf.atoms.end(); 74 ++atom) { 75 // Two data members for AtomType, name and mass 76 AtomType tempatomtype; 77 tempatomtype.name = atom->atom_type; 78 tempatomtype.mass = atom->mass; 79 // tempatomtype.symbolName = charmmDefaults->getNonbondedData(atom->atom_type).mySymbolName; 80 tempatomtype.symbolName = atomTypeToSymbolName(atom->atom_type); 81 tempatomtype.charge = atom->charge; 82 // Now check if this already exists (same name) 83 if (atomLookUpTable.find(tempatomtype.name) == atomLookUpTable.end()) { 84 atomLookUpTable[tempatomtype.name] = topo->atomTypes.size(); 85 topo->atomTypes.push_back(tempatomtype); 86 } 87 88 Atom tempatom; 89 // First, we need to find the index. (an integer corresponding 90 // to the type of the atom 91 tempatom.type = atomLookUpTable[tempatomtype.name]; 92 // Now, the scaled charge. This is straightforward. 93 tempatom.scaledCharge = (atom->charge)*Constant::SQRTCOULOMBCONSTANT; 94 tempatom.scaledMass = atom->mass; 95 96 // Now we need the size of the group for heavy atom ordering 97 // We need to parse the name for any H's then any numbers following 98 // First, if the atom is an H then this is 0 99 if (atom->atom_type == "H"){ 100 tempatom.hvyAtom = 0; 101 } 102 else{ 103 // Otherwise, we need to parse.. 104 // Initialize to 1 105 tempatom.hvyAtom = 1; 106 for (unsigned int pos = 0; pos < atom->atom_type.size(); ++pos){ 107 if (atom->atom_type[pos] == 'H'){ 108 string number = ""; 109 while (isdigit(atom->atom_type[++pos])) { 110 number += atom->atom_type[pos]; 111 } 112 if (number == "") // never entered loop, default is 1 113 number = "1"; 114 tempatom.hvyAtom += atoi(number.c_str()); 115 } 116 } 117 } 118 // C/C++ starts at 0, where PSF/PDB at 1 119 tempatom.atomNum = atom->number-1; 120 // Also the molecule - using residue sequence for now 121 topo->atoms.push_back(tempatom); 122 123 124 } 125 126 // calculate the # of degrees of freedom, if there are any bond constraints 127 // they will be subtracted later by ModifierShake 128 topo->degreesOfFreedom = 3 * topo->atoms.size() - 3; 129 130 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 131 // Get the bonds 132 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 133 134 // First create look-up-table 135 map<string,vector<PAR::Bond>::const_iterator> bondLookUpTable; 136 for (vector<PAR::Bond>::const_iterator bond = par.bonds.begin(); 137 bond != par.bonds.end(); 138 ++bond){ 139 bondLookUpTable[bond->atom1+","+bond->atom2] = bond; 140 //report << (*bond)<< ", " << bond->atom1 << ", " << bond->atom2 <<endr; 141 } 142 143 144 // Find the parameters from PAR 145 int ignoredBonds = 0; 146 for (vector<PSF::Bond>::const_iterator bond = psf.bonds.begin(); 147 bond != psf.bonds.end(); 148 ++bond){ 149 150 // store the ID numbers of the bonded atoms 151 int atom1 = bond->atom1-1; 152 int atom2 = bond->atom2-1; 153 154 // store the type names of the bonded atoms 155 string bond1(topo->atomTypes[topo->atoms[atom1].type].name); 156 string bond2(topo->atomTypes[topo->atoms[atom2].type].name); 157 158 map<string,vector<PAR::Bond>::const_iterator>::const_iterator currentbond = bondLookUpTable.find(bond1+","+bond2); 159 if(currentbond == bondLookUpTable.end()){ 160 currentbond = bondLookUpTable.find(bond2+","+bond1); 161 } 162 163 // if we still have not found this bond type in the PAR object, report an error 164 if(currentbond == bondLookUpTable.end()){ 165 report << error << "Could not find bond \'"<<bond1<<"\'-\'"<<bond2<<"\' ("<<bond->atom1 <<","<< bond->atom2<<")"<<std::endl; 166 for (map<string,vector<PAR::Bond>::const_iterator>::const_iterator i = bondLookUpTable.begin(); 167 i != bondLookUpTable.end(); 168 ++i){ 169 report << plain << i->first<<std::endl; 170 } 171 report << endr; 172 } 173 174 // if we have found this bond type then copy the bond parameters 175 // into the topology 176 Bond tempbond; 177 tempbond.springConstant = currentbond->second->forceConstant; 178 tempbond.restLength = currentbond->second->distance; 179 tempbond.atom1 = atom1; 180 tempbond.atom2 = atom2; 181 topo->bonds.push_back(tempbond); 182 183 // populate the vector of bonds maintained at each atom 184 //report << std::endl <<"Size of Bonds Vector "<<topo->bonds.size() <<endr; 185 topo->atoms[atom1].mybonds.push_back((topo->bonds.size())-1); 186 topo->atoms[atom2].mybonds.push_back((topo->bonds.size())-1); 187 // output to screen for testing purposes 188 //for (int j = 0; j < topo->atoms[atom1].mybonds.size(); j++){ 189 //report <<"Atom " << atom1 << " bond index = "<< topo->atoms[atom1].mybonds[j] <<endr; 190 //} 191 //for (int k = 0; k < topo->atoms[atom2].mybonds.size(); k++){ 192 //report <<"Atom " << atom2 << " bond index = "<< topo->atoms[atom2].mybonds[k] <<endr; 193 //} 194 195 if(tempbond.springConstant == 0.0) 196 ++ignoredBonds; 197 } 198 199 if(ignoredBonds > 0) 200 report << hint << "Systems contains "<<ignoredBonds<<" bonds with zero force constants."<<endr; 201 202 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 203 // Get the angles 204 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 205 206 // First create look-up-table 207 map<string,vector<PAR::Angle>::const_iterator> angleLookUpTable; 208 for (vector<PAR::Angle>::const_iterator angle = par.angles.begin(); 209 angle != par.angles.end(); 210 ++angle){ 211 angleLookUpTable[angle->atom1+","+angle->atom2+","+angle->atom3] = angle; 212 //report << (*angle)<<endr; 213 } 214 215 // Find the parameters from PAR 216 int ignoredAngles = 0; 217 218 // loop over the angle list in the PSF object 219 for (vector<PSF::Angle>::const_iterator angle = psf.angles.begin(); 220 angle != psf.angles.end(); ++angle){ 221 222 // store the ID numbers of the atoms in this angle 223 int atom1 = angle->atom1-1; 224 int atom2 = angle->atom2-1; 225 int atom3 = angle->atom3-1; 226 227 // store the type names of the atoms in this angle 228 string angle1(topo->atomTypes[topo->atoms[atom1].type].name); 229 string angle2(topo->atomTypes[topo->atoms[atom2].type].name); 230 string angle3(topo->atomTypes[topo->atoms[atom3].type].name); 231 232 map<string,vector<PAR::Angle>::const_iterator>::const_iterator currentangle = angleLookUpTable.find(angle1+","+angle2+","+angle3); 233 if(currentangle == angleLookUpTable.end()) 234 currentangle = angleLookUpTable.find(angle3+","+angle2+","+angle1); 235 236 // if we still have not found this angle type in the PAR object, report an error 237 if(currentangle == angleLookUpTable.end()) 238 report << error << "Could not find angle \'"<<angle1<<"\'-\'"<<angle2<<"\'-\'"<<angle3<<"\'."<<endr; 239 240 // if we have found this angle type then copy the angle parameters 241 // into the topology 242 Angle tempangle; 243 tempangle.atom1 = atom1; 244 tempangle.atom2 = atom2; 245 tempangle.atom3 = atom3; 246 tempangle.forceConstant = currentangle->second->forceConstant; 247 tempangle.restAngle = dtor(currentangle->second->angleval); 248 if (currentangle->second->ub_flag){ // do we want defaults for these 249 tempangle.ureyBradleyConstant = currentangle->second->k_ub; 250 tempangle.ureyBradleyRestLength = currentangle->second->r_ub; 251 } 252 // no Urey-Bradley term specified 253 else{ 254 tempangle.ureyBradleyConstant = 0.0; 255 tempangle.ureyBradleyRestLength = 0.0; 256 } 257 topo->angles.push_back(tempangle); 258 if(tempangle.forceConstant == 0.0) 259 ++ignoredAngles; 260 } 261 262 if(ignoredAngles > 0) 263 report << hint << "Systems contains "<<ignoredAngles<<" angles with zero force constants."<<endr; 264 265 266 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 267 // Get the dihedrals 268 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 269 // 270 // One change I made was to assume that a dihedral will only appear 271 // once in the .psf file regardless of it's multiplicity. The 272 // multiplicity should be handled in the .par file. 273 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 274 275 276 // First create look-up-table 277 map<string,vector<PAR::Dihedral>::const_iterator> dihedralLookUpTable; 278 for (vector<PAR::Dihedral>::const_iterator dihedral = par.dihedrals.begin(); 279 dihedral != par.dihedrals.end(); 280 ++dihedral){ 281 dihedralLookUpTable[dihedral->atom1+","+dihedral->atom2+","+dihedral->atom3+","+dihedral->atom4] = dihedral; 282 //report << (*dihedral)<<endr; 283 } 284 285 // Find the parameters from PAR 286 // loop over the dihedral list in the PSF object 287 for(vector<PSF::Dihedral>::const_iterator dihedral = psf.dihedrals.begin(); 288 dihedral != psf.dihedrals.end(); ++dihedral) { 289 290 // store the ID numbers of the atoms in this dihedral 291 int atom1 = dihedral->atom1 - 1; 292 int atom2 = dihedral->atom2 - 1; 293 int atom3 = dihedral->atom3 - 1; 294 int atom4 = dihedral->atom4 - 1; 295 296 // store the type names of the atoms in this dihedral 297 string dihedral1 = topo->atomTypes[topo->atoms[atom1].type].name; 298 string dihedral2 = topo->atomTypes[topo->atoms[atom2].type].name; 299 string dihedral3 = topo->atomTypes[topo->atoms[atom3].type].name; 300 string dihedral4 = topo->atomTypes[topo->atoms[atom4].type].name; 301 302 map<string,vector<PAR::Dihedral>::const_iterator>::const_iterator currentdihedral = 303 dihedralLookUpTable.find(dihedral1+","+dihedral2+","+dihedral3+","+dihedral4); 304 305 // if this dihedral type has not been found, try reversing the order of the atom types 306 if(currentdihedral == dihedralLookUpTable.end()) 307 currentdihedral = dihedralLookUpTable.find(dihedral4+","+dihedral3+","+dihedral2+","+dihedral1); 308 309 // Try wildcards if necessary 310 if(currentdihedral == dihedralLookUpTable.end()){ 311 currentdihedral = dihedralLookUpTable.find("X,"+dihedral2+","+dihedral3+",X"); 312 if(currentdihedral == dihedralLookUpTable.end()) 313 currentdihedral = dihedralLookUpTable.find("X,"+dihedral3+","+dihedral2+",X"); 314 } 315 316 // if we still have not found this dihedral type in the PAR object, report an error 317 if(currentdihedral == dihedralLookUpTable.end()) 318 report << error << "Could not find dihedral \'"<<dihedral1<<"\'-\'"<<dihedral2<<"\'-\'"<<dihedral3<<"\'-\'"<<dihedral4<<"\'."<<endr; 319 320 // if we have found this dihedral type then copy the 321 // dihedral parameters into the topology 322 Torsion torsion; 323 torsion.atom1 = atom1; 324 torsion.atom2 = atom2; 325 torsion.atom3 = atom3; 326 torsion.atom4 = atom4; 327 328 torsion.periodicity = currentdihedral->second->periodicity; 329 torsion.forceConstant = currentdihedral->second->forceConstant; 330 torsion.phaseShift = dtor(currentdihedral->second->phaseShift); 331 torsion.multiplicity = currentdihedral->second->multiplicity; 332 if(topo->dihedrals.empty() || 333 topo->dihedrals[topo->dihedrals.size()-1].atom1 != atom1 || 334 topo->dihedrals[topo->dihedrals.size()-1].atom2 != atom2 || 335 topo->dihedrals[topo->dihedrals.size()-1].atom3 != atom3 || 336 topo->dihedrals[topo->dihedrals.size()-1].atom4 != atom4){ 337 if(dihedralMultPSF){ 338 torsion.periodicity.resize(1); 339 torsion.forceConstant.resize(1); 340 torsion.phaseShift.resize(1); 341 torsion.multiplicity = 1; 342 } 343 topo->dihedrals.push_back(torsion); 344 } 345 else { 346 if(dihedralMultPSF){ 347 Torsion& tmp = topo->dihedrals[topo->dihedrals.size()-1]; 348 if(tmp.multiplicity > torsion.multiplicity ) 349 report << error<< "PSF multiplicity definition of dihedral (" 350 << dihedral1 <<"," 351 << dihedral2 <<"," 352 << dihedral3 <<"," 353 << dihedral4 <<") exceeded PAR definition."; 354 tmp.periodicity.push_back(torsion.periodicity[tmp.multiplicity]); 355 tmp.forceConstant.push_back(torsion.forceConstant[tmp.multiplicity]); 356 tmp.phaseShift.push_back(torsion.phaseShift[tmp.multiplicity]); 357 tmp.multiplicity++; 358 } 359 else { 360 report << error << "Unexpected PSF multiplicity definition of dihedral (" 361 << dihedral1 <<"," 362 << dihedral2 <<"," 363 << dihedral3 <<"," 364 << dihedral4 <<") occurred, use dihedral multiplicity PSF."; 365 } 366 } 367 } 368 369 370 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 371 // Get the impropers 372 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 373 // 374 // One change I made was to assume that a improper will only appear 375 // once in the .psf file regardless of it's multiplicity. The 376 // multiplicity should be handled in the .par file. 377 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 378 379 // No wildcard usage is allowed for bonds and angles. For dihedrals, 380 // two types are allowed; A - B - C - D (all four atoms specified) and 381 // X - A - B - X (only middle two atoms specified). Double dihedral 382 // specifications may be specified for the four atom type by listing a 383 // given set twice. When specifying this type in the topology file, specify 384 // a dihedral twice (with nothing intervening) and both forms will be used. 385 // 386 // There are five choices for wildcard usage for improper dihedrals; 387 // 1) A - B - C - D (all four atoms, double specification allowed) 388 // 2) A - X - X - B 389 // 3) X - A - B - C 390 // 4) X - A - B - X 391 // 5) X - X - A - B 392 // When classifying an improper dihedral, the first acceptable match (from 393 // the above order) is chosen. The match may be made in either direction 394 // ( A - B - C - D = D - C - B - A). 395 // 396 // The periodicity value for dihedrals and improper dihedral terms 397 // must be an integer. If it is positive, then a cosine functional form is used. 398 // Only positive values of 1,2,3,4,5 and 6 are allowed for the vector, parallel 399 // vector and cray routines. Slow and scalar routines can use any positive 400 // integer and thus dihedral constrains can be of any periodicity. 401 // Reference angle 0.0 and 180.0 degree correspond to minimum in staggered 402 // and eclipsed respectively. Any reference angle is allowed. The value 403 // 180 should be prefered over -180 since it is parsed faster and more 404 // accuratly. When the periodicity is given as zero, for OTHER THAN THE 405 // FIRST dihdral in a multiple dihedral set, then a the amplitude is a 406 // constant added to the energy. This is needed to effect the 407 // Ryckaert-Bellemans potential for hydrocarbons (see below). 408 409 410 411 412 // First create look-up-table 413 map<string,vector<PAR::Improper>::const_iterator> improperLookUpTable; 414 for (vector<PAR::Improper>::const_iterator improper = par.impropers.begin(); 415 improper != par.impropers.end(); 416 ++improper){ 417 improperLookUpTable[improper->atom1+","+improper->atom2+","+improper->atom3+","+improper->atom4] = improper; 418 //report << (*improper)<<endr; 419 } 420 421 // Find the parameters from PAR 422 // loop over the improper list in the PSF object 423 for(vector<PSF::Improper>::const_iterator improper = psf.impropers.begin(); 424 improper != psf.impropers.end(); ++improper) { 425 426 // store the ID numbers of the atoms in this improper 427 int atom1 = improper->atom1 - 1; 428 int atom2 = improper->atom2 - 1; 429 int atom3 = improper->atom3 - 1; 430 int atom4 = improper->atom4 - 1; 431 432 // store the type names of the atoms in this improper 433 string improper1 = topo->atomTypes[topo->atoms[atom1].type].name; 434 string improper2 = topo->atomTypes[topo->atoms[atom2].type].name; 435 string improper3 = topo->atomTypes[topo->atoms[atom3].type].name; 436 string improper4 = topo->atomTypes[topo->atoms[atom4].type].name; 437 438 439 map<string,vector<PAR::Improper>::const_iterator>::const_iterator currentimproper = 440 improperLookUpTable.find(improper1+","+improper2+","+improper3+","+improper4); 441 if(currentimproper == improperLookUpTable.end()) 442 currentimproper = improperLookUpTable.find(improper4+","+improper3+","+improper2+","+improper1); 443 444 // Try wildcards if necessary 445 // 2) A - X - X - B 446 if(currentimproper == improperLookUpTable.end()){ 447 currentimproper = improperLookUpTable.find(improper1+",X,X,"+improper4); 448 if(currentimproper == improperLookUpTable.end()) 449 currentimproper = improperLookUpTable.find(improper4+",X,X,"+improper1); 450 451 } 452 // 3) X - A - B - C 453 if(currentimproper == improperLookUpTable.end()){ 454 currentimproper = improperLookUpTable.find("X,"+improper2+","+improper3+","+improper4); 455 if(currentimproper == improperLookUpTable.end()) 456 currentimproper = improperLookUpTable.find(improper4+","+improper3+","+improper2+",X"); 457 } 458 459 // 4) X - A - B - X 460 if(currentimproper == improperLookUpTable.end()){ 461 currentimproper = improperLookUpTable.find("X,"+improper2+","+improper3+",X"); 462 if(currentimproper == improperLookUpTable.end()) 463 currentimproper = improperLookUpTable.find("X,"+improper3+","+improper2+",X"); 464 } 465 466 // 5) X - X - A - B 467 if(currentimproper == improperLookUpTable.end()){ 468 currentimproper = improperLookUpTable.find("X,X,"+improper3+","+improper4); 469 if(currentimproper == improperLookUpTable.end()) 470 currentimproper = improperLookUpTable.find(improper4+","+improper3+",X,X"); 471 } 472 473 // if we still have not found this improper type in the PAR object, report an error 474 if(currentimproper == improperLookUpTable.end()) 475 report << error << "Could not find improper."<<endr; 476 477 // if we have found this improper type then copy the 478 // improper parameters into the topology 479 Torsion torsion; 480 torsion.atom1 = atom1; 481 torsion.atom2 = atom2; 482 torsion.atom3 = atom3; 483 torsion.atom4 = atom4; 484 torsion.periodicity.push_back(currentimproper->second->periodicity); 485 torsion.forceConstant.push_back(currentimproper->second->forceConstant); 486 torsion.phaseShift.push_back(dtor(currentimproper->second->phaseShift)); 487 torsion.multiplicity = 1; 488 topo->impropers.push_back(torsion); 489 // report << plain<< (wildcard ? "#":"") 490 // << improper1 <<"," 491 // << improper2 <<"," 492 // << improper3 <<"," 493 // << improper4 <<"," 494 // << torsion.atom1 <<"," 495 // << torsion.atom2 <<"," 496 // << torsion.atom3 <<"," 497 // << torsion.atom4 <<"," 498 // << torsion.periodicity[0] <<"," 499 // << torsion.periodicity.size() <<"," 500 // << torsion.forceConstant[0] <<"," 501 // << torsion.phaseShift[0] <<"," 502 // << torsion.multiplicity << endr; 503 } 504 505 506 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 507 // LennardJonesParameters 508 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 509 510 // get some array sizes 511 unsigned int sizeAtomTypes = topo->atomTypes.size(); 512 unsigned int sizeNonbondeds = par.nonbondeds.size(); 513 unsigned int sizeNbfixs = par.nbfixs.size(); 514 515 topo->lennardJonesParameters.resize(sizeAtomTypes); 516 517 // Nonbonded 518 for(unsigned int i=0;i<sizeAtomTypes;++i){ 519 int ti = 0; 520 unsigned int bi = sizeNonbondeds; 521 522 for(unsigned int k=0;k<sizeNonbondeds;++k){ 523 int ok = equalWildcard(par.nonbondeds[k].atom,topo->atomTypes[i].name); 524 if(ok > ti){ 525 bi = k; 526 ti = ok; 527 } 528 if(ti > 1) 529 break; 530 } 531 532 if(ti<=0) 533 report << error <<"Could not find matching parameter nonbonded of atom \'"<<topo->atomTypes[i].name<<"\'."<<endr; 534 535 for(unsigned int j=i;j<sizeAtomTypes;++j){ 536 int tj = 0; 537 unsigned int bj = sizeNonbondeds; 538 for(unsigned int k=0;k<sizeNonbondeds;++k){ 539 int ok = equalWildcard(par.nonbondeds[k].atom,topo->atomTypes[j].name); 540 if(ok > tj){ 541 bj = k; 542 tj = ok; 543 } 544 if(tj > 1) 545 break; 546 } 547 548 if(tj<=0) 549 report << error <<"Could not find matching parameter nonbonded of atom \'"<<topo->atomTypes[j].name<<"\'."<<endr; 550 551 LennardJonesParameters paramsij; 552 553 //Charmm28 554 Real sigma_i = par.nonbondeds[bi].sigma; 555 Real sigma_j = par.nonbondeds[bj].sigma; 556 Real sigma14_i = par.nonbondeds[bi].sigma14; 557 Real sigma14_j = par.nonbondeds[bj].sigma14; 558 559 Real epsilon_i = par.nonbondeds[bi].epsilon; 560 Real epsilon_j = par.nonbondeds[bj].epsilon; 561 Real epsilon14_i = par.nonbondeds[bi].epsilon14; 562 Real epsilon14_j = par.nonbondeds[bj].epsilon14; 563 564 Real r_ij = sigma_i + sigma_j; 565 Real e_ij = sqrt(epsilon_i * epsilon_j); 566 Real r14_ij = sigma14_i + sigma14_j; 567 Real e14_ij = sqrt(epsilon14_i * epsilon14_j); 568 569 paramsij.A = power<12>(r_ij) * e_ij; 570 paramsij.B = 2 * power<6>(r_ij) * e_ij; 571 paramsij.A14 = power<12>(r14_ij) * e14_ij; 572 paramsij.B14 = 2 * power<6>(r14_ij) * e14_ij; 573 //report << topo->atomTypes[i].name<<","<<bi<<","<<topo->atomTypes[j].name<<","<<bj<<" # "<<paramsij.A<<","<<paramsij.B<<","<<paramsij.A14<<","<<paramsij.B14<<endr; 574 topo->lennardJonesParameters.set(i, j, paramsij); 575 } 576 577 } 578 579 // NbFix 580 for(unsigned int k=0;k<sizeNbfixs;++k){ 581 582 //report << par.nbfixs[k].atom1 << ","<<par.nbfixs[k].atom2<<endr; 583 584 int ti = 0; 585 int tj = 0; 586 unsigned int bi = sizeNbfixs; 587 unsigned int bj = sizeNbfixs; 588 589 for(unsigned int i=0;i<sizeAtomTypes;++i){ 590 int ok = equalWildcard(par.nbfixs[k].atom1,topo->atomTypes[i].name); 591 if(ok > ti){ 592 bi = i; 593 ti = ok; 594 } 595 if(ti > 1) 596 break; 597 } 598 599 if(ti <=0) // ??? 600 continue; 601 602 for(unsigned int j=0;j<sizeAtomTypes;++j){ 603 int ok = equalWildcard(par.nbfixs[k].atom2,topo->atomTypes[j].name); 604 if(ok > tj){ 605 bj = j; 606 tj = ok; 607 } 608 if(tj > 1) 609 break; 610 } 611 612 if(tj<=0) 613 report << error <<"Could not find matching parameter nbfix of atoms \'"<<par.nbfixs[k].atom1<<"\' - '"<<par.nbfixs[k].atom2<<"\'."<<endr; 614 615 LennardJonesParameters paramsij; 616 617 paramsij.A = par.nbfixs[k].a; 618 paramsij.B = par.nbfixs[k].b; 619 paramsij.A14 = par.nbfixs[k].a14; 620 paramsij.B14 = par.nbfixs[k].b14; 621 topo->lennardJonesParameters.set(bi, bj, paramsij); 622 623 //report << topo->atomTypes[bi].name<<","<<topo->atomTypes[bj].name<<","<<paramsij.A<<","<<paramsij.B<<","<<paramsij.A14<<","<<paramsij.B14<<endr; 624 } // end loop over NbFix types 625 626 // store the molecule information 627 buildMoleculeTable(topo); 628 buildExclusionTable(topo,topo->exclude); 629 } 630 631 632 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 633 // 634 // buildMoleculeTable 635 // 636 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 637 buildMoleculeTable(GenericTopology * topo)638 void buildMoleculeTable(GenericTopology *topo){ 639 640 // *** First we clear all molecules *** 641 topo->molecules.clear(); 642 643 const unsigned int numAtoms = topo->atoms.size(); 644 645 // *** Collecting all possible connections, building the graph *** 646 vector<vector<int> > graph(numAtoms,vector<int>()); 647 set<pair<int,int> > pairs; 648 // *** Bonds *** 649 for(unsigned int i=0;i<topo->bonds.size();++i){ 650 int a1 = topo->bonds[i].atom1; 651 int a2 = topo->bonds[i].atom2; 652 graph[a1].push_back(a2); 653 graph[a2].push_back(a1); 654 pairs.insert(pair<int,int>(std::min(a1,a2),std::max(a1,a2))); 655 } 656 unsigned int count = pairs.size(); 657 658 // *** Angles *** 659 for(unsigned int i=0;i<topo->angles.size();++i){ 660 int a1 = topo->angles[i].atom1; 661 int a2 = topo->angles[i].atom2; 662 int a3 = topo->angles[i].atom3; 663 graph[a1].push_back(a2); 664 graph[a1].push_back(a3); 665 graph[a2].push_back(a1); 666 graph[a2].push_back(a3); 667 graph[a3].push_back(a1); 668 graph[a3].push_back(a2); 669 pairs.insert(pair<int,int>(std::min(a1,a2),std::max(a1,a2))); 670 pairs.insert(pair<int,int>(std::min(a3,a2),std::max(a3,a2))); 671 } 672 673 if(count < pairs.size()) 674 report << hint << "Angles added "<<pairs.size()-count<<" new bond(s)." <<endr; 675 count = pairs.size(); 676 677 // *** Dihedrals *** 678 for(unsigned int i=0;i<topo->dihedrals.size();++i){ 679 int a1 = topo->dihedrals[i].atom1; 680 int a2 = topo->dihedrals[i].atom2; 681 int a3 = topo->dihedrals[i].atom3; 682 int a4 = topo->dihedrals[i].atom4; 683 graph[a1].push_back(a2); 684 graph[a1].push_back(a3); 685 graph[a1].push_back(a4); 686 graph[a2].push_back(a1); 687 graph[a2].push_back(a3); 688 graph[a2].push_back(a4); 689 graph[a3].push_back(a1); 690 graph[a3].push_back(a2); 691 graph[a3].push_back(a4); 692 graph[a4].push_back(a1); 693 graph[a4].push_back(a2); 694 graph[a4].push_back(a3); 695 pairs.insert(pair<int,int>(std::min(a1,a2),std::max(a1,a2))); 696 pairs.insert(pair<int,int>(std::min(a3,a2),std::max(a3,a2))); 697 pairs.insert(pair<int,int>(std::min(a3,a4),std::max(a3,a4))); 698 } 699 if(count < pairs.size()) 700 report << hint << "Dihedrals added "<<pairs.size()-count<<" new bond(s)." <<endr; 701 count = pairs.size(); 702 703 // *** Impropers *** 704 set<pair<int,int> > pairsAddImpropers; 705 // Impropers are defined over the bonds 1-2,1-3,1-4 or 4-1,4-3,4-2 706 // but MTorsionSystemForce computes distances betweeen 1-2,2-3,3-4 707 // we have to take care about these differences ... 708 for(unsigned int i=0;i<topo->impropers.size();++i){ 709 int a1 = topo->impropers[i].atom1; 710 int a2 = topo->impropers[i].atom2; 711 int a3 = topo->impropers[i].atom3; 712 int a4 = topo->impropers[i].atom4; 713 graph[a1].push_back(a2); 714 graph[a1].push_back(a3); 715 graph[a1].push_back(a4); 716 graph[a2].push_back(a1); 717 graph[a2].push_back(a3); 718 graph[a2].push_back(a4); 719 graph[a3].push_back(a1); 720 graph[a3].push_back(a2); 721 graph[a3].push_back(a4); 722 graph[a4].push_back(a1); 723 graph[a4].push_back(a2); 724 graph[a4].push_back(a3); 725 pair<int,int> p0(std::min(a1,a2),std::max(a1,a2)); 726 pair<int,int> p1(std::min(a1,a3),std::max(a1,a3)); 727 pair<int,int> p2(std::min(a1,a4),std::max(a1,a4)); 728 pair<int,int> p3(std::min(a2,a3),std::max(a2,a3)); 729 pair<int,int> p4(std::min(a2,a4),std::max(a2,a4)); 730 pair<int,int> p5(std::min(a3,a4),std::max(a3,a4)); 731 int j0 = 0; 732 int j1 = 0; 733 int j2 = 0; 734 int j3 = 0; 735 int j4 = 0; 736 int j5 = 0; 737 if(pairs.find(p0) != pairs.end()) j0++; 738 if(pairs.find(p1) != pairs.end()) j1++; 739 if(pairs.find(p2) != pairs.end()) j2++; 740 if(pairs.find(p3) != pairs.end()) j3++; 741 if(pairs.find(p4) != pairs.end()) j4++; 742 if(pairs.find(p5) != pairs.end()) j5++; 743 if(j0+j1+j2+j3+j4+j5 < 3){ 744 pairs.insert(p0); 745 pairs.insert(p1); 746 pairs.insert(p2); 747 } 748 pairsAddImpropers.insert(p0); 749 pairsAddImpropers.insert(p3); 750 pairsAddImpropers.insert(p5); 751 } 752 if(count < pairs.size()) 753 report << hint << "Impropers added "<<pairs.size()-count<<" new bond(s)." <<endr; 754 755 // Now add the improper pairs 756 for(set<pair<int,int> >::const_iterator i= pairsAddImpropers.begin();i != pairsAddImpropers.end();++i) 757 pairs.insert(*i); 758 759 760 count = pairs.size(); 761 //report << hint << count << endr; 762 // To keep track which atoms already have been added 763 // to molecules. 764 vector<bool> unused(numAtoms,true); 765 766 767 // Recursively finding the atoms beloning to a molecule 768 for(unsigned int i=0;i<numAtoms;++i){ 769 vector<int> v; 770 vector<PairInt> p; 771 findNextNeighbor(i,v,p,unused,graph,pairs); 772 if(!v.empty()){ 773 std::sort(v.begin(),v.end()); 774 // add this atom list to the molecules array 775 Molecule mol; 776 mol.atoms = v; 777 for(unsigned int j=0;j<p.size();++j) 778 if(p[j].first > p[j].second) 779 swap(p[j].first,p[j].second); 780 std::sort(p.begin(),p.end()); 781 mol.pairs = p; 782 topo->molecules.push_back(mol); 783 } 784 } 785 786 // Uncomment to sort descending after size() 787 // std::sort(topo->molecules.begin(),topo->molecules.end(),cmpSize); 788 789 // Look up table for atoms 790 const string h("H"); 791 const string o("O"); 792 for(unsigned int i=0;i<topo->molecules.size();++i){ 793 Real mass = 0.0; 794 const vector<int>& mol = topo->molecules[i].atoms; 795 for(unsigned int j=0;j<mol.size();++j){ 796 int k = mol[j]; 797 topo->atoms[k].molecule = i; 798 mass += topo->atoms[k].scaledMass; 799 } 800 topo->molecules[i].mass = mass; 801 topo->molecules[i].water = (mol.size() == 3 && 802 ((topo->atomTypes[topo->atoms[mol[0]].type].symbolName == h && 803 topo->atomTypes[topo->atoms[mol[1]].type].symbolName == h && 804 topo->atomTypes[topo->atoms[mol[2]].type].symbolName == o) || 805 (topo->atomTypes[topo->atoms[mol[0]].type].symbolName == h && 806 topo->atomTypes[topo->atoms[mol[1]].type].symbolName == o && 807 topo->atomTypes[topo->atoms[mol[2]].type].symbolName == h) || 808 (topo->atomTypes[topo->atoms[mol[0]].type].symbolName == o && 809 topo->atomTypes[topo->atoms[mol[1]].type].symbolName == h && 810 topo->atomTypes[topo->atoms[mol[2]].type].symbolName == h))); 811 812 } 813 814 #if defined(DEBUG_PRINT_MOLECULETABLE) 815 report<< plain << endl << "[buildMoleculeTable]: molecule table printout:" << endl; 816 for(int i=0;i<topo->molecules.size();++i){ 817 for(int j=0;j<topo->molecules[i].size();++j){ 818 report << topo->molecules[i][j]<<" "; 819 } 820 report << endl; 821 } 822 report << endr; 823 #endif 824 } 825 826 827 //____________________________________________________________findNextNeighbor findNextNeighbor(int a,vector<int> & v,vector<PairInt> & p,vector<bool> & unused,const vector<vector<int>> & graph,set<PairInt> & pairs)828 void findNextNeighbor(int a, vector<int>& v, vector<PairInt>& p, vector<bool>& unused, 829 const vector<vector<int> >& graph, set<PairInt>& pairs){ 830 if(unused[a]){ 831 v.push_back(a); 832 unused[a] = false; 833 for(unsigned int i=0;i<graph[a].size();++i){ 834 set<PairInt>::iterator itr = pairs.find(PairInt(std::min(a,graph[a][i]),std::max(a,graph[a][i]))); 835 if(itr != pairs.end()){ 836 p.push_back(PairInt(a,graph[a][i])); 837 pairs.erase(itr); 838 } 839 findNextNeighbor(graph[a][i],v,p,unused,graph,pairs); 840 } 841 } 842 } 843 844 //_____________________________________________________________________ cmpSize 845 // bool cmpSize(const vector<int>& m1, const vector<int>& m2){ 846 // return (m1.size() > m2.size()); 847 // } 848 849 850 851 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 852 // 853 // buildExclusionTable 854 // 855 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 856 buildExclusionTable(GenericTopology * topo,const ExclusionType & exclusionType)857 void buildExclusionTable(GenericTopology* topo, const ExclusionType& exclusionType) { 858 859 if(!exclusionType.valid()) 860 report << error <<"[buildExclusionTable()] Exclusion type not defined/valid."<<endr; 861 862 topo->exclude = exclusionType; 863 864 const int numBonds = topo->bonds.size(); 865 const int numAtoms = topo->atoms.size(); 866 int a1,a2,a3; 867 868 // *** Reformatting bond list *** 869 vector< set<int> > bondsParsed(numAtoms); 870 for (int i = 0; i < numBonds; ++i) { 871 a1=topo->bonds[i].atom1; 872 a2=topo->bonds[i].atom2; 873 bondsParsed[a1].insert(a2); 874 bondsParsed[a2].insert(a1); 875 } 876 877 // *** Building exclusion list *** 878 topo->exclusions.resize(numAtoms); 879 topo->exclusions.clear(); 880 list<ExclusionPair> one3s; 881 int exclusionsInserted = 0; 882 if(exclusionType!=ExclusionType::NONE) { 883 for(a2 = 0; a2 < numAtoms; ++a2) { 884 for(set<int>::iterator i=bondsParsed[a2].begin();i!=bondsParsed[a2].end();++i) { 885 a1=*i; 886 if(a1<a2) { 887 topo->exclusions.add(a1,a2,EXCLUSION_FULL); 888 ++exclusionsInserted; 889 } 890 if(exclusionType==ExclusionType::ONE3||exclusionType==ExclusionType::ONE4||exclusionType==ExclusionType::ONE4MODIFIED) { 891 for(set<int>::iterator j=i; j!=bondsParsed[a2].end();++j) { 892 a3=*j; 893 topo->exclusions.add(a1,a3,EXCLUSION_FULL); 894 ++exclusionsInserted; 895 one3s.push_front(ExclusionPair(a1,a3)); 896 one3s.push_front(ExclusionPair(a3,a1)); 897 } 898 } 899 } 900 } 901 if(exclusionType==ExclusionType::ONE4||exclusionType==ExclusionType::ONE4MODIFIED) { 902 ExclusionClass currentType = EXCLUSION_NONE; 903 if(exclusionType==ExclusionType::ONE4) 904 currentType=EXCLUSION_FULL; 905 if(exclusionType==ExclusionType::ONE4MODIFIED) 906 currentType=EXCLUSION_MODIFIED; 907 int curSrc=0; 908 for(list<ExclusionPair>::iterator i=one3s.begin(); i!=one3s.end(); ++i,++curSrc) { 909 a1=i->a1; 910 a2=i->a2; 911 for(set<int>::iterator j=bondsParsed[a2].lower_bound(a1);j!=bondsParsed[a2].end();++j) { 912 a3=*j; 913 if(!topo->exclusions.check(a1,a3)) { 914 topo->exclusions.add(a1,a3,currentType); 915 ++exclusionsInserted; 916 } 917 } 918 } 919 } 920 } 921 topo->exclusions.cleanTemporaries(); 922 } 923 924 925 } 926