# | 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 | ||
17 | #ifdef IS_MPI | |
18 | #include "mpiBASS.h" | |
# | Line 147 | Line 148 | void SimSetup::createSim(void){ | |
148 | // make the output filenames | |
149 | ||
150 | makeOutNames(); | |
150 | – | |
151 | – | if (globals->haveMinimizer()) |
152 | – | // make minimizer |
153 | – | makeMinimizer(); |
154 | – | else |
155 | – | // make the integrator |
156 | – | makeIntegrator(); |
151 | ||
152 | #ifdef IS_MPI | |
153 | mpiSim->mpiRefresh(); | |
# | Line 162 | Line 156 | void SimSetup::createSim(void){ | |
156 | // initialize the Fortran | |
157 | ||
158 | initFortran(); | |
159 | + | |
160 | + | if (globals->haveMinimizer()) |
161 | + | // make minimizer |
162 | + | makeMinimizer(); |
163 | + | else |
164 | + | // make the integrator |
165 | + | makeIntegrator(); |
166 | + | |
167 | } | |
168 | ||
169 | ||
170 | void SimSetup::makeMolecules(void){ | |
171 | int i, j, k; | |
172 | < | int exI, exJ, exK, exL, slI; |
172 | > | int exI, exJ, exK, exL, slI, slJ; |
173 | int tempI, tempJ, tempK, tempL; | |
174 | < | int molI; |
175 | < | int stampID, atomOffset, rbOffset; |
174 | > | int molI, globalID; |
175 | > | int stampID, atomOffset, rbOffset, groupOffset; |
176 | molInit molInfo; | |
177 | DirectionalAtom* dAtom; | |
178 | RigidBody* myRB; | |
# | Line 182 | Line 184 | void SimSetup::makeMolecules(void){ | |
184 | BendStamp* currentBend; | |
185 | TorsionStamp* currentTorsion; | |
186 | RigidBodyStamp* currentRigidBody; | |
187 | + | CutoffGroupStamp* currentCutoffGroup; |
188 | + | CutoffGroup* myCutoffGroup; |
189 | + | int nCutoffGroups;// number of cutoff group of a molecule defined in mdl file |
190 | + | set<int> cutoffAtomSet; //atoms belong to cutoffgroup defined at mdl file |
191 | ||
192 | bond_pair* theBonds; | |
193 | bend_set* theBends; | |
# | Line 190 | Line 196 | void SimSetup::makeMolecules(void){ | |
196 | set<int> skipList; | |
197 | ||
198 | double phi, theta, psi; | |
199 | + | char* molName; |
200 | + | char rbName[100]; |
201 | ||
202 | + | //ConstraintPair* consPair; //constraint pair |
203 | + | //ConstraintElement* consElement1; //first element of constraint pair |
204 | + | //ConstraintElement* consElement2; //second element of constraint pair |
205 | + | //int whichRigidBody; |
206 | + | //int consAtomIndex; //index of constraint atom in rigid body's atom array |
207 | + | //vector<pair<int, int> > jointAtoms; |
208 | //init the forceField paramters | |
209 | ||
210 | the_ff->readParams(); | |
# | Line 201 | Line 215 | void SimSetup::makeMolecules(void){ | |
215 | ||
216 | for (k = 0; k < nInfo; k++){ | |
217 | the_ff->setSimInfo(&(info[k])); | |
218 | + | |
219 | + | #ifdef IS_MPI |
220 | + | info[k].globalGroupMembership = new int[mpiSim->getNAtomsGlobal()]; |
221 | + | for (i = 0; i < mpiSim->getNAtomsGlobal(); i++) |
222 | + | info[k].globalGroupMembership[i] = 0; |
223 | + | #else |
224 | + | info[k].globalGroupMembership = new int[info[k].n_atoms]; |
225 | + | for (i = 0; i < info[k].n_atoms; i++) |
226 | + | info[k].globalGroupMembership[i] = 0; |
227 | + | #endif |
228 | ||
229 | atomOffset = 0; | |
230 | + | groupOffset = 0; |
231 | ||
232 | for (i = 0; i < info[k].n_mol; i++){ | |
233 | stampID = info[k].molecules[i].getStampID(); | |
234 | + | molName = comp_stamps[stampID]->getID(); |
235 | ||
236 | molInfo.nAtoms = comp_stamps[stampID]->getNAtoms(); | |
237 | molInfo.nBonds = comp_stamps[stampID]->getNBonds(); | |
238 | molInfo.nBends = comp_stamps[stampID]->getNBends(); | |
239 | molInfo.nTorsions = comp_stamps[stampID]->getNTorsions(); | |
240 | molInfo.nRigidBodies = comp_stamps[stampID]->getNRigidBodies(); | |
241 | + | |
242 | + | nCutoffGroups = comp_stamps[stampID]->getNCutoffGroups(); |
243 | ||
244 | molInfo.myAtoms = &(info[k].atoms[atomOffset]); | |
245 | ||
246 | if (molInfo.nBonds > 0) | |
247 | < | molInfo.myBonds = new (Bond *) [molInfo.nBonds]; |
247 | > | molInfo.myBonds = new Bond*[molInfo.nBonds]; |
248 | else | |
249 | molInfo.myBonds = NULL; | |
250 | ||
251 | if (molInfo.nBends > 0) | |
252 | < | molInfo.myBends = new (Bend *) [molInfo.nBends]; |
252 | > | molInfo.myBends = new Bend*[molInfo.nBends]; |
253 | else | |
254 | molInfo.myBends = NULL; | |
255 | ||
256 | if (molInfo.nTorsions > 0) | |
257 | < | molInfo.myTorsions = new (Torsion *) [molInfo.nTorsions]; |
257 | > | molInfo.myTorsions = new Torsion *[molInfo.nTorsions]; |
258 | else | |
259 | molInfo.myTorsions = NULL; | |
260 | ||
# | Line 259 | Line 287 | void SimSetup::makeMolecules(void){ | |
287 | else{ | |
288 | ||
289 | molInfo.myAtoms[j] = new Atom((j + atomOffset), info[k].getConfiguration()); | |
290 | + | |
291 | } | |
292 | ||
293 | molInfo.myAtoms[j]->setType(currentAtom->getType()); | |
265 | – | |
294 | #ifdef IS_MPI | |
295 | < | |
268 | < | molInfo.myAtoms[j]->setGlobalIndex(globalIndex[j + atomOffset]); |
269 | < | |
295 | > | molInfo.myAtoms[j]->setGlobalIndex(globalAtomIndex[j + atomOffset]); |
296 | #endif // is_mpi | |
297 | } | |
298 | ||
# | Line 406 | Line 432 | void SimSetup::makeMolecules(void){ | |
432 | info[k].excludes->addPair(exK, exL); | |
433 | } | |
434 | ||
435 | + | |
436 | + | molInfo.myRigidBodies.clear(); |
437 | + | |
438 | for (j = 0; j < molInfo.nRigidBodies; j++){ | |
439 | ||
440 | currentRigidBody = comp_stamps[stampID]->getRigidBody(j); | |
# | Line 414 | Line 443 | void SimSetup::makeMolecules(void){ | |
443 | // Create the Rigid Body: | |
444 | ||
445 | myRB = new RigidBody(); | |
446 | + | |
447 | + | sprintf(rbName,"%s_RB_%d", molName, j); |
448 | + | myRB->setType(rbName); |
449 | ||
450 | for (rb1 = 0; rb1 < nMembers; rb1++) { | |
451 | ||
# | Line 454 | Line 486 | void SimSetup::makeMolecules(void){ | |
486 | // used for the exclude list: | |
487 | ||
488 | #ifdef IS_MPI | |
489 | < | exI = info[k].atoms[tempI]->getGlobalIndex() + 1; |
490 | < | exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1; |
489 | > | exI = molInfo.myAtoms[tempI]->getGlobalIndex() + 1; |
490 | > | exJ = molInfo.myAtoms[tempJ]->getGlobalIndex() + 1; |
491 | #else | |
492 | < | exI = tempI + 1; |
493 | < | exJ = tempJ + 1; |
492 | > | exI = molInfo.myAtoms[tempI]->getIndex() + 1; |
493 | > | exJ = molInfo.myAtoms[tempJ]->getIndex() + 1; |
494 | #endif | |
495 | ||
496 | info[k].excludes->addPair(exI, exJ); | |
497 | ||
498 | } | |
499 | } | |
500 | + | |
501 | + | molInfo.myRigidBodies.push_back(myRB); |
502 | + | info[k].rigidBodies.push_back(myRB); |
503 | } | |
504 | ||
505 | + | |
506 | + | //create cutoff group for molecule |
507 | + | |
508 | + | cutoffAtomSet.clear(); |
509 | + | molInfo.myCutoffGroups.clear(); |
510 | + | |
511 | + | for (j = 0; j < nCutoffGroups; j++){ |
512 | + | |
513 | + | currentCutoffGroup = comp_stamps[stampID]->getCutoffGroup(j); |
514 | + | nMembers = currentCutoffGroup->getNMembers(); |
515 | + | |
516 | + | myCutoffGroup = new CutoffGroup(); |
517 | + | |
518 | + | #ifdef IS_MPI |
519 | + | myCutoffGroup->setGlobalIndex(globalGroupIndex[groupOffset]); |
520 | + | #else |
521 | + | myCutoffGroup->setGlobalIndex(groupOffset); |
522 | + | #endif |
523 | + | |
524 | + | for (int cg = 0; cg < nMembers; cg++) { |
525 | + | |
526 | + | // molI is atom numbering inside this molecule |
527 | + | molI = currentCutoffGroup->getMember(cg); |
528 | + | |
529 | + | // tempI is atom numbering on local processor |
530 | + | tempI = molI + atomOffset; |
531 | + | |
532 | + | #ifdef IS_MPI |
533 | + | globalID = info[k].atoms[tempI]->getGlobalIndex(); |
534 | + | info[k].globalGroupMembership[globalID] = globalGroupIndex[groupOffset]; |
535 | + | #else |
536 | + | globalID = info[k].atoms[tempI]->getIndex(); |
537 | + | info[k].globalGroupMembership[globalID] = groupOffset; |
538 | + | #endif |
539 | + | myCutoffGroup->addAtom(info[k].atoms[tempI]); |
540 | + | cutoffAtomSet.insert(tempI); |
541 | + | } |
542 | + | |
543 | + | molInfo.myCutoffGroups.push_back(myCutoffGroup); |
544 | + | groupOffset++; |
545 | + | |
546 | + | }//end for (j = 0; j < molInfo.nCutoffGroups; j++) |
547 | + | |
548 | + | |
549 | + | // create a cutoff group for every atom in current molecule which |
550 | + | // does not belong to cutoffgroup defined at mdl file |
551 | + | |
552 | + | for(j = 0; j < molInfo.nAtoms; j++){ |
553 | + | |
554 | + | if(cutoffAtomSet.find(molInfo.myAtoms[j]->getIndex()) == cutoffAtomSet.end()){ |
555 | + | myCutoffGroup = new CutoffGroup(); |
556 | + | myCutoffGroup->addAtom(molInfo.myAtoms[j]); |
557 | + | |
558 | + | #ifdef IS_MPI |
559 | + | myCutoffGroup->setGlobalIndex(globalGroupIndex[groupOffset]); |
560 | + | globalID = info[k].atoms[atomOffset + j]->getGlobalIndex(); |
561 | + | info[k].globalGroupMembership[globalID] = globalGroupIndex[groupOffset]; |
562 | + | #else |
563 | + | myCutoffGroup->setGlobalIndex(groupOffset); |
564 | + | globalID = info[k].atoms[atomOffset + j]->getIndex(); |
565 | + | info[k].globalGroupMembership[globalID] = groupOffset; |
566 | + | #endif |
567 | + | molInfo.myCutoffGroups.push_back(myCutoffGroup); |
568 | + | groupOffset++; |
569 | + | } |
570 | + | } |
571 | + | |
572 | + | // After this is all set up, scan through the atoms to |
573 | + | // see if they can be added to the integrableObjects: |
574 | + | |
575 | + | molInfo.myIntegrableObjects.clear(); |
576 | + | |
577 | + | |
578 | + | for (j = 0; j < molInfo.nAtoms; j++){ |
579 | + | |
580 | + | #ifdef IS_MPI |
581 | + | slJ = molInfo.myAtoms[j]->getGlobalIndex(); |
582 | + | #else |
583 | + | slJ = j+atomOffset; |
584 | + | #endif |
585 | + | |
586 | + | // if they aren't on the skip list, then they can be integrated |
587 | + | |
588 | + | if (skipList.find(slJ) == skipList.end()) { |
589 | + | mySD = (StuntDouble *) molInfo.myAtoms[j]; |
590 | + | info[k].integrableObjects.push_back(mySD); |
591 | + | molInfo.myIntegrableObjects.push_back(mySD); |
592 | + | } |
593 | + | } |
594 | + | |
595 | + | // all rigid bodies are integrated: |
596 | + | |
597 | + | for (j = 0; j < molInfo.nRigidBodies; j++) { |
598 | + | mySD = (StuntDouble *) molInfo.myRigidBodies[j]; |
599 | + | info[k].integrableObjects.push_back(mySD); |
600 | + | molInfo.myIntegrableObjects.push_back(mySD); |
601 | + | } |
602 | + | |
603 | + | |
604 | + | /* |
605 | + | |
606 | + | //creat ConstraintPair. |
607 | + | molInfo.myConstraintPair.clear(); |
608 | + | |
609 | + | for (j = 0; j < molInfo.nBonds; j++){ |
610 | + | |
611 | + | //if both atoms are in the same rigid body, just skip it |
612 | + | currentBond = comp_stamps[stampID]->getBond(j); |
613 | + | if(!comp_stamps[stampID]->isBondInSameRigidBody(currentBond)){ |
614 | + | |
615 | + | tempI = currentBond->getA() + atomOffset; |
616 | + | if( comp_stamps[stampID]->isAtomInRigidBody(currentBond->getA(), whichRigidBody, consAtomIndex)) |
617 | + | consElement1 = new ConstraintRigidBody(molInfo.myRigidBodies[whichRigidBody], consAtomIndex); |
618 | + | else |
619 | + | consElement1 = new ConstraintAtom(info[k].atoms[tempI]); |
620 | + | |
621 | + | tempJ = currentBond->getB() + atomOffset; |
622 | + | if(comp_stamps[stampID]->isAtomInRigidBody(currentBond->getB(), whichRigidBody, consAtomIndex)) |
623 | + | consElement2 = new ConstraintRigidBody(molInfo.myRigidBodies[whichRigidBody], consAtomIndex); |
624 | + | else |
625 | + | consElement2 = new ConstraintAtom(info[k].atoms[tempJ]); |
626 | + | |
627 | + | consPair = new DistanceConstraintPair(consElement1, consElement2); |
628 | + | molInfo.myConstraintPairs.push_back(consPair); |
629 | + | } |
630 | + | } |
631 | + | |
632 | + | //loop over rigid bodies, if two rigid bodies share same joint, creat a HingeConstraintPair |
633 | + | for (int rb1 = 0; rb1 < molInfo.nRigidBodies -1 ; rb1++){ |
634 | + | for (int rb2 = rb1 + 1; rb2 < molInfo.nRigidBodies ; rb2++){ |
635 | + | |
636 | + | jointAtoms = comp_stamps[stampID]->getJointAtoms(rb1, rb2); |
637 | + | |
638 | + | for(size_t m = 0; m < jointAtoms.size(); m++){ |
639 | + | consElement1 = new ConstraintRigidBody(molInfo.myRigidBodies[rb1], jointAtoms[m].first); |
640 | + | consElement2 = new ConstraintRigidBody(molInfo.myRigidBodies[rb2], jointAtoms[m].second); |
641 | + | |
642 | + | consPair = new JointConstraintPair(consElement1, consElement2); |
643 | + | molInfo.myConstraintPairs.push_back(consPair); |
644 | + | } |
645 | + | |
646 | + | } |
647 | + | } |
648 | + | |
649 | + | */ |
650 | // send the arrays off to the forceField for init. | |
651 | ||
652 | the_ff->initializeAtoms(molInfo.nAtoms, molInfo.myAtoms); | |
# | Line 476 | Line 656 | void SimSetup::makeMolecules(void){ | |
656 | theTorsions); | |
657 | ||
658 | info[k].molecules[i].initialize(molInfo); | |
659 | < | |
660 | < | |
659 | > | |
660 | > | |
661 | atomOffset += molInfo.nAtoms; | |
662 | delete[] theBonds; | |
663 | delete[] theBends; | |
664 | delete[] theTorsions; | |
665 | } | |
666 | ||
487 | – | // build up the integrableObjects vector: |
667 | ||
489 | – | for (i = 0; i < info[k].n_atoms; i++) { |
490 | – | |
491 | – | #ifdef IS_MPI |
492 | – | slI = info[k].atoms[i]->getGlobalIndex(); |
493 | – | #else |
494 | – | slI = i; |
495 | – | #endif |
668 | ||
669 | < | if (skipList.find(slI) == skipList.end()) { |
670 | < | mySD = (StuntDouble *) info[k].atoms[i]; |
671 | < | info[k].integrableObjects.push_back(mySD); |
672 | < | } |
673 | < | } |
674 | < | for (i = 0; i < info[k].rigidBodies.size(); i++) { |
675 | < | mySD = (StuntDouble *) info[k].rigidBodies[i]; |
676 | < | info[k].integrableObjects.push_back(mySD); |
677 | < | } |
669 | > | #ifdef IS_MPI |
670 | > | // Since the globalGroupMembership has been zero filled and we've only |
671 | > | // poked values into the atoms we know, we can do an Allreduce |
672 | > | // to get the full globalGroupMembership array (We think). |
673 | > | // This would be prettier if we could use MPI_IN_PLACE like the MPI-2 |
674 | > | // docs said we could. |
675 | > | |
676 | > | int* ggMjunk = new int[mpiSim->getNAtomsGlobal()]; |
677 | > | |
678 | > | MPI_Allreduce(info[k].globalGroupMembership, |
679 | > | ggMjunk, |
680 | > | mpiSim->getNAtomsGlobal(), |
681 | > | MPI_INT, MPI_SUM, MPI_COMM_WORLD); |
682 | > | |
683 | > | for (i = 0; i < mpiSim->getNAtomsGlobal(); i++) |
684 | > | info[k].globalGroupMembership[i] = ggMjunk[i]; |
685 | > | |
686 | > | delete[] ggMjunk; |
687 | ||
688 | + | #endif |
689 | + | |
690 | + | |
691 | + | |
692 | } | |
693 | ||
694 | #ifdef IS_MPI | |
695 | sprintf(checkPointMsg, "all molecules initialized succesfully"); | |
696 | MPIcheckPoint(); | |
697 | #endif // is_mpi | |
513 | – | |
514 | – | // clean up the forcefield |
515 | – | |
516 | – | if (!globals->haveLJrcut()){ |
698 | ||
518 | – | the_ff->calcRcut(); |
519 | – | |
520 | – | } else { |
521 | – | |
522 | – | the_ff->setRcut( globals->getLJrcut() ); |
523 | – | } |
524 | – | |
525 | – | the_ff->cleanMe(); |
699 | } | |
700 | ||
701 | void SimSetup::initFromBass(void){ | |
# | Line 809 | Line 982 | void SimSetup::gatherInfo(void){ | |
982 | } | |
983 | ||
984 | //check whether sample time, status time, thermal time and reset time are divisble by dt | |
985 | < | if (!isDivisible(globals->getSampleTime(), globals->getDt())){ |
985 | > | if (globals->haveSampleTime() && !isDivisible(globals->getSampleTime(), globals->getDt())){ |
986 | sprintf(painCave.errMsg, | |
987 | "Sample time is not divisible by dt.\n" | |
988 | "\tThis will result in samples that are not uniformly\n" | |
# | Line 819 | Line 992 | void SimSetup::gatherInfo(void){ | |
992 | simError(); | |
993 | } | |
994 | ||
995 | < | if (globals->haveStatusTime() && !isDivisible(globals->getSampleTime(), globals->getDt())){ |
995 | > | if (globals->haveStatusTime() && !isDivisible(globals->getStatusTime(), globals->getDt())){ |
996 | sprintf(painCave.errMsg, | |
997 | "Status time is not divisible by dt.\n" | |
998 | "\tThis will result in status reports that are not uniformly\n" | |
# | Line 855 | Line 1028 | void SimSetup::gatherInfo(void){ | |
1028 | if (globals->haveSampleTime()){ | |
1029 | info[i].sampleTime = globals->getSampleTime(); | |
1030 | info[i].statusTime = info[i].sampleTime; | |
858 | – | info[i].thermalTime = info[i].sampleTime; |
1031 | } | |
1032 | else{ | |
1033 | info[i].sampleTime = globals->getRunTime(); | |
1034 | info[i].statusTime = info[i].sampleTime; | |
863 | – | info[i].thermalTime = info[i].sampleTime; |
1035 | } | |
1036 | ||
1037 | if (globals->haveStatusTime()){ | |
# | Line 869 | Line 1040 | void SimSetup::gatherInfo(void){ | |
1040 | ||
1041 | if (globals->haveThermalTime()){ | |
1042 | info[i].thermalTime = globals->getThermalTime(); | |
1043 | + | } else { |
1044 | + | info[i].thermalTime = globals->getRunTime(); |
1045 | } | |
1046 | ||
1047 | info[i].resetIntegrator = 0; | |
# | Line 886 | Line 1059 | void SimSetup::gatherInfo(void){ | |
1059 | ||
1060 | info[i].useInitXSstate = globals->getUseInitXSstate(); | |
1061 | info[i].orthoTolerance = globals->getOrthoBoxTolerance(); | |
1062 | < | |
1062 | > | |
1063 | > | // check for thermodynamic integration |
1064 | > | if (globals->getUseSolidThermInt() && !globals->getUseLiquidThermInt()) { |
1065 | > | if (globals->haveThermIntLambda() && globals->haveThermIntK()) { |
1066 | > | info[i].useSolidThermInt = globals->getUseSolidThermInt(); |
1067 | > | info[i].thermIntLambda = globals->getThermIntLambda(); |
1068 | > | info[i].thermIntK = globals->getThermIntK(); |
1069 | > | |
1070 | > | Restraints *myRestraint = new Restraints(tot_nmol, info[i].thermIntLambda, info[i].thermIntK); |
1071 | > | info[i].restraint = myRestraint; |
1072 | > | } |
1073 | > | else { |
1074 | > | sprintf(painCave.errMsg, |
1075 | > | "SimSetup Error:\n" |
1076 | > | "\tKeyword useSolidThermInt was set to 'true' but\n" |
1077 | > | "\tthermodynamicIntegrationLambda (and/or\n" |
1078 | > | "\tthermodynamicIntegrationK) was not specified.\n" |
1079 | > | "\tPlease provide a lambda value and k value in your .bass file.\n"); |
1080 | > | painCave.isFatal = 1; |
1081 | > | simError(); |
1082 | > | } |
1083 | > | } |
1084 | > | else if(globals->getUseLiquidThermInt()) { |
1085 | > | if (globals->getUseSolidThermInt()) { |
1086 | > | sprintf( painCave.errMsg, |
1087 | > | "SimSetup Warning: It appears that you have both solid and\n" |
1088 | > | "\tliquid thermodynamic integration activated in your .bass\n" |
1089 | > | "\tfile. To avoid confusion, specify only one technique in\n" |
1090 | > | "\tyour .bass file. Liquid-state thermodynamic integration\n" |
1091 | > | "\twill be assumed for the current simulation. If this is not\n" |
1092 | > | "\twhat you desire, set useSolidThermInt to 'true' and\n" |
1093 | > | "\tuseLiquidThermInt to 'false' in your .bass file.\n"); |
1094 | > | painCave.isFatal = 0; |
1095 | > | simError(); |
1096 | > | } |
1097 | > | if (globals->haveThermIntLambda() && globals->haveThermIntK()) { |
1098 | > | info[i].useLiquidThermInt = globals->getUseLiquidThermInt(); |
1099 | > | info[i].thermIntLambda = globals->getThermIntLambda(); |
1100 | > | info[i].thermIntK = globals->getThermIntK(); |
1101 | > | } |
1102 | > | else { |
1103 | > | sprintf(painCave.errMsg, |
1104 | > | "SimSetup Error:\n" |
1105 | > | "\tKeyword useLiquidThermInt was set to 'true' but\n" |
1106 | > | "\tthermodynamicIntegrationLambda (and/or\n" |
1107 | > | "\tthermodynamicIntegrationK) was not specified.\n" |
1108 | > | "\tPlease provide a lambda value and k value in your .bass file.\n"); |
1109 | > | painCave.isFatal = 1; |
1110 | > | simError(); |
1111 | > | } |
1112 | > | } |
1113 | > | else if(globals->haveThermIntLambda() || globals->haveThermIntK()){ |
1114 | > | sprintf(painCave.errMsg, |
1115 | > | "SimSetup Warning: If you want to use Thermodynamic\n" |
1116 | > | "\tIntegration, set useSolidThermInt or useLiquidThermInt to\n" |
1117 | > | "\t'true' in your .bass file. These keywords are set to\n" |
1118 | > | "\t'false' by default, so your lambda and/or k values are\n" |
1119 | > | "\tbeing ignored.\n"); |
1120 | > | painCave.isFatal = 0; |
1121 | > | simError(); |
1122 | > | } |
1123 | } | |
1124 | ||
1125 | //setup seed for random number generator | |
# | Line 939 | Line 1172 | void SimSetup::finalInfoCheck(void){ | |
1172 | void SimSetup::finalInfoCheck(void){ | |
1173 | int index; | |
1174 | int usesDipoles; | |
1175 | + | int usesCharges; |
1176 | int i; | |
1177 | ||
1178 | for (i = 0; i < nInfo; i++){ | |
# | Line 950 | Line 1184 | void SimSetup::finalInfoCheck(void){ | |
1184 | usesDipoles = (info[i].atoms[index])->hasDipole(); | |
1185 | index++; | |
1186 | } | |
1187 | < | |
1187 | > | index = 0; |
1188 | > | usesCharges = 0; |
1189 | > | while ((index < info[i].n_atoms) && !usesCharges){ |
1190 | > | usesCharges= (info[i].atoms[index])->hasCharge(); |
1191 | > | index++; |
1192 | > | } |
1193 | #ifdef IS_MPI | |
1194 | int myUse = usesDipoles; | |
1195 | MPI_Allreduce(&myUse, &usesDipoles, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | |
1196 | #endif //is_mpi | |
1197 | ||
1198 | < | double theEcr, theEst; |
1198 | > | double theRcut, theRsw; |
1199 | ||
1200 | + | if (globals->haveRcut()) { |
1201 | + | theRcut = globals->getRcut(); |
1202 | + | |
1203 | + | if (globals->haveRsw()) |
1204 | + | theRsw = globals->getRsw(); |
1205 | + | else |
1206 | + | theRsw = theRcut; |
1207 | + | |
1208 | + | info[i].setDefaultRcut(theRcut, theRsw); |
1209 | + | |
1210 | + | } else { |
1211 | + | |
1212 | + | the_ff->calcRcut(); |
1213 | + | theRcut = info[i].getRcut(); |
1214 | + | |
1215 | + | if (globals->haveRsw()) |
1216 | + | theRsw = globals->getRsw(); |
1217 | + | else |
1218 | + | theRsw = theRcut; |
1219 | + | |
1220 | + | info[i].setDefaultRcut(theRcut, theRsw); |
1221 | + | } |
1222 | + | |
1223 | if (globals->getUseRF()){ | |
1224 | info[i].useReactionField = 1; | |
1225 | < | |
1226 | < | if (!globals->haveECR()){ |
1225 | > | |
1226 | > | if (!globals->haveRcut()){ |
1227 | sprintf(painCave.errMsg, | |
1228 | < | "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n" |
1228 | > | "SimSetup Warning: No value was set for the cutoffRadius.\n" |
1229 | "\tOOPSE will use a default value of 15.0 angstroms" | |
1230 | < | "\tfor the electrostaticCutoffRadius.\n"); |
1230 | > | "\tfor the cutoffRadius.\n"); |
1231 | painCave.isFatal = 0; | |
1232 | simError(); | |
1233 | < | theEcr = 15.0; |
1233 | > | theRcut = 15.0; |
1234 | } | |
1235 | else{ | |
1236 | < | theEcr = globals->getECR(); |
1236 | > | theRcut = globals->getRcut(); |
1237 | } | |
1238 | ||
1239 | < | if (!globals->haveEST()){ |
1239 | > | if (!globals->haveRsw()){ |
1240 | sprintf(painCave.errMsg, | |
1241 | < | "SimSetup Warning: No value was set for electrostaticSkinThickness.\n" |
1241 | > | "SimSetup Warning: No value was set for switchingRadius.\n" |
1242 | "\tOOPSE will use a default value of\n" | |
1243 | < | "\t0.05 * electrostaticCutoffRadius\n" |
982 | < | "\tfor the electrostaticSkinThickness\n"); |
1243 | > | "\t0.95 * cutoffRadius for the switchingRadius\n"); |
1244 | painCave.isFatal = 0; | |
1245 | simError(); | |
1246 | < | theEst = 0.05 * theEcr; |
1246 | > | theRsw = 0.95 * theRcut; |
1247 | } | |
1248 | else{ | |
1249 | < | theEst = globals->getEST(); |
1249 | > | theRsw = globals->getRsw(); |
1250 | } | |
1251 | ||
1252 | < | info[i].setDefaultEcr(theEcr, theEst); |
1252 | > | info[i].setDefaultRcut(theRcut, theRsw); |
1253 | ||
1254 | if (!globals->haveDielectric()){ | |
1255 | sprintf(painCave.errMsg, | |
# | Line 1001 | Line 1262 | void SimSetup::finalInfoCheck(void){ | |
1262 | info[i].dielectric = globals->getDielectric(); | |
1263 | } | |
1264 | else{ | |
1265 | < | if (usesDipoles){ |
1266 | < | if (!globals->haveECR()){ |
1265 | > | if (usesDipoles || usesCharges){ |
1266 | > | |
1267 | > | if (!globals->haveRcut()){ |
1268 | sprintf(painCave.errMsg, | |
1269 | < | "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n" |
1269 | > | "SimSetup Warning: No value was set for the cutoffRadius.\n" |
1270 | "\tOOPSE will use a default value of 15.0 angstroms" | |
1271 | < | "\tfor the electrostaticCutoffRadius.\n"); |
1272 | < | painCave.isFatal = 0; |
1273 | < | simError(); |
1274 | < | theEcr = 15.0; |
1275 | < | } |
1276 | < | else{ |
1277 | < | theEcr = globals->getECR(); |
1271 | > | "\tfor the cutoffRadius.\n"); |
1272 | > | painCave.isFatal = 0; |
1273 | > | simError(); |
1274 | > | theRcut = 15.0; |
1275 | > | } |
1276 | > | else{ |
1277 | > | theRcut = globals->getRcut(); |
1278 | } | |
1279 | < | |
1280 | < | if (!globals->haveEST()){ |
1279 | > | |
1280 | > | if (!globals->haveRsw()){ |
1281 | sprintf(painCave.errMsg, | |
1282 | < | "SimSetup Warning: No value was set for electrostaticSkinThickness.\n" |
1282 | > | "SimSetup Warning: No value was set for switchingRadius.\n" |
1283 | "\tOOPSE will use a default value of\n" | |
1284 | < | "\t0.05 * electrostaticCutoffRadius\n" |
1023 | < | "\tfor the electrostaticSkinThickness\n"); |
1284 | > | "\t0.95 * cutoffRadius for the switchingRadius\n"); |
1285 | painCave.isFatal = 0; | |
1286 | simError(); | |
1287 | < | theEst = 0.05 * theEcr; |
1287 | > | theRsw = 0.95 * theRcut; |
1288 | } | |
1289 | else{ | |
1290 | < | theEst = globals->getEST(); |
1290 | > | theRsw = globals->getRsw(); |
1291 | } | |
1292 | + | |
1293 | + | info[i].setDefaultRcut(theRcut, theRsw); |
1294 | ||
1032 | – | info[i].setDefaultEcr(theEcr, theEst); |
1295 | } | |
1296 | } | |
1297 | } | |
# | Line 1037 | Line 1299 | void SimSetup::finalInfoCheck(void){ | |
1299 | strcpy(checkPointMsg, "post processing checks out"); | |
1300 | MPIcheckPoint(); | |
1301 | #endif // is_mpi | |
1302 | + | |
1303 | + | // clean up the forcefield |
1304 | + | the_ff->cleanMe(); |
1305 | } | |
1306 | ||
1307 | void SimSetup::initSystemCoords(void){ | |
# | Line 1167 | Line 1432 | void SimSetup::makeOutNames(void){ | |
1432 | } | |
1433 | } | |
1434 | ||
1435 | + | strcpy(info[k].rawPotName, inFileName); |
1436 | + | nameLength = strlen(info[k].rawPotName); |
1437 | + | endTest = &(info[k].rawPotName[nameLength - 5]); |
1438 | + | if (!strcmp(endTest, ".bass")){ |
1439 | + | strcpy(endTest, ".raw"); |
1440 | + | } |
1441 | + | else if (!strcmp(endTest, ".BASS")){ |
1442 | + | strcpy(endTest, ".raw"); |
1443 | + | } |
1444 | + | else{ |
1445 | + | endTest = &(info[k].rawPotName[nameLength - 4]); |
1446 | + | if (!strcmp(endTest, ".bss")){ |
1447 | + | strcpy(endTest, ".raw"); |
1448 | + | } |
1449 | + | else if (!strcmp(endTest, ".mdl")){ |
1450 | + | strcpy(endTest, ".raw"); |
1451 | + | } |
1452 | + | else{ |
1453 | + | strcat(info[k].rawPotName, ".raw"); |
1454 | + | } |
1455 | + | } |
1456 | + | |
1457 | #ifdef IS_MPI | |
1458 | ||
1459 | } | |
# | Line 1251 | Line 1538 | void SimSetup::compList(void){ | |
1538 | LinkedMolStamp* headStamp = new LinkedMolStamp(); | |
1539 | LinkedMolStamp* currentStamp = NULL; | |
1540 | comp_stamps = new MoleculeStamp * [n_components]; | |
1541 | + | bool haveCutoffGroups; |
1542 | ||
1543 | + | haveCutoffGroups = false; |
1544 | + | |
1545 | // make an array of molecule stamps that match the components used. | |
1546 | // also extract the used stamps out into a separate linked list | |
1547 | ||
# | Line 1286 | Line 1576 | void SimSetup::compList(void){ | |
1576 | headStamp->add(currentStamp); | |
1577 | comp_stamps[i] = headStamp->match(id); | |
1578 | } | |
1579 | + | |
1580 | + | if(comp_stamps[i]->getNCutoffGroups() > 0) |
1581 | + | haveCutoffGroups = true; |
1582 | } | |
1583 | + | |
1584 | + | for (i = 0; i < nInfo; i++) |
1585 | + | info[i].haveCutoffGroups = haveCutoffGroups; |
1586 | ||
1587 | #ifdef IS_MPI | |
1588 | strcpy(checkPointMsg, "Component stamps successfully extracted\n"); | |
# | Line 1295 | Line 1591 | void SimSetup::calcSysValues(void){ | |
1591 | } | |
1592 | ||
1593 | void SimSetup::calcSysValues(void){ | |
1594 | < | int i; |
1594 | > | int i, j; |
1595 | > | int ncutgroups, atomsingroups, ngroupsinstamp; |
1596 | ||
1597 | int* molMembershipArray; | |
1598 | + | CutoffGroupStamp* cg; |
1599 | ||
1600 | tot_atoms = 0; | |
1601 | tot_bonds = 0; | |
1602 | tot_bends = 0; | |
1603 | tot_torsions = 0; | |
1604 | tot_rigid = 0; | |
1605 | + | tot_groups = 0; |
1606 | for (i = 0; i < n_components; i++){ | |
1607 | tot_atoms += components_nmol[i] * comp_stamps[i]->getNAtoms(); | |
1608 | tot_bonds += components_nmol[i] * comp_stamps[i]->getNBonds(); | |
1609 | tot_bends += components_nmol[i] * comp_stamps[i]->getNBends(); | |
1610 | tot_torsions += components_nmol[i] * comp_stamps[i]->getNTorsions(); | |
1611 | tot_rigid += components_nmol[i] * comp_stamps[i]->getNRigidBodies(); | |
1612 | + | |
1613 | + | ncutgroups = comp_stamps[i]->getNCutoffGroups(); |
1614 | + | atomsingroups = 0; |
1615 | + | for (j=0; j < ncutgroups; j++) { |
1616 | + | cg = comp_stamps[i]->getCutoffGroup(j); |
1617 | + | atomsingroups += cg->getNMembers(); |
1618 | + | } |
1619 | + | ngroupsinstamp = comp_stamps[i]->getNAtoms() - atomsingroups + ncutgroups; |
1620 | + | tot_groups += components_nmol[i] * ngroupsinstamp; |
1621 | } | |
1622 | ||
1623 | tot_SRI = tot_bonds + tot_bends + tot_torsions; | |
# | Line 1322 | Line 1630 | void SimSetup::calcSysValues(void){ | |
1630 | info[i].n_torsions = tot_torsions; | |
1631 | info[i].n_SRI = tot_SRI; | |
1632 | info[i].n_mol = tot_nmol; | |
1633 | < | |
1633 | > | info[i].ngroup = tot_groups; |
1634 | info[i].molMembershipArray = molMembershipArray; | |
1635 | } | |
1636 | } | |
# | Line 1333 | Line 1641 | void SimSetup::mpiMolDivide(void){ | |
1641 | int i, j, k; | |
1642 | int localMol, allMol; | |
1643 | int local_atoms, local_bonds, local_bends, local_torsions, local_SRI; | |
1644 | < | int local_rigid; |
1644 | > | int local_rigid, local_groups; |
1645 | > | vector<int> globalMolIndex; |
1646 | > | int ncutgroups, atomsingroups, ngroupsinstamp; |
1647 | > | CutoffGroupStamp* cg; |
1648 | ||
1649 | mpiSim = new mpiSimulation(info); | |
1650 | ||
1651 | < | globalIndex = mpiSim->divideLabor(); |
1651 | > | mpiSim->divideLabor(); |
1652 | > | globalAtomIndex = mpiSim->getGlobalAtomIndex(); |
1653 | > | globalGroupIndex = mpiSim->getGlobalGroupIndex(); |
1654 | > | //globalMolIndex = mpiSim->getGlobalMolIndex(); |
1655 | ||
1656 | // set up the local variables | |
1657 | ||
# | Line 1351 | Line 1665 | void SimSetup::mpiMolDivide(void){ | |
1665 | local_bends = 0; | |
1666 | local_torsions = 0; | |
1667 | local_rigid = 0; | |
1668 | < | globalAtomIndex = 0; |
1668 | > | local_groups = 0; |
1669 | > | globalAtomCounter = 0; |
1670 | ||
1671 | for (i = 0; i < n_components; i++){ | |
1672 | for (j = 0; j < components_nmol[i]; j++){ | |
# | Line 1361 | Line 1676 | void SimSetup::mpiMolDivide(void){ | |
1676 | local_bends += comp_stamps[i]->getNBends(); | |
1677 | local_torsions += comp_stamps[i]->getNTorsions(); | |
1678 | local_rigid += comp_stamps[i]->getNRigidBodies(); | |
1679 | + | |
1680 | + | ncutgroups = comp_stamps[i]->getNCutoffGroups(); |
1681 | + | atomsingroups = 0; |
1682 | + | for (k=0; k < ncutgroups; k++) { |
1683 | + | cg = comp_stamps[i]->getCutoffGroup(k); |
1684 | + | atomsingroups += cg->getNMembers(); |
1685 | + | } |
1686 | + | ngroupsinstamp = comp_stamps[i]->getNAtoms() - atomsingroups + |
1687 | + | ncutgroups; |
1688 | + | local_groups += ngroupsinstamp; |
1689 | + | |
1690 | localMol++; | |
1691 | } | |
1692 | for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){ | |
1693 | < | info[0].molMembershipArray[globalAtomIndex] = allMol; |
1694 | < | globalAtomIndex++; |
1693 | > | info[0].molMembershipArray[globalAtomCounter] = allMol; |
1694 | > | globalAtomCounter++; |
1695 | } | |
1696 | ||
1697 | allMol++; | |
# | Line 1373 | Line 1699 | void SimSetup::mpiMolDivide(void){ | |
1699 | } | |
1700 | local_SRI = local_bonds + local_bends + local_torsions; | |
1701 | ||
1702 | < | info[0].n_atoms = mpiSim->getMyNlocal(); |
1702 | > | info[0].n_atoms = mpiSim->getNAtomsLocal(); |
1703 | ||
1378 | – | |
1704 | if (local_atoms != info[0].n_atoms){ | |
1705 | sprintf(painCave.errMsg, | |
1706 | "SimSetup error: mpiSim's localAtom (%d) and SimSetup's\n" | |
# | Line 1385 | Line 1710 | void SimSetup::mpiMolDivide(void){ | |
1710 | simError(); | |
1711 | } | |
1712 | ||
1713 | + | info[0].ngroup = mpiSim->getNGroupsLocal(); |
1714 | + | if (local_groups != info[0].ngroup){ |
1715 | + | sprintf(painCave.errMsg, |
1716 | + | "SimSetup error: mpiSim's localGroups (%d) and SimSetup's\n" |
1717 | + | "\tlocalGroups (%d) are not equal.\n", |
1718 | + | info[0].ngroup, local_groups); |
1719 | + | painCave.isFatal = 1; |
1720 | + | simError(); |
1721 | + | } |
1722 | + | |
1723 | info[0].n_bonds = local_bonds; | |
1724 | info[0].n_bends = local_bends; | |
1725 | info[0].n_torsions = local_torsions; | |
# | Line 1421 | Line 1756 | void SimSetup::makeSysArrays(void){ | |
1756 | ||
1757 | ||
1758 | molIndex = 0; | |
1759 | < | for (i = 0; i < mpiSim->getTotNmol(); i++){ |
1759 | > | for (i = 0; i < mpiSim->getNMolGlobal(); i++){ |
1760 | if (mol2proc[i] == worldRank){ | |
1761 | the_molecules[molIndex].setStampID(molCompType[i]); | |
1762 | the_molecules[molIndex].setMyIndex(molIndex); | |
# | Line 1433 | Line 1768 | void SimSetup::makeSysArrays(void){ | |
1768 | #else // is_mpi | |
1769 | ||
1770 | molIndex = 0; | |
1771 | < | globalAtomIndex = 0; |
1771 | > | globalAtomCounter = 0; |
1772 | for (i = 0; i < n_components; i++){ | |
1773 | for (j = 0; j < components_nmol[i]; j++){ | |
1774 | the_molecules[molIndex].setStampID(i); | |
1775 | the_molecules[molIndex].setMyIndex(molIndex); | |
1776 | the_molecules[molIndex].setGlobalIndex(molIndex); | |
1777 | for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){ | |
1778 | < | info[l].molMembershipArray[globalAtomIndex] = molIndex; |
1779 | < | globalAtomIndex++; |
1778 | > | info[l].molMembershipArray[globalAtomCounter] = molIndex; |
1779 | > | globalAtomCounter++; |
1780 | } | |
1781 | molIndex++; | |
1782 | } | |
# | Line 1458 | Line 1793 | void SimSetup::makeSysArrays(void){ | |
1793 | info[l].atoms = the_atoms; | |
1794 | info[l].molecules = the_molecules; | |
1795 | info[l].nGlobalExcludes = 0; | |
1796 | < | |
1796 | > | |
1797 | the_ff->setSimInfo(info); | |
1798 | } | |
1799 | } |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |