openmm / openmm

OpenMM is a toolkit for molecular simulation using high performance GPU code.
1.46k stars 511 forks source link

Implementation of Configurational-Bias Monte Carlo #4450

Closed ruyeli closed 3 weeks ago

ruyeli commented 6 months ago

Hi,

I am currently working on implementing a Configurational-Bias Monte Carlo(CBMC) simulation using OpenMM with empirical force field like opls, where my goal is to insert additional water molecules into an existing water box, which initially contains 256 water molecules.

The CBMC approach I'm employing involves generating multiple candidate positions for the insertion of an additional water molecule (e.g., positions A, B, C, D, E, as illustrated in the figure below). For each candidate position, I need to calculate the energy contribution of the newly added molecule to determine its Boltzmann weight and decide the most favorable insertion position.

water2

Currently, my workflow involves calculating the total energy of the original 256-water box and then recalculating the total energy for the entire system (257 waters) for each candidate insertion. The energy difference gives the contribution of the newly added molecule. However, this method requires 20 full system energy calculations for 20 candidate positions, leading to significant computational inefficiency.

I am seeking advice or suggestions on a more efficient approach within OpenMM that would allow me to directly calculate the energy contribution of just the newly added water molecule, without the need to recalculate the energy of the entire system for each candidate position. Such a method would greatly enhance the efficiency of the CBMC simulations I am conducting.

Thank you in advance for your insights and assistance.

swails commented 6 months ago

A few observations here.

An OpenMM System object has a fixed number of particles, so you would need to use a separate System (and corresponding Context) object for each particle configuration you want to be accessible to your simulation. Or you'd need to repeatedly update particle parameters for the "new" molecules to eliminate their nonbonded potential terms (which is an expensive operation, particularly on GPUs). This is a consequence of how OpenMM (and many other MD packages) are written -- your use case is not well supported in the way OpenMM models molecular systems.

I suspect it's possible to implement your application with the existing interface (you may have to make some approximations, or introduce some small extra terms that hopefully yield minimal error).

I am seeking advice or suggestions on a more efficient approach within OpenMM that would allow me to directly calculate the energy contribution of just the newly added water molecule, without the need to recalculate the energy of the entire system for each candidate position.

This is more challenging at a more conceptual level. The premise of the statement "energy contribution of just the newly added water molecule" is that the potential energy function is (largely) pairwise-decomposable. What I mean by this is that, given 3 water molecules (call them A, B, and C), the potential energy is an additive function of interactions between pairs of particles. Said specifically, if you calculate the potential energy of each 3 pairs (A-B, A-C, B-C), the total energy is simply the sum. If this is the case, you can find the energy of any subsystem (A-B, B-C, or A-C) if you simply store the 3 pairwise interactions without having to recompute them.

However, most of the commonly used potential energy functions are not actually pairwise decomposable (that is, adding water molecule C will actually change the interaction between A and B). While the traditional classical force fields define a pairwise nonbonded potential, the way long-range electrostatics are handled by Ewald-based methods (like PME) are not pairwise decomposable, and continuum-based implicit solvation models are also not pairwise decomposable because the dielectric boundaries are nonlocal functions over all atoms in the system.

You might be able to make an approximate model by leveraging CustomNonbondedForce with defined interaction groups defined between "added" waters and the rest of the system, either picking a long-range electrostatic model that is pairwise decomposable (like reaction field), or assuming that the contribution of a single water molecule on the reciprocal space calculation in PME is likely to be small, and it's acceptable to treat the "added" particle entirely in direct space (an assumption whose errors would likely need to be measured).


What you're trying to do sounds pretty challenging from an implementation perspective, and would probably require a reasonably deep understanding of how OpenMM works at a low (implementation) level. You'll probably need to make use of the following features:

Good luck!