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 841 by mmeineke, Wed Oct 29 17:55:28 2003 UTC vs.
Revision 984 by gezelter, Mon Jan 26 21:52:56 2004 UTC

# Line 30 | Line 30 | using namespace std;
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;
# Line 103 | Line 127 | void SimSetup::createSim(void){
127  
128    sysObjectsCreation();
129  
130 +  // check on the post processing info
131 +
132 +  finalInfoCheck();
133 +
134    // initialize the system coordinates
135  
136    if ( !initSuspend ){
# Line 112 | Line 140 | void SimSetup::createSim(void){
140        info[0].currentTime = 0.0;
141    }  
142  
115 –  // check on the post processing info
116 –
117 –  finalInfoCheck();
118 –
143    // make the output filenames
144  
145    makeOutNames();
# Line 150 | Line 174 | void SimSetup::makeMolecules(void){
174    bend_set* theBends;
175    torsion_set* theTorsions;
176  
153 –
177    //init the forceField paramters
178  
179    the_ff->readParams();
# Line 158 | 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 194 | 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 613 | Line 663 | void SimSetup::gatherInfo(void){
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 646 | 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 667 | 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 699 | Line 790 | void SimSetup::gatherInfo(void){
790      
791      if (globals->haveTempSet())
792        info[i].setTemp = globals->getTempSet();
793 +
794 +    // check for the extended State init
795 +
796 +    info[i].useInitXSstate = globals->getUseInitXSstate();
797 +    info[i].orthoTolerance = globals->getOrthoBoxTolerance();
798      
799    }
800    
# Line 776 | 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;
785 <        smallest = info[i].boxL[0];
786 <        if (info[i].boxL[1] <= smallest)
787 <          smallest = info[i].boxL[1];
788 <        if (info[i].boxL[2] <= smallest)
789 <          smallest = info[i].boxL[2];
790 <        theEcr = 0.5 * smallest;
880 >        theEcr = 15.0;
881        }
882        else{
883          theEcr = globals->getECR();
# Line 795 | 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 809 | Line 901 | void SimSetup::finalInfoCheck(void){
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 820 | 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;
829 <          smallest = info[i].boxL[0];
830 <          if (info[i].boxL[1] <= smallest)
831 <            smallest = info[i].boxL[1];
832 <          if (info[i].boxL[2] <= smallest)
833 <            smallest = info[i].boxL[2];
834 <          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 849 | Line 937 | void SimSetup::finalInfoCheck(void){
937          else{
938            theEst = globals->getEST();
939          }
940 <
940 >        
941          info[i].setDefaultEcr(theEcr, theEst);
942        }
943      }
944    }
857 –
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 888 | Line 975 | void SimSetup::initSystemCoords(void){
975      delete fileInit;
976    }
977    else{
978 < #ifdef IS_MPI
892 <
978 >    
979      // no init from bass
980 <
980 >    
981      sprintf(painCave.errMsg,
982 <            "Cannot intialize a parallel simulation without an initial configuration file.\n");
982 >            "Cannot intialize a simulation without an initial configuration file.\n");
983      painCave.isFatal = 1;;
984      simError();
985 <
900 < #else
901 <
902 <    initFromBass();
903 <
904 <
905 < #endif
985 >    
986    }
987  
988   #ifdef IS_MPI
# Line 1198 | 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 1341 | 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          }
# Line 1364 | 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 1374 | 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 1384 | 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          }
# Line 1407 | 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());
1508 +
1509          else{
1510            sprintf(painCave.errMsg,
1511                    "SimSetup error: If you use an NPT\n"
1512 <                  "    ensemble, you must set tauBarostat.\n");
1512 >                  "\tensemble, you must set tauBarostat.\n");
1513            painCave.isFatal = 1;
1514            simError();
1515          }
# Line 1450 | Line 1532 | void SimSetup::makeIntegrator(void){
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          }    
# Line 1460 | Line 1542 | void SimSetup::makeIntegrator(void){
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          }
# Line 1470 | Line 1552 | void SimSetup::makeIntegrator(void){
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          }
# Line 1523 | 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 1539 | 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 1558 | Line 1641 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1641    }
1642    else{
1643      sprintf(painCave.errMsg,
1644 <            "ZConstraint Warning: User does not set force Subtraction policy, "
1645 <            "PolicyByMass 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("BYMASS");
# Line 1603 | Line 1687 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1687    //check the uniqueness of index  
1688    if(!zconsParaData->isIndexUnique()){
1689      sprintf(painCave.errMsg,
1690 <            "ZConstraint Error: molIndex is not unique\n");
1690 >            "ZConstraint Error: molIndex is not unique!\n");
1691      painCave.isFatal = 1;
1692      simError();
1693    }

Diff Legend

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