[gmx-users] TI calculation in vacuum gives number which is not zero for small molecule.
Berk Hess
gmx3 at hotmail.com
Tue Aug 3 09:36:49 CEST 2010
Hi,
You use reaction-field (even generalized rf, which I wouldn't use).
Reaction-field also has a correction for excluded pairs.
So rf, correctly, shows an effect of decoupling the molecule
from the surrounding dielectric.
PS using mdrun -seppot will show in the log file that the contribution
to dhdl is from rf-excl.
Berk
> Date: Tue, 3 Aug 2010 13:26:32 +1000
> From: itamar.kass at monash.edu
> To: gmx-users at gromacs.org
> Subject: [gmx-users] TI calculation in vacuum gives number which is not zero for small molecule.
>
> Shalom all,
>
> I wish to calculate the free energy of solvation of acetic acid using
> thermodynamics integration. So I built the topology:
>
> ; Include forcefield parameters
> #include "ffG53a6.itp"
>
> [ moleculetype ]
> ; Name nrexcl
> Protein 3
>
> [ atoms ]
> ; nr type resnr residue atom cgnr charge mass
> typeB chargeB massB
> 1 CH3 1 ASPH CB 2 0 15.035
> DUM 0 15.035
> 2 C 1 ASPH CG 3 0.33 12.011
> DUM 0 12.011
> 3 O 1 ASPH OD1 3 -0.45 15.9994
> DUM 0 15.9994
> 4 OA 1 ASPH OD2 3 -0.288 15.9994
> DUM 0 15.9994
> 5 H 1 ASPH HD2 3 0.408 1.008
> DUM 0 1.008
> [ bonds ]
> ; ai aj funct c0 c1 c2 c3
> 1 2 2 gb_27
> 2 3 2 gb_5
> 2 4 2 gb_13
> 4 5 2 gb_1
>
> [ pairs ]
> ; ai aj funct c0 c1 c2 c3
> 1 5 1
>
> [ angles ]
> ; ai aj ak funct c0 c1
> c2 c3
> 1 2 3 2 ga_30
> 1 2 4 2 ga_19
> 3 2 4 2 ga_33
> 2 4 5 2 ga_12
>
> [ dihedrals ]
> ; ai aj ak al funct c0 c1
> c2 c3 c4 c5
> 1 2 4 5 1 gd_12
>
> [ dihedrals ]
> ; ai aj ak al funct c0 c1
> c2 c3
> 1 4 3 2 2 gi_1
>
> and put the molecule in an empty box at the size of 5.32 nm^3 and
> simulated using the below parameters:
>
> integrator = sd
> ld_seed = -1
> dt = 0.002 ; ps !
> nsteps = 550000 ; total 1.1ns.
>
> ; OPTIONS FOR BONDS
> ; Constrain control
> constraints = all-bonds
> ; Do not constrain the start configuration
> continuation = no
> ; Type of constraint algorithm
> constraint-algorithm = lincs
>
> ; OUTPUT CONTROL OPTIONS
> ; Output frequency for coords (x), velocities (v) and forces (f)
> nstxout = 10000
> nstvout = 10000
> nstfout = 10000
> ; Output frequency and precision for xtc file
> nstxtcout = 500
> xtc-precision = 1000
> ; Energy monitoring
> energygrps = protein
> nstenergy = 500
>
> ; NEIGHBORSEARCHING PARAMETERS
> ; nblist update frequency
> nstlist = 5
> ; ns algorithm (simple or grid)
> ns-type = Grid
> ; nblist cut-off
> rlist = 0.8
>
> ; OPTIONS FOR ELECTROSTATICS AND VDW
> ; Method for doing electrostatics
> coulombtype = Generalized-Reaction-Field
> rcoulomb = 1.4
> epsilon_rf = 62
> ; Method for doing Van der Waals
> vdw-type = Cut-off
> ; cut-off lengths
> rvdw-switch = 0
> rvdw = 1.4
> ; Apply long range dispersion corrections for Energy and Pressure
> DispCorr = no
>
> ; OPTIONS FOR WEAK COUPLING ALGORITHMS
> ; Temperature coupling
> tcoupl = no
> ; Groups to couple separately, time constant (ps) and reference
> temperature (K)
> tc-grps = system
> tau-t = 0.2
> ref-t = 298
> ; Pressure coupling
> Pcoupl = no
> ; Generate velocites is on at 290K - do not get velocity from gro file.
> gen_vel = yes
> gen_temp = 290
> gen-seed = -1
>
> ; Free energy control stuff
> free-energy = yes
> init-lambda = 0 ; Topology A (lambda=0) to topology B
> (lambda=1).
> delta-lambda = 0
> sc-alpha = 0.5 ; soft core potantial
> sc-power = 1
> sc-sigma = 0.3
>
> ; Center of mass control
> nstcomm = 1000
> ; Periodic boundary conditions
> pbc = xyz
> ; Mode for center of mass motion removal
> comm_mode = Linear
> ; Groups for center of mass motion removal
> comm-grps = system
>
> The FE value out of the in vacuum simulation should be 0, as I don't
> have any interactions larger then 1-4, but what I get is -0.0956224
> kJ/mol. Some more technical data - I have done 40 (delta_lambda =
> 0.025), for each lambda value I had calculated the dg/dl for the last
> 1ns (out of 1.1ns). Any idea why this happen?
>
> All the best,
> Itamar
>
> --
>
>
> "In theory, there is no difference between theory and practice. But, in practice, there is." - Jan L.A. van de Snepscheut
>
> ===========================================
> | Itamar Kass, Ph.D.
> | Postdoctoral Research Fellow
> |
> | Department of Biochemistry and Molecular Biology
> | Building 77 Clayton Campus
> | Wellington Road
> | Monash University,
> | Victoria 3800
> | Australia
> |
> | Tel: +61 3 9902 9376
> | Fax: +61 3 9902 9500
> | E-mail: Itamar.Kass at monash.edu
> ============================================
>
> --
> gmx-users mailing list gmx-users at gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at http://www.gromacs.org/search before posting!
> 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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20100803/10bf702b/attachment.html>
More information about the gromacs.org_gmx-users
mailing list