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 984 by gezelter, Mon Jan 26 21:52:56 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  
26 +
27   #define FF_DUFF 0
28   #define FF_LJ   1
29   #define FF_EAM  2
30  
31   using namespace std;
32  
33 + /**
34 + * Check whether dividend is divisble by divisor or not
35 + */
36 + bool isDivisible(double dividend, double divisor){
37 +  double tolerance = 0.000001;
38 +  double quotient;
39 +  double diff;
40 +  int intQuotient;
41 +  
42 +  quotient = dividend / divisor;
43 +
44 +  if (quotient < 0)
45 +    quotient = -quotient;
46 +
47 +  intQuotient = int (quotient + tolerance);
48 +
49 +  diff = fabs(fabs(dividend) - intQuotient  * fabs(divisor));
50 +
51 +  if (diff <= tolerance)
52 +    return true;
53 +  else
54 +    return false;  
55 + }
56 +
57   SimSetup::SimSetup(){
58 +  
59 +  initSuspend = false;
60    isInfoArray = 0;
61    nInfo = 1;
62  
# Line 54 | Line 79 | void SimSetup::setSimInfo(SimInfo* the_info, int theNi
79    info = the_info;
80    nInfo = theNinfo;
81    isInfoArray = 1;
82 +  initSuspend = true;
83   }
84  
85  
# Line 92 | Line 118 | void SimSetup::createSim(void){
118   #endif // is_mpi
119  
120   void SimSetup::createSim(void){
95 –  int i, j, k, globalAtomIndex;
121  
122    // gather all of the information from the Bass file
123  
# Line 108 | Line 133 | void SimSetup::createSim(void){
133  
134    // initialize the system coordinates
135  
136 <  if (!isInfoArray){
136 >  if ( !initSuspend ){
137      initSystemCoords();
138 +
139 +    if( !(globals->getUseInitTime()) )
140 +      info[0].currentTime = 0.0;
141    }  
142  
143    // make the output filenames
# Line 131 | Line 159 | void SimSetup::makeMolecules(void){
159  
160  
161   void SimSetup::makeMolecules(void){
162 <  int k, l;
162 >  int k;
163    int i, j, exI, exJ, tempEx, stampID, atomOffset, excludeOffset;
164    molInit molInfo;
165    DirectionalAtom* dAtom;
# Line 146 | Line 174 | void SimSetup::makeMolecules(void){
174    bend_set* theBends;
175    torsion_set* theTorsions;
176  
149 –
177    //init the forceField paramters
178  
179    the_ff->readParams();
# Line 154 | Line 181 | void SimSetup::makeMolecules(void){
181  
182    // init the atoms
183  
184 +  double phi, theta, psi;
185 +  double sux, suy, suz;
186 +  double Axx, Axy, Axz, Ayx, Ayy, Ayz, Azx, Azy, Azz;
187    double ux, uy, uz, u, uSqr;
188  
189    for (k = 0; k < nInfo; k++){
# Line 190 | Line 220 | void SimSetup::makeMolecules(void){
220            info[k].n_oriented++;
221            molInfo.myAtoms[j] = dAtom;
222  
223 <          ux = currentAtom->getOrntX();
224 <          uy = currentAtom->getOrntY();
225 <          uz = currentAtom->getOrntZ();
223 >          // Directional Atoms have standard unit vectors which are oriented
224 >          // in space using the three Euler angles.  We assume the standard
225 >          // unit vector was originally along the z axis below.
226 >
227 >          phi = currentAtom->getEulerPhi() * M_PI / 180.0;
228 >          theta = currentAtom->getEulerTheta() * M_PI / 180.0;
229 >          psi = currentAtom->getEulerPsi()* M_PI / 180.0;
230 >            
231 >          Axx = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
232 >          Axy = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
233 >          Axz = sin(theta) * sin(psi);
234 >          
235 >          Ayx = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
236 >          Ayy = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
237 >          Ayz = sin(theta) * cos(psi);
238 >          
239 >          Azx = sin(phi) * sin(theta);
240 >          Azy = -cos(phi) * sin(theta);
241 >          Azz = cos(theta);
242  
243 +          sux = 0.0;
244 +          suy = 0.0;
245 +          suz = 1.0;
246 +
247 +          ux = (Axx * sux) + (Ayx * suy) + (Azx * suz);
248 +          uy = (Axy * sux) + (Ayy * suy) + (Azy * suz);
249 +          uz = (Axz * sux) + (Ayz * suy) + (Azz * suz);
250 +
251            uSqr = (ux * ux) + (uy * uy) + (uz * uz);
252  
253            u = sqrt(uSqr);
# Line 553 | Line 607 | void SimSetup::gatherInfo(void){
607  
608  
609   void SimSetup::gatherInfo(void){
610 <  int i, j, k;
610 >  int i;
611  
612    ensembleCase = -1;
613    ffCase = -1;
# Line 603 | Line 657 | void SimSetup::gatherInfo(void){
657    }
658    else if (!strcasecmp(ensemble, "NPTf")){
659      ensembleCase = NPTf_ENS;
606 –  }
607 –  else if (!strcasecmp(ensemble, "NPTim")){
608 –    ensembleCase = NPTim_ENS;
660    }
661 <  else if (!strcasecmp(ensemble, "NPTfm")){
662 <    ensembleCase = NPTfm_ENS;
661 >  else if (!strcasecmp(ensemble, "NPTxyz")){
662 >    ensembleCase = NPTxyz_ENS;
663    }
664    else{
665      sprintf(painCave.errMsg,
666 <            "SimSetup Warning. Unrecognized Ensemble -> %s, "
667 <            "reverting to NVE for this simulation.\n",
666 >            "SimSetup Warning. Unrecognized Ensemble -> %s \n"
667 >            "\treverting to NVE for this simulation.\n",
668              ensemble);
669           painCave.isFatal = 0;
670           simError();
# Line 645 | Line 696 | void SimSetup::gatherInfo(void){
696        if (!the_components[i]->haveNMol()){
697          // we have a problem
698          sprintf(painCave.errMsg,
699 <                "SimSetup Error. No global NMol or component NMol"
700 <                " given. Cannot calculate the number of atoms.\n");
699 >                "SimSetup Error. No global NMol or component NMol given.\n"
700 >                "\tCannot calculate the number of atoms.\n");
701          painCave.isFatal = 1;
702          simError();
703        }
# Line 666 | Line 717 | void SimSetup::gatherInfo(void){
717      simError();
718    }
719  
720 +  //check whether sample time, status time, thermal time and reset time are divisble by dt
721 +  if (!isDivisible(globals->getSampleTime(), globals->getDt())){
722 +    sprintf(painCave.errMsg,
723 +            "Sample time is not divisible by dt.\n"
724 +            "\tThis will result in samples that are not uniformly\n"
725 +            "\tdistributed in time.  If this is a problem, change\n"
726 +            "\tyour sampleTime variable.\n");
727 +    painCave.isFatal = 0;
728 +    simError();    
729 +  }
730 +
731 +  if (globals->haveStatusTime() && !isDivisible(globals->getSampleTime(), globals->getDt())){
732 +    sprintf(painCave.errMsg,
733 +            "Status time is not divisible by dt.\n"
734 +            "\tThis will result in status reports that are not uniformly\n"
735 +            "\tdistributed in time.  If this is a problem, change \n"
736 +            "\tyour statusTime variable.\n");
737 +    painCave.isFatal = 0;
738 +    simError();    
739 +  }
740 +
741 +  if (globals->haveThermalTime() && !isDivisible(globals->getThermalTime(), globals->getDt())){
742 +    sprintf(painCave.errMsg,
743 +            "Thermal time is not divisible by dt.\n"
744 +            "\tThis will result in thermalizations that are not uniformly\n"
745 +            "\tdistributed in time.  If this is a problem, change \n"
746 +            "\tyour thermalTime variable.\n");
747 +    painCave.isFatal = 0;
748 +    simError();    
749 +  }  
750 +
751 +  if (globals->haveResetTime() && !isDivisible(globals->getResetTime(), globals->getDt())){
752 +    sprintf(painCave.errMsg,
753 +            "Reset time is not divisible by dt.\n"
754 +            "\tThis will result in integrator resets that are not uniformly\n"
755 +            "\tdistributed in time.  If this is a problem, change\n"
756 +            "\tyour resetTime variable.\n");
757 +    painCave.isFatal = 0;
758 +    simError();    
759 +  }
760 +
761    // set the status, sample, and thermal kick times
762  
763    for (i = 0; i < nInfo; i++){
# Line 688 | Line 780 | void SimSetup::gatherInfo(void){
780        info[i].thermalTime = globals->getThermalTime();
781      }
782  
783 <    // check for the temperature set flag
783 >    info[i].resetIntegrator = 0;
784 >    if( globals->haveResetTime() ){
785 >      info[i].resetTime = globals->getResetTime();
786 >      info[i].resetIntegrator = 1;
787 >    }
788  
789 +    // check for the temperature set flag
790 +    
791      if (globals->haveTempSet())
792        info[i].setTemp = globals->getTempSet();
793  
794 <    // get some of the tricky things that may still be in the globals
794 >    // check for the extended State init
795  
796 <    double boxVector[3];
797 <    if (globals->haveBox()){
798 <      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 <    }
796 >    info[i].useInitXSstate = globals->getUseInitXSstate();
797 >    info[i].orthoTolerance = globals->getOrthoBoxTolerance();
798 >    
799    }
800 <
800 >  
801    //setup seed for random number generator
802    int seedValue;
803  
# Line 815 | Line 872 | void SimSetup::finalInfoCheck(void){
872  
873        if (!globals->haveECR()){
874          sprintf(painCave.errMsg,
875 <                "SimSetup Warning: using default value of 1/2 the smallest "
876 <                "box length for the electrostaticCutoffRadius.\n"
877 <                "I hope you have a very fast processor!\n");
875 >                "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n"
876 >                "\tOOPSE will use a default value of 15.0 angstroms"
877 >                "\tfor the electrostaticCutoffRadius.\n");
878          painCave.isFatal = 0;
879          simError();
880 <        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;
880 >        theEcr = 15.0;
881        }
882        else{
883          theEcr = globals->getECR();
# Line 834 | Line 885 | void SimSetup::finalInfoCheck(void){
885  
886        if (!globals->haveEST()){
887          sprintf(painCave.errMsg,
888 <                "SimSetup Warning: using default value of 0.05 * the "
889 <                "electrostaticCutoffRadius for the electrostaticSkinThickness\n");
888 >                "SimSetup Warning: No value was set for electrostaticSkinThickness.\n"
889 >                "\tOOPSE will use a default value of\n"
890 >                "\t0.05 * electrostaticCutoffRadius\n"
891 >                "\tfor the electrostaticSkinThickness\n");
892          painCave.isFatal = 0;
893          simError();
894          theEst = 0.05 * theEcr;
# Line 844 | Line 897 | void SimSetup::finalInfoCheck(void){
897          theEst = globals->getEST();
898        }
899  
900 <      info[i].setEcr(theEcr, theEst);
900 >      info[i].setDefaultEcr(theEcr, theEst);
901  
902        if (!globals->haveDielectric()){
903          sprintf(painCave.errMsg,
904 <                "SimSetup Error: You are trying to use Reaction Field without"
905 <                "setting a dielectric constant!\n");
904 >                "SimSetup Error: No Dielectric constant was set.\n"
905 >                "\tYou are trying to use Reaction Field without"
906 >                "\tsetting a dielectric constant!\n");
907          painCave.isFatal = 1;
908          simError();
909        }
# Line 859 | Line 913 | void SimSetup::finalInfoCheck(void){
913        if (usesDipoles){
914          if (!globals->haveECR()){
915            sprintf(painCave.errMsg,
916 <                  "SimSetup Warning: using default value of 1/2 the smallest "
917 <                  "box length for the electrostaticCutoffRadius.\n"
918 <                  "I hope you have a very fast processor!\n");
919 <          painCave.isFatal = 0;
920 <          simError();
921 <          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;
916 >                  "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n"
917 >                  "\tOOPSE will use a default value of 15.0 angstroms"
918 >                  "\tfor the electrostaticCutoffRadius.\n");
919 >          painCave.isFatal = 0;
920 >          simError();
921 >          theEcr = 15.0;
922          }
923          else{
924            theEcr = globals->getECR();
925          }
926 <
926 >        
927          if (!globals->haveEST()){
928            sprintf(painCave.errMsg,
929 <                  "SimSetup Warning: using default value of 0.05 * the "
930 <                  "electrostaticCutoffRadius for the "
931 <                  "electrostaticSkinThickness\n");
929 >                  "SimSetup Warning: No value was set for electrostaticSkinThickness.\n"
930 >                  "\tOOPSE will use a default value of\n"
931 >                  "\t0.05 * electrostaticCutoffRadius\n"
932 >                  "\tfor the electrostaticSkinThickness\n");
933            painCave.isFatal = 0;
934            simError();
935            theEst = 0.05 * theEcr;
# Line 888 | Line 937 | void SimSetup::finalInfoCheck(void){
937          else{
938            theEst = globals->getEST();
939          }
940 <
941 <        info[i].setEcr(theEcr, theEst);
940 >        
941 >        info[i].setDefaultEcr(theEcr, theEst);
942        }
943      }
944    }
896 –
945   #ifdef IS_MPI
946    strcpy(checkPointMsg, "post processing checks out");
947    MPIcheckPoint();
948   #endif // is_mpi
949   }
950 <
950 >  
951   void SimSetup::initSystemCoords(void){
952    int i;
953  
# Line 916 | Line 964 | void SimSetup::initSystemCoords(void){
964      if (worldRank == 0){
965   #endif //is_mpi
966        inName = globals->getInitialConfig();
919 –      double* tempDouble = new double[1000000];
967        fileInit = new InitializeFromFile(inName);
968   #ifdef IS_MPI
969      }
# Line 928 | Line 975 | void SimSetup::initSystemCoords(void){
975      delete fileInit;
976    }
977    else{
978 < #ifdef IS_MPI
932 <
978 >    
979      // no init from bass
980 <
980 >    
981      sprintf(painCave.errMsg,
982 <            "Cannot intialize a parallel simulation without an initial configuration file.\n");
983 <    painCave.isFatal;
982 >            "Cannot intialize a simulation without an initial configuration file.\n");
983 >    painCave.isFatal = 1;;
984      simError();
985 <
940 < #else
941 <
942 <    initFromBass();
943 <
944 <
945 < #endif
985 >    
986    }
987  
988   #ifdef IS_MPI
# Line 1160 | Line 1200 | void SimSetup::calcSysValues(void){
1200   }
1201  
1202   void SimSetup::calcSysValues(void){
1203 <  int i, j, k;
1203 >  int i;
1204  
1205    int* molMembershipArray;
1206  
# Line 1238 | Line 1278 | void SimSetup::mpiMolDivide(void){
1278  
1279    if (local_atoms != info[0].n_atoms){
1280      sprintf(painCave.errMsg,
1281 <            "SimSetup error: mpiSim's localAtom (%d) and SimSetup's"
1282 <            " localAtom (%d) are not equal.\n",
1281 >            "SimSetup error: mpiSim's localAtom (%d) and SimSetup's\n"
1282 >            "\tlocalAtom (%d) are not equal.\n",
1283              info[0].n_atoms, local_atoms);
1284      painCave.isFatal = 1;
1285      simError();
# Line 1259 | Line 1299 | void SimSetup::makeSysArrays(void){
1299  
1300  
1301   void SimSetup::makeSysArrays(void){
1302 <  int i, j, k, l;
1302 >
1303 > #ifndef IS_MPI
1304 >  int k, j;
1305 > #endif // is_mpi
1306 >  int i, l;
1307  
1308    Atom** the_atoms;
1309    Molecule* the_molecules;
# Line 1342 | Line 1386 | void SimSetup::makeIntegrator(void){
1386   void SimSetup::makeIntegrator(void){
1387    int k;
1388  
1389 +  NVE<RealIntegrator>* myNVE = NULL;
1390    NVT<RealIntegrator>* myNVT = NULL;
1391 <  NPTi<RealIntegrator>* myNPTi = NULL;
1392 <  NPTf<RealIntegrator>* myNPTf = NULL;
1393 <  NPTim<RealIntegrator>* myNPTim = NULL;
1394 <  NPTfm<RealIntegrator>* myNPTfm = NULL;
1350 <
1391 >  NPTi<NPT<RealIntegrator> >* myNPTi = NULL;
1392 >  NPTf<NPT<RealIntegrator> >* myNPTf = NULL;
1393 >  NPTxyz<NPT<RealIntegrator> >* myNPTxyz = NULL;
1394 >  
1395    for (k = 0; k < nInfo; k++){
1396      switch (ensembleCase){
1397        case NVE_ENS:
1398          if (globals->haveZconstraints()){
1399            setupZConstraint(info[k]);
1400 <          new ZConstraint<NVE<RealIntegrator> >(&(info[k]), the_ff);
1400 >          myNVE = new ZConstraint<NVE<RealIntegrator> >(&(info[k]), the_ff);
1401          }
1402 <        else
1403 <          new NVE<RealIntegrator>(&(info[k]), the_ff);
1402 >        else{
1403 >          myNVE = new NVE<RealIntegrator>(&(info[k]), the_ff);
1404 >        }
1405 >        
1406 >        info->the_integrator = myNVE;
1407          break;
1408  
1409        case NVT_ENS:
# Line 1374 | Line 1421 | void SimSetup::makeIntegrator(void){
1421          else{
1422            sprintf(painCave.errMsg,
1423                    "SimSetup error: If you use the NVT\n"
1424 <                  "    ensemble, you must set tauThermostat.\n");
1424 >                  "\tensemble, you must set tauThermostat.\n");
1425            painCave.isFatal = 1;
1426            simError();
1427          }
1428 +
1429 +        info->the_integrator = myNVT;
1430          break;
1431  
1432        case NPTi_ENS:
1433          if (globals->haveZconstraints()){
1434            setupZConstraint(info[k]);
1435 <          myNPTi = new ZConstraint<NPTi<RealIntegrator> >(&(info[k]), the_ff);
1435 >          myNPTi = new ZConstraint<NPTi<NPT <RealIntegrator> > >(&(info[k]), the_ff);
1436          }
1437          else
1438 <          myNPTi = new NPTi<RealIntegrator>(&(info[k]), the_ff);
1438 >          myNPTi = new NPTi<NPT<RealIntegrator> >(&(info[k]), the_ff);
1439  
1440          myNPTi->setTargetTemp(globals->getTargetTemp());
1441  
# Line 1395 | Line 1444 | void SimSetup::makeIntegrator(void){
1444          else{
1445            sprintf(painCave.errMsg,
1446                    "SimSetup error: If you use a constant pressure\n"
1447 <                  "    ensemble, you must set targetPressure in the BASS file.\n");
1447 >                  "\tensemble, you must set targetPressure in the BASS file.\n");
1448            painCave.isFatal = 1;
1449            simError();
1450          }
# Line 1405 | Line 1454 | void SimSetup::makeIntegrator(void){
1454          else{
1455            sprintf(painCave.errMsg,
1456                    "SimSetup error: If you use an NPT\n"
1457 <                  "    ensemble, you must set tauThermostat.\n");
1457 >                  "\tensemble, you must set tauThermostat.\n");
1458            painCave.isFatal = 1;
1459            simError();
1460          }
# Line 1415 | Line 1464 | void SimSetup::makeIntegrator(void){
1464          else{
1465            sprintf(painCave.errMsg,
1466                    "SimSetup error: If you use an NPT\n"
1467 <                  "    ensemble, you must set tauBarostat.\n");
1467 >                  "\tensemble, you must set tauBarostat.\n");
1468            painCave.isFatal = 1;
1469            simError();
1470          }
1471 +
1472 +        info->the_integrator = myNPTi;
1473          break;
1474  
1475        case NPTf_ENS:
1476          if (globals->haveZconstraints()){
1477            setupZConstraint(info[k]);
1478 <          myNPTf = new ZConstraint<NPTf<RealIntegrator> >(&(info[k]), the_ff);
1478 >          myNPTf = new ZConstraint<NPTf<NPT <RealIntegrator> > >(&(info[k]), the_ff);
1479          }
1480          else
1481 <          myNPTf = new NPTf<RealIntegrator>(&(info[k]), the_ff);
1481 >          myNPTf = new NPTf<NPT <RealIntegrator> >(&(info[k]), the_ff);
1482  
1483          myNPTf->setTargetTemp(globals->getTargetTemp());
1484  
# Line 1436 | Line 1487 | void SimSetup::makeIntegrator(void){
1487          else{
1488            sprintf(painCave.errMsg,
1489                    "SimSetup error: If you use a constant pressure\n"
1490 <                  "    ensemble, you must set targetPressure in the BASS file.\n");
1490 >                  "\tensemble, you must set targetPressure in the BASS file.\n");
1491            painCave.isFatal = 1;
1492            simError();
1493          }    
1494  
1495          if (globals->haveTauThermostat())
1496            myNPTf->setTauThermostat(globals->getTauThermostat());
1497 +
1498          else{
1499            sprintf(painCave.errMsg,
1500                    "SimSetup error: If you use an NPT\n"
1501 <                  "    ensemble, you must set tauThermostat.\n");
1501 >                  "\tensemble, you must set tauThermostat.\n");
1502            painCave.isFatal = 1;
1503            simError();
1504          }
1505  
1506          if (globals->haveTauBarostat())
1507            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;
1508  
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());
1509          else{
1510            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,
1511                    "SimSetup error: If you use an NPT\n"
1512 <                  "    ensemble, you must set tauThermostat.\n");
1512 >                  "\tensemble, you must set tauBarostat.\n");
1513            painCave.isFatal = 1;
1514            simError();
1515          }
1516  
1517 <        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 <        }
1517 >        info->the_integrator = myNPTf;
1518          break;
1519  
1520 <      case NPTfm_ENS:
1520 >      case NPTxyz_ENS:
1521          if (globals->haveZconstraints()){
1522            setupZConstraint(info[k]);
1523 <          myNPTfm = new ZConstraint<NPTfm<RealIntegrator> >(&(info[k]), the_ff);
1523 >          myNPTxyz = new ZConstraint<NPTxyz<NPT <RealIntegrator> > >(&(info[k]), the_ff);
1524          }
1525          else
1526 <          myNPTfm = new NPTfm<RealIntegrator>(&(info[k]), the_ff);
1526 >          myNPTxyz = new NPTxyz<NPT <RealIntegrator> >(&(info[k]), the_ff);
1527  
1528 <        myNPTfm->setTargetTemp(globals->getTargetTemp());
1528 >        myNPTxyz->setTargetTemp(globals->getTargetTemp());
1529  
1530          if (globals->haveTargetPressure())
1531 <          myNPTfm->setTargetPressure(globals->getTargetPressure());
1531 >          myNPTxyz->setTargetPressure(globals->getTargetPressure());
1532          else{
1533            sprintf(painCave.errMsg,
1534                    "SimSetup error: If you use a constant pressure\n"
1535 <                  "    ensemble, you must set targetPressure in the BASS file.\n");
1535 >                  "\tensemble, you must set targetPressure in the BASS file.\n");
1536            painCave.isFatal = 1;
1537            simError();
1538 <        }
1538 >        }    
1539  
1540          if (globals->haveTauThermostat())
1541 <          myNPTfm->setTauThermostat(globals->getTauThermostat());
1541 >          myNPTxyz->setTauThermostat(globals->getTauThermostat());
1542          else{
1543            sprintf(painCave.errMsg,
1544                    "SimSetup error: If you use an NPT\n"
1545 <                  "    ensemble, you must set tauThermostat.\n");
1545 >                  "\tensemble, you must set tauThermostat.\n");
1546            painCave.isFatal = 1;
1547            simError();
1548          }
1549  
1550          if (globals->haveTauBarostat())
1551 <          myNPTfm->setTauBarostat(globals->getTauBarostat());
1551 >          myNPTxyz->setTauBarostat(globals->getTauBarostat());
1552          else{
1553            sprintf(painCave.errMsg,
1554                    "SimSetup error: If you use an NPT\n"
1555 <                  "    ensemble, you must set tauBarostat.\n");
1555 >                  "\tensemble, you must set tauBarostat.\n");
1556            painCave.isFatal = 1;
1557            simError();
1558          }
1559 +
1560 +        info->the_integrator = myNPTxyz;
1561          break;
1562  
1563        default:
# Line 1589 | Line 1605 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1605    }
1606    else{
1607      sprintf(painCave.errMsg,
1608 <            "ZConstraint error: If you use an ZConstraint\n"
1609 <            " , you must set sample time.\n");
1608 >            "ZConstraint error: If you use a ZConstraint,\n"
1609 >            "\tyou must set zconsTime.\n");
1610      painCave.isFatal = 1;
1611      simError();
1612    }
# Line 1605 | Line 1621 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1621    else{
1622      double defaultZConsTol = 0.01;
1623      sprintf(painCave.errMsg,
1624 <            "ZConstraint Waring: Tolerance for z-constraint methodl is not specified\n"
1625 <            " , default value %f is used.\n",
1624 >            "ZConstraint Warning: Tolerance for z-constraint method is not specified.\n"
1625 >            "\tOOPSE will use a default value of %f.\n"
1626 >            "\tTo set the tolerance, use the zconsTol variable.\n",
1627              defaultZConsTol);
1628      painCave.isFatal = 0;
1629      simError();      
# Line 1615 | Line 1632 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1632    }
1633    theInfo.addProperty(zconsTol);
1634  
1635 <  //set Force Substraction Policy
1635 >  //set Force Subtraction Policy
1636    StringData* zconsForcePolicy = new StringData();
1637    zconsForcePolicy->setID(ZCONSFORCEPOLICY_ID);
1638  
# Line 1624 | Line 1641 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1641    }
1642    else{
1643      sprintf(painCave.errMsg,
1644 <            "ZConstraint Warning: User does not set force substraction policy, "
1645 <            "average force substraction policy is used\n");
1644 >            "ZConstraint Warning: No force subtraction policy was set.\n"
1645 >            "\tOOPSE will use PolicyByMass.\n"
1646 >            "\tTo set the policy, use the zconsForcePolicy variable.\n");
1647      painCave.isFatal = 0;
1648      simError();
1649 <    zconsForcePolicy->setData("BYNUMBER");
1649 >    zconsForcePolicy->setData("BYMASS");
1650    }
1651  
1652    theInfo.addProperty(zconsForcePolicy);
# Line 1666 | Line 1684 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1684      zconsParaData->addItem(tempParaItem);
1685    }
1686  
1687 +  //check the uniqueness of index  
1688 +  if(!zconsParaData->isIndexUnique()){
1689 +    sprintf(painCave.errMsg,
1690 +            "ZConstraint Error: molIndex is not unique!\n");
1691 +    painCave.isFatal = 1;
1692 +    simError();
1693 +  }
1694 +
1695    //sort the parameters by index of molecules
1696    zconsParaData->sortByIndex();
1697 <
1697 >  
1698    //push data into siminfo, therefore, we can retrieve later
1699    theInfo.addProperty(zconsParaData);
1700   }

Diff Legend

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