ash icon indicating copy to clipboard operation
ash copied to clipboard

Mechanical embedding in QMMM

Open quanp opened this issue 2 months ago • 8 comments

Hi @RagnarB83 ,

I just want to ask whether the mechanical embedding in QMMM MD is fully working.

A simple single point calculation gives me this

QM energy:             -1564.050764951681
MM energy:               -69.454757244004
Subtractive correction energy:        0.000000000000
QM/MM energy:          -1633.505522195685

I am a bit confused that the Subtractive correction energy value is 0.

Here is the input

from ash import *

pdbfile="warmup_MD_cycle2_lastframe_imaged.pdb"
frag = Fragment(pdbfile=pdbfile, charge=0, mult=1)
qmatoms = list(range(0,61))+[12151]

#Define QM, MM and QM/MM Theory
qm_theory = FairchemTheory(model_file="esen_sm_conserving_all.pt", device="cuda")
mm_theory = OpenMMTheory(xmlfiles=["LIG.xml", "amber/tip3p_standard.xml"], pdbfile=pdbfile, periodic=True, rigidwater=False, platform="CPU", autoconstraints=None)
qm_mm_theory = QMMMTheory(qm_theory=qm_theory, mm_theory = mm_theory, qmatoms=qmatoms, fragment=frag, qm_charge=1, qm_mult=1, embedding="mech")

Singlepoint(theory=qm_mm_theory, fragment=frag)

Thank you very much

quanp avatar Nov 12 '25 15:11 quanp

Hi quanp

Yes, it should be. The subtractive correction energy is an energy contribution only used for a special case in the QM/MM module (that has not yet been implemented actually). Set to 0.0 by default. The QM/MM module in ASH assumes an additive energy expression in general. For mechanical embedding, all the QM-MM coupling terms are calculated by the MM program.

RagnarB83 avatar Nov 13 '25 08:11 RagnarB83

Thanks you

I am still a bit confused. If I do the Singlepoint calculation with mm_theory only, I get an exact energy as the MM energy in the QMMM calculation (-69.454757244004 au). Does it mean the QMMM coupling terms are 0?

quanp avatar Nov 17 '25 06:11 quanp

Hi Quan Sorry for the delay, I will be looking into this today.

RagnarB83 avatar Nov 18 '25 12:11 RagnarB83

  1. So, one thing to be aware of is that if you define a QMMMTheory object using an OpenMMTheory object then the MM-system is modified (during the initialization of the QMMMTheory object) to account for QM/MM specific interactions. So defining both theories in the same script means that the MM-energy from OpenMMTheory will be modified. This could explain why you get the same exact MM energy from a Singlepoint calculation, if both theories were defined in the same script.

  2. I did spot a bug in the way the QM-MM interactions are treated in mechanical embedding. Had not spotted it before because mechanical embedding is relatively new in ASH and has not been used at all in my group. The fix has been added to the NEW branch. I still have not fully confirmed the correctness of the mechanical embedding against another code. Will be done hopefully today and tomorrow. Until then, be aware.

Thanks for bringing it up.

RagnarB83 avatar Nov 20 '25 10:11 RagnarB83

Thank you very much. Indeed, it was my mistake to keep the qm_mm_theory line while doing pure MM.

Another question: Is doing electrostatic QMMM calculations in combination with FairChem valid? Obviously, FairChem can not provide QM energy with polarization by MM point charges.

Another note, since I am testing the Fairchem interface. I notice that FAIRChemCalculator (and model) is created every time .run() is called, which is very inefficient for MD or geometry optimization loops. It takes a few seconds for FairChem to load the model. Should it be better to initialize the Fairchem model once in init(), and reuse it in run()?

Best regards

quanp avatar Nov 20 '25 15:11 quanp

Indeed, the electrostatic embedding QM/MM approach is not compatible with ML/MM. Mechanical embedding is hence required and some sensible charges defined for the ML region required. ONIOM is another applicable approach.

Good point about the FairChem interface. It was hastily written and meant to be revisited. I will try to fix the loading problems soon.

RagnarB83 avatar Nov 20 '25 15:11 RagnarB83

Thanks My goal is to run MD, but it seems that ONIOMTheory and OpenMM_MD are not compatible. print("Problem: ONIOMTheory containing an OpenMMTheory is currently not supported yet. Complain to developer")

quanp avatar Nov 20 '25 16:11 quanp

The Fairchem interface has been improved. Should be no significant interface overhead present any more.

I will check out the ONIOM-with-OpenMM MD limitation, had forgotten about that.

RagnarB83 avatar Nov 21 '25 05:11 RagnarB83

Hi Quan

Update: I can confirm that the mechanical embedding QM/MM implementation in ASH is fully correct on the NEW branch since Nov 20th. Confirmed by testing against another code.

The limitation of MD not being possible with ONIOM containing an OpenMMTheory, will be addressed in the next few days.

RagnarB83 avatar Nov 27 '25 19:11 RagnarB83

The "ONIOMTheory-for-MD" limitation (when ONIOM Theoryincludes an OpenMMTheory) has been lifted. The NEW branch now contains a fix for this.

There is a small limitation, however, in that when defining ONIOM regions, Region 1 definition must include a whole residue compatible with the forcefield defined by OpenMM. This should only really be a problem when using ONIOM for systems where there is a covalent QM-MM boundary. It's a tough problem to solve, so this will be a limitation for the time being. But not an issue if region 1 includes e.g. a whole molecule/ligand (and a forcefield available for the whole molecule).

ONIOM functionality will be improved in future versions of ASH.

RagnarB83 avatar Nov 28 '25 15:11 RagnarB83

Closing this as things should have been fixed.

RagnarB83 avatar Nov 30 '25 12:11 RagnarB83