Free Energy Calculation of Tryptophan

User discussions
1

Hello everyone, i’m a beginnner in gromacs, and this is my first post in this forum and im looking forward to connect with a new community here ! I installed the GROMACS version: 2025.2 with normal settings i would say.I wanted to do a free energy calculation of Tryptophan guided by the great Tutorials from Lemkul Lab.

So finally i got the point where i got results for my simulation and it looks like this:

Final results in kJ/mol:

point 0 - 1, DG 3.05 +/- 0.02
point 1 - 2, DG 3.01 +/- 0.03
point 2 - 3, DG 2.96 +/- 0.03
point 3 - 4, DG 2.91 +/- 0.02
point 4 - 5, DG 2.71 +/- 0.05
point 5 - 6, DG 2.50 +/- 0.04
point 6 - 7, DG 2.42 +/- 0.07
point 7 - 8, DG 2.27 +/- 0.05
point 8 - 9, DG 2.00 +/- 0.10
point 9 - 10, DG 1.61 +/- 0.08
point 10 - 11, DG 1.23 +/- 0.06
point 11 - 12, DG 0.74 +/- 0.10
point 12 - 13, DG -0.18 +/- 0.15
point 13 - 14, DG -1.25 +/- 0.24
point 14 - 15, DG -3.01 +/- 0.24
point 15 - 16, DG -6.43 +/- 0.05
point 16 - 17, DG -7.72 +/- 0.10
point 17 - 18, DG -6.02 +/- 0.10
point 18 - 19, DG -3.55 +/- 0.06
point 19 - 20, DG -1.15 +/- 0.03

total 0 - 20, DG -1.88 +/- 0.50

I used the charmm forcefield with just a Trp residue in a waterbox. The research by sandler et al and other researchers says it has to be at least -25 KJ/mol. I created my input files (by solvation builder), so my topol.top files looks like this:

; Include forcefield parameters
#include “toppar/forcefield.itp”
#include “toppar/PROA.itp”
#include “toppar/POT.itp”
#include “toppar/CLA.itp”
#include “toppar/TIP3.itp”

[ system ]
; Name
Title

[ molecules ]
; Compound #mols
PROA 1
POT 2
CLA 2
TIP3 825

Has anyone an idea what might went wrong or where i have to check ?

2

Iam answering myself, maybe someone can learn from that. I was looking at my mdp settings and i noticed, that i just did init_lambda settings for vdw like in the Lemkul Tutorial. Now that i don’t have a system like methane anymore there should be more interactions. so iam trying it again with the lambda settings for coulomb interactions like this:

; Vectors of lambda specified here
; Each combination is an index that is retrieved from init_lambda_state for each simulation
; init_lambda_state 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
vdw_lambdas = 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00
coul_lambdas =0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00
; We are not transforming any bonded or restrained interactions
bonded_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 0.00 0.00 0.00 0.00 0.00 0.00 0.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 0.00 0.00 0.00 0.00 0.00 0.00 0.00
; Masses are not changing (particle identities are the same at lambda = 0 and lambda = 1)
mass_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 0.00 0.00 0.00 0.00 0.00 0.00 0.00
; Not doing simulated temperting here
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 0.00 0.00 0.00 0.00 0.00 0.00 0.00