Pmx RBFE gives oddly large deltaG

pmx
1

Hello,
I’m new to pmx still and have been trying to run RBFE for a macrocycle host and small molecule guests: glycine anhydride morphing into glycolide (hence NH → O).
I’d like to believe the ligandHybrid step worked as the merged.itp file has:
[ atoms ]
; nr type resnr residue atom cgnr charge mass typeB chargeB massB
(…)
3 n 1 UNK N3 3 -0.574900 14.0100 os -0.395900 16.0000
4 n 1 UNK N4 4 -0.574900 14.0100 os -0.395900 16.0000
(…)
9 hn 1 UNK H9 9 0.333500 1.0080 DUM_hn 0.000000 1.0080
10 hn 1 UNK H10 10 0.333500 1.0080 DUM_hn 0.000000 1.0080
(…)

I then followed the standard md simulation steps, solvated with spc water model, minimised, equilibrated with temp (with single t coupling group) and pressure (Parinello-Rahman, isotropic), and run the production run over 100 ps. I kept LINCS constraints with order 6 (I wasn’t sure of the effect of changing this) and with FEP:

; Free energy control stuff
free-energy = yes
init-lambda = 0
delta-lambda = 4e-6
sc-function = gapsys
sc-alpha = 0.3
sc-sigma = 0.25
sc-power = 1
sc-coul = yes

The minimisation was a bit temperamental, requiring a few steep and then stochastic descent steps but the remainder run fine with no excessive variations in temperature or pressure.
I run the entire setup in forward and reverse directions:
init-lambda = 0/1

I was expecting the free energy to not be overly large and yet the pmx analyse produced:
Number of forward (0->1) trajectories: 1
Number of reverse (1->0) trajectories: 1
Temperature : 298.15 K

         Crooks Gaussian Intersection

CGI: Forward Gauss mean = 120.48 kJ/mol std = 0.00 kJ/mol
CGI: Reverse Gauss mean = 120.00 kJ/mol std = 0.00 kJ/mol

Gaussians too close for intersection calculation
→ Taking difference of mean values
CGI: dG = 120.24 kJ/mol
CGI: Std Err (bootstrap:parametric) = 0.00 kJ/mol
CGI: Std Err (bootstrap) = 0.00 kJ/mol
Forward: gaussian quality = nan
—> KS-Test Ok
Reverse: gaussian quality = nan
—> KS-Test Ok

The remainder of estimators give the same dG value.
Firstly, is the proximity of Gaussians a problem?
Secondly, what could be the source of such a large value?

Great thank you in advance,
Katarzyna Zator

2

Several observations:

  1. you only ran one transition in the forward and one in the reverse direction. This will not suffice. I recommend to increase the number of transitions to the order of ~100. See chapter 6 in this book chapter for the description of the non-eq protocol: http://www3.mpibpc.mpg.de/groups/de_groot/pdf/Gapsys_MolMod_2015.pdf
  2. so far you have calculated only dG of one leg in the thermodynamic cycle. To get ddG you also need the other leg.

Vytas