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.k...@monash.edu
============================================
--
gmx-users mailing list gmx-users@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-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/mailing_lists/users.php
