Extra point Sites on OPC water

[pramdhan]

Hello,

I am trying to input my parameters from AMBER into ASH. I used a ff19sb/opc water combination, however, when it attempts to be read by ash, it does not recognize the extra points from EPW placed. Is this a common issue? Should I delete it and continue forward without EPW?

[RagnarB83]

4-point water models are indeed problematic because of the extra site. We've not really used them before in our studies but I'm happy to try to have them supported.
I've just made some small changes to ASH Fragment and OpenMMTheory so that the virtual sites do not cause a crash anymore.

This basic ASH script:

from ash import *

frag = Fragment(pdbfile="water_with_vp.pdb")
theory = OpenMMTheory(pdbfile="water_with_vp.pdb", xmlfiles=["opc.xml"], autoconstraints=None, rigidwater=False)
Singlepoint(theory=theory, fragment=frag)

where PDB-file is:

REMARK   1 CREATED WITH OPENMM 8.3, 2026-06-10
HETATM    1  O   MOL A   1       0.225  -0.836  -1.476  1.00  0.00           O
HETATM    2  H1  MOL A   1       0.225  -0.065  -0.907  1.00  0.00           H
HETATM    3  H2  MOL A   1       0.225  -1.607  -0.907  1.00  0.00           H
HETATM    4 M    MOL A   1       0.225  -0.836  -1.308  1.00  0.00          EP
TER       5      MOL A   1
CONECT    1    2    3
CONECT    2    1
CONECT    3    1
END

does seem to work now.

Check out the NEW branch if you want to try it out.

[pramdhan]

Cool, thanks for doing that! Do you recommend uploading my qm/mm snapshot from amber as a PDB? I'm trying to use frag=Fragment(amber_prmtopfile=prmtopfile, amber_inpcrdfile=inpcrdfile), however it gets stuck with

ASH Fragment creation
Reading Amber INPCRD file
Reading coordinates from Amber INPCRD file: 'wt-clop-rc.ncrst' and PRMTOP file: '2c19-clop.parm7' into fragment.
Numatoms: 55557
Box dimensions read: [82.3731987, 82.3731987, 82.3731987, 109.471219, 109.471219, 109.471219]
Traceback (most recent call last):
File "/blue/lic/pramdhan1/benchmark_dft/clopidogrel/wt-rc/test.py", line 8, in
frag=Fragment(amber_prmtopfile=prmtopfile, amber_inpcrdfile=inpcrdfile)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/blue/lic/pramdhan1/.conda/envs/ash_amber/lib/python3.12/site-packages/ash/modules/module_coords.py", line 239, in init
self.read_amberfile(inpcrdfile=amber_inpcrdfile, prmtopfile=amber_prmtopfile, conncalc=conncalc)
File "/blue/lic/pramdhan1/.conda/envs/ash_amber/lib/python3.12/site-packages/ash/modules/module_coords.py", line 512, in read_amberfile
elems, coords, box_dims = read_ambercoordinates(prmtopfile=prmtopfile, inpcrdfile=inpcrdfile)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/blue/lic/pramdhan1/.conda/envs/ash_amber/lib/python3.12/site-packages/ash/modules/module_coords.py", line 2367, in read_ambercoordinates
elems += [reformat_element(int(i), isatomnum=True) for i in line.split()]
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/blue/lic/pramdhan1/.conda/envs/ash_amber/lib/python3.12/site-packages/ash/modules/module_coords.py", line 1199, in reformat_element
el_correct = ash.dictionaries_lists.element_dict_atnum[elem].symbol
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^^^^^^
KeyError: 0

[pramdhan]

When checking parmed, it appears that EPW is atom 0 so I'm unsure if this is just a naming issue where ash would be able to recognize it if changed? Not sure if you ran into an issue like this.

printDetails @0,1

The mask @0,1 matches 2 atoms:

ATOM RES RESNAME NAME TYPE At.# LJ Radius LJ Depth Mass Charge GB Radius GB Screen
1 1 LEU N N3 7 1.8240 0.1700 14.0100 0.1010 1.5500 0.7900
55557 12509 WAT EPW EP 0 0.0000 0.0000 0.0000 -1.3583 1.5000 0.8000

[RagnarB83]

The amber-parsing probably needs to be tweaked as well.
If you can send me the files I can probably fix this (ragnar.bjornsson@gmail.com).

[pramdhan]

I was playing around with the definitions but possible fix to it? It maintains EPW while being able to take in parm7 and ncrst files

==============================================================================

Virtual site patches (do not edit)

==============================================================================

import os
import sys
import parmed as pmd
from parmed.amber import AmberMask

import ash
from ash.dictionaries_lists import Element
import ash.modules.module_coords

ash.dictionaries_lists.element_dict_atnum[0] = Element("Dummy", "EP", 0)
ash.dictionaries_lists.element_dict_atname["ep"] = Element("Dummy", "EP", 0)
ash.dictionaries_lists.element_dict_atname["EP"] = Element("Dummy", "EP", 0)
ash.modules.module_coords.elematomnumbers["ep"] = 0
ash.modules.module_coords.elematomnumbers["EP"] = 0
ash.modules.module_coords.eldict_covrad["ep"] = 0.0
ash.modules.module_coords.eldict_covrad["EP"] = 0.0

from ash import *

def amber_mask_to_ash_indices(prmtopfile, inpcrdfile, mask):
"""Convert Amber mask to ASH 0-based indices, excluding virtual sites."""
struct = pmd.load_file(prmtopfile)
struct.load_rst7(inpcrdfile)
amask = AmberMask(struct, mask)
return sorted(
idx for idx, sel in enumerate(amask.Selection())
if sel == 1 and struct.atoms[idx].element_name.strip().upper() != "EP"
)

==============================================================================

USER SETTINGS - edit this section

==============================================================================

prmtopfile = "YOUR_PRMTOP.parm7" # Amber topology file
inpcrdfile = "YOUR_COORDS.ncrst" # Amber coordinate/restart file
qm_charge = 0 # Total charge of QM region
qm_mult = 1 # Spin multiplicity of QM region
numcores = 1 # Number of CPU cores for ORCA

qmmask = 'YOUR_QM_AMBER_MASK' # Amber mask for QM region
# e.g. ':LIG|:FE1|:RES@CA,CB'

actmask = 'YOUR_ACTIVE_AMBER_MASK' # Amber mask for active (mobile) region
# e.g. '(:LIG|:FE1)<:5.0'

orcasimpleinput = '! YOUR_METHOD YOUR_BASIS' # e.g. '! r2SCAN-3c' or '! XTB2'
orcablocks = """%maxcore 4000
%scf
MaxIter 1000
end
"""

==============================================================================

qmatoms = amber_mask_to_ash_indices(prmtopfile, inpcrdfile, qmmask)
actatoms = amber_mask_to_ash_indices(prmtopfile, inpcrdfile, actmask)
print(f"QM atoms ({len(qmatoms)}): {qmatoms}")
print(f"Active atoms ({len(actatoms)}): {actatoms}")

sys.stderr = open(os.devnull, 'w')
frag = Fragment(amber_prmtopfile=prmtopfile, amber_inpcrdfile=inpcrdfile)
sys.stderr = sys.stderr

omm = OpenMMTheory(
Amberfiles=True, amberprmtopfile=prmtopfile,
printlevel=1, periodic=True,
autoconstraints=None, rigidwater=False
)

orca = ORCATheory(
orcasimpleinput=orcasimpleinput,
orcablocks=orcablocks,
numcores=numcores
)

qmmm = QMMMTheory(
qm_theory=orca, mm_theory=omm,
fragment=frag, embedding="Elstat",
qmatoms=qmatoms, printlevel=2
)

Optimizer(
fragment=frag, theory=qmmm,
ActiveRegion=True, actatoms=actatoms,
maxiter=200, charge=qm_charge, mult=qm_mult
)

frag.write_xyzfile(xyzfilename="qmmm_optimized.xyz")
write_pdbfile(frag, outputname="qmmm_optimized.pdb", openmmobject=omm)

[RagnarB83]

I've checked in a fix to the NEW branch that hopefully takes care of virtual sites when present in Amber files. Worked for my simple test system so you can give it a go. If it does not work then I need your prmtop and coordinate files to further debug.

[pramdhan]

Sure thing, I'll email you.