Periodic boundary conditions with DrudeForce

[frgr2254]

Hi, I have a question regarding periodic boundary conditions and the DrudeForce. From the documentation for DrudeForce I understand that drudeforce.usesPeriodicBoundaryConditions() will return True or False depending on if the NonbondedForce uses PBC. However in my case drudeforce.usesPeriodicBoundaryConditions() returns False although nonbondedforce.usesPeriodicBoundaryConditions() returns True.

Since there is no drudeforce.setUsesPeriodicBoundaryConditions(True) how can I set drudeforce.usesPeriodicBoundaryConditions() to True?

I attach my code:

# Import openmm
from openmm.app import *
from openmm import *
from openmm.unit import *
from sys import stdout


# Read force field from XML file

forcefield = ForceField("forcefield.xml")


# Read positions and topology from Gromacs files, create a modeller object and add extra particles

gro = GromacsGroFile("slab.gro")

top = GromacsTopFile("topol.top", periodicBoxVectors=gro.getPeriodicBoxVectors(), includeDir='/usr/local/gromacs/share/gromacs/top') #gromacs top file

modeller = Modeller(top.topology,gro.positions)

modeller.addExtraParticles(forcefield)


# Create a system

system = forcefield.createSystem(modeller.topology, nonbondedMethod=PME, nonbondedCutoff=1.0*nanometer, constraints=HBonds)


# Activate PBC for bonds and angles

forces = system.getForces()

for force in forces:

    force_name = force.__class__.__name__

    if force_name == "HarmonicBondForce":

        force.setUsesPeriodicBoundaryConditions(True)

    elif force_name == "HarmonicAngleForce":

        force.setUsesPeriodicBoundaryConditions(True)


# Check options for periodic boundary conditions

for force in forces:

    force_name = force.__class__.__name__

    if force.usesPeriodicBoundaryConditions():
        
        print(f"{force_name} uses PBC")

    else: 
        print(f"{force_name} does not use PBC")

The output is

HarmonicBondForce uses PBC
NonbondedForce uses PBC
DrudeForce does not use PBC
CMMotionRemover does not use PBC
HarmonicAngleForce uses PBC

[peastman]

That documentation is a bit confusing. I'll fix it.

DrudeForce never uses periodic boundary conditions. It's a reasonable question whether it should have an option for it. The force connecting a Drude particle to its parent clearly should never use periodic boundary conditions. That just wouldn't make sense. But the interaction between screened pairs could make sense, if you're simulating an infinite molecule that wraps around.

[frgr2254]

Thank you very much for clarifying this. I apologize if I misunderstand how the program is working and if my suggestions do not make sense. I am developing a polarizable force field for metal oxides based on the Drude oscillator model. My goal is to simulate surfaces of these crystals (periodically repeated in the xy-plane) and study their interaction with water as well as biomolecules. In this situation I must specify force.setUsesPeriodicBoundaryConditions(True) for both HarmonicBondForce and HarmonicAngleForce in order for these interactions to work correctly across the boundary. I am of the opinion that the Thole-screened electrostatic interaction should handle periodic boundary conditions in exactly the same way as HarmonicBondForce and HarmonicAngleForce.

I will try to elaborate a bit what I mean. For simplicity consider a one-dimensional crystal with atoms A, B, C, ... each represented by a drude-core pair and A is bonded to B that is bonded to C and so on (see attached PNG figure). By A0, B0, C0, ... I denote particles in the central box and by A1, B1, C1, ... and A-1, B-1, C-1, ... their periodic images in adjacent boxes. In this situation I think that A0 (core and drude) should experience the Thole-interaction from F-1 and E-1 under periodic boundary conditions (in addition to its interaction with B0 and C0). B0 should experience the interaction from F-1 under periodic boundary conditions (in addition to its interaction with A0, C0 and D0). Analogously F0 should experience the interaction from A1 and B1 under periodic conditions (in addition to its interaction with E0 and D0). E0 should experience the interaction with A1 under periodic boundary conditions (in addition to its interaction with F0, D0 and C0).

[peastman]

That makes sense. I've marked this as a feature request. I should be able to add it soon.

[peastman]

The implementation is in #4523. If you're interested in testing it out, please let me know whether it works for you.