| 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 | 
  | 
 | 
| 103 | 
  | 
 | 
| 104 | 
  | 
  sysObjectsCreation(); | 
| 105 | 
  | 
 | 
| 105 | 
– | 
  // check on the post processing info | 
| 106 | 
– | 
 | 
| 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 | 
  | 
 | 
| 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, " | 
| 696 | 
  | 
    } | 
| 697 | 
  | 
 | 
| 698 | 
  | 
    // check for the temperature set flag | 
| 699 | 
< | 
 | 
| 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 | 
| 703 | 
> | 
    // check for the extended State init | 
| 704 | 
  | 
 | 
| 705 | 
< | 
    double boxVector[3]; | 
| 706 | 
< | 
    if (globals->haveBox()){ | 
| 707 | 
< | 
      boxVector[0] = globals->getBox(); | 
| 707 | 
< | 
      boxVector[1] = globals->getBox(); | 
| 708 | 
< | 
      boxVector[2] = globals->getBox(); | 
| 709 | 
< | 
 | 
| 710 | 
< | 
      info[i].setBox(boxVector); | 
| 711 | 
< | 
    } | 
| 712 | 
< | 
    else if (globals->haveDensity()){ | 
| 713 | 
< | 
      double vol; | 
| 714 | 
< | 
      vol = (double) tot_nmol / globals->getDensity(); | 
| 715 | 
< | 
      boxVector[0] = pow(vol, (1.0 / 3.0)); | 
| 716 | 
< | 
      boxVector[1] = boxVector[0]; | 
| 717 | 
< | 
      boxVector[2] = boxVector[0]; | 
| 718 | 
< | 
 | 
| 719 | 
< | 
      info[i].setBox(boxVector); | 
| 720 | 
< | 
    } | 
| 721 | 
< | 
    else{ | 
| 722 | 
< | 
      if (!globals->haveBoxX()){ | 
| 723 | 
< | 
        sprintf(painCave.errMsg, | 
| 724 | 
< | 
                "SimSetup error, no periodic BoxX size given.\n"); | 
| 725 | 
< | 
        painCave.isFatal = 1; | 
| 726 | 
< | 
        simError(); | 
| 727 | 
< | 
      } | 
| 728 | 
< | 
      boxVector[0] = globals->getBoxX(); | 
| 729 | 
< | 
 | 
| 730 | 
< | 
      if (!globals->haveBoxY()){ | 
| 731 | 
< | 
        sprintf(painCave.errMsg, | 
| 732 | 
< | 
                "SimSetup error, no periodic BoxY size given.\n"); | 
| 733 | 
< | 
        painCave.isFatal = 1; | 
| 734 | 
< | 
        simError(); | 
| 735 | 
< | 
      } | 
| 736 | 
< | 
      boxVector[1] = globals->getBoxY(); | 
| 737 | 
< | 
 | 
| 738 | 
< | 
      if (!globals->haveBoxZ()){ | 
| 739 | 
< | 
        sprintf(painCave.errMsg, | 
| 740 | 
< | 
                "SimSetup error, no periodic BoxZ size given.\n"); | 
| 741 | 
< | 
        painCave.isFatal = 1; | 
| 742 | 
< | 
        simError(); | 
| 743 | 
< | 
      } | 
| 744 | 
< | 
      boxVector[2] = globals->getBoxZ(); | 
| 745 | 
< | 
 | 
| 746 | 
< | 
      info[i].setBox(boxVector); | 
| 747 | 
< | 
    } | 
| 705 | 
> | 
    info[i].useInitXSstate = globals->getUseInitXSstate(); | 
| 706 | 
> | 
    info[i].orthoTolerance = globals->getOrthoBoxTolerance(); | 
| 707 | 
> | 
     | 
| 708 | 
  | 
  } | 
| 709 | 
< | 
 | 
| 709 | 
> | 
   | 
| 710 | 
  | 
  //setup seed for random number generator | 
| 711 | 
  | 
  int seedValue; | 
| 712 | 
  | 
 | 
| 810 | 
  | 
        theEst = globals->getEST(); | 
| 811 | 
  | 
      } | 
| 812 | 
  | 
 | 
| 813 | 
< | 
      info[i].setEcr(theEcr, theEst); | 
| 813 | 
> | 
      info[i].setDefaultEcr(theEcr, theEst); | 
| 814 | 
  | 
 | 
| 815 | 
  | 
      if (!globals->haveDielectric()){ | 
| 816 | 
  | 
        sprintf(painCave.errMsg, | 
| 855 | 
  | 
          theEst = globals->getEST(); | 
| 856 | 
  | 
        } | 
| 857 | 
  | 
 | 
| 858 | 
< | 
        info[i].setEcr(theEcr, theEst); | 
| 858 | 
> | 
        info[i].setDefaultEcr(theEcr, theEst); | 
| 859 | 
  | 
      } | 
| 860 | 
  | 
    } | 
| 861 | 
+ | 
 | 
| 862 | 
+ | 
    if( !initSuspend ) | 
| 863 | 
+ | 
      info[i].checkCutOffs(); | 
| 864 | 
  | 
  } | 
| 865 | 
  | 
 | 
| 866 | 
  | 
#ifdef IS_MPI | 
| 885 | 
  | 
    if (worldRank == 0){ | 
| 886 | 
  | 
#endif //is_mpi | 
| 887 | 
  | 
      inName = globals->getInitialConfig(); | 
| 925 | 
– | 
      double* tempDouble = new double[1000000]; | 
| 888 | 
  | 
      fileInit = new InitializeFromFile(inName); | 
| 889 | 
  | 
#ifdef IS_MPI | 
| 890 | 
  | 
    } | 
| 902 | 
  | 
 | 
| 903 | 
  | 
    sprintf(painCave.errMsg, | 
| 904 | 
  | 
            "Cannot intialize a parallel simulation without an initial configuration file.\n"); | 
| 905 | 
< | 
    painCave.isFatal; | 
| 905 | 
> | 
    painCave.isFatal = 1;; | 
| 906 | 
  | 
    simError(); | 
| 907 | 
  | 
 | 
| 908 | 
  | 
#else | 
| 1128 | 
  | 
} | 
| 1129 | 
  | 
 | 
| 1130 | 
  | 
void SimSetup::calcSysValues(void){ | 
| 1131 | 
< | 
  int i, j, k; | 
| 1131 | 
> | 
  int i; | 
| 1132 | 
  | 
 | 
| 1133 | 
  | 
  int* molMembershipArray; | 
| 1134 | 
  | 
 | 
| 1227 | 
  | 
 | 
| 1228 | 
  | 
 | 
| 1229 | 
  | 
void SimSetup::makeSysArrays(void){ | 
| 1230 | 
< | 
  int i, j, k, l; | 
| 1230 | 
> | 
  | 
| 1231 | 
> | 
#ifndef IS_MPI | 
| 1232 | 
> | 
  int k, j; | 
| 1233 | 
> | 
#endif // is_mpi | 
| 1234 | 
> | 
  int i, l; | 
| 1235 | 
  | 
 | 
| 1236 | 
  | 
  Atom** the_atoms; | 
| 1237 | 
  | 
  Molecule* the_molecules; | 
| 1314 | 
  | 
void SimSetup::makeIntegrator(void){ | 
| 1315 | 
  | 
  int k; | 
| 1316 | 
  | 
 | 
| 1317 | 
+ | 
  NVE<RealIntegrator>* myNVE = NULL; | 
| 1318 | 
  | 
  NVT<RealIntegrator>* myNVT = NULL; | 
| 1319 | 
  | 
  NPTi<NPT<RealIntegrator> >* myNPTi = NULL; | 
| 1320 | 
< | 
  NPTf<RealIntegrator>* myNPTf = NULL; | 
| 1321 | 
< | 
  NPTim<RealIntegrator>* myNPTim = NULL; | 
| 1355 | 
< | 
  NPTfm<RealIntegrator>* myNPTfm = NULL; | 
| 1320 | 
> | 
  NPTf<NPT<RealIntegrator> >* myNPTf = NULL; | 
| 1321 | 
> | 
  NPTxyz<NPT<RealIntegrator> >* myNPTxyz = NULL; | 
| 1322 | 
  | 
   | 
| 1323 | 
  | 
  for (k = 0; k < nInfo; k++){ | 
| 1324 | 
  | 
    switch (ensembleCase){ | 
| 1325 | 
  | 
      case NVE_ENS: | 
| 1326 | 
  | 
        if (globals->haveZconstraints()){ | 
| 1327 | 
  | 
          setupZConstraint(info[k]); | 
| 1328 | 
< | 
          new ZConstraint<NVE<RealIntegrator> >(&(info[k]), the_ff); | 
| 1328 | 
> | 
          myNVE = new ZConstraint<NVE<RealIntegrator> >(&(info[k]), the_ff); | 
| 1329 | 
  | 
        } | 
| 1330 | 
< | 
        else | 
| 1331 | 
< | 
          new NVE<RealIntegrator>(&(info[k]), the_ff); | 
| 1330 | 
> | 
        else{ | 
| 1331 | 
> | 
          myNVE = new NVE<RealIntegrator>(&(info[k]), the_ff); | 
| 1332 | 
> | 
        } | 
| 1333 | 
> | 
         | 
| 1334 | 
> | 
        info->the_integrator = myNVE; | 
| 1335 | 
  | 
        break; | 
| 1336 | 
  | 
 | 
| 1337 | 
  | 
      case NVT_ENS: | 
| 1353 | 
  | 
          painCave.isFatal = 1; | 
| 1354 | 
  | 
          simError(); | 
| 1355 | 
  | 
        } | 
| 1356 | 
+ | 
 | 
| 1357 | 
+ | 
        info->the_integrator = myNVT; | 
| 1358 | 
  | 
        break; | 
| 1359 | 
  | 
 | 
| 1360 | 
  | 
      case NPTi_ENS: | 
| 1396 | 
  | 
          painCave.isFatal = 1; | 
| 1397 | 
  | 
          simError(); | 
| 1398 | 
  | 
        } | 
| 1399 | 
+ | 
 | 
| 1400 | 
+ | 
        info->the_integrator = myNPTi; | 
| 1401 | 
  | 
        break; | 
| 1402 | 
  | 
 | 
| 1403 | 
  | 
      case NPTf_ENS: | 
| 1404 | 
  | 
        if (globals->haveZconstraints()){ | 
| 1405 | 
  | 
          setupZConstraint(info[k]); | 
| 1406 | 
< | 
          myNPTf = new ZConstraint<NPTf<RealIntegrator> >(&(info[k]), the_ff); | 
| 1406 | 
> | 
          myNPTf = new ZConstraint<NPTf<NPT <RealIntegrator> > >(&(info[k]), the_ff); | 
| 1407 | 
  | 
        } | 
| 1408 | 
  | 
        else | 
| 1409 | 
< | 
          myNPTf = new NPTf<RealIntegrator>(&(info[k]), the_ff); | 
| 1409 | 
> | 
          myNPTf = new NPTf<NPT <RealIntegrator> >(&(info[k]), the_ff); | 
| 1410 | 
  | 
 | 
| 1411 | 
  | 
        myNPTf->setTargetTemp(globals->getTargetTemp()); | 
| 1412 | 
  | 
 | 
| 1422 | 
  | 
 | 
| 1423 | 
  | 
        if (globals->haveTauThermostat()) | 
| 1424 | 
  | 
          myNPTf->setTauThermostat(globals->getTauThermostat()); | 
| 1425 | 
+ | 
 | 
| 1426 | 
  | 
        else{ | 
| 1427 | 
  | 
          sprintf(painCave.errMsg, | 
| 1428 | 
  | 
                  "SimSetup error: If you use an NPT\n" | 
| 1433 | 
  | 
 | 
| 1434 | 
  | 
        if (globals->haveTauBarostat()) | 
| 1435 | 
  | 
          myNPTf->setTauBarostat(globals->getTauBarostat()); | 
| 1436 | 
+ | 
 | 
| 1437 | 
  | 
        else{ | 
| 1438 | 
  | 
          sprintf(painCave.errMsg, | 
| 1439 | 
  | 
                  "SimSetup error: If you use an NPT\n" | 
| 1441 | 
  | 
          painCave.isFatal = 1; | 
| 1442 | 
  | 
          simError(); | 
| 1443 | 
  | 
        } | 
| 1444 | 
+ | 
 | 
| 1445 | 
+ | 
        info->the_integrator = myNPTf; | 
| 1446 | 
  | 
        break; | 
| 1447 | 
  | 
 | 
| 1448 | 
< | 
      case NPTim_ENS: | 
| 1448 | 
> | 
      case NPTxyz_ENS: | 
| 1449 | 
  | 
        if (globals->haveZconstraints()){ | 
| 1450 | 
  | 
          setupZConstraint(info[k]); | 
| 1451 | 
< | 
          myNPTim = new ZConstraint<NPTim<RealIntegrator> >(&(info[k]), the_ff); | 
| 1451 | 
> | 
          myNPTxyz = new ZConstraint<NPTxyz<NPT <RealIntegrator> > >(&(info[k]), the_ff); | 
| 1452 | 
  | 
        } | 
| 1453 | 
  | 
        else | 
| 1454 | 
< | 
          myNPTim = new NPTim<RealIntegrator>(&(info[k]), the_ff); | 
| 1454 | 
> | 
          myNPTxyz = new NPTxyz<NPT <RealIntegrator> >(&(info[k]), the_ff); | 
| 1455 | 
  | 
 | 
| 1456 | 
< | 
        myNPTim->setTargetTemp(globals->getTargetTemp()); | 
| 1456 | 
> | 
        myNPTxyz->setTargetTemp(globals->getTargetTemp()); | 
| 1457 | 
  | 
 | 
| 1458 | 
  | 
        if (globals->haveTargetPressure()) | 
| 1459 | 
< | 
          myNPTim->setTargetPressure(globals->getTargetPressure()); | 
| 1459 | 
> | 
          myNPTxyz->setTargetPressure(globals->getTargetPressure()); | 
| 1460 | 
  | 
        else{ | 
| 1461 | 
  | 
          sprintf(painCave.errMsg, | 
| 1462 | 
  | 
                  "SimSetup error: If you use a constant pressure\n" | 
| 1463 | 
  | 
                  "    ensemble, you must set targetPressure in the BASS file.\n"); | 
| 1464 | 
  | 
          painCave.isFatal = 1; | 
| 1465 | 
  | 
          simError(); | 
| 1466 | 
< | 
        } | 
| 1466 | 
> | 
        }     | 
| 1467 | 
  | 
 | 
| 1468 | 
  | 
        if (globals->haveTauThermostat()) | 
| 1469 | 
< | 
          myNPTim->setTauThermostat(globals->getTauThermostat()); | 
| 1469 | 
> | 
          myNPTxyz->setTauThermostat(globals->getTauThermostat()); | 
| 1470 | 
  | 
        else{ | 
| 1471 | 
  | 
          sprintf(painCave.errMsg, | 
| 1472 | 
  | 
                  "SimSetup error: If you use an NPT\n" | 
| 1476 | 
  | 
        } | 
| 1477 | 
  | 
 | 
| 1478 | 
  | 
        if (globals->haveTauBarostat()) | 
| 1479 | 
< | 
          myNPTim->setTauBarostat(globals->getTauBarostat()); | 
| 1479 | 
> | 
          myNPTxyz->setTauBarostat(globals->getTauBarostat()); | 
| 1480 | 
  | 
        else{ | 
| 1481 | 
  | 
          sprintf(painCave.errMsg, | 
| 1482 | 
  | 
                  "SimSetup error: If you use an NPT\n" | 
| 1484 | 
  | 
          painCave.isFatal = 1; | 
| 1485 | 
  | 
          simError(); | 
| 1486 | 
  | 
        } | 
| 1510 | 
– | 
        break; | 
| 1487 | 
  | 
 | 
| 1488 | 
< | 
      case NPTfm_ENS: | 
| 1513 | 
< | 
        if (globals->haveZconstraints()){ | 
| 1514 | 
< | 
          setupZConstraint(info[k]); | 
| 1515 | 
< | 
          myNPTfm = new ZConstraint<NPTfm<RealIntegrator> >(&(info[k]), the_ff); | 
| 1516 | 
< | 
        } | 
| 1517 | 
< | 
        else | 
| 1518 | 
< | 
          myNPTfm = new NPTfm<RealIntegrator>(&(info[k]), the_ff); | 
| 1519 | 
< | 
 | 
| 1520 | 
< | 
        myNPTfm->setTargetTemp(globals->getTargetTemp()); | 
| 1521 | 
< | 
 | 
| 1522 | 
< | 
        if (globals->haveTargetPressure()) | 
| 1523 | 
< | 
          myNPTfm->setTargetPressure(globals->getTargetPressure()); | 
| 1524 | 
< | 
        else{ | 
| 1525 | 
< | 
          sprintf(painCave.errMsg, | 
| 1526 | 
< | 
                  "SimSetup error: If you use a constant pressure\n" | 
| 1527 | 
< | 
                  "    ensemble, you must set targetPressure in the BASS file.\n"); | 
| 1528 | 
< | 
          painCave.isFatal = 1; | 
| 1529 | 
< | 
          simError(); | 
| 1530 | 
< | 
        } | 
| 1531 | 
< | 
 | 
| 1532 | 
< | 
        if (globals->haveTauThermostat()) | 
| 1533 | 
< | 
          myNPTfm->setTauThermostat(globals->getTauThermostat()); | 
| 1534 | 
< | 
        else{ | 
| 1535 | 
< | 
          sprintf(painCave.errMsg, | 
| 1536 | 
< | 
                  "SimSetup error: If you use an NPT\n" | 
| 1537 | 
< | 
                  "    ensemble, you must set tauThermostat.\n"); | 
| 1538 | 
< | 
          painCave.isFatal = 1; | 
| 1539 | 
< | 
          simError(); | 
| 1540 | 
< | 
        } | 
| 1541 | 
< | 
 | 
| 1542 | 
< | 
        if (globals->haveTauBarostat()) | 
| 1543 | 
< | 
          myNPTfm->setTauBarostat(globals->getTauBarostat()); | 
| 1544 | 
< | 
        else{ | 
| 1545 | 
< | 
          sprintf(painCave.errMsg, | 
| 1546 | 
< | 
                  "SimSetup error: If you use an NPT\n" | 
| 1547 | 
< | 
                  "    ensemble, you must set tauBarostat.\n"); | 
| 1548 | 
< | 
          painCave.isFatal = 1; | 
| 1549 | 
< | 
          simError(); | 
| 1550 | 
< | 
        } | 
| 1488 | 
> | 
        info->the_integrator = myNPTxyz; | 
| 1489 | 
  | 
        break; | 
| 1490 | 
  | 
 | 
| 1491 | 
  | 
      default: |