6
$\begingroup$

I'm trying to create a CV for angle theta. To make this calculation, the right atom selections need to be made. The function is given below:

def theta_(traj):
    angI = mdtraj_topology.select('resid ' + residA + dicACGT['A'])
    angII = mdtraj_topology.select('resid ' + residA + ' and name P OP1 OP2 or resid ' + str(int(residA) - 1) + ' and name P OP1 OP2')
                              
    subIIIa = mdtraj_topology.select('resid ' + str(int(residA) - 1) + dicACGT[neighbors[0]])
    subIIIb = mdtraj_topology.select('resid ' + str(int(residA) + 1) + dicACGT[neighbors[1]])
    subIIIc = mdtraj_topology.select('resid ' + str(int(residT) - 1) + dicACGT[neighbors[2]])
    subIIId = mdtraj_topology.select('resid ' + str(int(residT) + 1) + dicACGT[neighbors[3]])
                             
    angIII = np.concatenate((subIIIa, subIIIb, subIIIc, subIIId))
    outputI = get_center_of_mass(traj, angI)
    outputII = get_center_of_mass(traj, angII)
    outputIII = get_center_of_mass(traj, angIII)
    return(angIII)

    thta=[]
    for f in range(traj.n_frames):
        vec1 = np.subtract(outputI[f], outputII[f])
        vec2 = np.subtract(outputIII[f], outputII[f])
        thta.append(angle_between_vectors(vec1,vec2))
    #return(thta)

coms = paths.MDTrajFunctionCV('theta_',
                              theta_,
                              topology = engine.topology)

The atom selection (for variables angI, andII, subIIIa, subIIIb, subIIIc, subIIId), however, varies when I check with MDTraj trajectory which gives the right selection.

Selection with MDTraj trajectory

n = trajs[0].to_mdtraj()
theta_(n)

array([ 43,  44,  45,  46,  47,  48,  49,  50,  51,  52,  53,  54,  55,
        56, 107, 108, 109, 110, 111, 112, 114, 116, 120, 121, 202, 203,
       204, 205, 207, 211, 212, 213, 264, 265, 266, 267, 268, 269, 270,
       271, 272, 273, 274, 275, 276, 277])

Selection from OPS wrapper:

coms(trajs[0])

array([43, 44, 45, 46, 47, 48, 49, 50])

The angle calculation as a consequence fails as well when creating the CV. I think I might be missing an argument to the wrapper. How do I make the right selection?

$\endgroup$
6
$\begingroup$

I can't tell why you're getting that specific difference, but most likely it's related to the fact that the CV is using global variables. It's possible that something changed in those globals. OPS CVs shouldn't use global variables, and in fact, you'll get an error if you try to save a CV that uses global variables.

CVs shouldn't use globals

The theta_ function you give uses the global variables mdtraj_topology, residA, residT, dicACGT, and neighbors. It also uses the function get_center_of_mass (and, if the early return is removed, angle_between_vectors ). There are all defined in global scope, and the OPS CVs don't save their global scope. If you somehow changed the values in the globals, that would explain the error. You would also get an error warning you about this if you tried to save your CV.

For more background, here's an article on scope in Python. The function wrapped by an OPS CV will not know about anything defined in global scope.

In general, the way to handle this is to pass these as additional variables to your theta_ function, and freeze them with keyword arguments when creating your CV:

def theta_(traj, mdtraj_topology, residA, residT, dicACGT, 
           neighbors, get_center_of_mass, angle_between_vectors):
    ...  # fill in your stuff here

coms = paths.MDTrajFunctionCV(
    'theta_',
    theta_,
    topology=engine.topology,
    mdtraj_topology=mdtraj_topology,
    residA=residA,
    residT=residT,
    dicACGT=dicACGT,
    neighbors=neighbors,
    get_center_of_mass=get_center_of_mass,
    angle_between_vectors=angle_between_vectors
)

For the specific case of atom index selection, however, it's probably better to do that outside the CV (although we'll use the same idea of passing additional parameters as keywords).

Recommended approach: Do the selection outside the CV

Right now, you're doing the selection inside the CV. In most cases, that's not what you would want. Keep in mind that a given CV instance is already tied to a specific topology, so the atom indices you select will always be the same. To handle selections from different topologies, you can create multiple CV instances that wrap the same function.

So I'd recommend doing the selection outside the function, before creating your CV. This avoids problems with the different globals affecting your selections, because you provide a constant selection to your CV.

You could define your theta_ function as:

def theta_(traj, angI, angII, angIII, get_center_of_mass, angle_between_vectors):
    ...  # start after the selection


# fill these in, or create functions to extract them
angI = ...
angII = ...
angIII = ...

This also allows you to validate the indices you selected, so you're sure you have the right selection before you start your simulation.

Next you'll give angI, angII, and angIII as parameters that you freeze when you create the CV. Note that (just as above) you'll also need to pass functions that you'll use (such as get_center_of_mass) as parameters.

coms = paths.MDTrajFunctionCV(
    'theta_',
    theta_,
    topology=engine.topology,
    angI=angI,
    angII=angII,
    angIII=angIII,
    get_center_of_mass=get_center_of_mass,
    angle_between_vectors=angle_between_vectors
)
$\endgroup$
1
  • $\begingroup$ It works now when I make the selection outside the function. However, I'm not able to return the theta values at the end with the wrapper but the function works when used with an MDTraj trajectory. I still think there's something missing with how I call the data with the wrapper. The error I get says: 'AttributeError: 'list' object has no attribute 'shape'' $\endgroup$
    – rika-27
    Aug 21 at 9:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.