Minimum distance to a group of atoms as a collective variable

[lewis martin]

Hi all,
Is it possible to define a custom force as the minimum distance between one atom and a group of atoms, or alternatively a molecule COM and a group of atoms? The purpose is to use minimum distance to lipid headgroups as a collective variable for ligand diffusion across the bilayer, for example.

The starting point I came up with is a CustomCompoundBondForce, made up of many bonds between a ligand and an atom from each of the lipid headgroups. Then the energy expression might be

Code: Select all

 "min(distance(p1, p2), distance(p1,p3), .... ) "

and so on, where p1 is the ligand atom, and p's 2 through 50 are lipids.. Is that the best way?

Thanks,
Lewis

[Peter Eastman]

I don't really have a better solution. This will be very inefficient for large numbers of groups, but perhaps with 50 it won't be too terrible. If you want to use the molecule COM rather than a single atom you could do the same thing with CustomCentroidBondForce.

[lewis martin]

Thanks Peter.
I think to get a reasonably sized membrane I would need around 64 lipids per leaflet, so one restraint per headgroup, and then two leaflets, is 128 bonds. Surely that would be less extra distance calculation than adding even a single extra lipid.

Anyways Ill drop a gist here later with some example code if I figure it out.

[lewis martin]

OK I tested it out by adding forces that calculate minimum distance of a sodium ion to both the upper and lower leaflets of a CHARMMGUI lipid bilayer.

If anyone's interested, see a notebook here:
https://github.com/ljmartin/lipid_diffu ... mple.ipynb
and other files are here:
https://github.com/ljmartin/lipid_diffu ... ar_example

The addition of the two forces (one for each bilayer leaflet) reduces simulation speed to 80% of the same system without the forces, so I guess about 10% per force. More than I thought but certainly not a deal breaker.

Anyways, the goal was to reproduce the free energy surface seen in https://pubs.acs.org/doi/10.1021/acs.jctc.8b00882 , and this code + metadynamics should be enough to do that.
cheers
lewis