AWH very large free energy barriers

User discussions
, , ,
1

GROMACS version: 2023.3
GROMACS modification: No
I have been trying to run a Free energy simulation using AWH and no matter what I set some of the values too I keep getting the same error:

-------------------------------------------------------
Program:     gmx mdrun, version 2023.3
Source file: src/gromacs/applied_forces/awh/pointstate.h (line 355)
Function:    gmx::PointState::updateFreeEnergy(const gmx::BiasParams&, double)::<lambda()>

Assertion failed:
Condition: std::abs(freeEnergy_) < detail::c_largePositiveExponent
Very large free energy differences or badly normalized free energy in AWH
update.

For more information and tips for troubleshooting, please check the GROMACS
website at http://www.gromacs.org/Documentation/Errors
-------------------------------------------------------

This is the mdp file:

; Free energy control stuff
free_energy        = yes
init_lambda-state        = 0
delta_lambda       = 0
calc_lambda_neighbors    = -1  
couple-moltype     = Protein_chain_A
couple-lambda0     = vdw-q
couple-lambda1     = none
couple-intramol    = yes
fep-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
sc-alpha                 = 0.5
sc-coul                  = no       ; linear interpolation of Coulomb (none in this case)
sc-power                 = 1
sc-sigma                 = 0.3


; AWH stuff
awh = yes
awh-potential = umbrella
awh-nstout = 50000
awh-nbias = 1
awh-nstsample = 100
awh-nsamples-update = 100
awh1-error-init = 10
awh1-equilibrate-histogram = no
awh1-target = constant
awh1-growth = exp-linear
awh1-ndim = 1
awh1-dim1-start = 0.0
awh1-dim1-end = 1.0
awh1-dim1-coord-provider = fep-lambda
awh1-dim1-diffusion = 0.0001

Have tired making the awh1-dim1-diffusion and awh1-error-init options bigger and smaller, default, slightly altered default. I have tried setting awh-dim1-end to 0.75, 0.5, and 0.45 nothing as worked. I think I just need a professional opinion.

2

Does the molecule that you are decoupling have charges? If yes, than you need to set sc-coul=yes.

But the user should not get an assertion failure. Maybe we should change this to a normal error.

3

When you say charges you mean like when I create the topology file and it tells me whether or the not the systems charge is positive, negative, or neutral? Well in that case yes. If I did set sc-coul=yes and it works would I have to re-run the enitre simulation? For accuracy sake, because of the mistake?

4

No, I meant to ask if the molecule you are decoupling has atoms with partial charges or the charge of every atom in the molecule is zero.

5

Yeah it definitely has partial charges, it’s protein.

6

Then you need to use soft-core also for Coulomb interactions.

7

I filed an issue and created a fix that makes mdrun exit with a proper error message.

8

I tried just creating another TPR file with with that option set to yes and started the simulation again using it, and I just got the same error. So I’m going to completely restart the simulation with this option set to yes and see where that takes me. I’ll stay posted.

9

So I set sc-coul=yes and restarted the entire simulation but I’m still getting the same error, it was pushed back a little bit. Originally I was getting the error at around 15ns now I’m getting it at around17ns, but that could also be because I messed with the temperture. Got any other ideas?

10

@MagnusL mentioned that the same issue can occur when large molecules are decoupled with couple-intramol=yes and soft-core on. But I don’t see how that could cause issues. Magnus, can you explain?

11

With large molecules, such as proteins, and especially if there are formal charges, and even more so if they charges are balancing each other out, in the molecule, the free energy difference between fully interacting and fully eliminated (which is what you get with couple-intramol = yes) can be huge. But I don’t see any way around it.

12

Yes, I see now that decoupling a large molecule can exceed the limit on the free-energy of 700 kT. But I wonder if one can get reasonable convergence for such cases. With couple-intramol=yes sampling of the decoupled state can get difficult. With couple-intramol=no the cutoff for LJ limits the size of the molecule.

I have now uploaded an MR that changes the assertion failure to an exception with the message: “An AWH free energy difference is larger than 700 kT, which is not supported”.

13

Does that mean that free energy simualtions or alchemical simulations are more meant for ligands and small peptides rather that big protein structures?

14

We didn’t think that such large free energy differences were relevant when developing the AWH method. But with alchemical transformations one get can beyond 700 kT. It is not clear to me how practically useful such calculations are. If they are many use cases, we should look into moving this limit.

There are two separate aspects here. One is if it is scientifically useful to use alchemical transformations with such large free energy differences. Another is whether one can get converged results in reasonable time.

Note that you can move this limit a bit by dividing the lambda interval into parts. You can e.g. compute lambda 0-0.5 and 0.5-1 separately.