From 1771000486d609c8691bbc9450c978f5bd45994a Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Wed, 16 Nov 2022 12:07:15 +0000 Subject: [PATCH] Apply SOMD fixes. --- wrapper/Tools/OpenMMMD.py | 122 +++++++++++++++++++++++--------------- 1 file changed, 74 insertions(+), 48 deletions(-) diff --git a/wrapper/Tools/OpenMMMD.py b/wrapper/Tools/OpenMMMD.py index 74fa1ba3d..2cd14be96 100644 --- a/wrapper/Tools/OpenMMMD.py +++ b/wrapper/Tools/OpenMMMD.py @@ -307,10 +307,10 @@ def getSolute(system): # matching the perturbed_resnum.val. # Create the query string. - query = f"mol with resnum {perturbed_resnum.val}" + query = f"resnum {perturbed_resnum.val}" # Perform the search. - search = system.search(query) + search = system.search(query).molecules() # Make sure there is only one result. if len(search) != 1: @@ -858,6 +858,8 @@ def getDummies(molecule): from_dummies = None to_dummies = None + from_non_dummies = [] + to_non_dummies = [] for x in range(0, natoms): atom = atoms[x] @@ -866,13 +868,17 @@ def getDummies(molecule): from_dummies = molecule.selectAll(atom.index()) else: from_dummies += molecule.selectAll(atom.index()) - elif atom.property("final_ambertype") == "du": + else: + from_non_dummies.append(atom.index()) + if atom.property("final_ambertype") == "du": if to_dummies is None: to_dummies = molecule.selectAll(atom.index()) else: to_dummies += molecule.selectAll(atom.index()) + else: + to_non_dummies.append(atom.index()) - return to_dummies, from_dummies + return to_dummies, from_dummies, to_non_dummies, from_non_dummies def createSystemFreeEnergy(molecules): @@ -946,10 +952,21 @@ def createSystemFreeEnergy(molecules): solute_grp_ref_fromdummy = MoleculeGroup("solute_ref_fromdummy") solute_ref_hard = solute.selectAllAtoms() - solute_ref_todummy = solute_ref_hard.invert() - solute_ref_fromdummy = solute_ref_hard.invert() - - to_dummies, from_dummies = getDummies(solute) + solute_ref_todummy = solute.selectAllAtoms() + solute_ref_fromdummy = solute.selectAllAtoms() + + # N.B.: Currently Sire 2023 doesn't behave consistently with Sire < 2023 for + # selections with zero items. Previously it was possible to create an empty + # selector on a molecule by using .selectAllAtoms().invert(), then add items + # to that selector. For Sire 2023, an empty selector is Null, with no + # associated molecule. To work around this, we instead create an "all atom" + # selector, then remove unwanted atoms. An alternative approach is to use + # the .selection() operator on the solute, to create an AtomSelection. It is + # then possible to call .selectNone(), followed by repeated calls to .select() + # in order to add only the desired atoms. This can then be converted to a + # Selector_Atom_ by passing the solute and selection to the constructor. + + to_dummies, from_dummies, to_non_dummies, from_non_dummies = getDummies(solute) if to_dummies is not None: ndummies = to_dummies.count() @@ -958,7 +975,9 @@ def createSystemFreeEnergy(molecules): for x in range(0, ndummies): dummy_index = dummies[x].index() solute_ref_hard = solute_ref_hard.subtract(solute.select(dummy_index)) - solute_ref_todummy = solute_ref_todummy.add(solute.select(dummy_index)) + + for non_dummy in to_non_dummies: + solute_ref_todummy = solute_ref_todummy.subtract(solute.select(non_dummy)) if from_dummies is not None: ndummies = from_dummies.count() @@ -967,7 +986,9 @@ def createSystemFreeEnergy(molecules): for x in range(0, ndummies): dummy_index = dummies[x].index() solute_ref_hard = solute_ref_hard.subtract(solute.select(dummy_index)) - solute_ref_fromdummy = solute_ref_fromdummy.add(solute.select(dummy_index)) + + for non_dummy in from_non_dummies: + solute_ref_fromdummy = solute_ref_fromdummy.subtract(solute.select(non_dummy)) solute_grp_ref_hard.add(solute_ref_hard) solute_grp_ref_todummy.add(solute_ref_todummy) @@ -1357,7 +1378,7 @@ def generateDistanceRestraintsDict(system): # Find atom nearest to COG molecules = system.molecules() molnums = molecules.molNums() - solute = molecules.at(MolNum(1))[0].molecule() + solute = getSolute(system) nearestcog_atom = getAtomNearCOG( solute ) icoord = nearestcog_atom.property("coordinates") # Step 2) Find nearest 'CA' heavy atom in other solutes (skip water & ions) @@ -1377,39 +1398,43 @@ def generateDistanceRestraintsDict(system): if d < dmin: dmin = d closest = ca + + restraints = None + # Step 3) Compute position of 'mirror' CA. Find nearest CA atom to that point - jcoord = closest.property("coordinates") - mirror_coord = icoord-(jcoord-icoord) - dmin = 9999999.0 - mirror_closest = None - for molnum in molnums: - molecule = molecules.molecule(molnum)[0].molecule() - if molecule == solute: - continue - if molecule.residues()[0].name() == ResName("WAT"): - continue - #print (molecule) - ca_atoms = molecule.selectAll(AtomName("CA")) - for ca in ca_atoms: - jcoord = ca.property("coordinates") - d = Vector().distance(mirror_coord,jcoord) - if d < dmin: - dmin = d - mirror_closest = ca - #print (mirror_closest) - # Step 4) Setup restraint parameters - kl = 10.00 # kcal/mol/Angstrom^2 - Dl = 2.00 # Angstrom - i0 = nearestcog_atom.index().value() - i1 = closest.index().value() - i2 = mirror_closest.index().value() - r01 = Vector().distance(nearestcog_atom.property("coordinates"),closest.property("coordinates")) - r02 = Vector().distance(nearestcog_atom.property("coordinates"),mirror_closest.property("coordinates")) - restraints = { (i0, i1): (r01, kl, Dl), (i0,i2): (r02, kl, Dl) } - #print restraints - #distance_restraints_dict.val = restraints - #distance_restraints_dict - #import pdb; pdb.set_trace() + if closest: + jcoord = closest.property("coordinates") + mirror_coord = icoord-(jcoord-icoord) + dmin = 9999999.0 + mirror_closest = None + for molnum in molnums: + molecule = molecules.molecule(molnum)[0].molecule() + if molecule == solute: + continue + if molecule.residues()[0].name() == ResName("WAT"): + continue + #print (molecule) + ca_atoms = molecule.selectAll(AtomName("CA")) + for ca in ca_atoms: + jcoord = ca.property("coordinates") + d = Vector().distance(mirror_coord,jcoord) + if d < dmin: + dmin = d + mirror_closest = ca + #print (mirror_closest) + # Step 4) Setup restraint parameters + kl = 10.00 # kcal/mol/Angstrom^2 + Dl = 2.00 # Angstrom + i0 = nearestcog_atom.index().value() + i1 = closest.index().value() + i2 = mirror_closest.index().value() + r01 = Vector().distance(nearestcog_atom.property("coordinates"),closest.property("coordinates")) + r02 = Vector().distance(nearestcog_atom.property("coordinates"),mirror_closest.property("coordinates")) + restraints = { (i0, i1): (r01, kl, Dl), (i0,i2): (r02, kl, Dl) } + #print restraints + #distance_restraints_dict.val = restraints + #distance_restraints_dict + #import pdb; pdb.set_trace() return restraints @@ -1459,11 +1484,12 @@ def run(): if len(distance_restraints_dict.val) == 0: print ("Distance restraints have been activated, but none have been specified. Will autogenerate.") restraints = generateDistanceRestraintsDict(system) - # Save restraints - print ("Autogenerated distance restraints values: %s " % distance_restraints_dict) - stream = open("restraints.cfg",'w') - stream.write("distance restraints dictionary = %s\n" % restraints) - stream.close() + if restraints: + # Save restraints + print ("Autogenerated distance restraints values: %s " % distance_restraints_dict) + stream = open("restraints.cfg",'w') + stream.write("distance restraints dictionary = %s\n" % restraints) + stream.close() system = setupDistanceRestraints(system, restraints=restraints) if hydrogen_mass_repartitioning_factor.val > 1.0: