Mutating residue, free energy

User discussions
1

GROMACS version: Any version
GROMACS modification: No

Hi Everyone,

I was following this tutorial for the free energy calculation relative to a mutation: pmx_tutorial_template . I am wondering whether one can do the same, adopting a topology file with the same modifications needed in the tutorial itself (simply, it’s three commands to be performed, up to the ‘generate_hybrid_topology.py’ one), but without slow or fast growth methods.

In particular, something completely analogous to this one: Justin_Lemkul_tutorial . I am attaching the few lines of the mdp file I’d like to adopt:

; Free energy control stuff
free_energy              = yes
init_lambda_state        = 1
delta_lambda             = 0
calc_lambda_neighbors    = 1        ; only immediate neighboring windows
; init_lambda_state        0       1       2       3       4       5       6       7       8       9       10      11      12      13
vdw_lambdas              = 0.00    0.00922 0.04794 0.11505 0.20634 0.31608 0.43738 0.56262 0.68392 0.79366 0.88495 0.95206 0.99078 1.00
coul_lambdas             = 0.00    0.00922 0.04794 0.11505 0.20634 0.31608 0.43738 0.56262 0.68392 0.79366 0.88495 0.95206 0.99078 1.00
; We are not transforming any bonded or restrained interactions
bonded_lambdas             = 0.00    0.00922 0.04794 0.11505 0.20634 0.31608 0.43738 0.56262 0.68392 0.79366 0.88495 0.95206 0.99078 1.00
restraint_lambdas        = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
mass_lambdas             = 0.00    0.00922 0.04794 0.11505 0.20634 0.31608 0.43738 0.56262 0.68392 0.79366 0.88495 0.95206 0.99078 1.00
temperature_lambdas      = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
sc-alpha                 = 0 ; no need as we are not vanishing, but mutating
sc-sigma                = 0.3
nstdhdl                   = 10

Where you can infer it’s an MD at a certain fixed lambda value that will be carried out, I am just providing an example.

I am also wondering why I have a note for a W2F residue when using grompp with the former mdp file (both the grompp and mdrun command work), which is:

NOTE 1 [file topol.top, line 47]:
The bond in molecule-type Protein between atoms 846 DCD1 and 850 DCE1 has an estimated oscillational period of 7.1e-03 ps, which is less than 10 times the time step of 1.0e-03 ps. Maybe you forgot to change the constraints mdp option.

As you can see, I have dt = 0.001 and - being a carbon-carbon bond - I don’t think that putting:
constraints = all-bonds
is a good idea, but a way too large approximation. I wonder why a carbon-carbon bond raises this note for the phenylalanine to which the tryptophan is mutated, while all the other phenylalanine residues in the protein did never raise this note. I am sharing below the topol.top and relative .itp, .gro files, as well as the link to the force field adopted (which includes, of course, the mutations: it is the ‘mutff45.tar.gz’ one, amber99sbmut.ff directory), maybe this can help. Thank you for your time, I appreciate.

Bests,
Jacopo

top_itp_gro
force_field_with_mutations

2

I must admit I haven’t worked that much with mutations, but I think your mdp settings look OK.

I think the reason why you are getting the note about the bond oscillation period is that your dummy atoms have a mass of 1, so a bond between them will have a long oscillation period.

3

Dear Dr. Lundborg,

your answer about the note is reassuring to me, as the modifications are smoothly made thanks to pmx python package (not “by hand”). Maybe I should plug sc-alpha, too, as I am going from a W to a F, meaning that part of the atoms are vanishing from vdw-q to null, which I suppose could cause numerical inaccuracy. I appreciate your immediate response: thank you.

Sincerely,
Jacopo

4

Yes, I think it would be wise to use soft-core interactions in that case. And since you are turning on/off Coulomb and van der Waals interactions at the same, you should also you sc-coul.

1 Like