Solvation free energy calculation

User discussions
1

GROMACS version: 2024.2
GROMACS modification: No

Hi everyone,

good morning, thank you for your interest. I still have to get used to thermodynamic integration to compute the free energy of solvation.
I am writing to ask:

  1. Shoud sc-coul = yes if I want to decouple both vdw and q altogether, per each lambda run? In fact, from the manual I did not understand what does sc-coul mean.

  2. Can I compute the free energy of solvation via multiple lambda runs where I decouple separately with the following values of parameters?

couple_lambda0 = q
couple_lambda1 = none
vdw_lambda = 0.0 0.5 …… 0.0
coul lambdas = 0.0 0.0 …… 0.0
sc-coul = no (or yes, still?)

and then, at the same lambda, doing another run:

couple_lambda0 = vdw
couple_lambda1 = none
vdw_lambdas = 0.0 0.0 …… 0.0
coul_lambda = 0.0 0.5 …… 0.0
sc-coul = no (or yes, still?)

All the other parameters are not signaled since I know what to put. My doubts are on those above. By the way, is it more advisable in case 2) to put couple_lambda0 = vdwq and couple_lambda1 = vdw (or q, respectively) ? I know it is equivalent in theory, but I would like people pointing out which is the most accurate, and which the fastest (I think the one I extensively reported is the fastest, the second hereby mentioned the more accurate).

Much obliged to everyone answering, thank you for your time!

Sincerely,
Jacopo

2

With

couple_lambda0 = q
couple_lambda1 = none

you would definitely want to use soft-core for coulomb. But unless I’m completely mistaken sc-coul will not make any difference in this case. Just make sure that you have set the sc-alpha and the other soft-core parameters.

I might be wrong, but I think it’s only when using a setup like

couple_lambda0 = vdwq
couple_lambda1 = none

(or when using modifications in the topology) that sc-alpha makes a difference, i.e., if the soft-core parameters should be applied to both the vdw and q transformations.

I would recommend doing two sets of sequences of runs:

couple_lambda0 = vdwq
couple_lambda1 = vdw

and from the lambda 1 state output in that set of runs you run:

couple_lambda0 = vdw
couple_lambda1 = none

In theory, it is most efficient to run the respective lambda state simulations in each set in sequence, since you will then start each lambda state simulation from a state that is similar. But if you want to run very many simulations in parallel (better hardware utilisation) you will need to discard more time at the beginning of each simulation since it will take longer to reach equilibrium.

In the first of these sets you can set sc-alpha = 0 since no soft-core modifications are needed when there are still vdw interactions. In the second you should set the proper sc parameters. sc-coul should not be necessary in either of these cases.

3

Dear Dr. Magnus,

good morning, thank you for your interest. I’ve tried:

couple_lambda0 = q
couple_lambda1 = none
sc-coul = no (yes)

but the results are different in dhdl.xvg files, thus:

  • we should go for sc-coul = yes in that scenario;
  • also in the third code chunk you posted there is a difference, but naturally the sc-coul = no in that very case is perfectly fine.

Thank you for the insights, imagine that I would have engaged a vdwq —> q and then q —> none process, instead. In theory it should be equivalent but, given that in practice one has to apply soft core potentials, it is not. For instance, I understand the advantages you explained in your choice: both avoiding sc-coul = yes anywhere, as indeed you pointed out (in particular, this is implicitly stated when saying sc-alpha = 0 in your second-last sentence), and not employing soft core potentials at all in the first step, neither in VdW nor in Coulomb interactions.

Thanks again, wish you a good continuation.

Sincerely,
Jacopo