45 #if defined(USE_GETTIMEOFDAY) 50 gettimeofday(&t, NULL);
51 double sec = double(t.tv_sec);
52 double usec = double (t.tv_usec);
53 double totalSec = sec + 1e-6 * usec;
58 return CmiWallTimer();
62 #if(CMK_CCS_AVAILABLE && CMK_WEB_MODE) 63 extern "C" void CApplicationDepositNode0Data(
char *);
68 #define cbrt(x) pow(x,(double)(1.0/3.0)) 72 #define MIN_DEBUG_LEVEL 3 76 #define XXXBIGREAL 1.0e32 81 static char tmp_string[21];
82 sprintf(tmp_string,
"FepEnergy: %6d ",
X);
88 static char tmp_string[21];
89 sprintf(tmp_string,
"FepE_back: %6d ",
X);
95 static char tmp_string[21];
96 sprintf(tmp_string,
"FEP: %7d",
X);
102 static char tmp_string[21];
103 sprintf(tmp_string,
"TI: %7d",
X);
121 const char *myname,
int outputfreq)
122 : nslabs(numslabs), freq(outputfreq) {
123 name = strdup(myname);
124 nelements = 3*nslabs * numpartitions;
141 double inv_volume = 1.0 / lattice.
volume();
143 int arraysize = 3*nslabs;
144 for (i=0; i<nelements; i++) {
147 total[i % arraysize] += val;
150 if (!(step % freq)) {
154 iout <<
"PPROFILE" << name <<
": " << step <<
" ";
155 if (step == firsttimestep) {
157 for (i=0; i<nelements; i++)
158 iout << reduction->
item(i)*scalefac*inv_volume <<
" ";
162 for (i=0; i<nelements; i++)
178 computeChecksum(0), marginViolations(0), pairlistWarnings(0),
183 firstCTime(CmiTimer()),
215 int npairs = (ntypes * (ntypes+1))/2;
221 nslabs, npairs,
"BONDED", freq);
223 nslabs, npairs,
"NONBONDED", freq);
226 nslabs, ntypes,
"INTERNAL", freq);
230 nslabs, npairs,
"NONBONDED", freq);
252 #define AVGXY(T) T.xy = T.yx = 0; T.xx = T.yy = 0.5 * ( T.xx + T.yy );\ 253 T.xz = T.zx = T.yz = T.zy = 0.5 * ( T.xz + T.yz ); 283 for (
int i=0;i < n;i++) {
344 DebugM(4,
"Starting thread in controller on this=" <<
this <<
"\n");
346 threadWrapper->
thread = CthCreate((CthVoidFn)&(threadRun),(
void*)(
this),
CTRL_STK_SZ);
347 CthSetStrategyDefault(threadWrapper->
thread);
348 #if CMK_BLUEGENE_CHARM 355 CthAwaken(threadWrapper->
thread);
368 #ifdef NODEGROUP_FORCE_REGISTER 369 CProxy_PatchData cpdata(CkpvAccess(BOCclass_group).patchData);
370 PatchData *patchData = cpdata.ckLocalBranch();
372 switch ( scriptTask ) {
406 NAMD_die(
"Unable to revert, checkpoint was never called!");
434 NAMD_die(
"Unable to load checkpoint, requested key was never stored.");
451 NAMD_die(
"Unable to swap checkpoint, requested key was never stored.");
468 NAMD_die(
"Unable to free checkpoint, requested key was never stored.");
489 NAMD_bug(
"Unknown task in Controller::algorithm");
541 if ( zeroMomentum && dofull && ! (step %
slowFreq) )
550 sprintf(traceNote,
"s:%d", step);
551 traceUserSuppliedNote(traceNote);
563 for ( ++step ; step <= numberOfSteps; ++step )
574 if ( zeroMomentum && dofull && ! (step %
slowFreq) )
589 sprintf(traceNote,
"s:%d", step);
590 traceUserSuppliedNote(traceNote);
616 #ifdef MEASURE_NAMD_WITH_PAPI 618 int bstep =
simParams->papiMeasureStartStep;
619 int estep = bstep +
simParams->numPapiMeasureSteps;
621 papiMeasureBarrier(1, step);
624 papiMeasureBarrier(0, step);
647 #ifdef NODEGROUP_FORCE_REGISTER 670 if ( zeroMomentum && dofull && ! (step %
slowFreq) )
679 sprintf(traceNote,
"s:%d", step);
680 traceUserSuppliedNote(traceNote);
699 if ( zeroMomentum && dofull && ! (step %
slowFreq) )
714 sprintf(traceNote,
"s:%d", step);
715 traceUserSuppliedNote(traceNote);
725 printMinimizeEnergies(step); \ 726 outputExtendedSystem(step); \ 727 rebalanceLoad(step); \ 728 if ( step == numberOfSteps ) return; \ 732 if ( step == numberOfSteps ) { \ 733 if ( minVerbose ) { iout << "LINE MINIMIZER: RETURNING TO " << mid.x << " FROM " << last.x << "\n" << endi; } \ 734 if ( newDir || (mid.x-last.x) ) { \ 735 broadcast->minimizeCoefficient.publish(minSeq++,mid.x-last.x); \ 737 broadcast->minimizeCoefficient.publish(minSeq++,0.); \ 738 broadcast->minimizeCoefficient.publish(minSeq++,0.); \ 739 min_reduction->require(); \ 740 broadcast->minimizeCoefficient.publish(minSeq++,0.); \ 742 enqueueCollections(step); \ 744 } else if ( (X)-last.x ) { \ 745 broadcast->minimizeCoefficient.publish(minSeq++,(X)-last.x); \ 748 enqueueCollections(step); \ 750 last.u = min_energy; \ 751 last.dudx = -1. * min_f_dot_v; \ 752 last.noGradient = min_huge_count; \ 753 if ( minVerbose ) { \ 754 iout << "LINE MINIMIZER: POSITION " << last.x << " ENERGY " << last.u << " GRADIENT " << last.dudx; \ 755 if ( last.noGradient ) iout << " HUGECOUNT " << last.noGradient; \ 756 iout << "\n" << endi; \ 773 const BigReal goldenRatio = 0.5 * ( sqrt(5.0) - 1.0 );
780 int old_min_huge_count = 2000000000;
781 int steps_at_same_min_huge_count = 0;
782 for (
int i_slow = 0;
min_huge_count && i_slow < 100; ++i_slow ) {
784 if ( ++steps_at_same_min_huge_count > 10 )
break;
787 steps_at_same_min_huge_count = 0;
797 iout <<
"MINIMIZER STARTING CONJUGATE GRADIENT ALGORITHM\n" <<
endi;
800 int errorFactor = 10;
819 iout <<
"LINE MINIMIZER: POSITION " << last.
x <<
" ENERGY " << last.
u <<
" GRADIENT " << last.
dudx;
824 iout <<
"LINE MINIMIZER REDUCING GRADIENT FROM " <<
831 while ( last.
u < mid.
u ) {
832 if ( last.
dudx < mid.
dudx * (5.5 - x/maxstep)/5. ) {
833 x = 6 * maxstep;
break;
835 lo = mid; mid = last;
837 if ( x > 5.5 * maxstep )
break;
840 if ( x > 5.5 * maxstep ) {
841 iout <<
"MINIMIZER RESTARTING CONJUGATE GRADIENT ALGORITHM DUE TO POOR PROGRESS\n" <<
endi;
851 #define PRINT_BRACKET \ 852 iout << "LINE MINIMIZER BRACKET: DX " \ 853 << (mid.x-lo.x) << " " << (hi.x-mid.x) << \ 854 " DU " << (mid.u-lo.u) << " " << (hi.u-mid.u) << " DUDX " << \ 855 lo.dudx << " " << mid.dudx << " " << hi.dudx << " \n" << endi; 859 for ( itcnt = 10; itcnt > 0; --itcnt ) {
863 if ( ( mid.
x - lo.
x ) > ( hi.
x - mid.
x ) ) {
864 x = (1.0 - goldenRatio) * lo.
x + goldenRatio * mid.
x;
866 if ( last.
u <= mid.
u ) {
867 hi = mid; mid = last;
872 x = (1.0 - goldenRatio) * hi.
x + goldenRatio * mid.
x;
874 if ( last.
u <= mid.
u ) {
875 lo = mid; mid = last;
881 if ( mid.
dudx > 0. ) {
882 BigReal altxhi = 0.1 * lo.
x + 0.9 * mid.
x;
883 BigReal altxlo = 0.9 * lo.
x + 0.1 * mid.
x;
884 x = mid.
dudx*(mid.
x*mid.
x-lo.
x*lo.
x) + 2*mid.
x*(lo.
u-mid.
u);
886 if ( xdiv ) x /= xdiv;
else x = last.
x;
887 if ( x > altxhi ) x = altxhi;
888 if ( x < altxlo ) x = altxlo;
889 if ( x-last.
x == 0 )
break;
891 if ( last.
u <= mid.
u ) {
892 hi = mid; mid = last;
894 if ( lo.
dudx < 0. && last.
dudx > 0. ) progress = 0;
898 BigReal altxlo = 0.1 * hi.
x + 0.9 * mid.
x;
899 BigReal altxhi = 0.9 * hi.
x + 0.1 * mid.
x;
900 x = mid.
dudx*(mid.
x*mid.
x-hi.
x*hi.
x) + 2*mid.
x*(hi.
u-mid.
u);
902 if ( xdiv ) x /= xdiv;
else x = last.
x;
903 if ( x < altxlo ) x = altxlo;
904 if ( x > altxhi ) x = altxhi;
905 if ( x-last.
x == 0 )
break;
907 if ( last.
u <= mid.
u ) {
908 lo = mid; mid = last;
910 if ( hi.
dudx > 0. && last.
dudx < 0. ) progress = 0;
920 if ( fabs(last.
dudx) < tol )
break;
921 if ( lo.
x != mid.
x && (lo.
u-mid.
u) < tol * (mid.
x-lo.
x) )
break;
922 if ( mid.
x != hi.
x && (hi.
u-mid.
u) < tol * (hi.
x-mid.
x) )
break;
927 c = ( c > 1.5 ? 1.5 : c );
928 if ( atStart ) { c = 0; --atStart; }
931 if ( errorFactor < 100 ) errorFactor += 10;
934 iout <<
"MINIMIZER RESTARTING CONJUGATE GRADIENT ALGORITHM\n" <<
endi;
954 BigReal NGfac = 1.0/sqrt(3.0*NG + 1.0);
971 if (callNumber == 1) {
1002 calcPressure(step, 0, virial_normal, virial_nbond, virial_slow,
1003 intVirial_normal, intVirial_nbond, intVirial_slow,
1004 extForce_normal, extForce_nbond, extForce_slow);
1005 if (callNumber == 2) {
1035 for (
int i=n-1;i >= 0;i--) {
1039 }
else if (i == n-1) {
1049 if (callNumber == 1)
1071 for (
int i=0;i < n;i++) {
1075 }
else if (i == n-1) {
1109 BigReal W = (3.0*NG + 1.0)*kT0*tauV*tauV;
1144 BigReal weight = (totalEnergyNew - totalEnergyOld);
1145 weight += (pressureTarget * (volumeNew - volumeOld));
1146 weight -= (numAtom * kbT * log(volumeNew / volumeOld));
1147 weight -= (surfTensionTarget * (areaNew - areaOld));
1148 BigReal totalWeight = exp(-weight / kbT);
1152 printf(
"MC-data (accept): step: %d, Pe: %d, numAtom: %d\n", step, CkMyPe(), numAtom);
1153 printf(
" oldE: %f, newE: %f, deltaE: %f\n", totalEnergyOld, totalEnergyNew, totalEnergyNew - totalEnergyOld);
1154 printf(
" oldV: %f, newV: %f, deltaV: %f\n", volumeOld, volumeNew, volumeNew - volumeOld);
1155 printf(
" surfTen: %f, deltaA: %f \n", surfTensionTarget, areaNew - areaOld);
1156 printf(
" weight: %f, accepted: %d\n\n", totalWeight, accepted);
1172 if (accRate > 0.0) {
1184 iout <<
iWARN <<
"Small volume change in MC barostat!\n" <<
endi;
1185 iout <<
iWARN <<
"Maximum volume change is set to 0.001 A^3.\n" <<
endi;
1188 iout <<
iWARN <<
"Large volume change in MC barostat!\n" <<
endi;
1189 iout <<
iWARN <<
"Maximum volume change is set to " << maxVolume <<
" A^3.\n" <<
endi;
1197 iout <<
iWARN <<
"Small volume change in MC barostat!\n" <<
endi;
1198 iout <<
iWARN <<
"Maximum volume change is set to 0.001 A^3.\n" <<
endi;
1201 iout <<
iWARN <<
"Large volume change in MC barostat!\n" <<
endi;
1202 iout <<
iWARN <<
"Maximum volume change is set to " << maxVolume <<
" A^3.\n" <<
endi;
1210 iout <<
iWARN <<
"Small volume change in MC barostat!\n" <<
endi;
1211 iout <<
iWARN <<
"Maximum volume change is set to 0.001 A^3.\n" <<
endi;
1214 iout <<
iWARN <<
"Large volume change in MC barostat!\n" <<
endi;
1215 iout <<
iWARN <<
"Maximum volume change is set to " << maxVolume <<
" A^3.\n" <<
endi;
1227 " ACCEPTED: " <<
mc_totalAccept <<
" %ACCEPTANCE: " << acceptPerc;
1246 #ifdef NODEGROUP_FORCE_REGISTER 1247 publishedMCAcceptance.insert({step, accepted});
1270 bool constZ, constY, constX;
1276 constZ = constY = constX =
false;
1285 factor.
zz = (oldVolume + dV) / oldVolume;
1286 }
else if (constRatio) {
1290 if (probAxis < 1.0) {
1294 BigReal scale = sqrt((oldVolume + dV)/oldVolume);
1295 factor.
xx = factor.
yy = scale;
1297 }
else if (probAxis < 2.0 && !constZ) {
1301 BigReal scale = (oldVolume + dV)/oldVolume;
1310 if (probAxis < 1.0 && !constX) {
1314 factor.
xx = (oldVolume + dV) / oldVolume;
1316 }
else if (probAxis < 2.0 && !constY) {
1320 factor.
yy = (oldVolume + dV) / oldVolume;
1322 }
else if (probAxis < 3.0 && !constZ) {
1326 factor.
zz = (oldVolume + dV) / oldVolume;
1335 BigReal scale = pow((oldVolume + dV)/oldVolume, 1.0/3.0);
1336 factor.
xx = factor.
yy = factor.
zz = scale;
1345 #ifdef NODEGROUP_FORCE_REGISTER 1346 publishedRescaleFactors.insert({step, factor});
1353 #define LIMIT_SCALING(VAR,MIN,MAX,FLAG) {\ 1354 if ( VAR < (MIN) ) { VAR = (MIN); FLAG = 1; } \ 1355 if ( VAR > (MAX) ) { VAR = (MAX); FLAG = 1; } } 1379 " cell rescaling factor limited.\n" <<
endi;
1383 #ifdef NODEGROUP_FORCE_REGISTER 1384 publishedRescaleFactors.insert({step, factor});
1393 #undef LIMIT_SCALING 1407 #ifdef DEBUG_PRESSURE 1408 iout <<
iINFO <<
"entering langevinPiston1, strain rate: " << strainRate <<
"\n";
1414 BigReal f1 = exp( -0.5 * dt_long * gamma );
1415 BigReal f2 = sqrt( ( 1. - f1*f1 ) * kT / mass );
1431 #ifdef DEBUG_PRESSURE 1432 iout <<
iINFO <<
"applying langevin, strain rate: " << strainRate <<
"\n";
1446 strainRate += ( 0.5 * dt * cellDims *
state->lattice.
volume() / mass ) *
1456 #ifdef DEBUG_PRESSURE 1457 iout <<
iINFO <<
"integrating half step, strain rate: " << strainRate <<
"\n";
1466 factor.
xx = exp( dt_long * strainRate.xx );
1467 factor.
yy = exp( dt_long * strainRate.yy );
1469 factor.
xx = factor.
yy = 1;
1471 factor.
zz = exp( dt_long * strainRate.zz );
1472 #ifdef NODEGROUP_FORCE_REGISTER 1474 publishedRescaleFactors.insert({step, factor});
1481 #ifdef DEBUG_PRESSURE 1482 iout <<
iINFO <<
"rescaling by: " << factor <<
"\n";
1492 factor.
xx = exp( dt_long * strainRate.xx );
1493 factor.
yy = exp( dt_long * strainRate.yy );
1495 factor.
xx = factor.
yy = 1;
1497 factor.
zz = exp( dt_long * strainRate.zz );
1498 #ifdef NODEGROUP_FORCE_REGISTER 1500 publishedRescaleFactors.insert({step, factor});
1510 #ifdef DEBUG_PRESSURE 1522 factor.
xx = factor.
yy = 1;
1525 #ifdef NODEGROUP_FORCE_REGISTER 1527 publishedRescaleFactors.insert({step+1, factor});
1541 #ifdef DEBUG_PRESSURE 1542 iout <<
iINFO <<
"correcting strain rate for nbond, ";
1544 strainRate -= ( 0.5 * dt * cellDims *
state->lattice.
volume() / mass ) *
1546 #ifdef DEBUG_PRESSURE 1547 iout <<
"strain rate: " << strainRate <<
"\n";
1552 #ifdef DEBUG_PRESSURE 1553 iout <<
iINFO <<
"correcting strain rate for slow, ";
1555 strainRate -= ( 0.5 * dt * cellDims *
state->lattice.
volume() / mass ) *
1557 #ifdef DEBUG_PRESSURE 1558 iout <<
"strain rate: " << strainRate <<
"\n";
1586 #ifdef DEBUG_PRESSURE 1587 iout <<
iINFO <<
"correcting strain rate for nbond, ";
1589 strainRate += ( 0.5 * dt * cellDims *
state->lattice.
volume() / mass ) *
1591 #ifdef DEBUG_PRESSURE 1592 iout <<
"strain rate: " << strainRate <<
"\n";
1597 #ifdef DEBUG_PRESSURE 1598 iout <<
iINFO <<
"correcting strain rate for slow, ";
1600 strainRate += ( 0.5 * dt * cellDims *
state->lattice.
volume() / mass ) *
1602 #ifdef DEBUG_PRESSURE 1603 iout <<
"strain rate: " << strainRate <<
"\n";
1617 strainRate += ( 0.5 * dt * cellDims *
state->lattice.
volume() / mass ) *
1621 #ifdef DEBUG_PRESSURE 1622 iout <<
iINFO <<
"integrating half step, strain rate: " << strainRate <<
"\n";
1628 BigReal f1 = exp( -0.5 * dt_long * gamma );
1629 BigReal f2 = sqrt( ( 1. - f1*f1 ) * kT / mass );
1644 #ifdef DEBUG_PRESSURE 1645 iout <<
iINFO <<
"applying langevin, strain rate: " << strainRate <<
"\n";
1649 #ifdef DEBUG_PRESSURE 1650 iout <<
iINFO <<
"exiting langevinPiston2, strain rate: " << strainRate <<
"\n";
1665 if ( rescaleFreq > 0 ) {
1674 BigReal factor = sqrt(rescaleTemp/avgTemp);
1688 #ifdef NODEGROUP_FORCE_REGISTER 1700 #ifdef NODEGROUP_FORCE_REGISTER 1705 if ( momentum.
length2() == 0. )
1706 iout <<
iERROR <<
"Momentum is exactly zero; probably a bug.\n" <<
endi;
1708 NAMD_die(
"Total mass is zero in Controller::correctMomentum");
1710 momentum *= (-1./mass);
1725 iout <<
"FEP: RESETTING FOR NEW FEP WINDOW " 1726 <<
"LAMBDA SET TO " << alchLambda <<
" LAMBDA2 " << alchLambda2;
1727 if ( alchLambdaIDWS >= 0. ) {
1728 iout <<
" LAMBDA_IDWS " << alchLambdaIDWS;
1730 iout <<
"\nFEP: WINDOW TO HAVE " << alchEquilSteps
1731 <<
" STEPS OF EQUILIBRATION PRIOR TO FEP DATA COLLECTION.\n" 1732 <<
"FEP: USING CONSTANT TEMPERATURE OF " << alchTemp
1733 <<
" K FOR FEP CALCULATION\n" <<
endi;
1742 iout <<
"TI: RESETTING FOR NEW WINDOW " 1743 <<
"LAMBDA SET TO " << alchLambda
1744 <<
"\nTI: WINDOW TO HAVE " << alchEquilSteps
1745 <<
" STEPS OF EQUILIBRATION PRIOR TO TI DATA COLLECTION.\n" <<
endi;
1753 if ( ( reassignFreq > 0 ) && ! ( step % reassignFreq ) ) {
1760 if ( newTemp < simParams->reassignHold )
1763 iout <<
"REASSIGNING VELOCITIES AT STEP " << step
1764 <<
" TO " << newTemp <<
" KELVIN.\n" <<
endi;
1808 double coefficient = 1;
1832 static char tmp_string[50];
1833 static char format_string[50];
1834 const double maxnum = 99999999999.9999;
1835 if (
X > maxnum )
X = maxnum;
1836 if (
X < -maxnum )
X = -maxnum;
1838 int whole = (decimal <= 4 ? 14 : 10 + decimal);
1839 sprintf(format_string,
" %%%d.%df", whole, decimal);
1840 sprintf(tmp_string, format_string,
X);
1844 static char *
FORMAT(
const char *
X,
int decimal = 4)
1846 static char tmp_string[50];
1847 static char format_string[50];
1849 int width = (decimal <= 4 ? 14 : 10 + decimal);
1850 sprintf(format_string,
" %%%ds", width);
1851 sprintf(tmp_string, format_string,
X);
1857 static char tmp_string[21];
1858 sprintf(tmp_string,
"ENERGY: %7d",
X);
1899 int numFixedGroups = ( numFixedAtoms ? molecule->
numFixedGroups : 0 );
1902 if ( numFixedGroups ) numGroupDegFreedom -= 3 * numFixedGroups;
1906 numGroupDegFreedom -= 3;
1913 int numFixedRigidBonds =
1918 numDegFreedom -= ( numRigidBonds - numFixedRigidBonds - numLonepairs);
1920 BigReal groupKineticEnergyHalfstep, groupKineticEnergyCentered;
1921 #ifdef NODEGROUP_FORCE_REGISTER 1937 #ifdef NODEGROUP_FORCE_REGISTER 1947 BigReal groupTempHalfstep = 2.0 * groupKineticEnergyHalfstep
1949 BigReal groupTempCentered = 2.0 * groupKineticEnergyCentered
1961 #ifdef NODEGROUP_FORCE_REGISTER 1975 #ifdef NODEGROUP_FORCE_REGISTER 1980 #ifdef NODEGROUP_FORCE_REGISTER 1990 #ifdef NODEGROUP_FORCE_REGISTER 1995 #ifdef NODEGROUP_FORCE_REGISTER 2013 #ifdef NODEGROUP_FORCE_REGISTER 2027 BigReal drudeComKE, drudeBondKE;
2028 #ifdef NODEGROUP_FORCE_REGISTER 2040 #ifdef NODEGROUP_FORCE_REGISTER 2070 virial_normal, virial_nbond, virial_slow,
2071 intVirial_normal, intVirial_nbond, intVirial_slow,
2072 extForce_normal, extForce_nbond, extForce_slow);
2074 #ifdef DEBUG_PRESSURE 2103 const Tensor& virial_normal_in,
const Tensor& virial_nbond_in,
const Tensor& virial_slow_in,
2104 const Tensor& intVirial_normal,
const Tensor& intVirial_nbond,
const Tensor& intVirial_slow,
2105 const Vector& extForce_normal,
const Vector& extForce_nbond,
const Vector& extForce_slow) {
2107 Tensor virial_normal = virial_normal_in;
2108 Tensor virial_nbond = virial_nbond_in;
2109 Tensor virial_slow = virial_slow_in;
2112 Tensor tmp = virial_normal;
2113 fprintf(stderr,
"VIRIAL_NORMAL %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf\n",
2117 fprintf(stderr,
"VIRIAL_NBOND %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf\n",
2121 fprintf(stderr,
"VIRIAL_SLOW %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf\n",
2124 tmp = intVirial_normal;
2125 fprintf(stderr,
"INT_VIRIAL_NORMAL %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf\n",
2128 tmp = intVirial_nbond;
2129 fprintf(stderr,
"INT_VIRIAL_NBOND %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf\n",
2132 tmp = intVirial_slow;
2133 fprintf(stderr,
"INT_VIRIAL_SLOW %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf\n",
2136 Vector ext = extForce_normal;
2137 fprintf(stderr,
"EXT_FORCE_NORMAL %1.2lf %1.2lf %1.2lf\n",
2138 ext.
x, ext.
y, ext.
z);
2140 ext = extForce_nbond;
2141 fprintf(stderr,
"EXT_FORCE_NBOND %1.2lf %1.2lf %1.2lf\n",
2142 ext.
x, ext.
y, ext.
z);
2144 ext = extForce_slow;
2145 fprintf(stderr,
"EXT_FORCE_SLOW %1.2lf %1.2lf %1.2lf\n",
2146 ext.
x, ext.
y, ext.
z);
2149 fprintf(stderr,
"LATTICE %lf %lf %lf %lf %lf %lf %lf %lf %lf\n",
2154 fprintf(stderr,
"VOLUME %3.10lf\n",
state->lattice.
volume());
2160 Molecule *molecule = node->molecule;
2165 virial_normal -=
outer(extForce_normal,extPosition);
2166 virial_nbond -=
outer(extForce_nbond,extPosition);
2167 virial_slow -=
outer(extForce_slow,extPosition);
2169 if ( (volume=lattice.
volume()) != 0. )
2173 #ifdef MEM_OPT_VERSION 2174 NAMD_bug(
"LJcorrection not supported in memory optimized build.");
2269 #define AVGXY(T) T.xy = T.yx = 0; T.xx = T.yy = 0.5 * ( T.xx + T.yy );\ 2270 T.xz = T.zx = T.yz = T.zy = 0.5 * ( T.xz + T.yz ); 2297 else if(testV < *Vmin){
2302 delta = testV - *Vavg;
2304 *M2 += delta * (testV - (*Vavg));
2306 *sigmaV = sqrt(*M2/n);
2314 *k0 = (1-sigma0/sigmaV) * (Vmax-Vmin) / (Vavg-Vmin);
2322 *E = Vmin + (Vmax-Vmin)/(*k0);
2328 *k0 = (sigma0 / sigmaV) * (Vmax-Vmin) / (Vmax-Vavg);
2335 *k = *k0 / (Vmax - Vmin);
2345 *vir += vir_orig * (*factor);
2346 *dV += (*factor) * VE/2.;
2355 (
int step_n,
char type,
int V_n,
BigReal Vmax,
BigReal Vmin,
BigReal Vavg,
BigReal sigmaV,
BigReal M2,
BigReal E,
BigReal k,
bool write_topic,
bool lasttime){
2357 char timestepstr[20];
2361 const char *bsuffix;
2363 if(lasttime ||
simParams->restartFrequency == 0){
2372 baselen = strlen(filename);
2373 restart_name =
new char[baselen+26];
2375 strcpy(restart_name, filename);
2376 if (
simParams->restartSave && !lasttime) {
2377 sprintf(timestepstr,
".%d",step_n);
2378 strcat(restart_name, timestepstr);
2381 strcat(restart_name,
".gamd");
2386 rest = fopen(restart_name,
"w");
2388 fprintf(rest,
"# NAMD accelMDG restart file\n" 2389 "# D/T Vn Vmax Vmin Vavg sigmaV M2 E k\n" 2390 "%c %d %.17lg %.17lg %.17lg %.17lg %.17le %.17lg %.17lg\n",
2391 type, V_n-1, Vmax, Vmin, Vavg, sigmaV, M2, E, k);
2393 iout <<
"GAUSSIAN ACCELERATED MD RESTART DATA WRITTEN TO " << restart_name <<
"\n" <<
endi;
2396 iout <<
iWARN <<
"Cannot write accelMDG restart file\n" <<
endi;
2399 rest = fopen(restart_name,
"a");
2401 fprintf(rest,
"%c %d %.17lg %.17lg %.17lg %.17lg %.17le %.17lg %.17lg\n",
2402 type, V_n-1, Vmax, Vmin, Vavg, sigmaV, M2, E, k);
2406 iout <<
iWARN <<
"Cannot write accelMDG restart file\n" <<
endi;
2409 delete[] restart_name;
2443 static BigReal VmaxP, VminP, VavgP, M2P=0, sigmaVP;
2444 static BigReal VmaxD, VminD, VavgD, M2D=0, sigmaVD;
2447 static int V_n = 1, iEusedD, iEusedP;
2448 static char warnD, warnP;
2449 static int restfreq;
2468 BigReal amd_ljEnergy_Corr = 0.;
2485 volume = lattice.
volume();
2529 potentialEnergy = bondEnergy + angleEnergy + dihedralEnergy +
2530 improperEnergy + amd_electEnergy + amd_electEnergySlow + amd_ljEnergy +
2531 crosstermEnergy + boundaryEnergy + miscEnergy +
2532 amd_goTotalEnergy + amd_groLJEnergy + amd_groGaussEnergy;
2537 if(step == firststep) {
2539 warnD = warnP =
'\0';
2545 int dihe_n=0, tot_n=0;
2551 while(fgets(line, 256, rest)){
2553 dihe_n = sscanf(line+1,
" %d %la %la %la %la %la %la %la",
2554 &V_n, &VmaxD, &VminD, &VavgD, &sigmaVD, &M2D, &ED, &kD);
2556 else if(line[0] ==
'T'){
2557 tot_n = sscanf(line+1,
" %d %la %la %la %la %la %la %la",
2558 &V_n, &VmaxP, &VminP, &VavgP, &sigmaVP, &M2P, &EP, &kP);
2569 NAMD_die(
"Invalid dihedral potential energy format in accelMDG restart file");
2570 k0D = kD * (VmaxD - VminD);
2571 iout <<
"GAUSSIAN ACCELERATED MD READ FROM RESTART FILE: DIHED" 2572 <<
" Vmax " << VmaxD <<
" Vmin " << VminD
2573 <<
" Vavg " << VavgD <<
" sigmaV " << sigmaVD
2574 <<
" E " << ED <<
" k " << kD
2581 NAMD_die(
"Invalid total potential energy format in accelMDG restart file");
2582 k0P = kP * (VmaxP - VminP);
2583 iout <<
"GAUSSIAN ACCELERATED MD READ FROM RESTART FILE: TOTAL" 2584 <<
" Vmax " << VmaxP <<
" Vmin " << VminP
2585 <<
" Vavg " << VavgP <<
" sigmaV " << sigmaVP
2586 <<
" E " << EP <<
" k " << kP
2590 iEusedD = (ED == VmaxD) ? 1 : 2;
2591 iEusedP = (EP == VmaxP) ? 1 : 2;
2599 testV = dihedralEnergy + crosstermEnergy;
2602 if(step > firststep && step % restfreq == 0)
2603 write_accelMDG_rest_file(step,
'D', V_n, VmaxD, VminD, VavgD, sigmaVD, M2D, ED, kD,
2608 write_accelMDG_rest_file(step,
'D', V_n, VmaxD, VminD, VavgD, sigmaVD, M2D, ED, kD,
2612 write_accelMDG_rest_file(step,
'D', V_n, VmaxD, VminD, VavgD, sigmaVD, M2D, ED, kD,
2620 VmaxD = VminD = VavgD = testV;
2624 else if(step == ntcmdprep && ntcmdprep != 0){
2626 VmaxD = VminD = VavgD = testV;
2628 iout <<
"GAUSSIAN ACCELERATED MD: RESET DIHEDRAL POTENTIAL STATISTICS\n" <<
endi;
2631 else if(ntave > 0 && step % ntave == 0){
2633 if(testV > VmaxD) VmaxD = testV;
2634 if(testV < VminD) VminD = testV;
2638 iout <<
"GAUSSIAN ACCELERATED MD: RESET DIHEDRAL POTENTIAL AVG AND STD\n" <<
endi;
2643 &VmaxD, &VminD, &VavgD, &M2D, &sigmaVD) ;
2646 if(step == ntcmd - 1){
2648 VmaxD, VminD, VavgD, sigmaVD, &k0D, &kD, &ED, &iEusedD, &warnD);
2652 else if(step < nteb){
2654 &dV, &factor_dihe, &vir);
2659 VmaxD = VminD = VavgD = testV;
2661 iout <<
"GAUSSIAN ACCELERATED MD: RESET DIHEDRAL POTENTIAL STATISTICS\n" <<
endi;
2664 else if(ntave > 0 && step % ntave == 0){
2666 if(testV > VmaxD) VmaxD = testV;
2667 if(testV < VminD) VminD = testV;
2671 iout <<
"GAUSSIAN ACCELERATED MD: RESET DIHEDRAL POTENTIAL AVG AND STD\n" <<
endi;
2675 &VmaxD, &VminD, &VavgD, &M2D, &sigmaVD) ;
2678 if(step >= ntebprep)
2679 if((ntave > 0 && step % ntave == 0) || ntave <= 0)
2681 VmaxD, VminD, VavgD, sigmaVD, &k0D, &kD, &ED, &iEusedD, &warnD);
2686 &dV, &factor_dihe, &vir);
2692 Tensor vir_tot = vir_normal + vir_nbond + vir_slow;
2693 testV = potentialEnergy;
2695 testV -= dihedralEnergy + crosstermEnergy;
2696 vir_tot -= vir_dihe;
2700 if(step > firststep && step % restfreq == 0)
2701 write_accelMDG_rest_file(step,
'T', V_n, VmaxP, VminP, VavgP, sigmaVP, M2P, EP, kP,
2706 write_accelMDG_rest_file(step,
'T', V_n, VmaxP, VminP, VavgP, sigmaVP, M2P, EP, kP,
2710 write_accelMDG_rest_file(step,
'T', V_n, VmaxP, VminP, VavgP, sigmaVP, M2P, EP, kP,
2718 VmaxP = VminP = VavgP = testV;
2722 else if(step == ntcmdprep && ntcmdprep != 0){
2724 VmaxP = VminP = VavgP = testV;
2726 iout <<
"GAUSSIAN ACCELERATED MD: RESET TOTAL POTENTIAL STATISTICS\n" <<
endi;
2729 else if(ntave > 0 && step % ntave == 0){
2731 if(testV > VmaxP) VmaxP = testV;
2732 if(testV < VminP) VminP = testV;
2736 iout <<
"GAUSSIAN ACCELERATED MD: RESET TOTAL POTENTIAL AVG AND STD\n" <<
endi;
2741 &VmaxP, &VminP, &VavgP, &M2P, &sigmaVP);
2743 if(step == ntcmd - 1){
2745 VmaxP, VminP, VavgP, sigmaVP, &k0P, &kP, &EP, &iEusedP, &warnP);
2749 else if(step < nteb){
2751 &dV, &factor_tot, &vir);
2756 VmaxP = VminP = VavgP = testV;
2758 iout <<
"GAUSSIAN ACCELERATED MD: RESET TOTAL POTENTIAL STATISTICS\n" <<
endi;
2761 else if(ntave > 0 && step % ntave == 0){
2763 if(testV > VmaxP) VmaxP = testV;
2764 if(testV < VminP) VminP = testV;
2768 iout <<
"GAUSSIAN ACCELERATED MD: RESET TOTAL POTENTIAL AVG AND STD\n" <<
endi;
2772 &VmaxP, &VminP, &VavgP, &M2P, &sigmaVP);
2775 if(step >= ntebprep)
2776 if((ntave > 0 && step % ntave == 0) || ntave <= 0)
2778 VmaxP, VminP, VavgP, sigmaVP, &k0P, &kP, &EP, &iEusedP, &warnP);
2783 &dV, &factor_tot, &vir);
2790 if((ntcmdprep > 0 && step == ntcmdprep) ||
2791 (step < nteb && ntave > 0 && step % ntave == 0) ||
2804 testV = dihedralEnergy + crosstermEnergy;
2805 if ( testV < accelMDE ) {
2806 factor_dihe = accelMDalpha/(accelMDalpha + accelMDE - testV);
2807 factor_dihe *= factor_dihe;
2808 vir = vir_dihe * (factor_dihe - 1.0);
2809 dV = (accelMDE - testV)*(accelMDE - testV)/(accelMDalpha + accelMDE - testV);
2815 testV = dihedralEnergy + crosstermEnergy;
2816 if ( testV < accelMDE ) {
2817 factor_dihe = accelMDalpha/(accelMDalpha + accelMDE - testV);
2818 factor_dihe *= factor_dihe;
2819 vir = vir_dihe * (factor_dihe - 1.0) ;
2820 dV = (accelMDE - testV)*(accelMDE - testV)/(accelMDalpha + accelMDE - testV);
2823 testV = potentialEnergy - dihedralEnergy - crosstermEnergy;
2824 if ( testV < accelMDTE ) {
2825 factor_tot = accelMDTalpha/(accelMDTalpha + accelMDTE - testV);
2826 factor_tot *= factor_tot;
2827 vir += (vir_normal - vir_dihe + vir_nbond + vir_slow) * (factor_tot - 1.0);
2828 dV += (accelMDTE - testV)*(accelMDTE - testV)/(accelMDTalpha + accelMDTE - testV);
2834 testV = potentialEnergy;
2835 if ( testV < accelMDE ) {
2836 factor_tot = accelMDalpha/(accelMDalpha + accelMDE - testV);
2837 factor_tot *= factor_tot;
2838 vir = (vir_normal + vir_nbond + vir_slow) * (factor_tot - 1.0);
2839 dV = (accelMDE - testV)*(accelMDE - testV)/(accelMDalpha + accelMDE - testV);
2845 accelMDfactor[0]=factor_dihe;
2846 accelMDfactor[1]=factor_tot;
2852 if ( factor_tot < 0.001 ) {
2853 iout <<
iWARN <<
"accelMD is using a very high boost potential, simulation may become unstable!" 2857 if ( ! (step % accelMDOutFreq) ) {
2861 iout <<
"ACCELERATED MD: STEP " << step
2864 <<
" BOND " << bondEnergy
2865 <<
" ANGLE " << angleEnergy
2866 <<
" DIHED " << dihedralEnergy+crosstermEnergy
2867 <<
" IMPRP " << improperEnergy
2868 <<
" ELECT " << amd_electEnergy+amd_electEnergySlow
2869 <<
" VDW " << amd_ljEnergy
2870 <<
" POTENTIAL " << potentialEnergy;
2871 if(amd_ljEnergy_Corr)
2872 iout <<
" LJcorr " << amd_ljEnergy_Corr;
2876 iout <<
"GAUSSIAN ACCELERATED MD: DIHED iE " << iEusedD
2877 <<
" Vmax " << VmaxD <<
" Vmin " << VminD
2878 <<
" Vavg " << VavgD <<
" sigmaV " << sigmaVD
2879 <<
" E " << ED <<
" k0 " << k0D <<
" k " << kD
2882 iout <<
"GAUSSIAN ACCELERATED MD: DIHED !!! WARNING: k0 > 1, " 2883 <<
"SWITCHED TO iE = 1 MODE !!!\n" <<
endi;
2885 iout <<
"GAUSSIAN ACCELERATED MD: DIHED !!! WARNING: k0 <= 0, " 2886 <<
"SWITCHED TO iE = 1 MODE !!!\n" <<
endi;
2888 iout <<
"GAUSSIAN ACCELERATED MD: TOTAL iE " << iEusedP
2889 <<
" Vmax " << VmaxP <<
" Vmin " << VminP
2890 <<
" Vavg " << VavgP <<
" sigmaV " << sigmaVP
2891 <<
" E " << EP <<
" k0 " << k0P <<
" k " << kP
2894 iout <<
"GAUSSIAN ACCELERATED MD: TOTAL !!! WARNING: k0 > 1, " 2895 <<
"SWITCHED TO iE = 1 MODE !!!\n" <<
endi;
2897 iout <<
"GAUSSIAN ACCELERATED MD: TOTAL !!! WARNING: k0 <= 0, " 2898 <<
"SWITCHED TO iE = 1 MODE !!!\n" <<
endi;
2899 warnD = warnP =
'\0';
2909 if ( volume != 0. ) {
2910 p_normal = vir_normal/volume;
2911 p_nbond = vir_nbond/volume;
2912 p_slow = vir_slow/volume;
2915 iout <<
" accelMD Scaling Factor: " << accelMDfactor <<
"\n" 2916 <<
" accelMD PNORMAL " << trace(p_normal)*
PRESSUREFACTOR/3. <<
"\n" 2917 <<
" accelMD PNBOND " << trace(p_nbond)*
PRESSUREFACTOR/3. <<
"\n" 2927 iout <<
iINFO <<
"INITIALISING ADAPTIVE TEMPERING\n" <<
endi;
2932 iout <<
iINFO <<
"READING ADAPTIVE TEMPERING RESTART FILE\n" <<
endi;
2934 if (adaptTempRead) {
2939 adaptTempRead >> readInt;
2972 adaptTempRead.close();
2974 else NAMD_die(
"Could not open ADAPTIVE TEMPERING restart file.\n");
3014 iout <<
iINFO <<
"ADAPTIVE TEMPERING: ADAPTIVE BIN AVERAGING PARAMETER: " <<
adaptTempCg <<
"\n";
3015 iout <<
iINFO <<
"ADAPTIVE TEMPERING: SCALING CONSTANT FOR LANGEVIN TEMPERATURE UPDATES: " <<
adaptTempDt <<
"\n";
3021 NAMD_die(
"Error opening restart file");
3028 iout <<
"ADAPTEMP: WRITING RESTART FILE AT STEP " << step <<
"\n" <<
endi;
3060 if (
minimize || (step < simParams->adaptTempFirstStep ) ||
3070 <<
" adaptTempBeta: " << adaptTempBeta
3089 potEnergyVariance -= potEnergyAverage*potEnergyAverage;
3102 BigReal deltaBeta = 0.04*adaptTempBeta;
3109 betaMinus = adaptTempBeta - deltaBeta;
3110 betaPlus = adaptTempBeta + deltaBeta;
3111 BigReal betaMinus2 = betaMinus*betaMinus;
3112 BigReal betaPlus2 = betaPlus*betaPlus;
3130 DeltaE2AveSum += DeltaE2Ave;
3139 A0 /= (betaNp1-betaMinus);
3146 DeltaE2AveSum += DeltaE2Ave;
3157 A1 /= (betaPlus-betaNp1);
3166 if (DeltaE2Ave != 0) {
3167 aplus = (A0-A2)/(A0-A1);
3177 potEnergyAve0 /= (betaNp1-betaMinus);
3182 potEnergyAve1 /= (betaPlus-betaNp1);
3183 potEnergyAverage = aminus*potEnergyAve0;
3184 potEnergyAverage += aplus*potEnergyAve1;
3187 iout <<
"ADAPTEMP DEBUG:" <<
"\n" 3190 <<
" adaptTempBeta: " << adaptTempBeta <<
"\n" 3194 <<
" gammaAve: " << gammaAve <<
"\n" 3195 <<
" deltaBeta: " << deltaBeta <<
"\n" 3196 <<
" betaMinus: " << betaMinus <<
"\n" 3197 <<
" betaPlus: " << betaPlus <<
"\n" 3198 <<
" nMinus: " << nMinus <<
"\n" 3199 <<
" nPlus: " << nPlus <<
"\n" 3200 <<
" A0: " << A0 <<
"\n" 3201 <<
" A1: " << A1 <<
"\n" 3202 <<
" A2: " << A2 <<
"\n" 3203 <<
" a+: " << aplus <<
"\n" 3204 <<
" a-: " << aminus <<
"\n" 3210 iout <<
"ADAPTEMP DEBUG:" <<
"\n" 3213 <<
" adaptTempBeta: " << adaptTempBeta <<
"\n" 3217 <<
" gammaAve: " << gammaAve <<
"\n" 3218 <<
" deltaBeta: N/A\n" 3219 <<
" betaMinus: N/A\n" 3220 <<
" betaPlus: N/A\n" 3279 if (adaptTempTdiff > 0) {
3291 iout <<
"ADAPTEMP: " << step <<
" FRAC " << Frac <<
"\n";
3294 iout <<
"ADAPTEMP: Updating adaptTempDt to ";
3307 if (adaptTempTdiff > 0) {
3319 iout <<
"ADAPTEMP: " << step <<
" FRAC " << Frac <<
"\n";
3322 iout <<
"ADAPTEMP: Updating adaptTempDt to ";
3337 if ( ! (step % adaptTempOutFreq) ) {
3338 iout <<
"ADAPTEMP: STEP " << step
3340 <<
" ENERGY " << std::setprecision(10) << potentialEnergy
3341 <<
" ENERGYAVG " << std::setprecision(10) << potEnergyAverage
3342 <<
" ENERGYVAR " << std::setprecision(10) << potEnergyVariance;
3358 #ifdef NODEGROUP_FORCE_REGISTER 3364 #ifdef NODEGROUP_FORCE_REGISTER 3367 if ( ((
int)checksum) != molecule->
numAtoms && step != 0) {
3368 sprintf(errmsg,
"Bad global atom count! (%d vs %d)\n",
3369 (
int)checksum, molecule->
numAtoms);
3373 #ifdef NODEGROUP_FORCE_REGISTER 3378 #ifdef NODEGROUP_FORCE_REGISTER 3384 computesPartitioned = 0;
3385 }
else if ( (
int)checksum >
computeChecksum && ! computesPartitioned ) {
3386 computesPartitioned = 1;
3388 NAMD_bug(
"Bad global compute count!\n");
3392 #ifdef NODEGROUP_FORCE_REGISTER 3397 #ifdef NODEGROUP_FORCE_REGISTER 3400 if ( checksum_b && (((
int)checksum) != molecule->
numCalcBonds) ) {
3401 sprintf(errmsg,
"Bad global bond count! (%d vs %d)\n",
3403 if ( forgiving && (((
int)checksum) < molecule->
numCalcBonds) )
3408 #ifdef NODEGROUP_FORCE_REGISTER 3415 #ifdef NODEGROUP_FORCE_REGISTER 3418 if ( checksum_b && (((
int)checksum) != molecule->
numCalcAngles) ) {
3419 sprintf(errmsg,
"Bad global angle count! (%d vs %d)\n",
3421 if ( forgiving && (((
int)checksum) < molecule->
numCalcAngles) )
3425 #ifdef NODEGROUP_FORCE_REGISTER 3430 #ifdef NODEGROUP_FORCE_REGISTER 3434 sprintf(errmsg,
"Bad global dihedral count! (%d vs %d)\n",
3441 #ifdef NODEGROUP_FORCE_REGISTER 3446 #ifdef NODEGROUP_FORCE_REGISTER 3450 sprintf(errmsg,
"Bad global improper count! (%d vs %d)\n",
3457 #ifdef NODEGROUP_FORCE_REGISTER 3462 #ifdef NODEGROUP_FORCE_REGISTER 3465 if ( checksum_b && (((
int)checksum) != molecule->
numCalcTholes) ) {
3466 sprintf(errmsg,
"Bad global Thole count! (%d vs %d)\n",
3468 if ( forgiving && (((
int)checksum) < molecule->
numCalcTholes) )
3473 #ifdef NODEGROUP_FORCE_REGISTER 3478 #ifdef NODEGROUP_FORCE_REGISTER 3481 if ( checksum_b && (((
int)checksum) != molecule->
numCalcAnisos) ) {
3482 sprintf(errmsg,
"Bad global Aniso count! (%d vs %d)\n",
3484 if ( forgiving && (((
int)checksum) < molecule->
numCalcAnisos) )
3489 #ifdef NODEGROUP_FORCE_REGISTER 3494 #ifdef NODEGROUP_FORCE_REGISTER 3498 sprintf(errmsg,
"Bad global crossterm count! (%d vs %d)\n",
3505 #ifdef NODEGROUP_FORCE_REGISTER 3510 #ifdef NODEGROUP_FORCE_REGISTER 3516 iout <<
iWARN <<
"High global exclusion count (" 3517 << ((
int64)checksum) <<
" vs " 3519 <<
iWARN <<
"This warning is not unusual during minimization.\n" 3520 <<
iWARN <<
"Decreasing pairlistdist or cutoff that is too close to periodic cell size may avoid this.\n" <<
endi;
3523 sprintf(errmsg,
"High global exclusion count! (%lld vs %lld) System unstable or pairlistdist or cutoff too close to periodic cell size.\n",
3532 iout <<
iWARN <<
"Low global exclusion count! (" 3536 <<
iWARN <<
"This warning is not unusual during minimization.\n" 3537 <<
iWARN <<
"Increasing pairlistdist or cutoff may avoid this.\n";
3539 iout <<
" System unstable or pairlistdist or cutoff too small.\n";
3544 sprintf(errmsg,
"Low global exclusion count! (%lld vs %lld) System unstable or pairlistdist or cutoff too small.\n",
3550 #if defined(NAMD_CUDA) || defined(NAMD_HIP) 3551 #ifdef NODEGROUP_FORCE_REGISTER 3556 #ifdef NODEGROUP_FORCE_REGISTER 3562 iout <<
iWARN <<
"High global CUDA exclusion count (" 3563 << ((
int64)checksum) <<
" vs " 3565 <<
iWARN <<
"This warning is not unusual during minimization.\n" 3566 <<
iWARN <<
"Decreasing pairlistdist or cutoff that is too close to periodic cell size may avoid this.\n" <<
endi;
3569 sprintf(errmsg,
"High global CUDA exclusion count! (%lld vs %lld) System unstable or pairlistdist or cutoff too close to periodic cell size.\n",
3578 iout <<
iWARN <<
"Low global CUDA exclusion count! (" 3582 <<
iWARN <<
"This warning is not unusual during minimization.\n" 3583 <<
iWARN <<
"Increasing pairlistdist or cutoff may avoid this.\n";
3585 iout <<
" System unstable or pairlistdist or cutoff too small.\n";
3590 sprintf(errmsg,
"Low global CUDA exclusion count! (%lld vs %lld) System unstable or pairlistdist or cutoff too small.\n",
3597 #ifdef NODEGROUP_FORCE_REGISTER 3602 #ifdef NODEGROUP_FORCE_REGISTER 3606 iout <<
iERROR <<
"Margin is too small for " << ((int)checksum) <<
3607 " atoms during timestep " << step <<
".\n" <<
iERROR <<
3608 "Incorrect nonbonded forces and energies may be calculated!\n" <<
endi;
3612 #ifdef NODEGROUP_FORCE_REGISTER 3617 #ifdef NODEGROUP_FORCE_REGISTER 3622 "Pairlistdist is too small for " << ((int)checksum) <<
3623 " computes during timestep " << step <<
".\n" <<
endi;
3627 #ifdef NODEGROUP_FORCE_REGISTER 3632 #ifdef NODEGROUP_FORCE_REGISTER 3637 iout <<
iWARN <<
"Stray PME grid charges ignored!\n" <<
endi;
3638 else NAMD_bug(
"Stray PME grid charges detected!\n");
3647 const double endCTime = CmiTimer() - firstCTime;
3650 if ( (((
int)endWTime)>>6) != (((
int)startWTime)>>6 ) )
fflush_count = 2;
3652 const double elapsedW =
3654 const double elapsedC =
3657 const double remainingW = elapsedW * (
simParams->
N - step);
3658 const double remainingW_hours = remainingW / 3600;
3660 startWTime = endWTime;
3661 startCTime = endCTime;
3663 #if defined(NAMD_CUDA) || defined(NAMD_HIP) 3666 iout <<
iWARN <<
"Energy evaluation is expensive, increase computeEnergies to improve performance.\n" <<
endi;
3667 iout <<
iWARN <<
"If computeEnergies is not defined, its value defaults to the same as outputEnergies, and increasing outputEnergies would be helpful to improve the performance.\n" <<
endi;
3673 BigReal days = 1.0 / (24.0 * 60.0 * 60.0);
3674 BigReal nsPerDay = ns / (elapsedW * days);
3675 CmiPrintf(
"TIMING: %d CPU: %g, %g/step Wall: %g, %g/step" 3677 ", %g hours remaining, %f MB of memory in use.\n",
3678 step, endCTime, elapsedC, endWTime, elapsedW,
3683 CmiPrintf(
"TIMING: %d CPU: %g, %g/step Wall: %g, %g/step" 3684 ", %g hours remaining, %f MB of memory in use.\n",
3685 step, endCTime, elapsedC, endWTime, elapsedW,
3695 double ns_per_day = (
simParams->
dt / t_avg) * (60 * 60 * 24 * 1e-6);
3696 CmiPrintf(
"PERFORMANCE: %d averaging %g ns/day, %g sec/step with" 3697 " standard deviation %g\n", step, ns_per_day, t_avg, t_std);
3701 if (benchmarkTime > 0 && benchmarkTime <= endWTime) {
3704 sprintf(s,
"Exceeded benchmark time limit %g seconds",
3729 #ifdef DEBUG_MINIMIZE 3730 printf(
"%s, line %d\n", __FILE__, __LINE__);
3731 printf(
"Step %d:\n", step);
3755 BigReal totalPotentialEnergy = 0.0;
3756 BigReal bondEnergy, angleEnergy, dihedralEnergy;
3757 BigReal improperEnergy, crosstermEnergy;
3762 #ifdef NODEGROUP_FORCE_REGISTER 3780 #ifdef NODEGROUP_FORCE_REGISTER 3784 totalPotentialEnergy += bondEnergy + angleEnergy + dihedralEnergy +
3785 improperEnergy + crosstermEnergy + miscEnergy
3790 #ifdef NODEGROUP_FORCE_REGISTER 3799 #ifdef NODEGROUP_FORCE_REGISTER 3807 #ifdef NODEGROUP_FORCE_REGISTER 3813 #ifdef NODEGROUP_FORCE_REGISTER 3820 #ifdef MEM_OPT_VERSION 3821 NAMD_bug(
"LJcorrection not supported in memory optimized build.");
3830 printf(
"MC-data (TotalPotentialEnergy): step: %d, Pe: %d, Bond: %f, Angle: %f, Dih: %f, IMP: %f, Cross: %f\n",
3831 step, CkMyPe(), bondEnergy, angleEnergy, dihedralEnergy, improperEnergy, crosstermEnergy);
3832 printf(
"MC-data (TotalPotentialEnergy): step: %d, Pe: %d, LJ: %f, Real: %f, Recip: %f\n", step, CkMyPe(),
3834 printf(
"MC-data (TotalPotentialEnergy): step: %d, Pe: %d, MISC: %f, BC: %f, TOTAL: %f\n", step, CkMyPe(),
3835 miscEnergy, bcEnergy, totalPotentialEnergy);
3838 return totalPotentialEnergy;
3867 #ifdef NODEGROUP_FORCE_REGISTER 3885 #ifdef NODEGROUP_FORCE_REGISTER 3891 #ifdef NODEGROUP_FORCE_REGISTER 3939 #ifdef NODEGROUP_FORCE_REGISTER 3947 #ifdef NODEGROUP_FORCE_REGISTER 3967 #ifdef NODEGROUP_FORCE_REGISTER 3974 #ifdef MEM_OPT_VERSION 3975 NAMD_bug(
"LJcorrection not supported in memory optimized build.");
4004 #ifdef NODEGROUP_FORCE_REGISTER 4017 angularMomentum.
x =
reduction->
item(REDUCTION_ANGULAR_MOMENTUM_X);
4018 angularMomentum.
y =
reduction->
item(REDUCTION_ANGULAR_MOMENTUM_Y);
4019 angularMomentum.
z =
reduction->
item(REDUCTION_ANGULAR_MOMENTUM_Z);
4020 #ifdef NODEGROUP_FORCE_REGISTER 4025 potentialEnergy = (bondEnergy + angleEnergy + dihedralEnergy
4027 + crosstermEnergy + boundaryEnergy + miscEnergy +
goTotalEnergy 4031 if ( ! cudaIntegrator) {
4048 if (step <= simParams->firstTimestep &&
4056 iout <<
"MOMENTUM: " << step
4057 <<
" P: " << momentum
4058 <<
" L: " << angularMomentum
4063 if ( ! cudaIntegrator ) {
4069 if (cudaIntegrator) {
4105 iout <<
"PRESSURE: " << step <<
" " 4107 <<
"GPRESSURE: " << step <<
" " 4113 <<
"GPRESSAVG: " << step <<
" " 4119 iout <<
"PRESSURE: " << step <<
" " 4121 <<
"GPRESSURE: " << step <<
" " 4125 <<
"GPRESSAVG: " << step <<
" " 4141 memset(total, 0, arraysize*
sizeof(
BigReal));
4150 if (!(step % freq)) {
4154 iout <<
"PRESSUREPROFILE: " << step <<
" ";
4155 if (step == first) {
4157 for (
int i=0; i<arraysize; i++) {
4158 iout << total[i] * scalefac <<
" ";
4163 for (
int i=0; i<arraysize; i++)
4178 if ( benchPhase > 0 && benchPhase < 10 ) {
4179 #if defined(NAMD_CUDA) || defined(NAMD_HIP) 4181 iout <<
iWARN <<
"Energy evaluation is expensive, increase computeEnergies to improve performance.\n";
4182 iout <<
iWARN <<
"If computeEnergies is not defined, its value defaults to the same as outputEnergies, and increasing outputEnergies would be helpful to improve the performance.\n";
4186 if ( benchPhase < 4 )
iout <<
"Initial time: ";
4187 else iout <<
"Benchmark time: ";
4188 iout << CkNumPes() <<
" CPUs ";
4194 BigReal days = 1.0 / (24.0 * 60.0 * 60.0);
4197 BigReal nsPerDay = ns / (wallPerStep * days);
4198 iout << wallPerStep <<
" s/step " << nsPerDay <<
" ns/day ";
4201 BigReal daysPerNano = wallPerStep * days / ns;
4202 iout << wallPerStep <<
" s/step " << daysPerNano <<
" days/ns ";
4215 Vector pairVDWForce, pairElectForce;
4217 #ifdef NODEGROUP_FORCE_REGISTER 4225 #ifdef NODEGROUP_FORCE_REGISTER 4237 #define CALLBACKDATA(LABEL,VALUE) \ 4238 labels << (LABEL) << " "; values << (VALUE) << " "; 4239 #define CALLBACKLIST(LABEL,VALUE) \ 4240 labels << (LABEL) << " "; values << "{" << (VALUE) << "} "; 4243 std::ostringstream labels, values;
4246 values.precision(16);
4248 values << std::setprecision(16);
4271 labels <<
"PERIODIC "; values <<
"{" << lattice.
a_p() <<
" " 4272 << lattice.
b_p() <<
" " << lattice.
c_p() <<
"} ";
4298 CALLBACKLIST(
"BOND2", bondEnergy + angleEnergy + dihedralEnergy +
4305 labels <<
'\0'; values <<
'\0';
4306 state->callback_labelstring = labels.str();
4307 state->callback_valuestring = values.str();
4313 if ( ! cudaIntegrator) {
4325 " steps.\n" <<
endi;
4343 if (cudaIntegrator) {
4354 energies.
tstep = step;
4357 energies.
Epot = potentialEnergy;
4360 energies.
Ebond = bondEnergy;
4361 energies.
Eangle = angleEnergy;
4362 energies.
Edihe = dihedralEnergy + crosstermEnergy;
4363 energies.
Eimpr = improperEnergy;
4369 " margin violations detected since previous energy output.\n" <<
endi;
4376 iout <<
"ETITLE: TS";
4392 if (cudaIntegrator) {
4400 if ( volume != 0. ) {
4438 iout <<
"\nTITITLE: TS";
4452 iout <<
"\nFEPTITLE: TS";
4465 iout <<
"QMETITLE: TS";
4480 iout <<
FORMAT(dihedralEnergy+crosstermEnergy, precision);
4496 if (cudaIntegrator) {
4511 if (cudaIntegrator) {
4568 iout <<
FORMAT(bondEnergy + angleEnergy + dihedralEnergy