# | 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 | |
# | Line 32 | Line 31 | SimSetup::SimSetup(){ | |
31 | using namespace std; | |
32 | ||
33 | SimSetup::SimSetup(){ | |
34 | + | |
35 | + | initSuspend = false; |
36 | isInfoArray = 0; | |
37 | nInfo = 1; | |
38 | ||
# | Line 54 | Line 55 | void SimSetup::setSimInfo(SimInfo* the_info, int theNi | |
55 | info = the_info; | |
56 | nInfo = theNinfo; | |
57 | isInfoArray = 1; | |
58 | + | initSuspend = true; |
59 | } | |
60 | ||
61 | ||
# | Line 92 | Line 94 | void SimSetup::createSim(void){ | |
94 | #endif // is_mpi | |
95 | ||
96 | void SimSetup::createSim(void){ | |
95 | – | int i, j, k, globalAtomIndex; |
97 | ||
98 | // gather all of the information from the Bass file | |
99 | ||
# | Line 101 | Line 102 | void SimSetup::createSim(void){ | |
102 | // creation of complex system objects | |
103 | ||
104 | sysObjectsCreation(); | |
104 | – | |
105 | – | // check on the post processing info |
105 | ||
107 | – | finalInfoCheck(); |
108 | – | |
106 | // initialize the system coordinates | |
107 | ||
108 | < | if (!isInfoArray){ |
108 | > | if ( !initSuspend ){ |
109 | initSystemCoords(); | |
110 | + | |
111 | + | if( !(globals->getUseInitTime()) ) |
112 | + | info[0].currentTime = 0.0; |
113 | } | |
114 | ||
115 | + | // check on the post processing info |
116 | + | |
117 | + | finalInfoCheck(); |
118 | + | |
119 | // make the output filenames | |
120 | ||
121 | makeOutNames(); | |
# | Line 131 | Line 135 | void SimSetup::makeMolecules(void){ | |
135 | ||
136 | ||
137 | void SimSetup::makeMolecules(void){ | |
138 | < | int k, l; |
138 | > | int k; |
139 | int i, j, exI, exJ, tempEx, stampID, atomOffset, excludeOffset; | |
140 | molInit molInfo; | |
141 | DirectionalAtom* dAtom; | |
# | Line 553 | Line 557 | void SimSetup::gatherInfo(void){ | |
557 | ||
558 | ||
559 | void SimSetup::gatherInfo(void){ | |
560 | < | int i, j, k; |
560 | > | int i; |
561 | ||
562 | ensembleCase = -1; | |
563 | ffCase = -1; | |
# | Line 604 | Line 608 | void SimSetup::gatherInfo(void){ | |
608 | else if (!strcasecmp(ensemble, "NPTf")){ | |
609 | ensembleCase = NPTf_ENS; | |
610 | } | |
611 | < | else if (!strcasecmp(ensemble, "NPTim")){ |
612 | < | ensembleCase = NPTim_ENS; |
611 | > | else if (!strcasecmp(ensemble, "NPTxyz")){ |
612 | > | ensembleCase = NPTxyz_ENS; |
613 | } | |
610 | – | else if (!strcasecmp(ensemble, "NPTfm")){ |
611 | – | ensembleCase = NPTfm_ENS; |
612 | – | } |
614 | else{ | |
615 | sprintf(painCave.errMsg, | |
616 | "SimSetup Warning. Unrecognized Ensemble -> %s, " | |
# | Line 688 | Line 689 | void SimSetup::gatherInfo(void){ | |
689 | info[i].thermalTime = globals->getThermalTime(); | |
690 | } | |
691 | ||
692 | < | // check for the temperature set flag |
692 | > | info[i].resetIntegrator = 0; |
693 | > | if( globals->haveResetTime() ){ |
694 | > | info[i].resetTime = globals->getResetTime(); |
695 | > | info[i].resetIntegrator = 1; |
696 | > | } |
697 | ||
698 | + | // check for the temperature set flag |
699 | + | |
700 | if (globals->haveTempSet()) | |
701 | info[i].setTemp = globals->getTempSet(); | |
702 | < | |
703 | < | // get some of the tricky things that may still be in the globals |
704 | < | |
698 | < | double boxVector[3]; |
699 | < | if (globals->haveBox()){ |
700 | < | 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 | < | } |
742 | < | } |
743 | < | |
702 | > | |
703 | > | } |
704 | > | |
705 | //setup seed for random number generator | |
706 | int seedValue; | |
707 | ||
# | Line 844 | Line 805 | void SimSetup::finalInfoCheck(void){ | |
805 | theEst = globals->getEST(); | |
806 | } | |
807 | ||
808 | < | info[i].setEcr(theEcr, theEst); |
808 | > | info[i].setDefaultEcr(theEcr, theEst); |
809 | ||
810 | if (!globals->haveDielectric()){ | |
811 | sprintf(painCave.errMsg, | |
# | Line 889 | Line 850 | void SimSetup::finalInfoCheck(void){ | |
850 | theEst = globals->getEST(); | |
851 | } | |
852 | ||
853 | < | info[i].setEcr(theEcr, theEst); |
853 | > | info[i].setDefaultEcr(theEcr, theEst); |
854 | } | |
855 | } | |
856 | + | info[i].checkCutOffs(); |
857 | } | |
858 | ||
859 | #ifdef IS_MPI | |
# | Line 916 | Line 878 | void SimSetup::initSystemCoords(void){ | |
878 | if (worldRank == 0){ | |
879 | #endif //is_mpi | |
880 | inName = globals->getInitialConfig(); | |
919 | – | double* tempDouble = new double[1000000]; |
881 | fileInit = new InitializeFromFile(inName); | |
882 | #ifdef IS_MPI | |
883 | } | |
# | Line 934 | Line 895 | void SimSetup::initSystemCoords(void){ | |
895 | ||
896 | sprintf(painCave.errMsg, | |
897 | "Cannot intialize a parallel simulation without an initial configuration file.\n"); | |
898 | < | painCave.isFatal; |
898 | > | painCave.isFatal = 1;; |
899 | simError(); | |
900 | ||
901 | #else | |
# | Line 1160 | Line 1121 | void SimSetup::calcSysValues(void){ | |
1121 | } | |
1122 | ||
1123 | void SimSetup::calcSysValues(void){ | |
1124 | < | int i, j, k; |
1124 | > | int i; |
1125 | ||
1126 | int* molMembershipArray; | |
1127 | ||
# | Line 1259 | Line 1220 | void SimSetup::makeSysArrays(void){ | |
1220 | ||
1221 | ||
1222 | void SimSetup::makeSysArrays(void){ | |
1223 | < | int i, j, k, l; |
1223 | > | |
1224 | > | #ifndef IS_MPI |
1225 | > | int k, j; |
1226 | > | #endif // is_mpi |
1227 | > | int i, l; |
1228 | ||
1229 | Atom** the_atoms; | |
1230 | Molecule* the_molecules; | |
# | Line 1342 | Line 1307 | void SimSetup::makeIntegrator(void){ | |
1307 | void SimSetup::makeIntegrator(void){ | |
1308 | int k; | |
1309 | ||
1310 | + | NVE<RealIntegrator>* myNVE = NULL; |
1311 | NVT<RealIntegrator>* myNVT = NULL; | |
1312 | < | NPTi<RealIntegrator>* myNPTi = NULL; |
1313 | < | NPTf<RealIntegrator>* myNPTf = NULL; |
1314 | < | NPTim<RealIntegrator>* myNPTim = NULL; |
1315 | < | NPTfm<RealIntegrator>* myNPTfm = NULL; |
1350 | < | |
1312 | > | NPTi<NPT<RealIntegrator> >* myNPTi = NULL; |
1313 | > | NPTf<NPT<RealIntegrator> >* myNPTf = NULL; |
1314 | > | NPTxyz<NPT<RealIntegrator> >* myNPTxyz = NULL; |
1315 | > | |
1316 | for (k = 0; k < nInfo; k++){ | |
1317 | switch (ensembleCase){ | |
1318 | case NVE_ENS: | |
1319 | if (globals->haveZconstraints()){ | |
1320 | setupZConstraint(info[k]); | |
1321 | < | new ZConstraint<NVE<RealIntegrator> >(&(info[k]), the_ff); |
1321 | > | myNVE = new ZConstraint<NVE<RealIntegrator> >(&(info[k]), the_ff); |
1322 | } | |
1323 | < | else |
1324 | < | new NVE<RealIntegrator>(&(info[k]), the_ff); |
1323 | > | else{ |
1324 | > | myNVE = new NVE<RealIntegrator>(&(info[k]), the_ff); |
1325 | > | } |
1326 | > | |
1327 | > | info->the_integrator = myNVE; |
1328 | break; | |
1329 | ||
1330 | case NVT_ENS: | |
# | Line 1378 | Line 1346 | void SimSetup::makeIntegrator(void){ | |
1346 | painCave.isFatal = 1; | |
1347 | simError(); | |
1348 | } | |
1349 | + | |
1350 | + | info->the_integrator = myNVT; |
1351 | break; | |
1352 | ||
1353 | case NPTi_ENS: | |
1354 | if (globals->haveZconstraints()){ | |
1355 | setupZConstraint(info[k]); | |
1356 | < | myNPTi = new ZConstraint<NPTi<RealIntegrator> >(&(info[k]), the_ff); |
1356 | > | myNPTi = new ZConstraint<NPTi<NPT <RealIntegrator> > >(&(info[k]), the_ff); |
1357 | } | |
1358 | else | |
1359 | < | myNPTi = new NPTi<RealIntegrator>(&(info[k]), the_ff); |
1359 | > | myNPTi = new NPTi<NPT<RealIntegrator> >(&(info[k]), the_ff); |
1360 | ||
1361 | myNPTi->setTargetTemp(globals->getTargetTemp()); | |
1362 | ||
# | Line 1419 | Line 1389 | void SimSetup::makeIntegrator(void){ | |
1389 | painCave.isFatal = 1; | |
1390 | simError(); | |
1391 | } | |
1392 | + | |
1393 | + | info->the_integrator = myNPTi; |
1394 | break; | |
1395 | ||
1396 | case NPTf_ENS: | |
1397 | if (globals->haveZconstraints()){ | |
1398 | setupZConstraint(info[k]); | |
1399 | < | myNPTf = new ZConstraint<NPTf<RealIntegrator> >(&(info[k]), the_ff); |
1399 | > | myNPTf = new ZConstraint<NPTf<NPT <RealIntegrator> > >(&(info[k]), the_ff); |
1400 | } | |
1401 | else | |
1402 | < | myNPTf = new NPTf<RealIntegrator>(&(info[k]), the_ff); |
1402 | > | myNPTf = new NPTf<NPT <RealIntegrator> >(&(info[k]), the_ff); |
1403 | ||
1404 | myNPTf->setTargetTemp(globals->getTargetTemp()); | |
1405 | ||
# | Line 1443 | Line 1415 | void SimSetup::makeIntegrator(void){ | |
1415 | ||
1416 | if (globals->haveTauThermostat()) | |
1417 | myNPTf->setTauThermostat(globals->getTauThermostat()); | |
1418 | + | |
1419 | else{ | |
1420 | sprintf(painCave.errMsg, | |
1421 | "SimSetup error: If you use an NPT\n" | |
# | Line 1453 | Line 1426 | void SimSetup::makeIntegrator(void){ | |
1426 | ||
1427 | if (globals->haveTauBarostat()) | |
1428 | myNPTf->setTauBarostat(globals->getTauBarostat()); | |
1429 | + | |
1430 | else{ | |
1431 | sprintf(painCave.errMsg, | |
1432 | "SimSetup error: If you use an NPT\n" | |
# | Line 1460 | Line 1434 | void SimSetup::makeIntegrator(void){ | |
1434 | painCave.isFatal = 1; | |
1435 | simError(); | |
1436 | } | |
1437 | + | |
1438 | + | info->the_integrator = myNPTf; |
1439 | break; | |
1440 | ||
1441 | < | case NPTim_ENS: |
1441 | > | case NPTxyz_ENS: |
1442 | if (globals->haveZconstraints()){ | |
1443 | setupZConstraint(info[k]); | |
1444 | < | myNPTim = new ZConstraint<NPTim<RealIntegrator> >(&(info[k]), the_ff); |
1444 | > | myNPTxyz = new ZConstraint<NPTxyz<NPT <RealIntegrator> > >(&(info[k]), the_ff); |
1445 | } | |
1446 | else | |
1447 | < | myNPTim = new NPTim<RealIntegrator>(&(info[k]), the_ff); |
1447 | > | myNPTxyz = new NPTxyz<NPT <RealIntegrator> >(&(info[k]), the_ff); |
1448 | ||
1449 | < | myNPTim->setTargetTemp(globals->getTargetTemp()); |
1449 | > | myNPTxyz->setTargetTemp(globals->getTargetTemp()); |
1450 | ||
1451 | if (globals->haveTargetPressure()) | |
1452 | < | myNPTim->setTargetPressure(globals->getTargetPressure()); |
1452 | > | myNPTxyz->setTargetPressure(globals->getTargetPressure()); |
1453 | else{ | |
1454 | sprintf(painCave.errMsg, | |
1455 | "SimSetup error: If you use a constant pressure\n" | |
1456 | " ensemble, you must set targetPressure in the BASS file.\n"); | |
1457 | painCave.isFatal = 1; | |
1458 | simError(); | |
1459 | < | } |
1459 | > | } |
1460 | ||
1461 | if (globals->haveTauThermostat()) | |
1462 | < | myNPTim->setTauThermostat(globals->getTauThermostat()); |
1462 | > | myNPTxyz->setTauThermostat(globals->getTauThermostat()); |
1463 | else{ | |
1464 | sprintf(painCave.errMsg, | |
1465 | "SimSetup error: If you use an NPT\n" | |
# | Line 1493 | Line 1469 | void SimSetup::makeIntegrator(void){ | |
1469 | } | |
1470 | ||
1471 | if (globals->haveTauBarostat()) | |
1472 | < | myNPTim->setTauBarostat(globals->getTauBarostat()); |
1472 | > | myNPTxyz->setTauBarostat(globals->getTauBarostat()); |
1473 | else{ | |
1474 | sprintf(painCave.errMsg, | |
1475 | "SimSetup error: If you use an NPT\n" | |
1476 | " ensemble, you must set tauBarostat.\n"); | |
1501 | – | painCave.isFatal = 1; |
1502 | – | simError(); |
1503 | – | } |
1504 | – | break; |
1505 | – | |
1506 | – | case NPTfm_ENS: |
1507 | – | if (globals->haveZconstraints()){ |
1508 | – | setupZConstraint(info[k]); |
1509 | – | myNPTfm = new ZConstraint<NPTfm<RealIntegrator> >(&(info[k]), the_ff); |
1510 | – | } |
1511 | – | else |
1512 | – | myNPTfm = new NPTfm<RealIntegrator>(&(info[k]), the_ff); |
1513 | – | |
1514 | – | myNPTfm->setTargetTemp(globals->getTargetTemp()); |
1515 | – | |
1516 | – | if (globals->haveTargetPressure()) |
1517 | – | myNPTfm->setTargetPressure(globals->getTargetPressure()); |
1518 | – | else{ |
1519 | – | sprintf(painCave.errMsg, |
1520 | – | "SimSetup error: If you use a constant pressure\n" |
1521 | – | " ensemble, you must set targetPressure in the BASS file.\n"); |
1477 | painCave.isFatal = 1; | |
1478 | simError(); | |
1479 | } | |
1480 | ||
1481 | < | if (globals->haveTauThermostat()) |
1527 | < | myNPTfm->setTauThermostat(globals->getTauThermostat()); |
1528 | < | else{ |
1529 | < | sprintf(painCave.errMsg, |
1530 | < | "SimSetup error: If you use an NPT\n" |
1531 | < | " ensemble, you must set tauThermostat.\n"); |
1532 | < | painCave.isFatal = 1; |
1533 | < | simError(); |
1534 | < | } |
1535 | < | |
1536 | < | if (globals->haveTauBarostat()) |
1537 | < | myNPTfm->setTauBarostat(globals->getTauBarostat()); |
1538 | < | else{ |
1539 | < | sprintf(painCave.errMsg, |
1540 | < | "SimSetup error: If you use an NPT\n" |
1541 | < | " ensemble, you must set tauBarostat.\n"); |
1542 | < | painCave.isFatal = 1; |
1543 | < | simError(); |
1544 | < | } |
1481 | > | info->the_integrator = myNPTxyz; |
1482 | break; | |
1483 | ||
1484 | default: | |
# | Line 1615 | Line 1552 | void SimSetup::setupZConstraint(SimInfo& theInfo){ | |
1552 | } | |
1553 | theInfo.addProperty(zconsTol); | |
1554 | ||
1555 | < | //set Force Substraction Policy |
1555 | > | //set Force Subtraction Policy |
1556 | StringData* zconsForcePolicy = new StringData(); | |
1557 | zconsForcePolicy->setID(ZCONSFORCEPOLICY_ID); | |
1558 | ||
# | Line 1624 | Line 1561 | void SimSetup::setupZConstraint(SimInfo& theInfo){ | |
1561 | } | |
1562 | else{ | |
1563 | sprintf(painCave.errMsg, | |
1564 | < | "ZConstraint Warning: User does not set force substraction policy, " |
1565 | < | "average force substraction policy is used\n"); |
1564 | > | "ZConstraint Warning: User does not set force Subtraction policy, " |
1565 | > | "PolicyByMass is used\n"); |
1566 | painCave.isFatal = 0; | |
1567 | simError(); | |
1568 | < | zconsForcePolicy->setData("BYNUMBER"); |
1568 | > | zconsForcePolicy->setData("BYMASS"); |
1569 | } | |
1570 | ||
1571 | theInfo.addProperty(zconsForcePolicy); | |
# | Line 1666 | Line 1603 | void SimSetup::setupZConstraint(SimInfo& theInfo){ | |
1603 | zconsParaData->addItem(tempParaItem); | |
1604 | } | |
1605 | ||
1606 | + | //check the uniqueness of index |
1607 | + | if(!zconsParaData->isIndexUnique()){ |
1608 | + | sprintf(painCave.errMsg, |
1609 | + | "ZConstraint Error: molIndex is not unique\n"); |
1610 | + | painCave.isFatal = 1; |
1611 | + | simError(); |
1612 | + | } |
1613 | + | |
1614 | //sort the parameters by index of molecules | |
1615 | zconsParaData->sortByIndex(); | |
1616 | < | |
1616 | > | |
1617 | //push data into siminfo, therefore, we can retrieve later | |
1618 | theInfo.addProperty(zconsParaData); | |
1619 | } |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |