| 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 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 | 
  | 
 | 
| 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 | 
| 607 | 
  | 
  } | 
| 608 | 
  | 
  else if (!strcasecmp(ensemble, "NPTf")){ | 
| 609 | 
  | 
    ensembleCase = NPTf_ENS; | 
| 610 | 
+ | 
  } | 
| 611 | 
+ | 
  else if (!strcasecmp(ensemble, "NPTxyz")){ | 
| 612 | 
+ | 
    ensembleCase = NPTxyz_ENS; | 
| 613 | 
  | 
  } | 
| 614 | 
  | 
  else{ | 
| 615 | 
  | 
    sprintf(painCave.errMsg, | 
| 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 | 
| 694 | 
< | 
 | 
| 695 | 
< | 
    double boxVector[3]; | 
| 696 | 
< | 
    if (globals->haveBox()){ | 
| 697 | 
< | 
      boxVector[0] = globals->getBox(); | 
| 698 | 
< | 
      boxVector[1] = globals->getBox(); | 
| 699 | 
< | 
      boxVector[2] = globals->getBox(); | 
| 703 | 
> | 
    // check for the extended State init | 
| 704 | 
  | 
 | 
| 705 | 
< | 
      info[i].setBox(boxVector); | 
| 706 | 
< | 
    } | 
| 707 | 
< | 
    else if (globals->haveDensity()){ | 
| 704 | 
< | 
      double vol; | 
| 705 | 
< | 
      vol = (double) tot_nmol / globals->getDensity(); | 
| 706 | 
< | 
      boxVector[0] = pow(vol, (1.0 / 3.0)); | 
| 707 | 
< | 
      boxVector[1] = boxVector[0]; | 
| 708 | 
< | 
      boxVector[2] = boxVector[0]; | 
| 709 | 
< | 
 | 
| 710 | 
< | 
      info[i].setBox(boxVector); | 
| 711 | 
< | 
    } | 
| 712 | 
< | 
    else{ | 
| 713 | 
< | 
      if (!globals->haveBoxX()){ | 
| 714 | 
< | 
        sprintf(painCave.errMsg, | 
| 715 | 
< | 
                "SimSetup error, no periodic BoxX size given.\n"); | 
| 716 | 
< | 
        painCave.isFatal = 1; | 
| 717 | 
< | 
        simError(); | 
| 718 | 
< | 
      } | 
| 719 | 
< | 
      boxVector[0] = globals->getBoxX(); | 
| 720 | 
< | 
 | 
| 721 | 
< | 
      if (!globals->haveBoxY()){ | 
| 722 | 
< | 
        sprintf(painCave.errMsg, | 
| 723 | 
< | 
                "SimSetup error, no periodic BoxY size given.\n"); | 
| 724 | 
< | 
        painCave.isFatal = 1; | 
| 725 | 
< | 
        simError(); | 
| 726 | 
< | 
      } | 
| 727 | 
< | 
      boxVector[1] = globals->getBoxY(); | 
| 728 | 
< | 
 | 
| 729 | 
< | 
      if (!globals->haveBoxZ()){ | 
| 730 | 
< | 
        sprintf(painCave.errMsg, | 
| 731 | 
< | 
                "SimSetup error, no periodic BoxZ size given.\n"); | 
| 732 | 
< | 
        painCave.isFatal = 1; | 
| 733 | 
< | 
        simError(); | 
| 734 | 
< | 
      } | 
| 735 | 
< | 
      boxVector[2] = globals->getBoxZ(); | 
| 736 | 
< | 
 | 
| 737 | 
< | 
      info[i].setBox(boxVector); | 
| 738 | 
< | 
    } | 
| 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 | 
  | 
 | 
| 781 | 
  | 
 | 
| 782 | 
  | 
      if (!globals->haveECR()){ | 
| 783 | 
  | 
        sprintf(painCave.errMsg, | 
| 784 | 
< | 
                "SimSetup Warning: using default value of 1/2 the smallest " | 
| 785 | 
< | 
                "box length for the electrostaticCutoffRadius.\n" | 
| 817 | 
< | 
                "I hope you have a very fast processor!\n"); | 
| 784 | 
> | 
                "SimSetup Warning: using default value of 15.0 angstroms" | 
| 785 | 
> | 
                "box length for the electrostaticCutoffRadius.\n"); | 
| 786 | 
  | 
        painCave.isFatal = 0; | 
| 787 | 
  | 
        simError(); | 
| 788 | 
< | 
        double smallest; | 
| 821 | 
< | 
        smallest = info[i].boxL[0]; | 
| 822 | 
< | 
        if (info[i].boxL[1] <= smallest) | 
| 823 | 
< | 
          smallest = info[i].boxL[1]; | 
| 824 | 
< | 
        if (info[i].boxL[2] <= smallest) | 
| 825 | 
< | 
          smallest = info[i].boxL[2]; | 
| 826 | 
< | 
        theEcr = 0.5 * smallest; | 
| 788 | 
> | 
        theEcr = 15.0; | 
| 789 | 
  | 
      } | 
| 790 | 
  | 
      else{ | 
| 791 | 
  | 
        theEcr = globals->getECR(); | 
| 803 | 
  | 
        theEst = globals->getEST(); | 
| 804 | 
  | 
      } | 
| 805 | 
  | 
 | 
| 806 | 
< | 
      info[i].setEcr(theEcr, theEst); | 
| 806 | 
> | 
      info[i].setDefaultEcr(theEcr, theEst); | 
| 807 | 
  | 
 | 
| 808 | 
  | 
      if (!globals->haveDielectric()){ | 
| 809 | 
  | 
        sprintf(painCave.errMsg, | 
| 818 | 
  | 
      if (usesDipoles){ | 
| 819 | 
  | 
        if (!globals->haveECR()){ | 
| 820 | 
  | 
          sprintf(painCave.errMsg, | 
| 821 | 
< | 
                  "SimSetup Warning: using default value of 1/2 the smallest " | 
| 822 | 
< | 
                  "box length for the electrostaticCutoffRadius.\n" | 
| 823 | 
< | 
                  "I hope you have a very fast processor!\n"); | 
| 824 | 
< | 
          painCave.isFatal = 0; | 
| 825 | 
< | 
          simError(); | 
| 864 | 
< | 
          double smallest; | 
| 865 | 
< | 
          smallest = info[i].boxL[0]; | 
| 866 | 
< | 
          if (info[i].boxL[1] <= smallest) | 
| 867 | 
< | 
            smallest = info[i].boxL[1]; | 
| 868 | 
< | 
          if (info[i].boxL[2] <= smallest) | 
| 869 | 
< | 
            smallest = info[i].boxL[2]; | 
| 870 | 
< | 
          theEcr = 0.5 * smallest; | 
| 821 | 
> | 
                  "SimSetup Warning: using default value of 15.0 angstroms" | 
| 822 | 
> | 
                  "box length for the electrostaticCutoffRadius.\n"); | 
| 823 | 
> | 
          painCave.isFatal = 0; | 
| 824 | 
> | 
          simError(); | 
| 825 | 
> | 
          theEcr = 15.0; | 
| 826 | 
  | 
        } | 
| 827 | 
  | 
        else{ | 
| 828 | 
  | 
          theEcr = globals->getECR(); | 
| 829 | 
  | 
        } | 
| 830 | 
< | 
 | 
| 830 | 
> | 
         | 
| 831 | 
  | 
        if (!globals->haveEST()){ | 
| 832 | 
  | 
          sprintf(painCave.errMsg, | 
| 833 | 
  | 
                  "SimSetup Warning: using default value of 0.05 * the " | 
| 840 | 
  | 
        else{ | 
| 841 | 
  | 
          theEst = globals->getEST(); | 
| 842 | 
  | 
        } | 
| 843 | 
< | 
 | 
| 844 | 
< | 
        info[i].setEcr(theEcr, theEst); | 
| 843 | 
> | 
         | 
| 844 | 
> | 
        info[i].setDefaultEcr(theEcr, theEst); | 
| 845 | 
  | 
      } | 
| 846 | 
  | 
    } | 
| 847 | 
  | 
  } | 
| 893 | 
– | 
 | 
| 848 | 
  | 
#ifdef IS_MPI | 
| 849 | 
  | 
  strcpy(checkPointMsg, "post processing checks out"); | 
| 850 | 
  | 
  MPIcheckPoint(); | 
| 851 | 
  | 
#endif // is_mpi | 
| 852 | 
  | 
} | 
| 853 | 
< | 
 | 
| 853 | 
> | 
   | 
| 854 | 
  | 
void SimSetup::initSystemCoords(void){ | 
| 855 | 
  | 
  int i; | 
| 856 | 
  | 
 | 
| 878 | 
  | 
    delete fileInit; | 
| 879 | 
  | 
  } | 
| 880 | 
  | 
  else{ | 
| 881 | 
< | 
#ifdef IS_MPI  | 
| 928 | 
< | 
 | 
| 881 | 
> | 
     | 
| 882 | 
  | 
    // no init from bass | 
| 883 | 
< | 
 | 
| 883 | 
> | 
     | 
| 884 | 
  | 
    sprintf(painCave.errMsg, | 
| 885 | 
< | 
            "Cannot intialize a parallel simulation without an initial configuration file.\n"); | 
| 885 | 
> | 
            "Cannot intialize a simulation without an initial configuration file.\n"); | 
| 886 | 
  | 
    painCave.isFatal = 1;; | 
| 887 | 
  | 
    simError(); | 
| 888 | 
< | 
 | 
| 936 | 
< | 
#else | 
| 937 | 
< | 
 | 
| 938 | 
< | 
    initFromBass(); | 
| 939 | 
< | 
 | 
| 940 | 
< | 
 | 
| 941 | 
< | 
#endif | 
| 888 | 
> | 
     | 
| 889 | 
  | 
  } | 
| 890 | 
  | 
 | 
| 891 | 
  | 
#ifdef IS_MPI | 
| 1293 | 
  | 
  NVT<RealIntegrator>* myNVT = NULL; | 
| 1294 | 
  | 
  NPTi<NPT<RealIntegrator> >* myNPTi = NULL; | 
| 1295 | 
  | 
  NPTf<NPT<RealIntegrator> >* myNPTf = NULL; | 
| 1296 | 
+ | 
  NPTxyz<NPT<RealIntegrator> >* myNPTxyz = NULL; | 
| 1297 | 
  | 
   | 
| 1298 | 
  | 
  for (k = 0; k < nInfo; k++){ | 
| 1299 | 
  | 
    switch (ensembleCase){ | 
| 1397 | 
  | 
 | 
| 1398 | 
  | 
        if (globals->haveTauThermostat()) | 
| 1399 | 
  | 
          myNPTf->setTauThermostat(globals->getTauThermostat()); | 
| 1400 | 
+ | 
 | 
| 1401 | 
  | 
        else{ | 
| 1402 | 
  | 
          sprintf(painCave.errMsg, | 
| 1403 | 
  | 
                  "SimSetup error: If you use an NPT\n" | 
| 1408 | 
  | 
 | 
| 1409 | 
  | 
        if (globals->haveTauBarostat()) | 
| 1410 | 
  | 
          myNPTf->setTauBarostat(globals->getTauBarostat()); | 
| 1411 | 
+ | 
 | 
| 1412 | 
  | 
        else{ | 
| 1413 | 
  | 
          sprintf(painCave.errMsg, | 
| 1414 | 
  | 
                  "SimSetup error: If you use an NPT\n" | 
| 1418 | 
  | 
        } | 
| 1419 | 
  | 
 | 
| 1420 | 
  | 
        info->the_integrator = myNPTf; | 
| 1421 | 
+ | 
        break; | 
| 1422 | 
+ | 
 | 
| 1423 | 
+ | 
      case NPTxyz_ENS: | 
| 1424 | 
+ | 
        if (globals->haveZconstraints()){ | 
| 1425 | 
+ | 
          setupZConstraint(info[k]); | 
| 1426 | 
+ | 
          myNPTxyz = new ZConstraint<NPTxyz<NPT <RealIntegrator> > >(&(info[k]), the_ff); | 
| 1427 | 
+ | 
        } | 
| 1428 | 
+ | 
        else | 
| 1429 | 
+ | 
          myNPTxyz = new NPTxyz<NPT <RealIntegrator> >(&(info[k]), the_ff); | 
| 1430 | 
+ | 
 | 
| 1431 | 
+ | 
        myNPTxyz->setTargetTemp(globals->getTargetTemp()); | 
| 1432 | 
+ | 
 | 
| 1433 | 
+ | 
        if (globals->haveTargetPressure()) | 
| 1434 | 
+ | 
          myNPTxyz->setTargetPressure(globals->getTargetPressure()); | 
| 1435 | 
+ | 
        else{ | 
| 1436 | 
+ | 
          sprintf(painCave.errMsg, | 
| 1437 | 
+ | 
                  "SimSetup error: If you use a constant pressure\n" | 
| 1438 | 
+ | 
                  "    ensemble, you must set targetPressure in the BASS file.\n"); | 
| 1439 | 
+ | 
          painCave.isFatal = 1; | 
| 1440 | 
+ | 
          simError(); | 
| 1441 | 
+ | 
        }     | 
| 1442 | 
+ | 
 | 
| 1443 | 
+ | 
        if (globals->haveTauThermostat()) | 
| 1444 | 
+ | 
          myNPTxyz->setTauThermostat(globals->getTauThermostat()); | 
| 1445 | 
+ | 
        else{ | 
| 1446 | 
+ | 
          sprintf(painCave.errMsg, | 
| 1447 | 
+ | 
                  "SimSetup error: If you use an NPT\n" | 
| 1448 | 
+ | 
                  "    ensemble, you must set tauThermostat.\n"); | 
| 1449 | 
+ | 
          painCave.isFatal = 1; | 
| 1450 | 
+ | 
          simError(); | 
| 1451 | 
+ | 
        } | 
| 1452 | 
+ | 
 | 
| 1453 | 
+ | 
        if (globals->haveTauBarostat()) | 
| 1454 | 
+ | 
          myNPTxyz->setTauBarostat(globals->getTauBarostat()); | 
| 1455 | 
+ | 
        else{ | 
| 1456 | 
+ | 
          sprintf(painCave.errMsg, | 
| 1457 | 
+ | 
                  "SimSetup error: If you use an NPT\n" | 
| 1458 | 
+ | 
                  "    ensemble, you must set tauBarostat.\n"); | 
| 1459 | 
+ | 
          painCave.isFatal = 1; | 
| 1460 | 
+ | 
          simError(); | 
| 1461 | 
+ | 
        } | 
| 1462 | 
+ | 
 | 
| 1463 | 
+ | 
        info->the_integrator = myNPTxyz; | 
| 1464 | 
  | 
        break; | 
| 1465 | 
  | 
 | 
| 1466 | 
  | 
      default: |