Cannot perform free energy calculations due to the solvent (not the solute) possessing CBT bonds?

GROMACS version: 2020.5
GROMACS modification: No

Hi!

I’m attempting to calculate and compare the free energy of solvation of certain models of a solvent using a coarse-grained model. However, one of the solvent models has CBT bonds, causing GROMACS to return “Functiontype CBT dih does not support currentely [sic] free energy calculations”.

Given that the CBT bonded interactions exist within the solvent molecules and not the solute, there is no alteration in the hamiltonian through adjustment of the lambda value, which only affects the coupling of the solvent. It seems that, in theory, this should be possible, but is prevented nonetheless. Is anyone aware of a way around this problem?

Pasted below are the free energy parameters from the .mdp for a randomly picked lambda value; IPO is the solute to be decoupled, and TO is the solvent containing the CBT bonds (not mentioned here).

Thanks in advance!

; Free energy variables
free_energy = yes
init_lambda_state = 1
delta_lambda = 0
calc_lambda_neighbors = 1 ; only immediate neighboring windows
couple-moltype = IPO ; name of moleculetype to decouple
couple-lambda0 = vdw ; only van der Waals interactions
couple-lambda1 = none ; turn off everything, in this case only vdW
couple-intramol = no
; 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 21 22 23 24 25 26 27 28 29
vdw_lambdas = 0.00 0.05 0.10 0.15 0.20 0.25 0.27 0.29 0.31 0.33 0.35 0.37 0.38 0.39 0.41 0.43 0.45 0.47 0.49 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1
coul_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 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.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 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 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 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 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
; Options for the decoupling
sc-alpha = 1.3
sc-coul = no ; linear interpolation of Coulomb (none in this case)
sc-power = 1
sc-sigma = 0.47
nstdhdl = 10

I think this should be fixed in the 2022 release. This should be covered by the fix for: