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 733 by tim, Wed Aug 27 19:23:29 2003 UTC vs.
Revision 999 by chrisfen, Fri Jan 30 15:01:09 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"
# Line 22 | Line 21
21   #define NVT_ENS        1
22   #define NPTi_ENS       2
23   #define NPTf_ENS       3
24 < #define NPTim_ENS      4
26 < #define NPTfm_ENS      5
24 > #define NPTxyz_ENS     4
25  
28 – #define FF_DUFF 0
29 – #define FF_LJ   1
30 – #define FF_EAM  2
26  
27 + #define FF_DUFF  0
28 + #define FF_LJ    1
29 + #define FF_EAM   2
30 + #define FF_H2O 3
31 +
32   using namespace std;
33  
34 + /**
35 + * Check whether dividend is divisble by divisor or not
36 + */
37 + bool isDivisible(double dividend, double divisor){
38 +  double tolerance = 0.000001;
39 +  double quotient;
40 +  double diff;
41 +  int intQuotient;
42 +  
43 +  quotient = dividend / divisor;
44 +
45 +  if (quotient < 0)
46 +    quotient = -quotient;
47 +
48 +  intQuotient = int (quotient + tolerance);
49 +
50 +  diff = fabs(fabs(dividend) - intQuotient  * fabs(divisor));
51 +
52 +  if (diff <= tolerance)
53 +    return true;
54 +  else
55 +    return false;  
56 + }
57 +
58   SimSetup::SimSetup(){
59 +  
60 +  initSuspend = false;
61    isInfoArray = 0;
62    nInfo = 1;
63  
# Line 54 | Line 80 | void SimSetup::setSimInfo(SimInfo* the_info, int theNi
80    info = the_info;
81    nInfo = theNinfo;
82    isInfoArray = 1;
83 +  initSuspend = true;
84   }
85  
86  
# Line 92 | Line 119 | void SimSetup::createSim(void){
119   #endif // is_mpi
120  
121   void SimSetup::createSim(void){
95 –  int i, j, k, globalAtomIndex;
122  
123    // gather all of the information from the Bass file
124  
# Line 108 | Line 134 | void SimSetup::createSim(void){
134  
135    // initialize the system coordinates
136  
137 <  if (!isInfoArray){
137 >  if ( !initSuspend ){
138      initSystemCoords();
139 +
140 +    if( !(globals->getUseInitTime()) )
141 +      info[0].currentTime = 0.0;
142    }  
143  
144    // make the output filenames
# Line 131 | Line 160 | void SimSetup::makeMolecules(void){
160  
161  
162   void SimSetup::makeMolecules(void){
163 <  int k, l;
163 >  int k;
164    int i, j, exI, exJ, tempEx, stampID, atomOffset, excludeOffset;
165    molInit molInfo;
166    DirectionalAtom* dAtom;
# Line 146 | Line 175 | void SimSetup::makeMolecules(void){
175    bend_set* theBends;
176    torsion_set* theTorsions;
177  
149 –
178    //init the forceField paramters
179  
180    the_ff->readParams();
# Line 154 | Line 182 | void SimSetup::makeMolecules(void){
182  
183    // init the atoms
184  
185 +  double phi, theta, psi;
186 +  double sux, suy, suz;
187 +  double Axx, Axy, Axz, Ayx, Ayy, Ayz, Azx, Azy, Azz;
188    double ux, uy, uz, u, uSqr;
189  
190    for (k = 0; k < nInfo; k++){
# Line 190 | Line 221 | void SimSetup::makeMolecules(void){
221            info[k].n_oriented++;
222            molInfo.myAtoms[j] = dAtom;
223  
224 <          ux = currentAtom->getOrntX();
225 <          uy = currentAtom->getOrntY();
226 <          uz = currentAtom->getOrntZ();
224 >          // Directional Atoms have standard unit vectors which are oriented
225 >          // in space using the three Euler angles.  We assume the standard
226 >          // unit vector was originally along the z axis below.
227 >
228 >          phi = currentAtom->getEulerPhi() * M_PI / 180.0;
229 >          theta = currentAtom->getEulerTheta() * M_PI / 180.0;
230 >          psi = currentAtom->getEulerPsi()* M_PI / 180.0;
231 >            
232 >          Axx = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
233 >          Axy = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
234 >          Axz = sin(theta) * sin(psi);
235 >          
236 >          Ayx = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
237 >          Ayy = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
238 >          Ayz = sin(theta) * cos(psi);
239 >          
240 >          Azx = sin(phi) * sin(theta);
241 >          Azy = -cos(phi) * sin(theta);
242 >          Azz = cos(theta);
243  
244 +          sux = 0.0;
245 +          suy = 0.0;
246 +          suz = 1.0;
247 +
248 +          ux = (Axx * sux) + (Ayx * suy) + (Azx * suz);
249 +          uy = (Axy * sux) + (Ayy * suy) + (Azy * suz);
250 +          uz = (Axz * sux) + (Ayz * suy) + (Azz * suz);
251 +
252            uSqr = (ux * ux) + (uy * uy) + (uz * uz);
253  
254            u = sqrt(uSqr);
# Line 553 | Line 608 | void SimSetup::gatherInfo(void){
608  
609  
610   void SimSetup::gatherInfo(void){
611 <  int i, j, k;
611 >  int i;
612  
613    ensembleCase = -1;
614    ffCase = -1;
# Line 581 | Line 636 | void SimSetup::gatherInfo(void){
636    else if (!strcasecmp(force_field, "EAM")){
637      ffCase = FF_EAM;
638    }
639 +  else if (!strcasecmp(force_field, "WATER")){
640 +    ffCase = FF_H2O;
641 +  }
642    else{
643      sprintf(painCave.errMsg, "SimSetup Error. Unrecognized force field -> %s\n",
644              force_field);
# Line 604 | Line 662 | void SimSetup::gatherInfo(void){
662    else if (!strcasecmp(ensemble, "NPTf")){
663      ensembleCase = NPTf_ENS;
664    }
665 <  else if (!strcasecmp(ensemble, "NPTim")){
666 <    ensembleCase = NPTim_ENS;
665 >  else if (!strcasecmp(ensemble, "NPTxyz")){
666 >    ensembleCase = NPTxyz_ENS;
667    }
610 –  else if (!strcasecmp(ensemble, "NPTfm")){
611 –    ensembleCase = NPTfm_ENS;
612 –  }
668    else{
669      sprintf(painCave.errMsg,
670 <            "SimSetup Warning. Unrecognized Ensemble -> %s, "
671 <            "reverting to NVE for this simulation.\n",
670 >            "SimSetup Warning. Unrecognized Ensemble -> %s \n"
671 >            "\treverting to NVE for this simulation.\n",
672              ensemble);
673           painCave.isFatal = 0;
674           simError();
# Line 645 | Line 700 | void SimSetup::gatherInfo(void){
700        if (!the_components[i]->haveNMol()){
701          // we have a problem
702          sprintf(painCave.errMsg,
703 <                "SimSetup Error. No global NMol or component NMol"
704 <                " given. Cannot calculate the number of atoms.\n");
703 >                "SimSetup Error. No global NMol or component NMol given.\n"
704 >                "\tCannot calculate the number of atoms.\n");
705          painCave.isFatal = 1;
706          simError();
707        }
# Line 664 | Line 719 | void SimSetup::gatherInfo(void){
719              " Please give nMol in the components.\n");
720      painCave.isFatal = 1;
721      simError();
722 +  }
723 +
724 +  //check whether sample time, status time, thermal time and reset time are divisble by dt
725 +  if (!isDivisible(globals->getSampleTime(), globals->getDt())){
726 +    sprintf(painCave.errMsg,
727 +            "Sample time is not divisible by dt.\n"
728 +            "\tThis will result in samples that are not uniformly\n"
729 +            "\tdistributed in time.  If this is a problem, change\n"
730 +            "\tyour sampleTime variable.\n");
731 +    painCave.isFatal = 0;
732 +    simError();    
733    }
734  
735 +  if (globals->haveStatusTime() && !isDivisible(globals->getSampleTime(), globals->getDt())){
736 +    sprintf(painCave.errMsg,
737 +            "Status time is not divisible by dt.\n"
738 +            "\tThis will result in status reports that are not uniformly\n"
739 +            "\tdistributed in time.  If this is a problem, change \n"
740 +            "\tyour statusTime variable.\n");
741 +    painCave.isFatal = 0;
742 +    simError();    
743 +  }
744 +
745 +  if (globals->haveThermalTime() && !isDivisible(globals->getThermalTime(), globals->getDt())){
746 +    sprintf(painCave.errMsg,
747 +            "Thermal time is not divisible by dt.\n"
748 +            "\tThis will result in thermalizations that are not uniformly\n"
749 +            "\tdistributed in time.  If this is a problem, change \n"
750 +            "\tyour thermalTime variable.\n");
751 +    painCave.isFatal = 0;
752 +    simError();    
753 +  }  
754 +
755 +  if (globals->haveResetTime() && !isDivisible(globals->getResetTime(), globals->getDt())){
756 +    sprintf(painCave.errMsg,
757 +            "Reset time is not divisible by dt.\n"
758 +            "\tThis will result in integrator resets that are not uniformly\n"
759 +            "\tdistributed in time.  If this is a problem, change\n"
760 +            "\tyour resetTime variable.\n");
761 +    painCave.isFatal = 0;
762 +    simError();    
763 +  }
764 +
765    // set the status, sample, and thermal kick times
766  
767    for (i = 0; i < nInfo; i++){
# Line 688 | Line 784 | void SimSetup::gatherInfo(void){
784        info[i].thermalTime = globals->getThermalTime();
785      }
786  
787 <    // check for the temperature set flag
787 >    info[i].resetIntegrator = 0;
788 >    if( globals->haveResetTime() ){
789 >      info[i].resetTime = globals->getResetTime();
790 >      info[i].resetIntegrator = 1;
791 >    }
792  
793 +    // check for the temperature set flag
794 +    
795      if (globals->haveTempSet())
796        info[i].setTemp = globals->getTempSet();
797  
798 <    // get some of the tricky things that may still be in the globals
798 >    // check for the extended State init
799  
800 <    double boxVector[3];
801 <    if (globals->haveBox()){
802 <      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 <    }
800 >    info[i].useInitXSstate = globals->getUseInitXSstate();
801 >    info[i].orthoTolerance = globals->getOrthoBoxTolerance();
802 >    
803    }
804 <
804 >  
805    //setup seed for random number generator
806    int seedValue;
807  
# Line 782 | Line 843 | void SimSetup::gatherInfo(void){
843    }
844  
845   #ifdef IS_MPI
846 <  strcpy(checkPointMsg, "Succesfully gathered all information from Bass\n");
846 >  strcpy(checkPointMsg, "Successfully gathered all information from Bass\n");
847    MPIcheckPoint();
848   #endif // is_mpi
849   }
# Line 815 | Line 876 | void SimSetup::finalInfoCheck(void){
876  
877        if (!globals->haveECR()){
878          sprintf(painCave.errMsg,
879 <                "SimSetup Warning: using default value of 1/2 the smallest "
880 <                "box length for the electrostaticCutoffRadius.\n"
881 <                "I hope you have a very fast processor!\n");
879 >                "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n"
880 >                "\tOOPSE will use a default value of 15.0 angstroms"
881 >                "\tfor the electrostaticCutoffRadius.\n");
882          painCave.isFatal = 0;
883          simError();
884 <        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;
884 >        theEcr = 15.0;
885        }
886        else{
887          theEcr = globals->getECR();
# Line 834 | Line 889 | void SimSetup::finalInfoCheck(void){
889  
890        if (!globals->haveEST()){
891          sprintf(painCave.errMsg,
892 <                "SimSetup Warning: using default value of 0.05 * the "
893 <                "electrostaticCutoffRadius for the electrostaticSkinThickness\n");
892 >                "SimSetup Warning: No value was set for electrostaticSkinThickness.\n"
893 >                "\tOOPSE will use a default value of\n"
894 >                "\t0.05 * electrostaticCutoffRadius\n"
895 >                "\tfor the electrostaticSkinThickness\n");
896          painCave.isFatal = 0;
897          simError();
898          theEst = 0.05 * theEcr;
# Line 844 | Line 901 | void SimSetup::finalInfoCheck(void){
901          theEst = globals->getEST();
902        }
903  
904 <      info[i].setEcr(theEcr, theEst);
904 >      info[i].setDefaultEcr(theEcr, theEst);
905  
906        if (!globals->haveDielectric()){
907          sprintf(painCave.errMsg,
908 <                "SimSetup Error: You are trying to use Reaction Field without"
909 <                "setting a dielectric constant!\n");
908 >                "SimSetup Error: No Dielectric constant was set.\n"
909 >                "\tYou are trying to use Reaction Field without"
910 >                "\tsetting a dielectric constant!\n");
911          painCave.isFatal = 1;
912          simError();
913        }
# Line 859 | Line 917 | void SimSetup::finalInfoCheck(void){
917        if (usesDipoles){
918          if (!globals->haveECR()){
919            sprintf(painCave.errMsg,
920 <                  "SimSetup Warning: using default value of 1/2 the smallest "
921 <                  "box length for the electrostaticCutoffRadius.\n"
922 <                  "I hope you have a very fast processor!\n");
923 <          painCave.isFatal = 0;
924 <          simError();
925 <          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;
920 >                  "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n"
921 >                  "\tOOPSE will use a default value of 15.0 angstroms"
922 >                  "\tfor the electrostaticCutoffRadius.\n");
923 >          painCave.isFatal = 0;
924 >          simError();
925 >          theEcr = 15.0;
926          }
927          else{
928            theEcr = globals->getECR();
929          }
930 <
930 >        
931          if (!globals->haveEST()){
932            sprintf(painCave.errMsg,
933 <                  "SimSetup Warning: using default value of 0.05 * the "
934 <                  "electrostaticCutoffRadius for the "
935 <                  "electrostaticSkinThickness\n");
933 >                  "SimSetup Warning: No value was set for electrostaticSkinThickness.\n"
934 >                  "\tOOPSE will use a default value of\n"
935 >                  "\t0.05 * electrostaticCutoffRadius\n"
936 >                  "\tfor the electrostaticSkinThickness\n");
937            painCave.isFatal = 0;
938            simError();
939            theEst = 0.05 * theEcr;
# Line 888 | Line 941 | void SimSetup::finalInfoCheck(void){
941          else{
942            theEst = globals->getEST();
943          }
944 <
945 <        info[i].setEcr(theEcr, theEst);
944 >        
945 >        info[i].setDefaultEcr(theEcr, theEst);
946        }
947      }
948    }
896 –
949   #ifdef IS_MPI
950    strcpy(checkPointMsg, "post processing checks out");
951    MPIcheckPoint();
952   #endif // is_mpi
953   }
954 <
954 >  
955   void SimSetup::initSystemCoords(void){
956    int i;
957  
# Line 916 | Line 968 | void SimSetup::initSystemCoords(void){
968      if (worldRank == 0){
969   #endif //is_mpi
970        inName = globals->getInitialConfig();
919 –      double* tempDouble = new double[1000000];
971        fileInit = new InitializeFromFile(inName);
972   #ifdef IS_MPI
973      }
# Line 928 | Line 979 | void SimSetup::initSystemCoords(void){
979      delete fileInit;
980    }
981    else{
982 < #ifdef IS_MPI
932 <
982 >    
983      // no init from bass
984 <
984 >    
985      sprintf(painCave.errMsg,
986 <            "Cannot intialize a parallel simulation without an initial configuration file.\n");
987 <    painCave.isFatal;
986 >            "Cannot intialize a simulation without an initial configuration file.\n");
987 >    painCave.isFatal = 1;;
988      simError();
989 <
940 < #else
941 <
942 <    initFromBass();
943 <
944 <
945 < #endif
989 >    
990    }
991  
992   #ifdef IS_MPI
# Line 1096 | Line 1140 | void SimSetup::createFF(void){
1140        the_ff = new EAM_FF();
1141        break;
1142  
1143 +    case FF_H2O:
1144 +      the_ff = new WATER();
1145 +      break;
1146 +
1147      default:
1148        sprintf(painCave.errMsg,
1149                "SimSetup Error. Unrecognized force field in case statement.\n");
# Line 1160 | Line 1208 | void SimSetup::calcSysValues(void){
1208   }
1209  
1210   void SimSetup::calcSysValues(void){
1211 <  int i, j, k;
1211 >  int i;
1212  
1213    int* molMembershipArray;
1214  
# Line 1238 | Line 1286 | void SimSetup::mpiMolDivide(void){
1286  
1287    if (local_atoms != info[0].n_atoms){
1288      sprintf(painCave.errMsg,
1289 <            "SimSetup error: mpiSim's localAtom (%d) and SimSetup's"
1290 <            " localAtom (%d) are not equal.\n",
1289 >            "SimSetup error: mpiSim's localAtom (%d) and SimSetup's\n"
1290 >            "\tlocalAtom (%d) are not equal.\n",
1291              info[0].n_atoms, local_atoms);
1292      painCave.isFatal = 1;
1293      simError();
# Line 1259 | Line 1307 | void SimSetup::makeSysArrays(void){
1307  
1308  
1309   void SimSetup::makeSysArrays(void){
1310 <  int i, j, k, l;
1310 >
1311 > #ifndef IS_MPI
1312 >  int k, j;
1313 > #endif // is_mpi
1314 >  int i, l;
1315  
1316    Atom** the_atoms;
1317    Molecule* the_molecules;
# Line 1342 | Line 1394 | void SimSetup::makeIntegrator(void){
1394   void SimSetup::makeIntegrator(void){
1395    int k;
1396  
1397 +  NVE<RealIntegrator>* myNVE = NULL;
1398    NVT<RealIntegrator>* myNVT = NULL;
1399 <  NPTi<RealIntegrator>* myNPTi = NULL;
1400 <  NPTf<RealIntegrator>* myNPTf = NULL;
1401 <  NPTim<RealIntegrator>* myNPTim = NULL;
1349 <  NPTfm<RealIntegrator>* myNPTfm = NULL;
1399 >  NPTi<NPT<RealIntegrator> >* myNPTi = NULL;
1400 >  NPTf<NPT<RealIntegrator> >* myNPTf = NULL;
1401 >  NPTxyz<NPT<RealIntegrator> >* myNPTxyz = NULL;
1402    
1403    for (k = 0; k < nInfo; k++){
1404      switch (ensembleCase){
1405        case NVE_ENS:
1406          if (globals->haveZconstraints()){
1407            setupZConstraint(info[k]);
1408 <          new ZConstraint<NVE<RealIntegrator> >(&(info[k]), the_ff);
1408 >          myNVE = new ZConstraint<NVE<RealIntegrator> >(&(info[k]), the_ff);
1409          }
1410 <        else
1411 <          new NVE<RealIntegrator>(&(info[k]), the_ff);
1410 >        else{
1411 >          myNVE = new NVE<RealIntegrator>(&(info[k]), the_ff);
1412 >        }
1413 >        
1414 >        info->the_integrator = myNVE;
1415          break;
1416  
1417        case NVT_ENS:
# Line 1374 | Line 1429 | void SimSetup::makeIntegrator(void){
1429          else{
1430            sprintf(painCave.errMsg,
1431                    "SimSetup error: If you use the NVT\n"
1432 <                  "    ensemble, you must set tauThermostat.\n");
1432 >                  "\tensemble, you must set tauThermostat.\n");
1433            painCave.isFatal = 1;
1434            simError();
1435          }
1436 +
1437 +        info->the_integrator = myNVT;
1438          break;
1439  
1440        case NPTi_ENS:
1441          if (globals->haveZconstraints()){
1442            setupZConstraint(info[k]);
1443 <          myNPTi = new ZConstraint<NPTi<RealIntegrator> >(&(info[k]), the_ff);
1443 >          myNPTi = new ZConstraint<NPTi<NPT <RealIntegrator> > >(&(info[k]), the_ff);
1444          }
1445          else
1446 <          myNPTi = new NPTi<RealIntegrator>(&(info[k]), the_ff);
1446 >          myNPTi = new NPTi<NPT<RealIntegrator> >(&(info[k]), the_ff);
1447  
1448          myNPTi->setTargetTemp(globals->getTargetTemp());
1449  
# Line 1395 | Line 1452 | void SimSetup::makeIntegrator(void){
1452          else{
1453            sprintf(painCave.errMsg,
1454                    "SimSetup error: If you use a constant pressure\n"
1455 <                  "    ensemble, you must set targetPressure in the BASS file.\n");
1455 >                  "\tensemble, you must set targetPressure in the BASS file.\n");
1456            painCave.isFatal = 1;
1457            simError();
1458          }
# Line 1405 | Line 1462 | void SimSetup::makeIntegrator(void){
1462          else{
1463            sprintf(painCave.errMsg,
1464                    "SimSetup error: If you use an NPT\n"
1465 <                  "    ensemble, you must set tauThermostat.\n");
1465 >                  "\tensemble, you must set tauThermostat.\n");
1466            painCave.isFatal = 1;
1467            simError();
1468          }
# Line 1415 | Line 1472 | void SimSetup::makeIntegrator(void){
1472          else{
1473            sprintf(painCave.errMsg,
1474                    "SimSetup error: If you use an NPT\n"
1475 <                  "    ensemble, you must set tauBarostat.\n");
1475 >                  "\tensemble, you must set tauBarostat.\n");
1476            painCave.isFatal = 1;
1477            simError();
1478          }
1479 +
1480 +        info->the_integrator = myNPTi;
1481          break;
1482  
1483        case NPTf_ENS:
1484          if (globals->haveZconstraints()){
1485            setupZConstraint(info[k]);
1486 <          myNPTf = new ZConstraint<NPTf<RealIntegrator> >(&(info[k]), the_ff);
1486 >          myNPTf = new ZConstraint<NPTf<NPT <RealIntegrator> > >(&(info[k]), the_ff);
1487          }
1488          else
1489 <          myNPTf = new NPTf<RealIntegrator>(&(info[k]), the_ff);
1489 >          myNPTf = new NPTf<NPT <RealIntegrator> >(&(info[k]), the_ff);
1490  
1491          myNPTf->setTargetTemp(globals->getTargetTemp());
1492  
# Line 1436 | Line 1495 | void SimSetup::makeIntegrator(void){
1495          else{
1496            sprintf(painCave.errMsg,
1497                    "SimSetup error: If you use a constant pressure\n"
1498 <                  "    ensemble, you must set targetPressure in the BASS file.\n");
1498 >                  "\tensemble, you must set targetPressure in the BASS file.\n");
1499            painCave.isFatal = 1;
1500            simError();
1501          }    
1502  
1503          if (globals->haveTauThermostat())
1504            myNPTf->setTauThermostat(globals->getTauThermostat());
1505 +
1506          else{
1507            sprintf(painCave.errMsg,
1508                    "SimSetup error: If you use an NPT\n"
1509 <                  "    ensemble, you must set tauThermostat.\n");
1509 >                  "\tensemble, you must set tauThermostat.\n");
1510            painCave.isFatal = 1;
1511            simError();
1512          }
1513  
1514          if (globals->haveTauBarostat())
1515            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;
1516  
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());
1517          else{
1518            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,
1519                    "SimSetup error: If you use an NPT\n"
1520 <                  "    ensemble, you must set tauThermostat.\n");
1520 >                  "\tensemble, you must set tauBarostat.\n");
1521            painCave.isFatal = 1;
1522            simError();
1523          }
1524  
1525 <        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 <        }
1525 >        info->the_integrator = myNPTf;
1526          break;
1527  
1528 <      case NPTfm_ENS:
1528 >      case NPTxyz_ENS:
1529          if (globals->haveZconstraints()){
1530            setupZConstraint(info[k]);
1531 <          myNPTfm = new ZConstraint<NPTfm<RealIntegrator> >(&(info[k]), the_ff);
1531 >          myNPTxyz = new ZConstraint<NPTxyz<NPT <RealIntegrator> > >(&(info[k]), the_ff);
1532          }
1533          else
1534 <          myNPTfm = new NPTfm<RealIntegrator>(&(info[k]), the_ff);
1534 >          myNPTxyz = new NPTxyz<NPT <RealIntegrator> >(&(info[k]), the_ff);
1535  
1536 <        myNPTfm->setTargetTemp(globals->getTargetTemp());
1536 >        myNPTxyz->setTargetTemp(globals->getTargetTemp());
1537  
1538          if (globals->haveTargetPressure())
1539 <          myNPTfm->setTargetPressure(globals->getTargetPressure());
1539 >          myNPTxyz->setTargetPressure(globals->getTargetPressure());
1540          else{
1541            sprintf(painCave.errMsg,
1542                    "SimSetup error: If you use a constant pressure\n"
1543 <                  "    ensemble, you must set targetPressure in the BASS file.\n");
1543 >                  "\tensemble, you must set targetPressure in the BASS file.\n");
1544            painCave.isFatal = 1;
1545            simError();
1546 <        }
1546 >        }    
1547  
1548          if (globals->haveTauThermostat())
1549 <          myNPTfm->setTauThermostat(globals->getTauThermostat());
1549 >          myNPTxyz->setTauThermostat(globals->getTauThermostat());
1550          else{
1551            sprintf(painCave.errMsg,
1552                    "SimSetup error: If you use an NPT\n"
1553 <                  "    ensemble, you must set tauThermostat.\n");
1553 >                  "\tensemble, you must set tauThermostat.\n");
1554            painCave.isFatal = 1;
1555            simError();
1556          }
1557  
1558          if (globals->haveTauBarostat())
1559 <          myNPTfm->setTauBarostat(globals->getTauBarostat());
1559 >          myNPTxyz->setTauBarostat(globals->getTauBarostat());
1560          else{
1561            sprintf(painCave.errMsg,
1562                    "SimSetup error: If you use an NPT\n"
1563 <                  "    ensemble, you must set tauBarostat.\n");
1563 >                  "\tensemble, you must set tauBarostat.\n");
1564            painCave.isFatal = 1;
1565            simError();
1566          }
1567 +
1568 +        info->the_integrator = myNPTxyz;
1569          break;
1570  
1571        default:
# Line 1589 | Line 1613 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1613    }
1614    else{
1615      sprintf(painCave.errMsg,
1616 <            "ZConstraint error: If you use an ZConstraint\n"
1617 <            " , you must set sample time.\n");
1616 >            "ZConstraint error: If you use a ZConstraint,\n"
1617 >            "\tyou must set zconsTime.\n");
1618      painCave.isFatal = 1;
1619      simError();
1620    }
# Line 1605 | Line 1629 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1629    else{
1630      double defaultZConsTol = 0.01;
1631      sprintf(painCave.errMsg,
1632 <            "ZConstraint Waring: Tolerance for z-constraint methodl is not specified\n"
1633 <            " , default value %f is used.\n",
1632 >            "ZConstraint Warning: Tolerance for z-constraint method is not specified.\n"
1633 >            "\tOOPSE will use a default value of %f.\n"
1634 >            "\tTo set the tolerance, use the zconsTol variable.\n",
1635              defaultZConsTol);
1636      painCave.isFatal = 0;
1637      simError();      
# Line 1615 | Line 1640 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1640    }
1641    theInfo.addProperty(zconsTol);
1642  
1643 <  //set Force Substraction Policy
1643 >  //set Force Subtraction Policy
1644    StringData* zconsForcePolicy = new StringData();
1645    zconsForcePolicy->setID(ZCONSFORCEPOLICY_ID);
1646  
# Line 1624 | Line 1649 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1649    }
1650    else{
1651      sprintf(painCave.errMsg,
1652 <            "ZConstraint Warning: User does not set force substraction policy, "
1653 <            "average force substraction policy is used\n");
1652 >            "ZConstraint Warning: No force subtraction policy was set.\n"
1653 >            "\tOOPSE will use PolicyByMass.\n"
1654 >            "\tTo set the policy, use the zconsForcePolicy variable.\n");
1655      painCave.isFatal = 0;
1656      simError();
1657 <    zconsForcePolicy->setData("BYNUMBER");
1657 >    zconsForcePolicy->setData("BYMASS");
1658    }
1659  
1660    theInfo.addProperty(zconsForcePolicy);
# Line 1666 | Line 1692 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1692      zconsParaData->addItem(tempParaItem);
1693    }
1694  
1695 +  //check the uniqueness of index  
1696 +  if(!zconsParaData->isIndexUnique()){
1697 +    sprintf(painCave.errMsg,
1698 +            "ZConstraint Error: molIndex is not unique!\n");
1699 +    painCave.isFatal = 1;
1700 +    simError();
1701 +  }
1702 +
1703    //sort the parameters by index of molecules
1704    zconsParaData->sortByIndex();
1705 <
1705 >  
1706    //push data into siminfo, therefore, we can retrieve later
1707    theInfo.addProperty(zconsParaData);
1708   }

Diff Legend

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