diff --git a/ash/interfaces/interface_OpenMM.py b/ash/interfaces/interface_OpenMM.py index 38d68c073..617702e21 100644 --- a/ash/interfaces/interface_OpenMM.py +++ b/ash/interfaces/interface_OpenMM.py @@ -1911,6 +1911,32 @@ def delete_exceptions(self, atomlist): # print_time_rel(timeA, modulename="zero_nonbondedforce") # # self.create_simulation() + # Updating LJ interactions in OpenMM object. Used to set LJ sites to zero e.g. so that they do not contribute + # Can be used to get QM-MM LJ interaction energy + def update_LJ_epsilons(self, atomlist, epsilons): + import openmm + timeA = time.time() + print("Updating LJ interaction strengths in OpenMM object.") + assert len(atomlist) == len(epsilons) + for atomindex, newepsilon in zip(atomlist, epsilons): + print("atomindex:", atomindex) + print("newepsilon:", newepsilon) + charge, sigma, oldepsilon = self.nonbonded_force.getParticleParameters(atomindex) + print("charge:", charge) + print("oldepsilon", oldepsilon) + # Different depending on type of NonbondedForce + if isinstance(self.nonbonded_force, openmm.CustomNonbondedForce): + self.nonbonded_force.setParticleParameters(atomindex, [charge, sigma, newepsilon]) + # bla1,bla2,bla3 = self.nonbonded_force.getParticleParameters(i) + # print("bla1,bla2,bla3", bla1,bla2,bla3) + elif isinstance(self.nonbonded_force, openmm.NonbondedForce): + self.nonbonded_force.setParticleParameters(atomindex, charge, sigma, newepsilon) + bla1,bla2,bla3 = self.nonbonded_force.getParticleParameters(atomindex) + print("bla1,bla2,bla3", bla1,bla2,bla3) + + printdebug("done here") + print_time_rel(timeA, modulename="update_LJ_epsilons") + # Updating charges in OpenMM object. Used to set QM charges to 0 for example # Taking list of atom-indices and list of charges (usually zero) and setting new charge # Note: Exceptions also needs to be dealt with (see delete_exceptions) diff --git a/ash/tests/test_QM_MM_openmm.py b/ash/tests/test_QM_MM_openmm.py index fe27066ee..16bcad786 100644 --- a/ash/tests/test_QM_MM_openmm.py +++ b/ash/tests/test_QM_MM_openmm.py @@ -6,7 +6,7 @@ def test_qm_mm_pyscf_nonbondedtheory_MeOH_H2O(): # H2O...MeOH fragment defined. Reading XYZ file - H2O_MeOH = Fragment(xyzfile=f"{ashpath}/tests/xyzfile/h2o_MeOH.xyz") + H2O_MeOH = Fragment(xyzfile=f"{ashpath}/tests/xyzfiles/h2o_MeOH.xyz") # Specifying the QM atoms (3-8) by atom indices (MeOH). The other atoms (0,1,2) is the H2O and MM. # IMPORTANT: atom indices begin at 0. @@ -56,7 +56,7 @@ def test_qm_mm_pyscf_nonbondedtheory_MeOH_H2O(): def test_qm_mm_pyscf_openmm_MeOH_H2O(): # H2O...MeOH fragment defined. Reading XYZ file - H2O_MeOH = Fragment(xyzfile=f"{ashpath}/tests/xyzfile/h2o_MeOH.xyz") + H2O_MeOH = Fragment(xyzfile=f"{ashpath}/tests/xyzfiles/h2o_MeOH.xyz") # Write PDB-file for OpenMM (used for topology) H2O_MeOH.write_pdbfile_openmm(filename="h2o_MeOH.pdb", skip_connectivity=True)