# | Line 1 | Line 1 | |
---|---|---|
1 | #include <algorithm> | |
2 | < | #include <cstdlib> |
2 | > | #include <stdlib.h> |
3 | #include <iostream> | |
4 | < | #include <cmath> |
4 | > | #include <math.h> |
5 | #include <string> | |
6 | #include <sprng.h> | |
7 | – | |
7 | #include "SimSetup.hpp" | |
8 | #include "ReadWrite.hpp" | |
9 | #include "parse_me.h" | |
10 | #include "Integrator.hpp" | |
11 | #include "simError.h" | |
12 | + | //#include "ConjugateMinimizer.hpp" |
13 | + | #include "OOPSEMinimizer.hpp" |
14 | ||
15 | #ifdef IS_MPI | |
16 | #include "mpiBASS.h" | |
# | Line 22 | Line 23 | |
23 | #define NVT_ENS 1 | |
24 | #define NPTi_ENS 2 | |
25 | #define NPTf_ENS 3 | |
26 | + | #define NPTxyz_ENS 4 |
27 | ||
26 | – | #define FF_DUFF 0 |
27 | – | #define FF_LJ 1 |
28 | – | #define FF_EAM 2 |
28 | ||
29 | + | #define FF_DUFF 0 |
30 | + | #define FF_LJ 1 |
31 | + | #define FF_EAM 2 |
32 | + | #define FF_H2O 3 |
33 | + | |
34 | using namespace std; | |
35 | ||
36 | + | /** |
37 | + | * Check whether dividend is divisble by divisor or not |
38 | + | */ |
39 | + | bool isDivisible(double dividend, double divisor){ |
40 | + | double tolerance = 0.000001; |
41 | + | double quotient; |
42 | + | double diff; |
43 | + | int intQuotient; |
44 | + | |
45 | + | quotient = dividend / divisor; |
46 | + | |
47 | + | if (quotient < 0) |
48 | + | quotient = -quotient; |
49 | + | |
50 | + | intQuotient = int (quotient + tolerance); |
51 | + | |
52 | + | diff = fabs(fabs(dividend) - intQuotient * fabs(divisor)); |
53 | + | |
54 | + | if (diff <= tolerance) |
55 | + | return true; |
56 | + | else |
57 | + | return false; |
58 | + | } |
59 | + | |
60 | SimSetup::SimSetup(){ | |
61 | + | |
62 | + | initSuspend = false; |
63 | isInfoArray = 0; | |
64 | nInfo = 1; | |
65 | ||
# | Line 52 | Line 82 | void SimSetup::setSimInfo(SimInfo* the_info, int theNi | |
82 | info = the_info; | |
83 | nInfo = theNinfo; | |
84 | isInfoArray = 1; | |
85 | + | initSuspend = true; |
86 | } | |
87 | ||
88 | ||
# | Line 105 | Line 136 | void SimSetup::createSim(void){ | |
136 | ||
137 | // initialize the system coordinates | |
138 | ||
139 | < | if (!isInfoArray){ |
139 | > | if ( !initSuspend ){ |
140 | initSystemCoords(); | |
141 | ||
142 | if( !(globals->getUseInitTime()) ) | |
# | Line 116 | Line 147 | void SimSetup::createSim(void){ | |
147 | ||
148 | makeOutNames(); | |
149 | ||
150 | < | // make the integrator |
151 | < | |
152 | < | makeIntegrator(); |
153 | < | |
150 | > | if (globals->haveMinimizer()) |
151 | > | // make minimizer |
152 | > | makeMinimizer(); |
153 | > | else |
154 | > | // make the integrator |
155 | > | makeIntegrator(); |
156 | > | |
157 | #ifdef IS_MPI | |
158 | mpiSim->mpiRefresh(); | |
159 | #endif | |
# | Line 146 | Line 180 | void SimSetup::makeMolecules(void){ | |
180 | bend_set* theBends; | |
181 | torsion_set* theTorsions; | |
182 | ||
149 | – | |
183 | //init the forceField paramters | |
184 | ||
185 | the_ff->readParams(); | |
# | Line 154 | Line 187 | void SimSetup::makeMolecules(void){ | |
187 | ||
188 | // init the atoms | |
189 | ||
190 | + | double phi, theta, psi; |
191 | + | double sux, suy, suz; |
192 | + | double Axx, Axy, Axz, Ayx, Ayy, Ayz, Azx, Azy, Azz; |
193 | double ux, uy, uz, u, uSqr; | |
194 | ||
195 | for (k = 0; k < nInfo; k++){ | |
# | Line 190 | Line 226 | void SimSetup::makeMolecules(void){ | |
226 | info[k].n_oriented++; | |
227 | molInfo.myAtoms[j] = dAtom; | |
228 | ||
229 | < | ux = currentAtom->getOrntX(); |
230 | < | uy = currentAtom->getOrntY(); |
231 | < | uz = currentAtom->getOrntZ(); |
229 | > | // Directional Atoms have standard unit vectors which are oriented |
230 | > | // in space using the three Euler angles. We assume the standard |
231 | > | // unit vector was originally along the z axis below. |
232 | > | |
233 | > | phi = currentAtom->getEulerPhi() * M_PI / 180.0; |
234 | > | theta = currentAtom->getEulerTheta() * M_PI / 180.0; |
235 | > | psi = currentAtom->getEulerPsi()* M_PI / 180.0; |
236 | > | |
237 | > | Axx = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi)); |
238 | > | Axy = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi)); |
239 | > | Axz = sin(theta) * sin(psi); |
240 | > | |
241 | > | Ayx = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi)); |
242 | > | Ayy = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi)); |
243 | > | Ayz = sin(theta) * cos(psi); |
244 | > | |
245 | > | Azx = sin(phi) * sin(theta); |
246 | > | Azy = -cos(phi) * sin(theta); |
247 | > | Azz = cos(theta); |
248 | > | |
249 | > | sux = 0.0; |
250 | > | suy = 0.0; |
251 | > | suz = 1.0; |
252 | > | |
253 | > | ux = (Axx * sux) + (Ayx * suy) + (Azx * suz); |
254 | > | uy = (Axy * sux) + (Ayy * suy) + (Azy * suz); |
255 | > | uz = (Axz * sux) + (Ayz * suy) + (Azz * suz); |
256 | ||
257 | uSqr = (ux * ux) + (uy * uy) + (uz * uz); | |
258 | ||
# | Line 581 | Line 641 | void SimSetup::gatherInfo(void){ | |
641 | else if (!strcasecmp(force_field, "EAM")){ | |
642 | ffCase = FF_EAM; | |
643 | } | |
644 | + | else if (!strcasecmp(force_field, "WATER")){ |
645 | + | ffCase = FF_H2O; |
646 | + | } |
647 | else{ | |
648 | sprintf(painCave.errMsg, "SimSetup Error. Unrecognized force field -> %s\n", | |
649 | force_field); | |
# | Line 604 | Line 667 | void SimSetup::gatherInfo(void){ | |
667 | else if (!strcasecmp(ensemble, "NPTf")){ | |
668 | ensembleCase = NPTf_ENS; | |
669 | } | |
670 | + | else if (!strcasecmp(ensemble, "NPTxyz")){ |
671 | + | ensembleCase = NPTxyz_ENS; |
672 | + | } |
673 | else{ | |
674 | sprintf(painCave.errMsg, | |
675 | < | "SimSetup Warning. Unrecognized Ensemble -> %s, " |
676 | < | "reverting to NVE for this simulation.\n", |
675 | > | "SimSetup Warning. Unrecognized Ensemble -> %s \n" |
676 | > | "\treverting to NVE for this simulation.\n", |
677 | ensemble); | |
678 | painCave.isFatal = 0; | |
679 | simError(); | |
# | Line 639 | Line 705 | void SimSetup::gatherInfo(void){ | |
705 | if (!the_components[i]->haveNMol()){ | |
706 | // we have a problem | |
707 | sprintf(painCave.errMsg, | |
708 | < | "SimSetup Error. No global NMol or component NMol" |
709 | < | " given. Cannot calculate the number of atoms.\n"); |
708 | > | "SimSetup Error. No global NMol or component NMol given.\n" |
709 | > | "\tCannot calculate the number of atoms.\n"); |
710 | painCave.isFatal = 1; | |
711 | simError(); | |
712 | } | |
# | Line 658 | Line 724 | void SimSetup::gatherInfo(void){ | |
724 | " Please give nMol in the components.\n"); | |
725 | painCave.isFatal = 1; | |
726 | simError(); | |
727 | + | } |
728 | + | |
729 | + | //check whether sample time, status time, thermal time and reset time are divisble by dt |
730 | + | if (!isDivisible(globals->getSampleTime(), globals->getDt())){ |
731 | + | sprintf(painCave.errMsg, |
732 | + | "Sample time is not divisible by dt.\n" |
733 | + | "\tThis will result in samples that are not uniformly\n" |
734 | + | "\tdistributed in time. If this is a problem, change\n" |
735 | + | "\tyour sampleTime variable.\n"); |
736 | + | painCave.isFatal = 0; |
737 | + | simError(); |
738 | + | } |
739 | + | |
740 | + | if (globals->haveStatusTime() && !isDivisible(globals->getSampleTime(), globals->getDt())){ |
741 | + | sprintf(painCave.errMsg, |
742 | + | "Status time is not divisible by dt.\n" |
743 | + | "\tThis will result in status reports that are not uniformly\n" |
744 | + | "\tdistributed in time. If this is a problem, change \n" |
745 | + | "\tyour statusTime variable.\n"); |
746 | + | painCave.isFatal = 0; |
747 | + | simError(); |
748 | } | |
749 | ||
750 | + | if (globals->haveThermalTime() && !isDivisible(globals->getThermalTime(), globals->getDt())){ |
751 | + | sprintf(painCave.errMsg, |
752 | + | "Thermal time is not divisible by dt.\n" |
753 | + | "\tThis will result in thermalizations that are not uniformly\n" |
754 | + | "\tdistributed in time. If this is a problem, change \n" |
755 | + | "\tyour thermalTime variable.\n"); |
756 | + | painCave.isFatal = 0; |
757 | + | simError(); |
758 | + | } |
759 | + | |
760 | + | if (globals->haveResetTime() && !isDivisible(globals->getResetTime(), globals->getDt())){ |
761 | + | sprintf(painCave.errMsg, |
762 | + | "Reset time is not divisible by dt.\n" |
763 | + | "\tThis will result in integrator resets that are not uniformly\n" |
764 | + | "\tdistributed in time. If this is a problem, change\n" |
765 | + | "\tyour resetTime variable.\n"); |
766 | + | painCave.isFatal = 0; |
767 | + | simError(); |
768 | + | } |
769 | + | |
770 | // set the status, sample, and thermal kick times | |
771 | ||
772 | for (i = 0; i < nInfo; i++){ | |
# | Line 689 | Line 796 | void SimSetup::gatherInfo(void){ | |
796 | } | |
797 | ||
798 | // check for the temperature set flag | |
799 | < | |
799 | > | |
800 | if (globals->haveTempSet()) | |
801 | info[i].setTemp = globals->getTempSet(); | |
802 | ||
803 | < | // get some of the tricky things that may still be in the globals |
803 | > | // check for the extended State init |
804 | ||
805 | < | double boxVector[3]; |
806 | < | if (globals->haveBox()){ |
807 | < | boxVector[0] = globals->getBox(); |
701 | < | boxVector[1] = globals->getBox(); |
702 | < | boxVector[2] = globals->getBox(); |
703 | < | |
704 | < | info[i].setBox(boxVector); |
705 | < | } |
706 | < | else if (globals->haveDensity()){ |
707 | < | double vol; |
708 | < | vol = (double) tot_nmol / globals->getDensity(); |
709 | < | boxVector[0] = pow(vol, (1.0 / 3.0)); |
710 | < | boxVector[1] = boxVector[0]; |
711 | < | boxVector[2] = boxVector[0]; |
712 | < | |
713 | < | info[i].setBox(boxVector); |
714 | < | } |
715 | < | else{ |
716 | < | if (!globals->haveBoxX()){ |
717 | < | sprintf(painCave.errMsg, |
718 | < | "SimSetup error, no periodic BoxX size given.\n"); |
719 | < | painCave.isFatal = 1; |
720 | < | simError(); |
721 | < | } |
722 | < | boxVector[0] = globals->getBoxX(); |
723 | < | |
724 | < | if (!globals->haveBoxY()){ |
725 | < | sprintf(painCave.errMsg, |
726 | < | "SimSetup error, no periodic BoxY size given.\n"); |
727 | < | painCave.isFatal = 1; |
728 | < | simError(); |
729 | < | } |
730 | < | boxVector[1] = globals->getBoxY(); |
731 | < | |
732 | < | if (!globals->haveBoxZ()){ |
733 | < | sprintf(painCave.errMsg, |
734 | < | "SimSetup error, no periodic BoxZ size given.\n"); |
735 | < | painCave.isFatal = 1; |
736 | < | simError(); |
737 | < | } |
738 | < | boxVector[2] = globals->getBoxZ(); |
739 | < | |
740 | < | info[i].setBox(boxVector); |
741 | < | } |
805 | > | info[i].useInitXSstate = globals->getUseInitXSstate(); |
806 | > | info[i].orthoTolerance = globals->getOrthoBoxTolerance(); |
807 | > | |
808 | } | |
809 | < | |
809 | > | |
810 | //setup seed for random number generator | |
811 | int seedValue; | |
812 | ||
# | Line 780 | Line 846 | void SimSetup::gatherInfo(void){ | |
846 | for (int i = 0; i < nInfo; i++){ | |
847 | info[i].setSeed(seedValue); | |
848 | } | |
849 | < | |
849 | > | |
850 | #ifdef IS_MPI | |
851 | < | strcpy(checkPointMsg, "Succesfully gathered all information from Bass\n"); |
851 | > | strcpy(checkPointMsg, "Successfully gathered all information from Bass\n"); |
852 | MPIcheckPoint(); | |
853 | #endif // is_mpi | |
854 | } | |
# | Line 815 | Line 881 | void SimSetup::finalInfoCheck(void){ | |
881 | ||
882 | if (!globals->haveECR()){ | |
883 | sprintf(painCave.errMsg, | |
884 | < | "SimSetup Warning: using default value of 1/2 the smallest " |
885 | < | "box length for the electrostaticCutoffRadius.\n" |
886 | < | "I hope you have a very fast processor!\n"); |
884 | > | "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n" |
885 | > | "\tOOPSE will use a default value of 15.0 angstroms" |
886 | > | "\tfor the electrostaticCutoffRadius.\n"); |
887 | painCave.isFatal = 0; | |
888 | simError(); | |
889 | < | double smallest; |
824 | < | smallest = info[i].boxL[0]; |
825 | < | if (info[i].boxL[1] <= smallest) |
826 | < | smallest = info[i].boxL[1]; |
827 | < | if (info[i].boxL[2] <= smallest) |
828 | < | smallest = info[i].boxL[2]; |
829 | < | theEcr = 0.5 * smallest; |
889 | > | theEcr = 15.0; |
890 | } | |
891 | else{ | |
892 | theEcr = globals->getECR(); | |
# | Line 834 | Line 894 | void SimSetup::finalInfoCheck(void){ | |
894 | ||
895 | if (!globals->haveEST()){ | |
896 | sprintf(painCave.errMsg, | |
897 | < | "SimSetup Warning: using default value of 0.05 * the " |
898 | < | "electrostaticCutoffRadius for the electrostaticSkinThickness\n"); |
897 | > | "SimSetup Warning: No value was set for electrostaticSkinThickness.\n" |
898 | > | "\tOOPSE will use a default value of\n" |
899 | > | "\t0.05 * electrostaticCutoffRadius\n" |
900 | > | "\tfor the electrostaticSkinThickness\n"); |
901 | painCave.isFatal = 0; | |
902 | simError(); | |
903 | theEst = 0.05 * theEcr; | |
# | Line 844 | Line 906 | void SimSetup::finalInfoCheck(void){ | |
906 | theEst = globals->getEST(); | |
907 | } | |
908 | ||
909 | < | info[i].setEcr(theEcr, theEst); |
909 | > | info[i].setDefaultEcr(theEcr, theEst); |
910 | ||
911 | if (!globals->haveDielectric()){ | |
912 | sprintf(painCave.errMsg, | |
913 | < | "SimSetup Error: You are trying to use Reaction Field without" |
914 | < | "setting a dielectric constant!\n"); |
913 | > | "SimSetup Error: No Dielectric constant was set.\n" |
914 | > | "\tYou are trying to use Reaction Field without" |
915 | > | "\tsetting a dielectric constant!\n"); |
916 | painCave.isFatal = 1; | |
917 | simError(); | |
918 | } | |
# | Line 859 | Line 922 | void SimSetup::finalInfoCheck(void){ | |
922 | if (usesDipoles){ | |
923 | if (!globals->haveECR()){ | |
924 | sprintf(painCave.errMsg, | |
925 | < | "SimSetup Warning: using default value of 1/2 the smallest " |
926 | < | "box length for the electrostaticCutoffRadius.\n" |
927 | < | "I hope you have a very fast processor!\n"); |
928 | < | painCave.isFatal = 0; |
929 | < | simError(); |
930 | < | double smallest; |
868 | < | smallest = info[i].boxL[0]; |
869 | < | if (info[i].boxL[1] <= smallest) |
870 | < | smallest = info[i].boxL[1]; |
871 | < | if (info[i].boxL[2] <= smallest) |
872 | < | smallest = info[i].boxL[2]; |
873 | < | theEcr = 0.5 * smallest; |
925 | > | "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n" |
926 | > | "\tOOPSE will use a default value of 15.0 angstroms" |
927 | > | "\tfor the electrostaticCutoffRadius.\n"); |
928 | > | painCave.isFatal = 0; |
929 | > | simError(); |
930 | > | theEcr = 15.0; |
931 | } | |
932 | else{ | |
933 | theEcr = globals->getECR(); | |
934 | } | |
935 | < | |
935 | > | |
936 | if (!globals->haveEST()){ | |
937 | sprintf(painCave.errMsg, | |
938 | < | "SimSetup Warning: using default value of 0.05 * the " |
939 | < | "electrostaticCutoffRadius for the " |
940 | < | "electrostaticSkinThickness\n"); |
938 | > | "SimSetup Warning: No value was set for electrostaticSkinThickness.\n" |
939 | > | "\tOOPSE will use a default value of\n" |
940 | > | "\t0.05 * electrostaticCutoffRadius\n" |
941 | > | "\tfor the electrostaticSkinThickness\n"); |
942 | painCave.isFatal = 0; | |
943 | simError(); | |
944 | theEst = 0.05 * theEcr; | |
# | Line 888 | Line 946 | void SimSetup::finalInfoCheck(void){ | |
946 | else{ | |
947 | theEst = globals->getEST(); | |
948 | } | |
949 | < | |
950 | < | info[i].setEcr(theEcr, theEst); |
949 | > | |
950 | > | info[i].setDefaultEcr(theEcr, theEst); |
951 | } | |
952 | } | |
953 | } | |
896 | – | |
954 | #ifdef IS_MPI | |
955 | strcpy(checkPointMsg, "post processing checks out"); | |
956 | MPIcheckPoint(); | |
957 | #endif // is_mpi | |
958 | } | |
959 | < | |
959 | > | |
960 | void SimSetup::initSystemCoords(void){ | |
961 | int i; | |
962 | ||
# | Line 927 | Line 984 | void SimSetup::initSystemCoords(void){ | |
984 | delete fileInit; | |
985 | } | |
986 | else{ | |
987 | < | #ifdef IS_MPI |
931 | < | |
987 | > | |
988 | // no init from bass | |
989 | < | |
989 | > | |
990 | sprintf(painCave.errMsg, | |
991 | < | "Cannot intialize a parallel simulation without an initial configuration file.\n"); |
991 | > | "Cannot intialize a simulation without an initial configuration file.\n"); |
992 | painCave.isFatal = 1;; | |
993 | simError(); | |
994 | < | |
939 | < | #else |
940 | < | |
941 | < | initFromBass(); |
942 | < | |
943 | < | |
944 | < | #endif |
994 | > | |
995 | } | |
996 | ||
997 | #ifdef IS_MPI | |
# | Line 1093 | Line 1143 | void SimSetup::createFF(void){ | |
1143 | ||
1144 | case FF_EAM: | |
1145 | the_ff = new EAM_FF(); | |
1146 | + | break; |
1147 | + | |
1148 | + | case FF_H2O: |
1149 | + | the_ff = new WATER(); |
1150 | break; | |
1151 | ||
1152 | default: | |
# | Line 1237 | Line 1291 | void SimSetup::mpiMolDivide(void){ | |
1291 | ||
1292 | if (local_atoms != info[0].n_atoms){ | |
1293 | sprintf(painCave.errMsg, | |
1294 | < | "SimSetup error: mpiSim's localAtom (%d) and SimSetup's" |
1295 | < | " localAtom (%d) are not equal.\n", |
1294 | > | "SimSetup error: mpiSim's localAtom (%d) and SimSetup's\n" |
1295 | > | "\tlocalAtom (%d) are not equal.\n", |
1296 | info[0].n_atoms, local_atoms); | |
1297 | painCave.isFatal = 1; | |
1298 | simError(); | |
# | Line 1349 | Line 1403 | void SimSetup::makeIntegrator(void){ | |
1403 | NVT<RealIntegrator>* myNVT = NULL; | |
1404 | NPTi<NPT<RealIntegrator> >* myNPTi = NULL; | |
1405 | NPTf<NPT<RealIntegrator> >* myNPTf = NULL; | |
1406 | + | NPTxyz<NPT<RealIntegrator> >* myNPTxyz = NULL; |
1407 | ||
1408 | for (k = 0; k < nInfo; k++){ | |
1409 | switch (ensembleCase){ | |
# | Line 1379 | Line 1434 | void SimSetup::makeIntegrator(void){ | |
1434 | else{ | |
1435 | sprintf(painCave.errMsg, | |
1436 | "SimSetup error: If you use the NVT\n" | |
1437 | < | " ensemble, you must set tauThermostat.\n"); |
1437 | > | "\tensemble, you must set tauThermostat.\n"); |
1438 | painCave.isFatal = 1; | |
1439 | simError(); | |
1440 | } | |
# | Line 1402 | Line 1457 | void SimSetup::makeIntegrator(void){ | |
1457 | else{ | |
1458 | sprintf(painCave.errMsg, | |
1459 | "SimSetup error: If you use a constant pressure\n" | |
1460 | < | " ensemble, you must set targetPressure in the BASS file.\n"); |
1460 | > | "\tensemble, you must set targetPressure in the BASS file.\n"); |
1461 | painCave.isFatal = 1; | |
1462 | simError(); | |
1463 | } | |
# | Line 1412 | Line 1467 | void SimSetup::makeIntegrator(void){ | |
1467 | else{ | |
1468 | sprintf(painCave.errMsg, | |
1469 | "SimSetup error: If you use an NPT\n" | |
1470 | < | " ensemble, you must set tauThermostat.\n"); |
1470 | > | "\tensemble, you must set tauThermostat.\n"); |
1471 | painCave.isFatal = 1; | |
1472 | simError(); | |
1473 | } | |
# | Line 1422 | Line 1477 | void SimSetup::makeIntegrator(void){ | |
1477 | else{ | |
1478 | sprintf(painCave.errMsg, | |
1479 | "SimSetup error: If you use an NPT\n" | |
1480 | < | " ensemble, you must set tauBarostat.\n"); |
1480 | > | "\tensemble, you must set tauBarostat.\n"); |
1481 | painCave.isFatal = 1; | |
1482 | simError(); | |
1483 | } | |
# | Line 1445 | Line 1500 | void SimSetup::makeIntegrator(void){ | |
1500 | else{ | |
1501 | sprintf(painCave.errMsg, | |
1502 | "SimSetup error: If you use a constant pressure\n" | |
1503 | < | " ensemble, you must set targetPressure in the BASS file.\n"); |
1503 | > | "\tensemble, you must set targetPressure in the BASS file.\n"); |
1504 | painCave.isFatal = 1; | |
1505 | simError(); | |
1506 | } | |
1507 | ||
1508 | if (globals->haveTauThermostat()) | |
1509 | myNPTf->setTauThermostat(globals->getTauThermostat()); | |
1510 | + | |
1511 | else{ | |
1512 | sprintf(painCave.errMsg, | |
1513 | "SimSetup error: If you use an NPT\n" | |
1514 | < | " ensemble, you must set tauThermostat.\n"); |
1514 | > | "\tensemble, you must set tauThermostat.\n"); |
1515 | painCave.isFatal = 1; | |
1516 | simError(); | |
1517 | } | |
1518 | ||
1519 | if (globals->haveTauBarostat()) | |
1520 | myNPTf->setTauBarostat(globals->getTauBarostat()); | |
1521 | + | |
1522 | else{ | |
1523 | sprintf(painCave.errMsg, | |
1524 | "SimSetup error: If you use an NPT\n" | |
1525 | < | " ensemble, you must set tauBarostat.\n"); |
1525 | > | "\tensemble, you must set tauBarostat.\n"); |
1526 | painCave.isFatal = 1; | |
1527 | simError(); | |
1528 | } | |
1529 | ||
1530 | info->the_integrator = myNPTf; | |
1531 | + | break; |
1532 | + | |
1533 | + | case NPTxyz_ENS: |
1534 | + | if (globals->haveZconstraints()){ |
1535 | + | setupZConstraint(info[k]); |
1536 | + | myNPTxyz = new ZConstraint<NPTxyz<NPT <RealIntegrator> > >(&(info[k]), the_ff); |
1537 | + | } |
1538 | + | else |
1539 | + | myNPTxyz = new NPTxyz<NPT <RealIntegrator> >(&(info[k]), the_ff); |
1540 | + | |
1541 | + | myNPTxyz->setTargetTemp(globals->getTargetTemp()); |
1542 | + | |
1543 | + | if (globals->haveTargetPressure()) |
1544 | + | myNPTxyz->setTargetPressure(globals->getTargetPressure()); |
1545 | + | else{ |
1546 | + | sprintf(painCave.errMsg, |
1547 | + | "SimSetup error: If you use a constant pressure\n" |
1548 | + | "\tensemble, you must set targetPressure in the BASS file.\n"); |
1549 | + | painCave.isFatal = 1; |
1550 | + | simError(); |
1551 | + | } |
1552 | + | |
1553 | + | if (globals->haveTauThermostat()) |
1554 | + | myNPTxyz->setTauThermostat(globals->getTauThermostat()); |
1555 | + | else{ |
1556 | + | sprintf(painCave.errMsg, |
1557 | + | "SimSetup error: If you use an NPT\n" |
1558 | + | "\tensemble, you must set tauThermostat.\n"); |
1559 | + | painCave.isFatal = 1; |
1560 | + | simError(); |
1561 | + | } |
1562 | + | |
1563 | + | if (globals->haveTauBarostat()) |
1564 | + | myNPTxyz->setTauBarostat(globals->getTauBarostat()); |
1565 | + | else{ |
1566 | + | sprintf(painCave.errMsg, |
1567 | + | "SimSetup error: If you use an NPT\n" |
1568 | + | "\tensemble, you must set tauBarostat.\n"); |
1569 | + | painCave.isFatal = 1; |
1570 | + | simError(); |
1571 | + | } |
1572 | + | |
1573 | + | info->the_integrator = myNPTxyz; |
1574 | break; | |
1575 | ||
1576 | default: | |
# | Line 1518 | Line 1618 | void SimSetup::setupZConstraint(SimInfo& theInfo){ | |
1618 | } | |
1619 | else{ | |
1620 | sprintf(painCave.errMsg, | |
1621 | < | "ZConstraint error: If you use an ZConstraint\n" |
1622 | < | " , you must set sample time.\n"); |
1621 | > | "ZConstraint error: If you use a ZConstraint,\n" |
1622 | > | "\tyou must set zconsTime.\n"); |
1623 | painCave.isFatal = 1; | |
1624 | simError(); | |
1625 | } | |
# | Line 1534 | Line 1634 | void SimSetup::setupZConstraint(SimInfo& theInfo){ | |
1634 | else{ | |
1635 | double defaultZConsTol = 0.01; | |
1636 | sprintf(painCave.errMsg, | |
1637 | < | "ZConstraint Waring: Tolerance for z-constraint methodl is not specified\n" |
1638 | < | " , default value %f is used.\n", |
1637 | > | "ZConstraint Warning: Tolerance for z-constraint method is not specified.\n" |
1638 | > | "\tOOPSE will use a default value of %f.\n" |
1639 | > | "\tTo set the tolerance, use the zconsTol variable.\n", |
1640 | defaultZConsTol); | |
1641 | painCave.isFatal = 0; | |
1642 | simError(); | |
# | Line 1553 | Line 1654 | void SimSetup::setupZConstraint(SimInfo& theInfo){ | |
1654 | } | |
1655 | else{ | |
1656 | sprintf(painCave.errMsg, | |
1657 | < | "ZConstraint Warning: User does not set force Subtraction policy, " |
1658 | < | "PolicyByMass is used\n"); |
1657 | > | "ZConstraint Warning: No force subtraction policy was set.\n" |
1658 | > | "\tOOPSE will use PolicyByMass.\n" |
1659 | > | "\tTo set the policy, use the zconsForcePolicy variable.\n"); |
1660 | painCave.isFatal = 0; | |
1661 | simError(); | |
1662 | zconsForcePolicy->setData("BYMASS"); | |
# | Line 1562 | Line 1664 | void SimSetup::setupZConstraint(SimInfo& theInfo){ | |
1664 | ||
1665 | theInfo.addProperty(zconsForcePolicy); | |
1666 | ||
1667 | + | //set zcons gap |
1668 | + | DoubleData* zconsGap = new DoubleData(); |
1669 | + | zconsGap->setID(ZCONSGAP_ID); |
1670 | + | |
1671 | + | if (globals->haveZConsGap()){ |
1672 | + | zconsGap->setData(globals->getZconsGap()); |
1673 | + | theInfo.addProperty(zconsGap); |
1674 | + | } |
1675 | + | |
1676 | + | //set zcons fixtime |
1677 | + | DoubleData* zconsFixtime = new DoubleData(); |
1678 | + | zconsFixtime->setID(ZCONSFIXTIME_ID); |
1679 | + | |
1680 | + | if (globals->haveZConsFixTime()){ |
1681 | + | zconsFixtime->setData(globals->getZconsFixtime()); |
1682 | + | theInfo.addProperty(zconsFixtime); |
1683 | + | } |
1684 | + | |
1685 | + | |
1686 | //Determine the name of ouput file and add it into SimInfo's property list | |
1687 | //Be careful, do not use inFileName, since it is a pointer which | |
1688 | //point to a string at master node, and slave nodes do not contain that string | |
# | Line 1598 | Line 1719 | void SimSetup::setupZConstraint(SimInfo& theInfo){ | |
1719 | //check the uniqueness of index | |
1720 | if(!zconsParaData->isIndexUnique()){ | |
1721 | sprintf(painCave.errMsg, | |
1722 | < | "ZConstraint Error: molIndex is not unique\n"); |
1722 | > | "ZConstraint Error: molIndex is not unique!\n"); |
1723 | painCave.isFatal = 1; | |
1724 | simError(); | |
1725 | } | |
# | Line 1609 | Line 1730 | void SimSetup::setupZConstraint(SimInfo& theInfo){ | |
1730 | //push data into siminfo, therefore, we can retrieve later | |
1731 | theInfo.addProperty(zconsParaData); | |
1732 | } | |
1733 | + | |
1734 | + | void SimSetup::makeMinimizer(){ |
1735 | + | |
1736 | + | OOPSEMinimizer* myOOPSEMinimizer; |
1737 | + | MinimizerParameterSet* param; |
1738 | + | char minimizerName[100]; |
1739 | + | |
1740 | + | for (int i = 0; i < nInfo; i++){ |
1741 | + | |
1742 | + | //prepare parameter set for minimizer |
1743 | + | param = new MinimizerParameterSet(); |
1744 | + | param->setDefaultParameter(); |
1745 | + | |
1746 | + | if (globals->haveMinimizer()){ |
1747 | + | param->setFTol(globals->getMinFTol()); |
1748 | + | } |
1749 | + | |
1750 | + | if (globals->haveMinGTol()){ |
1751 | + | param->setGTol(globals->getMinGTol()); |
1752 | + | } |
1753 | + | |
1754 | + | if (globals->haveMinMaxIter()){ |
1755 | + | param->setMaxIteration(globals->getMinMaxIter()); |
1756 | + | } |
1757 | + | |
1758 | + | if (globals->haveMinWriteFrq()){ |
1759 | + | param->setMaxIteration(globals->getMinMaxIter()); |
1760 | + | } |
1761 | + | |
1762 | + | if (globals->haveMinWriteFrq()){ |
1763 | + | param->setWriteFrq(globals->getMinWriteFrq()); |
1764 | + | } |
1765 | + | |
1766 | + | if (globals->haveMinStepSize()){ |
1767 | + | param->setStepSize(globals->getMinStepSize()); |
1768 | + | } |
1769 | + | |
1770 | + | if (globals->haveMinLSMaxIter()){ |
1771 | + | param->setLineSearchMaxIteration(globals->getMinLSMaxIter()); |
1772 | + | } |
1773 | + | |
1774 | + | if (globals->haveMinLSTol()){ |
1775 | + | param->setLineSearchTol(globals->getMinLSTol()); |
1776 | + | } |
1777 | + | |
1778 | + | strcpy(minimizerName, globals->getMinimizer()); |
1779 | + | |
1780 | + | if (!strcasecmp(minimizerName, "CG")){ |
1781 | + | myOOPSEMinimizer = new PRCGMinimizer(&(info[i]), the_ff, param); |
1782 | + | } |
1783 | + | else if (!strcasecmp(minimizerName, "SD")){ |
1784 | + | //myOOPSEMinimizer = MinimizerFactory.creatMinimizer("", &(info[i]), the_ff, param); |
1785 | + | myOOPSEMinimizer = new SDMinimizer(&(info[i]), the_ff, param); |
1786 | + | } |
1787 | + | else{ |
1788 | + | sprintf(painCave.errMsg, |
1789 | + | "SimSetup error: Unrecognized Minimizer, use Conjugate Gradient \n"); |
1790 | + | painCave.isFatal = 0; |
1791 | + | simError(); |
1792 | + | |
1793 | + | myOOPSEMinimizer = new PRCGMinimizer(&(info[i]), the_ff, param); |
1794 | + | } |
1795 | + | info[i].the_integrator = myOOPSEMinimizer; |
1796 | + | |
1797 | + | //store the minimizer into simInfo |
1798 | + | info[i].the_minimizer = myOOPSEMinimizer; |
1799 | + | info[i].has_minimizer = true; |
1800 | + | } |
1801 | + | |
1802 | + | } |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |