| 183 |
|
BendStamp* currentBend; |
| 184 |
|
TorsionStamp* currentTorsion; |
| 185 |
|
RigidBodyStamp* currentRigidBody; |
| 186 |
+ |
CutoffGroupStamp* currentCutoffGroup; |
| 187 |
+ |
CutoffGroup* myCutoffGroup; |
| 188 |
+ |
int nCutoffGroups;// number of cutoff group of a molecule defined in mdl file |
| 189 |
+ |
set<int> cutoffAtomSet; //atoms belong to cutoffgroup defined at mdl file |
| 190 |
|
|
| 191 |
|
bond_pair* theBonds; |
| 192 |
|
bend_set* theBends; |
| 220 |
|
molInfo.nBends = comp_stamps[stampID]->getNBends(); |
| 221 |
|
molInfo.nTorsions = comp_stamps[stampID]->getNTorsions(); |
| 222 |
|
molInfo.nRigidBodies = comp_stamps[stampID]->getNRigidBodies(); |
| 223 |
+ |
|
| 224 |
+ |
nCutoffGroups = comp_stamps[stampID]->getNCutoffGroups(); |
| 225 |
|
|
| 226 |
|
molInfo.myAtoms = &(info[k].atoms[atomOffset]); |
| 227 |
|
|
| 228 |
|
if (molInfo.nBonds > 0) |
| 229 |
< |
molInfo.myBonds = new (Bond *) [molInfo.nBonds]; |
| 229 |
> |
molInfo.myBonds = new Bond*[molInfo.nBonds]; |
| 230 |
|
else |
| 231 |
|
molInfo.myBonds = NULL; |
| 232 |
|
|
| 233 |
|
if (molInfo.nBends > 0) |
| 234 |
< |
molInfo.myBends = new (Bend *) [molInfo.nBends]; |
| 234 |
> |
molInfo.myBends = new Bend*[molInfo.nBends]; |
| 235 |
|
else |
| 236 |
|
molInfo.myBends = NULL; |
| 237 |
|
|
| 238 |
|
if (molInfo.nTorsions > 0) |
| 239 |
< |
molInfo.myTorsions = new (Torsion *) [molInfo.nTorsions]; |
| 239 |
> |
molInfo.myTorsions = new Torsion *[molInfo.nTorsions]; |
| 240 |
|
else |
| 241 |
|
molInfo.myTorsions = NULL; |
| 242 |
|
|
| 486 |
|
info[k].rigidBodies.push_back(myRB); |
| 487 |
|
} |
| 488 |
|
|
| 489 |
+ |
|
| 490 |
+ |
//create cutoff group for molecule |
| 491 |
+ |
|
| 492 |
+ |
cutoffAtomSet.clear(); |
| 493 |
+ |
molInfo.myCutoffGroups.clear(); |
| 494 |
+ |
|
| 495 |
+ |
for (j = 0; j < nCutoffGroups; j++){ |
| 496 |
+ |
|
| 497 |
+ |
currentCutoffGroup = comp_stamps[stampID]->getCutoffGroup(j); |
| 498 |
+ |
nMembers = currentCutoffGroup->getNMembers(); |
| 499 |
+ |
|
| 500 |
+ |
myCutoffGroup = new CutoffGroup(); |
| 501 |
+ |
|
| 502 |
+ |
for (int cg = 0; cg < nMembers; cg++) { |
| 503 |
+ |
|
| 504 |
+ |
// molI is atom numbering inside this molecule |
| 505 |
+ |
molI = currentCutoffGroup->getMember(cg); |
| 506 |
+ |
|
| 507 |
+ |
// tempI is atom numbering on local processor |
| 508 |
+ |
tempI = molI + atomOffset; |
| 509 |
+ |
|
| 510 |
+ |
myCutoffGroup->addAtom(info[k].atoms[tempI]); |
| 511 |
+ |
|
| 512 |
+ |
cutoffAtomSet.insert(tempI); |
| 513 |
+ |
} |
| 514 |
+ |
|
| 515 |
+ |
molInfo.myCutoffGroups.push_back(myCutoffGroup); |
| 516 |
+ |
}//end for (j = 0; j < molInfo.nCutoffGroups; j++) |
| 517 |
+ |
|
| 518 |
+ |
//creat a cutoff group for every atom in current molecule which does not belong to cutoffgroup defined at mdl file |
| 519 |
|
|
| 520 |
+ |
for(j = 0; j < molInfo.nAtoms; j++){ |
| 521 |
+ |
|
| 522 |
+ |
if(cutoffAtomSet.find(molInfo.myAtoms[j]->getIndex()) == cutoffAtomSet.end()){ |
| 523 |
+ |
myCutoffGroup = new CutoffGroup(); |
| 524 |
+ |
myCutoffGroup->addAtom(molInfo.myAtoms[j]); |
| 525 |
+ |
molInfo.myCutoffGroups.push_back(myCutoffGroup); |
| 526 |
+ |
} |
| 527 |
+ |
|
| 528 |
+ |
} |
| 529 |
+ |
|
| 530 |
+ |
|
| 531 |
+ |
|
| 532 |
+ |
|
| 533 |
|
// After this is all set up, scan through the atoms to |
| 534 |
|
// see if they can be added to the integrableObjects: |
| 535 |
|
|
| 585 |
|
MPIcheckPoint(); |
| 586 |
|
#endif // is_mpi |
| 587 |
|
|
| 539 |
– |
// clean up the forcefield |
| 540 |
– |
|
| 541 |
– |
if (!globals->haveRcut()){ |
| 542 |
– |
|
| 543 |
– |
the_ff->calcRcut(); |
| 544 |
– |
|
| 545 |
– |
} else { |
| 546 |
– |
|
| 547 |
– |
the_ff->setRcut( globals->getRcut() ); |
| 548 |
– |
} |
| 549 |
– |
|
| 550 |
– |
the_ff->cleanMe(); |
| 588 |
|
} |
| 589 |
|
|
| 590 |
|
void SimSetup::initFromBass(void){ |
| 948 |
|
|
| 949 |
|
info[i].useInitXSstate = globals->getUseInitXSstate(); |
| 950 |
|
info[i].orthoTolerance = globals->getOrthoBoxTolerance(); |
| 951 |
< |
|
| 951 |
> |
|
| 952 |
> |
// check for thermodynamic integration |
| 953 |
> |
if (globals->getUseThermInt()) { |
| 954 |
> |
if (globals->haveThermIntLambda() && globals->haveThermIntK()) { |
| 955 |
> |
info[i].useThermInt = globals->getUseThermInt(); |
| 956 |
> |
info[i].thermIntLambda = globals->getThermIntLambda(); |
| 957 |
> |
info[i].thermIntK = globals->getThermIntK(); |
| 958 |
> |
|
| 959 |
> |
Restraints *myRestraint = new Restraints(tot_nmol, info[i].thermIntLambda, info[i].thermIntK); |
| 960 |
> |
info[i].restraint = myRestraint; |
| 961 |
> |
} |
| 962 |
> |
else { |
| 963 |
> |
sprintf(painCave.errMsg, |
| 964 |
> |
"SimSetup Error:\n" |
| 965 |
> |
"\tKeyword useThermInt was set to 'true' but\n" |
| 966 |
> |
"\tthermodynamicIntegrationLambda (and/or\n" |
| 967 |
> |
"\tthermodynamicIntegrationK) was not specified.\n" |
| 968 |
> |
"\tPlease provide a lambda value and k value in your .bass file.\n"); |
| 969 |
> |
painCave.isFatal = 1; |
| 970 |
> |
simError(); |
| 971 |
> |
} |
| 972 |
> |
} |
| 973 |
> |
else if(globals->haveThermIntLambda() || globals->haveThermIntK()){ |
| 974 |
> |
sprintf(painCave.errMsg, |
| 975 |
> |
"SimSetup Warning: If you want to use Thermodynamic\n" |
| 976 |
> |
"\tIntegration, set useThermInt to 'true' in your .bass file.\n" |
| 977 |
> |
"\tThe useThermInt keyword is 'false' by default, so your\n" |
| 978 |
> |
"\tlambda and/or k values are being ignored.\n"); |
| 979 |
> |
painCave.isFatal = 0; |
| 980 |
> |
simError(); |
| 981 |
> |
} |
| 982 |
|
} |
| 983 |
|
|
| 984 |
|
//setup seed for random number generator |
| 1056 |
|
|
| 1057 |
|
double theRcut, theRsw; |
| 1058 |
|
|
| 1059 |
+ |
if (globals->haveRcut()) { |
| 1060 |
+ |
theRcut = globals->getRcut(); |
| 1061 |
+ |
|
| 1062 |
+ |
if (globals->haveRsw()) |
| 1063 |
+ |
theRsw = globals->getRsw(); |
| 1064 |
+ |
else |
| 1065 |
+ |
theRsw = theRcut; |
| 1066 |
+ |
|
| 1067 |
+ |
info[i].setDefaultRcut(theRcut, theRsw); |
| 1068 |
+ |
|
| 1069 |
+ |
} else { |
| 1070 |
+ |
|
| 1071 |
+ |
the_ff->calcRcut(); |
| 1072 |
+ |
theRcut = info[i].getRcut(); |
| 1073 |
+ |
|
| 1074 |
+ |
if (globals->haveRsw()) |
| 1075 |
+ |
theRsw = globals->getRsw(); |
| 1076 |
+ |
else |
| 1077 |
+ |
theRsw = theRcut; |
| 1078 |
+ |
|
| 1079 |
+ |
info[i].setDefaultRcut(theRcut, theRsw); |
| 1080 |
+ |
} |
| 1081 |
+ |
|
| 1082 |
|
if (globals->getUseRF()){ |
| 1083 |
|
info[i].useReactionField = 1; |
| 1084 |
< |
|
| 1084 |
> |
|
| 1085 |
|
if (!globals->haveRcut()){ |
| 1086 |
|
sprintf(painCave.errMsg, |
| 1087 |
|
"SimSetup Warning: No value was set for the cutoffRadius.\n" |
| 1158 |
|
strcpy(checkPointMsg, "post processing checks out"); |
| 1159 |
|
MPIcheckPoint(); |
| 1160 |
|
#endif // is_mpi |
| 1161 |
+ |
|
| 1162 |
+ |
// clean up the forcefield |
| 1163 |
+ |
the_ff->cleanMe(); |
| 1164 |
|
} |
| 1165 |
|
|
| 1166 |
|
void SimSetup::initSystemCoords(void){ |
| 1291 |
|
} |
| 1292 |
|
} |
| 1293 |
|
|
| 1294 |
+ |
strcpy(info[k].rawPotName, inFileName); |
| 1295 |
+ |
nameLength = strlen(info[k].rawPotName); |
| 1296 |
+ |
endTest = &(info[k].rawPotName[nameLength - 5]); |
| 1297 |
+ |
if (!strcmp(endTest, ".bass")){ |
| 1298 |
+ |
strcpy(endTest, ".raw"); |
| 1299 |
+ |
} |
| 1300 |
+ |
else if (!strcmp(endTest, ".BASS")){ |
| 1301 |
+ |
strcpy(endTest, ".raw"); |
| 1302 |
+ |
} |
| 1303 |
+ |
else{ |
| 1304 |
+ |
endTest = &(info[k].rawPotName[nameLength - 4]); |
| 1305 |
+ |
if (!strcmp(endTest, ".bss")){ |
| 1306 |
+ |
strcpy(endTest, ".raw"); |
| 1307 |
+ |
} |
| 1308 |
+ |
else if (!strcmp(endTest, ".mdl")){ |
| 1309 |
+ |
strcpy(endTest, ".raw"); |
| 1310 |
+ |
} |
| 1311 |
+ |
else{ |
| 1312 |
+ |
strcat(info[k].rawPotName, ".raw"); |
| 1313 |
+ |
} |
| 1314 |
+ |
} |
| 1315 |
+ |
|
| 1316 |
|
#ifdef IS_MPI |
| 1317 |
|
|
| 1318 |
|
} |
| 1397 |
|
LinkedMolStamp* headStamp = new LinkedMolStamp(); |
| 1398 |
|
LinkedMolStamp* currentStamp = NULL; |
| 1399 |
|
comp_stamps = new MoleculeStamp * [n_components]; |
| 1400 |
+ |
bool haveCutoffGroups; |
| 1401 |
|
|
| 1402 |
+ |
haveCutoffGroups = false; |
| 1403 |
+ |
|
| 1404 |
|
// make an array of molecule stamps that match the components used. |
| 1405 |
|
// also extract the used stamps out into a separate linked list |
| 1406 |
|
|
| 1435 |
|
headStamp->add(currentStamp); |
| 1436 |
|
comp_stamps[i] = headStamp->match(id); |
| 1437 |
|
} |
| 1438 |
+ |
|
| 1439 |
+ |
if(comp_stamps[i]->getNCutoffGroups() > 0) |
| 1440 |
+ |
haveCutoffGroups = true; |
| 1441 |
|
} |
| 1442 |
+ |
|
| 1443 |
+ |
for (i = 0; i < nInfo; i++) |
| 1444 |
+ |
info[i].haveCutoffGroups = haveCutoffGroups; |
| 1445 |
|
|
| 1446 |
|
#ifdef IS_MPI |
| 1447 |
|
strcpy(checkPointMsg, "Component stamps successfully extracted\n"); |
| 1616 |
|
info[l].atoms = the_atoms; |
| 1617 |
|
info[l].molecules = the_molecules; |
| 1618 |
|
info[l].nGlobalExcludes = 0; |
| 1619 |
< |
|
| 1619 |
> |
|
| 1620 |
|
the_ff->setSimInfo(info); |
| 1621 |
|
} |
| 1622 |
|
} |
