[gmx-users] Running simulation differences

Gonzalez Fernandez, Cristina cristina.gonzalezfdez at unican.es
Wed Dec 26 09:50:32 CET 2018


Hi Justin,

Tha NAMD settings for vdw interactions are: 
Rvdw=1.2 nm
Smooth switching function between 1 and 1.2 nm
In this case, the CHARMM FF was used.

However, one of the paper authors repeated the simulation of the same lipid but using the GROMOS simulation package. In this case, he used the GROMOS96FF and obtained very close values to the NAMD simulation. He included the topology as supporting information and I took it, thus I think my mistake comes from the input to the energy function. 
In that case the settings are:
Pairlist generated every 5th,
Cut-off for short range pairlist=0.8
Cut-off for long range interactions=1.4

He also included the simulations settings. I include the settings related to vdw and coulomb interactions in GROMOS in case you could help me to identify what I am doing wrong:

# use grid based pairlist generation to speed up
PAIRLIST
#	algorithm: standard(0) (gromos96 like pairlist)
#		             grid(1) (XX grid pairlist)
#	SIZE:	   grid cell size (or auto = 0.5 * RCUTP)
#	TYPE:	   chargegoup(0) (chargegroup based cutoff)
#			         atomic(1) (atom based cutoff)
#
#	algorithm	  NSNB	RCUTP	RCUTL	  SIZE	TYPE
	        1	   5	  0.8	  1.4	   0.4	   0
END
# Longrange reaction field correction
NONBONDED
# NLRELE    APPAK      RCRF     EPSRF   NSLFEXCL
       1      0.0       1.4        61    1
# NSHAPE   ASHAPE    NA2CLC   TOLA2   EPSLS
       3       1.4        2   1e-10       0
# NKX    NKY   NKZ    KCUT
   10     10    10     100
# NGX   NGY   NGZ  NASORD  NFDORD   NALIAS  NSPORD
   32    32    32       3       2        3       4
# NQEVAL   FACCUR   NRDGRD   NWRGRD   NLRLJ    SLVDNS
  100000      1.6        0        0       0      33.3
END


With this information I set in GROMACS:
ns_type                = grid        ; search neighboring grid cells
nstlist                = 5        ; 20 fs, largely irrelevant with Verlet
rcoulomb            = 1.4        ; short-range electrostatic cutoff (in nm)
rvdw                = 1.4        ; short-range van der Waals cutoff (in nm)
rlist                   = 1.41
vdwtype                 = Cut-off
vdw_modifier            = None
verlet-buffer-drift     = -1


Could you help me to find what are the inputs of the wdw energy function that are not correct? 

Thank you for all your help,

C.
-----Mensaje original-----
De: gromacs.org_gmx-users-bounces at maillist.sys.kth.se <gromacs.org_gmx-users-bounces at maillist.sys.kth.se> En nombre de Justin Lemkul
Enviado el: viernes, 21 de diciembre de 2018 14:40
Para: Discussion list for GROMACS users <gmx-users at gromacs.org>
Asunto: Re: [gmx-users] Running simulation differences

On Thu, Dec 20, 2018 at 2:58 PM Gonzalez Fernandez, Cristina < cristina.gonzalezfdez at unican.es> wrote:

>
> Hi Justin,
>
>
> Sorry for the misunderstanding. I am simulating the interaction of a 
> lipid with water. I obtain practically the same value for the coulomb 
> interaction energy, but I am not able to get a value for the van der 
> waals energy close to the published one (although I have used the same 
> settings than the ones used in the article). Therefore, I checked tha 
> GMX manual and I tried to modified the settings that affect the van 
> der waals energy in order to get a better value: vdwtype, 
> vdw-modifier, rvdw and DispCorr (as I am using
> GROMOS96 FF I don't consider rvdw-switch).
>
>
What are your settings in NAMD? If you're getting a different result then either the input settings for the energy function are different or there is some difference in the system's topology.

-Justin


>
> This is the .mdp file:
>
>
> title                = Lipid A  simulation
> ; Run parameters
> integrator          = md        ; leap-frog integrator
> nsteps                = 5500000        ; 2 * 50000 = 100 ps
> dt                    = 0.002        ; 2 fs
> ; Output control
> nstxout                = 5000        ; save coordinates every 1.0 ps
> nstvout                 = 5000        ; save velocities every 1.0 ps
> nstenergy           = 5000        ; save energies every 1.0 ps
> nstlog                = 5000        ; update log file every 1.0 ps
> nstxtcout               = 5000
> xtc-grps                = System
> ; Bond parameters
> continuation            = yes        ; restarting after NVT
> constraint_algorithm    = lincs        ; holonomic constraints
> constraints            = all-bonds    ; all bonds (even heavy atom-H
> bonds) constrained
> lincs_iter            = 1            ; accuracy of LINCS
> lincs_order            = 4            ; also related to accuracy
> ; Nonbonded settings
> cutoff-scheme           = Verlet
> ns_type                = grid        ; search neighboring grid cells
> nstlist                = 5        ; 20 fs, largely irrelevant with Verlet
> rcoulomb            = 1.4        ; short-range electrostatic cutoff (in nm)
> rvdw                = 1.4        ; short-range van der Waals cutoff (in nm)
> rlist                   = 1.41
> vdwtype                 = Cut-off
> vdw_modifier            = None
> verlet-buffer-drift     = -1
> ; Electrostatics
> coulombtype            = Reaction-Field    ; Particle Mesh Ewald for
> long-range electrostatics
> epsilon_rf              = 61
> ; Temperature coupling is on
> tcoupl           = nose-hoover    ; modified Berendsen thermostat
> tc-grps                = System    ; two coupling groups - more accurate
> tau_t                = 0.2               ; time constant, in ps
> ref_t                = 298                ; reference temperature, one for
> each group, in K
> ; Pressure coupling is oN
> pcoupl                  = Parrinello-Rahman         ; Pressure coupling on
> in NPT
> pcoupltype              = isotropic                 ; uniform scaling of
> box vectors
> tau_p                   = 0.5                       ; time constant, in ps
> ref_p                   = 1.0                       ; reference pressure,
> in bar
> compressibility         = 4.5e-5                    ; isothermal
> compressibility of water, bar^-1
> ;refcoord_scaling        = com
> ; Periodic boundary conditions
> pbc                = xyz            ; 3-D PBC
> ;Dispersion correction
> ;DispCorr            = EnerPres    ; account for cut-off vdW scheme
> ; Velocity generation
> gen_vel                 = no        ; assign velocities from Maxwell
> distribution
> ;gen_temp            = 278        ; temperature for Maxwell distribution
> ;gen_seed            = -1        ; generate a random seed
>
> I would greatly appreciate if you could help me, because I don't know 
> what else I can do,
>
> C.
>
>
>
> ________________________________
> De: gromacs.org_gmx-users-bounces at maillist.sys.kth.se < 
> gromacs.org_gmx-users-bounces at maillist.sys.kth.se> en nombre de Justin 
> Lemkul <jalemkul at vt.edu>
> Enviado: jueves, 20 de diciembre de 2018 19:52:20
> Para: Discussion list for GROMACS users
> Asunto: Re: [gmx-users] Running simulation differences
>
> On Thu, Dec 20, 2018 at 1:01 PM Gonzalez Fernandez, Cristina < 
> cristina.gonzalezfdez at unican.es> wrote:
>
> > Hi Justin,
> >
> >
> > Do you think that the difference between a value of -280 ± 37 kj/mol 
> > obtained by Gromacs and -367 kj/mol obtained by NAMD can be acceptable?
> >
> >
> As I said before, no, a massive difference like that, outside the 
> error bars, is not acceptable.
>
> If you want more useful advice, you'll have to tell us exactly what 
> you're doing and trying to compare.
>
> -Justin
>
>
> >
> > Thank you so much for all your help
> >
> > C.
> >
> > ________________________________
> > De: gromacs.org_gmx-users-bounces at maillist.sys.kth.se < 
> > gromacs.org_gmx-users-bounces at maillist.sys.kth.se> en nombre de 
> > Justin Lemkul <jalemkul at vt.edu>
> > Enviado: martes, 18 de diciembre de 2018 15:34:07
> > Para: Discussion list for GROMACS users
> > Asunto: Re: [gmx-users] Running simulation differences
> >
> > On Tue, Dec 18, 2018 at 5:01 AM Gonzalez Fernandez, Cristina < 
> > cristina.gonzalezfdez at unican.es> wrote:
> >
> > > Hi again Justin,
> > >
> > > I've been thinking about the normal variability between MD 
> > > simulations
> > you
> > > commented. I am simulating the interaction of a lipid with water, 
> > > but
> the
> > > van der waals energy value that I obtain (-257.25  ± 35.86 kj/mol) 
> > > for
> me
> > > is considerably different from the published value (-361 kj/mol).
> > However,
> > > according to your email maybe these values can be considered similar.
> Do
> > > you think my value and the published one are similar? How can I 
> > > reduced
> > the
> > > error of the results because my value correspond to a long
> >
> >
> > I would not call those values similar at all. The two quantities you
> posted
> > before were within error of one another; those are similar. A 
> > difference
> of
> > over 100 kJ/mol, outside error, is enormous. You need to evaluate 
> > what you're doing because in that case, I suspect something is 
> > inconsistent between what is published and what you're doing.
> >
> >
> > > simulation (60 ns)? Could you recommend me a paper where is stated 
> > > the common deviations/errors that are acceptable to a simulation result?
> > >
> >
> > There is none to my knowledge, and certainly can't be generalized.
> > Different quantities have different tolerances. Different software 
> > may produce somewhat different results, as well.
> >
> > Also I wouldn't consider 60 ns long, by any means :)
> >
> > -Justin
> >
> > --
> >
> > ==========================================
> >
> > Justin A. Lemkul, Ph.D.
> >
> > Assistant Professor
> >
> > Office: 301 Fralin Hall
> > Lab: 303 Engel Hall
> >
> > Virginia Tech Department of Biochemistry
> > 340 West Campus Dr.
> > Blacksburg, VA 24061
> >
> > jalemkul at vt.edu | (540) 231-3129
> > http://www.thelemkullab.com
> >
> >
> > ==========================================
> > --
> > Gromacs Users mailing list
> >
> > * Please search the archive at
> > http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before 
> > posting!
> >
> > * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
> >
> > * For (un)subscribe requests visit
> > https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users 
> > or send a mail to gmx-users-request at gromacs.org.
> > --
> > Gromacs Users mailing list
> >
> > * Please search the archive at
> > http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before 
> > posting!
> >
> > * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
> >
> > * For (un)subscribe requests visit
> > https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users 
> > or send a mail to gmx-users-request at gromacs.org.
> >
>
>
> --
>
> ==========================================
>
> Justin A. Lemkul, Ph.D.
>
> Assistant Professor
>
> Office: 301 Fralin Hall
> Lab: 303 Engel Hall
>
> Virginia Tech Department of Biochemistry
> 340 West Campus Dr.
> Blacksburg, VA 24061
>
> jalemkul at vt.edu | (540) 231-3129
> http://www.thelemkullab.com
>
>
> ==========================================
> --
> Gromacs Users mailing list
>
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before 
> posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or 
> send a mail to gmx-users-request at gromacs.org.
> --
> Gromacs Users mailing list
>
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before 
> posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or 
> send a mail to gmx-users-request at gromacs.org.
>


-- 

==========================================

Justin A. Lemkul, Ph.D.

Assistant Professor

Office: 301 Fralin Hall
Lab: 303 Engel Hall

Virginia Tech Department of Biochemistry
340 West Campus Dr.
Blacksburg, VA 24061

jalemkul at vt.edu | (540) 231-3129
http://www.thelemkullab.com


==========================================
--
Gromacs Users mailing list

* Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-request at gromacs.org.


More information about the gromacs.org_gmx-users mailing list