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 722 by tim, Tue Aug 26 20:02:29 2003 UTC vs.
Revision 1032 by tim, Fri Feb 6 19:05:07 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  
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  
# 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 666 | Line 725 | void SimSetup::gatherInfo(void){
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;
1406 <  NPTfm<RealIntegrator>* myNPTfm = NULL;
1350 <
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 1615 | Line 1644 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1644    }
1645    theInfo.addProperty(zconsTol);
1646  
1647 <  //set Force Substraction Policy
1647 >  //set Force Subtraction Policy
1648    StringData* zconsForcePolicy = new StringData();
1649    zconsForcePolicy->setID(ZCONSFORCEPOLICY_ID);
1650  
# Line 1624 | Line 1653 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1653    }
1654    else{
1655      sprintf(painCave.errMsg,
1656 <            "ZConstraint Warning: User does not set force substraction policy, "
1657 <            "average force substraction policy 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("BYNUMBER");
1661 >    zconsForcePolicy->setData("BYMASS");
1662    }
1663  
1664    theInfo.addProperty(zconsForcePolicy);
# Line 1666 | Line 1696 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1696      zconsParaData->addItem(tempParaItem);
1697    }
1698  
1699 +  //check the uniqueness of index  
1700 +  if(!zconsParaData->isIndexUnique()){
1701 +    sprintf(painCave.errMsg,
1702 +            "ZConstraint Error: molIndex is not unique!\n");
1703 +    painCave.isFatal = 1;
1704 +    simError();
1705 +  }
1706 +
1707    //sort the parameters by index of molecules
1708    zconsParaData->sortByIndex();
1709 <
1709 >  
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 + }

Diff Legend

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