Here is what I did to compute average HBond distance data for use in a
network analysis project recently.
I did this using jupyter notebook, so I have pasted the commands from each
cell I ran below.
It should work as is provided you update the variables in the second block
of commands to point to your trajectory and topology files.
#Load needed packages
import pytraj as pt
import numpy as np
import pandas as pd
import collections
from collections import OrderedDict as od
import sys
import re
###UPDATE THIS BLOCK ACCORDINGLY###
#Setup paths to trajectory(ies) and topology
strucBaseDir='/home/wesley/work/igps_network_analysis/Analysis/structures/aligned_trajectories/Apo/'
trajDir=strucBaseDir+"/rep_1"
topFileName="step3_charmm2amber.psf"
trajFileName="all_production.aligned.dcd"
topFilePath=strucBaseDir+"/"+topFileName
trajFilePath=trajDir+"/"+trajFileName
iterTraj=pt.iterload(trajFilePath,top=topFilePath)
iterTraj
### ### ###
#compute hbond for loaded trajectory
#You can update this to match your system
#It is essentiall just a wrapper for the equivalent cpptraj command
hbond_data=pt.compute('hbond :1-454 dist 4.0 printatomnum
series',iterTraj,n_cores=4)
print hbond_data.keys()[1]
print hbond_data[hbond_data.keys()[1]]
#Find hbonds that existed for at least 1 frame...
hbList=[]
hbFreqs=[]
nEntries=len(hbond_data.keys())
print 'computing hbond frequencies |',
for iKey in np.arange(nEntries):
if iKey > 0:
if iKey%np.floor(nEntries/10)==0:
print '.',
sys.stdout.flush()
hbFreq=np.sum(hbond_data[hbond_data.keys()[iKey]])*1.0/iterTraj.n_frames
if hbFreq>0: #can set higher frequency threshold if desired
hbList.append(hbond_data.keys()[iKey])
hbFreqs.append(hbFreq)
print '|'
print 'Number of hbonds detected: '+str(len(hbList))
print hbList[0:9]
print hbFreqs[0:9]
#Use list of hbonds with non-zero frequency to compile a list
#of distance commands to be run
distCommands=[]
for hbString in hbList:
hbMasks=re.sub("[A-Z]*_",":",hbString).split("-")
distCommand="distance "+hbMasks[0]+" "+hbMasks[1]
distCommands.append(distCommand)
distCommands[0:5]
#Run the distance commands
dist_data=pt.compute(distCommands,iterTraj,n_cores=4)
print dist_data.keys()[0]+" = "+hbList[0]+" : ",
print dist_data[dist_data.keys()[0]]
#Compute average distances and compile into a pandas.DataFrame
#then export as an ascii data table (or csv format if desired)
averageDistances=pd.DataFrame({"HBond":hbList})
tempList=[]
for iKey in np.arange(len(hbList)):
tempList.append(np.mean(dist_data[dist_data.keys()[iKey]]))
averageDistances['Mean_Distance']=tempList
print averageDistances.head()
averageDistances.to_csv("average_hbond_distances.dat",sep=" ",index=False)
On Tue, Jul 24, 2018 at 10:16 AM, Sally Pias <sallypias.gmail.com> wrote:
> Hi Dan, Wesley, and Hai -
>
> Thank you all for your suggestions. The hbond command detected
> charge-dipole interactions between phosphate O acceptors and O-H or N-H
> donors. Nativecontacts also helped find headgroup phosphate-choline
> interactions
> -- achieved by setting mask1 and mask2 to the phosphate P and the choline
> N. One other functionality that would be useful to me would be to record
> the distance between atoms in the interactions detected. I want to know the
> average headgroup-headgroup and headgroup-organic molecule interaction
> distance in different bilayers.
>
> Best,
> Sally
>
> --
> Sally Pias
> Assistant Professor of Chemistry
> Faculty Adjunct, Department of Biology
> New Mexico Tech
> (575) 835-6204
> www.nmt.edu/~spias
> spias.nmt.edu
>
> On Fri, Jul 20, 2018 at 6:01 PM, Hai Nguyen <nhai.qn.gmail.com> wrote:
>
> > Just for clarity: pytraj works because it is cpptraj under the hood.
> >
> > Hai
> >
> > On Fri, Jul 20, 2018 at 7:57 PM Wesley Michael Botello-Smith <
> > wmsmith.uci.edu> wrote:
> >
> > > nativecontacts works... pytraj can do this as well. If you want to
> ensure
> > > it is for ion to ion or ion to dipole, you can accomplish this using
> > > appropriate selection masks. I think, if you use the convert the
> contact
> > > map into a matrix, you can use that to compute aggregate cluster sizes.
> > >
> > > On Fri, Jul 13, 2018, 3:19 PM Sally Pias <sallypias.gmail.com> wrote:
> > >
> > > > Hi all,
> > > >
> > > > I am looking for a way to count contacts between phospholipid
> > headgroups
> > > > and organic-lipid conjugate molecules embedded in a bilayer. Is
> there a
> > > > script or command similar to hbond that can calculate charge-charge
> and
> > > > ion-dipole contacts between nonsolvent residues (i.e., not water or
> > > > monatomic ions)? Aggregate contact counts, as well as number of
> > contacts
> > > > over time would both be useful. Ideally, we need to be able to
> identify
> > > the
> > > > residues and atoms involved.
> > > >
> > > > Thanks,
> > > > Sally
> > > >
> > > > --
> > > > Sally Pias
> > > > Associate Professor of Chemistry
> > > > Faculty Adjunct, Department of Biology
> > > > New Mexico Tech
> > > > _______________________________________________
> > > > AMBER mailing list
> > > > AMBER.ambermd.org
> > > > http://lists.ambermd.org/mailman/listinfo/amber
> > > >
> > > _______________________________________________
> > > AMBER mailing list
> > > AMBER.ambermd.org
> > > http://lists.ambermd.org/mailman/listinfo/amber
> > >
> > _______________________________________________
> > AMBER mailing list
> > AMBER.ambermd.org
> > http://lists.ambermd.org/mailman/listinfo/amber
> >
> _______________________________________________
> AMBER mailing list
> AMBER.ambermd.org
> http://lists.ambermd.org/mailman/listinfo/amber
>
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Tue Jul 24 2018 - 12:30:02 PDT
