我有个问题。我从这里使用脚本https://github.com/NMRLipids/MATCH/tree/master/scripts 我使用 Python 2.7（Anaconda，我安装了 MDAnalysis，numpy 等）我首先使用这个脚本
#!/bin/bash
# Gromacs 5.x version of ...
# bash wrapper for calculating Order Parameters of lipid bilayer
# python script/library calcOrderParameters.py
# meant for use with NMRlipids projects
#------------------------------------------------------------
# Made by J.Melcr, Last edit 2017/03/21
#------------------------------------------------------------
#
# Generation of non-water trajectory removed by O.H.S Ollila, 2018/08/23
#
#scriptdir='../../../../../scripts/'
traj_file_name="eq2_pbc_whole_290_300.xtc"
tpr_file_name="eq2.tpr"
conf_file_name="eq2.gro"
op_def_file="defFILE.def"
op_out_file="OrdPars_pal.dat"
#top="topol.top"
traj_pbc_file_name="eq2_pbc_whole_290_300_mol.xtc"
#f_conc=55430 # in mM/L
if ! [ -s $tpr_file_name ]
then
echo "We really need " $tpr_file_name " , but we can't find it!"
exit 1
fi
# remove PBC:
#! [ -s $traj_pbc_file_name ] && echo System | gmx trjconv -f $traj_file_name -s $tpr_file_name -o $traj_pbc_file_name -pbc mol
#CALCULATE ORDER PARAMETERS
python calcOrderParameters.py -i $op_def_file -t $conf_file_name -x $traj_pbc_file_name -o $op_out_file
然后我使用了这个脚本（因为 order parameters_calculate.sh 的最后一行是
python calcOrderParameters.py -i $op_def_file -t $conf_file_name -x $traj_pbc_file_name -o $op_out_file
这是 calcOrderParameters.py 的完整代码
"""
calculation of order parameters
from a MD trajectory
useful for example for lipid bilayers
meant for use with NMRlipids projects
------------------------------------------------------------
Made by J. Melcr, Last edit 2018/03/26
with contributions from:
H. Antila
------------------------------------------------------------
input: Order parameter definitions [order_parameter_file],
gro/tpr [topology_file], and xtc/trr/pbb file [trajectory_file] (or equivalents)
output: order parameters [output_file] (2 textfiles)
usage : python calcOrderParameters.py -i [order_parameter_file] -t [topology_file] -x [trajectory_file] -o [output_file]
--------------------------------------------------------
[order_parameter_file] should contain the definitions of order parameters in format:
OP_name1 Residue Carbon_name Hydrogen_name
OP_name2 Residue Carbon_name Hydrogen_name
.....
Example (CHARMM36):
beta1 POPC C12 H12A
beta2 POPC C12 H12B
alpha1 POPC C11 H11A
alpha2 POPC C11 H11B
g3_1 POPC C1 HA
g3_2 POPC C1 HB
g2_1 POPC C2 HS
g1_1 POPC C3 HX
g1_2 POPC C3 HY
"""
# coding: utf-8
import MDAnalysis as mda
import numpy as np
import math
import os, sys
from optparse import OptionParser
from collections import OrderedDict
bond_len_max=1.5 # in Angstroms, max distance between atoms for reasonable OP calculation (PBC and sanity check)
bond_len_max_sq=bond_len_max**2
class OrderParameter:
"""
Class for storing&manipulating
order parameter (OP) related metadata (definition, name, ...)
and OP trajectories
and methods to evaluate OPs.
OP definition consist of:
- name of the OP
- residue name
- involved atoms (exactly 2)
+ extra: mean, std.dev. & err. estimate
(using standard error of the means from individual residues)
of the OP (when reading-in an already calculated result)
Example (CHARMM36):
beta1 POPC C12 H12A
beta2 POPC C12 H12B
alpha1 POPC C11 H11A
alpha2 POPC C11 H11B
g3_1 POPC C1 HA
g3_2 POPC C1 HB
g2_1 POPC C2 HS
g1_1 POPC C3 HX
g1_2 POPC C3 HY
"""
def __init__(self, name, resname, atom_A_name, atom_B_name, *args):
"""
Initialization of an instance of this class.
it doesn't matter which atom comes first,
atom A or B, for OP calculation.
"""
self.name = name # name of the order parameter, a label
self.resname = resname # name of residue atoms are in
self.atAname = atom_A_name
self.atBname = atom_B_name
# variables for error estimate -- standard error of the mean (STEM)
self.avg = None # average/mean value from all residues
self.means = None # list of mean values from each individual residue
self.std = None # standard deviation (sqrt(variance))
self.stem = None # STandard Error of the Mean
# trajectory as list
self.traj = [] # for storing OPs
for field in self.__dict__:
if not isinstance(field, str):
raise UserWarning, "provided name >> {} << is not a string! \n \
Unexpected behaviour might occur.".format(field)
else:
if not field.strip():
raise RuntimeError, "provided name >> {} << is empty! \n \
Cannot use empty names for atoms and OP definitions.".format(field)
# extra optional arguments allow setting avg,std values -- suitable for reading-in results of this script
if len(args) == 2:
self.avg = args[0]
self.std = args[1]
elif len(args) == 3:
self.avg = args[0]
self.std = args[1]
self.stem = args[2]
else:
if len(args) != 0:
raise UserWarning, "Number of optional positional arguments is {len}, not 3, 2 or 0. Args: {args}\nWrong file format?".format(len=len(args), args=args)
def calc_OP(self, atoms):
"""
calculates Order Parameter according to equation
S = 1/2 * (3*cos(theta)^2 -1)
"""
vec = atoms[1].position - atoms[0].position
d2 = np.square(vec).sum()
if d2>bond_len_max_sq:
raise UserWarning, "Atomic distance for atoms \
{at1} and {at2} in residue no. {resnr} is suspiciously \
long: {d}!\nPBC removed???".format(at1=atoms[0].name, at2=atoms[1].name, resnr=atoms[0].resid, d=math.sqrt(d2))
cos2 = vec[2]**2/d2
S = 0.5*(3.0*cos2-1.0)
return S
def calc_angle(self, atoms, z_dim=45.0):
"""
calculates the angle between the vector and z-axis in degrees
no PBC check!
assuming a sim-box-centred membrane --> it's centre ~ z_dim/2
Warning: user has to make sure that correct z_dim is supplied,
otherwise - This is a bit DIRTY!!
-- this is taken care of in the main trajectory reader in this module
"""
vec = atoms[1].position - atoms[0].position
d = math.sqrt(np.square(vec).sum())
cos = vec[2]/d
# values for the bottom leaflet are inverted so that
# they have the same nomenclature as the top leaflet
cos *= math.copysign(1.0, atoms[0].position[2]-z_dim*0.5)
try:
angle = math.degrees(math.acos(cos))
except ValueError:
if abs(cos)>=1.0:
print "Cosine is too large = {} --> truncating it to +/-1.0".format(cos)
cos = math.copysign(1.0, cos)
angle = math.degrees(math.acos(cos))
return angle
@property
def get_avg_std_OP(self):
"""
Provides average and stddev of all OPs in self.traj
This method becomes deprecated after the introduction of error estimation
"""
# convert to numpy array
return (np.mean(self.traj), np.std(self.traj))
@property
def get_avg_std_stem_OP(self):
"""
Provides average, stddev and standard error of mean for all OPs in self.traj
"""
self.means = np.mean(self.traj, axis=0)
return ( np.mean(self.traj),
np.std(self.means),
np.std(self.means)/np.sqrt(len(self.means)) )
def read_trajs_calc_OPs(ordPars, top, trajs):
"""
procedure that
creates MDAnalysis (mda) Universe instance with topology top,
reads in trajectories trajs and then
goes through every frame and
evaluates each Order Parameter "S" from the list of OPs ordPars.
ordPars : list of OrderParameter class instances
each item in this list describes an Order parameter to be calculated in the trajectory
top : str
filename of a top file (e.g. conf.gro)
trajs : list of strings
filenames of trajectories
"""
# read-in topology and trajectory
mol = mda.Universe(top, trajs)
# make atom selections for each OP and store it as its attribute for later use with trajectory
for op in ordPars.values():
# selection = pairs of atoms, split-by residues
# this selection format preserves the order of the atoms (atA, atB) independent of their order in the topology
selection = mol.select_atoms("resname {rnm} and name {atA}".format(
rnm=op.resname, atA=op.atAname),
"resname {rnm} and name {atB}".format(
rnm=op.resname, atB=op.atBname)
).atoms.split("residue")
for res in selection:
# check if we have only 2 atoms (A & B) selected
if res.n_atoms != 2:
print res.resnames, res.resids
for atom in res.atoms:
print atom.name, atom.id
raise UserWarning, "Selection >> name {atA} {atB} << \
contains {nat} atoms, but should contain exactly 2!".format(
atA=op.atAname, atB=op.atBname, nat=res.n_atoms)
op.selection = selection
#Nres=len(op.selection)
#Nframes=len(mol.trajectory)
# go through trajectory frame-by-frame
# and calculate each OP from the list of OPs
# for each residue separately
for frame in mol.trajectory:
for op in ordPars.values():
# temporary list of order parameters for
# each individual residue for the given frame
temp_S = []
for residue in op.selection:
if "vec" in op.name:
S = op.calc_angle(residue, z_dim=frame.dimensions[2])
else:
S = op.calc_OP(residue)
temp_S.append(S)
# resulting S-trajectory will be a list of lists
# so that individual residues can be easily distinquished
op.traj.append(temp_S)
def parse_op_input(fname):
"""
parses input file with Order Parameter definitions
file format is as follows:
OP_name resname atom1 atom2 +extra: OP_mean OP_std
(flexible cols)
fname : string
input file name
returns : dictionary
with OrderParameters class instances
"""
# Using ordered dict since it preserves the read-in order. Might come in handy when comparing to experiments.
ordPars = OrderedDict()
try:
with open(fname,"r") as f:
for line in f.readlines():
if not line.startswith("#"):
items = line.split()
ordPars[items[0]] = OrderParameter(*items)
except:
raise RuntimeError, "Couldn't read input file >> {inpf} <<".format(inpf=opts.inp_fname)
return ordPars
#%%
if __name__ == "__main__":
# help message is automatically provided
# type=string, action=store is default
parser = OptionParser()
parser.add_option('-i', '--inp', dest='inp_fname', help='input (OP definitions) file name', default="Headgroup_Glycerol_OPs.def")
parser.add_option('-t', '--top', dest='top_fname', help='topology (gro, pdb) file name', default="last_frame_nonwat.gro")
parser.add_option('-x', '--traj', dest='traj_fname', help='beginning of trajectory (xtc, dcd) files names (will use all files beginning with this string).', default="traj")
parser.add_option('-o', '--out', dest='out_fname', help='output (OPs mean&std) file name', default="Headgroup_Glycerol_OPs.dat")
opts, args = parser.parse_args()
# dictionary for storing of OrderParameter class instances (name-wise, of course)
print "\nReading OP definitions ...\n"
ordPars = parse_op_input(opts.inp_fname)
# get all parts of trajectories
trajs = []
for file_name in os.listdir(os.getcwd()):
if file_name.startswith(opts.traj_fname):
trajs.append(file_name)
# read trajectory and calculate all OPs
print "Reading trajectories and calculating OPs ...\n"
read_trajs_calc_OPs(ordPars, opts.top_fname, trajs)
print "OP Name mean std err.est."
print "--------------------------------------------------------------------"
for op in ordPars.values():
try:
(op.avg, op.std, op.stem) = op.get_avg_std_stem_OP
print "{:10s} {: 2.4f} {: 2.4f} {: 2.4f}".format(op.name, op.avg, op.std, op.stem)
except:
print "{:s} -- problem calculating statistics ".format(op.name)
print "--------------------------------------------------------------------"
try:
with open(opts.out_fname,"w") as f:
f.write("# OP_name resname atom1 atom2 OP_mean OP_stddev OP_stem\n\
#--------------------------------------------------------------------\n")
for op in ordPars.values():
f.write( "{:20s} {:7s} {:5s} {:5s} {: 2.5f} {: 2.5f} {: 2.5f} \n".format(
op.name, op.resname, op.atAname, op.atBname,
op.avg, op.std, op.stem)
)
print "\nOrderParameters written to >> {fname} <<".format(fname=opts.out_fname)
except:
print "ERROR: Problems writing main output file."
# this single-line format may become soon deprecated, but
# it is the format that is used in NMRlipids projects for processing through awk+gnuplot
try:
conc_formatted_line = "conc {b1: 2.6f} {b1e: 2.6f} {b2: 2.6f} {b2e: 2.6f} {a1: 2.6f} {a1e: 2.6f} {a2: 2.6f} {a2e: 2.6f}".format(
b1=ordPars['beta1'].avg, b2=ordPars['beta2'].avg,
a1=ordPars['alpha1'].avg, a2=ordPars['alpha2'].avg,
b1e=ordPars['beta1'].stem, b2e=ordPars['beta2'].stem,
a1e=ordPars['alpha1'].stem, a2e=ordPars['alpha2'].stem)
print
print "Single line format:\nconc beta1 err beta2 err alpha1 err alpha2 err"
print conc_formatted_line
with open(opts.out_fname+".line","w") as f:
f.write(conc_formatted_line)
except:
print "ERROR: Problems writing the beta-alpha single line format file."
然后我得到一个错误
./order_parameters_calculate.sh
Reading OP definitions ...
Traceback (most recent call last):
File "calcOrderParameters.py", line 278, in <module>
ordPars = parse_op_input(opts.inp_fname)
File "calcOrderParameters.py", line 259, in parse_op_input
raise RuntimeError, "Couldn't read input file >> {inpf} <<".format(inpf=opts.inp_fname)
RuntimeError: Couldn't read input file >> defFILE.def <<
^C
我不知道为什么我有这个错误。我有 defFILE.def - 当我在终端中写 ls 时，我有
(mda_pyt27) jakub@jakub-Z370-HD3P:/media/jakub/WD/EnerPres_g2020_2/bez_zmian/POPC_300K_11_01_2021_bez_zmian_parametrow/calc_order$ ls
calcOrderParameters.py eq2_pbc_whole_290_300.xtc
defFILE.def eq2.tpr
eq2.gro order_parameters_calculate.sh
eq2_pbc_whole_290_300_mol.xtc
我的 defFile.def：
p1_1 POPC C12 H21
p1_2 POPC C12 H22
p2_1 POPC C29 H54
p2_2 POPC C29 H55
p3_1 POPC C30 H56
p3_2 POPC C30 H57
p4_1 POPC C31 H58
p4_2 POPC C31 H59
p5_1 POPC C32 H60
p5_2 POPC C32 H61
p6_1 POPC C33 H62
p6_2 POPC C33 H63
p7_1 POPC C34 H64
p7_2 POPC C34 H65
p8_1 POPC C35 H66
p8_2 POPC C35 H67
p9_1 POPC C36 H68
p9_2 POPC C36 H69
p10_1 POPC C37 H70
p10_2 POPC C37 H71
p11_1 POPC C38 H72
p11_2 POPC C38 H73
p12_1 POPC C39 H74
p12_2 POPC C39 H75
p13_1 POPC C40 H76
p13_2 POPC C40 H77
p14_1 POPC C41 H78
p14_2 POPC C41 H79
p15_1 POPC C42 H80
p15_2 POPC C42 H81
p15_3 POPC C42 H82
可能是什么问题？？？