| 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" | 
| 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 | 
| 31 | 
  | 
using namespace std; | 
| 32 | 
  | 
 | 
| 33 | 
  | 
SimSetup::SimSetup(){ | 
| 34 | 
+ | 
   | 
| 35 | 
+ | 
  initSuspend = false; | 
| 36 | 
  | 
  isInfoArray = 0; | 
| 37 | 
  | 
  nInfo = 1; | 
| 38 | 
  | 
 | 
| 55 | 
  | 
  info = the_info; | 
| 56 | 
  | 
  nInfo = theNinfo; | 
| 57 | 
  | 
  isInfoArray = 1; | 
| 58 | 
+ | 
  initSuspend = true; | 
| 59 | 
  | 
} | 
| 60 | 
  | 
 | 
| 61 | 
  | 
 | 
| 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 | 
  | 
 | 
| 109 | 
  | 
 | 
| 110 | 
  | 
  // initialize the system coordinates | 
| 111 | 
  | 
 | 
| 112 | 
< | 
  if (!isInfoArray){ | 
| 112 | 
> | 
  if ( !initSuspend ){ | 
| 113 | 
  | 
    initSystemCoords(); | 
| 114 | 
+ | 
 | 
| 115 | 
+ | 
    if( !(globals->getUseInitTime()) ) | 
| 116 | 
+ | 
      info[0].currentTime = 0.0; | 
| 117 | 
  | 
  }   | 
| 118 | 
  | 
 | 
| 119 | 
  | 
  // make the output filenames | 
| 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; | 
| 557 | 
  | 
 | 
| 558 | 
  | 
 | 
| 559 | 
  | 
void SimSetup::gatherInfo(void){ | 
| 560 | 
< | 
  int i, j, k; | 
| 560 | 
> | 
  int i; | 
| 561 | 
  | 
 | 
| 562 | 
  | 
  ensembleCase = -1; | 
| 563 | 
  | 
  ffCase = -1; | 
| 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, " | 
| 687 | 
  | 
 | 
| 688 | 
  | 
    if (globals->haveThermalTime()){ | 
| 689 | 
  | 
      info[i].thermalTime = globals->getThermalTime(); | 
| 690 | 
+ | 
    } | 
| 691 | 
+ | 
 | 
| 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 | 
| 923 | 
  | 
    if (worldRank == 0){ | 
| 924 | 
  | 
#endif //is_mpi | 
| 925 | 
  | 
      inName = globals->getInitialConfig(); | 
| 919 | 
– | 
      double* tempDouble = new double[1000000]; | 
| 926 | 
  | 
      fileInit = new InitializeFromFile(inName); | 
| 927 | 
  | 
#ifdef IS_MPI | 
| 928 | 
  | 
    } | 
| 940 | 
  | 
 | 
| 941 | 
  | 
    sprintf(painCave.errMsg, | 
| 942 | 
  | 
            "Cannot intialize a parallel simulation without an initial configuration file.\n"); | 
| 943 | 
< | 
    painCave.isFatal; | 
| 943 | 
> | 
    painCave.isFatal = 1;; | 
| 944 | 
  | 
    simError(); | 
| 945 | 
  | 
 | 
| 946 | 
  | 
#else | 
| 1166 | 
  | 
} | 
| 1167 | 
  | 
 | 
| 1168 | 
  | 
void SimSetup::calcSysValues(void){ | 
| 1169 | 
< | 
  int i, j, k; | 
| 1169 | 
> | 
  int i; | 
| 1170 | 
  | 
 | 
| 1171 | 
  | 
  int* molMembershipArray; | 
| 1172 | 
  | 
 | 
| 1265 | 
  | 
 | 
| 1266 | 
  | 
 | 
| 1267 | 
  | 
void SimSetup::makeSysArrays(void){ | 
| 1268 | 
< | 
  int i, j, k, l; | 
| 1268 | 
> | 
  | 
| 1269 | 
> | 
#ifndef IS_MPI | 
| 1270 | 
> | 
  int k, j; | 
| 1271 | 
> | 
#endif // is_mpi | 
| 1272 | 
> | 
  int i, l; | 
| 1273 | 
  | 
 | 
| 1274 | 
  | 
  Atom** the_atoms; | 
| 1275 | 
  | 
  Molecule* the_molecules; | 
| 1352 | 
  | 
void SimSetup::makeIntegrator(void){ | 
| 1353 | 
  | 
  int k; | 
| 1354 | 
  | 
 | 
| 1355 | 
+ | 
  NVE<RealIntegrator>* myNVE = NULL; | 
| 1356 | 
  | 
  NVT<RealIntegrator>* myNVT = NULL; | 
| 1357 | 
< | 
  NPTi<RealIntegrator>* myNPTi = NULL; | 
| 1358 | 
< | 
  NPTf<RealIntegrator>* myNPTf = NULL; | 
| 1359 | 
< | 
  NPTim<RealIntegrator>* myNPTim = NULL; | 
| 1349 | 
< | 
  NPTfm<RealIntegrator>* myNPTfm = NULL; | 
| 1357 | 
> | 
  NPTi<NPT<RealIntegrator> >* myNPTi = NULL; | 
| 1358 | 
> | 
  NPTf<NPT<RealIntegrator> >* myNPTf = NULL; | 
| 1359 | 
> | 
  NPTxyz<NPT<RealIntegrator> >* myNPTxyz = NULL; | 
| 1360 | 
  | 
   | 
| 1361 | 
  | 
  for (k = 0; k < nInfo; k++){ | 
| 1362 | 
  | 
    switch (ensembleCase){ | 
| 1363 | 
  | 
      case NVE_ENS: | 
| 1364 | 
  | 
        if (globals->haveZconstraints()){ | 
| 1365 | 
  | 
          setupZConstraint(info[k]); | 
| 1366 | 
< | 
          new ZConstraint<NVE<RealIntegrator> >(&(info[k]), the_ff); | 
| 1366 | 
> | 
          myNVE = new ZConstraint<NVE<RealIntegrator> >(&(info[k]), the_ff); | 
| 1367 | 
  | 
        } | 
| 1368 | 
< | 
        else | 
| 1369 | 
< | 
          new NVE<RealIntegrator>(&(info[k]), the_ff); | 
| 1368 | 
> | 
        else{ | 
| 1369 | 
> | 
          myNVE = new NVE<RealIntegrator>(&(info[k]), the_ff); | 
| 1370 | 
> | 
        } | 
| 1371 | 
> | 
         | 
| 1372 | 
> | 
        info->the_integrator = myNVE; | 
| 1373 | 
  | 
        break; | 
| 1374 | 
  | 
 | 
| 1375 | 
  | 
      case NVT_ENS: | 
| 1391 | 
  | 
          painCave.isFatal = 1; | 
| 1392 | 
  | 
          simError(); | 
| 1393 | 
  | 
        } | 
| 1394 | 
+ | 
 | 
| 1395 | 
+ | 
        info->the_integrator = myNVT; | 
| 1396 | 
  | 
        break; | 
| 1397 | 
  | 
 | 
| 1398 | 
  | 
      case NPTi_ENS: | 
| 1399 | 
  | 
        if (globals->haveZconstraints()){ | 
| 1400 | 
  | 
          setupZConstraint(info[k]); | 
| 1401 | 
< | 
          myNPTi = new ZConstraint<NPTi<RealIntegrator> >(&(info[k]), the_ff); | 
| 1401 | 
> | 
          myNPTi = new ZConstraint<NPTi<NPT <RealIntegrator> > >(&(info[k]), the_ff); | 
| 1402 | 
  | 
        } | 
| 1403 | 
  | 
        else | 
| 1404 | 
< | 
          myNPTi = new NPTi<RealIntegrator>(&(info[k]), the_ff); | 
| 1404 | 
> | 
          myNPTi = new NPTi<NPT<RealIntegrator> >(&(info[k]), the_ff); | 
| 1405 | 
  | 
 | 
| 1406 | 
  | 
        myNPTi->setTargetTemp(globals->getTargetTemp()); | 
| 1407 | 
  | 
 | 
| 1434 | 
  | 
          painCave.isFatal = 1; | 
| 1435 | 
  | 
          simError(); | 
| 1436 | 
  | 
        } | 
| 1437 | 
+ | 
 | 
| 1438 | 
+ | 
        info->the_integrator = myNPTi; | 
| 1439 | 
  | 
        break; | 
| 1440 | 
  | 
 | 
| 1441 | 
  | 
      case NPTf_ENS: | 
| 1442 | 
  | 
        if (globals->haveZconstraints()){ | 
| 1443 | 
  | 
          setupZConstraint(info[k]); | 
| 1444 | 
< | 
          myNPTf = new ZConstraint<NPTf<RealIntegrator> >(&(info[k]), the_ff); | 
| 1444 | 
> | 
          myNPTf = new ZConstraint<NPTf<NPT <RealIntegrator> > >(&(info[k]), the_ff); | 
| 1445 | 
  | 
        } | 
| 1446 | 
  | 
        else | 
| 1447 | 
< | 
          myNPTf = new NPTf<RealIntegrator>(&(info[k]), the_ff); | 
| 1447 | 
> | 
          myNPTf = new NPTf<NPT <RealIntegrator> >(&(info[k]), the_ff); | 
| 1448 | 
  | 
 | 
| 1449 | 
  | 
        myNPTf->setTargetTemp(globals->getTargetTemp()); | 
| 1450 | 
  | 
 | 
| 1474 | 
  | 
          sprintf(painCave.errMsg, | 
| 1475 | 
  | 
                  "SimSetup error: If you use an NPT\n" | 
| 1476 | 
  | 
                  "    ensemble, you must set tauBarostat.\n"); | 
| 1460 | 
– | 
          painCave.isFatal = 1; | 
| 1461 | 
– | 
          simError(); | 
| 1462 | 
– | 
        } | 
| 1463 | 
– | 
        break; | 
| 1464 | 
– | 
 | 
| 1465 | 
– | 
      case NPTim_ENS: | 
| 1466 | 
– | 
        if (globals->haveZconstraints()){ | 
| 1467 | 
– | 
          setupZConstraint(info[k]); | 
| 1468 | 
– | 
          myNPTim = new ZConstraint<NPTim<RealIntegrator> >(&(info[k]), the_ff); | 
| 1469 | 
– | 
        } | 
| 1470 | 
– | 
        else | 
| 1471 | 
– | 
          myNPTim = new NPTim<RealIntegrator>(&(info[k]), the_ff); | 
| 1472 | 
– | 
 | 
| 1473 | 
– | 
        myNPTim->setTargetTemp(globals->getTargetTemp()); | 
| 1474 | 
– | 
 | 
| 1475 | 
– | 
        if (globals->haveTargetPressure()) | 
| 1476 | 
– | 
          myNPTim->setTargetPressure(globals->getTargetPressure()); | 
| 1477 | 
– | 
        else{ | 
| 1478 | 
– | 
          sprintf(painCave.errMsg, | 
| 1479 | 
– | 
                  "SimSetup error: If you use a constant pressure\n" | 
| 1480 | 
– | 
                  "    ensemble, you must set targetPressure in the BASS file.\n"); | 
| 1477 | 
  | 
          painCave.isFatal = 1; | 
| 1478 | 
  | 
          simError(); | 
| 1479 | 
  | 
        } | 
| 1480 | 
  | 
 | 
| 1481 | 
< | 
        if (globals->haveTauThermostat()) | 
| 1486 | 
< | 
          myNPTim->setTauThermostat(globals->getTauThermostat()); | 
| 1487 | 
< | 
        else{ | 
| 1488 | 
< | 
          sprintf(painCave.errMsg, | 
| 1489 | 
< | 
                  "SimSetup error: If you use an NPT\n" | 
| 1490 | 
< | 
                  "    ensemble, you must set tauThermostat.\n"); | 
| 1491 | 
< | 
          painCave.isFatal = 1; | 
| 1492 | 
< | 
          simError(); | 
| 1493 | 
< | 
        } | 
| 1494 | 
< | 
 | 
| 1495 | 
< | 
        if (globals->haveTauBarostat()) | 
| 1496 | 
< | 
          myNPTim->setTauBarostat(globals->getTauBarostat()); | 
| 1497 | 
< | 
        else{ | 
| 1498 | 
< | 
          sprintf(painCave.errMsg, | 
| 1499 | 
< | 
                  "SimSetup error: If you use an NPT\n" | 
| 1500 | 
< | 
                  "    ensemble, you must set tauBarostat.\n"); | 
| 1501 | 
< | 
          painCave.isFatal = 1; | 
| 1502 | 
< | 
          simError(); | 
| 1503 | 
< | 
        } | 
| 1481 | 
> | 
        info->the_integrator = myNPTf; | 
| 1482 | 
  | 
        break; | 
| 1483 | 
  | 
 | 
| 1484 | 
< | 
      case NPTfm_ENS: | 
| 1484 | 
> | 
      case NPTxyz_ENS: | 
| 1485 | 
  | 
        if (globals->haveZconstraints()){ | 
| 1486 | 
  | 
          setupZConstraint(info[k]); | 
| 1487 | 
< | 
          myNPTfm = new ZConstraint<NPTfm<RealIntegrator> >(&(info[k]), the_ff); | 
| 1487 | 
> | 
          myNPTxyz = new ZConstraint<NPTxyz<NPT <RealIntegrator> > >(&(info[k]), the_ff); | 
| 1488 | 
  | 
        } | 
| 1489 | 
  | 
        else | 
| 1490 | 
< | 
          myNPTfm = new NPTfm<RealIntegrator>(&(info[k]), the_ff); | 
| 1490 | 
> | 
          myNPTxyz = new NPTxyz<NPT <RealIntegrator> >(&(info[k]), the_ff); | 
| 1491 | 
  | 
 | 
| 1492 | 
< | 
        myNPTfm->setTargetTemp(globals->getTargetTemp()); | 
| 1492 | 
> | 
        myNPTxyz->setTargetTemp(globals->getTargetTemp()); | 
| 1493 | 
  | 
 | 
| 1494 | 
  | 
        if (globals->haveTargetPressure()) | 
| 1495 | 
< | 
          myNPTfm->setTargetPressure(globals->getTargetPressure()); | 
| 1495 | 
> | 
          myNPTxyz->setTargetPressure(globals->getTargetPressure()); | 
| 1496 | 
  | 
        else{ | 
| 1497 | 
  | 
          sprintf(painCave.errMsg, | 
| 1498 | 
  | 
                  "SimSetup error: If you use a constant pressure\n" | 
| 1499 | 
  | 
                  "    ensemble, you must set targetPressure in the BASS file.\n"); | 
| 1500 | 
  | 
          painCave.isFatal = 1; | 
| 1501 | 
  | 
          simError(); | 
| 1502 | 
< | 
        } | 
| 1502 | 
> | 
        }     | 
| 1503 | 
  | 
 | 
| 1504 | 
  | 
        if (globals->haveTauThermostat()) | 
| 1505 | 
< | 
          myNPTfm->setTauThermostat(globals->getTauThermostat()); | 
| 1505 | 
> | 
          myNPTxyz->setTauThermostat(globals->getTauThermostat()); | 
| 1506 | 
  | 
        else{ | 
| 1507 | 
  | 
          sprintf(painCave.errMsg, | 
| 1508 | 
  | 
                  "SimSetup error: If you use an NPT\n" | 
| 1512 | 
  | 
        } | 
| 1513 | 
  | 
 | 
| 1514 | 
  | 
        if (globals->haveTauBarostat()) | 
| 1515 | 
< | 
          myNPTfm->setTauBarostat(globals->getTauBarostat()); | 
| 1515 | 
> | 
          myNPTxyz->setTauBarostat(globals->getTauBarostat()); | 
| 1516 | 
  | 
        else{ | 
| 1517 | 
  | 
          sprintf(painCave.errMsg, | 
| 1518 | 
  | 
                  "SimSetup error: If you use an NPT\n" | 
| 1520 | 
  | 
          painCave.isFatal = 1; | 
| 1521 | 
  | 
          simError(); | 
| 1522 | 
  | 
        } | 
| 1523 | 
+ | 
 | 
| 1524 | 
+ | 
        info->the_integrator = myNPTxyz; | 
| 1525 | 
  | 
        break; | 
| 1526 | 
  | 
 | 
| 1527 | 
  | 
      default: | 
| 1595 | 
  | 
  } | 
| 1596 | 
  | 
  theInfo.addProperty(zconsTol); | 
| 1597 | 
  | 
 | 
| 1598 | 
< | 
  //set Force Substraction Policy | 
| 1598 | 
> | 
  //set Force Subtraction Policy | 
| 1599 | 
  | 
  StringData* zconsForcePolicy = new StringData(); | 
| 1600 | 
  | 
  zconsForcePolicy->setID(ZCONSFORCEPOLICY_ID); | 
| 1601 | 
  | 
 | 
| 1604 | 
  | 
  } | 
| 1605 | 
  | 
  else{ | 
| 1606 | 
  | 
    sprintf(painCave.errMsg, | 
| 1607 | 
< | 
            "ZConstraint Warning: User does not set force substraction policy, " | 
| 1608 | 
< | 
            "average force substraction policy is used\n"); | 
| 1607 | 
> | 
            "ZConstraint Warning: User does not set force Subtraction policy, " | 
| 1608 | 
> | 
            "PolicyByMass is used\n"); | 
| 1609 | 
  | 
    painCave.isFatal = 0; | 
| 1610 | 
  | 
    simError();  | 
| 1611 | 
< | 
    zconsForcePolicy->setData("BYNUMBER"); | 
| 1611 | 
> | 
    zconsForcePolicy->setData("BYMASS"); | 
| 1612 | 
  | 
  } | 
| 1613 | 
  | 
 | 
| 1614 | 
  | 
  theInfo.addProperty(zconsForcePolicy); | 
| 1646 | 
  | 
    zconsParaData->addItem(tempParaItem); | 
| 1647 | 
  | 
  } | 
| 1648 | 
  | 
 | 
| 1649 | 
+ | 
  //check the uniqueness of index   | 
| 1650 | 
+ | 
  if(!zconsParaData->isIndexUnique()){ | 
| 1651 | 
+ | 
    sprintf(painCave.errMsg, | 
| 1652 | 
+ | 
            "ZConstraint Error: molIndex is not unique\n"); | 
| 1653 | 
+ | 
    painCave.isFatal = 1; | 
| 1654 | 
+ | 
    simError();  | 
| 1655 | 
+ | 
  } | 
| 1656 | 
+ | 
 | 
| 1657 | 
  | 
  //sort the parameters by index of molecules | 
| 1658 | 
  | 
  zconsParaData->sortByIndex(); | 
| 1659 | 
< | 
 | 
| 1659 | 
> | 
   | 
| 1660 | 
  | 
  //push data into siminfo, therefore, we can retrieve later | 
| 1661 | 
  | 
  theInfo.addProperty(zconsParaData); | 
| 1662 | 
  | 
} |