Merge pull request #6 from rinikerlab/CPU_MacOS

Publication updates and CPU compatibility

Author: katzberger

Analysis/Analyse_Speed_of_simulations.ipynb

@@ -301,6 +301,87 @@ def Data(self, Data):
         self._Data = Data
 
 
+class OpenFF_forcefield_OBC(OpenFF_forcefield):
+
+    def __init__(
+        self,
+        pdb_id,
+        solvent_model="TIP3P",
+        SA=None,
+        cache=None,
+        solvent_dielectric=78.5,
+        rdkit_mol=None,
+        partial_charges=None,
+        forcefield="openff-2.0.0",
+        constraints=HBonds,
+    ):
+        super().__init__(
+            pdb_id,
+            solvent_model,
+            cache=cache,
+            rdkit_mol=rdkit_mol,
+            partial_charges=partial_charges,
+            forcefield=forcefield,
+            constraints=constraints,
+        )
+        self._SA = SA
+        self._solvent_dielectric = solvent_dielectric
+
+    def create_system(
+        self, topology, nonbondedMethod=PME, nonbondedCutoff=1 * nanometer
+    ):
+
+        system = self._openmm_forcefield.createSystem(
+            topology=topology, nonbondedMethod=NoCutoff, constraints=HBonds
+        )
+        charges = np.array(
+            [
+                system.getForces()[0].getParticleParameters(i)[0]._value
+                for i in range(topology._numAtoms)
+            ]
+        )
+
+        force = GBSAOBC2Force(
+            cutoff=None,
+            SA=self._SA,
+            soluteDielectric=1,
+            solventDielectric=self._solvent_dielectric,
+        )
+        obc_parameters = np.empty((topology.getNumAtoms(), 3))
+        obc_parameters[:, 0] = charges  # Charges
+        obc_parameters[:, 1:] = force.getStandardParameters(
+            topology
+        )  # GBNeck2 parameters
+
+        self._Data = obc_parameters
+        # Add Particles and finalize force
+        force.addParticles(obc_parameters)
+        force.finalize()
+
+        # Create System and add force
+        system.addForce(force)
+
+        return system
+
+    def __str__(self):
+        if self._SA is None:
+            return "openff200_OBC_%.1f" % self._solvent_dielectric
+        else:
+            return "openff200_SAOBC_%.1f" % self._solvent_dielectric
+
+    @property
+    def water_model(self):
+        return "implicit"
+
+    @property
+    def Data(self):
+        return self._Data
+
+    @Data.setter
+    def Data(self, Data):
+        self._Data = Data
+
+
 class OpenFF_forcefield_SAGBNeck2(OpenFF_forcefield_GBNeck2):
 
     def __init__(self, pdb_id, solvent_model="TIP3P", SA="ACE", cache=None):

Analysis/Intra_molecular_HBond.ipynb

@@ -787,7 +787,7 @@ def _compute_energy_tensors(self, positions, solvent_model, solvent_dielectric):
             # Encode solvent
 
         # Get Graph
-        if positions.device != self._device:
+        if positions.device != torch.device("cuda"):
             positions = positions.float().to(self._device)
 
         # Build Graph

ForceField/Forcefield.py

@@ -2,14 +2,14 @@
 
 ## Publications
 
-[3] ChemRxiv. 2024; DOI: [https://doi.org/10.26434/chemrxiv-2024-1hb0b](https://doi.org/10.26434/chemrxiv-2024-1hb0b)
+[3] _J. Am. Chem. Soc._, **2025**, DOI: [https://doi.org/10.1021/jacs.4c17622](https://doi.org/10.1021/jacs.4c17622)
 
-[2] Chem. Sci., 2024, DOI: [https://doi.org/10.1039/D4SC02432J](https://doi.org/10.1039/D4SC02432J)
+[2] _Chem. Sci._, **2024**, DOI: [https://doi.org/10.1039/D4SC02432J](https://doi.org/10.1039/D4SC02432J)
 
-[1] J. Chem. Phys. 158, 204101 (2023), DOI: [https://doi.org/10.1063/5.0147027](https://doi.org/10.1063/5.0147027)
+[1] _J. Chem. Phys._, **2023**, DOI: [https://doi.org/10.1063/5.0147027](https://doi.org/10.1063/5.0147027)
 
 ## Abstract
-[3] Understanding and manipulating the conformational behavior of a molecule in different solvent environments is of great interest in the fields of drug discovery and organic synthesis. Molecular dynamics (MD) simulations with solvent molecules explicitly present are the gold standard to compute such conformational ensembles (within the accuracy of the underlying force field), complementing experimental findings and supporting their interpretation. However, conventional methods often face challenges related to computational cost (explicit solvent) or accuracy (implicit solvent). Here, we showcase how our graph neural network (GNN)-based implicit solvent (GNNIS) approach can be used to rapidly compute small molecule conformational ensembles in 39 common organic solvents with high accuracy compared to explicit-solvent simulations. We validate this approach using nuclear magnetic resonance (NMR) measurements, thus identifying the conformers contributing most to the experimental observable. The method allows the time required to accurately predict conformational ensembles to be reduced from days to minutes while achieving results within one kBT of the experimental values.
+[3] Understanding and manipulating the conformational behavior of a molecule in different solvent environments is of great interest in the fields of drug discovery and organic synthesis. Molecular dynamics (MD) simulations with solvent molecules explicitly present are the gold standard to compute such conformational ensembles (within the accuracy of the underlying force field), complementing experimental findings and supporting their interpretation. However, conventional methods often face challenges related to computational cost (explicit solvent) or accuracy (implicit solvent). Here, we showcase how our graph neural network (GNN)-based implicit solvent (GNNIS) approach can be used to rapidly compute small molecule conformational ensembles in 39 common organic solvents reproducing explicit-solvent simulations with high accuracy. We validate this approach using nuclear magnetic resonance (NMR) measurements, thus identifying the conformers contributing most to the experimental observable. The method allows the time required to accurately predict conformational ensembles to be reduced from days to minutes while achieving results within one kBT of the experimental values.
 
 [2] The dynamical behavior of small molecules in their environment can be studied with classical molecular dynamics (MD) simulations to gain deeper insight on an atomic level and thus complement and rationalize the interpretation of experimental findings. Such approaches are of great value in various areas of research, e.g., in the development of new therapeutics. The accurate description of solvation effects in such simulations is thereby key and has in consequence been an active field of research since the introduction of MD. So far, the most accurate approaches involve computationally expensive explicit solvent simulations, while widely applied models using an implicit solvent description suffer from reduced accuracy. Recently, machine learning (ML) approaches that provide a probabilistic representation of solvation effects have been proposed as potential alternatives. However, the associated computational costs and minimal or lack of transferability render them unusable in practice. Here, we report the first example of a transferable ML-based implicit solvent model trained on a diverse set of 3 000 000 molecular structures that can be applied to organic small molecules for simulations in water. Extensive testing against reference calculations demonstrated that the model delivers on par accuracy with explicit solvent simulations while providing an up to 18-fold increase in sampling rate.
 
@@ -26,7 +26,7 @@
 
 ## Installation
 
-First clone this repository to your work station and install the environmentusing conda or mamba:
+First clone this repository to your work station and install the environment using conda or mamba:
 
 ```bash
 mamba env create -f environment.yml
@@ -57,7 +57,7 @@ from rdkit.Chem import AllChem
 
 mol = Chem.MolFromSmiles('COCCO')
 mol = Chem.AddHs(mol)
-AllChem.EmbedMultipleConfs(mol, numConfs=128)
+AllChem.EmbedMultipleConfs(mol, numConfs=128, useExpTorsionAnglePrefs = False)
 
 minimized_mol, energies = minimize_mol(mol,"DMSO")
 entropies, free_energies = calculate_entropy(minimized_mol,"DMSO")
@@ -68,6 +68,7 @@ In addition an example workflow is provided in the [ExampleConformationEnsemble.
 ## Reproducibility
 
 This section is intended to provide a step-by-step guide to reproduce the results of the paper.
+For the exact versions of the packages used in this study please refer to the [reproduction_environment.yml](reproduction_environment.yml) file.
 
 ### Data Set Generation
 

MachineLearning/GNN_Models.py

@@ -450,6 +450,7 @@ def create_vac_sim(
         save_name=save_name,
         openff_forcefield=forcefield,
         constraints=constraints,
+        rdkit_mol=mol,
     )
     vac_sim.forcefield = OpenFF_forcefield_vacuum(
         pdb_id,
@@ -1164,7 +1165,7 @@ def get_gnn_sim(
             constraints=constraints,
         )
     else:
-        model_dict = torch.load(model_path)["model"]
+        model_dict = torch.load(model_path, map_location="cpu")["model"]
         gnn_sim = create_gnn_sim(
             smiles,
             cache=cache,

README.md

@@ -38,6 +38,7 @@
     OpenFF_forcefield_vacuum,
     OpenFF_forcefield_SAGBNeck2,
     OpenFF_TIP5P_forcefield,
+    OpenFF_forcefield_OBC,
 )
 import yaml
 from openmm import LangevinMiddleIntegrator, MonteCarloBarostat
@@ -49,6 +50,7 @@
 forcefield_dict = DefaultDict(lambda: OpenFF_forcefield)
 forcefield_dict["SAGBNeck2"] = OpenFF_forcefield_SAGBNeck2
 forcefield_dict["GBNeck2"] = OpenFF_forcefield_GBNeck2
+forcefield_dict["OBC"] = OpenFF_forcefield_OBC
 forcefield_dict["vac"] = OpenFF_forcefield_vacuum
 forcefield_dict["TIP5P"] = OpenFF_TIP5P_forcefield
 
@@ -59,7 +61,7 @@
 
 def solvent_model_dict(x):
 
-    if x in ["GBNeck2", "vac", "SAGBNeck2"]:
+    if x in ["GBNeck2", "vac", "SAGBNeck2", "OBC"]:
         return "v"
     if x in solvent_dict.keys():
         return solvent_dict[x]["SMILES"]
@@ -97,6 +99,18 @@ def solvent_model_dict(x):
         solvent_dielectric=args.solvent_dielectric,
         cache="run_caches/" + save_name + ".cache",
     )
+elif args.solvent == "OBC":
+    sim.forcefield = forcefield_dict[args.solvent](
+        pdb_id,
+        solvent_dielectric=args.solvent_dielectric,
+        cache="run_caches/" + save_name + ".cache",
+    )
+elif args.solvent == "SAGBNeck2":
+    sim.forcefield = forcefield_dict[args.solvent](
+        pdb_id,
+        solvent_dielectric=args.solvent_dielectric,
+        cache="run_caches/" + save_name + ".cache",
+    )
 else:
     sim.forcefield = forcefield_dict[args.solvent](
         pdb_id,