Coordinates to trajectory

[kmb413]

I remember we were talking about taking an array of coordinates and using the pt.iterload() function to load them, given the parmfile. But I forgot how. Would we have to write the coords to a .rst7 file and then load?

Thanks!

[hainm]

Would we have to write the coords to a .rst7 file and then load?

no.

you can combine any file formats you want (as long as pytraj/cpptraj supports http://amber-md.github.io/pytraj/latest/read_and_write.html#supported-file-formats)

fnames = ['myname_0.nc', 'myname_1.rst7', 'myname_2.pdb']
parmfile = 'my.parm7'
traj = pt.iterload(fnames, parmfile)
pt.rmsd(traj, ref=..., mask=...)

[hainm]

Note that you can use either load or iterload method. load will load all data to memory (good for small trajectory) but iterload does not (good for large trajectory).

[kmb413]

Right, but we're trying to go from pose to trajectory, of which I have access to only the coordinates. If we can avoid dumping the pose to a file (like a .pdb), that would be ideal.

[hainm]

example here

import pytraj as pt

topology = pt.load_topology('your.parm7')
# xyz is 3D coordinates, shape=(n_snapshots, n_atoms, 3)
traj = pt.Trajectory(xyz=xyz, top=topology)
pt.rmsd(traj, ref=..., mask=...)

You need to find a way to get 3D coordinates from rosetta files.

[hainm]

This is example showing how to connect pytraj to 3rd party programs

http://amber-md.github.io/pytraj/latest/tutorials/interface_with_other_packages.html

Should be doable for pyrosetta too.

[hainm]

note: it's better to convert any coordinates to numpy 3D array.

[hainm]

@kmb413 any progress?

[kmb413]

Getting there. Found an easy way to get all of the coordinates with one function, now I just need to figure out how the coordinates and parmfile are related. I assume the order in which the coordinates appear in the array are important for the parmfile. I don't think Rosetta and Amber list them in the same order, but I still have to check.

[hainm]

I assume the order in which the coordinates appear in the array are important for the parmfile.

Yes, it's important.

I don't think Rosetta and Amber list them in the same order, but I still have to check.

let me know if you need help or have any question. But the atom order in parmfile should match to the original pdb file.

[hainm]

@kmb413 if we can internally convert coords from rosetta to pytraj, we can use Jason's script to minimize, then evaluate energy

https://github.com/ParmEd/ParmEd/issues/243

(adapted from his code)

import pytraj as pt
import sander
from scipy.optimize import minimize
import numpy as np

def target(coords):
    """ Target function to minimize -- just the energy """
    sander.set_positions(coords)
    e, f = sander.energy_forces()
    return e.tot, -np.array(f) # forces are -gradient

def minimize_trajectory(traj, prmtop, igb=8, **min_options):
    '''

    Parameters
    ----------
    traj : pytraj.Trajectory (mutable)
    '''
    inp = sander.gas_input(igb)
    for index, frame in enumerate(traj):
        with sander.setup(prmtop, frame.xyz, frame.box, inp):
            res = minimize(target, frame.xyz.flatten(), **min_options)

        # update frame.xyz will update coordinates for traj too.
        frame.xyz[:] = res.x.reshape(traj.n_atoms, 3)
    return traj

def main():
    traj0 = pt.datafiles.load_tz2()
    parmfile = traj0.top.filename
    traj = pt.Trajectory(xyz=traj0[:2].xyz, top=parmfile)

    energy_dict = pt.energy_decomposition(traj, prmtop=parmfile, 
                              igb=8, dtype='dataframe')
    print(energy_dict[['tot', 'gb']])

    progress = False
    options = dict(maxiter=100,
                   disp=progress)
    min_options = dict(method='L-BFGS-B', jac=True, tol=1e-4)
    min_options['options'] = options

    traj = minimize_trajectory(traj, prmtop=parmfile, igb=8, **min_options)

    # if the minimization is time consuming, 
    # we can just save all snapshots to a single file 
    # (rather getting thousands of .rst7, .out, ...)
    # traj.save('my_traj.nc', overwrite=True)
    energy_dict = pt.energy_decomposition(traj, prmtop=parmfile, 
                         igb=8, dtype='dataframe')
    print(energy_dict[['tot', 'gb']])

if __name__ == '__main__':
    main()

updated: new code

(note: note validate the script yet)

[kmb413]

Cool! I think I'm close.

On Tue, Feb 23, 2016 at 6:22 PM, Hai Nguyen notifications@github.com wrote:

@kmb413 https://github.com/kmb413 if we can internally convert coords from rosetta to pytraj, we can use Jason's script to minimize, then evaluate energy

ParmEd/ParmEd#243 https://github.com/ParmEd/ParmEd/issues/243

(adapted from his a code)

import pytraj as ptimport sanderfrom scipy.optimize import minimize def target(coords): """ Target function to minimize -- just the energy """ sander.set_positions(coords) e, f = sander.energy_forces() return e.tot, -np.array(f) # forces are -gradient

parmfile = 'something.parm7' traj = pt.Trajectory(xyz=rosetta_coords, top=parmfile)

inp = sander.gas_input(8)for index, frame in enumerate(traj): with sander.setup(parmfile, frame.xyz, frame.box, inp): res = minimize(target, frame.xyz.flatten(), method='L-BFGS-B', jac=True, tol=1e-8, options=dict(maxiter=10000, disp=True))

traj.xyz[index] = res.x.reshape(traj.n_atoms, 3)

energy_dict = pt.energy_decomposition(traj, prmtop=parmfile, igb=8)

(note: note validate the script yet)

— Reply to this email directly or view it on GitHub https://github.com/kmb413/DecoyDiscrimination/issues/11#issuecomment-187962082 .

[kmb413]

Well, looks like since we want to do a pure python Rosetta/Amber interface, we're gonna use ParmEd to make the parmfile from the OpenMM structure (created from the pose), like in Jason's example.

So I've figured out how to make a topology from a parmed object, but how do I use this to make a trajectory? I've tried digging around the doc and couldn't find an example, and I've tried pt.iterload( pose-coordinates, topology, igb=8 ) and it gives

File "/scratch/kmb413/amber_git/amber/lib/python2.7/site-packages/pytraj/io.py", line 221, in iterload
     return load_traj(*args, **kwd)
File "/scratch/kmb413/amber_git/amber/lib/python2.7/site-packages/pytraj/io.py", line 250, in load_traj
     ts._load(filename)
File "/scratch/kmb413/amber_git/amber/lib/python2.7/site-packages/pytraj/trajectory_iterator.py", line 227, in _load
     raise ValueError("filename must a string or a list of string")

Can you point me to an example, maybe, of pytraj+parmed stuff?

[hainm]

Please try:

fn = 'my.parm7'

save parmed object to parm7

parm.save(fn)

load to pytraj's topology

top = pt.load_topology(fn)

construct pytraj Trajectory

traj = pt.Trajectory(xyz=some_coords, top=top)

compute energy

data = pt.energy_decomposition(traj, prmtop=fn, igb=8)

But please make sure that you can reproduce energy from regular workflow (use tleap)

Hai

On Mar 1, 2016, at 3:32 PM, Kristin Blacklock notifications@github.com wrote:

Well, looks like since we want to do a pure python Rosetta/Amber interface, we're gonna use ParmEd to make the parmfile from the OpenMM structure (created from the pose), like in Jason's example.

So I've figured out how to make a topology from a parmed object, but how do I use this to make a trajectory? I've tried digging around the doc and couldn't find an example, and I've tried pt.iterload( pose-coordinates, topology, igb=8 ) and it gives

File "/scratch/kmb413/amber_git/amber/lib/python2.7/site-packages/pytraj/io.py", line 221, in iterload return load_traj(_args, *_kwd) File "/scratch/kmb413/amber_git/amber/lib/python2.7/site-packages/pytraj/io.py", line 250, in load_traj ts._load(filename) File "/scratch/kmb413/amber_git/amber/lib/python2.7/site-packages/pytraj/trajectory_iterator.py", line 227, in _load raise ValueError("filename must a string or a list of string") Can you point me to an example, maybe, of pytraj+parmed stuff?

— Reply to this email directly or view it on GitHub.

[kmb413]

Oh, interesting, so we still need to save it as a parm7. That didn't quite work and I think it has something to do with the structure object for parmed not being parameterized yet <Structure 123 atoms; 12 residues; 122 bonds; NOT parametrized> so I will find out how to do that and then try this. Thanks! Have a nice conference.

[hainm]

No, that does not work because you still need to do further step like speciyfing force field, radii, ...

[hainm]

We need to save to parm7 because sander needs it.