# | Line 10 | Line 10 | |
---|---|---|
10 | #include "Integrator.hpp" | |
11 | #include "simError.h" | |
12 | #include "RigidBody.hpp" | |
13 | – | //#include "ConjugateMinimizer.hpp" |
13 | #include "OOPSEMinimizer.hpp" | |
14 | + | #include "ConstraintElement.hpp" |
15 | + | #include "ConstraintPair.hpp" |
16 | + | #include "ConstraintManager.hpp" |
17 | ||
18 | #ifdef IS_MPI | |
19 | #include "mpiBASS.h" | |
# | Line 156 | Line 158 | void SimSetup::createSim(void){ | |
158 | ||
159 | initFortran(); | |
160 | ||
161 | + | //creat constraint manager |
162 | + | for(int i = 0; i < nInfo; i++) |
163 | + | info[i].consMan = new ConstraintManager(&info[i]); |
164 | + | |
165 | if (globals->haveMinimizer()) | |
166 | // make minimizer | |
167 | makeMinimizer(); | |
# | Line 170 | Line 176 | void SimSetup::makeMolecules(void){ | |
176 | int i, j, k; | |
177 | int exI, exJ, exK, exL, slI, slJ; | |
178 | int tempI, tempJ, tempK, tempL; | |
179 | < | int molI; |
180 | < | int stampID, atomOffset, rbOffset; |
179 | > | int molI, globalID; |
180 | > | int stampID, atomOffset, rbOffset, groupOffset; |
181 | molInit molInfo; | |
182 | DirectionalAtom* dAtom; | |
183 | RigidBody* myRB; | |
# | Line 183 | Line 189 | void SimSetup::makeMolecules(void){ | |
189 | BendStamp* currentBend; | |
190 | TorsionStamp* currentTorsion; | |
191 | RigidBodyStamp* currentRigidBody; | |
192 | + | CutoffGroupStamp* currentCutoffGroup; |
193 | + | CutoffGroup* myCutoffGroup; |
194 | + | int nCutoffGroups;// number of cutoff group of a molecule defined in mdl file |
195 | + | set<int> cutoffAtomSet; //atoms belong to cutoffgroup defined at mdl file |
196 | ||
197 | bond_pair* theBonds; | |
198 | bend_set* theBends; | |
# | Line 194 | Line 204 | void SimSetup::makeMolecules(void){ | |
204 | char* molName; | |
205 | char rbName[100]; | |
206 | ||
207 | + | ConstraintPair* consPair; //constraint pair |
208 | + | ConstraintElement* consElement1; //first element of constraint pair |
209 | + | ConstraintElement* consElement2; //second element of constraint pair |
210 | + | int whichRigidBody; |
211 | + | int consAtomIndex; //index of constraint atom in rigid body's atom array |
212 | + | vector<pair<int, int> > jointAtoms; |
213 | + | double bondLength2; |
214 | //init the forceField paramters | |
215 | ||
216 | the_ff->readParams(); | |
# | Line 204 | Line 221 | void SimSetup::makeMolecules(void){ | |
221 | ||
222 | for (k = 0; k < nInfo; k++){ | |
223 | the_ff->setSimInfo(&(info[k])); | |
224 | + | |
225 | + | #ifdef IS_MPI |
226 | + | info[k].globalGroupMembership = new int[mpiSim->getNAtomsGlobal()]; |
227 | + | for (i = 0; i < mpiSim->getNAtomsGlobal(); i++) |
228 | + | info[k].globalGroupMembership[i] = 0; |
229 | + | #else |
230 | + | info[k].globalGroupMembership = new int[info[k].n_atoms]; |
231 | + | for (i = 0; i < info[k].n_atoms; i++) |
232 | + | info[k].globalGroupMembership[i] = 0; |
233 | + | #endif |
234 | ||
235 | atomOffset = 0; | |
236 | + | groupOffset = 0; |
237 | ||
238 | for (i = 0; i < info[k].n_mol; i++){ | |
239 | stampID = info[k].molecules[i].getStampID(); | |
# | Line 216 | Line 244 | void SimSetup::makeMolecules(void){ | |
244 | molInfo.nBends = comp_stamps[stampID]->getNBends(); | |
245 | molInfo.nTorsions = comp_stamps[stampID]->getNTorsions(); | |
246 | molInfo.nRigidBodies = comp_stamps[stampID]->getNRigidBodies(); | |
247 | + | |
248 | + | nCutoffGroups = comp_stamps[stampID]->getNCutoffGroups(); |
249 | ||
250 | molInfo.myAtoms = &(info[k].atoms[atomOffset]); | |
251 | ||
252 | if (molInfo.nBonds > 0) | |
253 | < | molInfo.myBonds = new (Bond *) [molInfo.nBonds]; |
253 | > | molInfo.myBonds = new Bond*[molInfo.nBonds]; |
254 | else | |
255 | molInfo.myBonds = NULL; | |
256 | ||
257 | if (molInfo.nBends > 0) | |
258 | < | molInfo.myBends = new (Bend *) [molInfo.nBends]; |
258 | > | molInfo.myBends = new Bend*[molInfo.nBends]; |
259 | else | |
260 | molInfo.myBends = NULL; | |
261 | ||
262 | if (molInfo.nTorsions > 0) | |
263 | < | molInfo.myTorsions = new (Torsion *) [molInfo.nTorsions]; |
263 | > | molInfo.myTorsions = new Torsion *[molInfo.nTorsions]; |
264 | else | |
265 | molInfo.myTorsions = NULL; | |
266 | ||
# | Line 268 | Line 298 | void SimSetup::makeMolecules(void){ | |
298 | ||
299 | molInfo.myAtoms[j]->setType(currentAtom->getType()); | |
300 | #ifdef IS_MPI | |
271 | – | |
301 | molInfo.myAtoms[j]->setGlobalIndex(globalAtomIndex[j + atomOffset]); | |
273 | – | |
302 | #endif // is_mpi | |
303 | } | |
304 | ||
# | Line 479 | Line 507 | void SimSetup::makeMolecules(void){ | |
507 | molInfo.myRigidBodies.push_back(myRB); | |
508 | info[k].rigidBodies.push_back(myRB); | |
509 | } | |
510 | + | |
511 | + | |
512 | + | //create cutoff group for molecule |
513 | + | |
514 | + | cutoffAtomSet.clear(); |
515 | + | molInfo.myCutoffGroups.clear(); |
516 | + | |
517 | + | for (j = 0; j < nCutoffGroups; j++){ |
518 | + | |
519 | + | currentCutoffGroup = comp_stamps[stampID]->getCutoffGroup(j); |
520 | + | nMembers = currentCutoffGroup->getNMembers(); |
521 | + | |
522 | + | myCutoffGroup = new CutoffGroup(); |
523 | + | |
524 | + | #ifdef IS_MPI |
525 | + | myCutoffGroup->setGlobalIndex(globalGroupIndex[groupOffset]); |
526 | + | #else |
527 | + | myCutoffGroup->setGlobalIndex(groupOffset); |
528 | + | #endif |
529 | + | |
530 | + | for (int cg = 0; cg < nMembers; cg++) { |
531 | + | |
532 | + | // molI is atom numbering inside this molecule |
533 | + | molI = currentCutoffGroup->getMember(cg); |
534 | + | |
535 | + | // tempI is atom numbering on local processor |
536 | + | tempI = molI + atomOffset; |
537 | + | |
538 | + | #ifdef IS_MPI |
539 | + | globalID = info[k].atoms[tempI]->getGlobalIndex(); |
540 | + | info[k].globalGroupMembership[globalID] = globalGroupIndex[groupOffset]; |
541 | + | #else |
542 | + | globalID = info[k].atoms[tempI]->getIndex(); |
543 | + | info[k].globalGroupMembership[globalID] = groupOffset; |
544 | + | #endif |
545 | + | myCutoffGroup->addAtom(info[k].atoms[tempI]); |
546 | + | cutoffAtomSet.insert(tempI); |
547 | + | } |
548 | + | |
549 | + | molInfo.myCutoffGroups.push_back(myCutoffGroup); |
550 | + | groupOffset++; |
551 | + | |
552 | + | }//end for (j = 0; j < molInfo.nCutoffGroups; j++) |
553 | ||
554 | + | |
555 | + | // create a cutoff group for every atom in current molecule which |
556 | + | // does not belong to cutoffgroup defined at mdl file |
557 | + | |
558 | + | for(j = 0; j < molInfo.nAtoms; j++){ |
559 | + | |
560 | + | if(cutoffAtomSet.find(molInfo.myAtoms[j]->getIndex()) == cutoffAtomSet.end()){ |
561 | + | myCutoffGroup = new CutoffGroup(); |
562 | + | myCutoffGroup->addAtom(molInfo.myAtoms[j]); |
563 | + | |
564 | + | #ifdef IS_MPI |
565 | + | myCutoffGroup->setGlobalIndex(globalGroupIndex[groupOffset]); |
566 | + | globalID = info[k].atoms[atomOffset + j]->getGlobalIndex(); |
567 | + | info[k].globalGroupMembership[globalID] = globalGroupIndex[groupOffset]; |
568 | + | #else |
569 | + | myCutoffGroup->setGlobalIndex(groupOffset); |
570 | + | globalID = info[k].atoms[atomOffset + j]->getIndex(); |
571 | + | info[k].globalGroupMembership[globalID] = groupOffset; |
572 | + | #endif |
573 | + | molInfo.myCutoffGroups.push_back(myCutoffGroup); |
574 | + | groupOffset++; |
575 | + | } |
576 | + | } |
577 | ||
578 | // After this is all set up, scan through the atoms to | |
579 | // see if they can be added to the integrableObjects: | |
# | Line 511 | Line 605 | void SimSetup::makeMolecules(void){ | |
605 | info[k].integrableObjects.push_back(mySD); | |
606 | molInfo.myIntegrableObjects.push_back(mySD); | |
607 | } | |
608 | < | |
515 | < | |
608 | > | |
609 | // send the arrays off to the forceField for init. | |
610 | ||
611 | the_ff->initializeAtoms(molInfo.nAtoms, molInfo.myAtoms); | |
# | Line 521 | Line 614 | void SimSetup::makeMolecules(void){ | |
614 | the_ff->initializeTorsions(molInfo.nTorsions, molInfo.myTorsions, | |
615 | theTorsions); | |
616 | ||
524 | – | info[k].molecules[i].initialize(molInfo); |
617 | ||
618 | + | //creat ConstraintPair. |
619 | + | molInfo.myConstraintPairs.clear(); |
620 | + | |
621 | + | for (j = 0; j < molInfo.nBonds; j++){ |
622 | ||
623 | + | //if bond is constrained bond, add it into constraint pair |
624 | + | if(molInfo.myBonds[j]->is_constrained()){ |
625 | + | |
626 | + | //if both atoms are in the same rigid body, just skip it |
627 | + | currentBond = comp_stamps[stampID]->getBond(j); |
628 | + | |
629 | + | if(!comp_stamps[stampID]->isBondInSameRigidBody(currentBond)){ |
630 | + | |
631 | + | tempI = currentBond->getA() + atomOffset; |
632 | + | if( comp_stamps[stampID]->isAtomInRigidBody(currentBond->getA(), whichRigidBody, consAtomIndex)) |
633 | + | consElement1 = new ConstraintRigidBody(molInfo.myRigidBodies[whichRigidBody], consAtomIndex); |
634 | + | else |
635 | + | consElement1 = new ConstraintAtom(info[k].atoms[tempI]); |
636 | + | |
637 | + | tempJ = currentBond->getB() + atomOffset; |
638 | + | if(comp_stamps[stampID]->isAtomInRigidBody(currentBond->getB(), whichRigidBody, consAtomIndex)) |
639 | + | consElement2 = new ConstraintRigidBody(molInfo.myRigidBodies[whichRigidBody], consAtomIndex); |
640 | + | else |
641 | + | consElement2 = new ConstraintAtom(info[k].atoms[tempJ]); |
642 | + | |
643 | + | bondLength2 = molInfo.myBonds[j]->get_constraint()->get_dsqr(); |
644 | + | consPair = new DistanceConstraintPair(consElement1, consElement2, bondLength2); |
645 | + | |
646 | + | molInfo.myConstraintPairs.push_back(consPair); |
647 | + | } |
648 | + | }//end if(molInfo.myBonds[j]->is_constrained()) |
649 | + | } |
650 | + | |
651 | + | //loop over rigid bodies, if two rigid bodies share same joint, creat a JointConstraintPair |
652 | + | for (int rb1 = 0; rb1 < molInfo.nRigidBodies -1 ; rb1++){ |
653 | + | for (int rb2 = rb1 + 1; rb2 < molInfo.nRigidBodies ; rb2++){ |
654 | + | |
655 | + | jointAtoms = comp_stamps[stampID]->getJointAtoms(rb1, rb2); |
656 | + | |
657 | + | for(size_t m = 0; m < jointAtoms.size(); m++){ |
658 | + | consElement1 = new ConstraintRigidBody(molInfo.myRigidBodies[rb1], jointAtoms[m].first); |
659 | + | consElement2 = new ConstraintRigidBody(molInfo.myRigidBodies[rb2], jointAtoms[m].second); |
660 | + | |
661 | + | consPair = new JointConstraintPair(consElement1, consElement2); |
662 | + | molInfo.myConstraintPairs.push_back(consPair); |
663 | + | } |
664 | + | |
665 | + | } |
666 | + | } |
667 | + | |
668 | + | |
669 | + | info[k].molecules[i].initialize(molInfo); |
670 | + | |
671 | + | |
672 | atomOffset += molInfo.nAtoms; | |
673 | delete[] theBonds; | |
674 | delete[] theBends; | |
675 | delete[] theTorsions; | |
676 | < | } |
532 | < | } |
676 | > | } |
677 | ||
534 | – | #ifdef IS_MPI |
535 | – | sprintf(checkPointMsg, "all molecules initialized succesfully"); |
536 | – | MPIcheckPoint(); |
537 | – | #endif // is_mpi |
678 | ||
539 | – | // clean up the forcefield |
679 | ||
680 | < | if (!globals->haveLJrcut()){ |
680 | > | #ifdef IS_MPI |
681 | > | // Since the globalGroupMembership has been zero filled and we've only |
682 | > | // poked values into the atoms we know, we can do an Allreduce |
683 | > | // to get the full globalGroupMembership array (We think). |
684 | > | // This would be prettier if we could use MPI_IN_PLACE like the MPI-2 |
685 | > | // docs said we could. |
686 | ||
687 | < | the_ff->calcRcut(); |
687 | > | int* ggMjunk = new int[mpiSim->getNAtomsGlobal()]; |
688 | ||
689 | < | } else { |
690 | < | |
691 | < | the_ff->setRcut( globals->getLJrcut() ); |
689 | > | MPI_Allreduce(info[k].globalGroupMembership, |
690 | > | ggMjunk, |
691 | > | mpiSim->getNAtomsGlobal(), |
692 | > | MPI_INT, MPI_SUM, MPI_COMM_WORLD); |
693 | > | |
694 | > | for (i = 0; i < mpiSim->getNAtomsGlobal(); i++) |
695 | > | info[k].globalGroupMembership[i] = ggMjunk[i]; |
696 | > | |
697 | > | delete[] ggMjunk; |
698 | > | |
699 | > | #endif |
700 | > | |
701 | > | |
702 | > | |
703 | } | |
704 | ||
705 | < | the_ff->cleanMe(); |
705 | > | #ifdef IS_MPI |
706 | > | sprintf(checkPointMsg, "all molecules initialized succesfully"); |
707 | > | MPIcheckPoint(); |
708 | > | #endif // is_mpi |
709 | > | |
710 | } | |
711 | ||
712 | void SimSetup::initFromBass(void){ | |
# | Line 757 | Line 916 | void SimSetup::gatherInfo(void){ | |
916 | painCave.isFatal = 1; | |
917 | simError(); | |
918 | } | |
919 | + | if (globals->haveForceFieldVariant()) { |
920 | + | strcpy(forcefield_variant, globals->getForceFieldVariant()); |
921 | + | has_forcefield_variant = 1; |
922 | + | } |
923 | + | |
924 | + | // get the ensemble |
925 | ||
761 | – | // get the ensemble |
762 | – | |
926 | strcpy(ensemble, globals->getEnsemble()); | |
927 | ||
928 | if (!strcasecmp(ensemble, "NVE")){ | |
# | Line 911 | Line 1074 | void SimSetup::gatherInfo(void){ | |
1074 | ||
1075 | info[i].useInitXSstate = globals->getUseInitXSstate(); | |
1076 | info[i].orthoTolerance = globals->getOrthoBoxTolerance(); | |
1077 | < | |
1077 | > | |
1078 | > | // check for thermodynamic integration |
1079 | > | if (globals->getUseSolidThermInt() && !globals->getUseLiquidThermInt()) { |
1080 | > | if (globals->haveThermIntLambda() && globals->haveThermIntK()) { |
1081 | > | info[i].useSolidThermInt = globals->getUseSolidThermInt(); |
1082 | > | info[i].thermIntLambda = globals->getThermIntLambda(); |
1083 | > | info[i].thermIntK = globals->getThermIntK(); |
1084 | > | |
1085 | > | Restraints *myRestraint = new Restraints(tot_nmol, info[i].thermIntLambda, info[i].thermIntK); |
1086 | > | info[i].restraint = myRestraint; |
1087 | > | } |
1088 | > | else { |
1089 | > | sprintf(painCave.errMsg, |
1090 | > | "SimSetup Error:\n" |
1091 | > | "\tKeyword useSolidThermInt was set to 'true' but\n" |
1092 | > | "\tthermodynamicIntegrationLambda (and/or\n" |
1093 | > | "\tthermodynamicIntegrationK) was not specified.\n" |
1094 | > | "\tPlease provide a lambda value and k value in your .bass file.\n"); |
1095 | > | painCave.isFatal = 1; |
1096 | > | simError(); |
1097 | > | } |
1098 | > | } |
1099 | > | else if(globals->getUseLiquidThermInt()) { |
1100 | > | if (globals->getUseSolidThermInt()) { |
1101 | > | sprintf( painCave.errMsg, |
1102 | > | "SimSetup Warning: It appears that you have both solid and\n" |
1103 | > | "\tliquid thermodynamic integration activated in your .bass\n" |
1104 | > | "\tfile. To avoid confusion, specify only one technique in\n" |
1105 | > | "\tyour .bass file. Liquid-state thermodynamic integration\n" |
1106 | > | "\twill be assumed for the current simulation. If this is not\n" |
1107 | > | "\twhat you desire, set useSolidThermInt to 'true' and\n" |
1108 | > | "\tuseLiquidThermInt to 'false' in your .bass file.\n"); |
1109 | > | painCave.isFatal = 0; |
1110 | > | simError(); |
1111 | > | } |
1112 | > | if (globals->haveThermIntLambda() && globals->haveThermIntK()) { |
1113 | > | info[i].useLiquidThermInt = globals->getUseLiquidThermInt(); |
1114 | > | info[i].thermIntLambda = globals->getThermIntLambda(); |
1115 | > | info[i].thermIntK = globals->getThermIntK(); |
1116 | > | } |
1117 | > | else { |
1118 | > | sprintf(painCave.errMsg, |
1119 | > | "SimSetup Error:\n" |
1120 | > | "\tKeyword useLiquidThermInt was set to 'true' but\n" |
1121 | > | "\tthermodynamicIntegrationLambda (and/or\n" |
1122 | > | "\tthermodynamicIntegrationK) was not specified.\n" |
1123 | > | "\tPlease provide a lambda value and k value in your .bass file.\n"); |
1124 | > | painCave.isFatal = 1; |
1125 | > | simError(); |
1126 | > | } |
1127 | > | } |
1128 | > | else if(globals->haveThermIntLambda() || globals->haveThermIntK()){ |
1129 | > | sprintf(painCave.errMsg, |
1130 | > | "SimSetup Warning: If you want to use Thermodynamic\n" |
1131 | > | "\tIntegration, set useSolidThermInt or useLiquidThermInt to\n" |
1132 | > | "\t'true' in your .bass file. These keywords are set to\n" |
1133 | > | "\t'false' by default, so your lambda and/or k values are\n" |
1134 | > | "\tbeing ignored.\n"); |
1135 | > | painCave.isFatal = 0; |
1136 | > | simError(); |
1137 | > | } |
1138 | } | |
1139 | ||
1140 | //setup seed for random number generator | |
# | Line 987 | Line 1210 | void SimSetup::finalInfoCheck(void){ | |
1210 | MPI_Allreduce(&myUse, &usesDipoles, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | |
1211 | #endif //is_mpi | |
1212 | ||
1213 | < | double theEcr, theEst; |
1213 | > | double theRcut, theRsw; |
1214 | ||
1215 | + | if (globals->haveRcut()) { |
1216 | + | theRcut = globals->getRcut(); |
1217 | + | |
1218 | + | if (globals->haveRsw()) |
1219 | + | theRsw = globals->getRsw(); |
1220 | + | else |
1221 | + | theRsw = theRcut; |
1222 | + | |
1223 | + | info[i].setDefaultRcut(theRcut, theRsw); |
1224 | + | |
1225 | + | } else { |
1226 | + | |
1227 | + | the_ff->calcRcut(); |
1228 | + | theRcut = info[i].getRcut(); |
1229 | + | |
1230 | + | if (globals->haveRsw()) |
1231 | + | theRsw = globals->getRsw(); |
1232 | + | else |
1233 | + | theRsw = theRcut; |
1234 | + | |
1235 | + | info[i].setDefaultRcut(theRcut, theRsw); |
1236 | + | } |
1237 | + | |
1238 | if (globals->getUseRF()){ | |
1239 | info[i].useReactionField = 1; | |
1240 | < | |
1241 | < | if (!globals->haveECR()){ |
1240 | > | |
1241 | > | if (!globals->haveRcut()){ |
1242 | sprintf(painCave.errMsg, | |
1243 | < | "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n" |
1243 | > | "SimSetup Warning: No value was set for the cutoffRadius.\n" |
1244 | "\tOOPSE will use a default value of 15.0 angstroms" | |
1245 | < | "\tfor the electrostaticCutoffRadius.\n"); |
1245 | > | "\tfor the cutoffRadius.\n"); |
1246 | painCave.isFatal = 0; | |
1247 | simError(); | |
1248 | < | theEcr = 15.0; |
1248 | > | theRcut = 15.0; |
1249 | } | |
1250 | else{ | |
1251 | < | theEcr = globals->getECR(); |
1251 | > | theRcut = globals->getRcut(); |
1252 | } | |
1253 | ||
1254 | < | if (!globals->haveEST()){ |
1254 | > | if (!globals->haveRsw()){ |
1255 | sprintf(painCave.errMsg, | |
1256 | < | "SimSetup Warning: No value was set for electrostaticSkinThickness.\n" |
1256 | > | "SimSetup Warning: No value was set for switchingRadius.\n" |
1257 | "\tOOPSE will use a default value of\n" | |
1258 | < | "\t0.05 * electrostaticCutoffRadius\n" |
1013 | < | "\tfor the electrostaticSkinThickness\n"); |
1258 | > | "\t0.95 * cutoffRadius for the switchingRadius\n"); |
1259 | painCave.isFatal = 0; | |
1260 | simError(); | |
1261 | < | theEst = 0.05 * theEcr; |
1261 | > | theRsw = 0.95 * theRcut; |
1262 | } | |
1263 | else{ | |
1264 | < | theEst = globals->getEST(); |
1264 | > | theRsw = globals->getRsw(); |
1265 | } | |
1266 | ||
1267 | < | info[i].setDefaultEcr(theEcr, theEst); |
1267 | > | info[i].setDefaultRcut(theRcut, theRsw); |
1268 | ||
1269 | if (!globals->haveDielectric()){ | |
1270 | sprintf(painCave.errMsg, | |
# | Line 1033 | Line 1278 | void SimSetup::finalInfoCheck(void){ | |
1278 | } | |
1279 | else{ | |
1280 | if (usesDipoles || usesCharges){ | |
1281 | < | if (!globals->haveECR()){ |
1281 | > | |
1282 | > | if (!globals->haveRcut()){ |
1283 | sprintf(painCave.errMsg, | |
1284 | < | "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n" |
1284 | > | "SimSetup Warning: No value was set for the cutoffRadius.\n" |
1285 | "\tOOPSE will use a default value of 15.0 angstroms" | |
1286 | < | "\tfor the electrostaticCutoffRadius.\n"); |
1287 | < | painCave.isFatal = 0; |
1288 | < | simError(); |
1289 | < | theEcr = 15.0; |
1290 | < | } |
1286 | > | "\tfor the cutoffRadius.\n"); |
1287 | > | painCave.isFatal = 0; |
1288 | > | simError(); |
1289 | > | theRcut = 15.0; |
1290 | > | } |
1291 | else{ | |
1292 | < | theEcr = globals->getECR(); |
1292 | > | theRcut = globals->getRcut(); |
1293 | } | |
1294 | < | |
1295 | < | if (!globals->haveEST()){ |
1294 | > | |
1295 | > | if (!globals->haveRsw()){ |
1296 | sprintf(painCave.errMsg, | |
1297 | < | "SimSetup Warning: No value was set for electrostaticSkinThickness.\n" |
1297 | > | "SimSetup Warning: No value was set for switchingRadius.\n" |
1298 | "\tOOPSE will use a default value of\n" | |
1299 | < | "\t0.05 * electrostaticCutoffRadius\n" |
1054 | < | "\tfor the electrostaticSkinThickness\n"); |
1299 | > | "\t0.95 * cutoffRadius for the switchingRadius\n"); |
1300 | painCave.isFatal = 0; | |
1301 | simError(); | |
1302 | < | theEst = 0.05 * theEcr; |
1302 | > | theRsw = 0.95 * theRcut; |
1303 | } | |
1304 | else{ | |
1305 | < | theEst = globals->getEST(); |
1305 | > | theRsw = globals->getRsw(); |
1306 | } | |
1307 | + | |
1308 | + | info[i].setDefaultRcut(theRcut, theRsw); |
1309 | ||
1063 | – | info[i].setDefaultEcr(theEcr, theEst); |
1310 | } | |
1311 | } | |
1312 | } | |
# | Line 1068 | Line 1314 | void SimSetup::finalInfoCheck(void){ | |
1314 | strcpy(checkPointMsg, "post processing checks out"); | |
1315 | MPIcheckPoint(); | |
1316 | #endif // is_mpi | |
1317 | + | |
1318 | + | // clean up the forcefield |
1319 | + | the_ff->cleanMe(); |
1320 | } | |
1321 | ||
1322 | void SimSetup::initSystemCoords(void){ | |
# | Line 1195 | Line 1444 | void SimSetup::makeOutNames(void){ | |
1444 | } | |
1445 | else{ | |
1446 | strcat(info[k].statusName, ".stat"); | |
1447 | + | } |
1448 | + | } |
1449 | + | |
1450 | + | strcpy(info[k].rawPotName, inFileName); |
1451 | + | nameLength = strlen(info[k].rawPotName); |
1452 | + | endTest = &(info[k].rawPotName[nameLength - 5]); |
1453 | + | if (!strcmp(endTest, ".bass")){ |
1454 | + | strcpy(endTest, ".raw"); |
1455 | + | } |
1456 | + | else if (!strcmp(endTest, ".BASS")){ |
1457 | + | strcpy(endTest, ".raw"); |
1458 | + | } |
1459 | + | else{ |
1460 | + | endTest = &(info[k].rawPotName[nameLength - 4]); |
1461 | + | if (!strcmp(endTest, ".bss")){ |
1462 | + | strcpy(endTest, ".raw"); |
1463 | } | |
1464 | + | else if (!strcmp(endTest, ".mdl")){ |
1465 | + | strcpy(endTest, ".raw"); |
1466 | + | } |
1467 | + | else{ |
1468 | + | strcat(info[k].rawPotName, ".raw"); |
1469 | + | } |
1470 | } | |
1471 | ||
1472 | #ifdef IS_MPI | |
# | Line 1247 | Line 1518 | void SimSetup::createFF(void){ | |
1518 | void SimSetup::createFF(void){ | |
1519 | switch (ffCase){ | |
1520 | case FF_DUFF: | |
1521 | < | the_ff = new DUFF(); |
1521 | > | the_ff = new DUFF(); |
1522 | break; | |
1523 | ||
1524 | case FF_LJ: | |
# | Line 1255 | Line 1526 | void SimSetup::createFF(void){ | |
1526 | break; | |
1527 | ||
1528 | case FF_EAM: | |
1529 | < | the_ff = new EAM_FF(); |
1529 | > | if (has_forcefield_variant) |
1530 | > | the_ff = new EAM_FF(forcefield_variant); |
1531 | > | else |
1532 | > | the_ff = new EAM_FF(); |
1533 | break; | |
1534 | ||
1535 | case FF_H2O: | |
# | Line 1269 | Line 1543 | void SimSetup::createFF(void){ | |
1543 | simError(); | |
1544 | } | |
1545 | ||
1546 | + | |
1547 | #ifdef IS_MPI | |
1548 | strcpy(checkPointMsg, "ForceField creation successful"); | |
1549 | MPIcheckPoint(); | |
# | Line 1282 | Line 1557 | void SimSetup::compList(void){ | |
1557 | LinkedMolStamp* headStamp = new LinkedMolStamp(); | |
1558 | LinkedMolStamp* currentStamp = NULL; | |
1559 | comp_stamps = new MoleculeStamp * [n_components]; | |
1560 | + | bool haveCutoffGroups; |
1561 | ||
1562 | + | haveCutoffGroups = false; |
1563 | + | |
1564 | // make an array of molecule stamps that match the components used. | |
1565 | // also extract the used stamps out into a separate linked list | |
1566 | ||
# | Line 1317 | Line 1595 | void SimSetup::compList(void){ | |
1595 | headStamp->add(currentStamp); | |
1596 | comp_stamps[i] = headStamp->match(id); | |
1597 | } | |
1598 | + | |
1599 | + | if(comp_stamps[i]->getNCutoffGroups() > 0) |
1600 | + | haveCutoffGroups = true; |
1601 | } | |
1602 | + | |
1603 | + | for (i = 0; i < nInfo; i++) |
1604 | + | info[i].haveCutoffGroups = haveCutoffGroups; |
1605 | ||
1606 | #ifdef IS_MPI | |
1607 | strcpy(checkPointMsg, "Component stamps successfully extracted\n"); | |
# | Line 1326 | Line 1610 | void SimSetup::calcSysValues(void){ | |
1610 | } | |
1611 | ||
1612 | void SimSetup::calcSysValues(void){ | |
1613 | < | int i; |
1613 | > | int i, j; |
1614 | > | int ncutgroups, atomsingroups, ngroupsinstamp; |
1615 | ||
1616 | int* molMembershipArray; | |
1617 | + | CutoffGroupStamp* cg; |
1618 | ||
1619 | tot_atoms = 0; | |
1620 | tot_bonds = 0; | |
1621 | tot_bends = 0; | |
1622 | tot_torsions = 0; | |
1623 | tot_rigid = 0; | |
1624 | + | tot_groups = 0; |
1625 | for (i = 0; i < n_components; i++){ | |
1626 | tot_atoms += components_nmol[i] * comp_stamps[i]->getNAtoms(); | |
1627 | tot_bonds += components_nmol[i] * comp_stamps[i]->getNBonds(); | |
1628 | tot_bends += components_nmol[i] * comp_stamps[i]->getNBends(); | |
1629 | tot_torsions += components_nmol[i] * comp_stamps[i]->getNTorsions(); | |
1630 | tot_rigid += components_nmol[i] * comp_stamps[i]->getNRigidBodies(); | |
1631 | + | |
1632 | + | ncutgroups = comp_stamps[i]->getNCutoffGroups(); |
1633 | + | atomsingroups = 0; |
1634 | + | for (j=0; j < ncutgroups; j++) { |
1635 | + | cg = comp_stamps[i]->getCutoffGroup(j); |
1636 | + | atomsingroups += cg->getNMembers(); |
1637 | + | } |
1638 | + | ngroupsinstamp = comp_stamps[i]->getNAtoms() - atomsingroups + ncutgroups; |
1639 | + | tot_groups += components_nmol[i] * ngroupsinstamp; |
1640 | } | |
1641 | ||
1642 | tot_SRI = tot_bonds + tot_bends + tot_torsions; | |
# | Line 1353 | Line 1649 | void SimSetup::calcSysValues(void){ | |
1649 | info[i].n_torsions = tot_torsions; | |
1650 | info[i].n_SRI = tot_SRI; | |
1651 | info[i].n_mol = tot_nmol; | |
1652 | < | |
1652 | > | info[i].ngroup = tot_groups; |
1653 | info[i].molMembershipArray = molMembershipArray; | |
1654 | } | |
1655 | } | |
# | Line 1364 | Line 1660 | void SimSetup::mpiMolDivide(void){ | |
1660 | int i, j, k; | |
1661 | int localMol, allMol; | |
1662 | int local_atoms, local_bonds, local_bends, local_torsions, local_SRI; | |
1663 | < | int local_rigid; |
1663 | > | int local_rigid, local_groups; |
1664 | vector<int> globalMolIndex; | |
1665 | + | int ncutgroups, atomsingroups, ngroupsinstamp; |
1666 | + | CutoffGroupStamp* cg; |
1667 | ||
1668 | mpiSim = new mpiSimulation(info); | |
1669 | ||
1670 | mpiSim->divideLabor(); | |
1671 | globalAtomIndex = mpiSim->getGlobalAtomIndex(); | |
1672 | + | globalGroupIndex = mpiSim->getGlobalGroupIndex(); |
1673 | //globalMolIndex = mpiSim->getGlobalMolIndex(); | |
1674 | ||
1675 | // set up the local variables | |
# | Line 1385 | Line 1684 | void SimSetup::mpiMolDivide(void){ | |
1684 | local_bends = 0; | |
1685 | local_torsions = 0; | |
1686 | local_rigid = 0; | |
1687 | + | local_groups = 0; |
1688 | globalAtomCounter = 0; | |
1689 | ||
1690 | for (i = 0; i < n_components; i++){ | |
# | Line 1395 | Line 1695 | void SimSetup::mpiMolDivide(void){ | |
1695 | local_bends += comp_stamps[i]->getNBends(); | |
1696 | local_torsions += comp_stamps[i]->getNTorsions(); | |
1697 | local_rigid += comp_stamps[i]->getNRigidBodies(); | |
1698 | + | |
1699 | + | ncutgroups = comp_stamps[i]->getNCutoffGroups(); |
1700 | + | atomsingroups = 0; |
1701 | + | for (k=0; k < ncutgroups; k++) { |
1702 | + | cg = comp_stamps[i]->getCutoffGroup(k); |
1703 | + | atomsingroups += cg->getNMembers(); |
1704 | + | } |
1705 | + | ngroupsinstamp = comp_stamps[i]->getNAtoms() - atomsingroups + |
1706 | + | ncutgroups; |
1707 | + | local_groups += ngroupsinstamp; |
1708 | + | |
1709 | localMol++; | |
1710 | } | |
1711 | for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){ | |
# | Line 1407 | Line 1718 | void SimSetup::mpiMolDivide(void){ | |
1718 | } | |
1719 | local_SRI = local_bonds + local_bends + local_torsions; | |
1720 | ||
1721 | < | info[0].n_atoms = mpiSim->getMyNlocal(); |
1721 | > | info[0].n_atoms = mpiSim->getNAtomsLocal(); |
1722 | ||
1412 | – | |
1723 | if (local_atoms != info[0].n_atoms){ | |
1724 | sprintf(painCave.errMsg, | |
1725 | "SimSetup error: mpiSim's localAtom (%d) and SimSetup's\n" | |
# | Line 1419 | Line 1729 | void SimSetup::mpiMolDivide(void){ | |
1729 | simError(); | |
1730 | } | |
1731 | ||
1732 | + | info[0].ngroup = mpiSim->getNGroupsLocal(); |
1733 | + | if (local_groups != info[0].ngroup){ |
1734 | + | sprintf(painCave.errMsg, |
1735 | + | "SimSetup error: mpiSim's localGroups (%d) and SimSetup's\n" |
1736 | + | "\tlocalGroups (%d) are not equal.\n", |
1737 | + | info[0].ngroup, local_groups); |
1738 | + | painCave.isFatal = 1; |
1739 | + | simError(); |
1740 | + | } |
1741 | + | |
1742 | info[0].n_bonds = local_bonds; | |
1743 | info[0].n_bends = local_bends; | |
1744 | info[0].n_torsions = local_torsions; | |
# | Line 1455 | Line 1775 | void SimSetup::makeSysArrays(void){ | |
1775 | ||
1776 | ||
1777 | molIndex = 0; | |
1778 | < | for (i = 0; i < mpiSim->getTotNmol(); i++){ |
1778 | > | for (i = 0; i < mpiSim->getNMolGlobal(); i++){ |
1779 | if (mol2proc[i] == worldRank){ | |
1780 | the_molecules[molIndex].setStampID(molCompType[i]); | |
1781 | the_molecules[molIndex].setMyIndex(molIndex); | |
# | Line 1492 | Line 1812 | void SimSetup::makeSysArrays(void){ | |
1812 | info[l].atoms = the_atoms; | |
1813 | info[l].molecules = the_molecules; | |
1814 | info[l].nGlobalExcludes = 0; | |
1815 | < | |
1815 | > | |
1816 | the_ff->setSimInfo(info); | |
1817 | } | |
1818 | } | |
# | Line 1500 | Line 1820 | void SimSetup::makeIntegrator(void){ | |
1820 | void SimSetup::makeIntegrator(void){ | |
1821 | int k; | |
1822 | ||
1823 | < | NVE<RealIntegrator>* myNVE = NULL; |
1824 | < | NVT<RealIntegrator>* myNVT = NULL; |
1825 | < | NPTi<NPT<RealIntegrator> >* myNPTi = NULL; |
1826 | < | NPTf<NPT<RealIntegrator> >* myNPTf = NULL; |
1827 | < | NPTxyz<NPT<RealIntegrator> >* myNPTxyz = NULL; |
1823 | > | NVE<Integrator<BaseIntegrator> >* myNVE = NULL; |
1824 | > | NVT<Integrator<BaseIntegrator> >* myNVT = NULL; |
1825 | > | NPTi<NPT<Integrator<BaseIntegrator> > >* myNPTi = NULL; |
1826 | > | NPTf<NPT<Integrator<BaseIntegrator> > >* myNPTf = NULL; |
1827 | > | NPTxyz<NPT<Integrator<BaseIntegrator> > >* myNPTxyz = NULL; |
1828 | ||
1829 | for (k = 0; k < nInfo; k++){ | |
1830 | switch (ensembleCase){ | |
# | Line 1514 | Line 1834 | void SimSetup::makeIntegrator(void){ | |
1834 | myNVE = new ZConstraint<NVE<RealIntegrator> >(&(info[k]), the_ff); | |
1835 | } | |
1836 | else{ | |
1837 | < | myNVE = new NVE<RealIntegrator>(&(info[k]), the_ff); |
1837 | > | if (globals->haveQuaternion()){ |
1838 | > | if (globals->getUseQuaternion()) |
1839 | > | info->the_integrator = new NVE<SQSIntegrator<RealIntegrator> >(&(info[k]), the_ff); |
1840 | > | } |
1841 | > | else |
1842 | > | info->the_integrator = new NVE<RealIntegrator>(&(info[k]), the_ff); |
1843 | > | break; |
1844 | > | |
1845 | > | //myNVE = new NVE<RealIntegrator>(&(info[k]), the_ff); |
1846 | } | |
1847 | ||
1848 | info->the_integrator = myNVE; |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |