[gmx-users] Re: Free energy calculation question

[Fabian Casteblanco]

Hello Justin,

You mean that only for vdW decoupling, you would need to use soft-core
potentials?  I had soft -core potentials on for decoupling the
electrostatic interactions (see below).   What would I use in its
place?  Thanks again for your help!


;Production MD
;------------------------------------------------------------------------------------
include		=-I/mphase/users2/fabian/CGenff
title		=CGenFF Lov/Eth Solution MD run

;PARAMETERS - describing what to do, when to stop and what to save
;------------------------------------------------------------------------------------

;Run parameters
integrator		=sd		;leap-frog integrator
ld_seed			=-1
nsteps			=500000		;2*500000=1000ps, 1 ns
dt			=0.002		;2 fs
nstcomm                 = 100		;*** - frequency for center of mass
motion removal

;Output control
nstxout			=1000		;save coordinates every 2ps
nstvout			=1000		;save velocities every 2 ps
nstenergy		=1000		;save energies every 2 ps
nstlog			=1000		;update log file every 2 ps
nstxtcout		=1000		;xtc compressed trajectory output every 2 ps
xtc-precision           =1000		;*** - precision to write to xtc trajectory

;Bond Parameters
continuation		=yes		;Restarting after NPT
constraint_algorithm	=lincs	    	;holonomic constraints
constraints		=all-bonds  	;all bonds (even heavy atom-H bonds) constr;ained
lincs_iter		=1	    	;accuracy of LINCS
lincs_order		=12	    	;also related to accuracy

;Neighborhood searching
ns_type			=grid		;search neighboring grid cells
nstlist			=5		;#steps.  5*0.002 ps = 5* 2 fs = 10 fs - Frequency to
update the neighbor list (and the long-range forces, when 					;using
twin-range cut-off’s). When this is 0, the neighbor list is made only
once.
rlist			=1.1		;short-range neighborlist cutoff (in nm)
rcoulomb		=1.1		;short-range electrostatic cutoff (in nm)
pbc	  		=xyz		; 3-D PBC

;Electrostatics	
coulombtype		=PME		;Particle Mesh Ewald for long-range electrostat;ics
pme_order		=4		;cubic interpolation
fourierspacing		=0.16		;grid spacing for FFT

; van der Waals
vdwtype	    		=Shift     	;Van der Waals for CHARMM
rvdw_switch 		=0.8
rvdw	    		=1.0    	;Short-range Van der Waals cut-off

;Dispersion correction
DispCorr    =EnerPres	;account for cut-off vdW scheme

;Temperature coupling is on
tcoupl	     		=V-rescale	;modified Berendsen thermostat
tc-grps	     		=SYSTEM  	;two coupling groups - more accurate
tau_t	     		=0.1		;time constant, in ps
ref_t	     		=298 		;reference temperature, on for each group, in K

;Pressure coupling is on
pcoupl	  =Parrinello-Rahman	;Pressure coupling on in NPT
pcoupltype			=isotropic	   ;uniform scaling of box vect;ors
tau_p				=2.0		   ;time constant, in ps
ref_p				=1.0		   ;reference pressure, in bar
compressibility			=4.5e-5		   ;isothermal compr of H2O, ba;r^(-1)

; Free energy control stuff
free_energy             = yes		;*** - Indicates we are doing a free
energy calculation, and that interpolation between the A and B states
of the 						;chosen molecule (defined below) will occur.
init_lambda             = 0.0		;*** - Value of λ
delta_lambda            = 0		;*** - The value of λ can be incremented
by some amount per timestep (i.e., δλ/δt) in a technique called "slow
						;growth." This method can have significant errors associated
with it, and thus we will make no time-dependent 						;changes to our
λ values.
foreign_lambda          = 0.05		;*** - Additional values of λ for
which ΔH will be written to dhdl.xvg (with frequency nstdhdl). The
						;configurations generated in the trajectory at λ = init_lambda
will have ΔH calculated for these same 						;configurations at all
values of λ = foreign_lambda.
sc-alpha                = 0.5		;*** - the soft-core parameter, a value
of 0 results in linear interpolation of the LJ and Coulomb
					;interactions
sc-power                = 1.0		;*** - the power for lambda in the
soft-core function, only the values 1 and 2 are supported
sc-sigma                = 0.3 		;*** - the soft-core sigma for
particles which have a C6 or C12 parameter equal to zero or a sigma
					;smaller than sc_sigma
couple-moltype          = LOV	  	;*** - name of moleculetype to decouple
couple-lambda0          = vdw-q		;*** - all interactions are on at lambda=0
couple-lambda1          = vdw     	;*** - only vdw interactions are on
at lambda=1
couple-intramol         = no		;*** - Do not decouple intramolecular
interactions. That is, the λ factor is applied to only solute-solvent
						;nonbonded interactions and not solute-solute nonbonded
interactions.
nstdhdl                 = 100		;*** - The frequency with which ∂H/∂λ
and ΔH are written to dhdl.xvg. A value of 100 would probably suffice,
since 						;the resulting values will be highly correlated and the
files will get very large. You may wish to increase this 						;value
to 100 for your own work.

;Velocity generation
gen_vel	  =no	;Velocity generation is off

;END

On Wed, Aug 31, 2011 at 11:24 AM, Fabian Casteblanco
<fabian.casteblanco at gmail.com> wrote:
> Hello Justin,
>
> I'm calculating the free energy of a drug in an alcohol solvent.  I
> have a question referring to your free energy tutorial.  You mentioned
> that decoupling of electrostatic interactions is linear and decoupling
> of vdW can vary.  Is this true for your case of methanol in water or
> for all cases?  When I ran my system of drug in ethanol solvent, I got
> a non linear dG for both electrostatic and vdW.  Also, is there no
> need to find dG of cav ( the free energy required to form the solute
> cavity within the solvent) ?  I have attached some graphs.
>
> Thanks for your help.  Your tutorial was extremely useful.
>
> --
> Best regards,
>
> Fabian F. Casteblanco
> Rutgers University --
> Chemical Engineering PhD Student
> C: +908 917 0723
> E:  fabian.casteblanco at gmail.com
>



-- 
Best regards,

Fabian F. Casteblanco
Rutgers University --
Chemical Engineering PhD Student
C: +908 917 0723
E:  fabian.casteblanco at gmail.com