# | 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 147 | Line 149 | void SimSetup::createSim(void){ | |
149 | // make the output filenames | |
150 | ||
151 | makeOutNames(); | |
150 | – | |
151 | – | if (globals->haveMinimizer()) |
152 | – | // make minimizer |
153 | – | makeMinimizer(); |
154 | – | else |
155 | – | // make the integrator |
156 | – | makeIntegrator(); |
152 | ||
153 | #ifdef IS_MPI | |
154 | mpiSim->mpiRefresh(); | |
# | Line 162 | Line 157 | void SimSetup::createSim(void){ | |
157 | // initialize the Fortran | |
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(); |
168 | + | else |
169 | + | // make the integrator |
170 | + | makeIntegrator(); |
171 | + | |
172 | } | |
173 | ||
174 | ||
# | Line 169 | 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 182 | 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 193 | 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 222 | void SimSetup::makeMolecules(void){ | |
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 215 | 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 262 | Line 293 | void SimSetup::makeMolecules(void){ | |
293 | else{ | |
294 | ||
295 | molInfo.myAtoms[j] = new Atom((j + atomOffset), info[k].getConfiguration()); | |
296 | + | |
297 | } | |
298 | ||
299 | molInfo.myAtoms[j]->setType(currentAtom->getType()); | |
268 | – | |
300 | #ifdef IS_MPI | |
270 | – | |
301 | molInfo.myAtoms[j]->setGlobalIndex(globalAtomIndex[j + atomOffset]); | |
272 | – | |
302 | #endif // is_mpi | |
303 | } | |
304 | ||
# | Line 409 | Line 438 | void SimSetup::makeMolecules(void){ | |
438 | info[k].excludes->addPair(exK, exL); | |
439 | } | |
440 | ||
441 | + | |
442 | + | molInfo.myRigidBodies.clear(); |
443 | + | |
444 | for (j = 0; j < molInfo.nRigidBodies; j++){ | |
445 | ||
446 | currentRigidBody = comp_stamps[stampID]->getRigidBody(j); | |
# | Line 418 | Line 450 | void SimSetup::makeMolecules(void){ | |
450 | ||
451 | myRB = new RigidBody(); | |
452 | ||
453 | < | sprintf(rbName,"%s_RB_%s", molName, j); |
453 | > | sprintf(rbName,"%s_RB_%d", molName, j); |
454 | myRB->setType(rbName); | |
455 | ||
456 | for (rb1 = 0; rb1 < nMembers; rb1++) { | |
# | Line 460 | Line 492 | void SimSetup::makeMolecules(void){ | |
492 | // used for the exclude list: | |
493 | ||
494 | #ifdef IS_MPI | |
495 | < | exI = info[k].atoms[tempI]->getGlobalIndex() + 1; |
496 | < | exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1; |
495 | > | exI = molInfo.myAtoms[tempI]->getGlobalIndex() + 1; |
496 | > | exJ = molInfo.myAtoms[tempJ]->getGlobalIndex() + 1; |
497 | #else | |
498 | < | exI = tempI + 1; |
499 | < | exJ = tempJ + 1; |
498 | > | exI = molInfo.myAtoms[tempI]->getIndex() + 1; |
499 | > | exJ = molInfo.myAtoms[tempJ]->getIndex() + 1; |
500 | #endif | |
501 | ||
502 | info[k].excludes->addPair(exI, exJ); | |
# | Line 476 | Line 508 | void SimSetup::makeMolecules(void){ | |
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: | |
580 | ||
581 | + | molInfo.myIntegrableObjects.clear(); |
582 | + | |
583 | + | |
584 | for (j = 0; j < molInfo.nAtoms; j++){ | |
585 | ||
586 | #ifdef IS_MPI | |
# | Line 504 | Line 605 | void SimSetup::makeMolecules(void){ | |
605 | info[k].integrableObjects.push_back(mySD); | |
606 | molInfo.myIntegrableObjects.push_back(mySD); | |
607 | } | |
608 | < | |
508 | < | |
608 | > | |
609 | // send the arrays off to the forceField for init. | |
610 | ||
611 | the_ff->initializeAtoms(molInfo.nAtoms, molInfo.myAtoms); | |
# | Line 514 | Line 614 | void SimSetup::makeMolecules(void){ | |
614 | the_ff->initializeTorsions(molInfo.nTorsions, molInfo.myTorsions, | |
615 | theTorsions); | |
616 | ||
517 | – | 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 | < | } |
525 | < | } |
676 | > | } |
677 | ||
527 | – | #ifdef IS_MPI |
528 | – | sprintf(checkPointMsg, "all molecules initialized succesfully"); |
529 | – | MPIcheckPoint(); |
530 | – | #endif // is_mpi |
678 | ||
532 | – | // 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 { |
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 | < | the_ff->setRcut( globals->getLJrcut() ); |
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 827 | Line 993 | void SimSetup::gatherInfo(void){ | |
993 | } | |
994 | ||
995 | //check whether sample time, status time, thermal time and reset time are divisble by dt | |
996 | < | if (!isDivisible(globals->getSampleTime(), globals->getDt())){ |
996 | > | if (globals->haveSampleTime() && !isDivisible(globals->getSampleTime(), globals->getDt())){ |
997 | sprintf(painCave.errMsg, | |
998 | "Sample time is not divisible by dt.\n" | |
999 | "\tThis will result in samples that are not uniformly\n" | |
# | Line 837 | Line 1003 | void SimSetup::gatherInfo(void){ | |
1003 | simError(); | |
1004 | } | |
1005 | ||
1006 | < | if (globals->haveStatusTime() && !isDivisible(globals->getSampleTime(), globals->getDt())){ |
1006 | > | if (globals->haveStatusTime() && !isDivisible(globals->getStatusTime(), globals->getDt())){ |
1007 | sprintf(painCave.errMsg, | |
1008 | "Status time is not divisible by dt.\n" | |
1009 | "\tThis will result in status reports that are not uniformly\n" | |
# | Line 873 | Line 1039 | void SimSetup::gatherInfo(void){ | |
1039 | if (globals->haveSampleTime()){ | |
1040 | info[i].sampleTime = globals->getSampleTime(); | |
1041 | info[i].statusTime = info[i].sampleTime; | |
876 | – | info[i].thermalTime = info[i].sampleTime; |
1042 | } | |
1043 | else{ | |
1044 | info[i].sampleTime = globals->getRunTime(); | |
1045 | info[i].statusTime = info[i].sampleTime; | |
881 | – | info[i].thermalTime = info[i].sampleTime; |
1046 | } | |
1047 | ||
1048 | if (globals->haveStatusTime()){ | |
# | Line 887 | Line 1051 | void SimSetup::gatherInfo(void){ | |
1051 | ||
1052 | if (globals->haveThermalTime()){ | |
1053 | info[i].thermalTime = globals->getThermalTime(); | |
1054 | + | } else { |
1055 | + | info[i].thermalTime = globals->getRunTime(); |
1056 | } | |
1057 | ||
1058 | info[i].resetIntegrator = 0; | |
# | Line 904 | Line 1070 | void SimSetup::gatherInfo(void){ | |
1070 | ||
1071 | info[i].useInitXSstate = globals->getUseInitXSstate(); | |
1072 | info[i].orthoTolerance = globals->getOrthoBoxTolerance(); | |
1073 | < | |
1073 | > | |
1074 | > | // check for thermodynamic integration |
1075 | > | if (globals->getUseSolidThermInt() && !globals->getUseLiquidThermInt()) { |
1076 | > | if (globals->haveThermIntLambda() && globals->haveThermIntK()) { |
1077 | > | info[i].useSolidThermInt = globals->getUseSolidThermInt(); |
1078 | > | info[i].thermIntLambda = globals->getThermIntLambda(); |
1079 | > | info[i].thermIntK = globals->getThermIntK(); |
1080 | > | |
1081 | > | Restraints *myRestraint = new Restraints(tot_nmol, info[i].thermIntLambda, info[i].thermIntK); |
1082 | > | info[i].restraint = myRestraint; |
1083 | > | } |
1084 | > | else { |
1085 | > | sprintf(painCave.errMsg, |
1086 | > | "SimSetup Error:\n" |
1087 | > | "\tKeyword useSolidThermInt was set to 'true' but\n" |
1088 | > | "\tthermodynamicIntegrationLambda (and/or\n" |
1089 | > | "\tthermodynamicIntegrationK) was not specified.\n" |
1090 | > | "\tPlease provide a lambda value and k value in your .bass file.\n"); |
1091 | > | painCave.isFatal = 1; |
1092 | > | simError(); |
1093 | > | } |
1094 | > | } |
1095 | > | else if(globals->getUseLiquidThermInt()) { |
1096 | > | if (globals->getUseSolidThermInt()) { |
1097 | > | sprintf( painCave.errMsg, |
1098 | > | "SimSetup Warning: It appears that you have both solid and\n" |
1099 | > | "\tliquid thermodynamic integration activated in your .bass\n" |
1100 | > | "\tfile. To avoid confusion, specify only one technique in\n" |
1101 | > | "\tyour .bass file. Liquid-state thermodynamic integration\n" |
1102 | > | "\twill be assumed for the current simulation. If this is not\n" |
1103 | > | "\twhat you desire, set useSolidThermInt to 'true' and\n" |
1104 | > | "\tuseLiquidThermInt to 'false' in your .bass file.\n"); |
1105 | > | painCave.isFatal = 0; |
1106 | > | simError(); |
1107 | > | } |
1108 | > | if (globals->haveThermIntLambda() && globals->haveThermIntK()) { |
1109 | > | info[i].useLiquidThermInt = globals->getUseLiquidThermInt(); |
1110 | > | info[i].thermIntLambda = globals->getThermIntLambda(); |
1111 | > | info[i].thermIntK = globals->getThermIntK(); |
1112 | > | } |
1113 | > | else { |
1114 | > | sprintf(painCave.errMsg, |
1115 | > | "SimSetup Error:\n" |
1116 | > | "\tKeyword useLiquidThermInt was set to 'true' but\n" |
1117 | > | "\tthermodynamicIntegrationLambda (and/or\n" |
1118 | > | "\tthermodynamicIntegrationK) was not specified.\n" |
1119 | > | "\tPlease provide a lambda value and k value in your .bass file.\n"); |
1120 | > | painCave.isFatal = 1; |
1121 | > | simError(); |
1122 | > | } |
1123 | > | } |
1124 | > | else if(globals->haveThermIntLambda() || globals->haveThermIntK()){ |
1125 | > | sprintf(painCave.errMsg, |
1126 | > | "SimSetup Warning: If you want to use Thermodynamic\n" |
1127 | > | "\tIntegration, set useSolidThermInt or useLiquidThermInt to\n" |
1128 | > | "\t'true' in your .bass file. These keywords are set to\n" |
1129 | > | "\t'false' by default, so your lambda and/or k values are\n" |
1130 | > | "\tbeing ignored.\n"); |
1131 | > | painCave.isFatal = 0; |
1132 | > | simError(); |
1133 | > | } |
1134 | } | |
1135 | ||
1136 | //setup seed for random number generator | |
# | Line 957 | Line 1183 | void SimSetup::finalInfoCheck(void){ | |
1183 | void SimSetup::finalInfoCheck(void){ | |
1184 | int index; | |
1185 | int usesDipoles; | |
1186 | + | int usesCharges; |
1187 | int i; | |
1188 | ||
1189 | for (i = 0; i < nInfo; i++){ | |
# | Line 968 | Line 1195 | void SimSetup::finalInfoCheck(void){ | |
1195 | usesDipoles = (info[i].atoms[index])->hasDipole(); | |
1196 | index++; | |
1197 | } | |
1198 | < | |
1198 | > | index = 0; |
1199 | > | usesCharges = 0; |
1200 | > | while ((index < info[i].n_atoms) && !usesCharges){ |
1201 | > | usesCharges= (info[i].atoms[index])->hasCharge(); |
1202 | > | index++; |
1203 | > | } |
1204 | #ifdef IS_MPI | |
1205 | int myUse = usesDipoles; | |
1206 | MPI_Allreduce(&myUse, &usesDipoles, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | |
1207 | #endif //is_mpi | |
1208 | ||
1209 | < | double theEcr, theEst; |
1209 | > | double theRcut, theRsw; |
1210 | ||
1211 | + | if (globals->haveRcut()) { |
1212 | + | theRcut = globals->getRcut(); |
1213 | + | |
1214 | + | if (globals->haveRsw()) |
1215 | + | theRsw = globals->getRsw(); |
1216 | + | else |
1217 | + | theRsw = theRcut; |
1218 | + | |
1219 | + | info[i].setDefaultRcut(theRcut, theRsw); |
1220 | + | |
1221 | + | } else { |
1222 | + | |
1223 | + | the_ff->calcRcut(); |
1224 | + | theRcut = info[i].getRcut(); |
1225 | + | |
1226 | + | if (globals->haveRsw()) |
1227 | + | theRsw = globals->getRsw(); |
1228 | + | else |
1229 | + | theRsw = theRcut; |
1230 | + | |
1231 | + | info[i].setDefaultRcut(theRcut, theRsw); |
1232 | + | } |
1233 | + | |
1234 | if (globals->getUseRF()){ | |
1235 | info[i].useReactionField = 1; | |
1236 | < | |
1237 | < | if (!globals->haveECR()){ |
1236 | > | |
1237 | > | if (!globals->haveRcut()){ |
1238 | sprintf(painCave.errMsg, | |
1239 | < | "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n" |
1239 | > | "SimSetup Warning: No value was set for the cutoffRadius.\n" |
1240 | "\tOOPSE will use a default value of 15.0 angstroms" | |
1241 | < | "\tfor the electrostaticCutoffRadius.\n"); |
1241 | > | "\tfor the cutoffRadius.\n"); |
1242 | painCave.isFatal = 0; | |
1243 | simError(); | |
1244 | < | theEcr = 15.0; |
1244 | > | theRcut = 15.0; |
1245 | } | |
1246 | else{ | |
1247 | < | theEcr = globals->getECR(); |
1247 | > | theRcut = globals->getRcut(); |
1248 | } | |
1249 | ||
1250 | < | if (!globals->haveEST()){ |
1250 | > | if (!globals->haveRsw()){ |
1251 | sprintf(painCave.errMsg, | |
1252 | < | "SimSetup Warning: No value was set for electrostaticSkinThickness.\n" |
1252 | > | "SimSetup Warning: No value was set for switchingRadius.\n" |
1253 | "\tOOPSE will use a default value of\n" | |
1254 | < | "\t0.05 * electrostaticCutoffRadius\n" |
1000 | < | "\tfor the electrostaticSkinThickness\n"); |
1254 | > | "\t0.95 * cutoffRadius for the switchingRadius\n"); |
1255 | painCave.isFatal = 0; | |
1256 | simError(); | |
1257 | < | theEst = 0.05 * theEcr; |
1257 | > | theRsw = 0.95 * theRcut; |
1258 | } | |
1259 | else{ | |
1260 | < | theEst = globals->getEST(); |
1260 | > | theRsw = globals->getRsw(); |
1261 | } | |
1262 | ||
1263 | < | info[i].setDefaultEcr(theEcr, theEst); |
1263 | > | info[i].setDefaultRcut(theRcut, theRsw); |
1264 | ||
1265 | if (!globals->haveDielectric()){ | |
1266 | sprintf(painCave.errMsg, | |
# | Line 1019 | Line 1273 | void SimSetup::finalInfoCheck(void){ | |
1273 | info[i].dielectric = globals->getDielectric(); | |
1274 | } | |
1275 | else{ | |
1276 | < | if (usesDipoles){ |
1277 | < | if (!globals->haveECR()){ |
1276 | > | if (usesDipoles || usesCharges){ |
1277 | > | |
1278 | > | if (!globals->haveRcut()){ |
1279 | sprintf(painCave.errMsg, | |
1280 | < | "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n" |
1280 | > | "SimSetup Warning: No value was set for the cutoffRadius.\n" |
1281 | "\tOOPSE will use a default value of 15.0 angstroms" | |
1282 | < | "\tfor the electrostaticCutoffRadius.\n"); |
1283 | < | painCave.isFatal = 0; |
1284 | < | simError(); |
1285 | < | theEcr = 15.0; |
1286 | < | } |
1287 | < | else{ |
1288 | < | theEcr = globals->getECR(); |
1282 | > | "\tfor the cutoffRadius.\n"); |
1283 | > | painCave.isFatal = 0; |
1284 | > | simError(); |
1285 | > | theRcut = 15.0; |
1286 | > | } |
1287 | > | else{ |
1288 | > | theRcut = globals->getRcut(); |
1289 | } | |
1290 | < | |
1291 | < | if (!globals->haveEST()){ |
1290 | > | |
1291 | > | if (!globals->haveRsw()){ |
1292 | sprintf(painCave.errMsg, | |
1293 | < | "SimSetup Warning: No value was set for electrostaticSkinThickness.\n" |
1293 | > | "SimSetup Warning: No value was set for switchingRadius.\n" |
1294 | "\tOOPSE will use a default value of\n" | |
1295 | < | "\t0.05 * electrostaticCutoffRadius\n" |
1041 | < | "\tfor the electrostaticSkinThickness\n"); |
1295 | > | "\t0.95 * cutoffRadius for the switchingRadius\n"); |
1296 | painCave.isFatal = 0; | |
1297 | simError(); | |
1298 | < | theEst = 0.05 * theEcr; |
1298 | > | theRsw = 0.95 * theRcut; |
1299 | } | |
1300 | else{ | |
1301 | < | theEst = globals->getEST(); |
1301 | > | theRsw = globals->getRsw(); |
1302 | } | |
1303 | + | |
1304 | + | info[i].setDefaultRcut(theRcut, theRsw); |
1305 | ||
1050 | – | info[i].setDefaultEcr(theEcr, theEst); |
1306 | } | |
1307 | } | |
1308 | } | |
# | Line 1055 | Line 1310 | void SimSetup::finalInfoCheck(void){ | |
1310 | strcpy(checkPointMsg, "post processing checks out"); | |
1311 | MPIcheckPoint(); | |
1312 | #endif // is_mpi | |
1313 | + | |
1314 | + | // clean up the forcefield |
1315 | + | the_ff->cleanMe(); |
1316 | } | |
1317 | ||
1318 | void SimSetup::initSystemCoords(void){ | |
# | Line 1185 | Line 1443 | void SimSetup::makeOutNames(void){ | |
1443 | } | |
1444 | } | |
1445 | ||
1446 | + | strcpy(info[k].rawPotName, inFileName); |
1447 | + | nameLength = strlen(info[k].rawPotName); |
1448 | + | endTest = &(info[k].rawPotName[nameLength - 5]); |
1449 | + | if (!strcmp(endTest, ".bass")){ |
1450 | + | strcpy(endTest, ".raw"); |
1451 | + | } |
1452 | + | else if (!strcmp(endTest, ".BASS")){ |
1453 | + | strcpy(endTest, ".raw"); |
1454 | + | } |
1455 | + | else{ |
1456 | + | endTest = &(info[k].rawPotName[nameLength - 4]); |
1457 | + | if (!strcmp(endTest, ".bss")){ |
1458 | + | strcpy(endTest, ".raw"); |
1459 | + | } |
1460 | + | else if (!strcmp(endTest, ".mdl")){ |
1461 | + | strcpy(endTest, ".raw"); |
1462 | + | } |
1463 | + | else{ |
1464 | + | strcat(info[k].rawPotName, ".raw"); |
1465 | + | } |
1466 | + | } |
1467 | + | |
1468 | #ifdef IS_MPI | |
1469 | ||
1470 | } | |
# | Line 1269 | Line 1549 | void SimSetup::compList(void){ | |
1549 | LinkedMolStamp* headStamp = new LinkedMolStamp(); | |
1550 | LinkedMolStamp* currentStamp = NULL; | |
1551 | comp_stamps = new MoleculeStamp * [n_components]; | |
1552 | + | bool haveCutoffGroups; |
1553 | ||
1554 | + | haveCutoffGroups = false; |
1555 | + | |
1556 | // make an array of molecule stamps that match the components used. | |
1557 | // also extract the used stamps out into a separate linked list | |
1558 | ||
# | Line 1304 | Line 1587 | void SimSetup::compList(void){ | |
1587 | headStamp->add(currentStamp); | |
1588 | comp_stamps[i] = headStamp->match(id); | |
1589 | } | |
1590 | + | |
1591 | + | if(comp_stamps[i]->getNCutoffGroups() > 0) |
1592 | + | haveCutoffGroups = true; |
1593 | } | |
1594 | + | |
1595 | + | for (i = 0; i < nInfo; i++) |
1596 | + | info[i].haveCutoffGroups = haveCutoffGroups; |
1597 | ||
1598 | #ifdef IS_MPI | |
1599 | strcpy(checkPointMsg, "Component stamps successfully extracted\n"); | |
# | Line 1313 | Line 1602 | void SimSetup::calcSysValues(void){ | |
1602 | } | |
1603 | ||
1604 | void SimSetup::calcSysValues(void){ | |
1605 | < | int i; |
1605 | > | int i, j; |
1606 | > | int ncutgroups, atomsingroups, ngroupsinstamp; |
1607 | ||
1608 | int* molMembershipArray; | |
1609 | + | CutoffGroupStamp* cg; |
1610 | ||
1611 | tot_atoms = 0; | |
1612 | tot_bonds = 0; | |
1613 | tot_bends = 0; | |
1614 | tot_torsions = 0; | |
1615 | tot_rigid = 0; | |
1616 | + | tot_groups = 0; |
1617 | for (i = 0; i < n_components; i++){ | |
1618 | tot_atoms += components_nmol[i] * comp_stamps[i]->getNAtoms(); | |
1619 | tot_bonds += components_nmol[i] * comp_stamps[i]->getNBonds(); | |
1620 | tot_bends += components_nmol[i] * comp_stamps[i]->getNBends(); | |
1621 | tot_torsions += components_nmol[i] * comp_stamps[i]->getNTorsions(); | |
1622 | tot_rigid += components_nmol[i] * comp_stamps[i]->getNRigidBodies(); | |
1623 | + | |
1624 | + | ncutgroups = comp_stamps[i]->getNCutoffGroups(); |
1625 | + | atomsingroups = 0; |
1626 | + | for (j=0; j < ncutgroups; j++) { |
1627 | + | cg = comp_stamps[i]->getCutoffGroup(j); |
1628 | + | atomsingroups += cg->getNMembers(); |
1629 | + | } |
1630 | + | ngroupsinstamp = comp_stamps[i]->getNAtoms() - atomsingroups + ncutgroups; |
1631 | + | tot_groups += components_nmol[i] * ngroupsinstamp; |
1632 | } | |
1633 | ||
1634 | tot_SRI = tot_bonds + tot_bends + tot_torsions; | |
# | Line 1340 | Line 1641 | void SimSetup::calcSysValues(void){ | |
1641 | info[i].n_torsions = tot_torsions; | |
1642 | info[i].n_SRI = tot_SRI; | |
1643 | info[i].n_mol = tot_nmol; | |
1644 | < | |
1644 | > | info[i].ngroup = tot_groups; |
1645 | info[i].molMembershipArray = molMembershipArray; | |
1646 | } | |
1647 | } | |
# | Line 1351 | Line 1652 | void SimSetup::mpiMolDivide(void){ | |
1652 | int i, j, k; | |
1653 | int localMol, allMol; | |
1654 | int local_atoms, local_bonds, local_bends, local_torsions, local_SRI; | |
1655 | < | int local_rigid; |
1355 | < | vector<int> globalAtomIndex; |
1655 | > | int local_rigid, local_groups; |
1656 | vector<int> globalMolIndex; | |
1657 | + | int ncutgroups, atomsingroups, ngroupsinstamp; |
1658 | + | CutoffGroupStamp* cg; |
1659 | ||
1660 | mpiSim = new mpiSimulation(info); | |
1661 | ||
1662 | mpiSim->divideLabor(); | |
1663 | globalAtomIndex = mpiSim->getGlobalAtomIndex(); | |
1664 | < | globalMolIndex = mpiSim->getGlobalMolIndex(); |
1664 | > | globalGroupIndex = mpiSim->getGlobalGroupIndex(); |
1665 | > | //globalMolIndex = mpiSim->getGlobalMolIndex(); |
1666 | ||
1667 | // set up the local variables | |
1668 | ||
# | Line 1373 | Line 1676 | void SimSetup::mpiMolDivide(void){ | |
1676 | local_bends = 0; | |
1677 | local_torsions = 0; | |
1678 | local_rigid = 0; | |
1679 | + | local_groups = 0; |
1680 | globalAtomCounter = 0; | |
1681 | ||
1682 | for (i = 0; i < n_components; i++){ | |
# | Line 1383 | Line 1687 | void SimSetup::mpiMolDivide(void){ | |
1687 | local_bends += comp_stamps[i]->getNBends(); | |
1688 | local_torsions += comp_stamps[i]->getNTorsions(); | |
1689 | local_rigid += comp_stamps[i]->getNRigidBodies(); | |
1690 | + | |
1691 | + | ncutgroups = comp_stamps[i]->getNCutoffGroups(); |
1692 | + | atomsingroups = 0; |
1693 | + | for (k=0; k < ncutgroups; k++) { |
1694 | + | cg = comp_stamps[i]->getCutoffGroup(k); |
1695 | + | atomsingroups += cg->getNMembers(); |
1696 | + | } |
1697 | + | ngroupsinstamp = comp_stamps[i]->getNAtoms() - atomsingroups + |
1698 | + | ncutgroups; |
1699 | + | local_groups += ngroupsinstamp; |
1700 | + | |
1701 | localMol++; | |
1702 | } | |
1703 | for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){ | |
# | Line 1395 | Line 1710 | void SimSetup::mpiMolDivide(void){ | |
1710 | } | |
1711 | local_SRI = local_bonds + local_bends + local_torsions; | |
1712 | ||
1713 | < | info[0].n_atoms = mpiSim->getMyNlocal(); |
1713 | > | info[0].n_atoms = mpiSim->getNAtomsLocal(); |
1714 | ||
1400 | – | |
1715 | if (local_atoms != info[0].n_atoms){ | |
1716 | sprintf(painCave.errMsg, | |
1717 | "SimSetup error: mpiSim's localAtom (%d) and SimSetup's\n" | |
# | Line 1407 | Line 1721 | void SimSetup::mpiMolDivide(void){ | |
1721 | simError(); | |
1722 | } | |
1723 | ||
1724 | + | info[0].ngroup = mpiSim->getNGroupsLocal(); |
1725 | + | if (local_groups != info[0].ngroup){ |
1726 | + | sprintf(painCave.errMsg, |
1727 | + | "SimSetup error: mpiSim's localGroups (%d) and SimSetup's\n" |
1728 | + | "\tlocalGroups (%d) are not equal.\n", |
1729 | + | info[0].ngroup, local_groups); |
1730 | + | painCave.isFatal = 1; |
1731 | + | simError(); |
1732 | + | } |
1733 | + | |
1734 | info[0].n_bonds = local_bonds; | |
1735 | info[0].n_bends = local_bends; | |
1736 | info[0].n_torsions = local_torsions; | |
# | Line 1443 | Line 1767 | void SimSetup::makeSysArrays(void){ | |
1767 | ||
1768 | ||
1769 | molIndex = 0; | |
1770 | < | for (i = 0; i < mpiSim->getTotNmol(); i++){ |
1770 | > | for (i = 0; i < mpiSim->getNMolGlobal(); i++){ |
1771 | if (mol2proc[i] == worldRank){ | |
1772 | the_molecules[molIndex].setStampID(molCompType[i]); | |
1773 | the_molecules[molIndex].setMyIndex(molIndex); | |
# | Line 1480 | Line 1804 | void SimSetup::makeSysArrays(void){ | |
1804 | info[l].atoms = the_atoms; | |
1805 | info[l].molecules = the_molecules; | |
1806 | info[l].nGlobalExcludes = 0; | |
1807 | < | |
1807 | > | |
1808 | the_ff->setSimInfo(info); | |
1809 | } | |
1810 | } |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |