| 29 |
|
#define FF_EAM 2 |
| 30 |
|
|
| 31 |
|
using namespace std; |
| 32 |
+ |
|
| 33 |
+ |
/** |
| 34 |
+ |
* Check whether dividend is divisble by divisor or not |
| 35 |
+ |
*/ |
| 36 |
+ |
bool isDivisible(double dividend, double divisor){ |
| 37 |
+ |
double tolerance = 0.000001; |
| 38 |
+ |
double quotient; |
| 39 |
+ |
double diff; |
| 40 |
+ |
int intQuotient; |
| 41 |
+ |
|
| 42 |
+ |
quotient = dividend / divisor; |
| 43 |
+ |
|
| 44 |
+ |
if (quotient < 0) |
| 45 |
+ |
quotient = -quotient; |
| 46 |
+ |
|
| 47 |
+ |
intQuotient = int (quotient + tolerance); |
| 48 |
+ |
|
| 49 |
+ |
diff = fabs(fabs(dividend) - intQuotient * fabs(divisor)); |
| 50 |
+ |
|
| 51 |
+ |
if (diff <= tolerance) |
| 52 |
+ |
return true; |
| 53 |
+ |
else |
| 54 |
+ |
return false; |
| 55 |
+ |
} |
| 56 |
|
|
| 57 |
|
SimSetup::SimSetup(){ |
| 58 |
|
|
| 689 |
|
" Please give nMol in the components.\n"); |
| 690 |
|
painCave.isFatal = 1; |
| 691 |
|
simError(); |
| 692 |
+ |
} |
| 693 |
+ |
|
| 694 |
+ |
//check whether sample time, status time, thermal time and reset time are divisble by dt |
| 695 |
+ |
if (!isDivisible(globals->getSampleTime(), globals->getDt())){ |
| 696 |
+ |
sprintf(painCave.errMsg, |
| 697 |
+ |
"Sample time is not divisible by dt \n"); |
| 698 |
+ |
painCave.isFatal = 0; |
| 699 |
+ |
simError(); |
| 700 |
+ |
} |
| 701 |
+ |
|
| 702 |
+ |
if (globals->haveStatusTime() && !isDivisible(globals->getSampleTime(), globals->getDt())){ |
| 703 |
+ |
sprintf(painCave.errMsg, |
| 704 |
+ |
"Status time is not divisible by dt\n"); |
| 705 |
+ |
painCave.isFatal = 0; |
| 706 |
+ |
simError(); |
| 707 |
|
} |
| 708 |
|
|
| 709 |
+ |
if (globals->haveThermalTime() && !isDivisible(globals->getThermalTime(), globals->getDt())){ |
| 710 |
+ |
sprintf(painCave.errMsg, |
| 711 |
+ |
"Thermal time is not divisible by dt\n"); |
| 712 |
+ |
painCave.isFatal = 0; |
| 713 |
+ |
simError(); |
| 714 |
+ |
} |
| 715 |
+ |
|
| 716 |
+ |
if (globals->haveResetTime() && !isDivisible(globals->getResetTime(), globals->getDt())){ |
| 717 |
+ |
sprintf(painCave.errMsg, |
| 718 |
+ |
"Reset time is not divisible by dt\n"); |
| 719 |
+ |
painCave.isFatal = 0; |
| 720 |
+ |
simError(); |
| 721 |
+ |
} |
| 722 |
+ |
|
| 723 |
|
// set the status, sample, and thermal kick times |
| 724 |
|
|
| 725 |
|
for (i = 0; i < nInfo; i++){ |
| 749 |
|
} |
| 750 |
|
|
| 751 |
|
// check for the temperature set flag |
| 752 |
< |
|
| 752 |
> |
|
| 753 |
|
if (globals->haveTempSet()) |
| 754 |
|
info[i].setTemp = globals->getTempSet(); |
| 755 |
|
|
| 756 |
< |
// get some of the tricky things that may still be in the globals |
| 756 |
> |
// check for the extended State init |
| 757 |
|
|
| 758 |
< |
double boxVector[3]; |
| 759 |
< |
if (globals->haveBox()){ |
| 760 |
< |
boxVector[0] = globals->getBox(); |
| 708 |
< |
boxVector[1] = globals->getBox(); |
| 709 |
< |
boxVector[2] = globals->getBox(); |
| 710 |
< |
|
| 711 |
< |
info[i].setBox(boxVector); |
| 712 |
< |
} |
| 713 |
< |
else if (globals->haveDensity()){ |
| 714 |
< |
double vol; |
| 715 |
< |
vol = (double) tot_nmol / globals->getDensity(); |
| 716 |
< |
boxVector[0] = pow(vol, (1.0 / 3.0)); |
| 717 |
< |
boxVector[1] = boxVector[0]; |
| 718 |
< |
boxVector[2] = boxVector[0]; |
| 719 |
< |
|
| 720 |
< |
info[i].setBox(boxVector); |
| 721 |
< |
} |
| 722 |
< |
else{ |
| 723 |
< |
if (!globals->haveBoxX()){ |
| 724 |
< |
sprintf(painCave.errMsg, |
| 725 |
< |
"SimSetup error, no periodic BoxX size given.\n"); |
| 726 |
< |
painCave.isFatal = 1; |
| 727 |
< |
simError(); |
| 728 |
< |
} |
| 729 |
< |
boxVector[0] = globals->getBoxX(); |
| 730 |
< |
|
| 731 |
< |
if (!globals->haveBoxY()){ |
| 732 |
< |
sprintf(painCave.errMsg, |
| 733 |
< |
"SimSetup error, no periodic BoxY size given.\n"); |
| 734 |
< |
painCave.isFatal = 1; |
| 735 |
< |
simError(); |
| 736 |
< |
} |
| 737 |
< |
boxVector[1] = globals->getBoxY(); |
| 738 |
< |
|
| 739 |
< |
if (!globals->haveBoxZ()){ |
| 740 |
< |
sprintf(painCave.errMsg, |
| 741 |
< |
"SimSetup error, no periodic BoxZ size given.\n"); |
| 742 |
< |
painCave.isFatal = 1; |
| 743 |
< |
simError(); |
| 744 |
< |
} |
| 745 |
< |
boxVector[2] = globals->getBoxZ(); |
| 746 |
< |
|
| 747 |
< |
info[i].setBox(boxVector); |
| 748 |
< |
} |
| 758 |
> |
info[i].useInitXSstate = globals->getUseInitXSstate(); |
| 759 |
> |
info[i].orthoTolerance = globals->getOrthoBoxTolerance(); |
| 760 |
> |
|
| 761 |
|
} |
| 762 |
< |
|
| 762 |
> |
|
| 763 |
|
//setup seed for random number generator |
| 764 |
|
int seedValue; |
| 765 |
|
|
| 834 |
|
|
| 835 |
|
if (!globals->haveECR()){ |
| 836 |
|
sprintf(painCave.errMsg, |
| 837 |
< |
"SimSetup Warning: using default value of 1/2 the smallest " |
| 838 |
< |
"box length for the electrostaticCutoffRadius.\n" |
| 827 |
< |
"I hope you have a very fast processor!\n"); |
| 837 |
> |
"SimSetup Warning: using default value of 15.0 angstroms" |
| 838 |
> |
"box length for the electrostaticCutoffRadius.\n"); |
| 839 |
|
painCave.isFatal = 0; |
| 840 |
|
simError(); |
| 841 |
< |
double smallest; |
| 831 |
< |
smallest = info[i].boxL[0]; |
| 832 |
< |
if (info[i].boxL[1] <= smallest) |
| 833 |
< |
smallest = info[i].boxL[1]; |
| 834 |
< |
if (info[i].boxL[2] <= smallest) |
| 835 |
< |
smallest = info[i].boxL[2]; |
| 836 |
< |
theEcr = 0.5 * smallest; |
| 841 |
> |
theEcr = 15.0; |
| 842 |
|
} |
| 843 |
|
else{ |
| 844 |
|
theEcr = globals->getECR(); |
| 856 |
|
theEst = globals->getEST(); |
| 857 |
|
} |
| 858 |
|
|
| 859 |
< |
info[i].setEcr(theEcr, theEst); |
| 859 |
> |
info[i].setDefaultEcr(theEcr, theEst); |
| 860 |
|
|
| 861 |
|
if (!globals->haveDielectric()){ |
| 862 |
|
sprintf(painCave.errMsg, |
| 871 |
|
if (usesDipoles){ |
| 872 |
|
if (!globals->haveECR()){ |
| 873 |
|
sprintf(painCave.errMsg, |
| 874 |
< |
"SimSetup Warning: using default value of 1/2 the smallest " |
| 875 |
< |
"box length for the electrostaticCutoffRadius.\n" |
| 876 |
< |
"I hope you have a very fast processor!\n"); |
| 877 |
< |
painCave.isFatal = 0; |
| 878 |
< |
simError(); |
| 874 |
< |
double smallest; |
| 875 |
< |
smallest = info[i].boxL[0]; |
| 876 |
< |
if (info[i].boxL[1] <= smallest) |
| 877 |
< |
smallest = info[i].boxL[1]; |
| 878 |
< |
if (info[i].boxL[2] <= smallest) |
| 879 |
< |
smallest = info[i].boxL[2]; |
| 880 |
< |
theEcr = 0.5 * smallest; |
| 874 |
> |
"SimSetup Warning: using default value of 15.0 angstroms" |
| 875 |
> |
"box length for the electrostaticCutoffRadius.\n"); |
| 876 |
> |
painCave.isFatal = 0; |
| 877 |
> |
simError(); |
| 878 |
> |
theEcr = 15.0; |
| 879 |
|
} |
| 880 |
|
else{ |
| 881 |
|
theEcr = globals->getECR(); |
| 882 |
|
} |
| 883 |
< |
|
| 883 |
> |
|
| 884 |
|
if (!globals->haveEST()){ |
| 885 |
|
sprintf(painCave.errMsg, |
| 886 |
|
"SimSetup Warning: using default value of 0.05 * the " |
| 893 |
|
else{ |
| 894 |
|
theEst = globals->getEST(); |
| 895 |
|
} |
| 896 |
< |
|
| 897 |
< |
info[i].setEcr(theEcr, theEst); |
| 896 |
> |
|
| 897 |
> |
info[i].setDefaultEcr(theEcr, theEst); |
| 898 |
|
} |
| 899 |
|
} |
| 900 |
|
} |
| 903 |
– |
|
| 901 |
|
#ifdef IS_MPI |
| 902 |
|
strcpy(checkPointMsg, "post processing checks out"); |
| 903 |
|
MPIcheckPoint(); |
| 904 |
|
#endif // is_mpi |
| 905 |
|
} |
| 906 |
< |
|
| 906 |
> |
|
| 907 |
|
void SimSetup::initSystemCoords(void){ |
| 908 |
|
int i; |
| 909 |
|
|
| 931 |
|
delete fileInit; |
| 932 |
|
} |
| 933 |
|
else{ |
| 934 |
< |
#ifdef IS_MPI |
| 938 |
< |
|
| 934 |
> |
|
| 935 |
|
// no init from bass |
| 936 |
< |
|
| 936 |
> |
|
| 937 |
|
sprintf(painCave.errMsg, |
| 938 |
< |
"Cannot intialize a parallel simulation without an initial configuration file.\n"); |
| 938 |
> |
"Cannot intialize a simulation without an initial configuration file.\n"); |
| 939 |
|
painCave.isFatal = 1;; |
| 940 |
|
simError(); |
| 941 |
< |
|
| 946 |
< |
#else |
| 947 |
< |
|
| 948 |
< |
initFromBass(); |
| 949 |
< |
|
| 950 |
< |
|
| 951 |
< |
#endif |
| 941 |
> |
|
| 942 |
|
} |
| 943 |
|
|
| 944 |
|
#ifdef IS_MPI |
| 1450 |
|
|
| 1451 |
|
if (globals->haveTauThermostat()) |
| 1452 |
|
myNPTf->setTauThermostat(globals->getTauThermostat()); |
| 1453 |
+ |
|
| 1454 |
|
else{ |
| 1455 |
|
sprintf(painCave.errMsg, |
| 1456 |
|
"SimSetup error: If you use an NPT\n" |
| 1461 |
|
|
| 1462 |
|
if (globals->haveTauBarostat()) |
| 1463 |
|
myNPTf->setTauBarostat(globals->getTauBarostat()); |
| 1464 |
+ |
|
| 1465 |
|
else{ |
| 1466 |
|
sprintf(painCave.errMsg, |
| 1467 |
|
"SimSetup error: If you use an NPT\n" |