Re: Problems with energy conservation in NVE dynamics with 2fs timestep

From: Peter Freddolino (petefred_at_ks.uiuc.edu)
Date: Thu Mar 08 2007 - 07:03:46 CST

Hi Anna,
while I don't generally use rigidbonds, I can suggest a couple other
factors that may be destabilizing your energy:
-there isn't much of a cushion between your pairlistdist and cutoff,
which can lead to problems if atoms move quickly. This is a higher risk
with longer timesteps
-You have a pretty sparse PME grid, less than one point per angstrom,
which could lead to some drift

Also, you're using prime numbers of pme grid points, which, while it
won't affect accuracy, will slow things down; you're best off using
numbers of grid points with small integer factors.

Peter

ahenneck_at_ix.urz.uni-heidelberg.de wrote:
> Dear NAMD users
>
> I observe a total energy increase of about 70kcal/ns in my NVE
> dynamics simulations (TOTAL, TOTAL2 and TOTAL3 are increasing,
> potential and kinetic energy are stable).
> This does not occur with a timestep of 1 fs, however since I fix all
> bonds to hydrogen atoms, a 2fs timestep should work as well.
>
> If anybody has run stable NVE simulations using the CHARMM22 force
> field, I'd be happy to get the simulation parameters optimized to this
> force field or a suggestion on how to better conserve the total energy.
>
> Cheers
> Anna
>
> My setup:
> - short (11 residue) aspartic-acid-rich peptide in explicit water
> (periodic boundary conditions)
> - The system was heated slowly. Equilibration was done in three steps,
> gradually releasing constraints from the peptide and using velocity
> rescaling every 500 timesteps
>
> #############################################################
> ## SIMULATION PARAMETERS ##
> #############################################################
>
> # Input
> paraTypeCharmm on
> parameters /home/par_all22_prot.inp
> # temperature 300
>
> # Force-Field Parameters
> exclude scaled1-4
> 1-4scaling 1.0
> cutoff 10.
> switching on
> switchdist 9.
> pairlistdist 11.5
>
>
> # Integrator Parameters
> timestep 2.0
> rigidBonds all
> nonbondedFreq 1
> fullElectFrequency 2
> stepspercycle 10
>
>
> # Constant Temperature Control
> langevin off
> #langevinDamping 1
> #langevinTemp 300
> #langevinHydrogen off
>
>
> # Periodic Boundary Conditions
> cellBasisVector1 70. 0. 0.
> cellBasisVector2 0. 70. 0.
> cellBasisVector3 0. 0 70.
> cellOrigin 0. 0. 0.
>
> # Image center all molecules
> wrapAll on
> wrapNearest on
>
> # PME (for full-system periodic electrostatics)
> PME yes
> PMEGridSizeX 61
> PMEGridSizeY 61
> PMEGridSizeZ 61
> # Order default 4 (cubic):
> PMEInterpOrder 6

This archive was generated by hypermail 2.1.6 : Wed Feb 29 2012 - 15:44:27 CST