Skip to content

Commit

Permalink
Merging in the changes from [michellab/Sire/pull/424](michellab/Sire#424
Browse files Browse the repository at this point in the history 
)

and [michellab/Sire/pull/425](michellab/Sire#425).

Thanks to [@fjclark](https://github.com/fjclark) for submitting. Sorry for the lack
of authorship - I had to copy these files in manually across to the openbiosim repo.
  • Loading branch information
@chryswoods
chryswoods committed Jan 17, 2023
1 parent c4a4c30 commit 29439dc
Show file tree
Hide file tree
Showing 3 changed files with 69 additions and 38 deletions.
72 changes: 50 additions & 22 deletions wrapper/Tools/LJcutoff.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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, k_boltz
from Sire.Vol import PeriodicBox

# Python dependencies
#
Expand All @@ -23,7 +33,7 @@
except:
pass

import numpy as np
import numpy as _np
#try:
# import mdtraj
#except ImportError:
Expand Down Expand Up @@ -170,27 +180,42 @@ 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")]
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
# 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 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()
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
Expand Down Expand Up @@ -320,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):
Expand Down Expand Up @@ -359,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%
Expand All @@ -385,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
Expand All @@ -408,7 +434,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())
Expand All @@ -420,10 +446,12 @@ 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)
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
32 changes: 16 additions & 16 deletions wrapper/Tools/OpenMMMD.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -898,12 +898,12 @@ class 'Sire.Base._Base.Properties': The properties required to
set up the Boresch distance restraint
"""

prop = Properties()
prop = Sire.Base.Properties()

prop.setProperty("AtomNum0", VariantProperty(boresch_dict['anchor_points']['l1']))
prop.setProperty("AtomNum1", VariantProperty(boresch_dict['anchor_points']['r1']))
prop.setProperty("equil_val", VariantProperty(boresch_dict['equilibrium_values']['r0']))
prop.setProperty("force_const", VariantProperty(boresch_dict['force_constants']['kr']))
prop.setProperty("AtomNum0", Sire.Base.VariantProperty(boresch_dict['anchor_points']['l1']))
prop.setProperty("AtomNum1", Sire.Base.VariantProperty(boresch_dict['anchor_points']['r1']))
prop.setProperty("equil_val", Sire.Base.VariantProperty(boresch_dict['equilibrium_values']['r0']))
prop.setProperty("force_const", Sire.Base.VariantProperty(boresch_dict['force_constants']['kr']))

return prop

Expand All @@ -918,7 +918,7 @@ class 'Sire.Base._Base.Properties': The properties required to
set up the Boresch angle restraints
"""

prop = Properties()
prop = Sire.Base.Properties()

angle_anchor_dict = {"thetaA":["r2", "r1", "l1"], "thetaB":["r1", "l1", "l2"]}

Expand All @@ -927,15 +927,15 @@ class 'Sire.Base._Base.Properties': The properties required to
if boresch_dict["force_constants"][f"k{angle}"] != 0:
for j in range(3):
prop.setProperty(f"AtomNum{j}-{i}",
VariantProperty(boresch_dict['anchor_points'][angle_anchor_dict[angle][j]]))
Sire.Base.VariantProperty(boresch_dict['anchor_points'][angle_anchor_dict[angle][j]]))
prop.setProperty(f"equil_val-{i}",
VariantProperty(boresch_dict["equilibrium_values"][f"{angle}0"]))
Sire.Base.VariantProperty(boresch_dict["equilibrium_values"][f"{angle}0"]))
prop.setProperty(f"force_const-{i}",
VariantProperty(boresch_dict["force_constants"][f"k{angle}"]))
Sire.Base.VariantProperty(boresch_dict["force_constants"][f"k{angle}"]))

i += 1

prop.setProperty("n_boresch_angle_restraints", VariantProperty(i));
prop.setProperty("n_boresch_angle_restraints", Sire.Base.VariantProperty(i));

return prop

Expand All @@ -950,7 +950,7 @@ class 'Sire.Base._Base.Properties': The properties required to
set up the Boresch dihedral restraints
"""

prop = Properties()
prop = Sire.Base.Properties()

dihedral_anchor_dict = {"phiA":["r3", "r2", "r1", "l1"], "phiB":["r2", "r1", "l1", "l2"],
"phiC":["r1", "l1", "l2", "l3"]}
Expand All @@ -960,15 +960,15 @@ class 'Sire.Base._Base.Properties': The properties required to
if boresch_dict["force_constants"][f"k{dihedral}"] != 0:
for j in range(4):
prop.setProperty(f"AtomNum{j}-{i}",
VariantProperty(boresch_dict['anchor_points'][dihedral_anchor_dict[dihedral][j]]))
Sire.Base.VariantProperty(boresch_dict['anchor_points'][dihedral_anchor_dict[dihedral][j]]))
prop.setProperty(f"equil_val-{i}",
VariantProperty(boresch_dict["equilibrium_values"][f"{dihedral}0"]))
Sire.Base.VariantProperty(boresch_dict["equilibrium_values"][f"{dihedral}0"]))
prop.setProperty(f"force_const-{i}",
VariantProperty(boresch_dict["force_constants"][f"k{dihedral}"]))
Sire.Base.VariantProperty(boresch_dict["force_constants"][f"k{dihedral}"]))

i += 1

prop.setProperty("n_boresch_dihedral_restraints", VariantProperty(i));
prop.setProperty("n_boresch_dihedral_restraints", Sire.Base.VariantProperty(i));

return prop

Expand Down Expand Up @@ -1056,7 +1056,7 @@ def saveTurnOnRestraintsModeProperty(system):
the distance or Boresch restraint information is also stored."""
solute = getSolute(system)
solute = solute.edit().setProperty("turn_on_restraints_mode",
VariantProperty(turn_on_restraints_mode.val)).commit()
Sire.Base.VariantProperty(turn_on_restraints_mode.val)).commit()
system.update(solute)

return(system)
Expand Down
3 changes: 3 additions & 0 deletions wrapper/Tools/__init__.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

0 comments on commit 29439dc

Please sign in to comment.