Request for Assistance: Detailed MM Energy Components in QM/MM Single-Point Energy Calculations

[0xBytefuzz]

When calculating single-point energy using the QM/MM method, I found that some frames have abnormal values in the MM part, while the QM part is reasonable upon inspection. Visual inspection did not reveal any obviously unreasonable areas in the structure. I would like to know how to print out the detailed components of the MM energy. Any tips you can provide would be very helpful. Thank you very much.

[RagnarB83]

Because of the way OpenMM does things, getting the MM forcefield energy decomposed into different contributions is unfortunately a bit of a hazzle. You can see some discussion on this here: https://github.com/openmm/openmm/issues/387

In ASH there is an option to OpenMMTheory : do_energy_decomposition=True that does work for CHARMM forcefields to get a decomposition automatically but only if the system was prepared using CHARMM files (not CHARMM via xml-files).

For other setups you would have to do something like deleting all forces except one to get the contribution of each force. See script below. This is not perfect but can give some insight. If the issue has to do with Coulomb interaction vs. Lennard-Jones interaction then this option will not help because this is all grouped together within the NonbondedForce. For that one could maybe do something like setting all LJ interactions to zero by setting all epsilons to zero in the forcefield file or the pointcharges.

from ash import *
forcefielddir="./"
psffile=forcefielddir+"step3_pbcsetup.psf"
topfile=forcefielddir+"top_all36_prot.rtf"
prmfile=forcefielddir+"par_all36_prot.prm"
xyzfile=forcefielddir+"coordinates.xyz"

frag = Fragment(xyzfile=xyzfile)

omm = OpenMMTheory(CHARMMfiles=True, psffile=psffile, charmmtopfile=topfile, autoconstraints=None, rigidwater=False,
    charmmprmfile=prmfile, periodic=True, periodic_cell_dimensions=[80, 80, 80, 90, 90, 90],
    dispersion_correction=False, periodic_nonbonded_cutoff=12, switching_function_distance=10,
    PMEparameters=[1.0/0.34, 90, 90, 90], do_energy_decomposition=True)

result = Singlepoint(theory=omm, fragment=frag)
energy_dict={}
energy_dict["Complete"] = result.energy

#Looping over all forces, deleting all forces except 1 particular force
for i in range(0, len(omm.system.getForces())):
    omm = OpenMMTheory(CHARMMfiles=True, psffile=psffile, charmmtopfile=topfile, autoconstraints=None, rigidwater=False,
        charmmprmfile=prmfile, periodic=True, periodic_cell_dimensions=[80, 80, 80, 90, 90, 90],
        dispersion_correction=False, periodic_nonbonded_cutoff=12, switching_function_distance=10,
        PMEparameters=[1.0/0.34, 90, 90, 90])

#Get forces in object
    forces = omm.system.getForces()
    forcenames = [f.getName() for f in forces]
    chosenforcename= forces[i].getName()
    print("forces:", forces)
    print("forcenames:", forcenames)
    print("Now looping over forcenames to delete")
    for fname in forcenames:
        if fname != chosenforcename:
            omm.remove_force_by_name(fname)
    # Remaining
    forces = omm.system.getForces()
    print("Forces remaining:", forces)
     # Singlepoint calculation
    result = Singlepoint(theory=omm, fragment=frag)
    energy_dict[f"System - {chosenforcename}"] = result.energy

print("energy_dict:", energy_dict)
print("")
print(" Forcename-cont   Energy")
print("----------------------------")
for name, e in energy_dict.items():
    print("{:10} {:13.10f}".format(name,e))

[RagnarB83]

Oh and I had realized that OpenMMTheory was missing options to remove forces so I just added this to the code: remove_force and remove_force_by_name

To use the script above, you have to checkout the latest ASH version in the NEW branch.

[0xBytefuzz]

Thank you very much for your response. Over the past few days, I've been busy with other matters and couldn't reply promptly. Although I have resolved the issue with the abnormal MM single-point energies, I am still looking forward to updates to ASH.

When calculating single-point energies, I used Amber's prmtop and rst files because this protein contains a non-standard residue (MIO), which I parameterized using antechamber. My simulations ran well, but when calculating single-point energies with ASH, the MM energies for some frames were positive. A visual inspection did not reveal any obvious structural errors. After I corrected for periodicity and centered the protein, the single-point energies calculated by ASH were normal. Based on my experience, the energy fluctuations were reasonable, and there were no sudden positive values.

Therefore, I suggest:

Are there any potential issues with periodic correction when ASH reads TOP files Should the documentation include a reminder to check for periodicity issues when reading intermediate frames from simulations Lastly, I wish you all the best.

[RagnarB83]

Ah I see, it was the periodicity. This will depend on what the input data looked like. If the frame coordinates were not centered before then the QM/MM will not work since the QM-part will normally not be periodic.

Documentation probably needs some updates there.