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 1035 by tim, Fri Feb 6 21:37: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  
14   #ifdef IS_MPI
15   #include "mpiBASS.h"
# Line 24 | Line 25
25   #define NPTxyz_ENS     4
26  
27  
28 < #define FF_DUFF 0
29 < #define FF_LJ   1
30 < #define FF_EAM  2
28 > #define FF_DUFF  0
29 > #define FF_LJ    1
30 > #define FF_EAM   2
31 > #define FF_H2O 3
32  
33   using namespace std;
34  
# Line 144 | Line 146 | void SimSetup::createSim(void){
146  
147    makeOutNames();
148  
149 <  // make the integrator
150 <
151 <  makeIntegrator();
152 <
149 >  if (globals->haveMinimizer())
150 >    // make minimizer
151 >    makeMinimizer();
152 >  else
153 >    // make the integrator
154 >    makeIntegrator();
155 >  
156   #ifdef IS_MPI
157    mpiSim->mpiRefresh();
158   #endif
# Line 174 | Line 179 | void SimSetup::makeMolecules(void){
179    bend_set* theBends;
180    torsion_set* theTorsions;
181  
177 –
182    //init the forceField paramters
183  
184    the_ff->readParams();
# Line 182 | Line 186 | void SimSetup::makeMolecules(void){
186  
187    // init the atoms
188  
189 +  double phi, theta, psi;
190 +  double sux, suy, suz;
191 +  double Axx, Axy, Axz, Ayx, Ayy, Ayz, Azx, Azy, Azz;
192    double ux, uy, uz, u, uSqr;
193  
194    for (k = 0; k < nInfo; k++){
# Line 218 | Line 225 | void SimSetup::makeMolecules(void){
225            info[k].n_oriented++;
226            molInfo.myAtoms[j] = dAtom;
227  
228 <          ux = currentAtom->getOrntX();
229 <          uy = currentAtom->getOrntY();
230 <          uz = currentAtom->getOrntZ();
228 >          // Directional Atoms have standard unit vectors which are oriented
229 >          // in space using the three Euler angles.  We assume the standard
230 >          // unit vector was originally along the z axis below.
231 >
232 >          phi = currentAtom->getEulerPhi() * M_PI / 180.0;
233 >          theta = currentAtom->getEulerTheta() * M_PI / 180.0;
234 >          psi = currentAtom->getEulerPsi()* M_PI / 180.0;
235 >            
236 >          Axx = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
237 >          Axy = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
238 >          Axz = sin(theta) * sin(psi);
239 >          
240 >          Ayx = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
241 >          Ayy = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
242 >          Ayz = sin(theta) * cos(psi);
243 >          
244 >          Azx = sin(phi) * sin(theta);
245 >          Azy = -cos(phi) * sin(theta);
246 >          Azz = cos(theta);
247  
248 +          sux = 0.0;
249 +          suy = 0.0;
250 +          suz = 1.0;
251 +
252 +          ux = (Axx * sux) + (Ayx * suy) + (Azx * suz);
253 +          uy = (Axy * sux) + (Ayy * suy) + (Azy * suz);
254 +          uz = (Axz * sux) + (Ayz * suy) + (Azz * suz);
255 +
256            uSqr = (ux * ux) + (uy * uy) + (uz * uz);
257  
258            u = sqrt(uSqr);
# Line 609 | Line 640 | void SimSetup::gatherInfo(void){
640    else if (!strcasecmp(force_field, "EAM")){
641      ffCase = FF_EAM;
642    }
643 +  else if (!strcasecmp(force_field, "WATER")){
644 +    ffCase = FF_H2O;
645 +  }
646    else{
647      sprintf(painCave.errMsg, "SimSetup Error. Unrecognized force field -> %s\n",
648              force_field);
# Line 637 | Line 671 | void SimSetup::gatherInfo(void){
671    }
672    else{
673      sprintf(painCave.errMsg,
674 <            "SimSetup Warning. Unrecognized Ensemble -> %s, "
675 <            "reverting to NVE for this simulation.\n",
674 >            "SimSetup Warning. Unrecognized Ensemble -> %s \n"
675 >            "\treverting to NVE for this simulation.\n",
676              ensemble);
677           painCave.isFatal = 0;
678           simError();
# Line 670 | Line 704 | void SimSetup::gatherInfo(void){
704        if (!the_components[i]->haveNMol()){
705          // we have a problem
706          sprintf(painCave.errMsg,
707 <                "SimSetup Error. No global NMol or component NMol"
708 <                " given. Cannot calculate the number of atoms.\n");
707 >                "SimSetup Error. No global NMol or component NMol given.\n"
708 >                "\tCannot calculate the number of atoms.\n");
709          painCave.isFatal = 1;
710          simError();
711        }
# Line 694 | Line 728 | void SimSetup::gatherInfo(void){
728    //check whether sample time, status time, thermal time and reset time are divisble by dt
729    if (!isDivisible(globals->getSampleTime(), globals->getDt())){
730      sprintf(painCave.errMsg,
731 <              "Sample time is not divisible by dt \n");
731 >            "Sample time is not divisible by dt.\n"
732 >            "\tThis will result in samples that are not uniformly\n"
733 >            "\tdistributed in time.  If this is a problem, change\n"
734 >            "\tyour sampleTime variable.\n");
735      painCave.isFatal = 0;
736      simError();    
737    }
738  
739    if (globals->haveStatusTime() && !isDivisible(globals->getSampleTime(), globals->getDt())){
740      sprintf(painCave.errMsg,
741 <              "Status time is not divisible by dt\n");
741 >            "Status time is not divisible by dt.\n"
742 >            "\tThis will result in status reports that are not uniformly\n"
743 >            "\tdistributed in time.  If this is a problem, change \n"
744 >            "\tyour statusTime variable.\n");
745      painCave.isFatal = 0;
746      simError();    
747    }
748  
749    if (globals->haveThermalTime() && !isDivisible(globals->getThermalTime(), globals->getDt())){
750      sprintf(painCave.errMsg,
751 <              "Thermal time is not divisible by dt\n");
751 >            "Thermal time is not divisible by dt.\n"
752 >            "\tThis will result in thermalizations that are not uniformly\n"
753 >            "\tdistributed in time.  If this is a problem, change \n"
754 >            "\tyour thermalTime variable.\n");
755      painCave.isFatal = 0;
756      simError();    
757    }  
758  
759    if (globals->haveResetTime() && !isDivisible(globals->getResetTime(), globals->getDt())){
760      sprintf(painCave.errMsg,
761 <              "Reset time is not divisible by dt\n");
761 >            "Reset time is not divisible by dt.\n"
762 >            "\tThis will result in integrator resets that are not uniformly\n"
763 >            "\tdistributed in time.  If this is a problem, change\n"
764 >            "\tyour resetTime variable.\n");
765      painCave.isFatal = 0;
766      simError();    
767    }
# Line 799 | Line 845 | void SimSetup::gatherInfo(void){
845    for (int i = 0; i < nInfo; i++){
846      info[i].setSeed(seedValue);
847    }
848 <
848 >  
849   #ifdef IS_MPI
850 <  strcpy(checkPointMsg, "Succesfully gathered all information from Bass\n");
850 >  strcpy(checkPointMsg, "Successfully gathered all information from Bass\n");
851    MPIcheckPoint();
852   #endif // is_mpi
853   }
# Line 834 | Line 880 | void SimSetup::finalInfoCheck(void){
880  
881        if (!globals->haveECR()){
882          sprintf(painCave.errMsg,
883 <                "SimSetup Warning: using default value of 15.0 angstroms"
884 <                "box length for the electrostaticCutoffRadius.\n");
883 >                "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n"
884 >                "\tOOPSE will use a default value of 15.0 angstroms"
885 >                "\tfor the electrostaticCutoffRadius.\n");
886          painCave.isFatal = 0;
887          simError();
888          theEcr = 15.0;
# Line 846 | Line 893 | void SimSetup::finalInfoCheck(void){
893  
894        if (!globals->haveEST()){
895          sprintf(painCave.errMsg,
896 <                "SimSetup Warning: using default value of 0.05 * the "
897 <                "electrostaticCutoffRadius for the electrostaticSkinThickness\n");
896 >                "SimSetup Warning: No value was set for electrostaticSkinThickness.\n"
897 >                "\tOOPSE will use a default value of\n"
898 >                "\t0.05 * electrostaticCutoffRadius\n"
899 >                "\tfor the electrostaticSkinThickness\n");
900          painCave.isFatal = 0;
901          simError();
902          theEst = 0.05 * theEcr;
# Line 860 | Line 909 | void SimSetup::finalInfoCheck(void){
909  
910        if (!globals->haveDielectric()){
911          sprintf(painCave.errMsg,
912 <                "SimSetup Error: You are trying to use Reaction Field without"
913 <                "setting a dielectric constant!\n");
912 >                "SimSetup Error: No Dielectric constant was set.\n"
913 >                "\tYou are trying to use Reaction Field without"
914 >                "\tsetting a dielectric constant!\n");
915          painCave.isFatal = 1;
916          simError();
917        }
# Line 871 | Line 921 | void SimSetup::finalInfoCheck(void){
921        if (usesDipoles){
922          if (!globals->haveECR()){
923            sprintf(painCave.errMsg,
924 <                  "SimSetup Warning: using default value of 15.0 angstroms"
925 <                  "box length for the electrostaticCutoffRadius.\n");
924 >                  "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n"
925 >                  "\tOOPSE will use a default value of 15.0 angstroms"
926 >                  "\tfor the electrostaticCutoffRadius.\n");
927            painCave.isFatal = 0;
928            simError();
929            theEcr = 15.0;
# Line 883 | Line 934 | void SimSetup::finalInfoCheck(void){
934          
935          if (!globals->haveEST()){
936            sprintf(painCave.errMsg,
937 <                  "SimSetup Warning: using default value of 0.05 * the "
938 <                  "electrostaticCutoffRadius for the "
939 <                  "electrostaticSkinThickness\n");
937 >                  "SimSetup Warning: No value was set for electrostaticSkinThickness.\n"
938 >                  "\tOOPSE will use a default value of\n"
939 >                  "\t0.05 * electrostaticCutoffRadius\n"
940 >                  "\tfor the electrostaticSkinThickness\n");
941            painCave.isFatal = 0;
942            simError();
943            theEst = 0.05 * theEcr;
# Line 1092 | Line 1144 | void SimSetup::createFF(void){
1144        the_ff = new EAM_FF();
1145        break;
1146  
1147 +    case FF_H2O:
1148 +      the_ff = new WATER();
1149 +      break;
1150 +
1151      default:
1152        sprintf(painCave.errMsg,
1153                "SimSetup Error. Unrecognized force field in case statement.\n");
# Line 1234 | Line 1290 | void SimSetup::mpiMolDivide(void){
1290  
1291    if (local_atoms != info[0].n_atoms){
1292      sprintf(painCave.errMsg,
1293 <            "SimSetup error: mpiSim's localAtom (%d) and SimSetup's"
1294 <            " localAtom (%d) are not equal.\n",
1293 >            "SimSetup error: mpiSim's localAtom (%d) and SimSetup's\n"
1294 >            "\tlocalAtom (%d) are not equal.\n",
1295              info[0].n_atoms, local_atoms);
1296      painCave.isFatal = 1;
1297      simError();
# Line 1377 | Line 1433 | void SimSetup::makeIntegrator(void){
1433          else{
1434            sprintf(painCave.errMsg,
1435                    "SimSetup error: If you use the NVT\n"
1436 <                  "    ensemble, you must set tauThermostat.\n");
1436 >                  "\tensemble, you must set tauThermostat.\n");
1437            painCave.isFatal = 1;
1438            simError();
1439          }
# Line 1400 | Line 1456 | void SimSetup::makeIntegrator(void){
1456          else{
1457            sprintf(painCave.errMsg,
1458                    "SimSetup error: If you use a constant pressure\n"
1459 <                  "    ensemble, you must set targetPressure in the BASS file.\n");
1459 >                  "\tensemble, you must set targetPressure in the BASS file.\n");
1460            painCave.isFatal = 1;
1461            simError();
1462          }
# Line 1410 | Line 1466 | void SimSetup::makeIntegrator(void){
1466          else{
1467            sprintf(painCave.errMsg,
1468                    "SimSetup error: If you use an NPT\n"
1469 <                  "    ensemble, you must set tauThermostat.\n");
1469 >                  "\tensemble, you must set tauThermostat.\n");
1470            painCave.isFatal = 1;
1471            simError();
1472          }
# Line 1420 | Line 1476 | void SimSetup::makeIntegrator(void){
1476          else{
1477            sprintf(painCave.errMsg,
1478                    "SimSetup error: If you use an NPT\n"
1479 <                  "    ensemble, you must set tauBarostat.\n");
1479 >                  "\tensemble, you must set tauBarostat.\n");
1480            painCave.isFatal = 1;
1481            simError();
1482          }
# Line 1443 | Line 1499 | void SimSetup::makeIntegrator(void){
1499          else{
1500            sprintf(painCave.errMsg,
1501                    "SimSetup error: If you use a constant pressure\n"
1502 <                  "    ensemble, you must set targetPressure in the BASS file.\n");
1502 >                  "\tensemble, you must set targetPressure in the BASS file.\n");
1503            painCave.isFatal = 1;
1504            simError();
1505          }    
# Line 1454 | Line 1510 | void SimSetup::makeIntegrator(void){
1510          else{
1511            sprintf(painCave.errMsg,
1512                    "SimSetup error: If you use an NPT\n"
1513 <                  "    ensemble, you must set tauThermostat.\n");
1513 >                  "\tensemble, you must set tauThermostat.\n");
1514            painCave.isFatal = 1;
1515            simError();
1516          }
# Line 1465 | Line 1521 | void SimSetup::makeIntegrator(void){
1521          else{
1522            sprintf(painCave.errMsg,
1523                    "SimSetup error: If you use an NPT\n"
1524 <                  "    ensemble, you must set tauBarostat.\n");
1524 >                  "\tensemble, you must set tauBarostat.\n");
1525            painCave.isFatal = 1;
1526            simError();
1527          }
# Line 1488 | Line 1544 | void SimSetup::makeIntegrator(void){
1544          else{
1545            sprintf(painCave.errMsg,
1546                    "SimSetup error: If you use a constant pressure\n"
1547 <                  "    ensemble, you must set targetPressure in the BASS file.\n");
1547 >                  "\tensemble, you must set targetPressure in the BASS file.\n");
1548            painCave.isFatal = 1;
1549            simError();
1550          }    
# Line 1498 | Line 1554 | void SimSetup::makeIntegrator(void){
1554          else{
1555            sprintf(painCave.errMsg,
1556                    "SimSetup error: If you use an NPT\n"
1557 <                  "    ensemble, you must set tauThermostat.\n");
1557 >                  "\tensemble, you must set tauThermostat.\n");
1558            painCave.isFatal = 1;
1559            simError();
1560          }
# Line 1508 | Line 1564 | void SimSetup::makeIntegrator(void){
1564          else{
1565            sprintf(painCave.errMsg,
1566                    "SimSetup error: If you use an NPT\n"
1567 <                  "    ensemble, you must set tauBarostat.\n");
1567 >                  "\tensemble, you must set tauBarostat.\n");
1568            painCave.isFatal = 1;
1569            simError();
1570          }
# Line 1561 | Line 1617 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1617    }
1618    else{
1619      sprintf(painCave.errMsg,
1620 <            "ZConstraint error: If you use an ZConstraint\n"
1621 <            " , you must set sample time.\n");
1620 >            "ZConstraint error: If you use a ZConstraint,\n"
1621 >            "\tyou must set zconsTime.\n");
1622      painCave.isFatal = 1;
1623      simError();
1624    }
# Line 1577 | Line 1633 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1633    else{
1634      double defaultZConsTol = 0.01;
1635      sprintf(painCave.errMsg,
1636 <            "ZConstraint Waring: Tolerance for z-constraint methodl is not specified\n"
1637 <            " , default value %f is used.\n",
1636 >            "ZConstraint Warning: Tolerance for z-constraint method is not specified.\n"
1637 >            "\tOOPSE will use a default value of %f.\n"
1638 >            "\tTo set the tolerance, use the zconsTol variable.\n",
1639              defaultZConsTol);
1640      painCave.isFatal = 0;
1641      simError();      
# Line 1596 | Line 1653 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1653    }
1654    else{
1655      sprintf(painCave.errMsg,
1656 <            "ZConstraint Warning: User does not set force Subtraction policy, "
1657 <            "PolicyByMass is used\n");
1656 >            "ZConstraint Warning: No force subtraction policy was set.\n"
1657 >            "\tOOPSE will use PolicyByMass.\n"
1658 >            "\tTo set the policy, use the zconsForcePolicy variable.\n");
1659      painCave.isFatal = 0;
1660      simError();
1661      zconsForcePolicy->setData("BYMASS");
# Line 1641 | Line 1699 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1699    //check the uniqueness of index  
1700    if(!zconsParaData->isIndexUnique()){
1701      sprintf(painCave.errMsg,
1702 <            "ZConstraint Error: molIndex is not unique\n");
1702 >            "ZConstraint Error: molIndex is not unique!\n");
1703      painCave.isFatal = 1;
1704      simError();
1705    }
# Line 1651 | Line 1709 | void SimSetup::setupZConstraint(SimInfo& theInfo){
1709    
1710    //push data into siminfo, therefore, we can retrieve later
1711    theInfo.addProperty(zconsParaData);
1712 + }
1713 +
1714 + void SimSetup::makeMinimizer(){
1715 +
1716 +  OOPSEMinimizerBase* myOOPSEMinimizerBase;
1717 +  ObjFunctor1 * objFunc;
1718 +  OutputFunctor* outputFunc;
1719 +  ConcreteNLModel1* nlp;
1720 +  MinimizerParameterSet* param;
1721 +  ConjugateMinimizerBase* minimizer;
1722 +  int dim;
1723 +  
1724 +  for (int i = 0; i < nInfo; i++){
1725 +    //creat
1726 +    myOOPSEMinimizerBase = new OOPSEMinimizerBase(&(info[i]), the_ff);
1727 +
1728 +     info[i].the_integrator = myOOPSEMinimizerBase;
1729 +    //creat the object functor;
1730 +    objFunc = (ObjFunctor1*) new ClassMemObjFunctor1<OOPSEMinimizerBase>
1731 +                                              (myOOPSEMinimizerBase, &OOPSEMinimizerBase::calcGradient);
1732 +
1733 +    //creat output functor;
1734 +    outputFunc =  new ClassMemOutputFunctor<OOPSEMinimizerBase>
1735 +                               (myOOPSEMinimizerBase, &OOPSEMinimizerBase::output);
1736 +
1737 +    //creat nonlinear model
1738 +    dim = myOOPSEMinimizerBase->getDim();    
1739 +    nlp = new ConcreteNLModel1(dim, objFunc);
1740 +
1741 +    nlp->setX(myOOPSEMinimizerBase->getCoor());
1742 +
1743 +    //prepare parameter set for minimizer
1744 +    param = new MinimizerParameterSet();
1745 +    param->setDefaultParameter();
1746 +
1747 +    if (globals->haveMinimizer()){
1748 +      param->setFTol(globals->getMinFTol());
1749 +    }
1750 +
1751 +    if (globals->haveMinGTol()){
1752 +      param->setGTol(globals->getMinGTol());
1753 +    }
1754 +
1755 +    if (globals->haveMinMaxIter()){
1756 +      param->setMaxIteration(globals->getMinMaxIter());
1757 +    }
1758 +
1759 +    if (globals->haveMinWriteFrq()){
1760 +      param->setMaxIteration(globals->getMinMaxIter());
1761 +    }
1762 +
1763 +    if (globals->haveMinWriteFrq()){
1764 +      param->setWriteFrq(globals->getMinWriteFrq());
1765 +    }
1766 +    
1767 +    if (globals->haveMinResetFrq()){
1768 +      param->setResetFrq(globals->getMinResetFrq());
1769 +    }
1770 +
1771 +    if (globals->haveMinLSMaxIter()){
1772 +      param->setLineSearchMaxIteration(globals->getMinLSMaxIter());
1773 +    }    
1774 +
1775 +    if (globals->haveMinLSTol()){
1776 +      param->setLineSearchTol(globals->getMinLSTol());
1777 +    }    
1778 +    
1779 +     //creat the minimizer
1780 +     minimizer = new PRCGMinimizer(nlp, param);
1781 +     minimizer->setLineSearchStrategy(nlp, GoldenSection);
1782 +     minimizer->setOutputFunctor(outputFunc);
1783 +
1784 +     //store the minimizer into simInfo
1785 +     info[i].the_minimizer = minimizer;
1786 +     info[i].has_minimizer = true;
1787 +  }
1788 +
1789   }

Diff Legend

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