Appropriate forcefield/residue type for graphene?

[Andrew MacBride]

I'm doing some graphene simulations in water, and am having some setup problems. In this instance I'm starting from a SMILES string and using RDKit to add appropriate hydrogens, then converting to a PDB file. As a result, the PDB residue field is unknown and set to UNL. I tried using the SMIRNOFFTemplateGenerator from openmmforcefields, but have found that the system setup (forcefield.createSystem) fails. (Details: AmberTools antechamber kept allocating memory until it reached 70G RSS, at which point I killed it.)

Setup: OpenMM 8.0.0 installed via conda
openmmforcefields 0.13.0
RDKit 2023.03.3
Ubuntu 22.04 LTS, 64G RAM, 16G swap, 16 threads

Input SMILES (4x6 graphene fragment): C1=CC2=CC3=CC4=C5C6=C3C7=C2C(=C1)C8=C9C7=C1C6=C2C3=C6C1=C1C9=C(C=CC1=CC6=CC1=C3C3=C6C(=C1)C=C1C=CC=C7C1=C6C1=C6C3=C2C5=C2C6=C3C(=CC2=C4)C=CC2=C3C1=C7C=C2)C=C

System size: 708 atoms including solvent

Code: Select all

alog("creating simulation system...")
alog("length of pdb.positions: " + str(len(pdb.positions)))

system = forcefield.createSystem(
    pdb.topology,
    nonbondedMethod=PME,
    nonbondedCutoff=1*unit.nanometer,
)
alog("system complete")

Is there a better residue type to use for graphene or similar systems? Is there a prebuilt forcefield (AMBER, CHARMM, etc.) that should be selected instead of making a custom FF via SMIRNOFFTemplateGenerator?

Any tips/tricks/RTFMs are welcome -- much appreciated. (For example, would a different choice for nonbondedMethod be more appropriate?)

Thanks,
Andrew

[Peter Eastman]

It's probably trying to compute the atomic charges using AM1-BCC, which is going to be pretty expensive. I haven't worked with graphene, but I think there are force fields specifically for simulating it. A web search for "graphene molecular force field" turned up lots of hits, for example

https://link.springer.com/article/10.10 ... 018-2115-5
https://pubs.aip.org/aip/jap/article/11 ... perties-of
https://www.mdpi.com/2073-4360/6/9/2404

A specially designed force field is likely to be more accurate than a generic one like OpenFF, especially since the latter is really designed for small molecules.

[Andrew MacBride]

Great, thanks for the information.

I've found that discoverability of force fields is a bit rough -- if anyone has a specific suggestion of a FF I can use out-of-the-box with OpenMM (or convert from another system and drop in), I'd appreciate it. At this point, perfect accuracy isn't as important as just getting it to work (with this small model system.)

[Peter Eastman]

You could always install openmm-xtb and use GFN-FF. It's slower than a lot of other force field, and I can't speak for the accuracy, but it's very easy to get any system running.

This paper uses OpenMM to simulate graphene. I'm not sure if their code is available.

[Andrew MacBride]

Thank you for the tip -- I've had some success in using xtb directly, and now I'd like to use it from within OpenMM. The only example code I've found was the TestXtbForce class in the test suite for openmm-xtb, but I wasn't sure how to use the low-level XtbForce with higher-level OpenMM classes/methods such as ForceField.createSystem(), Simulation, etc.

Is there a simple example floating around somewhere that shows how to use XtbForce, or should I assume that I'll have to construct the System by individually adding particles, as in the unit test?

[Peter Eastman]

Just create the System directly.

Code: Select all

system = System()
for atom in topology.atoms():
    system.addParticle(atom.element.mass)
indices = list(range(system.getNumParticles()))
numbers = [atom.element.atomic_number for atom in topology.atoms()]
system.addForce(XtbForce(method, charge, multiplicity, periodic, indices, numbers))

We could probably make the xtb plugin install a force field to automate this. I'll open a feature request for that.