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 962 by tim, Mon Jan 19 18:36:21 2004 UTC vs.
Revision 1093 by tim, Wed Mar 17 14:22:59 2004 UTC

# Line 9 | Line 9
9   #include "parse_me.h"
10   #include "Integrator.hpp"
11   #include "simError.h"
12 + //#include "ConjugateMinimizer.hpp"
13 + #include "OOPSEMinimizer.hpp"
14  
15   #ifdef IS_MPI
16   #include "mpiBASS.h"
# Line 24 | Line 26
26   #define NPTxyz_ENS     4
27  
28  
29 < #define FF_DUFF 0
30 < #define FF_LJ   1
31 < #define FF_EAM  2
29 > #define FF_DUFF  0
30 > #define FF_LJ    1
31 > #define FF_EAM   2
32 > #define FF_H2O   3
33  
34   using namespace std;
35  
# Line 144 | Line 147 | void SimSetup::createSim(void){
147  
148    makeOutNames();
149  
150 <  // make the integrator
151 <
152 <  makeIntegrator();
153 <
150 >  if (globals->haveMinimizer())
151 >    // make minimizer
152 >    makeMinimizer();
153 >  else
154 >    // make the integrator
155 >    makeIntegrator();
156 >  
157   #ifdef IS_MPI
158    mpiSim->mpiRefresh();
159   #endif
# Line 174 | Line 180 | void SimSetup::makeMolecules(void){
180    bend_set* theBends;
181    torsion_set* theTorsions;
182  
177 –
183    //init the forceField paramters
184  
185    the_ff->readParams();
# Line 182 | Line 187 | void SimSetup::makeMolecules(void){
187  
188    // init the atoms
189  
190 +  double phi, theta, psi;
191 +  double sux, suy, suz;
192 +  double Axx, Axy, Axz, Ayx, Ayy, Ayz, Azx, Azy, Azz;
193    double ux, uy, uz, u, uSqr;
194  
195    for (k = 0; k < nInfo; k++){
# Line 218 | Line 226 | void SimSetup::makeMolecules(void){
226            info[k].n_oriented++;
227            molInfo.myAtoms[j] = dAtom;
228  
229 <          ux = currentAtom->getOrntX();
230 <          uy = currentAtom->getOrntY();
231 <          uz = currentAtom->getOrntZ();
229 >          // Directional Atoms have standard unit vectors which are oriented
230 >          // in space using the three Euler angles.  We assume the standard
231 >          // unit vector was originally along the z axis below.
232 >
233 >          phi = currentAtom->getEulerPhi() * M_PI / 180.0;
234 >          theta = currentAtom->getEulerTheta() * M_PI / 180.0;
235 >          psi = currentAtom->getEulerPsi()* M_PI / 180.0;
236 >            
237 >          Axx = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
238 >          Axy = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
239 >          Axz = sin(theta) * sin(psi);
240 >          
241 >          Ayx = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
242 >          Ayy = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
243 >          Ayz = sin(theta) * cos(psi);
244 >          
245 >          Azx = sin(phi) * sin(theta);
246 >          Azy = -cos(phi) * sin(theta);
247 >          Azz = cos(theta);
248 >
249 >          sux = 0.0;
250 >          suy = 0.0;
251 >          suz = 1.0;
252 >
253 >          ux = (Axx * sux) + (Ayx * suy) + (Azx * suz);
254 >          uy = (Axy * sux) + (Ayy * suy) + (Azy * suz);
255 >          uz = (Axz * sux) + (Ayz * suy) + (Azz * suz);
256  
257            uSqr = (ux * ux) + (uy * uy) + (uz * uz);
258  
# Line 609 | Line 641 | void SimSetup::gatherInfo(void){
641    else if (!strcasecmp(force_field, "EAM")){
642      ffCase = FF_EAM;
643    }
644 +  else if (!strcasecmp(force_field, "WATER")){
645 +    ffCase = FF_H2O;
646 +  }
647    else{
648      sprintf(painCave.errMsg, "SimSetup Error. Unrecognized force field -> %s\n",
649              force_field);
# Line 637 | Line 672 | void SimSetup::gatherInfo(void){
672    }
673    else{
674      sprintf(painCave.errMsg,
675 <            "SimSetup Warning. Unrecognized Ensemble -> %s, "
676 <            "reverting to NVE for this simulation.\n",
675 >            "SimSetup Warning. Unrecognized Ensemble -> %s \n"
676 >            "\treverting to NVE for this simulation.\n",
677              ensemble);
678           painCave.isFatal = 0;
679           simError();
# Line 670 | Line 705 | void SimSetup::gatherInfo(void){
705        if (!the_components[i]->haveNMol()){
706          // we have a problem
707          sprintf(painCave.errMsg,
708 <                "SimSetup Error. No global NMol or component NMol"
709 <                " given. Cannot calculate the number of atoms.\n");
708 >                "SimSetup Error. No global NMol or component NMol given.\n"
709 >                "\tCannot calculate the number of atoms.\n");
710          painCave.isFatal = 1;
711          simError();
712        }
# Line 694 | Line 729 | void SimSetup::gatherInfo(void){
729    //check whether sample time, status time, thermal time and reset time are divisble by dt
730    if (!isDivisible(globals->getSampleTime(), globals->getDt())){
731      sprintf(painCave.errMsg,
732 <              "Sample time is not divisible by dt \n");
732 >            "Sample time is not divisible by dt.\n"
733 >            "\tThis will result in samples that are not uniformly\n"
734 >            "\tdistributed in time.  If this is a problem, change\n"
735 >            "\tyour sampleTime variable.\n");
736      painCave.isFatal = 0;
737      simError();    
738    }
739  
740    if (globals->haveStatusTime() && !isDivisible(globals->getSampleTime(), globals->getDt())){
741      sprintf(painCave.errMsg,
742 <              "Status time is not divisible by dt\n");
742 >            "Status time is not divisible by dt.\n"
743 >            "\tThis will result in status reports that are not uniformly\n"
744 >            "\tdistributed in time.  If this is a problem, change \n"
745 >            "\tyour statusTime variable.\n");
746      painCave.isFatal = 0;
747      simError();    
748    }
749  
750    if (globals->haveThermalTime() && !isDivisible(globals->getThermalTime(), globals->getDt())){
751      sprintf(painCave.errMsg,
752 <              "Thermal time is not divisible by dt\n");
752 >            "Thermal time is not divisible by dt.\n"
753 >            "\tThis will result in thermalizations that are not uniformly\n"
754 >            "\tdistributed in time.  If this is a problem, change \n"
755 >            "\tyour thermalTime variable.\n");
756      painCave.isFatal = 0;
757      simError();    
758    }  
759  
760    if (globals->haveResetTime() && !isDivisible(globals->getResetTime(), globals->getDt())){
761      sprintf(painCave.errMsg,
762 <              "Reset time is not divisible by dt\n");
762 >            "Reset time is not divisible by dt.\n"
763 >            "\tThis will result in integrator resets that are not uniformly\n"
764 >            "\tdistributed in time.  If this is a problem, change\n"
765 >            "\tyour resetTime variable.\n");
766      painCave.isFatal = 0;
767      simError();    
768    }
# Line 799 | Line 846 | void SimSetup::gatherInfo(void){
846    for (int i = 0; i < nInfo; i++){
847      info[i].setSeed(seedValue);
848    }
849 <
849 >  
850   #ifdef IS_MPI
851 <  strcpy(checkPointMsg, "Succesfully gathered all information from Bass\n");
851 >  strcpy(checkPointMsg, "Successfully gathered all information from Bass\n");
852    MPIcheckPoint();
853   #endif // is_mpi
854   }
# Line 834 | Line 881 | void SimSetup::finalInfoCheck(void){
881  
882        if (!globals->haveECR()){
883          sprintf(painCave.errMsg,
884 <                "SimSetup Warning: using default value of 15.0 angstroms"
885 <                "box length for the electrostaticCutoffRadius.\n");
884 >                "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n"
885 >                "\tOOPSE will use a default value of 15.0 angstroms"
886 >                "\tfor the electrostaticCutoffRadius.\n");
887          painCave.isFatal = 0;
888          simError();
889          theEcr = 15.0;
# Line 846 | Line 894 | void SimSetup::finalInfoCheck(void){
894  
895        if (!globals->haveEST()){
896          sprintf(painCave.errMsg,
897 <                "SimSetup Warning: using default value of 0.05 * the "
898 <                "electrostaticCutoffRadius for the electrostaticSkinThickness\n");
897 >                "SimSetup Warning: No value was set for electrostaticSkinThickness.\n"
898 >                "\tOOPSE will use a default value of\n"
899 >                "\t0.05 * electrostaticCutoffRadius\n"
900 >                "\tfor the electrostaticSkinThickness\n");
901          painCave.isFatal = 0;
902          simError();
903          theEst = 0.05 * theEcr;
# Line 860 | Line 910 | void SimSetup::finalInfoCheck(void){
910  
911        if (!globals->haveDielectric()){
912          sprintf(painCave.errMsg,
913 <                "SimSetup Error: You are trying to use Reaction Field without"
914 <                "setting a dielectric constant!\n");
913 >                "SimSetup Error: No Dielectric constant was set.\n"
914 >                "\tYou are trying to use Reaction Field without"
915 >                "\tsetting a dielectric constant!\n");
916          painCave.isFatal = 1;
917          simError();
918        }
# Line 871 | Line 922 | void SimSetup::finalInfoCheck(void){
922        if (usesDipoles){
923          if (!globals->haveECR()){
924            sprintf(painCave.errMsg,
925 <                  "SimSetup Warning: using default value of 15.0 angstroms"
926 <                  "box length for the electrostaticCutoffRadius.\n");
925 >                  "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n"
926 >                  "\tOOPSE will use a default value of 15.0 angstroms"
927 >                  "\tfor the electrostaticCutoffRadius.\n");
928            painCave.isFatal = 0;
929            simError();
930            theEcr = 15.0;
# Line 883 | Line 935 | void SimSetup::finalInfoCheck(void){
935          
936          if (!globals->haveEST()){
937            sprintf(painCave.errMsg,
938 <                  "SimSetup Warning: using default value of 0.05 * the "
939 <                  "electrostaticCutoffRadius for the "
940 <                  "electrostaticSkinThickness\n");
938 >                  "SimSetup Warning: No value was set for electrostaticSkinThickness.\n"
939 >                  "\tOOPSE will use a default value of\n"
940 >                  "\t0.05 * electrostaticCutoffRadius\n"
941 >                  "\tfor the electrostaticSkinThickness\n");
942            painCave.isFatal = 0;
943            simError();
944            theEst = 0.05 * theEcr;
# Line 1092 | Line 1145 | void SimSetup::createFF(void){
1145        the_ff = new EAM_FF();
1146        break;
1147  
1148 +    case FF_H2O:
1149 +      the_ff = new WATER();
1150 +      break;
1151 +
1152      default:
1153        sprintf(painCave.errMsg,
1154                "SimSetup Error. Unrecognized force field in case statement.\n");
# Line 1234 | Line 1291 | void SimSetup::mpiMolDivide(void){
1291  
1292    if (local_atoms != info[0].n_atoms){
1293      sprintf(painCave.errMsg,
1294 <            "SimSetup error: mpiSim's localAtom (%d) and SimSetup's"
1295 <            " localAtom (%d) are not equal.\n",
1294 >            "SimSetup error: mpiSim's localAtom (%d) and SimSetup's\n"
1295 >            "\tlocalAtom (%d) are not equal.\n",
1296              info[0].n_atoms, local_atoms);
1297      painCave.isFatal = 1;
1298      simError();
# Line 1377 | Line 1434 | void SimSetup::makeIntegrator(void){
1434          else{
1435            sprintf(painCave.errMsg,
1436                    "SimSetup error: If you use the NVT\n"
1437 <                  "    ensemble, you must set tauThermostat.\n");
1437 >                  "\tensemble, you must set tauThermostat.\n");
1438            painCave.isFatal = 1;
1439            simError();
1440          }
# Line 1400 | Line 1457 | void SimSetup::makeIntegrator(void){
1457          else{
1458            sprintf(painCave.errMsg,
1459                    "SimSetup error: If you use a constant pressure\n"
1460 <                  "    ensemble, you must set targetPressure in the BASS file.\n");
1460 >                  "\tensemble, you must set targetPressure in the BASS file.\n");
1461            painCave.isFatal = 1;
1462            simError();
1463          }
# Line 1410 | Line 1467 | void SimSetup::makeIntegrator(void){
1467          else{
1468            sprintf(painCave.errMsg,
1469                    "SimSetup error: If you use an NPT\n"
1470 <                  "    ensemble, you must set tauThermostat.\n");
1470 >                  "\tensemble, you must set tauThermostat.\n");
1471            painCave.isFatal = 1;
1472            simError();
1473          }
# Line 1420 | Line 1477 | void SimSetup::makeIntegrator(void){
1477          else{
1478            sprintf(painCave.errMsg,
1479                    "SimSetup error: If you use an NPT\n"
1480 <                  "    ensemble, you must set tauBarostat.\n");
1480 >                  "\tensemble, you must set tauBarostat.\n");
1481            painCave.isFatal = 1;
1482            simError();
1483          }
# Line 1443 | Line 1500 | void SimSetup::makeIntegrator(void){
1500          else{
1501            sprintf(painCave.errMsg,
1502                    "SimSetup error: If you use a constant pressure\n"
1503 <                  "    ensemble, you must set targetPressure in the BASS file.\n");
1503 >                  "\tensemble, you must set targetPressure in the BASS file.\n");
1504            painCave.isFatal = 1;
1505            simError();
1506          }    
# Line 1454 | Line 1511 | void SimSetup::makeIntegrator(void){
1511          else{
1512            sprintf(painCave.errMsg,
1513                    "SimSetup error: If you use an NPT\n"
1514 <                  "    ensemble, you must set tauThermostat.\n");
1514 >                  "\tensemble, you must set tauThermostat.\n");
1515            painCave.isFatal = 1;
1516            simError();
1517          }
# Line 1465 | Line 1522 | void SimSetup::makeIntegrator(void){
1522          else{
1523            sprintf(painCave.errMsg,
1524                    "SimSetup error: If you use an NPT\n"
1525 <                  "    ensemble, you must set tauBarostat.\n");
1525 >                  "\tensemble, you must set tauBarostat.\n");
1526            painCave.isFatal = 1;
1527            simError();
1528          }
# Line 1488 | Line 1545 | void SimSetup::makeIntegrator(void){
1545          else{
1546            sprintf(painCave.errMsg,
1547                    "SimSetup error: If you use a constant pressure\n"
1548 <                  "    ensemble, you must set targetPressure in the BASS file.\n");
1548 >                  "\tensemble, you must set targetPressure in the BASS file.\n");
1549            painCave.isFatal = 1;
1550            simError();
1551          }    
# Line 1498 | Line 1555 | void SimSetup::makeIntegrator(void){
1555          else{
1556            sprintf(painCave.errMsg,
1557                    "SimSetup error: If you use an NPT\n"
1558 <                  "    ensemble, you must set tauThermostat.\n");
1558 >                  "\tensemble, you must set tauThermostat.\n");
1559            painCave.isFatal = 1;
1560            simError();
1561          }
# Line 1508 | Line 1565 | void SimSetup::makeIntegrator(void){
1565          else{
1566            sprintf(painCave.errMsg,
1567                    "SimSetup error: If you use an NPT\n"
1568 <                  "    ensemble, you must set tauBarostat.\n");
1568 >                  "\tensemble, you must set tauBarostat.\n");
1569            painCave.isFatal = 1;
1570            simError();
1571          }
# Line 1561 | Line 1618 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1618    }
1619    else{
1620      sprintf(painCave.errMsg,
1621 <            "ZConstraint error: If you use an ZConstraint\n"
1622 <            " , you must set sample time.\n");
1621 >            "ZConstraint error: If you use a ZConstraint,\n"
1622 >            "\tyou must set zconsTime.\n");
1623      painCave.isFatal = 1;
1624      simError();
1625    }
# Line 1577 | Line 1634 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1634    else{
1635      double defaultZConsTol = 0.01;
1636      sprintf(painCave.errMsg,
1637 <            "ZConstraint Waring: Tolerance for z-constraint methodl is not specified\n"
1638 <            " , default value %f is used.\n",
1637 >            "ZConstraint Warning: Tolerance for z-constraint method is not specified.\n"
1638 >            "\tOOPSE will use a default value of %f.\n"
1639 >            "\tTo set the tolerance, use the zconsTol variable.\n",
1640              defaultZConsTol);
1641      painCave.isFatal = 0;
1642      simError();      
# Line 1596 | Line 1654 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1654    }
1655    else{
1656      sprintf(painCave.errMsg,
1657 <            "ZConstraint Warning: User does not set force Subtraction policy, "
1658 <            "PolicyByMass is used\n");
1657 >            "ZConstraint Warning: No force subtraction policy was set.\n"
1658 >            "\tOOPSE will use PolicyByMass.\n"
1659 >            "\tTo set the policy, use the zconsForcePolicy variable.\n");
1660      painCave.isFatal = 0;
1661      simError();
1662      zconsForcePolicy->setData("BYMASS");
# Line 1605 | Line 1664 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1664  
1665    theInfo.addProperty(zconsForcePolicy);
1666  
1667 +  //set zcons gap
1668 +  DoubleData* zconsGap = new DoubleData();
1669 +  zconsGap->setID(ZCONSGAP_ID);
1670 +
1671 +  if (globals->haveZConsGap()){
1672 +    zconsGap->setData(globals->getZconsGap());
1673 +    theInfo.addProperty(zconsGap);  
1674 +  }
1675 +
1676 +  //set zcons fixtime
1677 +  DoubleData* zconsFixtime = new DoubleData();
1678 +  zconsFixtime->setID(ZCONSFIXTIME_ID);
1679 +
1680 +  if (globals->haveZConsFixTime()){
1681 +    zconsFixtime->setData(globals->getZconsFixtime());
1682 +    theInfo.addProperty(zconsFixtime);  
1683 +  }
1684 +
1685 +  //set zconsUsingSMD
1686 +  IntData* zconsUsingSMD = new IntData();
1687 +  zconsUsingSMD->setID(ZCONSUSINGSMD_ID);
1688 +
1689 +  if (globals->haveZConsUsingSMD()){
1690 +    zconsUsingSMD->setData(globals->getZconsUsingSMD());
1691 +    theInfo.addProperty(zconsUsingSMD);  
1692 +  }
1693 +
1694    //Determine the name of ouput file and add it into SimInfo's property list
1695    //Be careful, do not use inFileName, since it is a pointer which
1696    //point to a string at master node, and slave nodes do not contain that string
# Line 1634 | Line 1720 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1720      tempParaItem.zPos = zconStamp[i]->getZpos();
1721      tempParaItem.zconsIndex = zconStamp[i]->getMolIndex();
1722      tempParaItem.kRatio = zconStamp[i]->getKratio();
1723 <
1723 >    tempParaItem.havingCantVel = zconStamp[i]->haveCantVel();
1724 >    tempParaItem.cantVel = zconStamp[i]->getCantVel();    
1725      zconsParaData->addItem(tempParaItem);
1726    }
1727  
1728    //check the uniqueness of index  
1729    if(!zconsParaData->isIndexUnique()){
1730      sprintf(painCave.errMsg,
1731 <            "ZConstraint Error: molIndex is not unique\n");
1731 >            "ZConstraint Error: molIndex is not unique!\n");
1732      painCave.isFatal = 1;
1733      simError();
1734    }
# Line 1651 | Line 1738 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1738    
1739    //push data into siminfo, therefore, we can retrieve later
1740    theInfo.addProperty(zconsParaData);
1741 + }
1742 +
1743 + void SimSetup::makeMinimizer(){
1744 +
1745 +  OOPSEMinimizer* myOOPSEMinimizer;
1746 +  MinimizerParameterSet* param;
1747 +  char minimizerName[100];
1748 +  
1749 +  for (int i = 0; i < nInfo; i++){
1750 +    
1751 +    //prepare parameter set for minimizer
1752 +    param = new MinimizerParameterSet();
1753 +    param->setDefaultParameter();
1754 +
1755 +    if (globals->haveMinimizer()){
1756 +      param->setFTol(globals->getMinFTol());
1757 +    }
1758 +
1759 +    if (globals->haveMinGTol()){
1760 +      param->setGTol(globals->getMinGTol());
1761 +    }
1762 +
1763 +    if (globals->haveMinMaxIter()){
1764 +      param->setMaxIteration(globals->getMinMaxIter());
1765 +    }
1766 +
1767 +    if (globals->haveMinWriteFrq()){
1768 +      param->setMaxIteration(globals->getMinMaxIter());
1769 +    }
1770 +
1771 +    if (globals->haveMinWriteFrq()){
1772 +      param->setWriteFrq(globals->getMinWriteFrq());
1773 +    }
1774 +    
1775 +    if (globals->haveMinStepSize()){
1776 +      param->setStepSize(globals->getMinStepSize());
1777 +    }
1778 +
1779 +    if (globals->haveMinLSMaxIter()){
1780 +      param->setLineSearchMaxIteration(globals->getMinLSMaxIter());
1781 +    }    
1782 +
1783 +    if (globals->haveMinLSTol()){
1784 +      param->setLineSearchTol(globals->getMinLSTol());
1785 +    }    
1786 +
1787 +    strcpy(minimizerName, globals->getMinimizer());
1788 +
1789 +    if (!strcasecmp(minimizerName, "CG")){
1790 +      myOOPSEMinimizer = new PRCGMinimizer(&(info[i]), the_ff, param);
1791 +    }
1792 +    else if (!strcasecmp(minimizerName, "SD")){
1793 +    //myOOPSEMinimizer = MinimizerFactory.creatMinimizer("", &(info[i]), the_ff, param);
1794 +      myOOPSEMinimizer = new SDMinimizer(&(info[i]), the_ff, param);
1795 +    }
1796 +    else{
1797 +          sprintf(painCave.errMsg,
1798 +                  "SimSetup error: Unrecognized Minimizer, use Conjugate Gradient \n");
1799 +          painCave.isFatal = 0;
1800 +          simError();
1801 +
1802 +      myOOPSEMinimizer = new PRCGMinimizer(&(info[i]), the_ff, param);          
1803 +    }
1804 +     info[i].the_integrator = myOOPSEMinimizer;
1805 +
1806 +     //store the minimizer into simInfo
1807 +     info[i].the_minimizer = myOOPSEMinimizer;
1808 +     info[i].has_minimizer = true;
1809 +  }
1810 +
1811   }

Diff Legend

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