[gmx-users] Re: NS in PME

Janne Hirvi janne.hirvi at joensuu.fi
Tue Apr 11 11:16:22 CEST 2006


Thanks for your response David!

I tried to simulate bulk water (1372 SPC/E molecules) with your parameters and
observed same kind of problems with energy conservation as with my earlier
parameters (rlist=rvdw=rcoulomb=1.2 & nstlist=10). In single precision NVE
simulation (500ps) drift in total energy was +1828.88kJ/mol with your
parameters and +1861.13kJ/mol with mine. So I am little bit confused and I
still have couple of questions: 

1) Are this large drifts normal or have I done something stupid? I mean have
you
observed this kind of artifacts as well? Atleast from the next mail I would
think that better energy conservation should be achieved with your parameters?

http://www.gromacs.org/pipermail/gmx-users/2002-December/003510.html

2) I still dont understand why parameter combination rlist>rcoulomb with PME
conserves total energy almost completely, but however I should use parameters
so
that rlist=rcoulomb. Can you comment on this because Berk may have noted your
request or then he have had no time to do that? Or is it so trivial that I
should find explanation from some trivial text book?


Thanks,

Janne



> Janne Hirvi wrote:
> > Hello David!
> > 
> > I read your interesting new article "The Origin of Layer Structure
> Artifacts in
> > Simulations of Liquid Water" and I would like to know especially what kind
> of
> > NS parameters you have used in the simulations with PME and are the
> following
> > parameters correct otherwise?
> > 
> > dt=0.002
> > nstlist=5
> > rlist=?
> > coulombtype=PME
> > rcoulomb=0.9
> > vdwtype=cut-off
> > rvdw=?    
> I'm attaching my mdp file.
> 
> We have just disabled the option to have rlist != rcoulomb with PME.
> 
> > 
> > I am interested especially about the rlist parameter because there seems to
> be
> > some energy conservation problems in my water simulations. At least I
> couldn't
> > achieve energy conservation with single precision complilation and with
> > double precision compilation energy is conserved only when switch(shift?)
> > function is used for vdw interactions (rlist>rvdw>rvdw_switch) and PME for
> > coulombic interactions (rlist>rcoulomb). If rlist is equal to rcoulomb I
> need
> > to update NS list at every time step but with preceding parameters energy
> is
> > conserved even when nstlist=10. 
> > 
> > So should I use rlist>rcoulomb with PME when nstlist>1 to take into
> account
> > movement of atoms/molecules(charge groups) near the real space cut-off
> limit
> > (rcoulomb) or is there something I don't understand? Atleast from the
> manual
> > someone could get the picture that rlist should be equal to rcoulomb and
> there
> > are also different kind of opinions on mailing list?
> 
> Maybe Berk can comment?
> 
> > 
> > 
> > Thanks for any comments in advance! 
> > 
> > Janne
> > 
> >
>
------------------------------------------------------------------------------
> > Janne Hirvi, MSc(Physical Chemistry), Researcher
> > University of Joensuu, Department of Chemistry, P.O.Box 111 80101 Joensuu,
> FI
> > Tel: +358 13 2514544 & +358 50 3474223
> > E-mail: Janne.Hirvi at joensuu.fi & hirvi at cc.joensuu.fi
> >
>
------------------------------------------------------------------------------
> > _______________________________________________
> > gmx-users mailing list    gmx-users at gromacs.org
> > http://www.gromacs.org/mailman/listinfo/gmx-users
> > Please don't post (un)subscribe requests to the list. Use the 
> > www interface or send it to gmx-users-request at gromacs.org.
> > Can't post? Read http://www.gromacs.org/mailing_lists/users.php
> 
> 
> -- 
> David.
> ________________________________________________________________________
> David van der Spoel, PhD, Assoc. Prof., Molecular Biophysics group,
> Dept. of Cell and Molecular Biology, Uppsala University.
> Husargatan 3, Box 596,  	75124 Uppsala, Sweden
> phone:	46 18 471 4205		fax: 46 18 511 755
> spoel at xray.bmc.uu.se	spoel at gromacs.org   http://folding.bmc.uu.se
> ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
> -------------- next part --------------
> ;
> ;	File 'mdout.mdp' was generated
> ;	By user: spoel (291)
> ;	On host: matisse
> ;	At date: Sun Aug 11 22:20:40 2002
> ;
> 
> ; VARIOUS PREPROCESSING OPTIONS = 
> title                    = 
> cpp                      = /lib/cpp
> 
> ; RUN CONTROL PARAMETERS = 
> integrator               = md
> ; start time and timestep in ps = 
> tinit                    = 0
> dt                       = 0.002
> nsteps                   = 1000000
> ; mode for center of mass motion removal = 
> comm-mode                = Linear
> ; number of steps for center of mass motion removal = 
> nstcomm                  = 1
> ; group(s) for center of mass motion removal = 
> comm-grps                = 
> 
> ; LANGEVIN DYNAMICS OPTIONS = 
> ; Temperature, friction coefficient (amu/ps) and random seed = 
> bd-temp                  = 300
> bd-fric                  = 0
> ld-seed                  = 1993
> 
> ; ENERGY MINIMIZATION OPTIONS = 
> ; Force tolerance and initial step-size = 
> emtol                    = 100
> emstep                   = 0.01
> ; Max number of iterations in relax_shells = 
> niter                    = 20
> ; Step size (1/ps^2) for minimization of flexible constraints = 
> fcstep                   = 0
> ; Frequency of steepest descents steps when doing CG = 
> nstcgsteep               = 1000
> 
> ; OUTPUT CONTROL OPTIONS = 
> ; Output frequency for coords (x), velocities (v) and forces (f) = 
> nstxout                  = 0
> nstvout                  = 0
> nstfout                  = 0
> ; Output frequency for energies to log file and energy file = 
> nstlog                   = 1000
> nstenergy                = 100
> ; Output frequency and precision for xtc file = 
> nstxtcout                = 250
> xtc-precision            = 1000
> ; This selects the subset of atoms for the xtc file. You can = 
> ; select multiple groups. By default all atoms will be written. = 
> xtc-grps                 = 
> ; Selection of energy groups = 
> energygrps               = 
> 
> ; NEIGHBORSEARCHING PARAMETERS = 
> ; nblist update frequency = 
> nstlist                  = 5
> ; ns algorithm (simple or grid) = 
> ns-type                  = Grid
> ; Periodic boundary conditions: xyz or no = 
> pbc                      = xyz
> ; nblist cut-off         = 
> domain-decomposition     = no
> 
> ; OPTIONS FOR ELECTROSTATICS AND VDW = 
> ; Method for doing electrostatics = 
> coulombtype              = PME
> rcoulomb-switch          = 0
> ; Dielectric constant (DC) for cut-off or DC of reaction field = 
> epsilon-r                = 1
> ; Method for doing Van der Waals = 
> vdw-type                 = Cut-off
> ; cut-off lengths        = 
> rvdw-switch              = 0
> ; Apply long range dispersion corrections for Energy and Pressure = 
> DispCorr                 = Ener
> ; Spacing for the PME/PPPM FFT grid = 
> fourierspacing           = 0.12
> ; FFT grid size, when a value is 0 fourierspacing will be used = 
> fourier_nx               = 0
> fourier_ny               = 0
> fourier_nz               = 0
> ; EWALD/PME/PPPM parameters = 
> pme_order                = 4
> ewald_rtol               = 1e-05
> ewald_geometry           = 3d
> epsilon_surface          = 0
> optimize_fft             = no
> 
> ; OPTIONS FOR WEAK COUPLING ALGORITHMS = 
> ; Temperature coupling   = 
> tcoupl                   = Berendsen
> ; Groups to couple separately = 
> tc-grps                  = System
> ; Time constant (ps) and reference temperature (K) = 
> tau-t                    = 0.1
> ref-t                    = 298.15 
> ; Pressure coupling      = 
> Pcoupl                   = Berendsen
> Pcoupltype               = Isotropic
> ; Time constant (ps), compressibility (1/bar) and reference P (bar) = 
> tau-p                    = 1.0
> compressibility          = 5e-5
> ref-p                    = 1
> 
> ; SIMULATED ANNEALING CONTROL = 
> annealing                = no
> ; Time at which temperature should be zero (ps) = 
> zero-temp_time           = 0
> 
> ; GENERATE VELOCITIES FOR STARTUP RUN = 
> gen-vel                  = yes
> gen-temp                 = 300
> gen-seed                 = 173529
> 
> ; OPTIONS FOR BONDS     = 
> constraints              = none
> ; Type of constraint algorithm = 
> constraint-algorithm     = shake
> ; Do not constrain the start configuration = 
> unconstrained-start      = no
> ; Use successive overrelaxation to reduce the number of shake iterations = 
> Shake-SOR                = no
> ; Relative tolerance of shake = 
> shake-tol                = 1e-04
> ; Highest order in the expansion of the constraint coupling matrix = 
> lincs-order              = 4
> ; Lincs will write a warning to the stderr if in one step a bond = 
> ; rotates over more degrees than = 
> lincs-warnangle          = 30
> ; Convert harmonic bonds to morse potentials = 
> morse                    = no
> 
> ; ENERGY GROUP EXCLUSIONS = 
> ; Pairs of energy groups for which all non-bonded interactions are excluded =
> 
> energygrp_excl           = 
> 
> ; NMR refinement stuff  = 
> ; Distance restraints type: No, Simple or Ensemble = 
> disre                    = No
> ; Force weighting of pairs in one distance restraint: Conservative or Equal =
> 
> disre-weighting          = Conservative
> ; Use sqrt of the time averaged times the instantaneous violation = 
> disre-mixed              = no
> disre-fc                 = 1000
> disre-tau                = 0
> ; Output frequency for pair distances to energy file = 
> nstdisreout              = 100
> ; Orientation restraints: No or Yes = 
> orire                    = no
> ; Orientation restraints force constant and tau for time averaging = 
> orire-fc                 = 0
> orire-tau                = 0
> orire-fitgrp             = 
> ; Output frequency for trace(SD) to energy file = 
> nstorireout              = 100
> 
> ; Free energy control stuff = 
> free-energy              = no
> init-lambda              = 0
> delta-lambda             = 0
> sc-alpha                 = 0
> sc-sigma                 = 0.3
> 
> ; Non-equilibrium MD stuff = 
> acc-grps                 = 
> accelerate               = 
> freezegrps               = 
> freezedim                = 
> cos-acceleration         = 0
> 
> ; Electric fields       = 
> ; Format is number of terms (int) and for all terms an amplitude (real) = 
> ; and a phase angle (real) = 
> E-x                      = 
> E-xt                     = 
> E-y                      = 
> E-yt                     = 
> E-z                      = 
> E-zt                     = 
> 
> ; User defined thingies = 
> user1-grps               = 
> user2-grps               = 
> userint1                 = 0
> userint2                 = 0
> userint3                 = 0
> userint4                 = 0
> userreal1                = 0
> userreal2                = 0
> userreal3                = 0
> userreal4                = 0
> rlist = 0.9
> rcoulomb = 0.9
> rvdw = 0.9



More information about the gromacs.org_gmx-users mailing list