# | Line 30 | Line 30 | using namespace std; | |
---|---|---|
30 | ||
31 | using namespace std; | |
32 | ||
33 | + | /** |
34 | + | * Check whether dividend is divisble by divisor or not |
35 | + | */ |
36 | + | bool isDivisible(double dividend, double divisor){ |
37 | + | double tolerance = 0.000001; |
38 | + | double quotient; |
39 | + | double diff; |
40 | + | int intQuotient; |
41 | + | |
42 | + | quotient = dividend / divisor; |
43 | + | |
44 | + | if (quotient < 0) |
45 | + | quotient = -quotient; |
46 | + | |
47 | + | intQuotient = int (quotient + tolerance); |
48 | + | |
49 | + | diff = fabs(fabs(dividend) - intQuotient * fabs(divisor)); |
50 | + | |
51 | + | if (diff <= tolerance) |
52 | + | return true; |
53 | + | else |
54 | + | return false; |
55 | + | } |
56 | + | |
57 | SimSetup::SimSetup(){ | |
58 | ||
59 | initSuspend = false; | |
# | Line 150 | Line 174 | void SimSetup::makeMolecules(void){ | |
174 | bend_set* theBends; | |
175 | torsion_set* theTorsions; | |
176 | ||
153 | – | |
177 | //init the forceField paramters | |
178 | ||
179 | the_ff->readParams(); | |
# | Line 158 | Line 181 | void SimSetup::makeMolecules(void){ | |
181 | ||
182 | // init the atoms | |
183 | ||
184 | + | double phi, theta, psi; |
185 | + | double sux, suy, suz; |
186 | + | double Axx, Axy, Axz, Ayx, Ayy, Ayz, Azx, Azy, Azz; |
187 | double ux, uy, uz, u, uSqr; | |
188 | ||
189 | for (k = 0; k < nInfo; k++){ | |
# | Line 194 | Line 220 | void SimSetup::makeMolecules(void){ | |
220 | info[k].n_oriented++; | |
221 | molInfo.myAtoms[j] = dAtom; | |
222 | ||
223 | < | ux = currentAtom->getOrntX(); |
224 | < | uy = currentAtom->getOrntY(); |
225 | < | uz = currentAtom->getOrntZ(); |
223 | > | // Directional Atoms have standard unit vectors which are oriented |
224 | > | // in space using the three Euler angles. We assume the standard |
225 | > | // unit vector was originally along the z axis below. |
226 | ||
227 | + | phi = currentAtom->getEulerPhi() * M_PI / 180.0; |
228 | + | theta = currentAtom->getEulerTheta() * M_PI / 180.0; |
229 | + | psi = currentAtom->getEulerPsi()* M_PI / 180.0; |
230 | + | |
231 | + | Axx = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi)); |
232 | + | Axy = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi)); |
233 | + | Axz = sin(theta) * sin(psi); |
234 | + | |
235 | + | Ayx = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi)); |
236 | + | Ayy = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi)); |
237 | + | Ayz = sin(theta) * cos(psi); |
238 | + | |
239 | + | Azx = sin(phi) * sin(theta); |
240 | + | Azy = -cos(phi) * sin(theta); |
241 | + | Azz = cos(theta); |
242 | + | |
243 | + | sux = 0.0; |
244 | + | suy = 0.0; |
245 | + | suz = 1.0; |
246 | + | |
247 | + | ux = (Axx * sux) + (Ayx * suy) + (Azx * suz); |
248 | + | uy = (Axy * sux) + (Ayy * suy) + (Azy * suz); |
249 | + | uz = (Axz * sux) + (Ayz * suy) + (Azz * suz); |
250 | + | |
251 | uSqr = (ux * ux) + (uy * uy) + (uz * uz); | |
252 | ||
253 | u = sqrt(uSqr); | |
# | Line 613 | Line 663 | void SimSetup::gatherInfo(void){ | |
663 | } | |
664 | else{ | |
665 | sprintf(painCave.errMsg, | |
666 | < | "SimSetup Warning. Unrecognized Ensemble -> %s, " |
667 | < | "reverting to NVE for this simulation.\n", |
666 | > | "SimSetup Warning. Unrecognized Ensemble -> %s \n" |
667 | > | "\treverting to NVE for this simulation.\n", |
668 | ensemble); | |
669 | painCave.isFatal = 0; | |
670 | simError(); | |
# | Line 646 | Line 696 | void SimSetup::gatherInfo(void){ | |
696 | if (!the_components[i]->haveNMol()){ | |
697 | // we have a problem | |
698 | sprintf(painCave.errMsg, | |
699 | < | "SimSetup Error. No global NMol or component NMol" |
700 | < | " given. Cannot calculate the number of atoms.\n"); |
699 | > | "SimSetup Error. No global NMol or component NMol given.\n" |
700 | > | "\tCannot calculate the number of atoms.\n"); |
701 | painCave.isFatal = 1; | |
702 | simError(); | |
703 | } | |
# | Line 665 | Line 715 | void SimSetup::gatherInfo(void){ | |
715 | " Please give nMol in the components.\n"); | |
716 | painCave.isFatal = 1; | |
717 | simError(); | |
718 | + | } |
719 | + | |
720 | + | //check whether sample time, status time, thermal time and reset time are divisble by dt |
721 | + | if (!isDivisible(globals->getSampleTime(), globals->getDt())){ |
722 | + | sprintf(painCave.errMsg, |
723 | + | "Sample time is not divisible by dt.\n" |
724 | + | "\tThis will result in samples that are not uniformly\n" |
725 | + | "\tdistributed in time. If this is a problem, change\n" |
726 | + | "\tyour sampleTime variable.\n"); |
727 | + | painCave.isFatal = 0; |
728 | + | simError(); |
729 | } | |
730 | ||
731 | + | if (globals->haveStatusTime() && !isDivisible(globals->getSampleTime(), globals->getDt())){ |
732 | + | sprintf(painCave.errMsg, |
733 | + | "Status time is not divisible by dt.\n" |
734 | + | "\tThis will result in status reports that are not uniformly\n" |
735 | + | "\tdistributed in time. If this is a problem, change \n" |
736 | + | "\tyour statusTime variable.\n"); |
737 | + | painCave.isFatal = 0; |
738 | + | simError(); |
739 | + | } |
740 | + | |
741 | + | if (globals->haveThermalTime() && !isDivisible(globals->getThermalTime(), globals->getDt())){ |
742 | + | sprintf(painCave.errMsg, |
743 | + | "Thermal time is not divisible by dt.\n" |
744 | + | "\tThis will result in thermalizations that are not uniformly\n" |
745 | + | "\tdistributed in time. If this is a problem, change \n" |
746 | + | "\tyour thermalTime variable.\n"); |
747 | + | painCave.isFatal = 0; |
748 | + | simError(); |
749 | + | } |
750 | + | |
751 | + | if (globals->haveResetTime() && !isDivisible(globals->getResetTime(), globals->getDt())){ |
752 | + | sprintf(painCave.errMsg, |
753 | + | "Reset time is not divisible by dt.\n" |
754 | + | "\tThis will result in integrator resets that are not uniformly\n" |
755 | + | "\tdistributed in time. If this is a problem, change\n" |
756 | + | "\tyour resetTime variable.\n"); |
757 | + | painCave.isFatal = 0; |
758 | + | simError(); |
759 | + | } |
760 | + | |
761 | // set the status, sample, and thermal kick times | |
762 | ||
763 | for (i = 0; i < nInfo; i++){ | |
# | Line 696 | Line 787 | void SimSetup::gatherInfo(void){ | |
787 | } | |
788 | ||
789 | // check for the temperature set flag | |
790 | < | |
790 | > | |
791 | if (globals->haveTempSet()) | |
792 | info[i].setTemp = globals->getTempSet(); | |
793 | ||
794 | < | // get some of the tricky things that may still be in the globals |
794 | > | // check for the extended State init |
795 | ||
796 | < | double boxVector[3]; |
797 | < | if (globals->haveBox()){ |
798 | < | boxVector[0] = globals->getBox(); |
708 | < | boxVector[1] = globals->getBox(); |
709 | < | boxVector[2] = globals->getBox(); |
710 | < | |
711 | < | info[i].setBox(boxVector); |
712 | < | } |
713 | < | else if (globals->haveDensity()){ |
714 | < | double vol; |
715 | < | vol = (double) tot_nmol / globals->getDensity(); |
716 | < | boxVector[0] = pow(vol, (1.0 / 3.0)); |
717 | < | boxVector[1] = boxVector[0]; |
718 | < | boxVector[2] = boxVector[0]; |
719 | < | |
720 | < | info[i].setBox(boxVector); |
721 | < | } |
722 | < | else{ |
723 | < | if (!globals->haveBoxX()){ |
724 | < | sprintf(painCave.errMsg, |
725 | < | "SimSetup error, no periodic BoxX size given.\n"); |
726 | < | painCave.isFatal = 1; |
727 | < | simError(); |
728 | < | } |
729 | < | boxVector[0] = globals->getBoxX(); |
730 | < | |
731 | < | if (!globals->haveBoxY()){ |
732 | < | sprintf(painCave.errMsg, |
733 | < | "SimSetup error, no periodic BoxY size given.\n"); |
734 | < | painCave.isFatal = 1; |
735 | < | simError(); |
736 | < | } |
737 | < | boxVector[1] = globals->getBoxY(); |
738 | < | |
739 | < | if (!globals->haveBoxZ()){ |
740 | < | sprintf(painCave.errMsg, |
741 | < | "SimSetup error, no periodic BoxZ size given.\n"); |
742 | < | painCave.isFatal = 1; |
743 | < | simError(); |
744 | < | } |
745 | < | boxVector[2] = globals->getBoxZ(); |
746 | < | |
747 | < | info[i].setBox(boxVector); |
748 | < | } |
796 | > | info[i].useInitXSstate = globals->getUseInitXSstate(); |
797 | > | info[i].orthoTolerance = globals->getOrthoBoxTolerance(); |
798 | > | |
799 | } | |
800 | < | |
800 | > | |
801 | //setup seed for random number generator | |
802 | int seedValue; | |
803 | ||
# | Line 822 | Line 872 | void SimSetup::finalInfoCheck(void){ | |
872 | ||
873 | if (!globals->haveECR()){ | |
874 | sprintf(painCave.errMsg, | |
875 | < | "SimSetup Warning: using default value of 1/2 the smallest " |
876 | < | "box length for the electrostaticCutoffRadius.\n" |
877 | < | "I hope you have a very fast processor!\n"); |
875 | > | "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n" |
876 | > | "\tOOPSE will use a default value of 15.0 angstroms" |
877 | > | "\tfor the electrostaticCutoffRadius.\n"); |
878 | painCave.isFatal = 0; | |
879 | simError(); | |
880 | < | double smallest; |
831 | < | smallest = info[i].boxL[0]; |
832 | < | if (info[i].boxL[1] <= smallest) |
833 | < | smallest = info[i].boxL[1]; |
834 | < | if (info[i].boxL[2] <= smallest) |
835 | < | smallest = info[i].boxL[2]; |
836 | < | theEcr = 0.5 * smallest; |
880 | > | theEcr = 15.0; |
881 | } | |
882 | else{ | |
883 | theEcr = globals->getECR(); | |
# | Line 841 | Line 885 | void SimSetup::finalInfoCheck(void){ | |
885 | ||
886 | if (!globals->haveEST()){ | |
887 | sprintf(painCave.errMsg, | |
888 | < | "SimSetup Warning: using default value of 0.05 * the " |
889 | < | "electrostaticCutoffRadius for the electrostaticSkinThickness\n"); |
888 | > | "SimSetup Warning: No value was set for electrostaticSkinThickness.\n" |
889 | > | "\tOOPSE will use a default value of\n" |
890 | > | "\t0.05 * electrostaticCutoffRadius\n" |
891 | > | "\tfor the electrostaticSkinThickness\n"); |
892 | painCave.isFatal = 0; | |
893 | simError(); | |
894 | theEst = 0.05 * theEcr; | |
# | Line 851 | Line 897 | void SimSetup::finalInfoCheck(void){ | |
897 | theEst = globals->getEST(); | |
898 | } | |
899 | ||
900 | < | info[i].setEcr(theEcr, theEst); |
900 | > | info[i].setDefaultEcr(theEcr, theEst); |
901 | ||
902 | if (!globals->haveDielectric()){ | |
903 | sprintf(painCave.errMsg, | |
904 | < | "SimSetup Error: You are trying to use Reaction Field without" |
905 | < | "setting a dielectric constant!\n"); |
904 | > | "SimSetup Error: No Dielectric constant was set.\n" |
905 | > | "\tYou are trying to use Reaction Field without" |
906 | > | "\tsetting a dielectric constant!\n"); |
907 | painCave.isFatal = 1; | |
908 | simError(); | |
909 | } | |
# | Line 866 | Line 913 | void SimSetup::finalInfoCheck(void){ | |
913 | if (usesDipoles){ | |
914 | if (!globals->haveECR()){ | |
915 | sprintf(painCave.errMsg, | |
916 | < | "SimSetup Warning: using default value of 1/2 the smallest " |
917 | < | "box length for the electrostaticCutoffRadius.\n" |
918 | < | "I hope you have a very fast processor!\n"); |
919 | < | painCave.isFatal = 0; |
920 | < | simError(); |
921 | < | double smallest; |
875 | < | smallest = info[i].boxL[0]; |
876 | < | if (info[i].boxL[1] <= smallest) |
877 | < | smallest = info[i].boxL[1]; |
878 | < | if (info[i].boxL[2] <= smallest) |
879 | < | smallest = info[i].boxL[2]; |
880 | < | theEcr = 0.5 * smallest; |
916 | > | "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n" |
917 | > | "\tOOPSE will use a default value of 15.0 angstroms" |
918 | > | "\tfor the electrostaticCutoffRadius.\n"); |
919 | > | painCave.isFatal = 0; |
920 | > | simError(); |
921 | > | theEcr = 15.0; |
922 | } | |
923 | else{ | |
924 | theEcr = globals->getECR(); | |
925 | } | |
926 | < | |
926 | > | |
927 | if (!globals->haveEST()){ | |
928 | sprintf(painCave.errMsg, | |
929 | < | "SimSetup Warning: using default value of 0.05 * the " |
930 | < | "electrostaticCutoffRadius for the " |
931 | < | "electrostaticSkinThickness\n"); |
929 | > | "SimSetup Warning: No value was set for electrostaticSkinThickness.\n" |
930 | > | "\tOOPSE will use a default value of\n" |
931 | > | "\t0.05 * electrostaticCutoffRadius\n" |
932 | > | "\tfor the electrostaticSkinThickness\n"); |
933 | painCave.isFatal = 0; | |
934 | simError(); | |
935 | theEst = 0.05 * theEcr; | |
# | Line 895 | Line 937 | void SimSetup::finalInfoCheck(void){ | |
937 | else{ | |
938 | theEst = globals->getEST(); | |
939 | } | |
940 | < | |
941 | < | info[i].setEcr(theEcr, theEst); |
940 | > | |
941 | > | info[i].setDefaultEcr(theEcr, theEst); |
942 | } | |
943 | } | |
944 | } | |
903 | – | |
945 | #ifdef IS_MPI | |
946 | strcpy(checkPointMsg, "post processing checks out"); | |
947 | MPIcheckPoint(); | |
948 | #endif // is_mpi | |
949 | } | |
950 | < | |
950 | > | |
951 | void SimSetup::initSystemCoords(void){ | |
952 | int i; | |
953 | ||
# | Line 934 | Line 975 | void SimSetup::initSystemCoords(void){ | |
975 | delete fileInit; | |
976 | } | |
977 | else{ | |
978 | < | #ifdef IS_MPI |
938 | < | |
978 | > | |
979 | // no init from bass | |
980 | < | |
980 | > | |
981 | sprintf(painCave.errMsg, | |
982 | < | "Cannot intialize a parallel simulation without an initial configuration file.\n"); |
982 | > | "Cannot intialize a simulation without an initial configuration file.\n"); |
983 | painCave.isFatal = 1;; | |
984 | simError(); | |
985 | < | |
946 | < | #else |
947 | < | |
948 | < | initFromBass(); |
949 | < | |
950 | < | |
951 | < | #endif |
985 | > | |
986 | } | |
987 | ||
988 | #ifdef IS_MPI | |
# | Line 1244 | Line 1278 | void SimSetup::mpiMolDivide(void){ | |
1278 | ||
1279 | if (local_atoms != info[0].n_atoms){ | |
1280 | sprintf(painCave.errMsg, | |
1281 | < | "SimSetup error: mpiSim's localAtom (%d) and SimSetup's" |
1282 | < | " localAtom (%d) are not equal.\n", |
1281 | > | "SimSetup error: mpiSim's localAtom (%d) and SimSetup's\n" |
1282 | > | "\tlocalAtom (%d) are not equal.\n", |
1283 | info[0].n_atoms, local_atoms); | |
1284 | painCave.isFatal = 1; | |
1285 | simError(); | |
# | Line 1387 | Line 1421 | void SimSetup::makeIntegrator(void){ | |
1421 | else{ | |
1422 | sprintf(painCave.errMsg, | |
1423 | "SimSetup error: If you use the NVT\n" | |
1424 | < | " ensemble, you must set tauThermostat.\n"); |
1424 | > | "\tensemble, you must set tauThermostat.\n"); |
1425 | painCave.isFatal = 1; | |
1426 | simError(); | |
1427 | } | |
# | Line 1410 | Line 1444 | void SimSetup::makeIntegrator(void){ | |
1444 | else{ | |
1445 | sprintf(painCave.errMsg, | |
1446 | "SimSetup error: If you use a constant pressure\n" | |
1447 | < | " ensemble, you must set targetPressure in the BASS file.\n"); |
1447 | > | "\tensemble, you must set targetPressure in the BASS file.\n"); |
1448 | painCave.isFatal = 1; | |
1449 | simError(); | |
1450 | } | |
# | Line 1420 | Line 1454 | void SimSetup::makeIntegrator(void){ | |
1454 | else{ | |
1455 | sprintf(painCave.errMsg, | |
1456 | "SimSetup error: If you use an NPT\n" | |
1457 | < | " ensemble, you must set tauThermostat.\n"); |
1457 | > | "\tensemble, you must set tauThermostat.\n"); |
1458 | painCave.isFatal = 1; | |
1459 | simError(); | |
1460 | } | |
# | Line 1430 | Line 1464 | void SimSetup::makeIntegrator(void){ | |
1464 | else{ | |
1465 | sprintf(painCave.errMsg, | |
1466 | "SimSetup error: If you use an NPT\n" | |
1467 | < | " ensemble, you must set tauBarostat.\n"); |
1467 | > | "\tensemble, you must set tauBarostat.\n"); |
1468 | painCave.isFatal = 1; | |
1469 | simError(); | |
1470 | } | |
# | Line 1453 | Line 1487 | void SimSetup::makeIntegrator(void){ | |
1487 | else{ | |
1488 | sprintf(painCave.errMsg, | |
1489 | "SimSetup error: If you use a constant pressure\n" | |
1490 | < | " ensemble, you must set targetPressure in the BASS file.\n"); |
1490 | > | "\tensemble, you must set targetPressure in the BASS file.\n"); |
1491 | painCave.isFatal = 1; | |
1492 | simError(); | |
1493 | } | |
1494 | ||
1495 | if (globals->haveTauThermostat()) | |
1496 | myNPTf->setTauThermostat(globals->getTauThermostat()); | |
1497 | + | |
1498 | else{ | |
1499 | sprintf(painCave.errMsg, | |
1500 | "SimSetup error: If you use an NPT\n" | |
1501 | < | " ensemble, you must set tauThermostat.\n"); |
1501 | > | "\tensemble, you must set tauThermostat.\n"); |
1502 | painCave.isFatal = 1; | |
1503 | simError(); | |
1504 | } | |
1505 | ||
1506 | if (globals->haveTauBarostat()) | |
1507 | myNPTf->setTauBarostat(globals->getTauBarostat()); | |
1508 | + | |
1509 | else{ | |
1510 | sprintf(painCave.errMsg, | |
1511 | "SimSetup error: If you use an NPT\n" | |
1512 | < | " ensemble, you must set tauBarostat.\n"); |
1512 | > | "\tensemble, you must set tauBarostat.\n"); |
1513 | painCave.isFatal = 1; | |
1514 | simError(); | |
1515 | } | |
# | Line 1496 | Line 1532 | void SimSetup::makeIntegrator(void){ | |
1532 | else{ | |
1533 | sprintf(painCave.errMsg, | |
1534 | "SimSetup error: If you use a constant pressure\n" | |
1535 | < | " ensemble, you must set targetPressure in the BASS file.\n"); |
1535 | > | "\tensemble, you must set targetPressure in the BASS file.\n"); |
1536 | painCave.isFatal = 1; | |
1537 | simError(); | |
1538 | } | |
# | Line 1506 | Line 1542 | void SimSetup::makeIntegrator(void){ | |
1542 | else{ | |
1543 | sprintf(painCave.errMsg, | |
1544 | "SimSetup error: If you use an NPT\n" | |
1545 | < | " ensemble, you must set tauThermostat.\n"); |
1545 | > | "\tensemble, you must set tauThermostat.\n"); |
1546 | painCave.isFatal = 1; | |
1547 | simError(); | |
1548 | } | |
# | Line 1516 | Line 1552 | void SimSetup::makeIntegrator(void){ | |
1552 | else{ | |
1553 | sprintf(painCave.errMsg, | |
1554 | "SimSetup error: If you use an NPT\n" | |
1555 | < | " ensemble, you must set tauBarostat.\n"); |
1555 | > | "\tensemble, you must set tauBarostat.\n"); |
1556 | painCave.isFatal = 1; | |
1557 | simError(); | |
1558 | } | |
# | Line 1569 | Line 1605 | void SimSetup::setupZConstraint(SimInfo& theInfo){ | |
1605 | } | |
1606 | else{ | |
1607 | sprintf(painCave.errMsg, | |
1608 | < | "ZConstraint error: If you use an ZConstraint\n" |
1609 | < | " , you must set sample time.\n"); |
1608 | > | "ZConstraint error: If you use a ZConstraint,\n" |
1609 | > | "\tyou must set zconsTime.\n"); |
1610 | painCave.isFatal = 1; | |
1611 | simError(); | |
1612 | } | |
# | Line 1585 | Line 1621 | void SimSetup::setupZConstraint(SimInfo& theInfo){ | |
1621 | else{ | |
1622 | double defaultZConsTol = 0.01; | |
1623 | sprintf(painCave.errMsg, | |
1624 | < | "ZConstraint Waring: Tolerance for z-constraint methodl is not specified\n" |
1625 | < | " , default value %f is used.\n", |
1624 | > | "ZConstraint Warning: Tolerance for z-constraint method is not specified.\n" |
1625 | > | "\tOOPSE will use a default value of %f.\n" |
1626 | > | "\tTo set the tolerance, use the zconsTol variable.\n", |
1627 | defaultZConsTol); | |
1628 | painCave.isFatal = 0; | |
1629 | simError(); | |
# | Line 1604 | Line 1641 | void SimSetup::setupZConstraint(SimInfo& theInfo){ | |
1641 | } | |
1642 | else{ | |
1643 | sprintf(painCave.errMsg, | |
1644 | < | "ZConstraint Warning: User does not set force Subtraction policy, " |
1645 | < | "PolicyByMass is used\n"); |
1644 | > | "ZConstraint Warning: No force subtraction policy was set.\n" |
1645 | > | "\tOOPSE will use PolicyByMass.\n" |
1646 | > | "\tTo set the policy, use the zconsForcePolicy variable.\n"); |
1647 | painCave.isFatal = 0; | |
1648 | simError(); | |
1649 | zconsForcePolicy->setData("BYMASS"); | |
# | Line 1649 | Line 1687 | void SimSetup::setupZConstraint(SimInfo& theInfo){ | |
1687 | //check the uniqueness of index | |
1688 | if(!zconsParaData->isIndexUnique()){ | |
1689 | sprintf(painCave.errMsg, | |
1690 | < | "ZConstraint Error: molIndex is not unique\n"); |
1690 | > | "ZConstraint Error: molIndex is not unique!\n"); |
1691 | painCave.isFatal = 1; | |
1692 | simError(); | |
1693 | } |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |