Hydration Free Energy with alchemlyb

[igordiy]

Dear openmmtools team,

I am trying to calculate the hydration free energy of pyrene molecule in a water/DMSO mixture and obtain nonphysical results.

Here is a brief description of what was done:
Pyrene and DMSO have been parametrized with OFF2, water is modeled with TIP3P.
I run alchemical simulations with LangevinMiddleIntegrator, MC barostat at 303.15 K, and 1 atm pressure. Apply HBonds constraints, and restrain the COM of pyrene in the center of the box (box size is 4x4x4 nm). The nonbonded cutoff is 1.0 nm, and water is rigid. Long-range interactions are treated with PME.

For HFE I use the following protocol:

1. Minimization with unperturbed Hamiltonian
2. Equilibration with perturbed Hamiltonian for 1 ns (1 fs timestep)
3. Production run with perturbed Hamiltonian for 4 ns (2 fs timestep)
   I decouple fist the electrostatic interactions in 11 steps (from 0 to 1.0 in 0.1 steps), then with fully-decoupled electrostatics decouple the VdW interactions (also in 11 steps, from 0 to 1.0 in 0.1 steps).

Settings used for AbsoluteAlchemicalFactory:
"consistent_exceptions": false,
"disable_alchemical_dispersion_correction": true,
"alchemical_pme_treatment": "coulomb"
Other settings for this object are default

Settings for AlchemicalRegion:
"annihilate_electrostatics": false,
"annihilate_sterics": false
Other settings for this object are default

Alchemical system is created from system.
The soft-core parameters used are default (0.5, 1, 1, 6)

Since I do the postprocessing with MBAR, I evaluate the potential energy of a system with all lambda values.

Below is the snippet used to run the simulation.

`def run_simulation(self):

    if self.run_type=="production":
        
        assert all([i is not None for i in [self.readout_frequency, self.save_values_path, self.lambda_values]])

        n_iter = int(self.steps/self.readout_frequency)
        
        energies = np.zeros([n_iter, len(self.lambda_values)])

        if self.logger:
            message = f"""Starting the alchemical {self.run_type} run.
            The simulation length is {self.simulation.integrator.getStepSize()*self.steps}
            """
            self.logger.info(message)

        start_time = time.time()

        for iteration in range(n_iter):
        
            if self.alch_type == 'VDW':
                
                self.alch_state.lamda_electrostatics = 0.0
                self.alch_state.lambda_sterics = self.lambda_for_propagation

            elif self.alch_type == 'ELEC':
                
                self.alch_state.lamda_sterics = 1.0
                self.alch_state.lambda_electrostatics = self.lambda_for_propagation
   
            self.alch_state.apply_to_context(self.simulation.context)
        
            self.simulation.step(self.readout_frequency)
        
            for ilv, lambda_value in enumerate(self.lambda_values):
            
                if self.alch_type == 'VDW':
            
                    self.alch_state.lambda_sterics = lambda_value

                elif self.alch_type == 'ELEC':
    
                    self.alch_state.lambda_electrostatics = lambda_value
            
                self.alch_state.apply_to_context(self.simulation.context)
            
                energies[iteration, ilv] = self.simulation.context.getState(getEnergy=True).getPotentialEnergy()._value

        end_time = time.time()

        elapsed_time = end_time - start_time

        np.save(self.save_values_path, energies)

        if self.logger:

            message = f"""The alchemical {self.run_type} run is successfull!
            The total time needed for the run:  {elapsed_time} s.\n
            The reduced potential from alchemical calculation is saved to {self.save_values_path}."""

            self.logger.info(message)`

The results obtained for electrostatics decoupling with MBAR is more than 1000 kJ/mol between lambda 0 and 1, which unphysically high.
Full matrix after MBAR is given below:

The FF parameters yield reasonable bulk properties of a solvent, as well as reasonable binding free energy of a pyrene to some host molecule. So, I suppose the issue is not in the FF parameters. I have also checked the vacuum simulation of pyrene, the energies of lambdas are all identical then, meaning that the pyrene is correctly decoupled, and not annihilated.

I am running out of ideas what can be wrong with the pipeline, so that it gives so nonsense results. The system is neutral by the way.

Looking forward to your replies!

Best,
Igor

[igordiy]

Update:
I have tested direct-space mode instead of Coulomb mode, and it has solved the problem!
It estimated the SFE to be -11 kJ/mol, which is very close to the experiment.
Does anybody know what is wrong with the Coulomb mode?

Best,
Igor

[ijpulidos]

@FreeSoulIG Thanks for reporting this possible issue. It is hard to know what problem you are having without some minimal working example showing the problem. If you can provide one, that would be very much appreciated. On the other hand, what exactly did you change in your scripts that made it work for your case?