Unwrapping a Trajectory in loos

[krb0073]

HI all,
I am trying to make a tool to calculate the Mean squared displacement of a pure bilayer trajectory, but I cant find a Way to unwrap the trajectory. What would you recommend doing to unwrap the trajectory via loos ?

Kyle

[agrossfield]

Hi Kyle,

There isn't a tool that does it. I've thought about how to do it, though, and I'd be happy to walk you through how to do it. Assuming you just want the unwrapped centers of mass or centroids, what I'd do is something like

lipids = an array of AtomicGroups made by selecting the lipids and doing splitByMolecule()
centers = an array of GCoords that is dimensioned [num_lipids]

first = True
for frame in traj:
    for index in range(len(lipids)):
        centers[index] = lipids[index].centroid()
    
    if not first:
        first = False
        # if a lipid has moved more than 1/2 of the box, it must have jumped
        box = frame.periodicBox()
        half_box = 0.5 * box
        for index in range(len(lipids)):
            diff = prev_centers[index] - centers[index]
            fix = loos.GCoord()
            if diff.x() > half_box: fix[0] -= box.x()
            elif diff.x() < -half_box: fix[0] +=box.x()
            # repeat for y and z

        print out centers
        prev_centers = centers

This is a little simplistic, but for a sane MD trajectory of a membrane it ought to work fine.

[agrossfield]

Kyle -- if you get something working, let me know. I'd be happy to include it in a future version of loos. On the flip side, if you run into trouble, ping me here or by email and I'll help out.

Alan

[krb0073]

Thanks Alan,
Here is what i was able to get done in regards to writing out the unwrapped postions of the lipid centers.

import sys
import loos
import loos.pyloos
from math import floor
# model file
lipid_sel = sys.argv[1]
# this is the gneric slection for the lipid
traj = sys.argv[3]
mdl = sys.argv[2]
#m kae the loos system
system = loos.createSystem(mdl)
# load the trajecotry
traj = loos.pyloos.Trajectory(traj,system)
# select the atoms for the lipids
## spilt into the invidaul molecules of lipid
lipids = loos.selectAtoms(system,lipid_sel).splitByMolecule()
# boolean to check if this is frame is frame 0
frist = True
# store the unwrapped
unwraped = []
# loop inro the frame
## this is the order of the lipid resname-resid-segname
### this will be usefulll later
print("# resname-resid-segid")
print('#' ,*[f"{lipids[i][0].resname()}-{lipids[i][0].resid()}-{lipids[i][0].segid()} " for i in range(len(lipids))])
print(" #frame Lipid_1_center_x Lipid_1_center_y Lipid_1_center_z ...")
for frame in traj:
    # loop into the lipids
    centers = []
    for index in range(len(lipids)):
        centers.append(lipids[index].centroid()) # gemotric center of the lipid
        # postion 0 for frame 0 is always the starting postion
        if frist == True:
            unwraped = centers
    # in frame one nothing should have moved
    if frist != True:
        box = frame.periodicBox() # get a G cord of the PBC
        half_box = 0.5 * box  # a square box's cneter is at the halfb way between all the postions
        for index in range(len(lipids)):
            diff = centers[index] -  prev_centers[index]
            fix = loos.GCoord() # blank GCorrd() == (0,0,0)
            # assign the x,y,z cord of fix is these condtion are met
            ## x, y, z may not scale perct to be a exact cudde
            if diff.x() > half_box.x():
                fix[0] -= box.x()
            elif diff.x() < -box.x():
                fix[0] += box.x()

            if diff.y() > half_box.y():
                fix[1] -= box.y()
            elif diff.y() < -half_box.y():
                fix[1] += box.y()

            if diff.z() > half_box.z():
                fix[2] -= box.z()
            elif diff.z() < - half_box.z():
                fix[2] += box.z()
            # with the fixed vector made , to figure out where
            ## lipid is in relation to the box we can adjust the postion of
            ### the center
            centers[index] += fix
    all = []
    for i in centers:
        for j in i:
            all.append(j)
    print(traj.index(),*all)
    prev_centers = centers
    frist = False

I did a few of the lipids in the system by hand indpendantly to see if I got the same answer as the output from the code. Seemed to be okay. I am unsure if it would be better to write out a tracjotry of the lipids center as they progess? Let me know if you have any commnets, and thanks for all of the help.
Kyle

[agrossfield]

This looks pretty good. I see some typos in the comments and unfortunate variable names that should be fixed (e.g. frist). I noticed a couple of things that look odd, though -- you set up the variable unwraped (please rename to unwrapped), but then don't actually use it -- you use prev_centers like I did in my example.

On the question of writing out a trajectory vs just writing the centroids: I see arguments both ways. If you want to write a trajectory, you can see an example in simple_traj_transform.py (this writes XTC files, but replace XTCWriter with DCDWriter and you'll get a dcd). However, you will get one piece of weirdness -- you'll inherit the periodic box size from the underlying AtomicGroup, but the coordinates will be way outside the box. I don't know if there would be a problem with doing that, but you might be breaking some underlying assumptions in the way xtc does its compression. You could add it as an option using a command line flag. The nice thing about writing it as a text file is you can test using something like gnuplot

plot "foo.dat" u 2:3 w l, "foo.dat" u 5:6" w l

would plot the trajectories of the first two lipids, and allow you to look for discontinuities or other weird behavior.

If you want to make a pull request and add it to loos, there will be some other things you'll need to do:

1. switch it over to using program options. See one of the other programs, eg all_stacking.py, for what that should look like
2. while you're at it, add a comment at the top indicating what the program does, a license block, and a detailed fullhelp. Again, all_stacking is good template to work from
3. The first line of the output should be the command line arguments, prepended by a "#"

I'm sure there's some other stuff, but this is what occurs to mind.

[krb0073]

Alan,
I finally got the ode to work the way i wanted the full thing is here.

# @Author: Kyle Billings <kbillings>
# @Date:   2021-09-22T18:00:20-04:00
# @Email:  [email protected]
# @Filename: lipid_un.py
# @Last modified by:   kbillings
# @Last modified time: 2021-09-25T16:20:39-04:00


#!/usr/bin/env python3

"""
Track a base stacking through a trajectory.  Intended for use with
nucleic acids, to identify base stacking.
"""

"""
  This file is part of LOOS.
  LOOS (Lightweight Object-Oriented Structure library)
  Copyright (c) 2021 Alan Grossfield, Grossfield Lab
  Department of Biochemistry and Biophysics
  School of Medicine & Dentistry, University of Rochester
  This package (LOOS) is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation under version 3 of the License.
  This package is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
  You should have received a copy of the GNU General Public License
  along with this program.  If not, see <http://www.gnu.org/licenses/>.
"""
import sys
import loos
import loos.pyloos
import numpy as np
import argparse
# making the fullhelp for the program
def fullhelp():
    print("""lipid_unwrap.py takes the coordinates from a trajectory file to unwrap
    the periodic box(PB) for a selection .At the moment the only rectangular
    box based periodic box conditions are able to be used. The subset of atoms
    that you use could be any valid loos selection ,but the program is intended
    to work on lipids,tracking the geometric center of a molecule at the
    current frame and comparing that position to the previous one.
    If the difference between the two positions vary by more than half of the
    P.B  length (in x,y, or z) a correction factor is added to the current position
    to correct the positions it moves outside the boundaries.

     If desired, a pdb and DCD of the unwrapped trjectory can be. This is done by
     taking every molecule if the selection to (0,0,0) then translating it to the
     unwrapped. For making the trajecotry,the assumtion is made that the
     Lx,Ly or Lz of the P.B will not be greater than 1,000 Angstroums.

    Usage: lipid_unwrap.py [OPTIONS] selection_string system_file trajectory_file
        --  selection_string identifies the atoms to be analyzed. They will be
            split by molecule number calculate position for every lipid adjusted
            without the wrapping
        --  system_file is a PDB, PSF, prmtop, etc
        --  trajectory_file eg DCD, xtc. Its contents must precisely
            match the system file

    Options:
        --skip: # of residues to exclude from the front of each trajectory
        --stride: how to step through the trajectory (eg --stride 10 will read
              every 10th frame)
        --fullhelp: prints this message
        --output_traj: output a dcd of the unwrapped trajectory defult(--output_traj 0)
        --output_prefix: prefix to use for the output defult is output

    Exmaple

""")
def findFix():
    """ find Fix is the takes the center befor it is unwrapped then checks the lipids
    until to see if the new postion has moved more than half the box then the box
    lenght is added to the postion of that atom to negate the wrapping. If the
     corrdinate of fixed atom still is lager than half the box a factor is multiped to
      the box lenght added to the system untill the distacne is within reason"""
    fix = loos.GCoord() # blank GCorrd() == (0,0,0)
    a ,b ,c, d,e,f = 1, 1,1,1,1,1
    check = True
    while check:
        diff = centers[index] -  prev_centers[index]
        tracker = 0
        if diff.x() >  half_box.x():
            fix[0] -= box.x() * a
        elif diff.x() < -half_box.x():
            fix[0] += box.x() * b
            tracker =1
        if diff.y() > half_box.y():
            fix[1] -= box.y() * c
            tracker = 2
        elif diff.y() < -half_box.y():
            fix[1] += box.y() * d
            tracker = 3

        if diff.z() > half_box.z():
            fix[2] -=  box.z() * e
            tracker = 4
        elif diff.z() < -half_box.z():
            fix[2] += box.z() * f
            tracker = 5
        ans = centers[index] + fix
        if prev_centers[index].distance(ans) > half_box.x():
            print('yup')
            C = [x for x in centers[index]]
            P = [x for x in prev_centers[index]]
            # now find the exact thing that messed up
            print(ans)
            if tracker == 0:
                if ans.x() - prev_centers[index].x() < half_box.x():
                    a += 1
            elif tracker == 1:
                if ans.x() -prev_centers[index].x() > half_box.x():
                    b += 1
            elif tracker == 2:
                if ans.y() - prev_centers[index].y() < half_box.y():
                    c += 1
            elif tracker == 3:
                if ans.y() - prev_center[index].y() > half_box.y():
                    d += 1
            elif tracker == 4:
                if ans.z() - prev_centers[index].z() < half_box.z():
                    c += 1
            elif tracker == 5:
                if ans.z() - prev_center[index].y() > half_box.z():
                    d += 1
        else:
            print("nope")
            check = False
    return fix
class FullHelp(argparse.Action):
    def __init__(self, option_strings, dest, nargs=None, **kwargs):
        kwargs['nargs'] = 0
        super(FullHelp, self).__init__(option_strings, dest, **kwargs)
    def __call__(self, parser, namespace, values, option_string=None):
        fullhelp()
        parser.print_help()
        setattr(namespace,self.dest,True)
        parser.exit()

parser = argparse.ArgumentParser(description="Unwrap PBC lipids")
parser.add_argument('system_file', help="File describing the system")
parser.add_argument('selection_string',help="Selection string describing which residues to use")
parser.add_argument('traj_file',help="File contraing the trajecotry")
parser.add_argument('--fullhelp',help="Print detailed description of all options",action=FullHelp)
parser.add_argument('--skip',type=int,default=0,help='# of frame to skip')
parser.add_argument('--stride' , type=int,default=1,help="Read every nth frame")
parser.add_argument('--output_traj',type=int,default=0,help='produce an unwrapped trajectory')
parser.add_argument('--output_prefix',default='unwrapped_output',type=str,help='name of the tractory file to write (DCD format)')
args = parser.parse_args()
print(args)
pre = args.output_prefix
header = " ".join([f"{i}" for i in sys.argv])
system = loos.createSystem(args.system_file)
# load the trajectory
traj = loos.pyloos.Trajectory(args.traj_file,system,stride=args.stride , skip=args.skip)
# select the atoms for the lipids
## spilt into the invidaul molecules of lipid
lipids = loos.selectAtoms(system,args.selection_string).splitByMolecule()
# boolean to check if this is frame is frame 0
first = True
# loop into the frame
## this is the order of the lipid resname-resid-segname
### this will be usefull later
outtraj = loos.DCDWriter(pre +".dcd")
header += '\nresname-resid-segid'
big_list = " ".join([f"{lipids[i][0].resname()}-{lipids[i][0].resid()}-{lipids[i][0].segid()} " for i in range(len(lipids))])
header += f"\n{big_list}"
header += "\nframe Lipid_1_center_x Lipid_1_center_y Lipid_1_center_z ..."
for frame in traj:
    # loop into the lipids
    centers = []
    for index in range(len(lipids)):
        centers.append(lipids[index].centroid()) # geometric center of the lipid
        # position 0 for frame 0 is always the starting position
    # in frame one nothing should have moved
    if not first:
        box = frame.periodicBox() # get a G cord of the PBC
        half_box = 0.5 * box  # a square box's cneter is at the halfb way between all the postions
        for index in range(len(lipids)):
            fix =findFix()
            centers[index]+= fix
            # fix the werid jump issssues if the atoms move more than wanted
            # if the user has the flag on for the trajectory then this is True
    if args.output_traj:
        # loop through the list of center locations and
        moved_lipids = loos.AtomicGroup()
        for lipid , center in zip(lipids,centers):
            lipid.centerAtOrigin() # move tha lipid to the center
            # translate to the center
            lipid.translate(center)
            moved_lipids.append(lipid)
        # setting the pbc to be large
        moved_lipids.periodicBox(loos.GCoord(1000,1000,1000))
        if len(lipids) >= 5:
            moved_lipids.centerAtOrigin()
        outtraj.writeFrame(moved_lipids)
        if first == True:
            temp_lipids = loos.selectAtoms(system,args.selection_string)
            frame_copy = temp_lipids.copy()
            temp_lipids.periodicBox(loos.GCoord(1000,1000,1000))
            frame_copy.pruneBonds()
            frame_copy.renumber()
            pdb = loos.PDB.fromAtomicGroup(frame_copy)
            pdb.remarks().add(header)
            with open(f"{pre}.pdb","w") as out:
                out.write(str(pdb))
    all = [traj.index()]
    for i in centers:
        for j in i:
            all.append(j)
    all = np.array(all)
    if first:
        all_centers = all
    else:
        all_centers = np.row_stack((all_centers,all))
    prev_centers = centers
    first = False
np.savetxt(f'{pre}.txt',all_centers,header=header)

Thanks again for all of your help.
Kyle

[agrossfield]

I'll go through it more carefully sometime in the next few days. If you'd like this to become part of loos, you have 2 options:

1. Use GitHub to make a fork of loos, add your new program to it, and create a pull request.
2. Tell me you want me to do it, and I'll merge it in.

The former has the advantage of making your contribution public and trackable, and will get done when you want it to. The latter requires less work from you, since I'll do the remaining tweaking.

Let me know which you prefer.