# | 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 | ||
14 | #ifdef IS_MPI | |
15 | #include "mpiBASS.h" | |
# | Line 22 | Line 22 | |
22 | #define NVT_ENS 1 | |
23 | #define NPTi_ENS 2 | |
24 | #define NPTf_ENS 3 | |
25 | < | #define NPTim_ENS 4 |
26 | < | #define NPTfm_ENS 5 |
25 | > | #define NPTxyz_ENS 4 |
26 | ||
28 | – | #define FF_DUFF 0 |
29 | – | #define FF_LJ 1 |
30 | – | #define FF_EAM 2 |
27 | ||
28 | + | #define FF_DUFF 0 |
29 | + | #define FF_LJ 1 |
30 | + | #define FF_EAM 2 |
31 | + | #define FF_H2O 3 |
32 | + | |
33 | using namespace std; | |
34 | ||
35 | + | /** |
36 | + | * Check whether dividend is divisble by divisor or not |
37 | + | */ |
38 | + | bool isDivisible(double dividend, double divisor){ |
39 | + | double tolerance = 0.000001; |
40 | + | double quotient; |
41 | + | double diff; |
42 | + | int intQuotient; |
43 | + | |
44 | + | quotient = dividend / divisor; |
45 | + | |
46 | + | if (quotient < 0) |
47 | + | quotient = -quotient; |
48 | + | |
49 | + | intQuotient = int (quotient + tolerance); |
50 | + | |
51 | + | diff = fabs(fabs(dividend) - intQuotient * fabs(divisor)); |
52 | + | |
53 | + | if (diff <= tolerance) |
54 | + | return true; |
55 | + | else |
56 | + | return false; |
57 | + | } |
58 | + | |
59 | SimSetup::SimSetup(){ | |
60 | + | |
61 | + | initSuspend = false; |
62 | isInfoArray = 0; | |
63 | nInfo = 1; | |
64 | ||
# | Line 54 | Line 81 | void SimSetup::setSimInfo(SimInfo* the_info, int theNi | |
81 | info = the_info; | |
82 | nInfo = theNinfo; | |
83 | isInfoArray = 1; | |
84 | + | initSuspend = true; |
85 | } | |
86 | ||
87 | ||
# | Line 92 | Line 120 | void SimSetup::createSim(void){ | |
120 | #endif // is_mpi | |
121 | ||
122 | void SimSetup::createSim(void){ | |
95 | – | int i, j, k, globalAtomIndex; |
123 | ||
124 | // gather all of the information from the Bass file | |
125 | ||
# | Line 108 | Line 135 | void SimSetup::createSim(void){ | |
135 | ||
136 | // initialize the system coordinates | |
137 | ||
138 | < | if (!isInfoArray){ |
138 | > | if ( !initSuspend ){ |
139 | initSystemCoords(); | |
140 | + | |
141 | + | if( !(globals->getUseInitTime()) ) |
142 | + | info[0].currentTime = 0.0; |
143 | } | |
144 | ||
145 | // make the output filenames | |
146 | ||
147 | makeOutNames(); | |
148 | ||
149 | < | // make the integrator |
150 | < | |
151 | < | makeIntegrator(); |
152 | < | |
149 | > | if (globals->haveMinimizer()) |
150 | > | // make minimizer |
151 | > | makeMinimizer(); |
152 | > | else |
153 | > | // make the integrator |
154 | > | makeIntegrator(); |
155 | > | |
156 | #ifdef IS_MPI | |
157 | mpiSim->mpiRefresh(); | |
158 | #endif | |
# | Line 131 | Line 164 | void SimSetup::makeMolecules(void){ | |
164 | ||
165 | ||
166 | void SimSetup::makeMolecules(void){ | |
167 | < | int k, l; |
167 | > | int k; |
168 | int i, j, exI, exJ, tempEx, stampID, atomOffset, excludeOffset; | |
169 | molInit molInfo; | |
170 | DirectionalAtom* dAtom; | |
# | Line 146 | Line 179 | void SimSetup::makeMolecules(void){ | |
179 | bend_set* theBends; | |
180 | torsion_set* theTorsions; | |
181 | ||
149 | – | |
182 | //init the forceField paramters | |
183 | ||
184 | the_ff->readParams(); | |
# | Line 154 | Line 186 | void SimSetup::makeMolecules(void){ | |
186 | ||
187 | // init the atoms | |
188 | ||
189 | + | double phi, theta, psi; |
190 | + | double sux, suy, suz; |
191 | + | double Axx, Axy, Axz, Ayx, Ayy, Ayz, Azx, Azy, Azz; |
192 | double ux, uy, uz, u, uSqr; | |
193 | ||
194 | for (k = 0; k < nInfo; k++){ | |
# | Line 190 | Line 225 | void SimSetup::makeMolecules(void){ | |
225 | info[k].n_oriented++; | |
226 | molInfo.myAtoms[j] = dAtom; | |
227 | ||
228 | < | ux = currentAtom->getOrntX(); |
229 | < | uy = currentAtom->getOrntY(); |
230 | < | uz = currentAtom->getOrntZ(); |
228 | > | // Directional Atoms have standard unit vectors which are oriented |
229 | > | // in space using the three Euler angles. We assume the standard |
230 | > | // unit vector was originally along the z axis below. |
231 | ||
232 | + | phi = currentAtom->getEulerPhi() * M_PI / 180.0; |
233 | + | theta = currentAtom->getEulerTheta() * M_PI / 180.0; |
234 | + | psi = currentAtom->getEulerPsi()* M_PI / 180.0; |
235 | + | |
236 | + | Axx = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi)); |
237 | + | Axy = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi)); |
238 | + | Axz = sin(theta) * sin(psi); |
239 | + | |
240 | + | Ayx = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi)); |
241 | + | Ayy = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi)); |
242 | + | Ayz = sin(theta) * cos(psi); |
243 | + | |
244 | + | Azx = sin(phi) * sin(theta); |
245 | + | Azy = -cos(phi) * sin(theta); |
246 | + | Azz = cos(theta); |
247 | + | |
248 | + | sux = 0.0; |
249 | + | suy = 0.0; |
250 | + | suz = 1.0; |
251 | + | |
252 | + | ux = (Axx * sux) + (Ayx * suy) + (Azx * suz); |
253 | + | uy = (Axy * sux) + (Ayy * suy) + (Azy * suz); |
254 | + | uz = (Axz * sux) + (Ayz * suy) + (Azz * suz); |
255 | + | |
256 | uSqr = (ux * ux) + (uy * uy) + (uz * uz); | |
257 | ||
258 | u = sqrt(uSqr); | |
# | Line 553 | Line 612 | void SimSetup::gatherInfo(void){ | |
612 | ||
613 | ||
614 | void SimSetup::gatherInfo(void){ | |
615 | < | int i, j, k; |
615 | > | int i; |
616 | ||
617 | ensembleCase = -1; | |
618 | ffCase = -1; | |
# | Line 581 | Line 640 | void SimSetup::gatherInfo(void){ | |
640 | else if (!strcasecmp(force_field, "EAM")){ | |
641 | ffCase = FF_EAM; | |
642 | } | |
643 | + | else if (!strcasecmp(force_field, "WATER")){ |
644 | + | ffCase = FF_H2O; |
645 | + | } |
646 | else{ | |
647 | sprintf(painCave.errMsg, "SimSetup Error. Unrecognized force field -> %s\n", | |
648 | force_field); | |
# | Line 604 | Line 666 | void SimSetup::gatherInfo(void){ | |
666 | else if (!strcasecmp(ensemble, "NPTf")){ | |
667 | ensembleCase = NPTf_ENS; | |
668 | } | |
669 | < | else if (!strcasecmp(ensemble, "NPTim")){ |
670 | < | ensembleCase = NPTim_ENS; |
669 | > | else if (!strcasecmp(ensemble, "NPTxyz")){ |
670 | > | ensembleCase = NPTxyz_ENS; |
671 | } | |
610 | – | else if (!strcasecmp(ensemble, "NPTfm")){ |
611 | – | ensembleCase = NPTfm_ENS; |
612 | – | } |
672 | else{ | |
673 | sprintf(painCave.errMsg, | |
674 | < | "SimSetup Warning. Unrecognized Ensemble -> %s, " |
675 | < | "reverting to NVE for this simulation.\n", |
674 | > | "SimSetup Warning. Unrecognized Ensemble -> %s \n" |
675 | > | "\treverting to NVE for this simulation.\n", |
676 | ensemble); | |
677 | painCave.isFatal = 0; | |
678 | simError(); | |
# | Line 645 | Line 704 | void SimSetup::gatherInfo(void){ | |
704 | if (!the_components[i]->haveNMol()){ | |
705 | // we have a problem | |
706 | sprintf(painCave.errMsg, | |
707 | < | "SimSetup Error. No global NMol or component NMol" |
708 | < | " given. Cannot calculate the number of atoms.\n"); |
707 | > | "SimSetup Error. No global NMol or component NMol given.\n" |
708 | > | "\tCannot calculate the number of atoms.\n"); |
709 | painCave.isFatal = 1; | |
710 | simError(); | |
711 | } | |
# | Line 664 | Line 723 | void SimSetup::gatherInfo(void){ | |
723 | " Please give nMol in the components.\n"); | |
724 | painCave.isFatal = 1; | |
725 | simError(); | |
726 | + | } |
727 | + | |
728 | + | //check whether sample time, status time, thermal time and reset time are divisble by dt |
729 | + | if (!isDivisible(globals->getSampleTime(), globals->getDt())){ |
730 | + | sprintf(painCave.errMsg, |
731 | + | "Sample time is not divisible by dt.\n" |
732 | + | "\tThis will result in samples that are not uniformly\n" |
733 | + | "\tdistributed in time. If this is a problem, change\n" |
734 | + | "\tyour sampleTime variable.\n"); |
735 | + | painCave.isFatal = 0; |
736 | + | simError(); |
737 | + | } |
738 | + | |
739 | + | if (globals->haveStatusTime() && !isDivisible(globals->getSampleTime(), globals->getDt())){ |
740 | + | sprintf(painCave.errMsg, |
741 | + | "Status time is not divisible by dt.\n" |
742 | + | "\tThis will result in status reports that are not uniformly\n" |
743 | + | "\tdistributed in time. If this is a problem, change \n" |
744 | + | "\tyour statusTime variable.\n"); |
745 | + | painCave.isFatal = 0; |
746 | + | simError(); |
747 | } | |
748 | ||
749 | + | if (globals->haveThermalTime() && !isDivisible(globals->getThermalTime(), globals->getDt())){ |
750 | + | sprintf(painCave.errMsg, |
751 | + | "Thermal time is not divisible by dt.\n" |
752 | + | "\tThis will result in thermalizations that are not uniformly\n" |
753 | + | "\tdistributed in time. If this is a problem, change \n" |
754 | + | "\tyour thermalTime variable.\n"); |
755 | + | painCave.isFatal = 0; |
756 | + | simError(); |
757 | + | } |
758 | + | |
759 | + | if (globals->haveResetTime() && !isDivisible(globals->getResetTime(), globals->getDt())){ |
760 | + | sprintf(painCave.errMsg, |
761 | + | "Reset time is not divisible by dt.\n" |
762 | + | "\tThis will result in integrator resets that are not uniformly\n" |
763 | + | "\tdistributed in time. If this is a problem, change\n" |
764 | + | "\tyour resetTime variable.\n"); |
765 | + | painCave.isFatal = 0; |
766 | + | simError(); |
767 | + | } |
768 | + | |
769 | // set the status, sample, and thermal kick times | |
770 | ||
771 | for (i = 0; i < nInfo; i++){ | |
# | Line 688 | Line 788 | void SimSetup::gatherInfo(void){ | |
788 | info[i].thermalTime = globals->getThermalTime(); | |
789 | } | |
790 | ||
791 | < | // check for the temperature set flag |
791 | > | info[i].resetIntegrator = 0; |
792 | > | if( globals->haveResetTime() ){ |
793 | > | info[i].resetTime = globals->getResetTime(); |
794 | > | info[i].resetIntegrator = 1; |
795 | > | } |
796 | ||
797 | + | // check for the temperature set flag |
798 | + | |
799 | if (globals->haveTempSet()) | |
800 | info[i].setTemp = globals->getTempSet(); | |
801 | ||
802 | < | // get some of the tricky things that may still be in the globals |
802 | > | // check for the extended State init |
803 | ||
804 | < | double boxVector[3]; |
805 | < | if (globals->haveBox()){ |
806 | < | 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 | < | } |
804 | > | info[i].useInitXSstate = globals->getUseInitXSstate(); |
805 | > | info[i].orthoTolerance = globals->getOrthoBoxTolerance(); |
806 | > | |
807 | } | |
808 | < | |
808 | > | |
809 | //setup seed for random number generator | |
810 | int seedValue; | |
811 | ||
# | Line 780 | Line 845 | void SimSetup::gatherInfo(void){ | |
845 | for (int i = 0; i < nInfo; i++){ | |
846 | info[i].setSeed(seedValue); | |
847 | } | |
848 | < | |
848 | > | |
849 | #ifdef IS_MPI | |
850 | < | strcpy(checkPointMsg, "Succesfully gathered all information from Bass\n"); |
850 | > | strcpy(checkPointMsg, "Successfully gathered all information from Bass\n"); |
851 | MPIcheckPoint(); | |
852 | #endif // is_mpi | |
853 | } | |
# | Line 815 | Line 880 | void SimSetup::finalInfoCheck(void){ | |
880 | ||
881 | if (!globals->haveECR()){ | |
882 | sprintf(painCave.errMsg, | |
883 | < | "SimSetup Warning: using default value of 1/2 the smallest " |
884 | < | "box length for the electrostaticCutoffRadius.\n" |
885 | < | "I hope you have a very fast processor!\n"); |
883 | > | "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n" |
884 | > | "\tOOPSE will use a default value of 15.0 angstroms" |
885 | > | "\tfor the electrostaticCutoffRadius.\n"); |
886 | painCave.isFatal = 0; | |
887 | simError(); | |
888 | < | 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; |
888 | > | theEcr = 15.0; |
889 | } | |
890 | else{ | |
891 | theEcr = globals->getECR(); | |
# | Line 834 | Line 893 | void SimSetup::finalInfoCheck(void){ | |
893 | ||
894 | if (!globals->haveEST()){ | |
895 | sprintf(painCave.errMsg, | |
896 | < | "SimSetup Warning: using default value of 0.05 * the " |
897 | < | "electrostaticCutoffRadius for the electrostaticSkinThickness\n"); |
896 | > | "SimSetup Warning: No value was set for electrostaticSkinThickness.\n" |
897 | > | "\tOOPSE will use a default value of\n" |
898 | > | "\t0.05 * electrostaticCutoffRadius\n" |
899 | > | "\tfor the electrostaticSkinThickness\n"); |
900 | painCave.isFatal = 0; | |
901 | simError(); | |
902 | theEst = 0.05 * theEcr; | |
# | Line 844 | Line 905 | void SimSetup::finalInfoCheck(void){ | |
905 | theEst = globals->getEST(); | |
906 | } | |
907 | ||
908 | < | info[i].setEcr(theEcr, theEst); |
908 | > | info[i].setDefaultEcr(theEcr, theEst); |
909 | ||
910 | if (!globals->haveDielectric()){ | |
911 | sprintf(painCave.errMsg, | |
912 | < | "SimSetup Error: You are trying to use Reaction Field without" |
913 | < | "setting a dielectric constant!\n"); |
912 | > | "SimSetup Error: No Dielectric constant was set.\n" |
913 | > | "\tYou are trying to use Reaction Field without" |
914 | > | "\tsetting a dielectric constant!\n"); |
915 | painCave.isFatal = 1; | |
916 | simError(); | |
917 | } | |
# | Line 859 | Line 921 | void SimSetup::finalInfoCheck(void){ | |
921 | if (usesDipoles){ | |
922 | if (!globals->haveECR()){ | |
923 | sprintf(painCave.errMsg, | |
924 | < | "SimSetup Warning: using default value of 1/2 the smallest " |
925 | < | "box length for the electrostaticCutoffRadius.\n" |
926 | < | "I hope you have a very fast processor!\n"); |
927 | < | painCave.isFatal = 0; |
928 | < | simError(); |
929 | < | 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; |
924 | > | "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n" |
925 | > | "\tOOPSE will use a default value of 15.0 angstroms" |
926 | > | "\tfor the electrostaticCutoffRadius.\n"); |
927 | > | painCave.isFatal = 0; |
928 | > | simError(); |
929 | > | theEcr = 15.0; |
930 | } | |
931 | else{ | |
932 | theEcr = globals->getECR(); | |
933 | } | |
934 | < | |
934 | > | |
935 | if (!globals->haveEST()){ | |
936 | sprintf(painCave.errMsg, | |
937 | < | "SimSetup Warning: using default value of 0.05 * the " |
938 | < | "electrostaticCutoffRadius for the " |
939 | < | "electrostaticSkinThickness\n"); |
937 | > | "SimSetup Warning: No value was set for electrostaticSkinThickness.\n" |
938 | > | "\tOOPSE will use a default value of\n" |
939 | > | "\t0.05 * electrostaticCutoffRadius\n" |
940 | > | "\tfor the electrostaticSkinThickness\n"); |
941 | painCave.isFatal = 0; | |
942 | simError(); | |
943 | theEst = 0.05 * theEcr; | |
# | Line 888 | Line 945 | void SimSetup::finalInfoCheck(void){ | |
945 | else{ | |
946 | theEst = globals->getEST(); | |
947 | } | |
948 | < | |
949 | < | info[i].setEcr(theEcr, theEst); |
948 | > | |
949 | > | info[i].setDefaultEcr(theEcr, theEst); |
950 | } | |
951 | } | |
952 | } | |
896 | – | |
953 | #ifdef IS_MPI | |
954 | strcpy(checkPointMsg, "post processing checks out"); | |
955 | MPIcheckPoint(); | |
956 | #endif // is_mpi | |
957 | } | |
958 | < | |
958 | > | |
959 | void SimSetup::initSystemCoords(void){ | |
960 | int i; | |
961 | ||
# | Line 916 | Line 972 | void SimSetup::initSystemCoords(void){ | |
972 | if (worldRank == 0){ | |
973 | #endif //is_mpi | |
974 | inName = globals->getInitialConfig(); | |
919 | – | double* tempDouble = new double[1000000]; |
975 | fileInit = new InitializeFromFile(inName); | |
976 | #ifdef IS_MPI | |
977 | } | |
# | Line 928 | Line 983 | void SimSetup::initSystemCoords(void){ | |
983 | delete fileInit; | |
984 | } | |
985 | else{ | |
986 | < | #ifdef IS_MPI |
932 | < | |
986 | > | |
987 | // no init from bass | |
988 | < | |
988 | > | |
989 | sprintf(painCave.errMsg, | |
990 | < | "Cannot intialize a parallel simulation without an initial configuration file.\n"); |
991 | < | painCave.isFatal; |
990 | > | "Cannot intialize a simulation without an initial configuration file.\n"); |
991 | > | painCave.isFatal = 1;; |
992 | simError(); | |
993 | < | |
940 | < | #else |
941 | < | |
942 | < | initFromBass(); |
943 | < | |
944 | < | |
945 | < | #endif |
993 | > | |
994 | } | |
995 | ||
996 | #ifdef IS_MPI | |
# | Line 1096 | Line 1144 | void SimSetup::createFF(void){ | |
1144 | the_ff = new EAM_FF(); | |
1145 | break; | |
1146 | ||
1147 | + | case FF_H2O: |
1148 | + | the_ff = new WATER(); |
1149 | + | break; |
1150 | + | |
1151 | default: | |
1152 | sprintf(painCave.errMsg, | |
1153 | "SimSetup Error. Unrecognized force field in case statement.\n"); | |
# | Line 1160 | Line 1212 | void SimSetup::calcSysValues(void){ | |
1212 | } | |
1213 | ||
1214 | void SimSetup::calcSysValues(void){ | |
1215 | < | int i, j, k; |
1215 | > | int i; |
1216 | ||
1217 | int* molMembershipArray; | |
1218 | ||
# | Line 1238 | Line 1290 | void SimSetup::mpiMolDivide(void){ | |
1290 | ||
1291 | if (local_atoms != info[0].n_atoms){ | |
1292 | sprintf(painCave.errMsg, | |
1293 | < | "SimSetup error: mpiSim's localAtom (%d) and SimSetup's" |
1294 | < | " localAtom (%d) are not equal.\n", |
1293 | > | "SimSetup error: mpiSim's localAtom (%d) and SimSetup's\n" |
1294 | > | "\tlocalAtom (%d) are not equal.\n", |
1295 | info[0].n_atoms, local_atoms); | |
1296 | painCave.isFatal = 1; | |
1297 | simError(); | |
# | Line 1259 | Line 1311 | void SimSetup::makeSysArrays(void){ | |
1311 | ||
1312 | ||
1313 | void SimSetup::makeSysArrays(void){ | |
1314 | < | int i, j, k, l; |
1314 | > | |
1315 | > | #ifndef IS_MPI |
1316 | > | int k, j; |
1317 | > | #endif // is_mpi |
1318 | > | int i, l; |
1319 | ||
1320 | Atom** the_atoms; | |
1321 | Molecule* the_molecules; | |
# | Line 1342 | Line 1398 | void SimSetup::makeIntegrator(void){ | |
1398 | void SimSetup::makeIntegrator(void){ | |
1399 | int k; | |
1400 | ||
1401 | + | NVE<RealIntegrator>* myNVE = NULL; |
1402 | NVT<RealIntegrator>* myNVT = NULL; | |
1403 | < | NPTi<RealIntegrator>* myNPTi = NULL; |
1404 | < | NPTf<RealIntegrator>* myNPTf = NULL; |
1405 | < | NPTim<RealIntegrator>* myNPTim = NULL; |
1349 | < | NPTfm<RealIntegrator>* myNPTfm = NULL; |
1403 | > | NPTi<NPT<RealIntegrator> >* myNPTi = NULL; |
1404 | > | NPTf<NPT<RealIntegrator> >* myNPTf = NULL; |
1405 | > | NPTxyz<NPT<RealIntegrator> >* myNPTxyz = NULL; |
1406 | ||
1407 | for (k = 0; k < nInfo; k++){ | |
1408 | switch (ensembleCase){ | |
1409 | case NVE_ENS: | |
1410 | if (globals->haveZconstraints()){ | |
1411 | setupZConstraint(info[k]); | |
1412 | < | new ZConstraint<NVE<RealIntegrator> >(&(info[k]), the_ff); |
1412 | > | myNVE = new ZConstraint<NVE<RealIntegrator> >(&(info[k]), the_ff); |
1413 | } | |
1414 | < | else |
1415 | < | new NVE<RealIntegrator>(&(info[k]), the_ff); |
1414 | > | else{ |
1415 | > | myNVE = new NVE<RealIntegrator>(&(info[k]), the_ff); |
1416 | > | } |
1417 | > | |
1418 | > | info->the_integrator = myNVE; |
1419 | break; | |
1420 | ||
1421 | case NVT_ENS: | |
# | Line 1374 | Line 1433 | void SimSetup::makeIntegrator(void){ | |
1433 | else{ | |
1434 | sprintf(painCave.errMsg, | |
1435 | "SimSetup error: If you use the NVT\n" | |
1436 | < | " ensemble, you must set tauThermostat.\n"); |
1436 | > | "\tensemble, you must set tauThermostat.\n"); |
1437 | painCave.isFatal = 1; | |
1438 | simError(); | |
1439 | } | |
1440 | + | |
1441 | + | info->the_integrator = myNVT; |
1442 | break; | |
1443 | ||
1444 | case NPTi_ENS: | |
1445 | if (globals->haveZconstraints()){ | |
1446 | setupZConstraint(info[k]); | |
1447 | < | myNPTi = new ZConstraint<NPTi<RealIntegrator> >(&(info[k]), the_ff); |
1447 | > | myNPTi = new ZConstraint<NPTi<NPT <RealIntegrator> > >(&(info[k]), the_ff); |
1448 | } | |
1449 | else | |
1450 | < | myNPTi = new NPTi<RealIntegrator>(&(info[k]), the_ff); |
1450 | > | myNPTi = new NPTi<NPT<RealIntegrator> >(&(info[k]), the_ff); |
1451 | ||
1452 | myNPTi->setTargetTemp(globals->getTargetTemp()); | |
1453 | ||
# | Line 1395 | Line 1456 | void SimSetup::makeIntegrator(void){ | |
1456 | else{ | |
1457 | sprintf(painCave.errMsg, | |
1458 | "SimSetup error: If you use a constant pressure\n" | |
1459 | < | " ensemble, you must set targetPressure in the BASS file.\n"); |
1459 | > | "\tensemble, you must set targetPressure in the BASS file.\n"); |
1460 | painCave.isFatal = 1; | |
1461 | simError(); | |
1462 | } | |
# | Line 1405 | Line 1466 | void SimSetup::makeIntegrator(void){ | |
1466 | else{ | |
1467 | sprintf(painCave.errMsg, | |
1468 | "SimSetup error: If you use an NPT\n" | |
1469 | < | " ensemble, you must set tauThermostat.\n"); |
1469 | > | "\tensemble, you must set tauThermostat.\n"); |
1470 | painCave.isFatal = 1; | |
1471 | simError(); | |
1472 | } | |
# | Line 1415 | Line 1476 | void SimSetup::makeIntegrator(void){ | |
1476 | else{ | |
1477 | sprintf(painCave.errMsg, | |
1478 | "SimSetup error: If you use an NPT\n" | |
1479 | < | " ensemble, you must set tauBarostat.\n"); |
1479 | > | "\tensemble, you must set tauBarostat.\n"); |
1480 | painCave.isFatal = 1; | |
1481 | simError(); | |
1482 | } | |
1483 | + | |
1484 | + | info->the_integrator = myNPTi; |
1485 | break; | |
1486 | ||
1487 | case NPTf_ENS: | |
1488 | if (globals->haveZconstraints()){ | |
1489 | setupZConstraint(info[k]); | |
1490 | < | myNPTf = new ZConstraint<NPTf<RealIntegrator> >(&(info[k]), the_ff); |
1490 | > | myNPTf = new ZConstraint<NPTf<NPT <RealIntegrator> > >(&(info[k]), the_ff); |
1491 | } | |
1492 | else | |
1493 | < | myNPTf = new NPTf<RealIntegrator>(&(info[k]), the_ff); |
1493 | > | myNPTf = new NPTf<NPT <RealIntegrator> >(&(info[k]), the_ff); |
1494 | ||
1495 | myNPTf->setTargetTemp(globals->getTargetTemp()); | |
1496 | ||
# | Line 1436 | Line 1499 | void SimSetup::makeIntegrator(void){ | |
1499 | else{ | |
1500 | sprintf(painCave.errMsg, | |
1501 | "SimSetup error: If you use a constant pressure\n" | |
1502 | < | " ensemble, you must set targetPressure in the BASS file.\n"); |
1502 | > | "\tensemble, you must set targetPressure in the BASS file.\n"); |
1503 | painCave.isFatal = 1; | |
1504 | simError(); | |
1505 | } | |
1506 | ||
1507 | if (globals->haveTauThermostat()) | |
1508 | myNPTf->setTauThermostat(globals->getTauThermostat()); | |
1509 | + | |
1510 | else{ | |
1511 | sprintf(painCave.errMsg, | |
1512 | "SimSetup error: If you use an NPT\n" | |
1513 | < | " ensemble, you must set tauThermostat.\n"); |
1513 | > | "\tensemble, you must set tauThermostat.\n"); |
1514 | painCave.isFatal = 1; | |
1515 | simError(); | |
1516 | } | |
1517 | ||
1518 | if (globals->haveTauBarostat()) | |
1519 | myNPTf->setTauBarostat(globals->getTauBarostat()); | |
1456 | – | else{ |
1457 | – | sprintf(painCave.errMsg, |
1458 | – | "SimSetup error: If you use an NPT\n" |
1459 | – | " ensemble, you must set tauBarostat.\n"); |
1460 | – | painCave.isFatal = 1; |
1461 | – | simError(); |
1462 | – | } |
1463 | – | break; |
1520 | ||
1465 | – | case NPTim_ENS: |
1466 | – | if (globals->haveZconstraints()){ |
1467 | – | setupZConstraint(info[k]); |
1468 | – | myNPTim = new ZConstraint<NPTim<RealIntegrator> >(&(info[k]), the_ff); |
1469 | – | } |
1470 | – | else |
1471 | – | myNPTim = new NPTim<RealIntegrator>(&(info[k]), the_ff); |
1472 | – | |
1473 | – | myNPTim->setTargetTemp(globals->getTargetTemp()); |
1474 | – | |
1475 | – | if (globals->haveTargetPressure()) |
1476 | – | myNPTim->setTargetPressure(globals->getTargetPressure()); |
1521 | else{ | |
1522 | sprintf(painCave.errMsg, | |
1479 | – | "SimSetup error: If you use a constant pressure\n" |
1480 | – | " ensemble, you must set targetPressure in the BASS file.\n"); |
1481 | – | painCave.isFatal = 1; |
1482 | – | simError(); |
1483 | – | } |
1484 | – | |
1485 | – | if (globals->haveTauThermostat()) |
1486 | – | myNPTim->setTauThermostat(globals->getTauThermostat()); |
1487 | – | else{ |
1488 | – | sprintf(painCave.errMsg, |
1523 | "SimSetup error: If you use an NPT\n" | |
1524 | < | " ensemble, you must set tauThermostat.\n"); |
1524 | > | "\tensemble, you must set tauBarostat.\n"); |
1525 | painCave.isFatal = 1; | |
1526 | simError(); | |
1527 | } | |
1528 | ||
1529 | < | if (globals->haveTauBarostat()) |
1496 | < | myNPTim->setTauBarostat(globals->getTauBarostat()); |
1497 | < | else{ |
1498 | < | sprintf(painCave.errMsg, |
1499 | < | "SimSetup error: If you use an NPT\n" |
1500 | < | " ensemble, you must set tauBarostat.\n"); |
1501 | < | painCave.isFatal = 1; |
1502 | < | simError(); |
1503 | < | } |
1529 | > | info->the_integrator = myNPTf; |
1530 | break; | |
1531 | ||
1532 | < | case NPTfm_ENS: |
1532 | > | case NPTxyz_ENS: |
1533 | if (globals->haveZconstraints()){ | |
1534 | setupZConstraint(info[k]); | |
1535 | < | myNPTfm = new ZConstraint<NPTfm<RealIntegrator> >(&(info[k]), the_ff); |
1535 | > | myNPTxyz = new ZConstraint<NPTxyz<NPT <RealIntegrator> > >(&(info[k]), the_ff); |
1536 | } | |
1537 | else | |
1538 | < | myNPTfm = new NPTfm<RealIntegrator>(&(info[k]), the_ff); |
1538 | > | myNPTxyz = new NPTxyz<NPT <RealIntegrator> >(&(info[k]), the_ff); |
1539 | ||
1540 | < | myNPTfm->setTargetTemp(globals->getTargetTemp()); |
1540 | > | myNPTxyz->setTargetTemp(globals->getTargetTemp()); |
1541 | ||
1542 | if (globals->haveTargetPressure()) | |
1543 | < | myNPTfm->setTargetPressure(globals->getTargetPressure()); |
1543 | > | myNPTxyz->setTargetPressure(globals->getTargetPressure()); |
1544 | else{ | |
1545 | sprintf(painCave.errMsg, | |
1546 | "SimSetup error: If you use a constant pressure\n" | |
1547 | < | " ensemble, you must set targetPressure in the BASS file.\n"); |
1547 | > | "\tensemble, you must set targetPressure in the BASS file.\n"); |
1548 | painCave.isFatal = 1; | |
1549 | simError(); | |
1550 | < | } |
1550 | > | } |
1551 | ||
1552 | if (globals->haveTauThermostat()) | |
1553 | < | myNPTfm->setTauThermostat(globals->getTauThermostat()); |
1553 | > | myNPTxyz->setTauThermostat(globals->getTauThermostat()); |
1554 | else{ | |
1555 | sprintf(painCave.errMsg, | |
1556 | "SimSetup error: If you use an NPT\n" | |
1557 | < | " ensemble, you must set tauThermostat.\n"); |
1557 | > | "\tensemble, you must set tauThermostat.\n"); |
1558 | painCave.isFatal = 1; | |
1559 | simError(); | |
1560 | } | |
1561 | ||
1562 | if (globals->haveTauBarostat()) | |
1563 | < | myNPTfm->setTauBarostat(globals->getTauBarostat()); |
1563 | > | myNPTxyz->setTauBarostat(globals->getTauBarostat()); |
1564 | else{ | |
1565 | sprintf(painCave.errMsg, | |
1566 | "SimSetup error: If you use an NPT\n" | |
1567 | < | " ensemble, you must set tauBarostat.\n"); |
1567 | > | "\tensemble, you must set tauBarostat.\n"); |
1568 | painCave.isFatal = 1; | |
1569 | simError(); | |
1570 | } | |
1571 | + | |
1572 | + | info->the_integrator = myNPTxyz; |
1573 | break; | |
1574 | ||
1575 | default: | |
# | Line 1589 | Line 1617 | void SimSetup::setupZConstraint(SimInfo& theInfo){ | |
1617 | } | |
1618 | else{ | |
1619 | sprintf(painCave.errMsg, | |
1620 | < | "ZConstraint error: If you use an ZConstraint\n" |
1621 | < | " , you must set sample time.\n"); |
1620 | > | "ZConstraint error: If you use a ZConstraint,\n" |
1621 | > | "\tyou must set zconsTime.\n"); |
1622 | painCave.isFatal = 1; | |
1623 | simError(); | |
1624 | } | |
# | Line 1605 | Line 1633 | void SimSetup::setupZConstraint(SimInfo& theInfo){ | |
1633 | else{ | |
1634 | double defaultZConsTol = 0.01; | |
1635 | sprintf(painCave.errMsg, | |
1636 | < | "ZConstraint Waring: Tolerance for z-constraint methodl is not specified\n" |
1637 | < | " , default value %f is used.\n", |
1636 | > | "ZConstraint Warning: Tolerance for z-constraint method is not specified.\n" |
1637 | > | "\tOOPSE will use a default value of %f.\n" |
1638 | > | "\tTo set the tolerance, use the zconsTol variable.\n", |
1639 | defaultZConsTol); | |
1640 | painCave.isFatal = 0; | |
1641 | simError(); | |
# | Line 1624 | Line 1653 | void SimSetup::setupZConstraint(SimInfo& theInfo){ | |
1653 | } | |
1654 | else{ | |
1655 | sprintf(painCave.errMsg, | |
1656 | < | "ZConstraint Warning: User does not set force Subtraction policy, " |
1657 | < | "PolicyByMass is used\n"); |
1656 | > | "ZConstraint Warning: No force subtraction policy was set.\n" |
1657 | > | "\tOOPSE will use PolicyByMass.\n" |
1658 | > | "\tTo set the policy, use the zconsForcePolicy variable.\n"); |
1659 | painCave.isFatal = 0; | |
1660 | simError(); | |
1661 | zconsForcePolicy->setData("BYMASS"); | |
# | Line 1669 | Line 1699 | void SimSetup::setupZConstraint(SimInfo& theInfo){ | |
1699 | //check the uniqueness of index | |
1700 | if(!zconsParaData->isIndexUnique()){ | |
1701 | sprintf(painCave.errMsg, | |
1702 | < | "ZConstraint Error: molIndex is not unique\n"); |
1702 | > | "ZConstraint Error: molIndex is not unique!\n"); |
1703 | painCave.isFatal = 1; | |
1704 | simError(); | |
1705 | } | |
# | Line 1680 | Line 1710 | void SimSetup::setupZConstraint(SimInfo& theInfo){ | |
1710 | //push data into siminfo, therefore, we can retrieve later | |
1711 | theInfo.addProperty(zconsParaData); | |
1712 | } | |
1713 | + | |
1714 | + | void SimSetup::makeMinimizer(){ |
1715 | + | /* |
1716 | + | OOPSEMinimizerBase* myOOPSEMinimizerBase; |
1717 | + | ObjFunctor1 * objFunc; |
1718 | + | OutputFunctor* outputFunc; |
1719 | + | ConcreteNLModel1* nlp; |
1720 | + | MinimizerParameterSet* param; |
1721 | + | ConjugateMinimizerBase* minimizer; |
1722 | + | int dim; |
1723 | + | |
1724 | + | for (int i = 0; i < nInfo; i++){ |
1725 | + | //creat |
1726 | + | myOOPSEMinimizerBase = new OOPSEMinimizerBase(&(info[i]), the_ff); |
1727 | + | |
1728 | + | //creat the object functor; |
1729 | + | objFunc = (ObjFunctor1*) new ClassMemObjFunctor1<OOPSEMinimizerBase> |
1730 | + | (myOOPSEMinimizerBase, &OOPSEMinimizerBase::calcGradient); |
1731 | + | |
1732 | + | //creat output functor; |
1733 | + | outputFunc = new ClassMemOutputFunctor<OOPSEMinimizerBase> |
1734 | + | (myOOPSEMinimizerBase, &OOPSEMinimizerBase::output); |
1735 | + | |
1736 | + | //creat nonlinear model |
1737 | + | dim = myOOPSEMinimizerBase->getDim(); |
1738 | + | nlp = new ConcreteNLModel1(dim, objFunc); |
1739 | + | |
1740 | + | //prepare parameter set for minimizer |
1741 | + | param = new MinimizerParameterSet(); |
1742 | + | param->setDefaultParameter(); |
1743 | + | |
1744 | + | if (globals->haveMinimizer()){ |
1745 | + | param->setFTol(globals->getMinFTol()); |
1746 | + | } |
1747 | + | |
1748 | + | if (globals->haveMinGTol()){ |
1749 | + | param->setGTol(globals->getMinGTol()); |
1750 | + | } |
1751 | + | |
1752 | + | if (globals->haveMinMaxIter()){ |
1753 | + | param->setMaxIteration(globals->getMinMaxIter()); |
1754 | + | } |
1755 | + | |
1756 | + | if (globals->haveMinWriteFrq()){ |
1757 | + | param->setMaxIteration(globals->getMinMaxIter()); |
1758 | + | } |
1759 | + | |
1760 | + | if (globals->haveMinWriteFrq()){ |
1761 | + | param->setWriteFrq(globals->getMinWriteFrq()); |
1762 | + | } |
1763 | + | |
1764 | + | if (globals->haveMinResetFrq()){ |
1765 | + | param->setResetFrq(globals->getMinResetFrq()); |
1766 | + | } |
1767 | + | |
1768 | + | if (globals->haveMinLSMaxIter()){ |
1769 | + | param->setLineSearchMaxIteration(globals->getMinLSMaxIter()); |
1770 | + | } |
1771 | + | |
1772 | + | if (globals->haveMinLSTol()){ |
1773 | + | param->setLineSearchTol(globals->getMinLSTol()); |
1774 | + | } |
1775 | + | |
1776 | + | //creat the minimizer |
1777 | + | minimizer = new PRCGMinimizer(nlp, param); |
1778 | + | minimizer->setLineSearchStrategy(nlp, GoldenSection); |
1779 | + | minimizer->setOutputFunctor(outputFunc); |
1780 | + | |
1781 | + | //store the minimizer into simInfo |
1782 | + | info[i].the_minimizer = minimizer; |
1783 | + | info[i].has_minimizer = true; |
1784 | + | } |
1785 | + | */ |
1786 | + | } |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |