From a5cef0749c2e854cfb1097aff1fd44ba57b3edaf Mon Sep 17 00:00:00 2001 From: finlayclark Date: Thu, 12 Jan 2023 15:35:03 +0000 Subject: [PATCH 1/5] Fix imports for LJCutoff [ci skip] --- wrapper/Tools/LJcutoff.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/wrapper/Tools/LJcutoff.py b/wrapper/Tools/LJcutoff.py index ba0e0096a..5fdbd9b40 100644 --- a/wrapper/Tools/LJcutoff.py +++ b/wrapper/Tools/LJcutoff.py @@ -8,8 +8,18 @@ import os,sys, random import math +from Sire.Base import VariantProperty +from Sire.CAS import Symbol, Max +from Sire.FF import FFName +from Sire.IO import Amber +from Sire.Maths import FreeEnergyAverage +from Sire.MM import IntraCLJFF, InterCLJFF, IntraSoftCLJFF, IntraGroupSoftCLJFF, InterGroupSoftCLJFF, InterGroupCLJFF, CHARMMSwitchingFunction +from Sire.Mol import AtomIdx, MGIdx, MoleculeGroup, MGName, AtomCoords +from Sire.System import System, PropertyConstraint, PerturbationConstraint from Sire.Tools.OpenMMMD import * from Sire.Tools import Parameter, resolveParameters +from Sire.Units import gram, centimeter, angstrom, g_per_mol, kcal_per_mol, mod_electron +from Sire.Vol import PeriodicBox # Python dependencies # From f9b5cd77583c6af4761b8074b1c5b966c156b7bb Mon Sep 17 00:00:00 2001 From: finlayclark Date: Thu, 12 Jan 2023 15:47:40 +0000 Subject: [PATCH 2/5] Fix typo [ci skip] --- wrapper/Tools/LJcutoff.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/wrapper/Tools/LJcutoff.py b/wrapper/Tools/LJcutoff.py index 5fdbd9b40..59b2bfb64 100644 --- a/wrapper/Tools/LJcutoff.py +++ b/wrapper/Tools/LJcutoff.py @@ -418,7 +418,7 @@ def runLambda(): #print (system_shortc.energy()) #print (system_longc.energy()) system_shortc = updateSystemfromTraj(system_shortc, frames_xyz, cell_lengths, cell_angles) - system_longcc = updateSystemfromTraj(system_longc, frames_xyz, cell_lengths, cell_angles) + system_longc = updateSystemfromTraj(system_longc, frames_xyz, cell_lengths, cell_angles) #print (system_shortc.energy()) #print (system_longc.energy()) delta_nrg = (system_longc.energy()+E_lrc_full - system_shortc.energy()) From 366f0e9a06a38893534038062c0a630718679ef5 Mon Sep 17 00:00:00 2001 From: finlayclark Date: Thu, 12 Jan 2023 16:08:41 +0000 Subject: [PATCH 3/5] Fix search for solvent molecule so that protein can not be selected, and that all molecules are searched (first and last molecules could both be non-solvent) [ci skip] --- wrapper/Tools/LJcutoff.py | 34 ++++++++++++++++++++++++++-------- 1 file changed, 26 insertions(+), 8 deletions(-) diff --git a/wrapper/Tools/LJcutoff.py b/wrapper/Tools/LJcutoff.py index 59b2bfb64..e897468ac 100644 --- a/wrapper/Tools/LJcutoff.py +++ b/wrapper/Tools/LJcutoff.py @@ -180,6 +180,18 @@ def setupLJFF(system, space, cutoff=10* angstrom): def addAnalyticalLRC(system, cutoff, bulk_density): + molecules = system[MGName("molecules")] + molnums = molecules.molNums() + #updated = [] + for molnum in molnums: + mol = molecules.molecule(molnum)[0].molecule() + editmol = mol.edit() + for x in range(0,mol.nAtoms()): + editatom = editmol.atom(AtomIdx(x)) + editatom.setProperty("charge", 0 * mod_electron) + editmol = editatom.molecule() + mol = editmol.commit() + system.update(mol) # 1) Collect all solvent particles # 2) Average all solvent sigma/epsilon parameters solvent = system[MGName("solvent")] @@ -192,15 +204,21 @@ def addAnalyticalLRC(system, cutoff, bulk_density): avg_sigma = 0.0 * angstrom avg_epsilon = 0.0 * kcal_per_mol LJsites = 0 - # Is this always the right molecule? - mol = solvent_mols.first()[0].molecule() - if (mol.nAtoms() == 1): - print ("This does not seem to be a solvent molecule. Picking up another one...") - mol = solvent_mols.last().molecule() - if (mol.nAtoms() == 1): - print ("This also does not seem to be a solvent molecule ! Abort!") + + # Loop over all solvent molecules to find something which + # is not an ion or a protein + # TODO: Make this more robust - there could be other non-ligand + # molecules with n_atoms > 1 < 30 + n_solvent_mols = len(solvent_molnums) + for i, molnum in enumerate(solvent_molnums): + mol = solvent_mols.molecule(molnum)[0].molecule() + if (mol.nAtoms() != 1) and (mol.nAtoms() < 30): + print (f"Solvent molecule identified with {mol.nAtoms()} atoms") + break + if i == n_solvent_mols - 1: + print ("Unable to identify solvent molecule ! Abort!") sys.exit(-1) - #molecule(molnum).molecule() + atoms = mol.atoms() natoms = mol.nAtoms() solv_mol_mass = 0 * g_per_mol From 859a60b47d5709a0271384a5ff3e673cf39f4851 Mon Sep 17 00:00:00 2001 From: finlayclark Date: Fri, 13 Jan 2023 11:25:50 +0000 Subject: [PATCH 4/5] Subtract mean from free energies to avoid overflow during exponential averaging [ci skip] --- wrapper/Tools/LJcutoff.py | 28 +++++++++++++--------------- 1 file changed, 13 insertions(+), 15 deletions(-) diff --git a/wrapper/Tools/LJcutoff.py b/wrapper/Tools/LJcutoff.py index e897468ac..fc834f7e9 100644 --- a/wrapper/Tools/LJcutoff.py +++ b/wrapper/Tools/LJcutoff.py @@ -18,7 +18,7 @@ from Sire.System import System, PropertyConstraint, PerturbationConstraint from Sire.Tools.OpenMMMD import * from Sire.Tools import Parameter, resolveParameters -from Sire.Units import gram, centimeter, angstrom, g_per_mol, kcal_per_mol, mod_electron +from Sire.Units import gram, centimeter, angstrom, g_per_mol, kcal_per_mol, mod_electron, k_boltz from Sire.Vol import PeriodicBox # Python dependencies @@ -33,7 +33,7 @@ except: pass -import numpy as np +import numpy as _np #try: # import mdtraj #except ImportError: @@ -197,10 +197,7 @@ def addAnalyticalLRC(system, cutoff, bulk_density): solvent = system[MGName("solvent")] solvent_mols = solvent.molecules() solvent_molnums = solvent_mols.molNums() - # - # What if solvent contains more than one type of molecule? - # - #for molnum in solvent_molnums: + avg_sigma = 0.0 * angstrom avg_epsilon = 0.0 * kcal_per_mol LJsites = 0 @@ -208,7 +205,7 @@ def addAnalyticalLRC(system, cutoff, bulk_density): # Loop over all solvent molecules to find something which # is not an ion or a protein # TODO: Make this more robust - there could be other non-ligand - # molecules with n_atoms > 1 < 30 + # molecules with n_atoms > 1 < 30 or more than one type of solvent n_solvent_mols = len(solvent_molnums) for i, molnum in enumerate(solvent_molnums): mol = solvent_mols.molecule(molnum)[0].molecule() @@ -348,10 +345,16 @@ def updateSystemfromTraj(system, frame_xyz, cell_lengths, cell_angles): return system def getFreeEnergy(delta_nrgs): + # Subtract mean value from all energies + # to avoid overflow during exponentiation + delta_nrgs = _np.array([nrg.value() for nrg in delta_nrgs]) + mean_nrg = delta_nrgs.mean() + delta_nrgs -= mean_nrg free_nrg = FreeEnergyAverage(temperature.val) for nrg in delta_nrgs: - free_nrg.accumulate(nrg.value()) - deltaG = free_nrg.average() * kcal_per_mol + free_nrg.accumulate(nrg) + # Add back the mean + deltaG = (free_nrg.average() + mean_nrg)* kcal_per_mol return deltaG def resample(values): @@ -387,15 +390,11 @@ def runLambda(): # !!! NEED TO DISABLE CHANGE IN COULOMBIC CUTOFF !! #system = zeroCharges(system) - #import pdb; pdb.set_trace() - # THIS IS THE ONE WITH SHORT CUTOFF system_shortc = System() system_shortc.copy(system) - #import pdb; pdb.set_trace() system_shortc = setupLJFF(system_shortc, space, \ cutoff=cutoff_dist.val) - #import pdb; pdb.set_trace() # Determine longest cutoff that can be used. Take lowest space dimension, # and decrease by 5% @@ -413,7 +412,6 @@ def runLambda(): cutoff=long_cutoff) # NOW ADD ANALYTICAL CORRECTION TERM TO longc E_lrc_full = addAnalyticalLRC(system_longc, long_cutoff, bulk_rho.val) - #import pdb; pdb.set_trace() # Now loop over snapshots in dcd and accumulate energies start_frame = 1 end_frame = 1000000000 @@ -448,7 +446,7 @@ def runLambda(): deltaG = getFreeEnergy(delta_nrgs) #print (deltaG) nbootstrap = 100 - deltaG_bootstrap = np.zeros(nbootstrap) + deltaG_bootstrap = _np.zeros(nbootstrap) for x in range(0,nbootstrap): resampled_nrgs = resample(delta_nrgs) dG = getFreeEnergy(resampled_nrgs) From 0cdb133e07dbbfa5e734240800b048241e37d02f Mon Sep 17 00:00:00 2001 From: finlayclark Date: Fri, 13 Jan 2023 13:27:20 +0000 Subject: [PATCH 5/5] Add return value to runLambda and modify resolveParameters to allow testing of LJ correction [ci skip] --- wrapper/Tools/LJcutoff.py | 2 ++ wrapper/Tools/__init__.py | 3 +++ 2 files changed, 5 insertions(+) diff --git a/wrapper/Tools/LJcutoff.py b/wrapper/Tools/LJcutoff.py index fc834f7e9..a133b171a 100644 --- a/wrapper/Tools/LJcutoff.py +++ b/wrapper/Tools/LJcutoff.py @@ -453,3 +453,5 @@ def runLambda(): deltaG_bootstrap[x] = dG.value() dev = deltaG_bootstrap.std() print ("DG_LJ = %8.5f +/- %8.5f kcal/mol (1 sigma) " % (deltaG.value(), dev)) + + return deltaG, dev \ No newline at end of file diff --git a/wrapper/Tools/__init__.py b/wrapper/Tools/__init__.py index a0b368d1e..9ae9b262c 100644 --- a/wrapper/Tools/__init__.py +++ b/wrapper/Tools/__init__.py @@ -174,10 +174,13 @@ def inner(params = {}): try: retval = func() except: + retval = None sys.exc_info()[0] Parameter.pop() raise Parameter.pop() + return retval + return inner