diff --git a/corelib/src/libs/SireMove/openmmfrenergyst.cpp b/corelib/src/libs/SireMove/openmmfrenergyst.cpp index 26e07ec75..48109223d 100644 --- a/corelib/src/libs/SireMove/openmmfrenergyst.cpp +++ b/corelib/src/libs/SireMove/openmmfrenergyst.cpp @@ -331,6 +331,52 @@ QString OpenMMFrEnergyST::toString() const return QObject::tr("OpenMMFrEnergyST()"); } +/** Takes in a molecule, system, and link properties and adds the link bonds to the system.*/ +void OpenMMFrEnergyST::addLinkBonds(OpenMM::System *system, + const SireMol::Molecule &molecule, + const SireBase::Properties &linkprop, + const QHash &atom_num_to_openmm_idx, + OpenMM::CustomBondForce *custom_force, + const bool debug) +{ + std::vector custom_bond_link_par(3); + const auto nlinks = linkprop.property(QString("nbondlinks")).asA().toInt(); + + if (debug) + qDebug() << "Adding link bonds to system"; + qDebug() << "Custom force is " << custom_force; + qDebug() << "Number of constraint links = " << nlinks; + + for (int i = 0; i < nlinks; i++) + { + const auto atomnum0 = + linkprop.property(QString("AtomNum0(%1)").arg(i)).asA().toInt(); + const auto atomnum1 = + linkprop.property(QString("AtomNum1(%1)").arg(i)).asA().toInt(); + const auto reql = linkprop.property(QString("reql(%1)").arg(i)).asA().toDouble(); + const auto kl = linkprop.property(QString("kl(%1)").arg(i)).asA().toDouble(); + const auto dl = linkprop.property(QString("dl(%1)").arg(i)).asA().toDouble(); + + const int openmmindex0 = atom_num_to_openmm_idx[atomnum0]; + const int openmmindex1 = atom_num_to_openmm_idx[atomnum1]; + + custom_bond_link_par[0] = reql * OpenMM::NmPerAngstrom; // req + custom_bond_link_par[1] = + kl * (OpenMM::KJPerKcal * OpenMM::AngstromsPerNm * OpenMM::AngstromsPerNm); // k + custom_bond_link_par[2] = dl * OpenMM::NmPerAngstrom; // dl + + if (debug) + { + qDebug() << "atomnum0 = " << atomnum0 << " openmmindex0 =" << openmmindex0; + qDebug() << "atomnum1 = " << atomnum1 << " openmmindex1 =" << openmmindex1; + qDebug() << "Req = " << reql << " kl = " << kl << " dl = " << dl; + } + + custom_force->addBond(openmmindex0, openmmindex1, custom_bond_link_par); + } + system->addForce(custom_force); +} + /** * initialises the openMM Free energy single topology calculation * Initialise must be called before anything else happens. @@ -1240,21 +1286,21 @@ void OpenMMFrEnergyST::initialise() // JM 05/23 CustomExternalForce has been deprecated as it is not recommended for implementing positional restraints // see thread https://github.com/openmm/openmm/issues/2262 - //OpenMM::CustomExternalForce *positionalRestraints_openmm = NULL; + // OpenMM::CustomExternalForce *positionalRestraints_openmm = NULL; OpenMM::HarmonicBondForce *positionalRestraints_openmm = NULL; if (Restraint_flag == true) { - //positionalRestraints_openmm = new OpenMM::CustomExternalForce("k*d2;" - // "d2 = max(0.0, d1 - d^2);" - // "d1 = (x-xref)^2 + (y-yref)^2 + (z-zref)^2"); - //positionalRestraints_openmm->addPerParticleParameter("xref"); - //positionalRestraints_openmm->addPerParticleParameter("yref"); - //positionalRestraints_openmm->addPerParticleParameter("zref"); - //positionalRestraints_openmm->addPerParticleParameter("k"); - //positionalRestraints_openmm->addPerParticleParameter("d"); - - positionalRestraints_openmm = new OpenMM::HarmonicBondForce(); + // positionalRestraints_openmm = new OpenMM::CustomExternalForce("k*d2;" + // "d2 = max(0.0, d1 - d^2);" + // "d1 = (x-xref)^2 + (y-yref)^2 + (z-zref)^2"); + // positionalRestraints_openmm->addPerParticleParameter("xref"); + // positionalRestraints_openmm->addPerParticleParameter("yref"); + // positionalRestraints_openmm->addPerParticleParameter("zref"); + // positionalRestraints_openmm->addPerParticleParameter("k"); + // positionalRestraints_openmm->addPerParticleParameter("d"); + + positionalRestraints_openmm = new OpenMM::HarmonicBondForce(); positionalRestraints_openmm->setUsesPeriodicBoundaryConditions(true); system_openmm->addForce(positionalRestraints_openmm); @@ -1283,7 +1329,7 @@ void OpenMMFrEnergyST::initialise() /****************************************BOND LINK POTENTIAL*****************************/ /* FC 12/21 CustomBondForce now (OpenMM 7.4.0) allows application of PBC checks*/ - OpenMM::CustomBondForce *custom_link_bond = new OpenMM::CustomBondForce("delta(min(0, r_eff))*lamrest*kl*r_eff^2;" + OpenMM::CustomBondForce *custom_link_bond = new OpenMM::CustomBondForce("delta(min(0, r_eff))*(lamrest^5)*kl*r_eff^2;" "r_eff=abs(r-reql)-dl"); custom_link_bond->addPerBondParameter("reql"); custom_link_bond->addPerBondParameter("kl"); @@ -1297,6 +1343,16 @@ void OpenMMFrEnergyST::initialise() else custom_link_bond->addGlobalParameter("lamrest", 1); + /*************************************PERMANENT BOND LINK POTENTIAL*****************************/ + // As above, but unaffected by lamrest + + OpenMM::CustomBondForce *custom_permanent_link_bond = new OpenMM::CustomBondForce("delta(min(0, r_eff))*kl*r_eff^2;" + "r_eff=abs(r-reql)-dl"); + custom_permanent_link_bond->addPerBondParameter("reql"); + custom_permanent_link_bond->addPerBondParameter("kl"); + custom_permanent_link_bond->addPerBondParameter("dl"); + custom_permanent_link_bond->setUsesPeriodicBoundaryConditions(true); + /****************************************BORESCH DISTANCE POTENTIAL*****************************/ OpenMM::CustomBondForce *custom_boresch_dist_rest = @@ -1841,28 +1897,28 @@ void OpenMMFrEnergyST::initialise() { //qDebug() << "atomnum " << atomnum << " openmmindex " << openmmindex << " x " << xref << " y " // << yref << " z " << zref << " k " << k << " d " << d; - qDebug() << "atomnum " << atomnum << " atopenmmindex " << atopenmmindex << " k " << k; + qDebug() << "atomnum " << atomnum << " atopenmmindex " << atopenmmindex << " k " << k; qDebug() << "anchornum " << anchornum << " anchoropenmmindex " << anchoropenmmindex << " k " << k; } - //int posrestrdim = 5; - //std::vector params(posrestrdim); - - //params[0] = xref * OpenMM::NmPerAngstrom; - //params[1] = yref * OpenMM::NmPerAngstrom; - //params[2] = zref * OpenMM::NmPerAngstrom; - //params[3] = k * (OpenMM::KJPerKcal * OpenMM::AngstromsPerNm * OpenMM::AngstromsPerNm); - //params[4] = d * OpenMM::NmPerAngstrom; - - //positionalRestraints_openmm->addParticle(openmmindex, params); - //int dummyidx = system_openmm->addParticle(0); - //std::vector custom_non_bonded_params{ 0., 0., 0., 0., 1., 1., 0., 0., 0., 1.}; - //custom_force_field->addParticle(custom_non_bonded_params); - custom_force_field->addExclusion(anchoropenmmindex, atopenmmindex); - positionalRestraints_openmm->addBond(anchoropenmmindex, atopenmmindex, - 0 * OpenMM::NmPerAngstrom, - k * (OpenMM::KJPerKcal * OpenMM::AngstromsPerNm * OpenMM::AngstromsPerNm)); + // int posrestrdim = 5; + // std::vector params(posrestrdim); + + // params[0] = xref * OpenMM::NmPerAngstrom; + // params[1] = yref * OpenMM::NmPerAngstrom; + // params[2] = zref * OpenMM::NmPerAngstrom; + // params[3] = k * (OpenMM::KJPerKcal * OpenMM::AngstromsPerNm * OpenMM::AngstromsPerNm); + // params[4] = d * OpenMM::NmPerAngstrom; + + // positionalRestraints_openmm->addParticle(openmmindex, params); + // int dummyidx = system_openmm->addParticle(0); + // std::vector custom_non_bonded_params{ 0., 0., 0., 0., 1., 1., 0., 0., 0., 1.}; + // custom_force_field->addParticle(custom_non_bonded_params); + custom_force_field->addExclusion(anchoropenmmindex, atopenmmindex); + positionalRestraints_openmm->addBond(anchoropenmmindex, atopenmmindex, + 0 * OpenMM::NmPerAngstrom, + k * (OpenMM::KJPerKcal * OpenMM::AngstromsPerNm * OpenMM::AngstromsPerNm)); } } } // end of restraint flag @@ -3049,6 +3105,7 @@ void OpenMMFrEnergyST::initialise() perturbed_energies = perturbed_energies_tmp; // IMPORTANT: PERTURBED ENERGY TORSIONS ARE ADDED ABOVE + bool UseLink_flag = true; // Distance Restraint. All the information is stored in the solute. @@ -3059,50 +3116,28 @@ void OpenMMFrEnergyST::initialise() { Molecule molecule = moleculegroup.moleculeAt(i).molecule(); - if (molecule.hasProperty("linkbonds")) + if (molecule.hasProperty("linkbonds") || molecule.hasProperty("permanent_linkbonds")) { - std::vector custom_bond_link_par(3); - - const auto linkprop = molecule.property("linkbonds").asA(); - - const auto nlinks = linkprop.property(QString("nbondlinks")).asA().toInt(); - - if (Debug) - qDebug() << "Number of constraint links = " << nlinks; - - for (int i = 0; i < nlinks; i++) + if (molecule.hasProperty("linkbonds")) { - const auto atomnum0 = - linkprop.property(QString("AtomNum0(%1)").arg(i)).asA().toInt(); - const auto atomnum1 = - linkprop.property(QString("AtomNum1(%1)").arg(i)).asA().toInt(); - const auto reql = linkprop.property(QString("reql(%1)").arg(i)).asA().toDouble(); - const auto kl = linkprop.property(QString("kl(%1)").arg(i)).asA().toDouble(); - const auto dl = linkprop.property(QString("dl(%1)").arg(i)).asA().toDouble(); - - const int openmmindex0 = AtomNumToOpenMMIndex[atomnum0]; - const int openmmindex1 = AtomNumToOpenMMIndex[atomnum1]; - - custom_bond_link_par[0] = reql * OpenMM::NmPerAngstrom; // req - custom_bond_link_par[1] = - kl * (OpenMM::KJPerKcal * OpenMM::AngstromsPerNm * OpenMM::AngstromsPerNm); // k - custom_bond_link_par[2] = dl * OpenMM::NmPerAngstrom; // dl - - if (Debug) - { - qDebug() << "atomnum0 = " << atomnum0 << " openmmindex0 =" << openmmindex0; - qDebug() << "atomnum1 = " << atomnum1 << " openmmindex1 =" << openmmindex1; - qDebug() << "Req = " << reql << " kl = " << kl << " dl = " << dl; - } - - custom_link_bond->addBond(openmmindex0, openmmindex1, custom_bond_link_par); + // Set up the link bonds + const auto linkprop = molecule.property("linkbonds").asA(); + addLinkBonds(system_openmm, molecule, linkprop, AtomNumToOpenMMIndex, + custom_link_bond, Debug); } - system_openmm->addForce(custom_link_bond); + if (molecule.hasProperty("permanent_linkbonds")) + { + // Set up the permanent link bonds + const auto linkprop = molecule.property("permanent_linkbonds").asA(); + addLinkBonds(system_openmm, molecule, linkprop, AtomNumToOpenMMIndex, + custom_permanent_link_bond, Debug); + } // We've found the molecule, exit the outer loop. break; } + } // end of loop over molecules in system } // end of bond link flag @@ -3484,8 +3519,8 @@ void OpenMMFrEnergyST::createContext(IntegratorWorkspace &workspace, SireUnits:: OpenMM::Vec3(c[j].x() * (OpenMM::NmPerAngstrom), c[j].y() * (OpenMM::NmPerAngstrom), c[j].z() * (OpenMM::NmPerAngstrom)); - //if (m[j] == 0.0) - // qDebug() << "\nWARNING - THE MASS OF PARTICLE " << system_index << " is ZERO\n"; + // if (m[j] == 0.0) + // qDebug() << "\nWARNING - THE MASS OF PARTICLE " << system_index << " is ZERO\n"; if (m[j] > SireMaths::small) { diff --git a/corelib/src/libs/SireMove/openmmfrenergyst.h b/corelib/src/libs/SireMove/openmmfrenergyst.h index 6f5ea1cdf..fded8e087 100644 --- a/corelib/src/libs/SireMove/openmmfrenergyst.h +++ b/corelib/src/libs/SireMove/openmmfrenergyst.h @@ -198,6 +198,12 @@ namespace SireMove void setDebug(bool); private: + void addLinkBonds(OpenMM::System* system, + const SireMol::Molecule& molecule, + const SireBase::Properties& linkprop, + const QHash& atom_num_to_openmm_idx, + OpenMM::CustomBondForce* custom_force, + const bool debug); void createContext(IntegratorWorkspace &workspace, SireUnits::Dimension::Time timestep); void destroyContext(); void updateBoxDimensions(OpenMM::State &state_openmm, QVector> &buffered_dimensions, diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index 051ff054f..80e55db12 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -89,6 +89,14 @@ organisation on `GitHub `__. Restraints can be named, meaning that you can scale different restraints at different stages and by different values across the λ-coordinate. +* Added support for one or more "permanent" distance restraints. These are + distance restraints that are always applied, and are never scaled by λ. + This allows the release of all other distance restraints to a single + harmonic or flat-bottomed restraint. When the ligand is fully decoupled, + the free energy of release of the single remaining restraint can be + computed without simulation. See + for more details. + * Please add the changelog entry for your PR here. We will add the link to your PR during the code review :-) diff --git a/tests/somd/test_boresch_corr.py b/tests/somd/test_boresch_corr.py index 8392cc64e..2aaf7e604 100644 --- a/tests/somd/test_boresch_corr.py +++ b/tests/somd/test_boresch_corr.py @@ -8,6 +8,7 @@ try: import scipy + have_scipy = True except ImportError: have_scipy = False @@ -70,9 +71,7 @@ def test_boresch_analytical_correction(tmpdir): print(output) # Check that the correction is correct - data = output.split( - "Analytical correction for releasing Boresch restraints =" - ) + data = output.split("Analytical correction for releasing Boresch restraints =") assert len(data) > 1 @@ -102,9 +101,7 @@ def test_boresch_numerical_correction(tmpdir): print(output) - data = output.split( - "Numerical correction for releasing Boresch restraints =" - ) + data = output.split("Numerical correction for releasing Boresch restraints =") assert len(data) > 1 diff --git a/wrapper/Tools/OpenMMMD.py b/wrapper/Tools/OpenMMMD.py index 8786e7b0d..6114c80f3 100644 --- a/wrapper/Tools/OpenMMMD.py +++ b/wrapper/Tools/OpenMMMD.py @@ -26,9 +26,7 @@ openmm_dir = os.environ["OPENMM_PLUGIN_DIR"] except KeyError: # Set to the default location of the bundled OpenMM package. - openmm_dir = os.environ["OPENMM_PLUGIN_DIR"] = ( - Sire.Base.getLibDir() + "/plugins" - ) + openmm_dir = os.environ["OPENMM_PLUGIN_DIR"] = Sire.Base.getLibDir() + "/plugins" finally: if not os.path.isdir(openmm_dir): warnings.warn("OPENMM_PLUGIN_DIR is not accessible") @@ -80,9 +78,7 @@ "temperature", 25.0 * Sire.Units.celsius, """Simulation temperature""" ) -pressure = Parameter( - "pressure", 1.0 * Sire.Units.atm, """Simulation pressure""" -) +pressure = Parameter("pressure", 1.0 * Sire.Units.atm, """Simulation pressure""") topfile = Parameter( "topfile", @@ -107,8 +103,7 @@ restart_file = Parameter( "restart file", "sim_restart.s3", - "Filename of the restart file to use to save progress during the " - "simulation.", + "Filename of the restart file to use to save progress during the " "simulation.", ) dcd_root = Parameter( @@ -120,8 +115,7 @@ nmoves = Parameter( "nmoves", 1000, - "Number of Molecular Dynamics moves to be performed during the " - "simulation.", + "Number of Molecular Dynamics moves to be performed during the " "simulation.", ) debug_seed = Parameter( @@ -164,8 +158,7 @@ minimal_coordinate_saving = Parameter( "minimal coordinate saving", False, - "Reduce the number of coordiantes writing for states" - "with lambda in ]0,1[", + "Reduce the number of coordiantes writing for states" "with lambda in ]0,1[", ) time_to_skip = Parameter( @@ -356,9 +349,26 @@ "Dictionary of pair of atoms whose distance is restrained, and restraint " "parameters. Syntax is {(atom0,atom1):(reql, kl, Dl)} where atom0, atom1 " "are atomic indices. reql the equilibrium distance. Kl the force constant " - "of the restraint. D the flat bottom radius. WARNING: PBC distance checks " - "not implemented, avoid restraining pair of atoms that may diffuse out of " - "the box.", + "of the restraint. D the flat bottom radius.", +) + +use_permanent_distance_restraints = Parameter( + "use permanent distance restraints", + False, + """ + Whether or not to use permanent (not affected by turn-on-restraints mode) + distance restraints distances between pairs of atoms. + """, +) + +permanent_distance_restraints_dict = Parameter( + "permanent distance restraints dictionary", + {}, + "Dictionary of pair of atoms whose distance is restrained, and restraint " + "parameters. Syntax is {(atom0,atom1):(reql, kl, Dl)} where atom0, atom1 " + "are atomic indices. reql the equilibrium distance. Kl the force constant " + "of the restraint. D the flat bottom radius. Permanent restraints are not " + "affected by turn-on-restraints mode and are always at full strength.", ) turn_on_restraints_mode = Parameter( @@ -416,8 +426,7 @@ morphfile = Parameter( "morphfile", "MORPH.pert", - "Name of the morph file containing the perturbation to apply to the " - "system.", + "Name of the morph file containing the perturbation to apply to the " "system.", ) lambda_val = Parameter( @@ -528,9 +537,7 @@ def setupDCD(system): def writeSystemData(system, moves, Trajectory, block, softcore_lambda=False): if softcore_lambda: if block == ncycles.val or block == 1: - Trajectory.writeModel( - system[MGName("all")], system.property("space") - ) + Trajectory.writeModel(system[MGName("all")], system.property("space")) else: if block % ncycles_per_snap.val == 0: if buffered_coords_freq.val > 0: @@ -539,13 +546,9 @@ def writeSystemData(system, moves, Trajectory, block, softcore_lambda=False): for prop in sysprops: if prop.startswith("buffered_space"): dimensions[str(prop)] = system.property(prop) - Trajectory.writeBufferedModels( - system[MGName("all")], dimensions - ) + Trajectory.writeBufferedModels(system[MGName("all")], dimensions) else: - Trajectory.writeModel( - system[MGName("all")], system.property("space") - ) + Trajectory.writeModel(system[MGName("all")], system.property("space")) # Write a PDB coordinate file each cycle. pdb = Sire.IO.PDB2(system, {"use_atom_numbers": Sire.Base.wrap(True)}) @@ -596,9 +599,7 @@ def centerSolute(system, space): box_center = space.dimensions() / 2 # TriclincBox. except: - box_center = 0.5 * ( - space.vector0() + space.vector1() + space.vector2() - ) + box_center = 0.5 * (space.vector0() + space.vector1() + space.vector2()) else: box_center = Sire.Maths.Vector(0.0, 0.0, 0.0) @@ -683,14 +684,10 @@ def setupForcefields(system, space): inter_ions_nonbondedff.add(ions) - inter_ions_molecules_nonbondedff = Sire.MM.InterGroupCLJFF( - "ions:molecules" - ) + inter_ions_molecules_nonbondedff = Sire.MM.InterGroupCLJFF("ions:molecules") if cutoff_type.val != "nocutoff": inter_ions_molecules_nonbondedff.setUseReactionField(True) - inter_ions_molecules_nonbondedff.setReactionFieldDielectric( - rf_dielectric.val - ) + inter_ions_molecules_nonbondedff.setReactionFieldDielectric(rf_dielectric.val) inter_ions_molecules_nonbondedff.add(ions, MGIdx(0)) inter_ions_molecules_nonbondedff.add(molecules, MGIdx(1)) @@ -758,9 +755,7 @@ def setupForcefields(system, space): system.setProperty( "switchingFunction", Sire.MM.CHARMMSwitchingFunction(cutoff_dist.val) ) - system.setProperty( - "combiningRules", Sire.Base.VariantProperty(combining_rules.val) - ) + system.setProperty("combiningRules", Sire.Base.VariantProperty(combining_rules.val)) total_nrg = ( internonbondedff.components().total() @@ -873,9 +868,7 @@ def atomNumVectorListToProperty(list): i = 0 for value in list: - prop.setProperty( - "AtomNum(%d)" % i, Sire.Base.VariantProperty(value[0].value()) - ) + prop.setProperty("AtomNum(%d)" % i, Sire.Base.VariantProperty(value[0].value())) prop.setProperty("x(%d)" % i, Sire.Base.VariantProperty(value[1].x())) prop.setProperty("y(%d)" % i, Sire.Base.VariantProperty(value[1].y())) prop.setProperty("z(%d)" % i, Sire.Base.VariantProperty(value[1].z())) @@ -888,17 +881,26 @@ def atomNumVectorListToProperty(list): def linkbondVectorListToProperty(list): + """ + Converts a list of linkbond vectors to a Sire property. + + Parameters + ---------- + list : list + List of linkbond vectors. + + Returns + ------- + prop : Sire.Base.Properties + Sire property containing the linkbond vectors. + """ prop = Sire.Base.Properties() i = 0 for value in list: - prop.setProperty( - "AtomNum0(%d)" % i, Sire.Base.VariantProperty(value[0]) - ) - prop.setProperty( - "AtomNum1(%d)" % i, Sire.Base.VariantProperty(value[1]) - ) + prop.setProperty("AtomNum0(%d)" % i, Sire.Base.VariantProperty(value[0])) + prop.setProperty("AtomNum1(%d)" % i, Sire.Base.VariantProperty(value[1])) prop.setProperty("reql(%d)" % i, Sire.Base.VariantProperty(value[2])) prop.setProperty("kl(%d)" % i, Sire.Base.VariantProperty(value[3])) prop.setProperty("dl(%d)" % i, Sire.Base.VariantProperty(value[4])) @@ -908,6 +910,7 @@ def linkbondVectorListToProperty(list): return prop + def PositionalRestraintsToProperty(posrest_list): """Generates properties to store positional restraint parameters. Args: @@ -930,6 +933,7 @@ class 'Sire.Base._Base.Properties': The properties required to return prop + def boreschDistRestraintToProperty(boresch_dict): """Generates properties to store information needed to set up the single Boresch distance restraint. @@ -988,9 +992,7 @@ class 'Sire.Base._Base.Properties': The properties required to prop.setProperty( f"AtomNum{j}-{i}", Sire.Base.VariantProperty( - boresch_dict["anchor_points"][ - angle_anchor_dict[angle][j] - ] + boresch_dict["anchor_points"][angle_anchor_dict[angle][j]] ), ) prop.setProperty( @@ -1001,16 +1003,12 @@ class 'Sire.Base._Base.Properties': The properties required to ) prop.setProperty( f"force_const-{i}", - Sire.Base.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", Sire.Base.VariantProperty(i) - ) + prop.setProperty("n_boresch_angle_restraints", Sire.Base.VariantProperty(i)) return prop @@ -1041,9 +1039,7 @@ class 'Sire.Base._Base.Properties': The properties required to prop.setProperty( f"AtomNum{j}-{i}", Sire.Base.VariantProperty( - boresch_dict["anchor_points"][ - dihedral_anchor_dict[dihedral][j] - ] + boresch_dict["anchor_points"][dihedral_anchor_dict[dihedral][j]] ), ) prop.setProperty( @@ -1061,9 +1057,7 @@ class 'Sire.Base._Base.Properties': The properties required to i += 1 - prop.setProperty( - "n_boresch_dihedral_restraints", Sire.Base.VariantProperty(i) - ) + prop.setProperty("n_boresch_dihedral_restraints", Sire.Base.VariantProperty(i)) return prop @@ -1099,6 +1093,7 @@ def propertyToAtomNumVectorList(prop): return list + def setupPositionalRestraints(system): # Tag atoms in the system that are listed in the provided restrained atoms file # this is used to activate positional restraints during initialisation @@ -1106,7 +1101,9 @@ def setupPositionalRestraints(system): # parse restrained atoms file if restrained_atoms is None: - print ("""Error: if positional restraints are activated you must supply a keyword argument restrained atoms = /path/to/file.""") + print( + """Error: if positional restraints are activated you must supply a keyword argument restrained atoms = /path/to/file.""" + ) sys.exit(-1) restrained_dict = restrained_atoms.val keys = restrained_dict.keys() @@ -1127,7 +1124,9 @@ def setupPositionalRestraints(system): atnumber = at.number() atnumberval = atnumber.value() if atnumberval in keys: - restrainedAtoms.append((atnumberval, restrained_dict[atnumberval] , k_restraint.val)) + restrainedAtoms.append( + (atnumberval, restrained_dict[atnumberval], k_restraint.val) + ) if len(restrainedAtoms) > 0: mol = ( @@ -1142,6 +1141,7 @@ def setupPositionalRestraints(system): return system + def setupRestraints(system): molecules = system[MGName("all")].molecules() @@ -1206,14 +1206,38 @@ def saveTurnOnRestraintsModeProperty(system): return system -def setupDistanceRestraints(system, restraints=None): +def setupDistanceRestraints(system, restraints=None, permanent=False): + """ + Sets up distance restraints between pairs of atoms in the system. + + Parameters + ---------- + system : Sire System + The system to which the restraints will be added. + restraints : dict, optional, default=None + A dictionary of the form {(atom0,atom1):(reql, kl, Dl)} where atom0, atom1 + are atomic indices, reql is the equilibrium distance, Kl is the force constant + of the restraint and Dl is the flat bottom radius. If none, this will be read + from the cfg file. + permanent : bool, optional, default=False + If True, the restraints are 'permanent' and will not be scaled by lambda + in turn-on-restraints mode (they will always be at full strength). + + Returns + ------- + None + """ prop_list = [] molecules = system[MGName("all")].molecules() + # Get the provided restraints, or read them from the cfg file + # if none are provided if restraints is None: - # dic_items = list(distance_restraints_dict.val.items()) - dic_items = list(dict(distance_restraints_dict.val).items()) + if not permanent: + dic_items = list(dict(distance_restraints_dict.val).items()) + else: # Permanent restraints + dic_items = list(dict(permanent_distance_restraints_dict.val).items()) else: dic_items = list(restraints.items()) @@ -1257,10 +1281,12 @@ def setupDistanceRestraints(system, restraints=None): print(unique_prop_list) # The solute will store all the information related to the receptor-ligand restraints solute = getSolute(system) + property_name = "linkbonds" if not permanent else "permanent_linkbonds" solute = ( solute.edit() .setProperty( - "linkbonds", linkbondVectorListToProperty(unique_prop_list) + property_name, + linkbondVectorListToProperty(unique_prop_list), ) .commit() ) @@ -1360,9 +1386,7 @@ def setupBoreschRestraints(system): print(anchor + "=" + str(at)) if anchors_not_present: - print( - "Error! The following anchor points do not not exist in the system:" - ) + print("Error! The following anchor points do not not exist in the system:") for anchor in anchors_not_present: print(f"{anchor}: index {anchors_dict[anchor]-1}") sys.exit(-1) @@ -1412,13 +1436,10 @@ def freezeResidues(system): atnumber = at.number() if at.residue().name().value() in frozen_residues.val: print( - "Freezing %s %s %s " - % (at, atnumber, at.residue().name().value()) + "Freezing %s %s %s " % (at, atnumber, at.residue().name().value()) ) mol = ( - at.edit() - .setProperty("mass", 0.0 * Sire.Units.g_per_mol) - .molecule() + at.edit().setProperty("mass", 0.0 * Sire.Units.g_per_mol).molecule() ) system.update(mol) @@ -1522,8 +1543,7 @@ def repartitionMasses(system, hmassfactor=4.0): print( "WARNING! The mass repartitioning algorithm is not conserving " "atomic masses for molecule %s (total delta is %s). " - "Report bug to a Sire developer." - % (molnum, total_delta.value()), + "Report bug to a Sire developer." % (molnum, total_delta.value()), ) sys.exit(-1) @@ -1544,12 +1564,7 @@ def repartitionMasses(system, hmassfactor=4.0): ) sys.exit(-1) - mol = ( - mol.edit() - .atom(atidx) - .setProperty("mass", newmass)[0] - .molecule() - ) + mol = mol.edit().atom(atidx).setProperty("mass", newmass)[0].molecule() system.update(mol) # import pdb; pdb.set_trace() @@ -1567,26 +1582,54 @@ def getDummies(molecule): to_dummies = None from_non_dummies = [] to_non_dummies = [] + # For simulations in which restraints are modified while ligand is non-interacting + always_dummies = None for x in range(0, natoms): atom = atoms[x] - if atom.property("initial_ambertype") == "du": - if from_dummies is None: - from_dummies = molecule.selectAll(atom.index()) + # First, check for special case where we are going from a dummy to a + # dummy. + if ( + atom.property("initial_ambertype") == "du" + and atom.property("final_ambertype") == "du" + ): + if always_dummies is None: + always_dummies = molecule.selectAll(atom.index()) else: - from_dummies += molecule.selectAll(atom.index()) + always_dummies += molecule.selectAll(atom.index()) + # Otherwise, this is a standard stage in a free energy calculation + # where we are going from a dummy to a non-dummy or vice versa. else: - from_non_dummies.append(atom.index()) - if atom.property("final_ambertype") == "du": - if to_dummies is None: - to_dummies = molecule.selectAll(atom.index()) + if atom.property("initial_ambertype") == "du": + if from_dummies is None: + from_dummies = molecule.selectAll(atom.index()) + else: + from_dummies += molecule.selectAll(atom.index()) else: - to_dummies += molecule.selectAll(atom.index()) - else: - to_non_dummies.append(atom.index()) + 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()) + + # Check that if we have any "always_dummies" then we don't have any other types of dummies + # because we expect the entire ligand to be non-interacting (and the logic to handle this + # is not implemented). + if always_dummies is not None: + if from_dummies is not None or to_dummies is not None: + raise RuntimeError( + "Cannot have any dummy atoms in a ligand that is always non-interacting" + ) + if not (len(from_non_dummies) == 0 and len(to_non_dummies) == 0): + raise RuntimeError( + "Cannot have any non-dummy atoms in a ligand that is always non-interacting" + ) - return to_dummies, from_dummies, to_non_dummies, from_non_dummies + return to_dummies, from_dummies, to_non_dummies, from_non_dummies, always_dummies def createSystemFreeEnergy(molecules): @@ -1618,9 +1661,7 @@ def createSystemFreeEnergy(molecules): # FIXME: assumption that perturbed molecule is at a fixed location for molecule in moleculeList: - if molecule.residue(ResIdx(0)).number() == ResNum( - perturbed_resnum.val - ): + if molecule.residue(ResIdx(0)).number() == ResNum(perturbed_resnum.val): solute = molecule moleculeList.remove(molecule) break @@ -1657,12 +1698,13 @@ def createSystemFreeEnergy(molecules): .commit() ) - # We put atoms in three groups depending on what happens in the + # We put atoms in four groups depending on what happens in the # perturbation # non dummy to non dummy --> the hard group, use a normal intermolecular FF # non dummy to dummy --> the todummy group, uses SoftFF with alpha = Lambda # dummy to non dummy --> the fromdummy group, uses SoftFF with # alpha = 1 - Lambda + # dummy to dummy --> the alwaysdummy group, uses SoftFF with alpha = 1 # We start assuming all atoms are hard atoms. Then we call getDummies to # find which atoms # start/end as dummies and update the hard, todummy and fromdummy groups @@ -1672,10 +1714,12 @@ def createSystemFreeEnergy(molecules): solute_grp_ref_hard = MoleculeGroup("solute_ref_hard") solute_grp_ref_todummy = MoleculeGroup("solute_ref_todummy") solute_grp_ref_fromdummy = MoleculeGroup("solute_ref_fromdummy") + solute_grp_ref_alwaysdummy = MoleculeGroup("solute_ref_alwaysdummy") solute_ref_hard = solute.selectAllAtoms() solute_ref_todummy = solute.selectAllAtoms() solute_ref_fromdummy = solute.selectAllAtoms() + solute_ref_alwaysdummy = 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 @@ -1688,9 +1732,13 @@ def createSystemFreeEnergy(molecules): # 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 - ) + ( + to_dummies, + from_dummies, + to_non_dummies, + from_non_dummies, + always_dummies, + ) = getDummies(solute) if to_dummies is not None: ndummies = to_dummies.count() @@ -1698,14 +1746,10 @@ 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_hard = solute_ref_hard.subtract(solute.select(dummy_index)) for non_dummy in to_non_dummies: - solute_ref_todummy = solute_ref_todummy.subtract( - solute.select(non_dummy) - ) + solute_ref_todummy = solute_ref_todummy.subtract(solute.select(non_dummy)) if from_dummies is not None: ndummies = from_dummies.count() @@ -1713,18 +1757,34 @@ def createSystemFreeEnergy(molecules): for x in range(0, ndummies): dummy_index = dummies[x].index() - solute_ref_hard = solute_ref_hard.subtract( + solute_ref_hard = solute_ref_hard.subtract(solute.select(dummy_index)) + + for non_dummy in from_non_dummies: + solute_ref_fromdummy = solute_ref_fromdummy.subtract(solute.select(non_dummy)) + + if always_dummies is not None: + # Remove all the atoms from the hard group, the todummy group and the + # fromdummy group + ndummies = always_dummies.count() + dummies = always_dummies.atoms() + + 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.subtract(solute.select(dummy_index)) + solute_ref_fromdummy = solute_ref_fromdummy.subtract( solute.select(dummy_index) ) - for non_dummy in from_non_dummies: - solute_ref_fromdummy = solute_ref_fromdummy.subtract( - solute.select(non_dummy) + else: # We have no alwaysdummy atoms - remove all atoms from the always dummy group + solute_ref_alwaysdummy = solute_ref_alwaysdummy.subtract( + solute.selectAllAtoms() ) solute_grp_ref_hard.add(solute_ref_hard) solute_grp_ref_todummy.add(solute_ref_todummy) solute_grp_ref_fromdummy.add(solute_ref_fromdummy) + solute_grp_ref_alwaysdummy.add(solute_ref_alwaysdummy) solutes = MoleculeGroup("solutes") solutes.add(solute) @@ -1748,6 +1808,7 @@ def createSystemFreeEnergy(molecules): all.add(solute_grp_ref_hard) all.add(solute_grp_ref_todummy) all.add(solute_grp_ref_fromdummy) + all.add(solute_grp_ref_alwaysdummy) # Add these groups to the System system = Sire.System.System() @@ -1758,6 +1819,7 @@ def createSystemFreeEnergy(molecules): system.add(solute_grp_ref_hard) system.add(solute_grp_ref_todummy) system.add(solute_grp_ref_fromdummy) + system.add(solute_grp_ref_alwaysdummy) system.add(molecules) @@ -1791,6 +1853,7 @@ def setupForceFieldsFreeEnergy(system, space): solute_hard = system[MGName("solute_ref_hard")] solute_todummy = system[MGName("solute_ref_todummy")] solute_fromdummy = system[MGName("solute_ref_fromdummy")] + solute_alwaysdummy = system[MGName("solute_ref_alwaysdummy")] solvent = system[MGName("solvent")] @@ -1831,9 +1894,7 @@ def setupForceFieldsFreeEnergy(system, space): solute_todummy_intraclj.setReactionFieldDielectric(rf_dielectric.val) solute_todummy_intraclj.add(solute_todummy) - solute_fromdummy_intraclj = Sire.MM.IntraSoftCLJFF( - "solute_fromdummy_intraclj" - ) + solute_fromdummy_intraclj = Sire.MM.IntraSoftCLJFF("solute_fromdummy_intraclj") solute_fromdummy_intraclj.setShiftDelta(shift_delta.val) solute_fromdummy_intraclj.setCoulombPower(coulomb_power.val) if cutoff_type.val != "nocutoff": @@ -1841,6 +1902,14 @@ def setupForceFieldsFreeEnergy(system, space): solute_fromdummy_intraclj.setReactionFieldDielectric(rf_dielectric.val) solute_fromdummy_intraclj.add(solute_fromdummy) + solute_alwaysdummy_intraclj = Sire.MM.IntraSoftCLJFF("solute_alwaysdummy_intraclj") + solute_alwaysdummy_intraclj.setShiftDelta(shift_delta.val) + solute_alwaysdummy_intraclj.setCoulombPower(coulomb_power.val) + if cutoff_type.val != "nocutoff": + solute_alwaysdummy_intraclj.setUseReactionField(True) + solute_alwaysdummy_intraclj.setReactionFieldDielectric(rf_dielectric.val) + solute_alwaysdummy_intraclj.add(solute_alwaysdummy) + solute_hard_todummy_intraclj = Sire.MM.IntraGroupSoftCLJFF( "solute_hard:todummy_intraclj" ) @@ -1848,9 +1917,7 @@ def setupForceFieldsFreeEnergy(system, space): solute_hard_todummy_intraclj.setCoulombPower(coulomb_power.val) if cutoff_type.val != "nocutoff": solute_hard_todummy_intraclj.setUseReactionField(True) - solute_hard_todummy_intraclj.setReactionFieldDielectric( - rf_dielectric.val - ) + solute_hard_todummy_intraclj.setReactionFieldDielectric(rf_dielectric.val) solute_hard_todummy_intraclj.add(solute_hard, MGIdx(0)) solute_hard_todummy_intraclj.add(solute_todummy, MGIdx(1)) @@ -1861,9 +1928,7 @@ def setupForceFieldsFreeEnergy(system, space): solute_hard_fromdummy_intraclj.setCoulombPower(coulomb_power.val) if cutoff_type.val != "nocutoff": solute_hard_fromdummy_intraclj.setUseReactionField(True) - solute_hard_fromdummy_intraclj.setReactionFieldDielectric( - rf_dielectric.val - ) + solute_hard_fromdummy_intraclj.setReactionFieldDielectric(rf_dielectric.val) solute_hard_fromdummy_intraclj.add(solute_hard, MGIdx(0)) solute_hard_fromdummy_intraclj.add(solute_fromdummy, MGIdx(1)) @@ -1874,9 +1939,7 @@ def setupForceFieldsFreeEnergy(system, space): solute_todummy_fromdummy_intraclj.setCoulombPower(coulomb_power.val) if cutoff_type.val != "nocutoff": solute_todummy_fromdummy_intraclj.setUseReactionField(True) - solute_todummy_fromdummy_intraclj.setReactionFieldDielectric( - rf_dielectric.val - ) + solute_todummy_fromdummy_intraclj.setReactionFieldDielectric(rf_dielectric.val) solute_todummy_fromdummy_intraclj.add(solute_todummy, MGIdx(0)) solute_todummy_fromdummy_intraclj.add(solute_fromdummy, MGIdx(1)) @@ -1888,23 +1951,17 @@ def setupForceFieldsFreeEnergy(system, space): solute_hard_solventff.add(solute_hard, MGIdx(0)) solute_hard_solventff.add(solvent, MGIdx(1)) - solute_todummy_solventff = Sire.MM.InterGroupSoftCLJFF( - "solute_todummy:solvent" - ) + solute_todummy_solventff = Sire.MM.InterGroupSoftCLJFF("solute_todummy:solvent") if cutoff_type.val != "nocutoff": solute_todummy_solventff.setUseReactionField(True) solute_todummy_solventff.setReactionFieldDielectric(rf_dielectric.val) solute_todummy_solventff.add(solute_todummy, MGIdx(0)) solute_todummy_solventff.add(solvent, MGIdx(1)) - solute_fromdummy_solventff = Sire.MM.InterGroupSoftCLJFF( - "solute_fromdummy:solvent" - ) + solute_fromdummy_solventff = Sire.MM.InterGroupSoftCLJFF("solute_fromdummy:solvent") if cutoff_type.val != "nocutoff": solute_fromdummy_solventff.setUseReactionField(True) - solute_fromdummy_solventff.setReactionFieldDielectric( - rf_dielectric.val - ) + solute_fromdummy_solventff.setReactionFieldDielectric(rf_dielectric.val) solute_fromdummy_solventff.add(solute_fromdummy, MGIdx(0)) solute_fromdummy_solventff.add(solvent, MGIdx(1)) @@ -1914,6 +1971,7 @@ def setupForceFieldsFreeEnergy(system, space): solute_hard_intraclj, solute_todummy_intraclj, solute_fromdummy_intraclj, + solute_alwaysdummy_intraclj, solute_hard_todummy_intraclj, solute_hard_fromdummy_intraclj, solute_todummy_fromdummy_intraclj, @@ -1938,15 +1996,9 @@ def setupForceFieldsFreeEnergy(system, space): else: system.setProperty("switchingFunction", Sire.MM.NoCutoff()) - system.setProperty( - "combiningRules", Sire.Base.VariantProperty(combining_rules.val) - ) - system.setProperty( - "coulombPower", Sire.Base.VariantProperty(coulomb_power.val) - ) - system.setProperty( - "shiftDelta", Sire.Base.VariantProperty(shift_delta.val) - ) + system.setProperty("combiningRules", Sire.Base.VariantProperty(combining_rules.val)) + system.setProperty("coulombPower", Sire.Base.VariantProperty(coulomb_power.val)) + system.setProperty("shiftDelta", Sire.Base.VariantProperty(shift_delta.val)) # TOTAL total_nrg = ( @@ -1954,6 +2006,7 @@ def setupForceFieldsFreeEnergy(system, space): + solute_hard_intraclj.components().total() + solute_todummy_intraclj.components().total(0) + solute_fromdummy_intraclj.components().total(0) + + solute_alwaysdummy_intraclj.components().total(0) + solute_hard_todummy_intraclj.components().total(0) + solute_hard_fromdummy_intraclj.components().total(0) + solute_todummy_fromdummy_intraclj.components().total(0) @@ -1987,6 +2040,11 @@ def setupForceFieldsFreeEnergy(system, space): "alpha0", Sire.FF.FFName("solute_fromdummy_intraclj"), 1 - lam ) ) + system.add( + Sire.System.PropertyConstraint( + "alpha0", Sire.FF.FFName("solute_alwaysdummy_intraclj"), 1 + ) + ) system.add( Sire.System.PropertyConstraint( "alpha0", Sire.FF.FFName("solute_hard:todummy_intraclj"), lam @@ -2336,9 +2394,7 @@ def computeOpenMMEnergy(prmtop_filename, inpcrd_filename, cutoff): state = context.getState(getEnergy=True) - return state.getPotentialEnergy().value_in_unit( - units.kilocalorie / units.mole - ) + return state.getPotentialEnergy().value_in_unit(units.kilocalorie / units.mole) ### This is how a TIP3P specific template for Cl- looks like @@ -2485,15 +2541,9 @@ def run(): print("\n### Running Molecular Dynamics simulation on %s ###" % host) if verbose.val: - print( - "###================= Simulation Parameters=====================" - "###" - ) + print("###================= Simulation Parameters=====================" "###") Parameter.printAll() - print( - "###===========================================================" - "###\n" - ) + print("###===========================================================" "###\n") timer = Sire.Qt.QTime() timer.start() @@ -2544,12 +2594,13 @@ def run(): % distance_restraints_dict ) stream = open("restraints.cfg", "w") - stream.write( - "distance restraints dictionary = %s\n" % restraints - ) + stream.write("distance restraints dictionary = %s\n" % restraints) stream.close() system = setupDistanceRestraints(system, restraints=restraints) + if use_permanent_distance_restraints.val: + system = setupDistanceRestraints(system, permanent=True) + if use_boresch_restraints.val: print("Setting up Boresch restraints...") system = setupBoreschRestraints(system) @@ -2567,10 +2618,7 @@ def run(): system = setupForcefields(system, space) if debug_seed.val != 0: - print( - "Setting up the simulation with debugging seed %s" - % debug_seed.val - ) + print("Setting up the simulation with debugging seed %s" % debug_seed.val) moves = setupMoves(system, debug_seed.val, gpu.val) @@ -2584,9 +2632,7 @@ def run(): move0.setIntegrator(integrator) moves = Sire.Move.WeightedMoves() moves.add(move0) - print( - "Index GPU = %s " % moves.moves()[0].integrator().getDeviceIndex() - ) + print("Index GPU = %s " % moves.moves()[0].integrator().getDeviceIndex()) print( "Loaded a restart file on which we have performed %d moves." % moves.nMoves() @@ -2612,17 +2658,13 @@ def run(): mdmoves = moves.moves()[0] integrator = mdmoves.integrator() - print( - "###===========================================================###\n" - ) + print("###===========================================================###\n") # energy = computeOpenMMEnergy(topfile.val, crdfile.val, cutoff_dist.val) # print(f'OpenMM Energy (PME): {energy}\n') if minimise.val: - print( - "###=======================Minimisation========================###" - ) + print("###=======================Minimisation========================###") print("Running minimization.") if verbose.val: print("Energy before the minimization: " + str(system.energy())) @@ -2641,22 +2683,17 @@ def run(): print("Energy minimization done.") integrator.setConstraintType(constraint.val) print( - "###===========================================================" - "###\n", + "###===========================================================" "###\n", flush=True, ) if equilibrate.val: - print( - "###======================Equilibration========================###" - ) + print("###======================Equilibration========================###") print("Running equilibration.") # Here we anneal lambda (To be determined) if verbose.val: print("Equilibration timestep " + str(equil_timestep.val)) - print( - "Number of equilibration steps: " + str(equil_iterations.val) - ) + print("Number of equilibration steps: " + str(equil_iterations.val)) system = integrator.equilibrateSystem( system, equil_timestep.val, equil_iterations.val ) @@ -2665,18 +2702,14 @@ def run(): print("Energy after the equilibration: " + str(system.energy())) print("Equilibration done.\n") print( - "###===========================================================" - "###\n", + "###===========================================================" "###\n", flush=True, ) simtime = nmoves.val * ncycles.val * timestep.val print("###=======================somd run============================###") print("Starting somd run...") - print( - "%s moves %s cycles, %s simulation time" - % (nmoves.val, ncycles.val, simtime) - ) + print("%s moves %s cycles, %s simulation time" % (nmoves.val, ncycles.val, simtime)) s1 = timer.elapsed() / 1000.0 for i in range(cycle_start, cycle_end): @@ -2707,15 +2740,9 @@ def runFreeNrg(): ) if verbose.val: - print( - "###================= Simulation Parameters=====================" - "###" - ) + print("###================= Simulation Parameters=====================" "###") Parameter.printAll() - print( - "###===========================================================" - "###\n" - ) + print("###===========================================================" "###\n") timer = Sire.Qt.QTime() timer.start() @@ -2771,11 +2798,14 @@ def runFreeNrg(): % distance_restraints_dict ) stream = open("restraints.cfg", "w") - stream.write( - "distance restraints dictionary = %s\n" % restraints - ) + stream.write("distance restraints dictionary = %s\n" % restraints) stream.close() - system = setupDistanceRestraints(system, restraints=restraints) + system = setupDistanceRestraints( + system, restraints=restraints, permanent=False + ) + + if use_permanent_distance_restraints.val: + system = setupDistanceRestraints(system, permanent=True) if use_boresch_restraints.val: print("Setting up Boresch restraints...") @@ -2794,18 +2824,13 @@ def runFreeNrg(): system = setupForceFieldsFreeEnergy(system, space) if debug_seed.val != 0: - print( - "Setting up the simulation with debugging seed %s" - % debug_seed.val - ) + print("Setting up the simulation with debugging seed %s" % debug_seed.val) if charge_diff.val != 0: print("The difference in charge is", charge_diff.val) system = selectWatersForPerturbation(system, charge_diff.val) - moves = setupMovesFreeEnergy( - system, debug_seed.val, gpu.val, lambda_val.val - ) + moves = setupMovesFreeEnergy(system, debug_seed.val, gpu.val, lambda_val.val) print("Saving restart") Sire.Stream.save([system, moves], restart_file.val) @@ -2842,9 +2867,7 @@ def runFreeNrg(): "UTF-8", ) ) - outfile.write( - bytes("#Generating lambda is\t\t " + lam_str + "\n", "UTF-8") - ) + outfile.write(bytes("#Generating lambda is\t\t " + lam_str + "\n", "UTF-8")) outfile.write( bytes( "#Alchemical array is\t\t " + str(lambda_array.val) + "\n", @@ -2890,9 +2913,7 @@ def runFreeNrg(): moves.add(move0) cycle_start = int(moves.nMoves() / nmoves.val) cycle_end = cycle_start + ncycles.val - print( - "Index GPU = %s " % moves.moves()[0].integrator().getDeviceIndex() - ) + print("Index GPU = %s " % moves.moves()[0].integrator().getDeviceIndex()) print( "Loaded a restart file on which we have performed %d moves." % moves.nMoves() @@ -2903,8 +2924,7 @@ def runFreeNrg(): if cycle_start > maxcycles.val: print( "Maxinum number of MD cycles reached (%s). If you wish to extend the " - "simulation increase the value of the parameter maxcycle." - % maxcycles.val + "simulation increase the value of the parameter maxcycle." % maxcycles.val ) sys.exit(-1) @@ -2922,9 +2942,7 @@ def runFreeNrg(): mdmoves = moves.moves()[0] integrator = mdmoves.integrator() - print( - "###===========================================================###\n" - ) + print("###===========================================================###\n") print(f"Initial energy: {integrator.getPotentialEnergy(system)}") @@ -2933,9 +2951,7 @@ def runFreeNrg(): # f'({cutoff_type}): {energy:.2f} kcal mol-1\n') if minimise.val: - print( - "###=======================Minimisation========================###" - ) + print("###=======================Minimisation========================###") print("Running minimization.") print(f"Tolerance for minimization: {str(minimise_tol.val)}") print( @@ -2950,26 +2966,18 @@ def runFreeNrg(): system.mustNowRecalculateFromScratch() print( - "Energy after the minimization: " - f"{integrator.getPotentialEnergy(system)}" + "Energy after the minimization: " f"{integrator.getPotentialEnergy(system)}" ) print("Energy minimization done.") - print( - "###===========================================================" - "###\n" - ) + print("###===========================================================" "###\n") if equilibrate.val: - print( - "###======================Equilibration========================###" - ) + print("###======================Equilibration========================###") print("Running lambda equilibration to lambda=%s." % lambda_val.val) # Here we anneal lambda (To be determined) if verbose.val: print("Equilibration timestep " + str(equil_timestep.val)) - print( - "Number of equilibration steps: " + str(equil_iterations.val) - ) + print("Number of equilibration steps: " + str(equil_iterations.val)) system = integrator.annealSystemToLambda( system, equil_timestep.val, equil_iterations.val ) @@ -2979,17 +2987,11 @@ def runFreeNrg(): f"Energy after the annealing: {integrator.getPotentialEnergy(system)}" ) print("Lambda annealing done.\n") - print( - "###===========================================================" - "###\n" - ) + print("###===========================================================" "###\n") print("###====================somd-freenrg run=======================###") print("Starting somd-freenrg run...") - print( - "%s moves %s cycles, %s simulation time" - % (nmoves.val, ncycles.val, simtime) - ) + print("%s moves %s cycles, %s simulation time" % (nmoves.val, ncycles.val, simtime)) softcore_lambda = False if minimal_coordinate_saving.val: @@ -3045,9 +3047,7 @@ def runFreeNrg(): outgradients.close() outfile.close() print("Simulation took %d s " % (s2 - s1)) - print( - "###===========================================================###\n" - ) + print("###===========================================================###\n") if os.path.exists("gradients.s3"): siregrads = Sire.Stream.load("gradients.s3")