From 90a48fac1978eaf9988d4502503a074d828b3356 Mon Sep 17 00:00:00 2001 From: krb0073 <41389987+krb0073@users.noreply.github.com> Date: Mon, 27 Sep 2021 18:19:55 -0400 Subject: [PATCH 1/6] Added lipid_unwrap.py to PYLOOS --- Packages/PyLOOS/lipid_unwrap.py | 235 ++++++++++++++++++++++++++++++++ 1 file changed, 235 insertions(+) create mode 100644 Packages/PyLOOS/lipid_unwrap.py diff --git a/Packages/PyLOOS/lipid_unwrap.py b/Packages/PyLOOS/lipid_unwrap.py new file mode 100644 index 000000000..14fac4786 --- /dev/null +++ b/Packages/PyLOOS/lipid_unwrap.py @@ -0,0 +1,235 @@ +# @Author: Kyle Billings +# @Date: 2021-09-22T18:00:20-04:00 +# @Email: krb0073@mix.wvu.edu +# @Filename: lipid_un.py +# @Last modified by: kbillings +# @Last modified time: 2021-09-26T15:07:15-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 . +""" +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. + + The output txt file contrain the radial position of the atom comparing it to the + center of mass to the bilayer. + + 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: + python /media/bak12/Analysis/unwrap/lipid_unwrap.py --ouput_traj 1 --output_prefix='unwrapped_foo' + foo.psf 'resname == "POPC"' foo.dcd + + We will obatian a dcd of the unwrapped trajectory with the name of unwrapped_foo.dcd + a pdb of unwrapped_foo.pdb to be used to veiw the dcd, and unwrapped_foo.txt that + has the radial distacne of the ceneter of the lipid to the COM of bilayer + at that frame for the POPC lipids in the membrane""") +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 # scaling factors of each + check = True # checks if the condtion is true + while check: # while the ditstance propational to a jump + diff = centers[index] - prev_centers[index] # final minus statrs + tracker = 0 # tracks which if it enterhs + # looking for points that jumpedd + 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 + # the new postion is the wrapped postion plus the fix + ans = centers[index] + fix + # if the distacne from the old location to the new is large than the + ## half_box + if prev_centers[index].distance(ans) > half_box.x(): + # using the value of tracker we find which if was used to fix the + ## bad corrdinate + if tracker == 0: + # if the ans - previous is less than half box the jump + ## was from the postive side of the pbc to the negative + if ans.x() - prev_centers[index].x() < half_box.x(): + a += 1 + elif tracker == 1: + # is the diffrecne is grater than half_box the jump was from + ## negative to postive + 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: + # if there was no large jump we exit the loop + 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() +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 +if args.output_traj: # fixes the issues where an emptry traj is made when +## not wanting one out + 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_R Lipid_2_center_R Lipid_3_center_R ..." +for frame in traj: + # loop into the lipids + COM_R = loos.selectAtoms(system,args.selection_string).centerOfMass() + 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 + R_coord = [traj.index()] + 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 + R_coord.append(np.sqrt((centers[index].x()-COM_R.x())**2 +(centers[index].y()-COM_R.y())**2+(centers[index].z()-COM_R.z())**2)) + # 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()] + if first: + for index in range(len(lipids)): + R_coord.append(np.sqrt((centers[index].x()-COM_R.x())**2 +(centers[index].y()-COM_R.y())**2+(centers[index].z()-COM_R.z())**2)) + all_centers = np.array(R_coord) + else: + R_coord = np.array((R_coord)) + all_centers = np.row_stack((all_centers,R_coord)) + prev_centers = centers + first = False +np.savetxt(f'{pre}.txt',all_centers,header=header) From 48e9a46aa355cb6f7c7c57f4cbcb4f5a9d2eeae1 Mon Sep 17 00:00:00 2001 From: krb0073 <41389987+krb0073@users.noreply.github.com> Date: Sun, 3 Oct 2021 15:25:56 -0400 Subject: [PATCH 2/6] Update lipid_unwrap.py Added a top commnet for what the file does. Changed the a,b...f to one dictonary that was all the values of the tracker to the fix. added if __name__ in main . Spelling fixed(not well). Commneted on what i am doin in the trajectory making stuff. removed so silly thing about making the pdb . incresed the number of the new PBC Box --- Packages/PyLOOS/lipid_unwrap.py | 249 +++++++++++++++++--------------- 1 file changed, 132 insertions(+), 117 deletions(-) diff --git a/Packages/PyLOOS/lipid_unwrap.py b/Packages/PyLOOS/lipid_unwrap.py index 14fac4786..3e2798ce1 100644 --- a/Packages/PyLOOS/lipid_unwrap.py +++ b/Packages/PyLOOS/lipid_unwrap.py @@ -1,16 +1,9 @@ -# @Author: Kyle Billings -# @Date: 2021-09-22T18:00:20-04:00 -# @Email: krb0073@mix.wvu.edu -# @Filename: lipid_un.py -# @Last modified by: kbillings -# @Last modified time: 2021-09-26T15:07:15-04:00 - - #!/usr/bin/env python3 """ -Track a base stacking through a trajectory. Intended for use with -nucleic acids, to identify base stacking. +Unwraps a systems perdoic box of for each frame of a +given selection. This creates an output for calcuting the +means sqaured diffusion. """ """ @@ -29,29 +22,34 @@ You should have received a copy of the GNU General Public License along with this program. If not, see . """ +# making the fullhelp for the program + + + + 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 + print( + """lipid_unwrap.py takes the coordinates from a trajectory file to unwrap + the periodic box 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 + 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 + periodic box 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 + If desired, a pdb and DCD of the unwrapped trajectory 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. + unwrapped. For making the trajectory,the assumption is made that the + Lx,Ly or Lz of the periodic box will not be greater than 1000000 Angstroms. - The output txt file contrain the radial position of the atom comparing it to the + The output text file containing the radial position of the atom comparing it to the center of mass to the bilayer. Usage: lipid_unwrap.py [OPTIONS] selection_string system_file trajectory_file @@ -67,169 +65,186 @@ def fullhelp(): --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 + --output_traj: output a dcd of the unwrapped trajectory default(--output_traj 0) + --output_prefix: prefix to use for the output default is output - Exmaple: + Example: python /media/bak12/Analysis/unwrap/lipid_unwrap.py --ouput_traj 1 --output_prefix='unwrapped_foo' foo.psf 'resname == "POPC"' foo.dcd - We will obatian a dcd of the unwrapped trajectory with the name of unwrapped_foo.dcd - a pdb of unwrapped_foo.pdb to be used to veiw the dcd, and unwrapped_foo.txt that - has the radial distacne of the ceneter of the lipid to the COM of bilayer + We will obtain a dcd of the unwrapped trajectory with the name of unwrapped_foo.dcd + a pdb of unwrapped_foo.pdb to be used to view the dcd, and unwrapped_foo.txt that + has the radial distantness of the center of the lipid to the COM of bilayer at that frame for the POPC lipids in the membrane""") + + 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 # scaling factors of each - check = True # checks if the condtion is true - while check: # while the ditstance propational to a jump - diff = centers[index] - prev_centers[index] # final minus statrs - tracker = 0 # tracks which if it enterhs - # looking for points that jumpedd - if diff.x() > half_box.x(): - fix[0] -= box.x() * a + """ find Fix is the takes the center before it is unwrapped then checks the lipids + until to see if the new position has moved more than half the box then the box + length is added to the position of that atom to negate the wrapping. If the + corrdinate of fixed atom still is lager than half the box a factor is multiple to + the box length added to the system until the distance is within reason""" + fix = loos.GCoord() # blank GCorrd() == (0,0,0) + bad = {0:1,1:1,2:1,3:1,4:1,5:1} # scaling factors of each + check = True # checks if the condtion is true + while check: # while the ditstance propational to a jump + diff = centers[index] - prev_centers[index] # final minus statrs + tracker = 0 # tracks which if it enterhs + # looking for points that jumped + bad_coord = {} + if diff.x() > half_box.x(): + fix[0] -= box.x() * bad[tracker] + bad_coord[tracker] = bad[tracker] elif diff.x() < -half_box.x(): - fix[0] += box.x() * b - tracker =1 + fix[0] += box.x() * bad[tracker] + tracker = 1 + bad_coord[tracker] = bad[tracker] if diff.y() > half_box.y(): - fix[1] -= box.y() * c + fix[1] -= box.y() * bad[tracker] tracker = 2 + bad_coord[tracker] = bad[tracker] elif diff.y() < -half_box.y(): - fix[1] += box.y() * d + fix[1] += box.y() * bad[tracker] tracker = 3 + bad_coord[tracker] = bad[tracker] if diff.z() > half_box.z(): - fix[2] -= box.z() * e + fix[2] -= box.z() * bad[tracker] tracker = 4 + bad_coord[tracker] = bad[tracker] elif diff.z() < -half_box.z(): - fix[2] += box.z() * f + fix[2] += box.z() * bad[tracker] tracker = 5 - # the new postion is the wrapped postion plus the fix + bad_coord[tracker] = bad[tracker] + # the new position is the wrapped position plus the fix ans = centers[index] + fix - # if the distacne from the old location to the new is large than the - ## half_box if prev_centers[index].distance(ans) > half_box.x(): - # using the value of tracker we find which if was used to fix the - ## bad corrdinate - if tracker == 0: - # if the ans - previous is less than half box the jump - ## was from the postive side of the pbc to the negative - if ans.x() - prev_centers[index].x() < half_box.x(): - a += 1 - elif tracker == 1: - # is the diffrecne is grater than half_box the jump was from - ## negative to postive - 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 + # can assume here that the bad cord is empty + if len(bad_coord) != 0: + for i in bad_coord.keys(): + bad[i] += 1 else: - # if there was no large jump we exit the loop 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) + 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)') +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() 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) +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() +# 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 -if args.output_traj: # fixes the issues where an emptry traj is made when -## not wanting one out - outtraj = loos.DCDWriter(pre +".dcd") +# this is the order of the lipid resname-resid-segname +# this will be useful later +if args.output_traj: # fixes the issues where an emptry traj is made when + # not wanting one out + 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))]) +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_R Lipid_2_center_R Lipid_3_center_R ..." for frame in traj: # loop into the lipids - COM_R = loos.selectAtoms(system,args.selection_string).centerOfMass() + COM_R = loos.selectAtoms(system, args.selection_string).centerOfMass() centers = [] for index in range(len(lipids)): - centers.append(lipids[index].centroid()) # geometric center of the lipid + # geometric center of the lipid + centers.append(lipids[index].centroid()) # position 0 for frame 0 is always the starting position # in frame one nothing should have moved R_coord = [traj.index()] 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 + box = frame.periodicBox() # get a G cord of the PBC + half_box = 0.5 * box # a square box's center is at the half box way between all the positions for index in range(len(lipids)): - fix =findFix() - centers[index]+= fix - R_coord.append(np.sqrt((centers[index].x()-COM_R.x())**2 +(centers[index].y()-COM_R.y())**2+(centers[index].z()-COM_R.z())**2)) - # fix the werid jump issssues if the atoms move more than wanted + fix = findFix() + centers[index] += fix + R_coord.append(np.sqrt((centers[index].x() - COM_R.x())**2 + ( + centers[index].y() - COM_R.y())**2 + (centers[index].z() - COM_R.z())**2)) + # fix the weird jump issues 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 + moved_lipids = loos.AtomicGroup() # blank Atoms group + for lipid, center in zip(lipids, centers): + lipid.centerAtOrigin() # move that lipid to the center # translate to the center - lipid.translate(center) - moved_lipids.append(lipid) + lipid.translate(center) # then move the lipids to the center of + # the unwrapped point + moved_lipids.append(lipid) # add the lipid to the blank Atom loos # 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() + moved_lipids.periodicBox(loos.GCoord(1000000, 1000000, 1000000)) # set the + # PBC to 100000 100000 100000 + if len(lipids) >= 5: # if we have a smaller selection moving the + # everything to the center of mass of the bilayer will not show + # the movement that we are looking for the 5 is kinda arbitrary + moved_lipids.centerAtOrigin() # center the whole selection a the + # origin + outtraj.writeFrame(moved_lipids) # write the frame + if first: # if frame 0 + lipids.periodicBox(loos.GCoord(1000000, 1000000, 1000000)) + # change the box size 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)) + with open(f"{pre}.pdb", "w") as out: + out.write(str(pdb)) # write a pdb all = [traj.index()] if first: for index in range(len(lipids)): - R_coord.append(np.sqrt((centers[index].x()-COM_R.x())**2 +(centers[index].y()-COM_R.y())**2+(centers[index].z()-COM_R.z())**2)) + R_coord.append(np.sqrt((centers[index].x() - COM_R.x())**2 + ( + centers[index].y() - COM_R.y())**2 + (centers[index].z() - COM_R.z())**2)) all_centers = np.array(R_coord) else: R_coord = np.array((R_coord)) - all_centers = np.row_stack((all_centers,R_coord)) + all_centers = np.row_stack((all_centers, R_coord)) prev_centers = centers first = False -np.savetxt(f'{pre}.txt',all_centers,header=header) +np.savetxt(f'{pre}.txt', all_centers, header=header) + +if __name__ in 'main': + pass + From 4e4ccbbf728fec1ba5b7387211dc890e04c9467f Mon Sep 17 00:00:00 2001 From: krb0073 <41389987+krb0073@users.noreply.github.com> Date: Tue, 12 Oct 2021 13:13:55 -0400 Subject: [PATCH 3/6] Update lipid_unwrap.py Used the code given to me by Dr.Grossfeild That moved all no function calls to after the if name in main. editted if name in main to "__main__". Removed old commnet that did not apply anymore around line 215, Added an arugment to not recenter the unwrapped trajectory. Made a 2D array instead of buibling one from scratch over time. --- Packages/PyLOOS/lipid_unwrap.py | 293 +++++++++++++++++--------------- 1 file changed, 154 insertions(+), 139 deletions(-) diff --git a/Packages/PyLOOS/lipid_unwrap.py b/Packages/PyLOOS/lipid_unwrap.py index 3e2798ce1..5dbb6ab8b 100644 --- a/Packages/PyLOOS/lipid_unwrap.py +++ b/Packages/PyLOOS/lipid_unwrap.py @@ -22,10 +22,6 @@ You should have received a copy of the GNU General Public License along with this program. If not, see . """ -# making the fullhelp for the program - - - import sys import loos @@ -33,24 +29,23 @@ import numpy as np import argparse def fullhelp(): - print( - """lipid_unwrap.py takes the coordinates from a trajectory file to unwrap + print("""lipid_unwrap.py takes the coordinates from a trajectory file to unwrap the periodic box 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 - periodic box length (in x,y, or z) a correction factor is added to the current position - to correct the positions it moves outside the boundaries. + periodic box 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 trajectory 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 trajectory,the assumption is made that the - Lx,Ly or Lz of the periodic box will not be greater than 1000000 Angstroms. + If desired, a pdb and DCD of the unwrapped trajectory 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 trajectory, the assumption is made that the + Lx, Ly or Lz of the periodic box will not be greater than 1000000 Angstroms. - The output text file containing the radial position of the atom comparing it to the - center of mass to the bilayer. + The output text file containing the radial position of the atom comparing + it to the center of mass to the bilayer. Usage: lipid_unwrap.py [OPTIONS] selection_string system_file trajectory_file -- selection_string identifies the atoms to be analyzed. They will be @@ -69,22 +64,24 @@ def fullhelp(): --output_prefix: prefix to use for the output default is output Example: - python /media/bak12/Analysis/unwrap/lipid_unwrap.py --ouput_traj 1 --output_prefix='unwrapped_foo' - foo.psf 'resname == "POPC"' foo.dcd + python3 lipid_unwrap.py --ouput_traj 1 --output_prefix='unwrapped_foo' foo.psf 'resname == "POPC"' foo.dcd - We will obtain a dcd of the unwrapped trajectory with the name of unwrapped_foo.dcd - a pdb of unwrapped_foo.pdb to be used to view the dcd, and unwrapped_foo.txt that - has the radial distantness of the center of the lipid to the COM of bilayer - at that frame for the POPC lipids in the membrane""") + We will obtain a dcd of the unwrapped trajectory with the name of + unwrapped_foo.dcd a pdb of unwrapped_foo.pdb to be used to view the dcd, and + unwrapped_foo.txt that has the radial distantness of the center of the lipid + to the COM of bilayer at that frame for the POPC lipids in the membrane + """) def findFix(): - """ find Fix is the takes the center before it is unwrapped then checks the lipids + """ + find Fix is the takes the center before it is unwrapped then checks the lipids until to see if the new position has moved more than half the box then the box length is added to the position of that atom to negate the wrapping. If the - corrdinate of fixed atom still is lager than half the box a factor is multiple to - the box length added to the system until the distance is within reason""" - fix = loos.GCoord() # blank GCorrd() == (0,0,0) + corrdinate of fixed atom still is lager than half the box a factor is multiple to + the box length added to the system until the distance is within reason + """ + fix = loos.GCoord() bad = {0:1,1:1,2:1,3:1,4:1,5:1} # scaling factors of each check = True # checks if the condtion is true while check: # while the ditstance propational to a jump @@ -94,27 +91,28 @@ def findFix(): bad_coord = {} if diff.x() > half_box.x(): fix[0] -= box.x() * bad[tracker] - bad_coord[tracker] = bad[tracker] + bad_coord[tracker] = bad[tracker] elif diff.x() < -half_box.x(): - fix[0] += box.x() * bad[tracker] + fix[0] += box.x() * bad[tracker] tracker = 1 bad_coord[tracker] = bad[tracker] if diff.y() > half_box.y(): fix[1] -= box.y() * bad[tracker] tracker = 2 - bad_coord[tracker] = bad[tracker] + bad_coord[tracker] = bad[tracker] elif diff.y() < -half_box.y(): - fix[1] += box.y() * bad[tracker] + fix[1] += box.y() * bad[tracker] tracker = 3 bad_coord[tracker] = bad[tracker] if diff.z() > half_box.z(): fix[2] -= box.z() * bad[tracker] tracker = 4 - bad_coord[tracker] = bad[tracker] + bad_coord[tracker] = bad[tracker] elif diff.z() < -half_box.z(): - fix[2] += box.z() * bad[tracker] + fix[2] += box.z() * bad[tracker] tracker = 5 - bad_coord[tracker] = bad[tracker] + bad_coord[tracker] = bad[tracker] + # the new position is the wrapped position plus the fix ans = centers[index] + fix if prev_centers[index].distance(ans) > half_box.x(): @@ -125,6 +123,8 @@ def findFix(): else: check = False return fix + + class FullHelp(argparse.Action): def __init__(self, option_strings, dest, nargs=None, **kwargs): kwargs['nargs'] = 0 @@ -137,114 +137,129 @@ def __call__(self, parser, namespace, values, option_string=None): 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() -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 useful later -if args.output_traj: # fixes the issues where an emptry traj is made when - # not wanting one out - 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_R Lipid_2_center_R Lipid_3_center_R ..." -for frame in traj: - # loop into the lipids - COM_R = loos.selectAtoms(system, args.selection_string).centerOfMass() - centers = [] - for index in range(len(lipids)): - # geometric center of the lipid - centers.append(lipids[index].centroid()) - # position 0 for frame 0 is always the starting position - # in frame one nothing should have moved - R_coord = [traj.index()] - if not first: - box = frame.periodicBox() # get a G cord of the PBC - half_box = 0.5 * box # a square box's center is at the half box way between all the positions - for index in range(len(lipids)): - fix = findFix() - centers[index] += fix - R_coord.append(np.sqrt((centers[index].x() - COM_R.x())**2 + ( - centers[index].y() - COM_R.y())**2 + (centers[index].z() - COM_R.z())**2)) - # fix the weird jump issues if the atoms move more than wanted - # if the user has the flag on for the trajectory then this is True + +if __name__ in '__main__': + 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)') + parser.add_argument('--no_center', + default=1, + type=int, + help='Do not wrap the trajecotry default is to center the lipid selection') + args = parser.parse_args() + pre = args.output_prefix + header = " ".join([f"{i}" for i in sys.argv]) + + system = loos.createSystem(args.system_file) + traj = loos.pyloos.Trajectory(args.traj_file, + system, + stride=args.stride, + skip=args.skip) + + # select the atoms for the lipids + # split into the individual molecules of lipid + lipid_sel = loos.selectAtoms(system, args.selection_string) + lipids =lipid_sel.splitByMolecule() + # boolean to check if this is frame is frame 0 + first = True + # making a numpy array to save the coorinates + num_lipids = len(lipids) + num_frames= len(traj) + all_centers = np.zeros((num_lipids,num_frames)) + # loop into the frame if args.output_traj: - # loop through the list of center locations and - moved_lipids = loos.AtomicGroup() # blank Atoms group - for lipid, center in zip(lipids, centers): - lipid.centerAtOrigin() # move that lipid to the center - # translate to the center - lipid.translate(center) # then move the lipids to the center of - # the unwrapped point - moved_lipids.append(lipid) # add the lipid to the blank Atom loos - # setting the pbc to be large - moved_lipids.periodicBox(loos.GCoord(1000000, 1000000, 1000000)) # set the - # PBC to 100000 100000 100000 - if len(lipids) >= 5: # if we have a smaller selection moving the - # everything to the center of mass of the bilayer will not show - # the movement that we are looking for the 5 is kinda arbitrary - moved_lipids.centerAtOrigin() # center the whole selection a the - # origin - outtraj.writeFrame(moved_lipids) # write the frame - if first: # if frame 0 - lipids.periodicBox(loos.GCoord(1000000, 1000000, 1000000)) - # change the box size - 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)) # write a pdb - all = [traj.index()] - if first: + outtraj = loos.DCDWriter(pre + ".dcd") + + # this is the order of the lipid resname-resid-segname + 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 += "\nLipid_1_center_R Lipid_2_center_R Lipid_3_center_R ..." + for frame in traj: + # shift the lipids such that their total centroid is at the origin + COM_R = loos.selectAtoms(system, args.selection_string).centroid() + system.translate(-COM_R) + centers = [] for index in range(len(lipids)): - R_coord.append(np.sqrt((centers[index].x() - COM_R.x())**2 + ( - centers[index].y() - COM_R.y())**2 + (centers[index].z() - COM_R.z())**2)) - all_centers = np.array(R_coord) - else: - R_coord = np.array((R_coord)) - all_centers = np.row_stack((all_centers, R_coord)) - prev_centers = centers - first = False -np.savetxt(f'{pre}.txt', all_centers, header=header) - -if __name__ in 'main': - pass + # geometric center of the lipid + centers.append(lipids[index].centroid()) + # position 0 for frame 0 is always the starting position + # in frame one nothing should have moved + R_coord = [traj.index()] + if not first: + box = frame.periodicBox() # get a G cord of the PBC + half_box = 0.5 * box # a square box's center is at the half box way between all the positions + for index in range(len(lipids)): + fix = findFix() + centers[index] += fix + R_coord.append(centers[index].length()) + if args.output_traj: + # loop through the list of center locations and + moved_lipids = loos.AtomicGroup() # blank Atoms group + for lipid, center in zip(lipids, centers): + lipid.centerAtOrigin() # move that lipid to the center + # translate to the center + lipid.translate(center) # then move the lipids to the center of + # the unwrapped point + moved_lipids.append(lipid) # add the lipid to the blank Atom loos + + # setting the pbc to be large + moved_lipids.periodicBox(loos.GCoord(1000000, 1000000, 1000000)) # set the + if args.no_center: + moved_lipids.centerAtOrigin() + + outtraj.writeFrame(moved_lipids) + + # write a pdb + if first: + moved_lipids.periodicBox(loos.GCoord(1000000, 1000000, 1000000)) + # change the box size + frame_copy = lipid_sel.copy() + 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)) + # now loop into the array of R_coord input the values + for index in range(len(lipids)): + # find the lenght of the centers indexs + all_centers[index,traj.index()] = centers[index].length() + """ if first: + for index in range(len(lipids)): + R_coord.append(centers[index].length()) + all_centers = np.array(R_coord) + + # TODO: you're doing this every frame? That's going to be + # insanely slow and wasteful. You're much better off just + # allocating a num_frames x num_lipids 2D array right at the beginning + else: + R_coord = np.array((R_coord)) + all_centers = np.row_stack((all_centers, R_coord))""" + prev_centers = centers + first = False + np.savetxt(f'{pre}.txt', all_centers, header=header) From 1daa840714344d0a2040b086b585fdfde345cad5 Mon Sep 17 00:00:00 2001 From: krb0073 <41389987+krb0073@users.noreply.github.com> Date: Tue, 12 Oct 2021 13:48:26 -0400 Subject: [PATCH 4/6] Update lipid_unwrap.py added a commnet for the fullhelp for the no_center flag --- Packages/PyLOOS/lipid_unwrap.py | 1 + 1 file changed, 1 insertion(+) diff --git a/Packages/PyLOOS/lipid_unwrap.py b/Packages/PyLOOS/lipid_unwrap.py index 5dbb6ab8b..3f637faa4 100644 --- a/Packages/PyLOOS/lipid_unwrap.py +++ b/Packages/PyLOOS/lipid_unwrap.py @@ -62,6 +62,7 @@ def fullhelp(): --fullhelp: prints this message --output_traj: output a dcd of the unwrapped trajectory default(--output_traj 0) --output_prefix: prefix to use for the output default is output + --no_center: Variable to turn off the recentering of the trajectory default(--no_center 1) Example: python3 lipid_unwrap.py --ouput_traj 1 --output_prefix='unwrapped_foo' foo.psf 'resname == "POPC"' foo.dcd From 44fec01bae8258db4af1d1c8318baa99febdf4bc Mon Sep 17 00:00:00 2001 From: krb0073 <41389987+krb0073@users.noreply.github.com> Date: Tue, 12 Oct 2021 13:50:20 -0400 Subject: [PATCH 5/6] Update lipid_unwrap.py removed old way of getting the matrix for the corrianted of the lipid centers. --- Packages/PyLOOS/lipid_unwrap.py | 11 ----------- 1 file changed, 11 deletions(-) diff --git a/Packages/PyLOOS/lipid_unwrap.py b/Packages/PyLOOS/lipid_unwrap.py index 3f637faa4..049953f7c 100644 --- a/Packages/PyLOOS/lipid_unwrap.py +++ b/Packages/PyLOOS/lipid_unwrap.py @@ -250,17 +250,6 @@ def __call__(self, parser, namespace, values, option_string=None): for index in range(len(lipids)): # find the lenght of the centers indexs all_centers[index,traj.index()] = centers[index].length() - """ if first: - for index in range(len(lipids)): - R_coord.append(centers[index].length()) - all_centers = np.array(R_coord) - - # TODO: you're doing this every frame? That's going to be - # insanely slow and wasteful. You're much better off just - # allocating a num_frames x num_lipids 2D array right at the beginning - else: - R_coord = np.array((R_coord)) - all_centers = np.row_stack((all_centers, R_coord))""" prev_centers = centers first = False np.savetxt(f'{pre}.txt', all_centers, header=header) From c2aa0041bffe48f6a8bed384541291efd08e8247 Mon Sep 17 00:00:00 2001 From: krb0073 <41389987+krb0073@users.noreply.github.com> Date: Tue, 12 Oct 2021 20:29:41 -0400 Subject: [PATCH 6/6] Update lipid_unwrap.py Fixed whitespace error --- Packages/PyLOOS/lipid_unwrap.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Packages/PyLOOS/lipid_unwrap.py b/Packages/PyLOOS/lipid_unwrap.py index 049953f7c..a83becf21 100644 --- a/Packages/PyLOOS/lipid_unwrap.py +++ b/Packages/PyLOOS/lipid_unwrap.py @@ -247,9 +247,9 @@ def __call__(self, parser, namespace, values, option_string=None): with open(f"{pre}.pdb", "w") as out: out.write(str(pdb)) # now loop into the array of R_coord input the values - for index in range(len(lipids)): - # find the lenght of the centers indexs - all_centers[index,traj.index()] = centers[index].length() + for index in range(len(lipids)): + # find the lenght of the centers indexs + all_centers[index,traj.index()] = centers[index].length() prev_centers = centers first = False np.savetxt(f'{pre}.txt', all_centers, header=header)