ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/SimSetup.cpp
(Generate patch)

Comparing trunk/OOPSE/libmdtools/SimSetup.cpp (file contents):
Revision 736 by tim, Thu Aug 28 21:09:47 2003 UTC vs.
Revision 1066 by tim, Tue Feb 24 16:36:33 2004 UTC

# 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 666 | Line 726 | void SimSetup::gatherInfo(void){
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 688 | Line 789 | void SimSetup::gatherInfo(void){
789        info[i].thermalTime = globals->getThermalTime();
790      }
791  
792 <    // check for the temperature set flag
792 >    info[i].resetIntegrator = 0;
793 >    if( globals->haveResetTime() ){
794 >      info[i].resetTime = globals->getResetTime();
795 >      info[i].resetIntegrator = 1;
796 >    }
797  
798 +    // check for the temperature set flag
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 916 | Line 973 | void SimSetup::initSystemCoords(void){
973      if (worldRank == 0){
974   #endif //is_mpi
975        inName = globals->getInitialConfig();
919 –      double* tempDouble = new double[1000000];
976        fileInit = new InitializeFromFile(inName);
977   #ifdef IS_MPI
978      }
# Line 928 | Line 984 | void SimSetup::initSystemCoords(void){
984      delete fileInit;
985    }
986    else{
987 < #ifdef IS_MPI
932 <
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 <
940 < #else
941 <
942 <    initFromBass();
943 <
944 <
945 < #endif
994 >    
995    }
996  
997   #ifdef IS_MPI
# Line 1096 | Line 1145 | void SimSetup::createFF(void){
1145        the_ff = new EAM_FF();
1146        break;
1147  
1148 +    case FF_H2O:
1149 +      the_ff = new WATER();
1150 +      break;
1151 +
1152      default:
1153        sprintf(painCave.errMsg,
1154                "SimSetup Error. Unrecognized force field in case statement.\n");
# Line 1160 | 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 1238 | 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 1259 | 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 1342 | 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<RealIntegrator>* myNPTi = NULL;
1405 <  NPTf<RealIntegrator>* myNPTf = NULL;
1406 <  NPTim<RealIntegrator>* myNPTim = NULL;
1349 <  NPTfm<RealIntegrator>* myNPTfm = 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){
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 1374 | 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:
1446          if (globals->haveZconstraints()){
1447            setupZConstraint(info[k]);
1448 <          myNPTi = new ZConstraint<NPTi<RealIntegrator> >(&(info[k]), the_ff);
1448 >          myNPTi = new ZConstraint<NPTi<NPT <RealIntegrator> > >(&(info[k]), the_ff);
1449          }
1450          else
1451 <          myNPTi = new NPTi<RealIntegrator>(&(info[k]), the_ff);
1451 >          myNPTi = new NPTi<NPT<RealIntegrator> >(&(info[k]), the_ff);
1452  
1453          myNPTi->setTargetTemp(globals->getTargetTemp());
1454  
# Line 1395 | 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 1405 | 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 1415 | 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 1436 | 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());
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;
1521  
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());
1522          else{
1523            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,
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())
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 <        }
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 1589 | 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 1605 | 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 1615 | Line 1645 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1645    }
1646    theInfo.addProperty(zconsTol);
1647  
1648 <  //set Force Substraction Policy
1648 >  //set Force Subtraction Policy
1649    StringData* zconsForcePolicy = new StringData();
1650    zconsForcePolicy->setID(ZCONSFORCEPOLICY_ID);
1651  
# Line 1624 | Line 1654 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1654    }
1655    else{
1656      sprintf(painCave.errMsg,
1657 <            "ZConstraint Warning: User does not set force substraction 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 1669 | Line 1700 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1700    //check the uniqueness of index  
1701    if(!zconsParaData->isIndexUnique()){
1702      sprintf(painCave.errMsg,
1703 <            "ZConstraint Error: molIndex is not unique\n");
1703 >            "ZConstraint Error: molIndex is not unique!\n");
1704      painCave.isFatal = 1;
1705      simError();
1706    }
# Line 1680 | Line 1711 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1711    //push data into siminfo, therefore, we can retrieve later
1712    theInfo.addProperty(zconsParaData);
1713   }
1714 +
1715 + void SimSetup::makeMinimizer(){
1716 +
1717 +  OOPSEMinimizer* myOOPSEMinimizer;
1718 +  MinimizerParameterSet* param;
1719 +  char minimizerName[100];
1720 +  
1721 +  for (int i = 0; i < nInfo; i++){
1722 +    
1723 +    //prepare parameter set for minimizer
1724 +    param = new MinimizerParameterSet();
1725 +    param->setDefaultParameter();
1726 +
1727 +    if (globals->haveMinimizer()){
1728 +      param->setFTol(globals->getMinFTol());
1729 +    }
1730 +
1731 +    if (globals->haveMinGTol()){
1732 +      param->setGTol(globals->getMinGTol());
1733 +    }
1734 +
1735 +    if (globals->haveMinMaxIter()){
1736 +      param->setMaxIteration(globals->getMinMaxIter());
1737 +    }
1738 +
1739 +    if (globals->haveMinWriteFrq()){
1740 +      param->setMaxIteration(globals->getMinMaxIter());
1741 +    }
1742 +
1743 +    if (globals->haveMinWriteFrq()){
1744 +      param->setWriteFrq(globals->getMinWriteFrq());
1745 +    }
1746 +    
1747 +    if (globals->haveMinStepSize()){
1748 +      param->setStepSize(globals->getMinStepSize());
1749 +    }
1750 +
1751 +    if (globals->haveMinLSMaxIter()){
1752 +      param->setLineSearchMaxIteration(globals->getMinLSMaxIter());
1753 +    }    
1754 +
1755 +    if (globals->haveMinLSTol()){
1756 +      param->setLineSearchTol(globals->getMinLSTol());
1757 +    }    
1758 +
1759 +    strcpy(minimizerName, globals->getMinimizer());
1760 +
1761 +    if (!strcasecmp(minimizerName, "CG")){
1762 +      myOOPSEMinimizer = new PRCGMinimizer(&(info[i]), the_ff, param);
1763 +    }
1764 +    else if (!strcasecmp(minimizerName, "SD")){
1765 +    //myOOPSEMinimizer = MinimizerFactory.creatMinimizer("", &(info[i]), the_ff, param);
1766 +      myOOPSEMinimizer = new SDMinimizer(&(info[i]), the_ff, param);
1767 +    }
1768 +    else{
1769 +          sprintf(painCave.errMsg,
1770 +                  "SimSetup error: Unrecognized Minimizer, use Conjugate Gradient \n");
1771 +          painCave.isFatal = 0;
1772 +          simError();
1773 +
1774 +      myOOPSEMinimizer = new PRCGMinimizer(&(info[i]), the_ff, param);          
1775 +    }
1776 +     info[i].the_integrator = myOOPSEMinimizer;
1777 +
1778 +     //store the minimizer into simInfo
1779 +     info[i].the_minimizer = myOOPSEMinimizer;
1780 +     info[i].has_minimizer = true;
1781 +  }
1782 +
1783 + }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines