[gmx-developers] energy calculation fails in "-rerun" mode

Wed May 27 05:12:33 CEST 2009

The problem appears when I am trying to recalculate energies and dhdl values
using "-rerun" option. Recalculated values do not coincide with original
(from MD run) values. The difference is definitely larger then the
uncertainty due to trajectory rounding (trajectory precision is 3 digits).
Trajectory has been obtained for the solute (Cl+,CL-) immersed in the
polarizable solvent (benzene).  The solvent is modeled by the polarizable
Drude oscillator model. Shell particles are used for benzene molecules. The
"rerun" procedure results to a strange output.

1) Average values are not printed out at the end of the log-file. There is
only some message at the beginning: "Will not sum the energies at every
step, therefore the energy file does not contain exact averages and
fluctuations."
But I didn't turn on the option -nosum.

2) Output energies are not updated for each configuration of my trajectory,
but only once per 20 snapshots, i.e. 20 consequent configurations have
exactly the same energies.

The points 1) and 2) seems to be related. There is also question regarding
minimization of the shell particle positions in "-rerun" mode.

3) My trajectory already have minimized (during MD) shell positions. I only
need to calculate energies corresponding to a given configuration,
therefore, I don't need minimize the shell particles. How can I turn off the
minimization in "-rerun" mode? I tried to specify "niter=0", but it results
to a strange output (unphysically large forces which are increased on each
step). Below is my mdp- and log-files:

==================================================================
nstlog = 1 ; print to logfile
energygrps = IPR ; energy groups
nstenergy = 0 ; print energies
nstlist = 1 ; update pairlist
ns_type = grid ;simple ; pairlist method
emtol = 0.8 ;the convergency criterion for maximum force
niter =100
fcstep = 0.001 ; shell particle minimization timestep
lincs_order = 8
lincs_iter = 4
unconstrained_start = yes
free_energy = yes ; perform slow growth
init_lambda = 1.0
delta_lambda = 0 ; increment per time step
pbc = xyz ; periodic boundary conditions
optimize_fft = yes ; perform FFT optimization at start
coulombtype = PME; Cut-off
rlist = 1.21 ; cut-off for ns
rvdw = 1.21 ; cut-off for vdw
rcoulomb = 1.21 ; cut-off for coulomb
Tcoupl = no;
==================================================================

:-) G R O M A C S (-:

Getting the Right Output Means no Artefacts in Calculating Stuff

:-) VERSION 4.0.4 (-:

parameters of the run:
integrator = md
nsteps = 0
init_step = 0
ns_type = Grid
nstlist = 1
ndelta = 2
nstcomm = 1
comm_mode = Linear
nstlog = 1
nstxout = 100000
nstvout = 100000
nstfout = 0
nstenergy = 0
nstxtcout = 0
init_t = 0
delta_t = 0.001
xtcprec = 1000
nkx = 32
nky = 32
nkz = 32
pme_order = 4
ewald_rtol = 1e-05
ewald_geometry = 0
epsilon_surface = 0
optimize_fft = TRUE
ePBC = xyz
bPeriodicMols = FALSE
bContinuation = TRUE
bShakeSOR = FALSE
etc = No
epc = No
epctype = Isotropic
tau_p = 1
ref_p (3x3):
ref_p[ 0]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
ref_p[ 1]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
ref_p[ 2]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
compress (3x3):
compress[ 0]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
compress[ 1]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
compress[ 2]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
refcoord_scaling = No
posres_com (3):
posres_com[0]= 0.00000e+00
posres_com[1]= 0.00000e+00
posres_com[2]= 0.00000e+00
posres_comB (3):
posres_comB[0]= 0.00000e+00
posres_comB[1]= 0.00000e+00
posres_comB[2]= 0.00000e+00
andersen_seed = 815131
rlist = 1.21
rtpi = 0.05
coulombtype = PME
rcoulomb_switch = 0
rcoulomb = 1.21
vdwtype = Cut-off
rvdw_switch = 0
rvdw = 1.21
epsilon_r = 1
epsilon_rf = 1
tabext = 1
implicit_solvent = No
gb_algorithm = Still
gb_epsilon_solvent = 80
nstgbradii = 1
rgbradii = 2
gb_saltconc = 0
gb_obc_alpha = 1
gb_obc_beta = 0.8
gb_obc_gamma = 4.85
sa_surface_tension = 2.092
DispCorr = No
free_energy = yes
init_lambda = 1
sc_alpha = 0
sc_power = 0
sc_sigma = 0.3
delta_lambda = 0
nwall = 0
wall_type = 9-3
wall_atomtype[0] = -1
wall_atomtype[1] = -1
wall_density[0] = 0
wall_density[1] = 0
wall_ewald_zfac = 3
pull = no
disre = No
disre_weighting = Conservative
disre_mixed = FALSE
dr_fc = 1000
dr_tau = 0
nstdisreout = 100
orires_fc = 0
orires_tau = 0
nstorireout = 100
dihre-fc = 1000
em_stepsize = 0.01
em_tol = 0.8
niter = 20
fc_stepsize = 0.001
nstcgsteep = 1000
nbfgscorr = 10
ConstAlg = Lincs
shake_tol = 1e-04
lincs_order = 8
lincs_warnangle = 30
lincs_iter = 4
bd_fric = 0
ld_seed = 1993
cos_accel = 0
deform (3x3):
deform[ 0]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
deform[ 1]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
deform[ 2]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
userint1 = 0
userint2 = 0
userint3 = 0
userint4 = 0
userreal1 = 0
userreal2 = 0
userreal3 = 0
userreal4 = 0
grpopts:
nrdf: 12603
ref_t: 0
tau_t: 0
anneal: No
ann_npoints: 0
acc: 0 0 0
nfreeze: N N N
energygrp_flags[ 0]: 0 0
energygrp_flags[ 1]: 0 0
efield-x:
n = 0
efield-xt:
n = 0
efield-y:
n = 0
efield-yt:
n = 0
efield-z:
n = 0
efield-zt:
n = 0
bQMMM = FALSE
QMconstraints = 0
QMMMscheme = 0
scalefactor = 1
qm_opts:
ngQM = 0
Initializing Domain Decomposition on 8 nodes
Dynamic load balancing: no
Will sort the charge groups at every domain (re)decomposition
Initial maximum inter charge-group distances:
two-body bonded interactions: 0.317 nm, Bond, atoms 5967 5961
multi-body bonded interactions: 0.365 nm, Proper Dih., atoms 276 285
Minimum cell size due to bonded interactions: 0.401 nm
Using 0 separate PME nodes
Optimizing the DD grid for 8 cells with a minimum initial size of 0.401 nm
The maximum allowed number of cells is: X 9 Y 9 Z 9
Domain decomposition grid 8 x 1 x 1, separate PME nodes 0
Domain decomposition nodeid 0, coordinates 0 0 0
Table routines are used for coulomb: TRUE
Table routines are used for vdw: FALSE
Will do PME sum in reciprocal space.
++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++
U. Essman, L. Perela, M. L. Berkowitz, T. Darden, H. Lee and L. G. Pedersen
A smooth particle mesh Ewald method
J. Chem. Phys. 103 (1995) pp. 8577-8592
-------- -------- --- Thank You --- -------- --------
Using a Gaussian width (1/beta) of 0.387397 nm for Ewald
Cut-off's: NS: 1.21 Coulomb: 1.21 LJ: 1.21
System total charge, top. A: -0.000 top. B: -0.000
Generated table with 1105 data points for Ewald.
Tabscale = 500 points/nm
Generated table with 1105 data points for LJ6.
Tabscale = 500 points/nm
Generated table with 1105 data points for LJ12.
Tabscale = 500 points/nm
Configuring nonbonded kernels...
Testing x86_64 SSE support... present.
There are 2 atoms and 2 charges for free energy perturbation
Linking all bonded interactions to atoms
There are 25200 inter charge-group exclusions,
will use an extra communication step for exclusion forces for PME
The initial number of communication pulses is: X 3
The initial domain decomposition cell size is: X 0.47 nm
The maximum allowed distance for charge groups involved in interactions is:
non-bonded interactions 1.210 nm
two-body bonded interactions (-rdd) 1.210 nm
multi-body bonded interactions (-rdd) 0.472 nm
Making 1D domain decomposition grid 8 x 1 x 1, home cell index 0 0 0
Will not sum the energies at every step,
therefore the energy file does not contain exact averages and fluctuations.
WARNING:
We should not remove the COM motion every step with option -nosum,
setting nstcomm to nstlist (1)
Center of mass motion removal mode is Linear
We have the following groups for center of mass motion removal:
0: rest
There are: 4202 Atoms
There are: 2100 Shells
Charge group distribution at step 0: 263 256 267 260 279 263 261 253
Grid: 8 x 6 x 6 cells
Initial temperature: 292.107 K
Started mdrun on node 0 Wed May 27 05:14:03 2009
Charge group distribution at step 0: 263 256 267 260 279 263 261 253

Step Time Lambda
0 0.00000 1.00000

step 0: EM did not converge in 20 iterations, RMS force 3.237
Energies (kJ/mol)
Bond U-B Proper Dih. LJ (SR) Coulomb (SR)
6.64705e+03 4.09717e+03 3.66230e+03 4.94161e+03 9.54908e+02
Coul. recip. Polarization Thole Pol. Potential Kinetic En.
-6.16992e+02 2.29306e+02 2.62844e+01 1.99416e+04 1.53046e+04
Total Energy Temperature Pressure (bar) dVpot/dlambda dEkin/dlambda
3.52463e+04 2.92107e+02 3.48768e+01 -5.41851e+02 0.00000e+00

DD step 3999 load imb.: force 8.2%
Charge group distribution at step 4000: 262 256 267 260 279 263 261 254
Step Time Lambda
4000 2.00000 1.00000

step 4000: EM did not converge in 20 iterations, RMS force 3.440
Energies (kJ/mol)
Bond U-B Proper Dih. LJ (SR) Coulomb (SR)
6.64705e+03 4.09717e+03 3.66230e+03 4.94162e+03 9.55294e+02
Coul. recip. Polarization Thole Pol. Potential Kinetic En.
-6.17053e+02 2.30621e+02 2.52502e+01 1.99423e+04 1.53046e+04
Total Energy Temperature Pressure (bar) dVpot/dlambda dEkin/dlambda
3.52469e+04 2.92107e+02 3.47506e+01 -5.41880e+02 0.00000e+00
DD step 7999 load imb.: force 8.3%
Charge group distribution at step 8000: 262 256 267 260 279 263 261 254
Step Time Lambda
8000 4.00000 1.00000
step 8000: EM did not converge in 20 iterations, RMS force 3.440
Energies (kJ/mol)
Bond U-B Proper Dih. LJ (SR) Coulomb (SR)
6.64705e+03 4.09717e+03 3.66230e+03 4.94162e+03 9.55294e+02
Coul. recip. Polarization Thole Pol. Potential Kinetic En.
-6.17053e+02 2.30621e+02 2.52502e+01 1.99423e+04 1.53046e+04
Total Energy Temperature Pressure (bar) dVpot/dlambda dEkin/dlambda
3.52469e+04 2.92107e+02 3.47506e+01 -5.41880e+02 0.00000e+00

-//-//-//-//

Charge group distribution at step 100000: 262 256 267 260 279 263 261 254
Step Time Lambda
100000 50.00000 1.00000

step 100000: EM did not converge in 20 iterations, RMS force 3.440
Energies (kJ/mol)
Bond U-B Proper Dih. LJ (SR) Coulomb (SR)
6.64705e+03 4.09717e+03 3.66230e+03 4.94162e+03 9.55294e+02
Coul. recip. Polarization Thole Pol. Potential Kinetic En.
-6.17053e+02 2.30621e+02 2.52502e+01 1.99423e+04 1.53046e+04
Total Energy Temperature Pressure (bar) dVpot/dlambda dEkin/dlambda
3.52469e+04 2.92107e+02 3.47506e+01 -5.41880e+02 0.00000e+00

DD step 103999 load imb.: force 15.1%
Step Time Lambda
104000 52.00000 1.00000

step 104000: EM did not converge in 20 iterations, RMS force 4.451
Energies (kJ/mol)
Bond U-B Proper Dih. LJ (SR) Coulomb (SR)
6.64705e+03 4.09717e+03 3.66230e+03 4.94161e+03 9.56206e+02
Coul. recip. Polarization Thole Pol. Potential Kinetic En.
-6.17090e+02 2.31805e+02 2.42288e+01 1.99433e+04 1.53046e+04
Total Energy Temperature Pressure (bar) dVpot/dlambda dEkin/dlambda
3.52479e+04 2.92107e+02 3.45236e+01 -5.41907e+02 0.00000e+00

-//-//-//-//

Step Time Lambda
10000000 5000.00000 1.00000

step 10000000: EM did not converge in 20 iterations, RMS force 5.071
Energies (kJ/mol)
Bond U-B Proper Dih. LJ (SR) Coulomb (SR)
6.64705e+03 4.09717e+03 3.66230e+03 4.94161e+03 9.56047e+02
Coul. recip. Polarization Thole Pol. Potential Kinetic En.
-6.17083e+02 2.32676e+02 2.35112e+01 1.99433e+04 1.53046e+04
Total Energy Temperature Pressure (bar) dVpot/dlambda dEkin/dlambda
3.52479e+04 2.92107e+02 3.44392e+01 -5.41893e+02 0.00000e+00

Fraction of iterations that converged: 0.00 %
Average number of force evaluations per MD step: 0.01

M E G A - F L O P S A C C O U N T I N G
RF=Reaction-Field FE=Free Energy SCFE=Soft-Core/Free Energy
T=Tabulated W3=SPC/TIP3p W4=TIP4p (single or pairs)
NF=No Forces
Computing: M-Number M-Flops % Flops
-----------------------------------------------------------------------
LJ 420.168000 13865.544 0.2
Coul(T) 74944.363200 3147663.254 38.3
Coul(T) + LJ 60207.591360 3311417.525 40.3
Free energy innerloop 89.185660 13377.849 0.2
Outer nonbonded loop 1291.216280 12912.163 0.2
Calc Weights 945.678120 34044.412 0.4
Spread Q Bspline 40348.933120 80697.866 1.0
Gather F Bspline 40348.933120 484187.197 5.9
3D-FFT 98343.321600 786746.573 9.6
Solve PME 1741.496320 111455.764 1.4
NS-Pairs 1644.786945 34540.526 0.4
CG-CoM 15.767604 47.303 0.0
Bonds 210.084000 12394.956 0.2
Propers 420.168000 96218.472 1.2
Thole Pol. 210.084000 62184.864 0.8
Virial 333.233240 5998.198 0.1
Calc-Ekin 15.767604 425.725 0.0
-----------------------------------------------------------------------
Total 8208178.193 100.0
-----------------------------------------------------------------------

D O M A I N D E C O M P O S I T I O N S T A T I S T I C S
av. #atoms communicated per step for force: 2 x 16176.0
Average load imbalance: 8.6 %
Part of the total run time spent waiting due to load imbalance: 6.2 %
NOTE: 6.2 % performance was lost due to load imbalance
in the domain decomposition.
You might want to use dynamic load balancing (option -dlb.)

R E A L C Y C L E A N D T I M E A C C O U N T I N G
Computing: Nodes Number G-Cycles Seconds %
-----------------------------------------------------------------------
Domain decomp. 8 2501 114.105 49.0 1.1
Comm. coord. 8 50020 91.331 39.2 0.9
Neighbor search 8 2501 247.347 106.2 2.4
Force 8 50020 7503.545 3222.5 72.6
Wait + Comm. F 8 50020 176.649 75.9 1.7
PME mesh 8 50020 2045.186 878.3 19.8
Write traj. 8 101 1.491 0.6 0.0
Comm. energies 8 2501 15.381 6.6 0.1
Rest 8 143.303 61.5 1.4
-----------------------------------------------------------------------
Total 8 10338.337 4440.0 100.0
----------------------------------------------------------------------- 