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

Wed May 27 05:19:28 CEST 2009

Leontyev Igor wrote:
> The problem appears when I am trying to recalculate energies and dhdl 
> values
> using "-rerun" option.

There were bugfixes for -rerun between 4.0.4 and 4.0.5. Please check 
release notes for this sort of thing :-) I don't know if you might still 
have a real issue, but please update to 4.0.5 and try again.

> 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).

Unless you did neighbour-searching every step of your original 
trajectory, then you will likely get different results, since nstlist != 
0 implies neighbour searching every step for reruns.

Mark

> 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
> -----------------------------------------------------------------------
> _______________________________________________
> gmx-developers mailing list
> gmx-developers at gromacs.org
> http://www.gromacs.org/mailman/listinfo/gmx-developers
> Please don't post (un)subscribe requests to the list. Use the www 
> interface or send it to gmx-developers-request at gromacs.org.
> 