search

openforcefield / openff-toolkit

The Open Forcefield Toolkit provides implementations of the SMIRNOFF format, parameterization engine, and other tools. Documentation available at http://open-forcefield-toolkit.readthedocs.io
http://openforcefield.org
MIT License
313 stars 92 forks source link

Fluctuating charges #556

Open matteo-ambrosetti opened 4 years ago

matteo-ambrosetti commented 4 years ago

Hello everyone, is it possible to change the atomic charges on the fly (or is there an easy way to do it)?

j-wags commented 4 years ago

Hi @matteo-ambrosetti,

I'm not sure exactly what you're looking to do. Could you describe your question a bit more?

A few general ways that we interact with charges are

Let me know if either of these solve your problem, otherwise we can discuss other possibilities for your specific needs.

matteo-ambrosetti commented 4 years ago

Yes, sorry, it was a quite dry question...

I was referring to the fluctuating charge force field (see https://aip.scitation.org/doi/pdf/10.1063/1.468398). The main difference with respect to a "standard" MD is in the ability of charges to fluctuate as a response to the environment. I have already implemented a code to compute the equilibrium charges and I was thinking to interface it with OpenMM. Since the atomic charges should change according to the system's geometry during the dynamics, I need to update them every time step.

j-wags commented 4 years ago

Ah, I see. Right now the Open Force Field Toolkit does not support parameterizing systems with fluctuating charges. In the future, we may begin offering support for certain polarizability models that are supported by OpenMM. This isn't on our development roadmap yet, so I can not provide a date for when that will be ready. If you're interested in adding it, I'd be happy to guide you through its implementation, though it would probably take a few weeks.

matteo-ambrosetti commented 4 years ago

That would be great, I am interested in adding it! I have been using the fluctuating charge (FQ) force field in the context of QM/MM models so far (relying on fixed charge FFs for the sampling part). I think that being able to perform a FQ-MD would be highly beneficial to my research.

davidlmobley commented 4 years ago

@matteo-ambrosetti we have no objection in principle to including FlucQ models; a big picture goal is to enable our toolkit/suite to support all the different types of force fields people might be interested in using so that they all can be parameterized/compared on equal footing.

[Note: I updated the above; I left out the word "no" in "no objection" in the first post in error.]

I think critical milestones to achieving this would be things like:

  1. Ensuring the science you want is supported in underlying simulation packages (esp. OpenMM but also AMBER, GROMACS...)
  2. Working with our devs/group to define how this would appear in the SMIRNOFF spec
  3. Planning (w devs) how to implement into the toolkit
  4. Implementing/getting code review etc.

I hope that helps!

matteo-ambrosetti commented 4 years ago

Thanks @davidlmobley, I will keep your suggestion in mind. This is the first time I deal with a code like this, so suggestions and help are welcomed! I will start by checking if the science needed is supported (I know for sure that CHARMM already supports FQ).

j-wags commented 4 years ago

Sounds good, @matteo-ambrosetti.

A bit of introduction for @davidlmobley's item 1: Right now we natively produce parameterized systems in OpenMM's System format, mostly because it has the highest information content of any MM file format (this lets us convert to other formats via ParmEd easily, since we're almost always removing information, rather than needing to guess it). Architecturally, this means that our parameterization routine is performed by the ForceField.create_openmm_system function. So, the most important thing is that we be able to write fluctuating charges into an OpenMM System object. The first step to that is to ensure that we can implement the fluctuating charge model into an OpenMM Force object (perhaps a CustomNonbondedForce offers enough flexibility to implement fluctuating charges?).

The step after that will be seeing if it'll be possible to convert this fluctuating charge model from an OpenMM force into other packages (like CHARMM). We usually handle this using ParmEd, but since this will be a new functional form, I don't think that ParmEd supports this conversion yet. So we'll need to add this if we want to simulate outside OpenMM.

Anyway, I'd be interested to know if the fluctuating charges model can be computed in OpenMM -- Please update us as you look into that!

davidlmobley commented 4 years ago

Right, an addendum to Jeff's point is that it's also important to know whether the CHARMM, AMBER and OpenMM packages can all support the SAME type of FlucQ model that you want to implement, so that moving it into all packages would just be an issue of conversion. Otherwise, broad support gets more complicated.

That wouldn't mean you couldn't proceed, but it would just affect how one would go about it.

matteo-ambrosetti commented 4 years ago

CHARMM implementation

From what I got from literature the Fluctuating Charge (FQ) Force Field (FF) is implemented only in CHARMM. The principal differences with respect to a non-polarizable molecular dynamics are:

  1. Electrostatic interaction

    • 1-2, 1-3, 1-4 interactions are computed with a screened Coulomb potential which depends on both atomic electronegativity and chemical hardness (see eq. (3)).
    • Atoms in a molecule separated by more than three bonds interact by means of the standard Coulomb 1/r interaction, as do charge sites on different molecules.

  2. Dynamics

    • The charge sites are given small masses and the entire system is propagated with an extened Lagrangian (see eq. (5)).

GROMACS implementation

The only implementation of a polarizable FF is the one reported in the first article on ForceBalance, but the model implemented is the QTPIE which differs from FQ.

First step implementation

I don't know if it is a good idea or not, but I am trying to set up a first (wrong) attempt just to see if I have everything I need. The code is something like this:

for step in range(num_steps):
    # 1. Evolve the system of one step
    sim.step(1)
    # 2. Compute the charges of each atom based on the new geometry
    # 3. Assign them to the corresponding atom

I am experiencing some problem with the third step, the update of the charges. Indeed, if I increase the charges at each step nothing seems to change in the energy. I attach an example of the code.

import simtk.openmm as omm
import simtk.openmm.app as app
import simtk.unit as unit
import openforcefield.topology as off_top
from openforcefield.typing.engines.smirnoff import ForceField
import numpy as np

# Read input geometry
pdbfile = app.PDBFile("water_box.pdb")
omm_topology = pdbfile.topology

# Water
solvent = off_top.Molecule.from_smiles('O')

# System
all_system = [solvent]

# Create the Open Force Field Topology from an OpenMM Topology object.
off_topology = off_top.Topology.from_openmm(openmm_topology=omm_topology, unique_molecules=all_system)

# Load the OpenFF "Parsley" force field.
forcefield = ForceField('openff-1.1.0.offxml')

# Parametrize the topology and create an OpenMM System.
system = forcefield.create_openmm_system(off_topology, charge_from_molecules=all_system)

# Propagate the System with Langevin dynamics.
time_step   = 2*unit.femtoseconds  # simulation timestep
temperature = 300*unit.kelvin      # simulation temperature
friction    = 1/unit.picosecond    # collision rate
num_steps   = 10000                # number of integration steps to run

# Output options.
trj_freq  = 1 # int(num_steps/1000)  # number of steps per written trajectory frame
data_freq = 1 # int(num_steps/100)   # number of steps per written simulation statistics

# System setup.
nonbonded_method  = app.PME
nonbonded_cutoff  = 8.0*unit.angstroms
constraints       = app.HBonds #None, app.HBonds, app.AllBonds

# Define the platform to use; CUDA, OpenCL, CPU, or Reference. Or do not specify
# the platform to use the default (fastest) platform
platform = omm.Platform.getPlatformByName('CPU')

# Parametrize the topology and create an OpenMM System.
system = forcefield.create_openmm_system(off_topology, charge_from_molecules=all_system)

# Create the integrator to do Langevin dynamics
integrator = omm.LangevinIntegrator(temperature, friction, time_step)

# Create the Simulation object
sim = app.Simulation(omm_topology, system, integrator, platform)

# System positions
positions = pdbfile.getPositions()

# Set the particle positions
sim.context.setPositions(positions)

# Randomize the velocities from a Boltzmann distribution at a given temperature.
sim.context.setVelocitiesToTemperature(temperature)

# Minimize the energy
print('Minimizing energy...')
sim.minimizeEnergy(maxIterations=500)
print("Energy minimization complete.")

# Set up the reporters to report energies and coordinates every 100 steps
pdb_reporter = omm.app.PDBReporter('trajectory.pdb', trj_freq)
state_data_reporter = omm.app.StateDataReporter('data.csv', data_freq, step=True,
                                                 potentialEnergy=True, temperature=True,
                                                 density=True)
sim.reporters.append(pdb_reporter)
sim.reporters.append(state_data_reporter)

for step in range(500):
    # 1. Evolve the system of one step
    sim.step(1)
    print("Energy of the system : ", sim.context.getState(getEnergy=True).getPotentialEnergy())
    # 2. Compute the charges of each atom based on the new geometry
    # 3. Assign them to the corresponding atom
    # Compute the forces and check their value
    forces = sim.context.getState(getForces=True).getForces()
    for force in system.getForces():
        # Select the nonbonded forces
        if isinstance(force, omm.openmm.NonbondedForce):
            # Run over all the particles of the system
            for j in range(system.getNumParticles()):
                # Get the parameters associated to each particle
                charge, sigma, epsilon = force.getParticleParameters(j)
                print(charge, sigma, epsilon)
                if j > 2:
                    new_charge = charge + 10.0e+10*unit.elementary_charge
                else:
                    new_charge = charge - 10.0e+10*unit.elementary_charge
                force.setParticleParameters(j, new_charge, sigma, epsilon)
                charge, sigma, epsilon = force.getParticleParameters(j)
                print(charge, sigma, epsilon)
    print("New energy of the system : ", sim.context.getState(getEnergy=True).getPotentialEnergy(), "\n")

This is the geometry:

HEADER
TITLE     Built with Packmol through Packmol Memgen
REMARK   Packmol generated pdb file, Packmol Memgen estimated parameters
REMARK   Home-Page: http://www.ime.unicamp.br/~martinez/packmol
REMARK
CRYST1    20.00    20.00    20.00  90.00  90.00  90.00 P 1           1
HETATM    1 H1   HOH A   1       7.387  12.021  11.223
HETATM    2 H2   HOH A   1       8.249  10.734  11.739
HETATM    3 O    HOH A   1       7.877  11.624  12.000
HETATM    4 H1   HOH A   2      16.284   9.965  12.754
HETATM    5 H2   HOH A   2      14.923   9.105  13.028
HETATM    6 O    HOH A   2      15.376   9.737  12.401
END
davidlmobley commented 4 years ago

@matteo-ambrosetti I personally don't have the expertise with OpenMM to help you attempt this; I'd guess that the OpenMM issue tracker is the best place to follow this up, since you're dealing with whether or not you can get it to work in OpenMM.

If the answer is "yes", then you should come back over here and we can talk about how to get it working in OpenFF.

I suspect that @jchodera 's Yank package makes updates to charges on the fly, or it did at one point, so it may be worth looking there (and at their related work on Perses) to see how they do this. However, this is not in the context of fluctuating charge models but rather free energy calculations.