# | 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 NPTim_ENS 4 |
26 | < | #define NPTfm_ENS 5 |
26 | > | #define NPTxyz_ENS 4 |
27 | ||
28 | – | #define FF_DUFF 0 |
29 | – | #define FF_LJ 1 |
30 | – | #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 54 | 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 92 | Line 121 | void SimSetup::createSim(void){ | |
121 | #endif // is_mpi | |
122 | ||
123 | void SimSetup::createSim(void){ | |
95 | – | int i, j, k, globalAtomIndex; |
124 | ||
125 | // gather all of the information from the Bass file | |
126 | ||
# | Line 108 | 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()) ) |
143 | + | info[0].currentTime = 0.0; |
144 | } | |
145 | ||
146 | // make the output filenames | |
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 131 | Line 165 | void SimSetup::makeMolecules(void){ | |
165 | ||
166 | ||
167 | void SimSetup::makeMolecules(void){ | |
168 | < | int k, l; |
168 | > | int k; |
169 | int i, j, exI, exJ, tempEx, stampID, atomOffset, excludeOffset; | |
170 | molInit molInfo; | |
171 | DirectionalAtom* dAtom; | |
# | 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 | ||
259 | u = sqrt(uSqr); | |
# | Line 553 | Line 613 | void SimSetup::gatherInfo(void){ | |
613 | ||
614 | ||
615 | void SimSetup::gatherInfo(void){ | |
616 | < | int i, j, k; |
616 | > | int i; |
617 | ||
618 | ensembleCase = -1; | |
619 | ffCase = -1; | |
# | 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, "NPTim")){ |
671 | < | ensembleCase = NPTim_ENS; |
670 | > | else if (!strcasecmp(ensemble, "NPTxyz")){ |
671 | > | ensembleCase = NPTxyz_ENS; |
672 | } | |
610 | – | else if (!strcasecmp(ensemble, "NPTfm")){ |
611 | – | ensembleCase = NPTfm_ENS; |
612 | – | } |
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 645 | 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 664 | 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 695 | 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(); |
707 | < | boxVector[1] = globals->getBox(); |
708 | < | boxVector[2] = globals->getBox(); |
709 | < | |
710 | < | info[i].setBox(boxVector); |
711 | < | } |
712 | < | else if (globals->haveDensity()){ |
713 | < | double vol; |
714 | < | vol = (double) tot_nmol / globals->getDensity(); |
715 | < | boxVector[0] = pow(vol, (1.0 / 3.0)); |
716 | < | boxVector[1] = boxVector[0]; |
717 | < | boxVector[2] = boxVector[0]; |
718 | < | |
719 | < | info[i].setBox(boxVector); |
720 | < | } |
721 | < | else{ |
722 | < | if (!globals->haveBoxX()){ |
723 | < | sprintf(painCave.errMsg, |
724 | < | "SimSetup error, no periodic BoxX size given.\n"); |
725 | < | painCave.isFatal = 1; |
726 | < | simError(); |
727 | < | } |
728 | < | boxVector[0] = globals->getBoxX(); |
729 | < | |
730 | < | if (!globals->haveBoxY()){ |
731 | < | sprintf(painCave.errMsg, |
732 | < | "SimSetup error, no periodic BoxY size given.\n"); |
733 | < | painCave.isFatal = 1; |
734 | < | simError(); |
735 | < | } |
736 | < | boxVector[1] = globals->getBoxY(); |
737 | < | |
738 | < | if (!globals->haveBoxZ()){ |
739 | < | sprintf(painCave.errMsg, |
740 | < | "SimSetup error, no periodic BoxZ size given.\n"); |
741 | < | painCave.isFatal = 1; |
742 | < | simError(); |
743 | < | } |
744 | < | boxVector[2] = globals->getBoxZ(); |
745 | < | |
746 | < | info[i].setBox(boxVector); |
747 | < | } |
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 786 | 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 821 | 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; |
830 | < | smallest = info[i].boxL[0]; |
831 | < | if (info[i].boxL[1] <= smallest) |
832 | < | smallest = info[i].boxL[1]; |
833 | < | if (info[i].boxL[2] <= smallest) |
834 | < | smallest = info[i].boxL[2]; |
835 | < | theEcr = 0.5 * smallest; |
889 | > | theEcr = 15.0; |
890 | } | |
891 | else{ | |
892 | theEcr = globals->getECR(); | |
# | Line 840 | 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 850 | 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 865 | 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; |
874 | < | smallest = info[i].boxL[0]; |
875 | < | if (info[i].boxL[1] <= smallest) |
876 | < | smallest = info[i].boxL[1]; |
877 | < | if (info[i].boxL[2] <= smallest) |
878 | < | smallest = info[i].boxL[2]; |
879 | < | 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 894 | 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 | } | |
902 | – | |
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 922 | Line 973 | void SimSetup::initSystemCoords(void){ | |
973 | if (worldRank == 0){ | |
974 | #endif //is_mpi | |
975 | inName = globals->getInitialConfig(); | |
925 | – | double* tempDouble = new double[1000000]; |
976 | fileInit = new InitializeFromFile(inName); | |
977 | #ifdef IS_MPI | |
978 | } | |
# | Line 934 | Line 984 | void SimSetup::initSystemCoords(void){ | |
984 | delete fileInit; | |
985 | } | |
986 | else{ | |
987 | < | #ifdef IS_MPI |
938 | < | |
987 | > | |
988 | // no init from bass | |
989 | < | |
989 | > | |
990 | sprintf(painCave.errMsg, | |
991 | < | "Cannot intialize a parallel simulation without an initial configuration file.\n"); |
992 | < | painCave.isFatal; |
991 | > | "Cannot intialize a simulation without an initial configuration file.\n"); |
992 | > | painCave.isFatal = 1;; |
993 | simError(); | |
994 | < | |
946 | < | #else |
947 | < | |
948 | < | initFromBass(); |
949 | < | |
950 | < | |
951 | < | #endif |
994 | > | |
995 | } | |
996 | ||
997 | #ifdef IS_MPI | |
# | Line 1100 | 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 1166 | Line 1213 | void SimSetup::calcSysValues(void){ | |
1213 | } | |
1214 | ||
1215 | void SimSetup::calcSysValues(void){ | |
1216 | < | int i, j, k; |
1216 | > | int i; |
1217 | ||
1218 | int* molMembershipArray; | |
1219 | ||
# | Line 1244 | 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 1265 | Line 1312 | void SimSetup::makeSysArrays(void){ | |
1312 | ||
1313 | ||
1314 | void SimSetup::makeSysArrays(void){ | |
1315 | < | int i, j, k, l; |
1315 | > | |
1316 | > | #ifndef IS_MPI |
1317 | > | int k, j; |
1318 | > | #endif // is_mpi |
1319 | > | int i, l; |
1320 | ||
1321 | Atom** the_atoms; | |
1322 | Molecule* the_molecules; | |
# | Line 1348 | Line 1399 | void SimSetup::makeIntegrator(void){ | |
1399 | void SimSetup::makeIntegrator(void){ | |
1400 | int k; | |
1401 | ||
1402 | + | NVE<RealIntegrator>* myNVE = NULL; |
1403 | NVT<RealIntegrator>* myNVT = NULL; | |
1404 | NPTi<NPT<RealIntegrator> >* myNPTi = NULL; | |
1405 | < | NPTf<RealIntegrator>* myNPTf = NULL; |
1406 | < | NPTim<RealIntegrator>* myNPTim = NULL; |
1355 | < | NPTfm<RealIntegrator>* myNPTfm = NULL; |
1405 | > | NPTf<NPT<RealIntegrator> >* myNPTf = NULL; |
1406 | > | NPTxyz<NPT<RealIntegrator> >* myNPTxyz = NULL; |
1407 | ||
1408 | for (k = 0; k < nInfo; k++){ | |
1409 | switch (ensembleCase){ | |
1410 | case NVE_ENS: | |
1411 | if (globals->haveZconstraints()){ | |
1412 | setupZConstraint(info[k]); | |
1413 | < | new ZConstraint<NVE<RealIntegrator> >(&(info[k]), the_ff); |
1413 | > | myNVE = new ZConstraint<NVE<RealIntegrator> >(&(info[k]), the_ff); |
1414 | } | |
1415 | < | else |
1416 | < | new NVE<RealIntegrator>(&(info[k]), the_ff); |
1415 | > | else{ |
1416 | > | myNVE = new NVE<RealIntegrator>(&(info[k]), the_ff); |
1417 | > | } |
1418 | > | |
1419 | > | info->the_integrator = myNVE; |
1420 | break; | |
1421 | ||
1422 | case NVT_ENS: | |
# | Line 1380 | 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 | } | |
1441 | + | |
1442 | + | info->the_integrator = myNVT; |
1443 | break; | |
1444 | ||
1445 | case NPTi_ENS: | |
# | Line 1401 | 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 1411 | 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 1421 | 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 | } | |
1484 | + | |
1485 | + | info->the_integrator = myNPTi; |
1486 | break; | |
1487 | ||
1488 | case NPTf_ENS: | |
1489 | if (globals->haveZconstraints()){ | |
1490 | setupZConstraint(info[k]); | |
1491 | < | myNPTf = new ZConstraint<NPTf<RealIntegrator> >(&(info[k]), the_ff); |
1491 | > | myNPTf = new ZConstraint<NPTf<NPT <RealIntegrator> > >(&(info[k]), the_ff); |
1492 | } | |
1493 | else | |
1494 | < | myNPTf = new NPTf<RealIntegrator>(&(info[k]), the_ff); |
1494 | > | myNPTf = new NPTf<NPT <RealIntegrator> >(&(info[k]), the_ff); |
1495 | ||
1496 | myNPTf->setTargetTemp(globals->getTargetTemp()); | |
1497 | ||
# | Line 1442 | 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()); | |
1462 | – | else{ |
1463 | – | sprintf(painCave.errMsg, |
1464 | – | "SimSetup error: If you use an NPT\n" |
1465 | – | " ensemble, you must set tauBarostat.\n"); |
1466 | – | painCave.isFatal = 1; |
1467 | – | simError(); |
1468 | – | } |
1469 | – | break; |
1521 | ||
1471 | – | case NPTim_ENS: |
1472 | – | if (globals->haveZconstraints()){ |
1473 | – | setupZConstraint(info[k]); |
1474 | – | myNPTim = new ZConstraint<NPTim<RealIntegrator> >(&(info[k]), the_ff); |
1475 | – | } |
1476 | – | else |
1477 | – | myNPTim = new NPTim<RealIntegrator>(&(info[k]), the_ff); |
1478 | – | |
1479 | – | myNPTim->setTargetTemp(globals->getTargetTemp()); |
1480 | – | |
1481 | – | if (globals->haveTargetPressure()) |
1482 | – | myNPTim->setTargetPressure(globals->getTargetPressure()); |
1522 | else{ | |
1523 | sprintf(painCave.errMsg, | |
1485 | – | "SimSetup error: If you use a constant pressure\n" |
1486 | – | " ensemble, you must set targetPressure in the BASS file.\n"); |
1487 | – | painCave.isFatal = 1; |
1488 | – | simError(); |
1489 | – | } |
1490 | – | |
1491 | – | if (globals->haveTauThermostat()) |
1492 | – | myNPTim->setTauThermostat(globals->getTauThermostat()); |
1493 | – | else{ |
1494 | – | sprintf(painCave.errMsg, |
1524 | "SimSetup error: If you use an NPT\n" | |
1525 | < | " ensemble, you must set tauThermostat.\n"); |
1525 | > | "\tensemble, you must set tauBarostat.\n"); |
1526 | painCave.isFatal = 1; | |
1527 | simError(); | |
1528 | } | |
1529 | ||
1530 | < | if (globals->haveTauBarostat()) |
1502 | < | myNPTim->setTauBarostat(globals->getTauBarostat()); |
1503 | < | else{ |
1504 | < | sprintf(painCave.errMsg, |
1505 | < | "SimSetup error: If you use an NPT\n" |
1506 | < | " ensemble, you must set tauBarostat.\n"); |
1507 | < | painCave.isFatal = 1; |
1508 | < | simError(); |
1509 | < | } |
1530 | > | info->the_integrator = myNPTf; |
1531 | break; | |
1532 | ||
1533 | < | case NPTfm_ENS: |
1533 | > | case NPTxyz_ENS: |
1534 | if (globals->haveZconstraints()){ | |
1535 | setupZConstraint(info[k]); | |
1536 | < | myNPTfm = new ZConstraint<NPTfm<RealIntegrator> >(&(info[k]), the_ff); |
1536 | > | myNPTxyz = new ZConstraint<NPTxyz<NPT <RealIntegrator> > >(&(info[k]), the_ff); |
1537 | } | |
1538 | else | |
1539 | < | myNPTfm = new NPTfm<RealIntegrator>(&(info[k]), the_ff); |
1539 | > | myNPTxyz = new NPTxyz<NPT <RealIntegrator> >(&(info[k]), the_ff); |
1540 | ||
1541 | < | myNPTfm->setTargetTemp(globals->getTargetTemp()); |
1541 | > | myNPTxyz->setTargetTemp(globals->getTargetTemp()); |
1542 | ||
1543 | if (globals->haveTargetPressure()) | |
1544 | < | myNPTfm->setTargetPressure(globals->getTargetPressure()); |
1544 | > | myNPTxyz->setTargetPressure(globals->getTargetPressure()); |
1545 | else{ | |
1546 | sprintf(painCave.errMsg, | |
1547 | "SimSetup error: If you use a constant pressure\n" | |
1548 | < | " ensemble, you must set targetPressure in the BASS file.\n"); |
1548 | > | "\tensemble, you must set targetPressure in the BASS file.\n"); |
1549 | painCave.isFatal = 1; | |
1550 | simError(); | |
1551 | < | } |
1551 | > | } |
1552 | ||
1553 | if (globals->haveTauThermostat()) | |
1554 | < | myNPTfm->setTauThermostat(globals->getTauThermostat()); |
1554 | > | myNPTxyz->setTauThermostat(globals->getTauThermostat()); |
1555 | else{ | |
1556 | sprintf(painCave.errMsg, | |
1557 | "SimSetup error: If you use an NPT\n" | |
1558 | < | " ensemble, you must set tauThermostat.\n"); |
1558 | > | "\tensemble, you must set tauThermostat.\n"); |
1559 | painCave.isFatal = 1; | |
1560 | simError(); | |
1561 | } | |
1562 | ||
1563 | if (globals->haveTauBarostat()) | |
1564 | < | myNPTfm->setTauBarostat(globals->getTauBarostat()); |
1564 | > | myNPTxyz->setTauBarostat(globals->getTauBarostat()); |
1565 | else{ | |
1566 | sprintf(painCave.errMsg, | |
1567 | "SimSetup error: If you use an NPT\n" | |
1568 | < | " ensemble, you must set tauBarostat.\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 1595 | 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 1611 | 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 1630 | 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 1639 | 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 1675 | 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 1686 | 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 |