distance restraints

[Thomas Evangelidis]

Greetings,

Is it possible to apply NMR-like (flat-welled) distance restraints in OpenMM without using the Meld plugin. I just want to restrain some waters to be within 15 Angstroms from the center of mass (COM) of the ligand. I found this thread where they discuss a similar case and suggest Meld plugin:

https://github.com/pandegroup/openmm/issues/924

In case flat-welled restraints are not support directly by OpenMM, maybe this example from the manual can work:

Code: Select all

force = CustomExternalForce('100*max(0, r-2)^2; r=sqrt(x*x+y*y+z*z)')
system.addForce(force)
for i in range(system.getNumParticles()):
force.addParticle(i, [])

Provided that the residue name of the ligand is "LIG", could someone please show me how to modify this so that "r" is the distance from the COM of the ligand?

Thank you in advance.
Thomas

[Thomas Evangelidis]

OK, I did some progress. I pass through the command line the coordinates of the COM of the ligand and set up distance restraints for the waters as follows:

Code: Select all

# Set up the distance restraints for the waters
if args.ligand_COM:
    force = CustomExternalForce('1000*max(0, r-'+str(0.1*args.WATER_CUTOFF)+')^2; r=sqrt( (x-x0)*(x-x0) + (y-y0)*(y-y0) + (z-z0)*(z-z0) )')   # distances are in nm!
    water_mask = AmberMask(parm, ':WAT & @O')
    force.addPerParticleParameter('x0')
    force.addPerParticleParameter('y0')
    force.addPerParticleParameter('z0')
    COM_coords = args.ligand_COM.split()
    for idx in water_mask.Selected():
        force.addParticle(idx, [float(COM_coords[0])*u.angstroms, float(COM_coords[1])*u.angstroms, float(COM_coords[2])*u.angstroms])
    system.addForce(force)

where args.WATER_CUTOFF=15 and args.ligand_COM="-17.4308 2.3502 15.8220" for example. Is there any way to ask OpenMM to calculate the COM coordinates and update x0, y0, z0 in every frame?

[Peter Eastman]

The best way to do this is with CustomCentroidBondForce. You can use it to apply restraints between the centers of mass of groups of atoms. See http://docs.openmm.org/7.1.0/api-python ... dBondForce.

[Thomas Evangelidis]

Can I apply this conditional restraint with CustomCentroidBondForce, namely if the distance is greater than the cutoff apply a force, otherwise not?

[Peter Eastman]

Yes. You do it in exactly the same way: by using a max() function so the energy remains 0 until r gets large enough.