Solvation free energy calculation issues

[SHervoe-Hansen]

Hi OpenMMTools devs!

I have been trying to determine solvation free energies and more specifically octanol-water partitioning coefficient using the alchemical library from OpenMMTools but somehow manage to get the wrong result consistently - opposite sign of nearly correct magnitudes.

Am I doing something wrong in my approach outlined below? The code shown is for the solute (hexanol) in water with the lambda value being for the fully coupled state:

# Imports
import sys
import os
import numpy as np
import openmm as mm
from openmm import app
from openmm import unit as u
from openmmtools import alchemy

lambdas = {1: {'sterics': 0.0, 'electrostatics': 0.0}, 2: {'sterics': 0.1, 'electrostatics': 0.0}, 3: {'sterics': 0.2, 'electrostatics': 0.0}, 4: {'sterics': 0.3, 'electrostatics': 0.0}, 5: {'sterics': 0.4, 'electrostatics': 0.0}, 6: {'sterics': 0.5, 'electrostatics': 0.0}, 7: {'sterics': 0.6, 'electrostatics': 0.0}, 8: {'sterics': 0.7, 'electrostatics': 0.0}, 9: {'sterics': 0.8, 'electrostatics': 0.0}, 10: {'sterics': 0.9, 'electrostatics': 0.0}, 11: {'sterics': 1.0, 'electrostatics': 0.0}, 12: {'sterics': 1.0, 'electrostatics': 0.2}, 13: {'sterics': 1.0, 'electrostatics': 0.4}, 14: {'sterics': 1.0, 'electrostatics': 0.6}, 15: {'sterics': 1.0, 'electrostatics': 0.8}, 16: {'sterics': 1.0, 'electrostatics': 1.0}}

print('Loading initial configuration and toplogy')
pdb = app.pdbfile.PDBFile('Hexanol_in_water.pdb')
forcefield = app.ForceField('/home/users/eau/Octanol-Water_partitioning/Force_fields/1-hexanol.xml',
                            '/home/users/eau/Octanol-Water_partitioning/Force_fields/1-octanol.xml',
                            '/home/users/eau/Octanol-Water_partitioning/Force_fields/tip3p.xml')

for residue in pdb.topology.residues():
    if residue.name == 'HEX':
        hexanol_atoms = []
        for atom in residue.atoms():
            hexanol_atoms.append(atom.index)
        break

# Creating system
print('Creating OpenMM System')
system = forcefield.createSystem(pdb.topology, nonbondedMethod=app.PME, ewaldErrorTolerance=0.0005,
                                 nonbondedCutoff=1.2*u.nanometers, constraints=app.HBonds, rigidWater=True)

# Pressure-coupling by a Monte Carlo Barostat (NPT)
system.addForce(mm.MonteCarloBarostat(1*u.bar, 298.15*u.kelvin, 25))

# Enable alchemical system
alchemical_region = alchemy.AlchemicalRegion(alchemical_atoms=hexanol_atoms,
                                             annihilate_electrostatics=True,
                                             annihilate_sterics=False)
factory = alchemy.AbsoluteAlchemicalFactory()
alchemical_system = factory.create_alchemical_system(system, alchemical_region)

# Calculating total mass of system
total_mass = u.sum([system.getParticleMass(i) for i in range(system.getNumParticles())])

# Temperature-coupling by leap frog (BAOAB) Langevin integrator (NVT)
integrator = mm.LangevinMiddleIntegrator(298.15*u.kelvin, 1.0/u.picoseconds, 2.0*u.femtoseconds)

# Select platform
platform = mm.Platform.getPlatformByName('CUDA')
properties = {'CudaPrecision': 'mixed', 'CudaDeviceIndex': '0'}

# Create the Simulation object
sim = app.Simulation(pdb.topology, alchemical_system, integrator, platform, properties)

# Set the particle positions
sim.context.setPositions(pdb.positions)

# Set sampling alchemical state
alchemical_state = alchemy.AlchemicalState.from_system(sim.system)
alchemical_state.lambda_sterics = lambdas[16]['sterics']
alchemical_state.lambda_electrostatics = lambdas[16]['electrostatics']

# Minimize the energy
print('Minimizing energy')
sim.minimizeEnergy()

# Draw initial MB velocities
sim.context.setVelocitiesToTemperature(298.15*u.kelvin)

# Equlibrate simulation
print('Equilibrating...')
sim.step(5000000)  # 5000000*2 fs = 10 ns

# Set up the reporters
sim.reporters.append(app.StateDataReporter('output.dat', 500, totalSteps=25000000+5000000,
    time=True, potentialEnergy=True, kineticEnergy=True, temperature=True, volume=True, density=True,
    systemMass=total_mass, remainingTime=True, speed=True, separator='  '))

# Set up trajectory reporter
sim.reporters.append(app.XTCReporter('trajectory.xtc', reportInterval=500, append=False))

# Prepare collection of potential energies
alchemical_reduced_energies = np.zeros(shape=(50000, len(lambdas)))

# Run dynamics and calculate potential energies
print('Running dynamics!')
for alchemical_snapshot in range(50000):
    # Step 1: Sample at a given lambda value
    sim.step(500)

    # Step 2: Calculate energies given native and perturbed Hamiltonians
    for i, (_, l) in enumerate(lambdas.items()):
        alchemical_state.lambda_sterics = l['sterics']
        alchemical_state.lambda_electrostatics = l['electrostatics']
        alchemical_state.apply_to_context(sim.context)
        energy = sim.context.getState(getEnergy=True).getPotentialEnergy()
        alchemical_reduced_energies[alchemical_snapshot, i] = energy / (u.MOLAR_GAS_CONSTANT_R * 298.15*u.kelvin)

    # Step 3: Revert back to lambda value for sampling
    alchemical_state.lambda_sterics = lambdas[16]['sterics']
    alchemical_state.lambda_electrostatics = lambdas[16]['electrostatics']
    alchemical_state.apply_to_context(sim.context)

# Dump data to file
np.save('reduced_energies_lambda_16', alchemical_reduced_energies)
print('Done!')

The simulations are afterwards analyzed using the mbar python package using the following code:

def process_k(k, u_kln_water, u_kln_octanol):
    """
    Detects equilibration and subsamples correlated data for both water and octanol
    at the given index k.

    Parameters:
    -----------
    k : int
        The lambda index to process.
    u_kln_water : np.ndarray
        The reduced potential energy matrix for water.
    u_kln_octanol : np.ndarray
        The reduced potential energy matrix for octanol.

    Returns:
    --------
    u_kln_water_sub : np.ndarray
        Subsampled reduced energies for water at index k.
    N_k_water : int
        Number of uncorrelated samples for water at index k.
    u_kln_octanol_sub : np.ndarray
        Subsampled reduced energies for octanol at index k.
    N_k_octanol : int
        Number of uncorrelated samples for octanol at index k.
    """
    # Process water
    [nequil, g, Neff_max] = timeseries.detectEquilibration(u_kln_water[k, k, :])
    indices_water = timeseries.subsampleCorrelatedData(u_kln_water[k, k, :], g=g)
    N_k_water = len(indices_water)
    u_kln_water_sub = u_kln_water[k, :, indices_water].T
    
    # Process octanol
    [nequil, g, Neff_max] = timeseries.detectEquilibration(u_kln_octanol[k, k, :])
    indices_octanol = timeseries.subsampleCorrelatedData(u_kln_octanol[k, k, :], g=g)
    N_k_octanol = len(indices_octanol)
    u_kln_octanol_sub = u_kln_octanol[k, :, indices_octanol].T
    
    return u_kln_water_sub, N_k_water, u_kln_octanol_sub, N_k_octanol

# Loading solvation free energies
u_kln_water = np.zeros(shape=(len(lambdas), len(lambdas), NumAlchemicalSnapshots))
u_kln_octanol = np.zeros(shape=(len(lambdas), len(lambdas), NumAlchemicalSnapshots))

for l in lambdas:
    l = '{0:.0f}'.format(l)
    reduced_energies = np.load('Simulations/BAR/Water/{l}/reduced_energies_lambda_{l}.npy'.format(l=l))
    u_kln_water[int(l)-1,:,:] = reduced_energies.T
    
    reduced_energies = np.load('Simulations/BAR/1-Octanol/{l}/reduced_energies_lambda_{l}.npy'.format(l=l))
    u_kln_octanol[int(l)-1,:,:] = reduced_energies.T

# Determination of number of uncorrelated samples
N_k_water = np.zeros([len(lambdas)], np.int32)
N_k_octanol = np.zeros([len(lambdas)], np.int32)

# Parallelize the processing of equilibration and subsampling
with Pool() as pool:
    results = pool.starmap(process_k, [(k, u_kln_water, u_kln_octanol) for k in range(len(lambdas)-1)])

# Collect results from parallel processing
for k, (u_kln_water_sub, N_k_water_k, u_kln_octanol_sub, N_k_octanol_k) in enumerate(results):
    u_kln_water[k, :, 0:N_k_water_k] = u_kln_water_sub
    N_k_water[k] = N_k_water_k
    
    u_kln_octanol[k, :, 0:N_k_octanol_k] = u_kln_octanol_sub
    N_k_octanol[k] = N_k_octanol_k
    
# Compute free energy differences
mbar_water = MBAR(u_kln_water, N_k_water)
mbar_octanol = MBAR(u_kln_octanol, N_k_octanol)

dG_water = mbar_water.getFreeEnergyDifferences(compute_uncertainty=True)
dG_octanol = mbar_octanol.getFreeEnergyDifferences(compute_uncertainty=True)

print('Solvation free energy of hexanol into water: {:.2f} +/- {:.2f} RT'.format(dG_water[0][-1][0], dG_water[1][-1][0]))
print('Solvation free energy of hexanol into octanol:: {:.2f} +/- {:.2f} RT'.format(dG_octanol[0][-1][0], dG_octanol[1][-1][0]))

The current output of the current code is:

Solvation free energy of hexanol into water: 4.03 +/- 0.04 RT
Solvation free energy of hexanol into octanol:: 10.83 +/- 0.17 RT

Which is obviously not in agreement with experiment or literature.
Can anybody see any obvious mistake I am making?