46 #if(CMK_CCS_AVAILABLE && CMK_WEB_MODE)
47 extern "C" void CApplicationDepositNode0Data(
char *);
52 #define cbrt(x) pow(x,(double)(1.0/3.0))
56 #define MIN_DEBUG_LEVEL 3
60 #define XXXBIGREAL 1.0e32
75 const char *myname,
int outputfreq)
76 : nslabs(numslabs), freq(outputfreq) {
77 name = strdup(myname);
78 nelements = 3*nslabs * numpartitions;
81 average =
new BigReal[nelements];
95 double inv_volume = 1.0 / lattice.
volume();
97 int arraysize = 3*nslabs;
98 for (i=0; i<nelements; i++) {
101 total[i % arraysize] += val;
104 if (!(step % freq)) {
108 iout <<
"PPROFILE" << name <<
": " << step <<
" ";
109 if (step == firsttimestep) {
111 for (i=0; i<nelements; i++)
112 iout << reduction->
item(i)*scalefac*inv_volume <<
" ";
116 for (i=0; i<nelements; i++)
117 iout << average[i]*scalefac <<
" ";
121 memset(average, 0, nelements*
sizeof(
BigReal));
129 computeChecksum(0), marginViolations(0), pairlistWarnings(0),
134 firstCTime(CmiTimer()),
136 firstWTime(CmiWallTimer()),
165 int npairs = (ntypes * (ntypes+1))/2;
171 nslabs, npairs,
"BONDED", freq);
173 nslabs, npairs,
"NONBONDED", freq);
176 nslabs, ntypes,
"INTERNAL", freq);
180 nslabs, npairs,
"NONBONDED", freq);
202 #define AVGXY(T) T.xy = T.yx = 0; T.xx = T.yy = 0.5 * ( T.xx + T.yy );\
203 T.xz = T.zx = T.yz = T.zy = 0.5 * ( T.xz + T.yz );
220 for (
int i=0;i < n;i++) {
271 DebugM(4,
"Starting thread in controller on this=" <<
this <<
"\n");
272 thread = CthCreate((CthVoidFn)&(threadRun),(
void*)(
this),
CTRL_STK_SZ);
273 CthSetStrategyDefault(thread);
274 #if CMK_BLUEGENE_CHARM
287 switch ( scriptTask ) {
321 NAMD_die(
"Unable to revert, checkpoint was never called!");
349 NAMD_die(
"Unable to load checkpoint, requested key was never stored.");
366 NAMD_die(
"Unable to swap checkpoint, requested key was never stored.");
383 NAMD_die(
"Unable to free checkpoint, requested key was never stored.");
404 NAMD_bug(
"Unknown task in Controller::algorithm");
423 if (sig == SIGINT) gotsigint = 1;
454 if ( zeroMomentum && dofull && ! (step %
slowFreq) )
462 traceUserEvent(eventEndOfTimeStep);
463 sprintf(traceNote,
"s:%d", step);
464 traceUserSuppliedNote(traceNote);
476 for ( ++step ; step <= numberOfSteps; ++step )
488 if ( zeroMomentum && dofull && ! (step %
slowFreq) )
502 traceUserEvent(eventEndOfTimeStep);
503 sprintf(traceNote,
"s:%d", step);
504 traceUserSuppliedNote(traceNote);
530 #ifdef MEASURE_NAMD_WITH_PAPI
532 int bstep =
simParams->papiMeasureStartStep;
533 int estep = bstep +
simParams->numPapiMeasureSteps;
535 papiMeasureBarrier(1, step);
538 papiMeasureBarrier(0, step);
554 printMinimizeEnergies(step); \
555 outputExtendedSystem(step); \
556 rebalanceLoad(step); \
557 if ( step == numberOfSteps ) return; \
561 if ( step == numberOfSteps ) { \
562 if ( minVerbose ) { iout << "LINE MINIMIZER: RETURNING TO " << mid.x << " FROM " << last.x << "\n" << endi; } \
563 if ( newDir || (mid.x-last.x) ) { \
564 broadcast->minimizeCoefficient.publish(minSeq++,mid.x-last.x); \
566 broadcast->minimizeCoefficient.publish(minSeq++,0.); \
567 broadcast->minimizeCoefficient.publish(minSeq++,0.); \
568 min_reduction->require(); \
569 broadcast->minimizeCoefficient.publish(minSeq++,0.); \
571 enqueueCollections(step); \
573 } else if ( (X)-last.x ) { \
574 broadcast->minimizeCoefficient.publish(minSeq++,(X)-last.x); \
577 enqueueCollections(step); \
579 last.u = min_energy; \
580 last.dudx = -1. * min_f_dot_v; \
581 last.noGradient = min_huge_count; \
582 if ( minVerbose ) { \
583 iout << "LINE MINIMIZER: POSITION " << last.x << " ENERGY " << last.u << " GRADIENT " << last.dudx; \
584 if ( last.noGradient ) iout << " HUGECOUNT " << last.noGradient; \
585 iout << "\n" << endi; \
602 const BigReal goldenRatio = 0.5 * ( sqrt(5.0) - 1.0 );
609 int old_min_huge_count = 2000000000;
610 int steps_at_same_min_huge_count = 0;
611 for (
int i_slow = 0;
min_huge_count && i_slow < 100; ++i_slow ) {
613 if ( ++steps_at_same_min_huge_count > 10 )
break;
616 steps_at_same_min_huge_count = 0;
626 iout <<
"MINIMIZER STARTING CONJUGATE GRADIENT ALGORITHM\n" <<
endi;
629 int errorFactor = 10;
648 iout <<
"LINE MINIMIZER: POSITION " << last.
x <<
" ENERGY " << last.
u <<
" GRADIENT " << last.
dudx;
653 iout <<
"LINE MINIMIZER REDUCING GRADIENT FROM " <<
654 fabs(min_f_dot_v) <<
" TO " << tol <<
"\n" <<
endi;
660 while ( last.
u < mid.
u ) {
661 if ( last.
dudx < mid.
dudx * (5.5 - x/maxstep)/5. ) {
662 x = 6 * maxstep;
break;
664 lo = mid; mid = last;
666 if ( x > 5.5 * maxstep )
break;
669 if ( x > 5.5 * maxstep ) {
670 iout <<
"MINIMIZER RESTARTING CONJUGATE GRADIENT ALGORITHM DUE TO POOR PROGRESS\n" <<
endi;
680 #define PRINT_BRACKET \
681 iout << "LINE MINIMIZER BRACKET: DX " \
682 << (mid.x-lo.x) << " " << (hi.x-mid.x) << \
683 " DU " << (mid.u-lo.u) << " " << (hi.u-mid.u) << " DUDX " << \
684 lo.dudx << " " << mid.dudx << " " << hi.dudx << " \n" << endi;
688 for ( itcnt = 10; itcnt > 0; --itcnt ) {
692 if ( ( mid.
x - lo.
x ) > ( hi.
x - mid.
x ) ) {
693 x = (1.0 - goldenRatio) * lo.
x + goldenRatio * mid.
x;
695 if ( last.
u <= mid.
u ) {
696 hi = mid; mid = last;
701 x = (1.0 - goldenRatio) * hi.
x + goldenRatio * mid.
x;
703 if ( last.
u <= mid.
u ) {
704 lo = mid; mid = last;
710 if ( mid.
dudx > 0. ) {
711 BigReal altxhi = 0.1 * lo.
x + 0.9 * mid.
x;
712 BigReal altxlo = 0.9 * lo.
x + 0.1 * mid.
x;
713 x = mid.
dudx*(mid.
x*mid.
x-lo.
x*lo.
x) + 2*mid.
x*(lo.
u-mid.
u);
715 if ( xdiv ) x /= xdiv;
else x = last.
x;
716 if ( x > altxhi ) x = altxhi;
717 if ( x < altxlo ) x = altxlo;
718 if ( x-last.
x == 0 )
break;
720 if ( last.
u <= mid.
u ) {
721 hi = mid; mid = last;
723 if ( lo.
dudx < 0. && last.
dudx > 0. ) progress = 0;
727 BigReal altxlo = 0.1 * hi.
x + 0.9 * mid.
x;
728 BigReal altxhi = 0.9 * hi.
x + 0.1 * mid.
x;
729 x = mid.
dudx*(mid.
x*mid.
x-hi.
x*hi.
x) + 2*mid.
x*(hi.
u-mid.
u);
731 if ( xdiv ) x /= xdiv;
else x = last.
x;
732 if ( x < altxlo ) x = altxlo;
733 if ( x > altxhi ) x = altxhi;
734 if ( x-last.
x == 0 )
break;
736 if ( last.
u <= mid.
u ) {
737 lo = mid; mid = last;
739 if ( hi.
dudx > 0. && last.
dudx < 0. ) progress = 0;
749 if ( fabs(last.
dudx) < tol )
break;
750 if ( lo.
x != mid.
x && (lo.
u-mid.
u) < tol * (mid.
x-lo.
x) )
break;
751 if ( mid.
x != hi.
x && (hi.
u-mid.
u) < tol * (hi.
x-mid.
x) )
break;
756 c = ( c > 1.5 ? 1.5 : c );
757 if ( atStart ) { c = 0; --atStart; }
760 if ( errorFactor < 100 ) errorFactor += 10;
763 iout <<
"MINIMIZER RESTARTING CONJUGATE GRADIENT ALGORITHM\n" <<
endi;
783 BigReal NGfac = 1.0/sqrt(3.0*NG + 1.0);
794 BigReal velScale = exp(-s*NGfac*(1.0 + 1.0/NG)*multigratorXiT);
800 if (callNumber == 1) {
831 calcPressure(step, 0, virial_normal, virial_nbond, virial_slow,
832 intVirial_normal, intVirial_nbond, intVirial_slow,
833 extForce_normal, extForce_nbond, extForce_slow);
834 if (callNumber == 2) {
864 for (
int i=n-1;i >= 0;i--) {
868 }
else if (i == n-1) {
900 for (
int i=0;i < n;i++) {
904 }
else if (i == n-1) {
938 BigReal W = (3.0*NG + 1.0)*kT0*tauV*tauV;
966 #define LIMIT_SCALING(VAR,MIN,MAX,FLAG) {\
967 if ( VAR < (MIN) ) { VAR = (MIN); FLAG = 1; } \
968 if ( VAR > (MAX) ) { VAR = (MAX); FLAG = 1; } }
976 " cell rescaling factor limited.\n" <<
endi;
999 #ifdef DEBUG_PRESSURE
1000 iout <<
iINFO <<
"entering langevinPiston1, strain rate: " << strainRate <<
"\n";
1006 BigReal f1 = exp( -0.5 * dt_long * gamma );
1007 BigReal f2 = sqrt( ( 1. - f1*f1 ) * kT / mass );
1023 #ifdef DEBUG_PRESSURE
1024 iout <<
iINFO <<
"applying langevin, strain rate: " << strainRate <<
"\n";
1038 strainRate += ( 0.5 * dt * cellDims *
state->lattice.
volume() / mass ) *
1048 #ifdef DEBUG_PRESSURE
1049 iout <<
iINFO <<
"integrating half step, strain rate: " << strainRate <<
"\n";
1058 factor.
xx = exp( dt_long * strainRate.xx );
1059 factor.
yy = exp( dt_long * strainRate.yy );
1061 factor.
xx = factor.
yy = 1;
1063 factor.
zz = exp( dt_long * strainRate.zz );
1066 #ifdef DEBUG_PRESSURE
1067 iout <<
iINFO <<
"rescaling by: " << factor <<
"\n";
1077 factor.
xx = exp( dt_long * strainRate.xx );
1078 factor.
yy = exp( dt_long * strainRate.yy );
1080 factor.
xx = factor.
yy = 1;
1082 factor.
zz = exp( dt_long * strainRate.zz );
1088 #ifdef DEBUG_PRESSURE
1100 factor.
xx = factor.
yy = 1;
1112 #ifdef DEBUG_PRESSURE
1113 iout <<
iINFO <<
"correcting strain rate for nbond, ";
1115 strainRate -= ( 0.5 * dt * cellDims *
state->lattice.
volume() / mass ) *
1117 #ifdef DEBUG_PRESSURE
1118 iout <<
"strain rate: " << strainRate <<
"\n";
1123 #ifdef DEBUG_PRESSURE
1124 iout <<
iINFO <<
"correcting strain rate for slow, ";
1126 strainRate -= ( 0.5 * dt * cellDims *
state->lattice.
volume() / mass ) *
1128 #ifdef DEBUG_PRESSURE
1129 iout <<
"strain rate: " << strainRate <<
"\n";
1157 #ifdef DEBUG_PRESSURE
1158 iout <<
iINFO <<
"correcting strain rate for nbond, ";
1160 strainRate += ( 0.5 * dt * cellDims *
state->lattice.
volume() / mass ) *
1162 #ifdef DEBUG_PRESSURE
1163 iout <<
"strain rate: " << strainRate <<
"\n";
1168 #ifdef DEBUG_PRESSURE
1169 iout <<
iINFO <<
"correcting strain rate for slow, ";
1171 strainRate += ( 0.5 * dt * cellDims *
state->lattice.
volume() / mass ) *
1173 #ifdef DEBUG_PRESSURE
1174 iout <<
"strain rate: " << strainRate <<
"\n";
1188 strainRate += ( 0.5 * dt * cellDims *
state->lattice.
volume() / mass ) *
1192 #ifdef DEBUG_PRESSURE
1193 iout <<
iINFO <<
"integrating half step, strain rate: " << strainRate <<
"\n";
1196 if ( ! ( step % slowFreq ) )
1199 BigReal f1 = exp( -0.5 * dt_long * gamma );
1200 BigReal f2 = sqrt( ( 1. - f1*f1 ) * kT / mass );
1215 #ifdef DEBUG_PRESSURE
1216 iout <<
iINFO <<
"applying langevin, strain rate: " << strainRate <<
"\n";
1220 #ifdef DEBUG_PRESSURE
1221 iout <<
iINFO <<
"exiting langevinPiston2, strain rate: " << strainRate <<
"\n";
1236 if ( rescaleFreq > 0 ) {
1245 BigReal factor = sqrt(rescaleTemp/avgTemp);
1263 if ( momentum.
length2() == 0. )
1264 iout <<
iERROR <<
"Momentum is exactly zero; probably a bug.\n" <<
endi;
1266 NAMD_die(
"Total mass is zero in Controller::correctMomentum");
1268 momentum *= (-1./mass);
1283 iout <<
"FEP: RESETTING FOR NEW FEP WINDOW "
1284 <<
"LAMBDA SET TO " << alchLambda <<
" LAMBDA2 " << alchLambda2;
1285 if ( alchLambdaIDWS >= 0. ) {
1286 iout <<
" LAMBDA_IDWS " << alchLambdaIDWS;
1288 iout <<
"\nFEP: WINDOW TO HAVE " << alchEquilSteps
1289 <<
" STEPS OF EQUILIBRATION PRIOR TO FEP DATA COLLECTION.\n"
1290 <<
"FEP: USING CONSTANT TEMPERATURE OF " << alchTemp
1291 <<
" K FOR FEP CALCULATION\n" <<
endi;
1300 iout <<
"TI: RESETTING FOR NEW WINDOW "
1301 <<
"LAMBDA SET TO " << alchLambda
1302 <<
"\nTI: WINDOW TO HAVE " << alchEquilSteps
1303 <<
" STEPS OF EQUILIBRATION PRIOR TO TI DATA COLLECTION.\n" <<
endi;
1311 if ( ( reassignFreq > 0 ) && ! ( step % reassignFreq ) ) {
1318 if ( newTemp < simParams->reassignHold )
1321 iout <<
"REASSIGNING VELOCITIES AT STEP " << step
1322 <<
" TO " << newTemp <<
" KELVIN.\n" <<
endi;
1366 double coefficient = 1;
1390 static char tmp_string[25];
1391 const double maxnum = 99999999999.9999;
1392 if ( X > maxnum ) X = maxnum;
1393 if ( X < -maxnum ) X = -maxnum;
1394 sprintf(tmp_string,
" %14.4f",X);
1400 static char tmp_string[25];
1401 sprintf(tmp_string,
" %14s",X);
1407 static char tmp_string[21];
1408 sprintf(tmp_string,
"ENERGY: %7d",X);
1449 int numFixedGroups = ( numFixedAtoms ? molecule->
numFixedGroups : 0 );
1452 if ( numFixedGroups ) numGroupDegFreedom -= 3 * numFixedGroups;
1456 numGroupDegFreedom -= 3;
1463 int numFixedRigidBonds =
1468 numDegFreedom -= ( numRigidBonds - numFixedRigidBonds - numLonepairs);
1483 BigReal groupTempHalfstep = 2.0 * groupKineticEnergyHalfstep
1485 BigReal groupTempCentered = 2.0 * groupKineticEnergyCentered
1558 virial_normal, virial_nbond, virial_slow,
1559 intVirial_normal, intVirial_nbond, intVirial_slow,
1560 extForce_normal, extForce_nbond, extForce_slow);
1562 #ifdef DEBUG_PRESSURE
1591 const Tensor& virial_normal_in,
const Tensor& virial_nbond_in,
const Tensor& virial_slow_in,
1592 const Tensor& intVirial_normal,
const Tensor& intVirial_nbond,
const Tensor& intVirial_slow,
1593 const Vector& extForce_normal,
const Vector& extForce_nbond,
const Vector& extForce_slow) {
1595 Tensor virial_normal = virial_normal_in;
1596 Tensor virial_nbond = virial_nbond_in;
1597 Tensor virial_slow = virial_slow_in;
1611 virial_normal -=
outer(extForce_normal,extPosition);
1612 virial_nbond -=
outer(extForce_nbond,extPosition);
1613 virial_slow -=
outer(extForce_slow,extPosition);
1615 if ( (volume=lattice.
volume()) != 0. )
1619 #ifdef MEM_OPT_VERSION
1620 NAMD_bug(
"LJcorrection not supported in memory optimized build.");
1638 if ( minimize || ! ( step %
nbondFreq ) )
1644 if ( minimize || ! ( step %
slowFreq ) )
1687 #define AVGXY(T) T.xy = T.yx = 0; T.xx = T.yy = 0.5 * ( T.xx + T.yy );\
1688 T.xz = T.zx = T.yz = T.zy = 0.5 * ( T.xz + T.yz );
1715 else if(testV < *Vmin){
1720 delta = testV - *Vavg;
1722 *M2 += delta * (testV - (*Vavg));
1724 *sigmaV = sqrt(*M2/n);
1732 *k0 = (1-sigma0/sigmaV) * (Vmax-Vmin) / (Vavg-Vmin);
1740 *E = Vmin + (Vmax-Vmin)/(*k0);
1746 *k0 = (sigma0 / sigmaV) * (Vmax-Vmin) / (Vmax-Vavg);
1753 *k = *k0 / (Vmax - Vmin);
1763 *vir += vir_orig * (*factor);
1764 *dV += (*factor) * VE/2.;
1773 (
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){
1775 char timestepstr[20];
1779 const char *bsuffix;
1781 if(lasttime ||
simParams->restartFrequency == 0){
1790 baselen = strlen(filename);
1791 restart_name =
new char[baselen+26];
1793 strcpy(restart_name, filename);
1794 if (
simParams->restartSave && !lasttime) {
1795 sprintf(timestepstr,
".%d",step_n);
1796 strcat(restart_name, timestepstr);
1799 strcat(restart_name,
".gamd");
1804 rest = fopen(restart_name,
"w");
1806 fprintf(rest,
"# NAMD accelMDG restart file\n"
1807 "# D/T Vn Vmax Vmin Vavg sigmaV M2 E k\n"
1808 "%c %d %.17lg %.17lg %.17lg %.17lg %.17le %.17lg %.17lg\n",
1809 type, V_n-1, Vmax, Vmin, Vavg, sigmaV, M2, E, k);
1811 iout <<
"GAUSSIAN ACCELERATED MD RESTART DATA WRITTEN TO " << restart_name <<
"\n" <<
endi;
1814 iout <<
iWARN <<
"Cannot write accelMDG restart file\n" <<
endi;
1817 rest = fopen(restart_name,
"a");
1819 fprintf(rest,
"%c %d %.17lg %.17lg %.17lg %.17lg %.17le %.17lg %.17lg\n",
1820 type, V_n-1, Vmax, Vmin, Vavg, sigmaV, M2, E, k);
1824 iout <<
iWARN <<
"Cannot write accelMDG restart file\n" <<
endi;
1827 delete[] restart_name;
1861 static BigReal VmaxP, VminP, VavgP, M2P=0, sigmaVP;
1862 static BigReal VmaxD, VminD, VavgD, M2D=0, sigmaVD;
1865 static int V_n = 1, iEusedD, iEusedP;
1866 static char warnD, warnP;
1867 static int restfreq;
1886 BigReal amd_ljEnergy_Corr = 0.;
1903 volume = lattice.
volume();
1947 potentialEnergy = bondEnergy + angleEnergy + dihedralEnergy +
1948 improperEnergy + amd_electEnergy + amd_electEnergySlow + amd_ljEnergy +
1949 crosstermEnergy + boundaryEnergy + miscEnergy +
1950 amd_goTotalEnergy + amd_groLJEnergy + amd_groGaussEnergy;
1955 if(step == firststep) {
1957 warnD = warnP =
'\0';
1963 int dihe_n=0, tot_n=0;
1969 while(fgets(line, 256, rest)){
1971 dihe_n = sscanf(line+1,
" %d %la %la %la %la %la %la %la",
1972 &V_n, &VmaxD, &VminD, &VavgD, &sigmaVD, &M2D, &ED, &kD);
1974 else if(line[0] ==
'T'){
1975 tot_n = sscanf(line+1,
" %d %la %la %la %la %la %la %la",
1976 &V_n, &VmaxP, &VminP, &VavgP, &sigmaVP, &M2P, &EP, &kP);
1987 NAMD_die(
"Invalid dihedral potential energy format in accelMDG restart file");
1988 k0D = kD * (VmaxD - VminD);
1989 iout <<
"GAUSSIAN ACCELERATED MD READ FROM RESTART FILE: DIHED"
1990 <<
" Vmax " << VmaxD <<
" Vmin " << VminD
1991 <<
" Vavg " << VavgD <<
" sigmaV " << sigmaVD
1992 <<
" E " << ED <<
" k " << kD
1999 NAMD_die(
"Invalid total potential energy format in accelMDG restart file");
2000 k0P = kP * (VmaxP - VminP);
2001 iout <<
"GAUSSIAN ACCELERATED MD READ FROM RESTART FILE: TOTAL"
2002 <<
" Vmax " << VmaxP <<
" Vmin " << VminP
2003 <<
" Vavg " << VavgP <<
" sigmaV " << sigmaVP
2004 <<
" E " << EP <<
" k " << kP
2008 iEusedD = (ED == VmaxD) ? 1 : 2;
2009 iEusedP = (EP == VmaxP) ? 1 : 2;
2017 testV = dihedralEnergy + crosstermEnergy;
2020 if(step > firststep && step % restfreq == 0)
2021 write_accelMDG_rest_file(step,
'D', V_n, VmaxD, VminD, VavgD, sigmaVD, M2D, ED, kD,
2026 write_accelMDG_rest_file(step,
'D', V_n, VmaxD, VminD, VavgD, sigmaVD, M2D, ED, kD,
2030 write_accelMDG_rest_file(step,
'D', V_n, VmaxD, VminD, VavgD, sigmaVD, M2D, ED, kD,
2038 VmaxD = VminD = VavgD = testV;
2042 else if(step == ntcmdprep && ntcmdprep != 0){
2044 VmaxD = VminD = VavgD = testV;
2046 iout <<
"GAUSSIAN ACCELERATED MD: RESET DIHEDRAL POTENTIAL STATISTICS\n" <<
endi;
2049 else if(ntave > 0 && step % ntave == 0){
2051 if(testV > VmaxD) VmaxD = testV;
2052 if(testV < VminD) VminD = testV;
2056 iout <<
"GAUSSIAN ACCELERATED MD: RESET DIHEDRAL POTENTIAL AVG AND STD\n" <<
endi;
2061 &VmaxD, &VminD, &VavgD, &M2D, &sigmaVD) ;
2064 if(step == ntcmd - 1){
2066 VmaxD, VminD, VavgD, sigmaVD, &k0D, &kD, &ED, &iEusedD, &warnD);
2070 else if(step < nteb){
2072 &dV, &factor_dihe, &vir);
2077 VmaxD = VminD = VavgD = testV;
2079 iout <<
"GAUSSIAN ACCELERATED MD: RESET DIHEDRAL POTENTIAL STATISTICS\n" <<
endi;
2082 else if(ntave > 0 && step % ntave == 0){
2084 if(testV > VmaxD) VmaxD = testV;
2085 if(testV < VminD) VminD = testV;
2089 iout <<
"GAUSSIAN ACCELERATED MD: RESET DIHEDRAL POTENTIAL AVG AND STD\n" <<
endi;
2093 &VmaxD, &VminD, &VavgD, &M2D, &sigmaVD) ;
2096 if(step >= ntebprep)
2097 if((ntave > 0 && step % ntave == 0) || ntave <= 0)
2099 VmaxD, VminD, VavgD, sigmaVD, &k0D, &kD, &ED, &iEusedD, &warnD);
2104 &dV, &factor_dihe, &vir);
2110 Tensor vir_tot = vir_normal + vir_nbond + vir_slow;
2111 testV = potentialEnergy;
2113 testV -= dihedralEnergy + crosstermEnergy;
2114 vir_tot -= vir_dihe;
2118 if(step > firststep && step % restfreq == 0)
2119 write_accelMDG_rest_file(step,
'T', V_n, VmaxP, VminP, VavgP, sigmaVP, M2P, EP, kP,
2124 write_accelMDG_rest_file(step,
'T', V_n, VmaxP, VminP, VavgP, sigmaVP, M2P, EP, kP,
2128 write_accelMDG_rest_file(step,
'T', V_n, VmaxP, VminP, VavgP, sigmaVP, M2P, EP, kP,
2136 VmaxP = VminP = VavgP = testV;
2140 else if(step == ntcmdprep && ntcmdprep != 0){
2142 VmaxP = VminP = VavgP = testV;
2144 iout <<
"GAUSSIAN ACCELERATED MD: RESET TOTAL POTENTIAL STATISTICS\n" <<
endi;
2147 else if(ntave > 0 && step % ntave == 0){
2149 if(testV > VmaxP) VmaxP = testV;
2150 if(testV < VminP) VminP = testV;
2154 iout <<
"GAUSSIAN ACCELERATED MD: RESET TOTAL POTENTIAL AVG AND STD\n" <<
endi;
2159 &VmaxP, &VminP, &VavgP, &M2P, &sigmaVP);
2161 if(step == ntcmd - 1){
2163 VmaxP, VminP, VavgP, sigmaVP, &k0P, &kP, &EP, &iEusedP, &warnP);
2167 else if(step < nteb){
2169 &dV, &factor_tot, &vir);
2174 VmaxP = VminP = VavgP = testV;
2176 iout <<
"GAUSSIAN ACCELERATED MD: RESET TOTAL POTENTIAL STATISTICS\n" <<
endi;
2179 else if(ntave > 0 && step % ntave == 0){
2181 if(testV > VmaxP) VmaxP = testV;
2182 if(testV < VminP) VminP = testV;
2186 iout <<
"GAUSSIAN ACCELERATED MD: RESET TOTAL POTENTIAL AVG AND STD\n" <<
endi;
2190 &VmaxP, &VminP, &VavgP, &M2P, &sigmaVP);
2193 if(step >= ntebprep)
2194 if((ntave > 0 && step % ntave == 0) || ntave <= 0)
2196 VmaxP, VminP, VavgP, sigmaVP, &k0P, &kP, &EP, &iEusedP, &warnP);
2201 &dV, &factor_tot, &vir);
2208 if((ntcmdprep > 0 && step == ntcmdprep) ||
2209 (step < nteb && ntave > 0 && step % ntave == 0) ||
2222 testV = dihedralEnergy + crosstermEnergy;
2223 if ( testV < accelMDE ) {
2224 factor_dihe = accelMDalpha/(accelMDalpha + accelMDE - testV);
2225 factor_dihe *= factor_dihe;
2226 vir = vir_dihe * (factor_dihe - 1.0);
2227 dV = (accelMDE - testV)*(accelMDE - testV)/(accelMDalpha + accelMDE - testV);
2233 testV = dihedralEnergy + crosstermEnergy;
2234 if ( testV < accelMDE ) {
2235 factor_dihe = accelMDalpha/(accelMDalpha + accelMDE - testV);
2236 factor_dihe *= factor_dihe;
2237 vir = vir_dihe * (factor_dihe - 1.0) ;
2238 dV = (accelMDE - testV)*(accelMDE - testV)/(accelMDalpha + accelMDE - testV);
2241 testV = potentialEnergy - dihedralEnergy - crosstermEnergy;
2242 if ( testV < accelMDTE ) {
2243 factor_tot = accelMDTalpha/(accelMDTalpha + accelMDTE - testV);
2244 factor_tot *= factor_tot;
2245 vir += (vir_normal - vir_dihe + vir_nbond + vir_slow) * (factor_tot - 1.0);
2246 dV += (accelMDTE - testV)*(accelMDTE - testV)/(accelMDTalpha + accelMDTE - testV);
2252 testV = potentialEnergy;
2253 if ( testV < accelMDE ) {
2254 factor_tot = accelMDalpha/(accelMDalpha + accelMDE - testV);
2255 factor_tot *= factor_tot;
2256 vir = (vir_normal + vir_nbond + vir_slow) * (factor_tot - 1.0);
2257 dV = (accelMDE - testV)*(accelMDE - testV)/(accelMDalpha + accelMDE - testV);
2263 accelMDfactor[0]=factor_dihe;
2264 accelMDfactor[1]=factor_tot;
2270 if ( factor_tot < 0.001 ) {
2271 iout <<
iWARN <<
"accelMD is using a very high boost potential, simulation may become unstable!"
2275 if ( ! (step % accelMDOutFreq) ) {
2279 iout <<
"ACCELERATED MD: STEP " << step
2282 <<
" BOND " << bondEnergy
2283 <<
" ANGLE " << angleEnergy
2284 <<
" DIHED " << dihedralEnergy+crosstermEnergy
2285 <<
" IMPRP " << improperEnergy
2286 <<
" ELECT " << amd_electEnergy+amd_electEnergySlow
2287 <<
" VDW " << amd_ljEnergy
2288 <<
" POTENTIAL " << potentialEnergy;
2289 if(amd_ljEnergy_Corr)
2290 iout <<
" LJcorr " << amd_ljEnergy_Corr;
2294 iout <<
"GAUSSIAN ACCELERATED MD: DIHED iE " << iEusedD
2295 <<
" Vmax " << VmaxD <<
" Vmin " << VminD
2296 <<
" Vavg " << VavgD <<
" sigmaV " << sigmaVD
2297 <<
" E " << ED <<
" k0 " << k0D <<
" k " << kD
2300 iout <<
"GAUSSIAN ACCELERATED MD: DIHED !!! WARNING: k0 > 1, "
2301 <<
"SWITCHED TO iE = 1 MODE !!!\n" <<
endi;
2303 iout <<
"GAUSSIAN ACCELERATED MD: DIHED !!! WARNING: k0 <= 0, "
2304 <<
"SWITCHED TO iE = 1 MODE !!!\n" <<
endi;
2306 iout <<
"GAUSSIAN ACCELERATED MD: TOTAL iE " << iEusedP
2307 <<
" Vmax " << VmaxP <<
" Vmin " << VminP
2308 <<
" Vavg " << VavgP <<
" sigmaV " << sigmaVP
2309 <<
" E " << EP <<
" k0 " << k0P <<
" k " << kP
2312 iout <<
"GAUSSIAN ACCELERATED MD: TOTAL !!! WARNING: k0 > 1, "
2313 <<
"SWITCHED TO iE = 1 MODE !!!\n" <<
endi;
2315 iout <<
"GAUSSIAN ACCELERATED MD: TOTAL !!! WARNING: k0 <= 0, "
2316 <<
"SWITCHED TO iE = 1 MODE !!!\n" <<
endi;
2317 warnD = warnP =
'\0';
2327 if ( volume != 0. ) {
2328 p_normal = vir_normal/volume;
2329 p_nbond = vir_nbond/volume;
2330 p_slow = vir_slow/volume;
2333 iout <<
" accelMD Scaling Factor: " << accelMDfactor <<
"\n"
2334 <<
" accelMD PNORMAL " << trace(p_normal)*
PRESSUREFACTOR/3. <<
"\n"
2335 <<
" accelMD PNBOND " << trace(p_nbond)*
PRESSUREFACTOR/3. <<
"\n"
2345 iout <<
iINFO <<
"INITIALISING ADAPTIVE TEMPERING\n" <<
endi;
2350 iout <<
iINFO <<
"READING ADAPTIVE TEMPERING RESTART FILE\n" <<
endi;
2352 if (adaptTempRead) {
2357 adaptTempRead >> readInt;
2365 if (adaptTempBetaMin > adaptTempBetaMax){
2390 adaptTempRead.close();
2392 else NAMD_die(
"Could not open ADAPTIVE TEMPERING restart file.\n");
2432 iout <<
iINFO <<
"ADAPTIVE TEMPERING: ADAPTIVE BIN AVERAGING PARAMETER: " <<
adaptTempCg <<
"\n";
2433 iout <<
iINFO <<
"ADAPTIVE TEMPERING: SCALING CONSTANT FOR LANGEVIN TEMPERATURE UPDATES: " <<
adaptTempDt <<
"\n";
2439 NAMD_die(
"Error opening restart file");
2446 iout <<
"ADAPTEMP: WRITING RESTART FILE AT STEP " << step <<
"\n" <<
endi;
2478 if ( minimize || (step < simParams->adaptTempFirstStep ) ||
2488 <<
" adaptTempBeta: " << adaptTempBeta
2507 potEnergyVariance -= potEnergyAverage*potEnergyAverage;
2520 BigReal deltaBeta = 0.04*adaptTempBeta;
2526 if ( adaptTempBetaMax-adaptTempBeta < deltaBeta ) deltaBeta = adaptTempBetaMax-adaptTempBeta;
2527 betaMinus = adaptTempBeta - deltaBeta;
2528 betaPlus = adaptTempBeta + deltaBeta;
2529 BigReal betaMinus2 = betaMinus*betaMinus;
2530 BigReal betaPlus2 = betaPlus*betaPlus;
2531 nMinus = (int)floor((betaMinus-adaptTempBetaMin)/
adaptTempDBeta);
2548 DeltaE2AveSum += DeltaE2Ave;
2552 A0 += (adaptTempBetaMin+0.5*
adaptTempDBeta-betaMinus)*DeltaE2AveSum;
2557 A0 /= (betaNp1-betaMinus);
2563 DeltaE2Ave = adaptTempPotEnergyVar[j];
2564 DeltaE2AveSum += DeltaE2Ave;
2568 A1 += (adaptTempBetaMin+0.5*
adaptTempDBeta-betaPlus)*DeltaE2AveSum;
2570 if ( nPlus < adaptTempBins ) {
2573 A1 -= 0.5*betaj*adaptTempPotEnergyVar[nPlus];
2575 A1 /= (betaPlus-betaNp1);
2584 if (DeltaE2Ave != 0) {
2585 aplus = (A0-A2)/(A0-A1);
2595 potEnergyAve0 /= (betaNp1-betaMinus);
2597 if ( nPlus < adaptTempBins ) {
2600 potEnergyAve1 /= (betaPlus-betaNp1);
2601 potEnergyAverage = aminus*potEnergyAve0;
2602 potEnergyAverage += aplus*potEnergyAve1;
2605 iout <<
"ADAPTEMP DEBUG:" <<
"\n"
2608 <<
" adaptTempBeta: " << adaptTempBeta <<
"\n"
2610 <<
" betaMin: " << adaptTempBetaMin <<
"\n"
2611 <<
" betaMax: " << adaptTempBetaMax <<
"\n"
2612 <<
" gammaAve: " << gammaAve <<
"\n"
2613 <<
" deltaBeta: " << deltaBeta <<
"\n"
2614 <<
" betaMinus: " << betaMinus <<
"\n"
2615 <<
" betaPlus: " << betaPlus <<
"\n"
2616 <<
" nMinus: " << nMinus <<
"\n"
2617 <<
" nPlus: " << nPlus <<
"\n"
2618 <<
" A0: " << A0 <<
"\n"
2619 <<
" A1: " << A1 <<
"\n"
2620 <<
" A2: " << A2 <<
"\n"
2621 <<
" a+: " << aplus <<
"\n"
2622 <<
" a-: " << aminus <<
"\n"
2628 iout <<
"ADAPTEMP DEBUG:" <<
"\n"
2631 <<
" adaptTempBeta: " << adaptTempBeta <<
"\n"
2634 <<
" betaMax: " << adaptTempBetaMax <<
"\n"
2635 <<
" gammaAve: " << gammaAve <<
"\n"
2636 <<
" deltaBeta: N/A\n"
2637 <<
" betaMinus: N/A\n"
2638 <<
" betaPlus: N/A\n"
2678 else if ( dT < 1./adaptTempBetaMax ) {
2697 if (adaptTempTdiff > 0) {
2709 iout <<
"ADAPTEMP: " << step <<
" FRAC " << Frac <<
"\n";
2712 iout <<
"ADAPTEMP: Updating adaptTempDt to ";
2725 if (adaptTempTdiff > 0) {
2737 iout <<
"ADAPTEMP: " << step <<
" FRAC " << Frac <<
"\n";
2740 iout <<
"ADAPTEMP: Updating adaptTempDt to ";
2755 if ( ! (step % adaptTempOutFreq) ) {
2756 iout <<
"ADAPTEMP: STEP " << step
2758 <<
" ENERGY " << std::setprecision(10) << potentialEnergy
2759 <<
" ENERGYAVG " << std::setprecision(10) << potEnergyAverage
2760 <<
" ENERGYVAR " << std::setprecision(10) << potEnergyVariance;
2776 if ( ((
int)checksum) != molecule->
numAtoms ) {
2777 sprintf(errmsg,
"Bad global atom count! (%d vs %d)\n",
2778 (
int)checksum, molecule->
numAtoms);
2785 computesPartitioned = 0;
2786 }
else if ( (
int)checksum >
computeChecksum && ! computesPartitioned ) {
2787 computesPartitioned = 1;
2789 NAMD_bug(
"Bad global compute count!\n");
2795 if ( checksum_b && (((
int)checksum) != molecule->
numCalcBonds) ) {
2796 sprintf(errmsg,
"Bad global bond count! (%d vs %d)\n",
2798 if ( forgiving && (((
int)checksum) < molecule->
numCalcBonds) )
2804 if ( checksum_b && (((
int)checksum) != molecule->
numCalcAngles) ) {
2805 sprintf(errmsg,
"Bad global angle count! (%d vs %d)\n",
2807 if ( forgiving && (((
int)checksum) < molecule->
numCalcAngles) )
2814 sprintf(errmsg,
"Bad global dihedral count! (%d vs %d)\n",
2823 sprintf(errmsg,
"Bad global improper count! (%d vs %d)\n",
2831 if ( checksum_b && (((
int)checksum) != molecule->
numCalcTholes) ) {
2832 sprintf(errmsg,
"Bad global Thole count! (%d vs %d)\n",
2834 if ( forgiving && (((
int)checksum) < molecule->
numCalcTholes) )
2840 if ( checksum_b && (((
int)checksum) != molecule->
numCalcAnisos) ) {
2841 sprintf(errmsg,
"Bad global Aniso count! (%d vs %d)\n",
2843 if ( forgiving && (((
int)checksum) < molecule->
numCalcAnisos) )
2850 sprintf(errmsg,
"Bad global crossterm count! (%d vs %d)\n",
2861 iout <<
iWARN <<
"High global exclusion count ("
2862 << ((
int64)checksum) <<
" vs "
2864 <<
iWARN <<
"This warning is not unusual during minimization.\n"
2865 <<
iWARN <<
"Decreasing pairlistdist or cutoff that is too close to periodic cell size may avoid this.\n" <<
endi;
2868 sprintf(errmsg,
"High global exclusion count! (%lld vs %lld) System unstable or pairlistdist or cutoff too close to periodic cell size.\n",
2877 iout <<
iWARN <<
"Low global exclusion count! ("
2881 <<
iWARN <<
"This warning is not unusual during minimization.\n"
2882 <<
iWARN <<
"Increasing pairlistdist or cutoff may avoid this.\n";
2884 iout <<
" System unstable or pairlistdist or cutoff too small.\n";
2889 sprintf(errmsg,
"Low global exclusion count! (%lld vs %lld) System unstable or pairlistdist or cutoff too small.\n",
2900 iout <<
iWARN <<
"High global CUDA exclusion count ("
2901 << ((
int64)checksum) <<
" vs "
2903 <<
iWARN <<
"This warning is not unusual during minimization.\n"
2904 <<
iWARN <<
"Decreasing pairlistdist or cutoff that is too close to periodic cell size may avoid this.\n" <<
endi;
2907 sprintf(errmsg,
"High global CUDA exclusion count! (%lld vs %lld) System unstable or pairlistdist or cutoff too close to periodic cell size.\n",
2916 iout <<
iWARN <<
"Low global CUDA exclusion count! ("
2920 <<
iWARN <<
"This warning is not unusual during minimization.\n"
2921 <<
iWARN <<
"Increasing pairlistdist or cutoff may avoid this.\n";
2923 iout <<
" System unstable or pairlistdist or cutoff too small.\n";
2928 sprintf(errmsg,
"Low global CUDA exclusion count! (%lld vs %lld) System unstable or pairlistdist or cutoff too small.\n",
2937 iout <<
iERROR <<
"Margin is too small for " << ((int)checksum) <<
2938 " atoms during timestep " << step <<
".\n" <<
iERROR <<
2939 "Incorrect nonbonded forces and energies may be calculated!\n" <<
endi;
2946 "Pairlistdist is too small for " << ((int)checksum) <<
2947 " computes during timestep " << step <<
".\n" <<
endi;
2954 iout <<
iWARN <<
"Stray PME grid charges ignored!\n" <<
endi;
2955 else NAMD_bug(
"Stray PME grid charges detected!\n");
2963 const double endWTime = CmiWallTimer() - firstWTime;
2964 const double endCTime = CmiTimer() - firstCTime;
2967 if ( (((
int)endWTime)>>6) != (((
int)startWTime)>>6 ) )
fflush_count = 2;
2969 const double elapsedW =
2971 const double elapsedC =
2974 const double remainingW = elapsedW * (
simParams->
N - step);
2975 const double remainingW_hours = remainingW / 3600;
2977 startWTime = endWTime;
2978 startCTime = endCTime;
2983 iout <<
iWARN <<
"Energy evaluation is expensive, increase outputEnergies to improve performance.\n" <<
endi;
2987 CmiPrintf(
"TIMING: %d CPU: %g, %g/step Wall: %g, %g/step"
2988 ", %g hours remaining, %f MB of memory in use.\n",
2989 step, endCTime, elapsedC, endWTime, elapsedW,
3061 if ( minimize || ! ( step %
nbondFreq ) )
3088 if ( minimize || ! ( step %
slowFreq ) )
3102 #ifdef MEM_OPT_VERSION
3103 NAMD_bug(
"LJcorrection not supported in memory optimized build.");
3109 if (simParameters->
alchOn) {
3136 angularMomentum.
x =
reduction->
item(REDUCTION_ANGULAR_MOMENTUM_X);
3137 angularMomentum.
y =
reduction->
item(REDUCTION_ANGULAR_MOMENTUM_Y);
3138 angularMomentum.
z =
reduction->
item(REDUCTION_ANGULAR_MOMENTUM_Z);
3141 potentialEnergy = (bondEnergy + angleEnergy + dihedralEnergy
3143 + crosstermEnergy + boundaryEnergy + miscEnergy +
goTotalEnergy
3149 if ( !(step%slowFreq) ) {
3162 if (step <= simParameters->firstTimestep &&
3170 iout <<
"MOMENTUM: " << step
3171 <<
" P: " << momentum
3172 <<
" L: " << angularMomentum
3181 iout <<
"PRESSURE: " << step <<
" "
3183 <<
"GPRESSURE: " << step <<
" "
3187 <<
"GPRESSAVG: " << step <<
" "
3191 groupPressure_tavg = 0;
3202 memset(total, 0, arraysize*
sizeof(
BigReal));
3211 if (!(step % freq)) {
3215 iout <<
"PRESSUREPROFILE: " << step <<
" ";
3216 if (step == first) {
3218 for (
int i=0; i<arraysize; i++) {
3219 iout << total[i] * scalefac <<
" ";
3224 for (
int i=0; i<arraysize; i++)
3239 if ( benchPhase > 0 && benchPhase < 10 ) {
3242 iout <<
iWARN <<
"Energy evaluation is expensive, increase outputEnergies to improve performance.\n";
3246 if ( benchPhase < 4 )
iout <<
"Initial time: ";
3247 else iout <<
"Benchmark time: ";
3248 iout << CkNumPes() <<
" CPUs ";
3253 BigReal days = 1.0 / (24.0 * 60.0 * 60.0);
3254 BigReal daysPerNano = wallPerStep * days / ns;
3255 iout << wallPerStep <<
" s/step " << daysPerNano <<
" days/ns ";
3259 startBenchTime = CmiWallTimer();
3266 Vector pairVDWForce, pairElectForce;
3279 #define CALLBACKDATA(LABEL,VALUE) \
3280 labels << (LABEL) << " "; values << (VALUE) << " ";
3281 #define CALLBACKLIST(LABEL,VALUE) \
3282 labels << (LABEL) << " "; values << "{" << (VALUE) << "} ";
3284 std::ostringstream labels,
values;
3287 values.precision(16);
3289 values << std::setprecision(16);
3312 labels <<
"PERIODIC "; values <<
"{" << lattice.
a_p() <<
" "
3313 << lattice.
b_p() <<
" " << lattice.
c_p() <<
"} ";
3325 if (simParameters->
alchOn) {
3339 CALLBACKLIST(
"BOND2", bondEnergy + angleEnergy + dihedralEnergy +
3346 labels <<
'\0'; values <<
'\0';
3347 state->callback_labelstring = labels.str();
3348 state->callback_valuestring = values.str();
3365 " steps.\n" <<
endi;
3375 if ( ! minimize && step % simParameters->
outputEnergies )
return;
3378 if (simParameters->
IMDon && !(step % simParameters->
IMDfreq)) {
3380 energies.
tstep = step;
3383 energies.
Epot = potentialEnergy;
3386 energies.
Ebond = bondEnergy;
3387 energies.
Eangle = angleEnergy;
3388 energies.
Edihe = dihedralEnergy + crosstermEnergy;
3389 energies.
Eimpr = improperEnergy;
3395 " margin violations detected since previous energy output.\n" <<
endi;
3399 if ( (step % (10 * (minimize?1:simParameters->
outputEnergies) ) ) == 0 )
3401 iout <<
"ETITLE: TS";
3420 if ( volume != 0. ) {
3456 if (simParameters->
alchOn) {
3458 iout <<
"\nTITITLE: TS";
3472 iout <<
"\nFEPTITLE: TS";
3485 iout <<
"QMETITLE: TS";
3500 iout <<
FORMAT(dihedralEnergy+crosstermEnergy);
3555 if (simParameters->
alchOn) {
3576 iout <<
FORMAT(bondEnergy + angleEnergy + dihedralEnergy
3588 #if(CMK_CCS_AVAILABLE && CMK_WEB_MODE)
3591 (
int)(potentialEnergy),
3593 CApplicationDepositNode0Data(webout);
3597 iout <<
"PAIR INTERACTION:";
3598 iout <<
" STEP: " << step;
3599 iout <<
" VDW_FORCE: ";
3603 iout <<
" ELECT_FORCE: ";
3623 file <<
"#$LABELS step";
3624 if ( lattice.
a_p() ) file <<
" a_x a_y a_z";
3625 if ( lattice.
b_p() ) file <<
" b_x b_y b_z";
3626 if ( lattice.
c_p() ) file <<
" c_x c_y c_z";
3627 file <<
" o_x o_y o_z";
3629 file <<
" s_x s_y s_z s_u s_v s_w";
3638 if ( lattice.
a_p() ) file <<
" " << lattice.
a().
x <<
" " << lattice.
a().
y <<
" " << lattice.
a().
z;
3639 if ( lattice.
b_p() ) file <<
" " << lattice.
b().
x <<
" " << lattice.
b().
y <<
" " << lattice.
b().
z;
3640 if ( lattice.
c_p() ) file <<
" " << lattice.
c().
x <<
" " << lattice.
c().
y <<
" " << lattice.
c().
z;
3645 file <<
" " << strainRate.
x;
3646 file <<
" " << strainRate.
y;
3647 file <<
" " << strainRate.
z;
3648 file <<
" " << strainRate2.
x;
3649 file <<
" " << strainRate2.
y;
3650 file <<
" " << strainRate2.
z;
3674 if (alchEnsembleAvg && (stepInRun == 0 || stepInRun == alchEquilSteps)) {
3692 if (stepInRun == 0) {
3696 iout <<
"OPENING FEP ENERGY OUTPUT FILE\n" <<
endi;
3697 if(alchEnsembleAvg){
3700 <<
"vdW dE dE_avg Temp dG\n"
3702 <<
" l l+dl E(l+dl)-E(l)" << std::endl;
3708 <<
" l l+dl E(l+dl)-E(l)" << std::endl;
3712 fepFile <<
"#NEW FEP WINDOW: "
3713 <<
"LAMBDA SET TO " << alchLambda <<
" LAMBDA2 "
3721 if ((alchEquilSteps) && (stepInRun == alchEquilSteps)) {
3722 fepFile <<
"#" << alchEquilSteps <<
" STEPS OF EQUILIBRATION AT "
3724 <<
"#STARTING COLLECTION OF ENSEMBLE AVERAGE" << std::endl;
3730 if (alchEnsembleAvg && (step ==
simParams->
N)) {
3732 fepFile <<
"#Free energy change for lambda window [ " << alchLambda
3734 <<
" ; net change until now is " <<
fepSum << std::endl;
3747 if (stepInRun == 0 || stepInRun == alchEquilSteps) {
3769 if (stepInRun == 0 || stepInRun == alchEquilSteps
3793 if (stepInRun == 0) {
3803 iout <<
"OPENING TI ENERGY OUTPUT FILE\n" <<
endi;
3804 tiFile <<
"#TITITLE: TS";
3819 if (alchLambdaFreq > 0) {
3826 if (alchLambdaFreq > 0) {
3827 tiFile <<
"#ALCHEMICAL SWITCHING ACTIVE "
3829 <<
"\n#LAMBDA SCHEDULE: "
3831 <<
" Freq: " << alchLambdaFreq;
3841 tiFile <<
"#NEW TI WINDOW: "
3842 <<
"LAMBDA " << alchLambda
3843 <<
"\n#PARTITION 1 SCALING: BOND " << bond_lambda_1
3844 <<
" VDW " << vdw_lambda_1 <<
" ELEC " << elec_lambda_1
3845 <<
"\n#PARTITION 2 SCALING: BOND " << bond_lambda_2
3846 <<
" VDW " << vdw_lambda_2 <<
" ELEC " << elec_lambda_2;
3852 if ((alchEquilSteps) && (stepInRun == alchEquilSteps)) {
3853 tiFile <<
"#" << alchEquilSteps <<
" STEPS OF EQUILIBRATION AT "
3855 <<
"#STARTING COLLECTION OF ENSEMBLE AVERAGE" << std::endl;
3889 (elec_lambda_12 - elec_lambda_1)*
3893 (elec_lambda_22 - elec_lambda_2)*
3922 if(alchEnsembleAvg){
3927 if(alchEnsembleAvg){
3970 iout <<
"OPENING EXTENDED SYSTEM TRAJECTORY FILE\n" <<
endi;
3974 if ( errno == EINTR ) {
3975 CkPrintf(
"Warning: Interrupted system call opening XST trajectory file, retrying.\n");
3984 xstFile <<
"# NAMD extended system trajectory file" << std::endl;
3997 iout <<
"WRITING EXTENDED SYSTEM TO RESTART FILE AT STEP "
3998 << step <<
"\n" <<
endi;
4000 const char *bsuffix =
".old";
4003 char timestepstr[20];
4004 sprintf(timestepstr,
".%d",step);
4005 strcat(fname, timestepstr);
4008 strcat(fname,
".xsc");
4012 if ( errno == EINTR ) {
4013 CkPrintf(
"Warning: Interrupted system call opening XSC restart file, retrying.\n");
4015 xscFile.
open(fname);
4019 sprintf(err_msg,
"Error opening XSC restart file %s",fname);
4022 xscFile <<
"# NAMD extended system configuration restart file" << std::endl;
4027 sprintf(err_msg,
"Error writing XSC restart file %s",fname);
4039 iout <<
"WRITING EXTENDED SYSTEM TO OUTPUT FILE AT STEP "
4040 << realstep <<
"\n" <<
endi;
4041 static char fname[140];
4043 strcat(fname,
".xsc");
4047 if ( errno == EINTR ) {
4048 CkPrintf(
"Warning: Interrupted system call opening XSC output file, retrying.\n");
4050 xscFile.
open(fname);
4054 sprintf(err_msg,
"Error opening XSC output file %s",fname);
4057 xscFile <<
"# NAMD extended system configuration output file" << std::endl;
4062 sprintf(err_msg,
"Error writing XSC output file %s",fname);
4071 iout <<
"CLOSING EXTENDED SYSTEM TRAJECTORY FILE\n" <<
endi;
4088 NAMD_die(
"Unable to load checkpoint, requested key was never stored.");
4094 NAMD_die(
"Unable to swap checkpoint, requested key was never stored.");
4100 NAMD_die(
"Unable to free checkpoint, requested key was never stored.");
4113 CkpvAccess(_qd)->process();
4123 startBenchTime -= CmiWallTimer();
4126 startBenchTime += CmiWallTimer();
4134 broadcast->cycleBarrier.publish(step,1);
4135 CkPrintf(
"Cycle time at sync Wall: %f CPU %f\n",
4136 CmiWallTimer()-firstWTime,CmiTimer()-firstCTime);
4142 CkPrintf(
"Cycle time at trace sync (begin) Wall at step %d: %f CPU %f\n", step, CmiWallTimer()-firstWTime,CmiTimer()-firstCTime);
4143 CProxy_Node nd(CkpvAccess(BOCclass_group).node);
4144 nd.traceBarrier(turnOnTrace, step);
4150 CkPrintf(
"Cycle time at trace sync (end) Wall at step %d: %f CPU %f\n", step, CmiWallTimer()-firstWTime,CmiTimer()-firstCTime);
4154 #ifdef MEASURE_NAMD_WITH_PAPI
4155 void Controller::papiMeasureBarrier(
int turnOnMeasure,
int step){
4156 CkPrintf(
"Cycle time at PAPI measurement sync (begin) Wall at step %d: %f CPU %f\n", step, CmiWallTimer()-firstWTime,CmiTimer()-firstCTime);
4157 CProxy_Node nd(CkpvAccess(BOCclass_group).node);
4158 nd.papiMeasureBarrier(turnOnMeasure, step);
4162 void Controller::resumeAfterPapiMeasureBarrier(
int step){
4163 broadcast->papiMeasureBarrier.publish(step,1);
4164 CkPrintf(
"Cycle time at PAPI measurement sync (end) Wall at step %d: %f CPU %f\n", step, CmiWallTimer()-firstWTime,CmiTimer()-firstCTime);
Vector gaussian_vector(void)
Bool accelMDGresetVaftercmd
BigReal berendsenPressureCompressibility
void enqueueCollections(int)
Tensor controlPressure_slow
char scriptStringArg1[128]
BigReal berendsenPressureRelaxationTime
std::ostream & iINFO(std::ostream &s)
#define CALLBACKDATA(LABEL, VALUE)
void recvCheckpointReq(const char *key, int task, checkpoint &cp)
void cycleBarrier(int, int)
void rescaleVelocities(int)
void rescaleaccelMD(int step, int minimize=0)
static int velocityNeeded(int)
void write_accelMDG_rest_file(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)
BigReal getLambdaDelta(void)
void compareChecksums(int, int=0)
Tensor groupPressure_nbond
void NAMD_err(const char *err_msg)
std::map< std::string, checkpoint * > checkpoints
BigReal langevinPistonTemp
SimpleBroadcastObject< int > traceBarrier
void calcPressure(int step, int minimize, const Tensor &virial_normal_in, const Tensor &virial_nbond_in, const Tensor &virial_slow_in, const Tensor &intVirial_normal, const Tensor &intVirial_nbond, const Tensor &intVirial_slow, const Vector &extForce_normal, const Vector &extForce_nbond, const Vector &extForce_slow)
static Tensor diagonal(const Vector &v1)
static int coordinateNeeded(int)
void adaptTempWriteRestart(int step)
static __global__ void(const patch_pair *patch_pairs, const atom *atoms, const atom_param *atom_params, const int *vdw_types, unsigned int *plist, float4 *tmpforces, float4 *slow_tmpforces, float4 *forces, float4 *slow_forces, float *tmpvirials, float *slow_tmpvirials, float *virials, float *slow_virials, unsigned int *global_counters, int *force_ready_queue, const unsigned int *overflow_exclusions, const int npatches, const int block_begin, const int total_block_count, int *block_order, exclmask *exclmasks, const int lj_table_size, const float3 lata, const float3 latb, const float3 latc, const float cutoff2, const float plcutoff2, const int doSlow)
SimpleBroadcastObject< Vector > momentumCorrection
ScriptTcl * getScript(void)
virtual void algorithm(void)
static PatchMap * Object()
void calc_accelMDG_mean_std(BigReal testV, int step_n, BigReal *Vmax, BigReal *Vmin, BigReal *Vavg, BigReal *M2, BigReal *sigmaV)
void adaptTempInit(int step)
__global__ void const int const TileList *__restrict__ TileExcl *__restrict__ const int *__restrict__ const int const float2 *__restrict__ const int *__restrict__ const float3 const float3 const float3 const float4 *__restrict__ const float cudaTextureObject_t cudaTextureObject_t cudaTextureObject_t const int const float const PatchPairRecord *__restrict__ const int *__restrict__ const int2 *__restrict__ const unsigned int *__restrict__ unsigned int *__restrict__ int *__restrict__ int *__restrict__ TileListStat *__restrict__ const BoundingBox *__restrict__ float4 *__restrict__ float4 *__restrict__ float *__restrict__ float *__restrict__ float *__restrict__ float *__restrict__ float *__restrict__ float *__restrict__ float *__restrict__ float *__restrict__ const int numPatches
char adaptTempInFile[NAMD_FILENAME_BUFFER_SIZE]
PressureProfileReduction(int rtag, int numslabs, int numpartitions, const char *myname, int outputfreq)
SimParameters * simParameters
void printMinimizeEnergies(int)
BigReal multigratorPressureTarget
int accelMDGEquiPrepSteps
double stochRescaleCoefficient()
BigReal surfaceTensionTarget
void writeExtendedSystemLabels(ofstream_namd &file)
void enqueueVelocities(int seq)
std::vector< BigReal > multigratorOmega
BigReal sum_of_squared_gaussians(int64_t num_gaussians)
Return the sum of i.i.d. squared Gaussians.
BigReal electEnergySlow_f
std::vector< BigReal > multigratorNuT
int berendsenPressureFreq
std::ostream & iWARN(std::ostream &s)
SubmitReduction * willSubmit(int setID, int size=-1)
if(ComputeNonbondedUtil::goMethod==2)
SimpleBroadcastObject< BigReal > adaptTemperature
Tensor groupPressure_tavg
Bool langevinPistonBarrier
BigReal * adaptTempPotEnergyAveNum
static ReductionMgr * Object(void)
#define GET_TENSOR(O, R, A)
int num_fixed_atoms() const
BigReal recent_dEdl_bond_2
void enqueuePositions(int seq, Lattice &lattice)
SimpleBroadcastObject< BigReal > tcoupleCoefficient
void outputPatchComputeMaps(const char *filename, int tag)
static int forceNeeded(int)
~PressureProfileReduction()
void swap(Array< T > &s, Array< T > &t)
void split(int iStream, int numStreams)
RequireReduction * amd_reduction
char outputFilename[NAMD_FILENAME_BUFFER_SIZE]
void multigratorTemperature(int step, int callNumber)
int64_t num_deg_freedom(int isInitialReport=0) const
void gather_energies(IMDEnergies *energies)
SimpleBroadcastObject< BigReal > stochRescaleCoefficient
void langevinPiston1(int)
void outputTiEnergy(int step)
static char * FEPTITLE2(int X)
Tensor langevinPiston_strainRate
BigReal langevinPistonDecay
std::vector< BigReal > multigratorNu
BigReal adaptTempDTavenum
Lattice checkpoint_lattice
PressureProfileReduction * ppint
void rescale(Tensor factor)
ofstream_namd adaptTempRestartFile
BigReal berendsenPressureTarget
ControllerState checkpoint_state