| 4 | 
  | 
#include <cmath> | 
| 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 | 
| 91 | 
  | 
#endif // is_mpi | 
| 92 | 
  | 
 | 
| 93 | 
  | 
void SimSetup::createSim(void){ | 
| 95 | 
– | 
  int i, j, k, globalAtomIndex; | 
| 94 | 
  | 
 | 
| 95 | 
  | 
  // gather all of the information from the Bass file | 
| 96 | 
  | 
 | 
| 108 | 
  | 
 | 
| 109 | 
  | 
  if (!isInfoArray){ | 
| 110 | 
  | 
    initSystemCoords(); | 
| 111 | 
+ | 
 | 
| 112 | 
+ | 
    if( !(globals->getUseInitTime()) ) | 
| 113 | 
+ | 
      info[0].currentTime = 0.0; | 
| 114 | 
  | 
  }   | 
| 115 | 
  | 
 | 
| 116 | 
  | 
  // make the output filenames | 
| 132 | 
  | 
 | 
| 133 | 
  | 
 | 
| 134 | 
  | 
void SimSetup::makeMolecules(void){ | 
| 135 | 
< | 
  int k, l; | 
| 135 | 
> | 
  int k; | 
| 136 | 
  | 
  int i, j, exI, exJ, tempEx, stampID, atomOffset, excludeOffset; | 
| 137 | 
  | 
  molInit molInfo; | 
| 138 | 
  | 
  DirectionalAtom* dAtom; | 
| 554 | 
  | 
 | 
| 555 | 
  | 
 | 
| 556 | 
  | 
void SimSetup::gatherInfo(void){ | 
| 557 | 
< | 
  int i, j, k; | 
| 557 | 
> | 
  int i; | 
| 558 | 
  | 
 | 
| 559 | 
  | 
  ensembleCase = -1; | 
| 560 | 
  | 
  ffCase = -1; | 
| 604 | 
  | 
  } | 
| 605 | 
  | 
  else if (!strcasecmp(ensemble, "NPTf")){ | 
| 606 | 
  | 
    ensembleCase = NPTf_ENS; | 
| 606 | 
– | 
  } | 
| 607 | 
– | 
  else if (!strcasecmp(ensemble, "NPTim")){ | 
| 608 | 
– | 
    ensembleCase = NPTim_ENS; | 
| 607 | 
  | 
  } | 
| 608 | 
< | 
  else if (!strcasecmp(ensemble, "NPTfm")){ | 
| 609 | 
< | 
    ensembleCase = NPTfm_ENS; | 
| 608 | 
> | 
  else if (!strcasecmp(ensemble, "NPTxyz")){ | 
| 609 | 
> | 
    ensembleCase = NPTxyz_ENS; | 
| 610 | 
  | 
  } | 
| 611 | 
  | 
  else{ | 
| 612 | 
  | 
    sprintf(painCave.errMsg, | 
| 920 | 
  | 
    if (worldRank == 0){ | 
| 921 | 
  | 
#endif //is_mpi | 
| 922 | 
  | 
      inName = globals->getInitialConfig(); | 
| 925 | 
– | 
      double* tempDouble = new double[1000000]; | 
| 923 | 
  | 
      fileInit = new InitializeFromFile(inName); | 
| 924 | 
  | 
#ifdef IS_MPI | 
| 925 | 
  | 
    } | 
| 937 | 
  | 
 | 
| 938 | 
  | 
    sprintf(painCave.errMsg, | 
| 939 | 
  | 
            "Cannot intialize a parallel simulation without an initial configuration file.\n"); | 
| 940 | 
< | 
    painCave.isFatal; | 
| 940 | 
> | 
    painCave.isFatal = 1;; | 
| 941 | 
  | 
    simError(); | 
| 942 | 
  | 
 | 
| 943 | 
  | 
#else | 
| 1163 | 
  | 
} | 
| 1164 | 
  | 
 | 
| 1165 | 
  | 
void SimSetup::calcSysValues(void){ | 
| 1166 | 
< | 
  int i, j, k; | 
| 1166 | 
> | 
  int i; | 
| 1167 | 
  | 
 | 
| 1168 | 
  | 
  int* molMembershipArray; | 
| 1169 | 
  | 
 | 
| 1262 | 
  | 
 | 
| 1263 | 
  | 
 | 
| 1264 | 
  | 
void SimSetup::makeSysArrays(void){ | 
| 1265 | 
< | 
  int i, j, k, l; | 
| 1265 | 
> | 
  | 
| 1266 | 
> | 
#ifndef IS_MPI | 
| 1267 | 
> | 
  int k, j; | 
| 1268 | 
> | 
#endif // is_mpi | 
| 1269 | 
> | 
  int i, l; | 
| 1270 | 
  | 
 | 
| 1271 | 
  | 
  Atom** the_atoms; | 
| 1272 | 
  | 
  Molecule* the_molecules; | 
| 1349 | 
  | 
void SimSetup::makeIntegrator(void){ | 
| 1350 | 
  | 
  int k; | 
| 1351 | 
  | 
 | 
| 1352 | 
+ | 
  NVE<RealIntegrator>* myNVE = NULL; | 
| 1353 | 
  | 
  NVT<RealIntegrator>* myNVT = NULL; | 
| 1354 | 
< | 
  NPTi<RealIntegrator>* myNPTi = NULL; | 
| 1355 | 
< | 
  NPTf<RealIntegrator>* myNPTf = NULL; | 
| 1356 | 
< | 
  NPTim<RealIntegrator>* myNPTim = NULL; | 
| 1355 | 
< | 
  NPTfm<RealIntegrator>* myNPTfm = NULL; | 
| 1354 | 
> | 
  NPTi<NPT<RealIntegrator> >* myNPTi = NULL; | 
| 1355 | 
> | 
  NPTf<NPT<RealIntegrator> >* myNPTf = NULL; | 
| 1356 | 
> | 
  NPTxyz<NPT<RealIntegrator> >* myNPTxyz = NULL; | 
| 1357 | 
  | 
   | 
| 1358 | 
  | 
  for (k = 0; k < nInfo; k++){ | 
| 1359 | 
  | 
    switch (ensembleCase){ | 
| 1360 | 
  | 
      case NVE_ENS: | 
| 1361 | 
  | 
        if (globals->haveZconstraints()){ | 
| 1362 | 
  | 
          setupZConstraint(info[k]); | 
| 1363 | 
< | 
          new ZConstraint<NVE<RealIntegrator> >(&(info[k]), the_ff); | 
| 1363 | 
> | 
          myNVE = new ZConstraint<NVE<RealIntegrator> >(&(info[k]), the_ff); | 
| 1364 | 
  | 
        } | 
| 1365 | 
< | 
        else | 
| 1366 | 
< | 
          new NVE<RealIntegrator>(&(info[k]), the_ff); | 
| 1365 | 
> | 
        else{ | 
| 1366 | 
> | 
          myNVE = new NVE<RealIntegrator>(&(info[k]), the_ff); | 
| 1367 | 
> | 
        } | 
| 1368 | 
> | 
         | 
| 1369 | 
> | 
        info->the_integrator = myNVE; | 
| 1370 | 
  | 
        break; | 
| 1371 | 
  | 
 | 
| 1372 | 
  | 
      case NVT_ENS: | 
| 1388 | 
  | 
          painCave.isFatal = 1; | 
| 1389 | 
  | 
          simError(); | 
| 1390 | 
  | 
        } | 
| 1391 | 
+ | 
 | 
| 1392 | 
+ | 
        info->the_integrator = myNVT; | 
| 1393 | 
  | 
        break; | 
| 1394 | 
  | 
 | 
| 1395 | 
  | 
      case NPTi_ENS: | 
| 1396 | 
  | 
        if (globals->haveZconstraints()){ | 
| 1397 | 
  | 
          setupZConstraint(info[k]); | 
| 1398 | 
< | 
          myNPTi = new ZConstraint<NPTi<RealIntegrator> >(&(info[k]), the_ff); | 
| 1398 | 
> | 
          myNPTi = new ZConstraint<NPTi<NPT <RealIntegrator> > >(&(info[k]), the_ff); | 
| 1399 | 
  | 
        } | 
| 1400 | 
  | 
        else | 
| 1401 | 
< | 
          myNPTi = new NPTi<RealIntegrator>(&(info[k]), the_ff); | 
| 1401 | 
> | 
          myNPTi = new NPTi<NPT<RealIntegrator> >(&(info[k]), the_ff); | 
| 1402 | 
  | 
 | 
| 1403 | 
  | 
        myNPTi->setTargetTemp(globals->getTargetTemp()); | 
| 1404 | 
  | 
 | 
| 1431 | 
  | 
          painCave.isFatal = 1; | 
| 1432 | 
  | 
          simError(); | 
| 1433 | 
  | 
        } | 
| 1434 | 
+ | 
 | 
| 1435 | 
+ | 
        info->the_integrator = myNPTi; | 
| 1436 | 
  | 
        break; | 
| 1437 | 
  | 
 | 
| 1438 | 
  | 
      case NPTf_ENS: | 
| 1439 | 
  | 
        if (globals->haveZconstraints()){ | 
| 1440 | 
  | 
          setupZConstraint(info[k]); | 
| 1441 | 
< | 
          myNPTf = new ZConstraint<NPTf<RealIntegrator> >(&(info[k]), the_ff); | 
| 1441 | 
> | 
          myNPTf = new ZConstraint<NPTf<NPT <RealIntegrator> > >(&(info[k]), the_ff); | 
| 1442 | 
  | 
        } | 
| 1443 | 
  | 
        else | 
| 1444 | 
< | 
          myNPTf = new NPTf<RealIntegrator>(&(info[k]), the_ff); | 
| 1444 | 
> | 
          myNPTf = new NPTf<NPT <RealIntegrator> >(&(info[k]), the_ff); | 
| 1445 | 
  | 
 | 
| 1446 | 
  | 
        myNPTf->setTargetTemp(globals->getTargetTemp()); | 
| 1447 | 
  | 
 | 
| 1474 | 
  | 
          painCave.isFatal = 1; | 
| 1475 | 
  | 
          simError(); | 
| 1476 | 
  | 
        } | 
| 1469 | 
– | 
        break; | 
| 1477 | 
  | 
 | 
| 1478 | 
< | 
      case NPTim_ENS: | 
| 1472 | 
< | 
        if (globals->haveZconstraints()){ | 
| 1473 | 
< | 
          setupZConstraint(info[k]); | 
| 1474 | 
< | 
          myNPTim = new ZConstraint<NPTim<RealIntegrator> >(&(info[k]), the_ff); | 
| 1475 | 
< | 
        } | 
| 1476 | 
< | 
        else | 
| 1477 | 
< | 
          myNPTim = new NPTim<RealIntegrator>(&(info[k]), the_ff); | 
| 1478 | 
< | 
 | 
| 1479 | 
< | 
        myNPTim->setTargetTemp(globals->getTargetTemp()); | 
| 1480 | 
< | 
 | 
| 1481 | 
< | 
        if (globals->haveTargetPressure()) | 
| 1482 | 
< | 
          myNPTim->setTargetPressure(globals->getTargetPressure()); | 
| 1483 | 
< | 
        else{ | 
| 1484 | 
< | 
          sprintf(painCave.errMsg, | 
| 1485 | 
< | 
                  "SimSetup error: If you use a constant pressure\n" | 
| 1486 | 
< | 
                  "    ensemble, you must set targetPressure in the BASS file.\n"); | 
| 1487 | 
< | 
          painCave.isFatal = 1; | 
| 1488 | 
< | 
          simError(); | 
| 1489 | 
< | 
        } | 
| 1490 | 
< | 
 | 
| 1491 | 
< | 
        if (globals->haveTauThermostat()) | 
| 1492 | 
< | 
          myNPTim->setTauThermostat(globals->getTauThermostat()); | 
| 1493 | 
< | 
        else{ | 
| 1494 | 
< | 
          sprintf(painCave.errMsg, | 
| 1495 | 
< | 
                  "SimSetup error: If you use an NPT\n" | 
| 1496 | 
< | 
                  "    ensemble, you must set tauThermostat.\n"); | 
| 1497 | 
< | 
          painCave.isFatal = 1; | 
| 1498 | 
< | 
          simError(); | 
| 1499 | 
< | 
        } | 
| 1500 | 
< | 
 | 
| 1501 | 
< | 
        if (globals->haveTauBarostat()) | 
| 1502 | 
< | 
          myNPTim->setTauBarostat(globals->getTauBarostat()); | 
| 1503 | 
< | 
        else{ | 
| 1504 | 
< | 
          sprintf(painCave.errMsg, | 
| 1505 | 
< | 
                  "SimSetup error: If you use an NPT\n" | 
| 1506 | 
< | 
                  "    ensemble, you must set tauBarostat.\n"); | 
| 1507 | 
< | 
          painCave.isFatal = 1; | 
| 1508 | 
< | 
          simError(); | 
| 1509 | 
< | 
        } | 
| 1478 | 
> | 
        info->the_integrator = myNPTf; | 
| 1479 | 
  | 
        break; | 
| 1480 | 
  | 
 | 
| 1481 | 
< | 
      case NPTfm_ENS: | 
| 1481 | 
> | 
      case NPTxyz_ENS: | 
| 1482 | 
  | 
        if (globals->haveZconstraints()){ | 
| 1483 | 
  | 
          setupZConstraint(info[k]); | 
| 1484 | 
< | 
          myNPTfm = new ZConstraint<NPTfm<RealIntegrator> >(&(info[k]), the_ff); | 
| 1484 | 
> | 
          myNPTxyz = new ZConstraint<NPTxyz<NPT <RealIntegrator> > >(&(info[k]), the_ff); | 
| 1485 | 
  | 
        } | 
| 1486 | 
  | 
        else | 
| 1487 | 
< | 
          myNPTfm = new NPTfm<RealIntegrator>(&(info[k]), the_ff); | 
| 1487 | 
> | 
          myNPTxyz = new NPTxyz<NPT <RealIntegrator> >(&(info[k]), the_ff); | 
| 1488 | 
  | 
 | 
| 1489 | 
< | 
        myNPTfm->setTargetTemp(globals->getTargetTemp()); | 
| 1489 | 
> | 
        myNPTxyz->setTargetTemp(globals->getTargetTemp()); | 
| 1490 | 
  | 
 | 
| 1491 | 
  | 
        if (globals->haveTargetPressure()) | 
| 1492 | 
< | 
          myNPTfm->setTargetPressure(globals->getTargetPressure()); | 
| 1492 | 
> | 
          myNPTxyz->setTargetPressure(globals->getTargetPressure()); | 
| 1493 | 
  | 
        else{ | 
| 1494 | 
  | 
          sprintf(painCave.errMsg, | 
| 1495 | 
  | 
                  "SimSetup error: If you use a constant pressure\n" | 
| 1496 | 
  | 
                  "    ensemble, you must set targetPressure in the BASS file.\n"); | 
| 1497 | 
  | 
          painCave.isFatal = 1; | 
| 1498 | 
  | 
          simError(); | 
| 1499 | 
< | 
        } | 
| 1499 | 
> | 
        }     | 
| 1500 | 
  | 
 | 
| 1501 | 
  | 
        if (globals->haveTauThermostat()) | 
| 1502 | 
< | 
          myNPTfm->setTauThermostat(globals->getTauThermostat()); | 
| 1502 | 
> | 
          myNPTxyz->setTauThermostat(globals->getTauThermostat()); | 
| 1503 | 
  | 
        else{ | 
| 1504 | 
  | 
          sprintf(painCave.errMsg, | 
| 1505 | 
  | 
                  "SimSetup error: If you use an NPT\n" | 
| 1509 | 
  | 
        } | 
| 1510 | 
  | 
 | 
| 1511 | 
  | 
        if (globals->haveTauBarostat()) | 
| 1512 | 
< | 
          myNPTfm->setTauBarostat(globals->getTauBarostat()); | 
| 1512 | 
> | 
          myNPTxyz->setTauBarostat(globals->getTauBarostat()); | 
| 1513 | 
  | 
        else{ | 
| 1514 | 
  | 
          sprintf(painCave.errMsg, | 
| 1515 | 
  | 
                  "SimSetup error: If you use an NPT\n" | 
| 1517 | 
  | 
          painCave.isFatal = 1; | 
| 1518 | 
  | 
          simError(); | 
| 1519 | 
  | 
        } | 
| 1520 | 
+ | 
 | 
| 1521 | 
+ | 
        info->the_integrator = myNPTxyz; | 
| 1522 | 
  | 
        break; | 
| 1523 | 
  | 
 | 
| 1524 | 
  | 
      default: |