| 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 | 
  | 
 | 
| 121 | 
  | 
  makeOutNames(); | 
| 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, " | 
| 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 | 
< | 
 | 
| 696 | 
< | 
    // get some of the tricky things that may still be in the globals | 
| 697 | 
< | 
 | 
| 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 | 
< | 
    } | 
| 702 | 
> | 
     | 
| 703 | 
  | 
  } | 
| 704 | 
< | 
 | 
| 704 | 
> | 
   | 
| 705 | 
  | 
  //setup seed for random number generator | 
| 706 | 
  | 
  int seedValue; | 
| 707 | 
  | 
 | 
| 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, | 
| 850 | 
  | 
          theEst = globals->getEST(); | 
| 851 | 
  | 
        } | 
| 852 | 
  | 
 | 
| 853 | 
< | 
        info[i].setEcr(theEcr, theEst); | 
| 853 | 
> | 
        info[i].setDefaultEcr(theEcr, theEst); | 
| 854 | 
  | 
      } | 
| 855 | 
  | 
    } | 
| 856 | 
  | 
  } | 
| 877 | 
  | 
    if (worldRank == 0){ | 
| 878 | 
  | 
#endif //is_mpi | 
| 879 | 
  | 
      inName = globals->getInitialConfig(); | 
| 919 | 
– | 
      double* tempDouble = new double[1000000]; | 
| 880 | 
  | 
      fileInit = new InitializeFromFile(inName); | 
| 881 | 
  | 
#ifdef IS_MPI | 
| 882 | 
  | 
    } | 
| 894 | 
  | 
 | 
| 895 | 
  | 
    sprintf(painCave.errMsg, | 
| 896 | 
  | 
            "Cannot intialize a parallel simulation without an initial configuration file.\n"); | 
| 897 | 
< | 
    painCave.isFatal; | 
| 897 | 
> | 
    painCave.isFatal = 1;; | 
| 898 | 
  | 
    simError(); | 
| 899 | 
  | 
 | 
| 900 | 
  | 
#else | 
| 1120 | 
  | 
} | 
| 1121 | 
  | 
 | 
| 1122 | 
  | 
void SimSetup::calcSysValues(void){ | 
| 1123 | 
< | 
  int i, j, k; | 
| 1123 | 
> | 
  int i; | 
| 1124 | 
  | 
 | 
| 1125 | 
  | 
  int* molMembershipArray; | 
| 1126 | 
  | 
 | 
| 1219 | 
  | 
 | 
| 1220 | 
  | 
 | 
| 1221 | 
  | 
void SimSetup::makeSysArrays(void){ | 
| 1222 | 
< | 
  int i, j, k, l; | 
| 1222 | 
> | 
  | 
| 1223 | 
> | 
#ifndef IS_MPI | 
| 1224 | 
> | 
  int k, j; | 
| 1225 | 
> | 
#endif // is_mpi | 
| 1226 | 
> | 
  int i, l; | 
| 1227 | 
  | 
 | 
| 1228 | 
  | 
  Atom** the_atoms; | 
| 1229 | 
  | 
  Molecule* the_molecules; | 
| 1306 | 
  | 
void SimSetup::makeIntegrator(void){ | 
| 1307 | 
  | 
  int k; | 
| 1308 | 
  | 
 | 
| 1309 | 
+ | 
  NVE<RealIntegrator>* myNVE = NULL; | 
| 1310 | 
  | 
  NVT<RealIntegrator>* myNVT = NULL; | 
| 1311 | 
< | 
  NPTi<RealIntegrator>* myNPTi = NULL; | 
| 1312 | 
< | 
  NPTf<RealIntegrator>* myNPTf = NULL; | 
| 1313 | 
< | 
  NPTim<RealIntegrator>* myNPTim = NULL; | 
| 1349 | 
< | 
  NPTfm<RealIntegrator>* myNPTfm = NULL; | 
| 1311 | 
> | 
  NPTi<NPT<RealIntegrator> >* myNPTi = NULL; | 
| 1312 | 
> | 
  NPTf<NPT<RealIntegrator> >* myNPTf = NULL; | 
| 1313 | 
> | 
  NPTxyz<NPT<RealIntegrator> >* myNPTxyz = NULL; | 
| 1314 | 
  | 
   | 
| 1315 | 
  | 
  for (k = 0; k < nInfo; k++){ | 
| 1316 | 
  | 
    switch (ensembleCase){ | 
| 1317 | 
  | 
      case NVE_ENS: | 
| 1318 | 
  | 
        if (globals->haveZconstraints()){ | 
| 1319 | 
  | 
          setupZConstraint(info[k]); | 
| 1320 | 
< | 
          new ZConstraint<NVE<RealIntegrator> >(&(info[k]), the_ff); | 
| 1320 | 
> | 
          myNVE = new ZConstraint<NVE<RealIntegrator> >(&(info[k]), the_ff); | 
| 1321 | 
  | 
        } | 
| 1322 | 
< | 
        else | 
| 1323 | 
< | 
          new NVE<RealIntegrator>(&(info[k]), the_ff); | 
| 1322 | 
> | 
        else{ | 
| 1323 | 
> | 
          myNVE = new NVE<RealIntegrator>(&(info[k]), the_ff); | 
| 1324 | 
> | 
        } | 
| 1325 | 
> | 
         | 
| 1326 | 
> | 
        info->the_integrator = myNVE; | 
| 1327 | 
  | 
        break; | 
| 1328 | 
  | 
 | 
| 1329 | 
  | 
      case NVT_ENS: | 
| 1345 | 
  | 
          painCave.isFatal = 1; | 
| 1346 | 
  | 
          simError(); | 
| 1347 | 
  | 
        } | 
| 1348 | 
+ | 
 | 
| 1349 | 
+ | 
        info->the_integrator = myNVT; | 
| 1350 | 
  | 
        break; | 
| 1351 | 
  | 
 | 
| 1352 | 
  | 
      case NPTi_ENS: | 
| 1353 | 
  | 
        if (globals->haveZconstraints()){ | 
| 1354 | 
  | 
          setupZConstraint(info[k]); | 
| 1355 | 
< | 
          myNPTi = new ZConstraint<NPTi<RealIntegrator> >(&(info[k]), the_ff); | 
| 1355 | 
> | 
          myNPTi = new ZConstraint<NPTi<NPT <RealIntegrator> > >(&(info[k]), the_ff); | 
| 1356 | 
  | 
        } | 
| 1357 | 
  | 
        else | 
| 1358 | 
< | 
          myNPTi = new NPTi<RealIntegrator>(&(info[k]), the_ff); | 
| 1358 | 
> | 
          myNPTi = new NPTi<NPT<RealIntegrator> >(&(info[k]), the_ff); | 
| 1359 | 
  | 
 | 
| 1360 | 
  | 
        myNPTi->setTargetTemp(globals->getTargetTemp()); | 
| 1361 | 
  | 
 | 
| 1388 | 
  | 
          painCave.isFatal = 1; | 
| 1389 | 
  | 
          simError(); | 
| 1390 | 
  | 
        } | 
| 1391 | 
+ | 
 | 
| 1392 | 
+ | 
        info->the_integrator = myNPTi; | 
| 1393 | 
  | 
        break; | 
| 1394 | 
  | 
 | 
| 1395 | 
  | 
      case NPTf_ENS: | 
| 1396 | 
  | 
        if (globals->haveZconstraints()){ | 
| 1397 | 
  | 
          setupZConstraint(info[k]); | 
| 1398 | 
< | 
          myNPTf = new ZConstraint<NPTf<RealIntegrator> >(&(info[k]), the_ff); | 
| 1398 | 
> | 
          myNPTf = new ZConstraint<NPTf<NPT <RealIntegrator> > >(&(info[k]), the_ff); | 
| 1399 | 
  | 
        } | 
| 1400 | 
  | 
        else | 
| 1401 | 
< | 
          myNPTf = new NPTf<RealIntegrator>(&(info[k]), the_ff); | 
| 1401 | 
> | 
          myNPTf = new NPTf<NPT <RealIntegrator> >(&(info[k]), the_ff); | 
| 1402 | 
  | 
 | 
| 1403 | 
  | 
        myNPTf->setTargetTemp(globals->getTargetTemp()); | 
| 1404 | 
  | 
 | 
| 1431 | 
  | 
          painCave.isFatal = 1; | 
| 1432 | 
  | 
          simError(); | 
| 1433 | 
  | 
        } | 
| 1434 | 
+ | 
 | 
| 1435 | 
+ | 
        info->the_integrator = myNPTf; | 
| 1436 | 
  | 
        break; | 
| 1437 | 
  | 
 | 
| 1438 | 
< | 
      case NPTim_ENS: | 
| 1438 | 
> | 
      case NPTxyz_ENS: | 
| 1439 | 
  | 
        if (globals->haveZconstraints()){ | 
| 1440 | 
  | 
          setupZConstraint(info[k]); | 
| 1441 | 
< | 
          myNPTim = new ZConstraint<NPTim<RealIntegrator> >(&(info[k]), the_ff); | 
| 1441 | 
> | 
          myNPTxyz = new ZConstraint<NPTxyz<NPT <RealIntegrator> > >(&(info[k]), the_ff); | 
| 1442 | 
  | 
        } | 
| 1443 | 
  | 
        else | 
| 1444 | 
< | 
          myNPTim = new NPTim<RealIntegrator>(&(info[k]), the_ff); | 
| 1444 | 
> | 
          myNPTxyz = new NPTxyz<NPT <RealIntegrator> >(&(info[k]), the_ff); | 
| 1445 | 
  | 
 | 
| 1446 | 
< | 
        myNPTim->setTargetTemp(globals->getTargetTemp()); | 
| 1446 | 
> | 
        myNPTxyz->setTargetTemp(globals->getTargetTemp()); | 
| 1447 | 
  | 
 | 
| 1448 | 
  | 
        if (globals->haveTargetPressure()) | 
| 1449 | 
< | 
          myNPTim->setTargetPressure(globals->getTargetPressure()); | 
| 1449 | 
> | 
          myNPTxyz->setTargetPressure(globals->getTargetPressure()); | 
| 1450 | 
  | 
        else{ | 
| 1451 | 
  | 
          sprintf(painCave.errMsg, | 
| 1452 | 
  | 
                  "SimSetup error: If you use a constant pressure\n" | 
| 1453 | 
  | 
                  "    ensemble, you must set targetPressure in the BASS file.\n"); | 
| 1454 | 
  | 
          painCave.isFatal = 1; | 
| 1455 | 
  | 
          simError(); | 
| 1456 | 
< | 
        } | 
| 1456 | 
> | 
        }     | 
| 1457 | 
  | 
 | 
| 1458 | 
  | 
        if (globals->haveTauThermostat()) | 
| 1459 | 
< | 
          myNPTim->setTauThermostat(globals->getTauThermostat()); | 
| 1459 | 
> | 
          myNPTxyz->setTauThermostat(globals->getTauThermostat()); | 
| 1460 | 
  | 
        else{ | 
| 1461 | 
  | 
          sprintf(painCave.errMsg, | 
| 1462 | 
  | 
                  "SimSetup error: If you use an NPT\n" | 
| 1466 | 
  | 
        } | 
| 1467 | 
  | 
 | 
| 1468 | 
  | 
        if (globals->haveTauBarostat()) | 
| 1469 | 
< | 
          myNPTim->setTauBarostat(globals->getTauBarostat()); | 
| 1469 | 
> | 
          myNPTxyz->setTauBarostat(globals->getTauBarostat()); | 
| 1470 | 
  | 
        else{ | 
| 1471 | 
  | 
          sprintf(painCave.errMsg, | 
| 1472 | 
  | 
                  "SimSetup error: If you use an NPT\n" | 
| 1474 | 
  | 
          painCave.isFatal = 1; | 
| 1475 | 
  | 
          simError(); | 
| 1476 | 
  | 
        } | 
| 1504 | 
– | 
        break; | 
| 1477 | 
  | 
 | 
| 1478 | 
< | 
      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"); | 
| 1522 | 
< | 
          painCave.isFatal = 1; | 
| 1523 | 
< | 
          simError(); | 
| 1524 | 
< | 
        } | 
| 1525 | 
< | 
 | 
| 1526 | 
< | 
        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 | 
< | 
        } | 
| 1478 | 
> | 
        info->the_integrator = myNPTxyz; | 
| 1479 | 
  | 
        break; | 
| 1480 | 
  | 
 | 
| 1481 | 
  | 
      default: |