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.k...@monash.edu > To: gmx-users@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.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
-- 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
