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 812 by mmeineke, Wed Oct 22 21:17:32 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 31 | 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;
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 107 | 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()) )
# Line 148 | Line 174 | void SimSetup::makeMolecules(void){
174    bend_set* theBends;
175    torsion_set* theTorsions;
176  
151 –
177    //init the forceField paramters
178  
179    the_ff->readParams();
# Line 156 | 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 192 | 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 611 | 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 644 | 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 665 | 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 694 | Line 787 | void SimSetup::gatherInfo(void){
787      }
788  
789      // check for the temperature set flag
790 <
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
702 <
703 <    double boxVector[3];
704 <    if (globals->haveBox()){
705 <      boxVector[0] = globals->getBox();
706 <      boxVector[1] = globals->getBox();
707 <      boxVector[2] = globals->getBox();
708 <
709 <      info[i].setBox(boxVector);
710 <    }
711 <    else if (globals->haveDensity()){
712 <      double vol;
713 <      vol = (double) tot_nmol / globals->getDensity();
714 <      boxVector[0] = pow(vol, (1.0 / 3.0));
715 <      boxVector[1] = boxVector[0];
716 <      boxVector[2] = boxVector[0];
717 <
718 <      info[i].setBox(boxVector);
719 <    }
720 <    else{
721 <      if (!globals->haveBoxX()){
722 <        sprintf(painCave.errMsg,
723 <                "SimSetup error, no periodic BoxX size given.\n");
724 <        painCave.isFatal = 1;
725 <        simError();
726 <      }
727 <      boxVector[0] = globals->getBoxX();
794 >    // check for the extended State init
795  
796 <      if (!globals->haveBoxY()){
797 <        sprintf(painCave.errMsg,
798 <                "SimSetup error, no periodic BoxY size given.\n");
732 <        painCave.isFatal = 1;
733 <        simError();
734 <      }
735 <      boxVector[1] = globals->getBoxY();
736 <
737 <      if (!globals->haveBoxZ()){
738 <        sprintf(painCave.errMsg,
739 <                "SimSetup error, no periodic BoxZ size given.\n");
740 <        painCave.isFatal = 1;
741 <        simError();
742 <      }
743 <      boxVector[2] = globals->getBoxZ();
744 <
745 <      info[i].setBox(boxVector);
746 <    }
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 820 | 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;
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;
880 >        theEcr = 15.0;
881        }
882        else{
883          theEcr = globals->getECR();
# Line 839 | 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 849 | 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 864 | 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;
873 <          smallest = info[i].boxL[0];
874 <          if (info[i].boxL[1] <= smallest)
875 <            smallest = info[i].boxL[1];
876 <          if (info[i].boxL[2] <= smallest)
877 <            smallest = info[i].boxL[2];
878 <          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 893 | 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    }
901 –
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 932 | Line 975 | void SimSetup::initSystemCoords(void){
975      delete fileInit;
976    }
977    else{
978 < #ifdef IS_MPI
936 <
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 <
944 < #else
945 <
946 <    initFromBass();
947 <
948 <
949 < #endif
985 >    
986    }
987  
988   #ifdef IS_MPI
# Line 1242 | 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 1385 | 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 1408 | 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 1418 | 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 1428 | 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 1451 | 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 1494 | 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 1504 | 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 1514 | 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 1567 | 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 1583 | 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 1602 | 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 1647 | 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