Re: How to conservate energy during NVE simulations???

From: Marcos Sotomayor (sotomayo_at_ks.uiuc.edu)
Date: Thu Jun 16 2005 - 10:30:28 CDT

Hi Leonardo,

The energy drift observed in your simulations is probably caused by the
time steps you are using. If you want to have a very good energy
conservation, a *uniform* 1 femtosecond time step should work (you could
also try 2 fs with nonbondFreq 1, fullElectFrequency 1 and rigidBonds).

Marcos

PS: Good references about energy conservation and use of multiple timestep
algorithms are

[3] Nonlinear Resonance Artifacts in Molecular Dynamics Simulations. T.
Schlick, M. Mandziuk, R. D. Skeel and K. Srinivas. J. Comput. Phys.,
140:1-29, 1998
[4] Verlet-I/r-RESPA Is Limited by Nonlinear Instability. Q. Ma, J.
Izaguirre, and R. D. Skeel. SIAM J. Sci. Comput., 24(6), 1951-1973, 2003
[5] Difficulties with Multiple Timestepping and the Fast Multipole
Algorithm in Molecular Dynamics. T. Bishop and Robert D. Skeel and K.
Schulten. J. Comp. Chem., 18(14):1785-1791, 1997.
[6] Langevin Stabilization of Molecular Dynamics J. A. Izaguirre, D. P.
Catarello, J. M. Wozniak, and R. D. Skeel. J. Chem. Phys., 114(4),
2090-2098, 2000.

On Wed, 15 Jun 2005, [ISO-8859-1] Leonardo Sepulveda Durán wrote:

> Hello everybody
>
> I was trying to do a NVE simulation in NAMD. I used imput coordinate
> and velocity files from a NVT simulation which was equilibrated at
> desired temperature (498 K)
>
> Looking at TOTE3, it is obvious my system's energy underwent a slow
> drift (from -20235 to -20225 in 250 ps ). The .conf I used is below:
>
> structure ./ubq-wb.psf
> coordinates ./ubq-min.coor
>
> # For Continuing a run, otherwise comment all but firsttimestep
> set inputname ubq-heat
> binCoordinates $inputname.coor ;# coordinates for last run (binary))
> binVelocities $inputname.vel ;# velocities for last run (binary)
> extendedSystem $inputname.xsc ;# Cell dymentions for last run
> firsttimestep 41500 ;# Last step previous run
> numsteps 125000 ;# run stops when this
> step is reached
>
> # Force-Field Parameters
> paraTypeCharmm on
> parameters /home/RPMs/NAMD/toppar/par_all22_prot.inp
> exclude scaled1-4
> 1-4scaling 1.0
> switching on
>
> switchdist 10.; # cutoff-2
> cutoff 12.; # may be 10 with PME
> pairlistdist 14.; # cutoff+2
> margin 0.498;
>
> # Integrator Parameters
> timestep 2.0; # 2fs/step
> rigidBonds all; # needed for 2fs steps (SHAKE)
> nonbondedFreq 1
> fullElectFrequency 2
> stepspercycle 10
> outputPairlists 1000
>
> # Periodic Boundary Conditions
> cellBasisVector1 50.68 0. 0.
> cellBasisVector2 0. 52.37 0.
> cellBasisVector3 0. 0. 56.27
> cellOrigin 30.85 29.14 17.90
>
> wrapAll on
>
> # PME (for full-system periodic electrostatics)
> PME yes
> PMEGridSizeX 54
> PMEGridSizeY 54
> PMEGridSizeZ 60
>
> # Output
> outputName ubq-eq
> binaryOutput yes;
>
> restartfreq 500; # 500steps = every 1ps
> dcdfreq 500
> xstFreq 500
> outputEnergies 100
> outputPressure 100
> outputTiming 1000
>
> I was wondering if someone can give me an advice about how to get my
> simulations conservative. I am confused about some concepts and
> parameters...for example, in the manual, stepspercycle is explained
> as the number of timesteps beetween atom
> reasignments; I thought it was refering to pairlist reasignment, but
> later "pairlistpercycle" is defined as the number of times pairlists
> are regenerated in one cycle...so if I have
>
> stepspercycle = 20
> pairlistpercycle = 2
>
> is the same as using
>
> stepspercycle = 10
> pairlistpercycle = 1 ?????
>
> I really dont understand the difference. Also I put outputPairlists
> because it sound useful to avoid atoms to pass the pairlist boundary,
> but i am not sure if only give the pairlists violations. So I would
> like to ask if someone have done NVE simulations succesfully, and what
> parameters can be adecuate to do it.
>
> Thanks
>
> Leonardo
>
> PD: I attach a .pdf with temperature and energy from my simulation.
>

This archive was generated by hypermail 2.1.6 : Wed Feb 29 2012 - 15:40:51 CST