Permanent monopole energies: Nonbonded x AmoebaMultipoleForce

[gcarneiroq]

Dear all,

I am currently testing some possibilities of modeling a periodic bulk structure of quartz. The topology, in principle, makes use of bond across images.

I would like to evaluate the possibility to use the AMOEBA plugin to calculate the induced polarization interactions, keeping only the partial charges of the atoms (excluding dipoles and quadrupoles).

I confess I am having a hard time on understanding the scaling of the permanent electrostatic interactions. I performed two tests (polarizabilities are zeroed for the moment)

1. Ignoring bonds - turning off the bonded interactions the energy of nonbondedforce (without LJ) meet the result of the AmoebaMultipoleForce with both considering PME.

2. For trying to understand the scaling factors, I set the covalentmap(iatom, 2 ) and covalentmap(iatom, 3) to empty lists, in that I was hoping to recover full 1-5 and 1-4 permanent electrostatic interactions. However, I see a deviation of about 10% in the AmoebaMultipoleForce and the nonbondedforce. In this test I used the harmonicbondforce with PBC.

My question now is if there is some scaling involving the 1-3 interactions or just by setting the aforementioned terms of the covalentmaps is not enough to recover the full 1-4 and 1-5 permament electrostatics?

Thank you very much for any information.

Kind regards,

Guilherme

PS. Below is a part of the python script.

for force in system.getForces(): if isinstance(force, HarmonicBondForce): force.setUsesPeriodicBoundaryConditions(True) elif isinstance(force, HarmonicAngleForce): force.setUsesPeriodicBoundaryConditions(True)
force.setCutoffDistance(1.2nanometer) force.setNonbondedMethod(AmoebaMultipoleForce.PME) print(force.usesPeriodicBoundaryConditions()) for i in range(0,n_atoms): cov_map = force.getCovalentMaps(i) for j in cov_map[0]: if [i,j] and [j,i] not in exceptions_12: exceptions_12.append([i,j]) for j in cov_map[1]: if [i,j] and [j,i] not in exceptions_13: exceptions_13.append([i,j]) for j in cov_map[2]: if [i,j] and [j,i] not in exceptions_14: exceptions_14.append([i,j]) for j in cov_map[3]: if [i,j] and [j,i] not in exceptions_15: exceptions_15.append([i,j]) force.setCovalentMap(i,3,[]) force.setCovalentMap(i,2,[]) elif isinstance(force, CustomNonbondedForce): force.setNonbondedMethod(CustomNonbondedForce.CutoffPeriodic) force.setCutoffDistance(1.2nanometer) force.setUseLongRangeCorrection(True) else: pass

forcegroups = forcegroupify(system)

for i, f in enumerate(system.getForces()): f.setForceGroup(i)

integrator = LangevinIntegrator(

system temperature:

            temperature, 
            # collision rate
            1.0/picoseconds, 
            # timestep:
            0.5*femtoseconds)

integrator.setConstraintTolerance(0.00001) if npt:

Add a bartostat to NPT simulations

system.addForce(
        MonteCarloBarostat(pressure, temperature, 25))

platform = Platform.getPlatformByName('CUDA') properties={}

simulation = Simulation(model_topology, system, integrator, platform,properties)

simulation.context.setPositions(model_positions)

[gcarneiroq]

I performed another test, the energies of test2 met when no Long range corrections are applied. It might be the case that another parameter controls the scaling of the long range correction?

[peastman]

Can you edit your post and add lines containing ``` before and after the code? It's getting rendered as markdown, which makes it hard to interpret.

I don't understand what your code is meant to do. How did you create the System in the first place? Do you expect all the covalent maps to be correct at the beginning? You have variables called exceptions_12, exceptions_13, etc. which apparently already exist, and then you add more entries to them, but never use them for anything. What is that for? And your code checks the type of force, looking for HarmonicBondForce and HarmonicAngleForce. And then you start treating it as an AmoebaMultipoleForce without ever having checked for that class.

[gcarneiroq]

Hi, thank you for looking on that!

I apologize for the lousy way I wrote. Here is the part of the code in a more clean version:

for force in system.getForces():
    if isinstance(force, HarmonicBondForce):
        force.setUsesPeriodicBoundaryConditions(True)
    elif isinstance(force, HarmonicAngleForce):
        force.setUsesPeriodicBoundaryConditions(True)  
    elif isinstance(force, AmoebaMultipoleForce):
        force.setCutoffDistance(1.2*nanometer)
        print(force.usesPeriodicBoundaryConditions())
        for i in range(0,n_atoms):
            force.setCovalentMap(i,3,[])
            force.setCovalentMap(i,2,[])
            force.setCovalentMap(i,7,[])
        force.setNonbondedMethod(AmoebaMultipoleForce.PME)
    elif isinstance(force, CustomNonbondedForce):
        force.setNonbondedMethod(CustomNonbondedForce.CutoffPeriodic)
        force.setCutoffDistance(1.2*nanometer)
        force.setUseLongRangeCorrection(True)
    else:
        pass

I want to compare the electrostatics calculated by the Nonbondedforce and by the Amoebamultipoleforce (with only the partial charges defined). Later I will turn on the polarization.

The point is that I would like to have only the 1-2 and 1-3 interactions ignored by the by the Nonbondedforce and the amoebamultipoleforce.

So i changed the covalentmaps for 14 and 15 to empty lists.

By doing that I manage to reproduce the exact same electrostatics without PME, that is, using a plain cutoff.

In another test I ignored the bonds and angles, that is I use a topology without any connections. In that case the the results coincide with and without PME.

Regards

[peastman]

You're saying that if you don't have any bonds, they agree perfectly, but if you do have bonds they don't?

If so, presumably it's a difference in the treatment of bonded interactions. You didn't answer the question about how you created the System. What exceptions does it contain? Typically a force field will omit all 1-2 and 1-3 interactions, and reduce the strength of 1-4 interactions. Is that what yours does?

AMOEBA is different in that it treats bonded pairs differently in the polarization calculation, hence the different maps for Covalent and PolarizationCovalent. How do you have those set up?

[gcarneiroq]

Not only that. If I do define the bonds, but do not make use of PME the energies match as well.

The topology is defined in the attached xml file.

In the script here are the lines defining the system:

import mbuild as mb
import MDAnalysis as mda
from openmm import *
from openmm.app import *
from openmm.unit import *

mol_file = '/home/gui/Documents/simulations/openmm/quartz/build/quartz_pb.mol2'

q = mb.load(mol_file)
n_atoms = len(q.xyz)
q_u = mda.Universe(mol_file)

topology = Topology()
topology.setPeriodicBoxVectors((u,v,w))
chain = topology.addChain()
residue = topology.addResidue("QUA",chain)
element_Si = Element.getByAtomicNumber(14)
element_O = Element.getByAtomicNumber(8)

for i in range(0,n_atoms):
    if q[i].name == 'Si':
        globals()['atom'+str(i)] = topology.addAtom("Si"+str(i), element_Si, residue)
    else:
        globals()['atom'+str(i)] = topology.addAtom("O"+str(i), element_O, residue)

for i in q_u.bonds:
       topology.addBond(globals()['atom'+str(i[0].index)],globals()['atom'+str(i[1].index)])

forcefield = ForceField('./bonds.xml')
model = modeller.Modeller(topology,q.xyz)
model_topology = model.getTopology()
model_topology.setPeriodicBoxVectors((u,v,w))
model_positions = model.getPositions()

system = forcefield.createSystem(model_topology, nonbondedCutoff=1.2*nanometer,constraints='None',removeCMMotion=True,nonbondedMethod=PME,polarization='mutual',mutualInducedTargetEpsilon=0.00001,ewaldErrorTolerance=1e-5)

And then goes to the other piece of code I pasted before.

The bondcutoff in XML variable is 2 so I need to empty the 13, 14 and 17 covalent maps as written above.

In that case, the energies of Nonbonded and AmoebaMultipole matches when no PME is evoked.

Here is the big chunk of the xml file (I tried to upload but it didnt accept the format, and I ommitted all the atoms and bonds, but they are defined in the original file )

<ForceField>

<AtomTypes>
 <Type name="1" class="Si" element="Si" mass="28.0855"/>
 <Type name="2" class="O" element="O" mass="15.999"/>
</AtomTypes>

<Residues>
 <Residue name="QUA">
  <Atom name="Si0" type="1"/>
  <Atom name="O1" type="2"/>
  ....
 <Bond atomName1="Si0" atomName2="O3712"/>
...
</Residue>
</Residues>

<HarmonicBondForce>
<Bond class1="Si" class2="O" length="0.163662" k="261578.0375"/>
</HarmonicBondForce>

<HarmonicAngleForce>
<Angle class1="O" class2="Si" class3="O" angle="1.9111334947163243" k="681.206083333333"/>
<Angle class1="Si" class2="O" class3="Si" angle="2.2406152964667765" k="686.179"/>
</HarmonicAngleForce>

<AmoebaMultipoleForce direct11Scale="0.0"  direct12Scale="1.0"  direct13Scale="1.0"  direct14Scale="1.0"  mpole12Scale="0.0"  mpole13Scale="0.0"  mpole14Scale="0.4"  mpole15Scale="0.8"  mutual11Scale="1.0"  mutual12Scale="1.0"  mutual13Scale="1.0"  mutual14Scale="1.0"  polar12Scale="0.0"  polar13Scale="0.0"  polar14Intra="0.5"  polar14Scale="1.0"  polar15Scale="1.0" >
<Multipole type="1" kz="0" kx="0" c0="1.864536" d1="0" d2="0" d3="0" q11="0" q21="0" q22="0" q31="0" q32="0" q33="0"/>
<Multipole type="2" kz="0" kx="0" c0="-0.932268" d1="0" d2="0" d3="0" q11="0" q21="0" q22="0" q31="0" q32="0" q33="0"/>
 <Polarize type="1" polarizability="0.0" dampingFactor="1." thole="0.33" />
 <Polarize type="2" polarizability="0.0" dampingFactor="1." thole="0.33" />
</AmoebaMultipoleForce>

</ForceField>

[peastman]

Can you put everything needed to reproduce it (complete scripts and input files) into a zip file and post it?

[gcarneiroq]

Sure, no problem. I attached a zip files with different folders corresponding to each test. For the moment, I am only interested in matching the electrostatic energies from the Nonbonded with the Amoebamultipole, so I still need to deal with exceptions of the customnonbondedforce, once the electrostatics is settled.

The reason behind all that is that i want to simulate a system with a customnonbonded force with the 1-4 interactions fully on (both dispersion and electrostatics).

Thank you very much for taking time to dig into it. clean_scripts.zip

[peastman]

The files I should be looking at are the ones in the folder amoeba_bonds_pme, right?

[gcarneiroq]

Yes, precisely. This is the one i am interested and which I dont understand the energy value.

Em ter., 29 de out. de 2024 00:37, Peter Eastman @.***> escreveu:

The files I should be looking at are the ones in the folder amoeba_bonds_pme, right?

— Reply to this email directly, view it on GitHub https://github.com/openmm/openmm/issues/4705#issuecomment-2442872014, or unsubscribe https://github.com/notifications/unsubscribe-auth/AFNTMN6HX7EMTCOOXHOVYDTZ53DFDAVCNFSM6AAAAABQRNKYFGVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDINBSHA3TEMBRGQ . You are receiving this because you authored the thread.Message ID: @.***>

[peastman]

What is mbuild? Your script imports it. I'm not familiar with it.

[gcarneiroq]

Mbuild is a library within the mosdef hub for building systems and topology, this was the way i found to create topology files with periodic bonds across the boxes.

https://mbuild.mosdef.org/en/stable/

Perhaps if i save a pdbx file we could run without the need of regenerating the topology all the time.

Em ter., 29 de out. de 2024 00:46, Peter Eastman @.***> escreveu:

What is mbuild? Your script imports it. I'm not familiar with it.

— Reply to this email directly, view it on GitHub https://github.com/openmm/openmm/issues/4705#issuecomment-2442880979, or unsubscribe https://github.com/notifications/unsubscribe-auth/AFNTMN4IBCTWUCPTQJQDKX3Z53EGJAVCNFSM6AAAAABQRNKYFGVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDINBSHA4DAOJXHE . You are receiving this because you authored the thread.Message ID: @.***>

[peastman]

That would be great if you can do that.

[gcarneiroq]

amoeba_pme_bonds_pdbx.zip

Here it is! I tested here and gives the same result as previously building the topology using mbuild

[peastman]

Got it, thanks. And the question is why the energy of the AmoebaMultipoleForce is different from the energy of the NonbondedForce in coul_bonds_pme_with_14. Is that right?

[gcarneiroq]

Exactly

Em ter., 29 de out. de 2024 05:41, Peter Eastman @.***> escreveu:

Got it, thanks. And the question is why the energy of the AmoebaMultipoleForce is different from the energy of the NonbondedForce in coul_bonds_pme_with_14. Is that right?

— Reply to this email directly, view it on GitHub https://github.com/openmm/openmm/issues/4705#issuecomment-2443197291, or unsubscribe https://github.com/notifications/unsubscribe-auth/AFNTMN4IGUIB2SXCYMNDSY3Z54GV3AVCNFSM6AAAAABQRNKYFGVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDINBTGE4TOMRZGE . You are receiving this because you authored the thread.Message ID: @.***>

[peastman]

I think I understand what the difference is, but not why it's different. So I think this is a question for @jayponder.

Jay, the goal here is to get the AMOEBA code to reproduce a conventional point charge force field. If you set all the permanent dipoles and quadrupoles to 0, and likewise polarizability to 0 for every atom, the result should be identical to a standard point charge coulomb. And that works if you don't use PME, or alternatively if your system contains no bonds. But if you use PME and the system contains bonds, they don't agree.

The difference comes from the handling of excluded interactions, e.g. 1-2 interactions. They should be completely omitted, but there's no way to prevent the reciprocal space energy from including a contribution from them. Since that energy is incorrectly included in the reciprocal space calculation, we compensate by subtracting it back off in direct space.

Anyway, that's what we do for standard Coulomb interactions. But the AMOEBA implementation in OpenMM doesn't do that. And I searched through the Tinker code, and can't find any sign that it does it either. Should it be subtracting off the reciprocal space part of excluded interactions? If not, then why not?

[jayponder]

Hi Everyone, I'm not sure what the OpenMM plugin may be doing, but I believe the test Peter describes immediately above works correctly in the Tinker codes (from github.com/TinkerTools).

For years I've kept files in the /example directory of the Tinker distribution that implement exactly this kind of test. See the "formbig" test case, which is a large monoclinic chunk of crystalline formamide. This test uses PME and computes the energy for simple OPLS-AA atomic charges, and for AMOEBA-like "multipoles" where the charges are the same as OPLS-AA and all dipoles, quadrupoles and polarizabilities are zero. For this test the total simple atomic charge energy and total "multipole" energy are exactly the same- as they should be.

In contrast to what Peter suggests, I believe the Tinker code does remove the reciprocal part of the excluded interactions in the real space code. This is done via the mscale() array and scalek variable in the real space routines in "empole.f" and similar source files in Tinker. Perhaps, as Peter also suggests, the OpenMM plugin is not doing this removal.

I'm happy to provide more details and/or to help in any way I can.

[peastman]

Do you know if that behavior changed in Tinker at some point? Mark created the original OpenMM AMOEBA imlementation based on the Tinker code back in 2010. I looked in the Tinker source repository, but the earliest version of the code there is from 2014.

[jayponder]

Well, the various Tinker code bases have certainly been reorganized multiple times over the years :)

I was able to dig up an old version of "Tinker5" from February 2010, and it has a "scaling" mechanism in the real space multipole routines to remove or scale excluded interactions similar to the current Tinker8. Both this old Tinker5 and the Tinker8 CPU code do the removal as an integral part of the loop over all real space interactions. As I'm sure you are aware, it is also possible to do the exclusions in a separate small loop over just the excluded and scaled interactions. This latter method is likely preferred on GPUs.

I believe all versions of Tinker have always had some method for doing this correction for stuff included in reciprocal space. If you don't do so then the answers you get are simply wrong. As I say above, we've had a test for agreement between simple partial charge electrostatics and charge-only "polarizable multipoles" since the beginning days of AMOEBA- it was a very early sanity check...

[peastman]

Agreed. We do have test cases that create a variety of systems, including ones with both PME and bonds, and compare the forces and energies to results computed with Tinker. I need to dig into it further.

[peastman]

I computed the energy of a box of AMOEBA water with PME with both Tinker and OpenMM. They give perfect agreement for all components. Is that a sufficient test that it's handling this correctly, or are there other cases I should check?

[jayponder]

Hi Peter, AMOEBA completely neglects 1-2 and 1-3 permanent multipole interactions. But both 1-4 and 1-5 interactions are scaled, by 0.4 and 0.8 respectively. Similarly, there is partial exclusion ("scaling") of intramolecular 1-4 polarization. So it might be good to check something bigger than water- maybe a protein in a water box.