How to add a radius of gyration constraint?

[Wei Lu]

I want apply the constraint to a subset of atoms.
The energy E is defined as below.

Should I use CustomCVForce? I didn't find any example using this force, so I'm not sure how to write it.
We can also write the rg^2 as following

Then, use CustomBondForce, but I don't know how to include Rg0 into the energy computation.

Any help would be greatly appreciated. Thank you!

[Wei Lu]

We are able to solve this ourself. The code is copied below in case other people are interested.

Code: Select all

def rg_bias_term(oa, k_rg=4.184, rg0=0):
    nres, ca = oa.nres, oa.ca
    rg_square = CustomBondForce("1/normalization*r^2")
    rg_square.addGlobalParameter("normalization", nres*nres)
    for i in range(nres):
        for j in range(i+1, nres):
            rg_square.addBond(ca[i], ca[j], [])
    rg = CustomCVForce(f"{k_rg}*(rg_square^0.5-{rg0})^2")
    rg.addCollectiveVariable("rg_square", rg_square)
    rg.setForceGroup(27)
    return rg