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