| 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 |
|
|
| 94 |
|
#endif // is_mpi |
| 95 |
|
|
| 96 |
|
void SimSetup::createSim(void){ |
| 93 |
– |
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, "NPTxyz")){ |
| 612 |
+ |
ensembleCase = NPTxyz_ENS; |
| 613 |
+ |
} |
| 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(); |
| 699 |
< |
boxVector[1] = globals->getBox(); |
| 700 |
< |
boxVector[2] = globals->getBox(); |
| 701 |
< |
|
| 702 |
< |
info[i].setBox(boxVector); |
| 703 |
< |
} |
| 704 |
< |
else if (globals->haveDensity()){ |
| 705 |
< |
double vol; |
| 706 |
< |
vol = (double) tot_nmol / globals->getDensity(); |
| 707 |
< |
boxVector[0] = pow(vol, (1.0 / 3.0)); |
| 708 |
< |
boxVector[1] = boxVector[0]; |
| 709 |
< |
boxVector[2] = boxVector[0]; |
| 710 |
< |
|
| 711 |
< |
info[i].setBox(boxVector); |
| 712 |
< |
} |
| 713 |
< |
else{ |
| 714 |
< |
if (!globals->haveBoxX()){ |
| 715 |
< |
sprintf(painCave.errMsg, |
| 716 |
< |
"SimSetup error, no periodic BoxX size given.\n"); |
| 717 |
< |
painCave.isFatal = 1; |
| 718 |
< |
simError(); |
| 719 |
< |
} |
| 720 |
< |
boxVector[0] = globals->getBoxX(); |
| 721 |
< |
|
| 722 |
< |
if (!globals->haveBoxY()){ |
| 723 |
< |
sprintf(painCave.errMsg, |
| 724 |
< |
"SimSetup error, no periodic BoxY size given.\n"); |
| 725 |
< |
painCave.isFatal = 1; |
| 726 |
< |
simError(); |
| 727 |
< |
} |
| 728 |
< |
boxVector[1] = globals->getBoxY(); |
| 729 |
< |
|
| 730 |
< |
if (!globals->haveBoxZ()){ |
| 731 |
< |
sprintf(painCave.errMsg, |
| 732 |
< |
"SimSetup error, no periodic BoxZ size given.\n"); |
| 733 |
< |
painCave.isFatal = 1; |
| 734 |
< |
simError(); |
| 735 |
< |
} |
| 736 |
< |
boxVector[2] = globals->getBoxZ(); |
| 737 |
< |
|
| 738 |
< |
info[i].setBox(boxVector); |
| 739 |
< |
} |
| 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" |
| 818 |
< |
"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; |
| 822 |
< |
smallest = info[i].boxL[0]; |
| 823 |
< |
if (info[i].boxL[1] <= smallest) |
| 824 |
< |
smallest = info[i].boxL[1]; |
| 825 |
< |
if (info[i].boxL[2] <= smallest) |
| 826 |
< |
smallest = info[i].boxL[2]; |
| 827 |
< |
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(); |
| 865 |
< |
double smallest; |
| 866 |
< |
smallest = info[i].boxL[0]; |
| 867 |
< |
if (info[i].boxL[1] <= smallest) |
| 868 |
< |
smallest = info[i].boxL[1]; |
| 869 |
< |
if (info[i].boxL[2] <= smallest) |
| 870 |
< |
smallest = info[i].boxL[2]; |
| 871 |
< |
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 |
|
} |
| 894 |
– |
|
| 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 |
|
|
| 867 |
|
if (worldRank == 0){ |
| 868 |
|
#endif //is_mpi |
| 869 |
|
inName = globals->getInitialConfig(); |
| 917 |
– |
double* tempDouble = new double[1000000]; |
| 870 |
|
fileInit = new InitializeFromFile(inName); |
| 871 |
|
#ifdef IS_MPI |
| 872 |
|
} |
| 878 |
|
delete fileInit; |
| 879 |
|
} |
| 880 |
|
else{ |
| 881 |
< |
#ifdef IS_MPI |
| 930 |
< |
|
| 881 |
> |
|
| 882 |
|
// no init from bass |
| 883 |
< |
|
| 883 |
> |
|
| 884 |
|
sprintf(painCave.errMsg, |
| 885 |
< |
"Cannot intialize a parallel simulation without an initial configuration file.\n"); |
| 886 |
< |
painCave.isFatal; |
| 885 |
> |
"Cannot intialize a simulation without an initial configuration file.\n"); |
| 886 |
> |
painCave.isFatal = 1;; |
| 887 |
|
simError(); |
| 888 |
< |
|
| 938 |
< |
#else |
| 939 |
< |
|
| 940 |
< |
initFromBass(); |
| 941 |
< |
|
| 942 |
< |
|
| 943 |
< |
#endif |
| 888 |
> |
|
| 889 |
|
} |
| 890 |
|
|
| 891 |
|
#ifdef IS_MPI |
| 1103 |
|
} |
| 1104 |
|
|
| 1105 |
|
void SimSetup::calcSysValues(void){ |
| 1106 |
< |
int i, j, k; |
| 1106 |
> |
int i; |
| 1107 |
|
|
| 1108 |
|
int* molMembershipArray; |
| 1109 |
|
|
| 1202 |
|
|
| 1203 |
|
|
| 1204 |
|
void SimSetup::makeSysArrays(void){ |
| 1205 |
< |
int i, j, k, l; |
| 1205 |
> |
|
| 1206 |
> |
#ifndef IS_MPI |
| 1207 |
> |
int k, j; |
| 1208 |
> |
#endif // is_mpi |
| 1209 |
> |
int i, l; |
| 1210 |
|
|
| 1211 |
|
Atom** the_atoms; |
| 1212 |
|
Molecule* the_molecules; |
| 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" |
| 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: |
| 1467 |
|
sprintf(painCave.errMsg, |
| 1468 |
|
"SimSetup Error. Unrecognized ensemble in case statement.\n"); |