Hello,
i'm interested in change of free energy of protein + ligand complex upon
mutation. While there's some information and tutorials on decoupling of
small molecules in water, i found little information on set up for free
energy upon mutation, with a dated toluene -> p-cresol tutorial by
Gilles Pieffet and Alan E. Mark being best hit i came across. I'm having
problems with my topology as well simulation parameters.
My work-flow:
I take a pdb reference protein with ligand and mutate it (Modeller) to
represent a specific genotype. I then proceed to do EM with Gromacs on
the structure (using Antechamber to parametrize the ligand, for which
reason i then use amber99sb-ildn FF in the simulation). For the
resulting minimized structure i want to proceed with estimation of free
energy change upon L -> V mutation. In the topology, for B state of the
residue in question, i introduce 3 H atoms, change 6 H to dummies, do
other atom substitutions/charge adaptations to valine and update all the
following atom indices (topology snippet below).
====================
; residue 76 LEU rtp LEU q 0.0
1208 N 76 LEU N 1208 -0.4157 14.01 ; qtot
0.5843
1209 H 76 LEU H 1209 0.2719 1.008 ; qtot
0.8562
1210 CT 76 LEU CA 1210 -0.0518 12.01
CT -0.0875 12.01 ; qtot 0.7687
1211 H1 76 LEU HA 1211 0.0922 1.008
H1 0.0969 1.008 ; qtot 0.8656
1212 CT 76 LEU CB 1212 -0.1102 12.01
CT 0.2985 12.01 ; qtot 1.164
1213 HC 76 LEU HB1 1213 0.0457 1.008
HC -0.0297 1.008 ; qtot 1.134
1214 HC 76 LEU HB2 1214 0.0457 1.008
CT -0.3192 12.01 ; qtot 0.8152
1215 DUM 76 LEU DUM 1215 0 1.008
HC 0.0791 1.008 ; qtot 0.8943
1216 DUM 76 LEU DUM 1216 0 1.008
HC 0.0791 1.008 ; qtot 0.9734
1217 DUM 76 LEU DUM 1217 0 1.008
HC 0.0791 1.008 ; qtot 1.052
1218 CT 76 LEU CG 1218 0.3531 12.01
CT -0.3192 12.01 ; qtot 0.7333
1219 HC 76 LEU HG 1219 -0.0361 1.008
HC 0.0791 1.008 ; qtot 0.8124
1220 CT 76 LEU CD1 1220 -0.4121 12.01
HC 0.0791 1.008 ; qtot 0.8915
1221 HC 76 LEU HD11 1221 0.1 1.008
DUM 0 1.008 ; qtot 0.8915
1222 HC 76 LEU HD12 1222 0.1 1.008
DUM 0 1.008 ; qtot 0.8915
1223 HC 76 LEU HD13 1223 0.1 1.008
DUM 0 1.008 ; qtot 0.8915
1224 CT 76 LEU CD2 1224 -0.4121 12.01
HC 0.0791 1.008 ; qtot 0.9706
1225 HC 76 LEU HD21 1225 0.1 1.008
DUM 0 1.008 ; qtot 0.9706
1226 HC 76 LEU HD22 1226 0.1 1.008
DUM 0 1.008 ; qtot 0.9706
1227 HC 76 LEU HD23 1227 0.1 1.008
DUM 0 1.008 ; qtot 0.9706
1228 C 76 LEU C 1228 0.5973 12.01 ; qtot
1.568
1229 O 76 LEU O 1229 -0.5679 16 ; qtot 1
====================
Added bonds/angles/dihedrals snippets:
====================
[ bonds ]
; ai aj funct c0 c1 c2 c3
1214 1215 1
1214 1216 1
1214 1217 1
[ angles ]
; ai aj ak funct c0 c1 c2 c3
1212 1214 1215 1
1212 1214 1216 1
1212 1214 1217 1
1215 1214 1216 1
1215 1214 1217 1
1216 1214 1217 1
[ dihedrals ]
; ai aj ak al funct c0 c1 c2
c3 c4 c5
1210 1212 1214 1215 9
1210 1212 1214 1216 9
1210 1212 1214 1217 9
1213 1212 1214 1215 9
1213 1212 1214 1216 9
1213 1212 1214 1217 9
1221 1212 1214 1215 9
1221 1212 1214 1216 9
====================
In the structure file the following dummies are added in the vicinity of
HB2 (-> CT) (LEU residue snipplet):
====================
76LEU N 1208 3.836 3.954 3.263
76LEU H 1209 3.914 3.962 3.328
76LEU CA 1210 3.765 4.077 3.231
76LEU HA 1211 3.672 4.053 3.180
76LEU CB 1212 3.730 4.151 3.361
76LEU HB1 1213 3.667 4.235 3.336
76LEU HB2 1214 3.822 4.191 3.405
76LEU DUM 1215 3.910 4.191 3.379
76LEU DUM 1216 3.910 4.231 3.428
76LEU DUM 1217 3.910 4.131 3.451
76LEU CG 1218 3.661 4.064 3.468
76LEU HG 1219 3.729 3.988 3.506
76LEU CD1 1220 3.618 4.153 3.583
76LEU HD11 1221 3.566 4.094 3.658
76LEU HD12 1222 3.705 4.199 3.629
76LEU HD13 1223 3.551 4.230 3.547
76LEU CD2 1224 3.534 3.997 3.416
76LEU HD21 1225 3.484 3.946 3.498
76LEU HD22 1226 3.470 4.074 3.373
76LEU HD23 1227 3.560 3.924 3.340
76LEU C 1228 3.849 4.164 3.138
76LEU O 1229 3.968 4.180 3.162
====================
To account for the interactions with dummy atoms, i include my bonded
itp file:
====================
[ bondtypes ]
; i j func b0 kb
DUM H 1 0.1080 307105.6
DUM HC 1 0.1080 307105.6
DUM C 1 0.1080 307105.6
[ angletypes ]
; i j k func th0 cth
DUM HC DUM 1 109.500 292.880 ;
DUM HC CT 1 109.500 418.400 ;
CT DUM DUM 1 109.500 418.400 ;
[ dihedraltypes ]
; i j k l func
CT CT HC DUM 9 0.0 0.62760 3 ;
HC CT CT DUM 9 0.0 0.62760 3 ;
HC CT HC DUM 9 0.0 0.62760 3 ;
====================
I'm basically using em mdp file from J. Lemkul's Methane in Water
tutorial, with couple- parameters removed, as i'm using the B-state
parameters in topology instead (am i correct here?) and constraints
applied to all-bonds:
====================
; Run control
integrator = steep
nsteps = 5000
; EM criteria and other stuff
emtol = 100
emstep = 0.01
niter = 20
nbfgscorr = 10
; Output control
nstlog = 1
nstenergy = 1
; Neighborsearching and short-range nonbonded interactions
nstlist = 1
ns_type = grid
pbc = xyz
rlist = 1.0
; Electrostatics
coulombtype = PME
rcoulomb = 1.0
; van der Waals
vdw-type = switch
rvdw-switch = 0.8
rvdw = 0.9
; Apply long range dispersion corrections for Energy and Pressure
DispCorr = EnerPres
; Spacing for the PME/PPPM FFT grid
fourierspacing = 0.12
; EWALD/PME/PPPM parameters
pme_order = 6
ewald_rtol = 1e-06
epsilon_surface = 0
optimize_fft = no
; Temperature and pressure coupling are off during EM
tcoupl = no
pcoupl = no
; Free energy control stuff
free_energy = yes
init_lambda = 0.0
delta_lambda = 0
foreign_lambda = 0.05
sc-alpha = 0.5
sc-power = 1.0
sc-sigma = 0.3
nstdhdl = 10
; Generate velocities to start
gen_vel = no
; options for bonds
constraints = all-bonds
; Type of constraint algorithm
constraint-algorithm = lincs
; Do not constrain the starting configuration
continuation = no
; Highest order in the expansion of the constraint coupling matrix
lincs-order = 12
====================
Output of the "grompp -f em_steep.mdp -c ../em_stateB.gro -p
../Complex_stateB.top -o min_test.tpr":
====================
NOTE 1 [file em_steep.mdp]:
init-lambda is deprecated for setting lambda state (except for slow
growth). Use init-lambda-state instead.
NOTE 2 [file em_steep.mdp]:
With PME there is a minor soft core effect present at the cut-off,
proportional to (LJsigma/rcoulomb)^6. This could have a minor effect on
energy conservation, but usually other effects dominate. With a common
sigma value of 0.34 nm the fraction of the particle-particle potential at
the cut-off at lambda=0.5 is around 6.4e-05, while ewald-rtol is 1.0e-06.
NOTE 3 [file em_steep.mdp]:
The switch/shift interaction settings are just for compatibility; you
will get better performance from applying potential modifiers to your
interactions!
Generated 3240 of the 3240 non-bonded parameter combinations
Generating 1-4 interactions: fudge = 0.5
Generated 3240 of the 3240 1-4 parameter combinations
WARNING 1 [file Protein_Protein_chain_A_stateB.itp, line 2926]:
No default Bond types for perturbed atoms, using normal values
WARNING 2 [file Protein_Protein_chain_A_stateB.itp, line 13733]:
Some parameters for bonded interaction involving perturbed atoms are
specified explicitly in state A, but not B - copying A to B
ERROR 1 [file Protein_Protein_chain_A_stateB.itp, line 13762]:
Cannot automatically perturb a torsion with multiple terms to different
form.
Please specify perturbed parameters manually for this torsion in your
topology!
ERROR 2 [file Protein_Protein_chain_A_stateB.itp, line 13763]:
Cannot automatically perturb a torsion with multiple terms to different
form.
Please specify perturbed parameters manually for this torsion in your
topology!
Excluding 3 bonded neighbours molecule type 'Protein_chain_A'
Excluding 3 bonded neighbours molecule type 'Protein_chain_B'
Excluding 3 bonded neighbours molecule type 'SQV'
Excluding 2 bonded neighbours molecule type 'SOL'
Excluding 1 bonded neighbours molecule type 'CL'
There were 3 notes
There were 2 warnings
-------------------------------------------------------
Program grompp, VERSION 4.6.5
Source code file: /X/gromacs-4.6.5_src/src/kernel/grompp.c, line: 1610
Fatal error:
There were 2 errors in input file(s)
For more information and tips for troubleshooting, please check the GROMACS
website at http://www.gromacs.org/Documentation/Errors
-------------------------------------------------------
turning all bonds into constraints...
turning all bonds into constraints...
turning all bonds into constraints...
turning all bonds into constraints...
turning all bonds into constraints...
====================
Below given the lines 2926, 13733, 13762, 13763 from the topology file
from the dihedrals section:
====================
1210 1225 1
...
1228 1210 1212 1218 9 torsion_LEU_C_CA_CB_CG_mult1
...
1210 1212 1218 1220 9
1210 1212 1218 1224 9
====================
Regarding errors 1-2: These being groups of CT atoms, dihedral types for
which are defined in ambers ffbonded.itp, why is there an error or what
additional torsion parameters are expected?
Regarding warning 2: The specific updated torsion for LEU will not apply
for VAL in state B anymore - should then this specification just be
dropped from the topology for that particular position?
I would also be very pleased to hear comments on my approach for the
given system in general.
Thank you
tomas
--
Gromacs Users mailing list
* Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a
mail to [email protected].
