From 5339692b5d0d0c5b077ac517c3dc4cab83976fba Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Wed, 2 May 2018 11:40:53 +0200 Subject: [PATCH 01/23] dihedral distribution script Script to compute either the dihedral or the improper distributions of a system. The computed angle is signed, i.e. the range is between -180 and 180 degrees. --- scripts/dihedral_distribute.py | 1434 ++++++++++++++++++++++++++++++++ 1 file changed, 1434 insertions(+) create mode 100644 scripts/dihedral_distribute.py diff --git a/scripts/dihedral_distribute.py b/scripts/dihedral_distribute.py new file mode 100644 index 0000000..dc5864e --- /dev/null +++ b/scripts/dihedral_distribute.py @@ -0,0 +1,1434 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + Pizza.py/dihedral_distribute.py at 217539187969af996b81ad06628b8a905fb2cbc9 · evoyiatzis/Pizza.py + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+ Skip to content +
+ + + + + + + + + + + + +
+ +
+ +
+
+ + + +
+
+
+ + + + + + + + +
+
+ +
    +
  • +
    + +
    + + + +
    +
    +
    + + Notifications +
    + + + +
    +
    +
    +
    +
    • + +
  • + +
    +
    + + + +
    +
    + + + +
    + +
    • + +
  • +
    + +
    + +
    • +
+ +

+ + /Pizza.py + + + forked from slitvinov/Pizza.py + +

+ +
+ + + + +
+ +
+
+ + + Permalink + + + +
+ +
+ + +
+ +
+
+ + Switch branches/tags +
+ +
+
+ +
+
+ +
+
+ + + +
+
+ + +
+ +
Nothing to show
+
+ +
+
+
+ +
+ + Find file + + + Copy path + +
+ +
+ + + +
+ Fetching contributors… +
+ +
+ + Cannot retrieve contributors at this time +
+ + + +
+
+
+ +
+ Raw + Blame + History +
+ + + + + +
+ +
+ 160 lines (127 sloc) + + 4.62 KB +
+
+ + + +
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
#!/usr/bin/python
+
# Script: dihedral_distribute.py
# Purpose: binned signed dihedral distributions by dihedral type
# Syntax: dihedral_distribute.py datafile nbin theta_min theta_max outfile files ...
# datafile = lammps data file
# nbin = # of bins per dihedral type
# theta_min = min expected angle
# theta_max = max expected angle length
# outfile = file to write stats to
# files = series of dump files
# Example: dihedral_distribute.py pore.data 1000 110. 120. angles.out pore.dump.1
# Author: Evangelos Voyiatzis (TU Darmstadt)
+
# enable script to run from Python directly w/out Pizza.py
+
import sys
import math
from dump import dump
+
if not globals().has_key("argv"): argv = sys.argv
+
def PBC(distance,boxlength): # function to compute the minimum image
+
if abs(distance) > 0.5*xprd:
if distance < 0.0:
distance += boxlength
else:
distance -= boxlength
+
return distance
+
def InnerProduct(Vector1,Vector2): # function to compute the inner product of two vectors
+
return (Vector1[0]*Vector2[0] + Vector1[1]*Vector2[1] + Vector1[2]*Vector2[2])
+
def CrossProduct(Vector1,Vector2): # function to compute the Cross product of two vectors
+
Vector3 = [0] *3
Vector3[0] = Vector1[1]*Vector2[2] - Vector1[2]*Vector2[1]
Vector3[1] = Vector1[2]*Vector2[0] - Vector1[0]*Vector2[2]
Vector3[2] = Vector1[0]*Vector2[1] - Vector1[1]*Vector2[0]
+
return Vector3
+
# main script
+
if len(argv) < 8:
raise StandardError, \
"Syntax: dihedral_distribute.py dihedrals/impropers datafile nbin theta_min theta_max outfile files ..."
+
dt = data(argv[2])
nbins = int(argv[3])
theta_min = float(argv[4])
theta_max = float(argv[5])
outfile = argv[6]
files = ' '.join(argv[7:])
+
# get the angles from the data file
if == "dihedrals":
angle = dt.get("Dihedrals")
else if == "impropers":
angle = dt.get("Impropers")
else:
raise StandardError, "The second keyword is neither 'dihedrals' nor 'impropers' "
sys.exit()
+
nangles = len(angle)
atype = nangles * [0]
iatom = nangles * [0]
jatom = nangles * [0]
katom = nangles * [0]
latom = nangles * [0]
for i in xrange(nangles):
atype[i] = int(angle[i][1] - 1)
iatom[i] = int(angle[i][2] - 1)
jatom[i] = int(angle[i][3] - 1)
katom[i] = int(angle[i][4] - 1)
latom[i] = int(angle[i][5] - 1)
+
ntypes = 0
for i in xrange(nangles): ntypes = max(angle[i][1],ntypes)
ntypes = int(ntypes)
ncount = ntypes * [0]
bin = nbins * [0]
for i in xrange(nbins):
bin[i] = ntypes * [0]
+
# read snapshots one-at-a-time
+
d = dump(files,0)
d.map(1,"id",2,"type",3,"x",4,"y",5,"z")
+
while 1:
time = d.next()
if time == -1: break
box = (d.snaps[-1].xlo,d.snaps[-1].ylo,d.snaps[-1].zlo,
d.snaps[-1].xhi,d.snaps[-1].yhi,d.snaps[-1].zhi)
xprd = box[3] - box[0]
yprd = box[4] - box[1]
zprd = box[5] - box[2]
d.unscale()
d.sort()
x,y,z = d.vecs(time,"x","y","z")
+
VectorA = [0] * 3
VectorB = [0] * 3
VectorC = [0] * 3
for i in xrange(nangles):
# first vector
VectorA[0] = PBC(x[jatom[i]] - x[iatom[i]],xprd)
VectorA[1] = PBC(y[jatom[i]] - y[iatom[i]],yprd)
VectorA[1] = PBC(z[jatom[i]] - z[iatom[i]],zprd)
+
# second vector
VectorB[0] = PBC(x[katom[i]] - x[jatom[i]],xprd)
VectorB[1] = PBC(y[katom[i]] - y[jatom[i]],yprd)
VectorB[2] = PBC(z[katom[i]] - z[jatom[i]],zprd)
normVectorB = math.sqrt(InnerProduct(VectorB,VectorB))
+
# third vector
VectorC[0] = PBC(x[latom[i]] - x[katom[i]],xprd)
VectorC[1] = PBC(y[latom[i]] - y[katom[i]],yprd)
VectorC[2] = PBC(z[latom[i]] - z[katom[i]],zprd)
# compute the signed dihedral angle
HelpVector = CrossProduct(CrossProduct(VectorA,VectorB),CrossProduct(VectorB,VectorC))
Argument1 = InnerProduct(HelpVector,VectorB/normVectorB)
Argument2 = InnerProduct(CrossProduct(VectorA,VectorB),CrossProduct(VectorB,VectorC))
theta = 180.0*math.atan2(Argument1,Argument2)/math.pi
ibin = int(nbins*(theta - theta_min)/(theta_max - theta_min) + 0.5)
if ((ibin >= 0) and (ibin <= nbins-1)):
bin[ibin][atype[i]] += nbins
ncount[atype[i]] += 1
else:
print "Warning: dihedral outside specified range"
print "dihedral type:", atype[i]+1
print "dihedral number:", i
print time,
print
print "Printing normalized dihedral distributions to",outfile
fp = open(outfile,"w")
theta_range = theta_max - theta_min
for i in xrange(nbins):
print >>fp, theta_min + theta_range*float(i)/float(nbins),
for j in xrange(ntypes):
if (ncount[j] > 0):
print >>fp, float(bin[i][j])/float(ncount[j])/theta_range,
else:
print >>fp, 0.0,
print >>fp
fp.close()
+ + + +
+ +
+ + + + + +
+
+
+ +
+
+ +
+ + + + + + +
+ + + You can’t perform that action at this time. +
+ + + + + + + + + + +
+ + You signed in with another tab or window. Reload to refresh your session. + You signed out in another tab or window. Reload to refresh your session. +
+ + + + +
+ Press h to open a hovercard with more details. +
+ + + + + From c1c74356c5b25b25946ba328b485046a0c6f739e Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Wed, 2 May 2018 19:46:58 +0200 Subject: [PATCH 02/23] Delete dihedral_distribute.py --- scripts/dihedral_distribute.py | 1434 -------------------------------- 1 file changed, 1434 deletions(-) delete mode 100644 scripts/dihedral_distribute.py diff --git a/scripts/dihedral_distribute.py b/scripts/dihedral_distribute.py deleted file mode 100644 index dc5864e..0000000 --- a/scripts/dihedral_distribute.py +++ /dev/null @@ -1,1434 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Pizza.py/dihedral_distribute.py at 217539187969af996b81ad06628b8a905fb2cbc9 · evoyiatzis/Pizza.py - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- Skip to content -
- - - - - - - - - - - - -
- -
- -
-
- - - -
-
-
- - - - - - - - -
-
- -
    -
  • -
    - -
    - - - -
    -
    -
    - - Notifications -
    - - - -
    -
    -
    -
    -
    • - -
  • - -
    -
    - - - -
    -
    - - - -
    - -
    • - -
  • -
    - -
    - -
    • -
- -

- - /Pizza.py - - - forked from slitvinov/Pizza.py - -

- -
- - - - -
- -
-
- - - Permalink - - - -
- -
- - -
- -
-
- - Switch branches/tags -
- -
-
- -
-
- -
-
- - - -
-
- - -
- -
Nothing to show
-
- -
-
-
- -
- - Find file - - - Copy path - -
- -
- - - -
- Fetching contributors… -
- -
- - Cannot retrieve contributors at this time -
- - - -
-
-
- -
- Raw - Blame - History -
- - - - - -
- -
- 160 lines (127 sloc) - - 4.62 KB -
-
- - - -
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
#!/usr/bin/python
-
# Script: dihedral_distribute.py
# Purpose: binned signed dihedral distributions by dihedral type
# Syntax: dihedral_distribute.py datafile nbin theta_min theta_max outfile files ...
# datafile = lammps data file
# nbin = # of bins per dihedral type
# theta_min = min expected angle
# theta_max = max expected angle length
# outfile = file to write stats to
# files = series of dump files
# Example: dihedral_distribute.py pore.data 1000 110. 120. angles.out pore.dump.1
# Author: Evangelos Voyiatzis (TU Darmstadt)
-
# enable script to run from Python directly w/out Pizza.py
-
import sys
import math
from dump import dump
-
if not globals().has_key("argv"): argv = sys.argv
-
def PBC(distance,boxlength): # function to compute the minimum image
-
if abs(distance) > 0.5*xprd:
if distance < 0.0:
distance += boxlength
else:
distance -= boxlength
-
return distance
-
def InnerProduct(Vector1,Vector2): # function to compute the inner product of two vectors
-
return (Vector1[0]*Vector2[0] + Vector1[1]*Vector2[1] + Vector1[2]*Vector2[2])
-
def CrossProduct(Vector1,Vector2): # function to compute the Cross product of two vectors
-
Vector3 = [0] *3
Vector3[0] = Vector1[1]*Vector2[2] - Vector1[2]*Vector2[1]
Vector3[1] = Vector1[2]*Vector2[0] - Vector1[0]*Vector2[2]
Vector3[2] = Vector1[0]*Vector2[1] - Vector1[1]*Vector2[0]
-
return Vector3
-
# main script
-
if len(argv) < 8:
raise StandardError, \
"Syntax: dihedral_distribute.py dihedrals/impropers datafile nbin theta_min theta_max outfile files ..."
-
dt = data(argv[2])
nbins = int(argv[3])
theta_min = float(argv[4])
theta_max = float(argv[5])
outfile = argv[6]
files = ' '.join(argv[7:])
-
# get the angles from the data file
if == "dihedrals":
angle = dt.get("Dihedrals")
else if == "impropers":
angle = dt.get("Impropers")
else:
raise StandardError, "The second keyword is neither 'dihedrals' nor 'impropers' "
sys.exit()
-
nangles = len(angle)
atype = nangles * [0]
iatom = nangles * [0]
jatom = nangles * [0]
katom = nangles * [0]
latom = nangles * [0]
for i in xrange(nangles):
atype[i] = int(angle[i][1] - 1)
iatom[i] = int(angle[i][2] - 1)
jatom[i] = int(angle[i][3] - 1)
katom[i] = int(angle[i][4] - 1)
latom[i] = int(angle[i][5] - 1)
-
ntypes = 0
for i in xrange(nangles): ntypes = max(angle[i][1],ntypes)
ntypes = int(ntypes)
ncount = ntypes * [0]
bin = nbins * [0]
for i in xrange(nbins):
bin[i] = ntypes * [0]
-
# read snapshots one-at-a-time
-
d = dump(files,0)
d.map(1,"id",2,"type",3,"x",4,"y",5,"z")
-
while 1:
time = d.next()
if time == -1: break
box = (d.snaps[-1].xlo,d.snaps[-1].ylo,d.snaps[-1].zlo,
d.snaps[-1].xhi,d.snaps[-1].yhi,d.snaps[-1].zhi)
xprd = box[3] - box[0]
yprd = box[4] - box[1]
zprd = box[5] - box[2]
d.unscale()
d.sort()
x,y,z = d.vecs(time,"x","y","z")
-
VectorA = [0] * 3
VectorB = [0] * 3
VectorC = [0] * 3
for i in xrange(nangles):
# first vector
VectorA[0] = PBC(x[jatom[i]] - x[iatom[i]],xprd)
VectorA[1] = PBC(y[jatom[i]] - y[iatom[i]],yprd)
VectorA[1] = PBC(z[jatom[i]] - z[iatom[i]],zprd)
-
# second vector
VectorB[0] = PBC(x[katom[i]] - x[jatom[i]],xprd)
VectorB[1] = PBC(y[katom[i]] - y[jatom[i]],yprd)
VectorB[2] = PBC(z[katom[i]] - z[jatom[i]],zprd)
normVectorB = math.sqrt(InnerProduct(VectorB,VectorB))
-
# third vector
VectorC[0] = PBC(x[latom[i]] - x[katom[i]],xprd)
VectorC[1] = PBC(y[latom[i]] - y[katom[i]],yprd)
VectorC[2] = PBC(z[latom[i]] - z[katom[i]],zprd)
# compute the signed dihedral angle
HelpVector = CrossProduct(CrossProduct(VectorA,VectorB),CrossProduct(VectorB,VectorC))
Argument1 = InnerProduct(HelpVector,VectorB/normVectorB)
Argument2 = InnerProduct(CrossProduct(VectorA,VectorB),CrossProduct(VectorB,VectorC))
theta = 180.0*math.atan2(Argument1,Argument2)/math.pi
ibin = int(nbins*(theta - theta_min)/(theta_max - theta_min) + 0.5)
if ((ibin >= 0) and (ibin <= nbins-1)):
bin[ibin][atype[i]] += nbins
ncount[atype[i]] += 1
else:
print "Warning: dihedral outside specified range"
print "dihedral type:", atype[i]+1
print "dihedral number:", i
print time,
print
print "Printing normalized dihedral distributions to",outfile
fp = open(outfile,"w")
theta_range = theta_max - theta_min
for i in xrange(nbins):
print >>fp, theta_min + theta_range*float(i)/float(nbins),
for j in xrange(ntypes):
if (ncount[j] > 0):
print >>fp, float(bin[i][j])/float(ncount[j])/theta_range,
else:
print >>fp, 0.0,
print >>fp
fp.close()
- - - -
- -
- - - - - -
-
-
- -
-
- -
- - - - - - -
- - - You can’t perform that action at this time. -
- - - - - - - - - - -
- - You signed in with another tab or window. Reload to refresh your session. - You signed out in another tab or window. Reload to refresh your session. -
- - - - -
- Press h to open a hovercard with more details. -
- - - - - From f7b94b115ae4a753c9b0fd67e543496319f1069e Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Wed, 2 May 2018 19:48:29 +0200 Subject: [PATCH 03/23] script to compute the dihedral distribution compute the distribution of dihedral angles similar to bond and angle scripts. The dihedral angle is signed i.e. ranging from -180 to 180 degrees --- scripts/dihedral_distribute.py | 159 +++++++++++++++++++++++++++++++++ 1 file changed, 159 insertions(+) create mode 100644 scripts/dihedral_distribute.py diff --git a/scripts/dihedral_distribute.py b/scripts/dihedral_distribute.py new file mode 100644 index 0000000..0fee86d --- /dev/null +++ b/scripts/dihedral_distribute.py @@ -0,0 +1,159 @@ +#!/usr/bin/python + +# Script: dihedral_distribute.py +# Purpose: binned signed dihedral distributions by dihedral type +# Syntax: dihedral_distribute.py datafile nbin theta_min theta_max outfile files ... +# datafile = lammps data file +# nbin = # of bins per dihedral type +# theta_min = min expected angle +# theta_max = max expected angle length +# outfile = file to write stats to +# files = series of dump files +# Example: dihedral_distribute.py pore.data 1000 110. 120. angles.out pore.dump.1 +# Author: Evangelos Voyiatzis (TU Darmstadt) + +# enable script to run from Python directly w/out Pizza.py + +import sys +import math +from dump import dump + +if not globals().has_key("argv"): argv = sys.argv + +def PBC(distance,boxlength): # function to compute the minimum image + + if abs(distance) > 0.5*xprd: + if distance < 0.0: + distance += boxlength + else: + distance -= boxlength + + return distance + +def InnerProduct(Vector1,Vector2): # function to compute the inner product of two vectors + + return (Vector1[0]*Vector2[0] + Vector1[1]*Vector2[1] + Vector1[2]*Vector2[2]) + +def CrossProduct(Vector1,Vector2): # function to compute the Cross product of two vectors + + Vector3 = [0] *3 + Vector3[0] = Vector1[1]*Vector2[2] - Vector1[2]*Vector2[1] + Vector3[1] = Vector1[2]*Vector2[0] - Vector1[0]*Vector2[2] + Vector3[2] = Vector1[0]*Vector2[1] - Vector1[1]*Vector2[0] + + return Vector3 + +# main script + +if len(argv) < 8: + raise StandardError, \ + "Syntax: dihedral_distribute.py dihedrals/impropers datafile nbin theta_min theta_max outfile files ..." + +dt = data(argv[2]) +nbins = int(argv[3]) +theta_min = float(argv[4]) +theta_max = float(argv[5]) +outfile = argv[6] +files = ' '.join(argv[7:]) + +# get the angles from the data file +if == "dihedrals": + angle = dt.get("Dihedrals") +else if == "impropers": + angle = dt.get("Impropers") +else: + raise StandardError, "The second keyword is neither 'dihedrals' nor 'impropers' " + sys.exit() + +nangles = len(angle) +atype = nangles * [0] +iatom = nangles * [0] +jatom = nangles * [0] +katom = nangles * [0] +latom = nangles * [0] +for i in xrange(nangles): + atype[i] = int(angle[i][1] - 1) + iatom[i] = int(angle[i][2] - 1) + jatom[i] = int(angle[i][3] - 1) + katom[i] = int(angle[i][4] - 1) + latom[i] = int(angle[i][5] - 1) + +ntypes = 0 +for i in xrange(nangles): ntypes = max(angle[i][1],ntypes) +ntypes = int(ntypes) +ncount = ntypes * [0] +bin = nbins * [0] +for i in xrange(nbins): + bin[i] = ntypes * [0] + +# read snapshots one-at-a-time + +d = dump(files,0) +d.map(1,"id",2,"type",3,"x",4,"y",5,"z") + +while 1: + time = d.next() + if time == -1: break + + box = (d.snaps[-1].xlo,d.snaps[-1].ylo,d.snaps[-1].zlo, + d.snaps[-1].xhi,d.snaps[-1].yhi,d.snaps[-1].zhi) + + xprd = box[3] - box[0] + yprd = box[4] - box[1] + zprd = box[5] - box[2] + + d.unscale() + d.sort() + x,y,z = d.vecs(time,"x","y","z") + + VectorA = [0] * 3 + VectorB = [0] * 3 + VectorC = [0] * 3 + + for i in xrange(nangles): + # first vector + VectorA[0] = PBC(x[jatom[i]] - x[iatom[i]],xprd) + VectorA[1] = PBC(y[jatom[i]] - y[iatom[i]],yprd) + VectorA[1] = PBC(z[jatom[i]] - z[iatom[i]],zprd) + + # second vector + VectorB[0] = PBC(x[katom[i]] - x[jatom[i]],xprd) + VectorB[1] = PBC(y[katom[i]] - y[jatom[i]],yprd) + VectorB[2] = PBC(z[katom[i]] - z[jatom[i]],zprd) + normVectorB = math.sqrt(InnerProduct(VectorB,VectorB)) + + # third vector + VectorC[0] = PBC(x[latom[i]] - x[katom[i]],xprd) + VectorC[1] = PBC(y[latom[i]] - y[katom[i]],yprd) + VectorC[2] = PBC(z[latom[i]] - z[katom[i]],zprd) + + # compute the signed dihedral angle + HelpVector = CrossProduct(CrossProduct(VectorA,VectorB),CrossProduct(VectorB,VectorC)) + Argument1 = InnerProduct(HelpVector,VectorB/normVectorB) + Argument2 = InnerProduct(CrossProduct(VectorA,VectorB),CrossProduct(VectorB,VectorC)) + theta = 180.0*math.atan2(Argument1,Argument2)/math.pi + + ibin = int(nbins*(theta - theta_min)/(theta_max - theta_min) + 0.5) + if ((ibin >= 0) and (ibin <= nbins-1)): + bin[ibin][atype[i]] += nbins + ncount[atype[i]] += 1 + else: + print "Warning: dihedral outside specified range" + print "dihedral type:", atype[i]+1 + print "dihedral number:", i + print time, + +print +print "Printing normalized dihedral distributions to",outfile + +fp = open(outfile,"w") +theta_range = theta_max - theta_min +for i in xrange(nbins): + print >>fp, theta_min + theta_range*float(i)/float(nbins), + for j in xrange(ntypes): + if (ncount[j] > 0): + print >>fp, float(bin[i][j])/float(ncount[j])/theta_range, + else: + print >>fp, 0.0, + print >>fp +fp.close() From 1d252bbaefe170794d878d14dd6ba315662e52cb Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Wed, 9 May 2018 10:47:44 +0200 Subject: [PATCH 04/23] fixing missing import statement the data.py module is used without being imported. the added line of code is just the import command. --- scripts/angle_distribute.py | 1 + 1 file changed, 1 insertion(+) diff --git a/scripts/angle_distribute.py b/scripts/angle_distribute.py index 1f03398..e1988af 100644 --- a/scripts/angle_distribute.py +++ b/scripts/angle_distribute.py @@ -16,6 +16,7 @@ import sys from dump import dump +from data import data from math import sqrt,acos,atan if not globals().has_key("argv"): argv = sys.argv From 6e74b31ea132cd543ccba22e438eb21e43f174d7 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Wed, 9 May 2018 10:48:42 +0200 Subject: [PATCH 05/23] fixing missing import statement the data.py module is used without being imported. The added line takes care of that --- scripts/bond_distribute.py | 1 + 1 file changed, 1 insertion(+) diff --git a/scripts/bond_distribute.py b/scripts/bond_distribute.py index de2a600..a4c78ab 100644 --- a/scripts/bond_distribute.py +++ b/scripts/bond_distribute.py @@ -16,6 +16,7 @@ import sys from dump import dump +from data import data from math import sqrt if not globals().has_key("argv"): argv = sys.argv From 4b2e286deb4f1b2880449dd13ec7b68ac26df1ec Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Wed, 9 May 2018 10:49:46 +0200 Subject: [PATCH 06/23] fixing missing import statement The data.py module is used without being imported. The added line of code takes care of the import statement --- scripts/dihedral_distribute.py | 1 + 1 file changed, 1 insertion(+) diff --git a/scripts/dihedral_distribute.py b/scripts/dihedral_distribute.py index 0fee86d..21e0d57 100644 --- a/scripts/dihedral_distribute.py +++ b/scripts/dihedral_distribute.py @@ -16,6 +16,7 @@ import sys import math +from data import data from dump import dump if not globals().has_key("argv"): argv = sys.argv From ca0eac12ce00136a02c4c51456c0ffa5716f9ec7 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Wed, 9 May 2018 11:20:19 +0200 Subject: [PATCH 07/23] fixing bug in PBC function The variable xprd (the length in the x-direction) was replaced with boxlength which is the length of the simulation box in the desired cartesian coordinate --- scripts/dihedral_distribute.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/dihedral_distribute.py b/scripts/dihedral_distribute.py index 21e0d57..7461940 100644 --- a/scripts/dihedral_distribute.py +++ b/scripts/dihedral_distribute.py @@ -23,7 +23,7 @@ def PBC(distance,boxlength): # function to compute the minimum image - if abs(distance) > 0.5*xprd: + if abs(distance) > 0.5*boxlength: if distance < 0.0: distance += boxlength else: From 959accd5feac2ab31ec3549292ab498a6fe14745 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Wed, 9 May 2018 11:27:24 +0200 Subject: [PATCH 08/23] Update dihedral_distribute.py --- scripts/dihedral_distribute.py | 8 +------- 1 file changed, 1 insertion(+), 7 deletions(-) diff --git a/scripts/dihedral_distribute.py b/scripts/dihedral_distribute.py index 7461940..434f991 100644 --- a/scripts/dihedral_distribute.py +++ b/scripts/dihedral_distribute.py @@ -58,13 +58,7 @@ def CrossProduct(Vector1,Vector2): # function to compute the Cross product of tw files = ' '.join(argv[7:]) # get the angles from the data file -if == "dihedrals": - angle = dt.get("Dihedrals") -else if == "impropers": - angle = dt.get("Impropers") -else: - raise StandardError, "The second keyword is neither 'dihedrals' nor 'impropers' " - sys.exit() +angle = dt.get("Dihedrals") nangles = len(angle) atype = nangles * [0] From 8d955eff54e43d43e41449bd58a480bdf4cb4268 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Wed, 9 May 2018 11:28:32 +0200 Subject: [PATCH 09/23] Update dihedral_distribute.py --- scripts/dihedral_distribute.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/scripts/dihedral_distribute.py b/scripts/dihedral_distribute.py index 434f991..81311c2 100644 --- a/scripts/dihedral_distribute.py +++ b/scripts/dihedral_distribute.py @@ -46,9 +46,9 @@ def CrossProduct(Vector1,Vector2): # function to compute the Cross product of tw # main script -if len(argv) < 8: +if len(argv) < 7: raise StandardError, \ - "Syntax: dihedral_distribute.py dihedrals/impropers datafile nbin theta_min theta_max outfile files ..." + "Syntax: dihedral_distribute.py datafile nbin theta_min theta_max outfile files ..." dt = data(argv[2]) nbins = int(argv[3]) From 7f53913d66b6c0d32c26922c76f12bd0f62d9617 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Wed, 9 May 2018 11:36:06 +0200 Subject: [PATCH 10/23] Update dihedral_distribute.py --- scripts/dihedral_distribute.py | 18 ++++++++++++++++-- 1 file changed, 16 insertions(+), 2 deletions(-) diff --git a/scripts/dihedral_distribute.py b/scripts/dihedral_distribute.py index 81311c2..ff84fe6 100644 --- a/scripts/dihedral_distribute.py +++ b/scripts/dihedral_distribute.py @@ -46,10 +46,19 @@ def CrossProduct(Vector1,Vector2): # function to compute the Cross product of tw # main script -if len(argv) < 7: +if len(argv) < 8: raise StandardError, \ "Syntax: dihedral_distribute.py datafile nbin theta_min theta_max outfile files ..." +# get the angles from the data file +if argv[1] == "dihedrals": + angle = dt.get("Dihedrals") +elif argv[1] == "impropers": + angle = dt.get("Impropers") +else: + raise StandardError, "The second keyword is neither 'dihedrals' nor 'impropers' " + sys.exit() + dt = data(argv[2]) nbins = int(argv[3]) theta_min = float(argv[4]) @@ -104,6 +113,7 @@ def CrossProduct(Vector1,Vector2): # function to compute the Cross product of tw VectorA = [0] * 3 VectorB = [0] * 3 VectorC = [0] * 3 + VectorD = [0] * 3 for i in xrange(nangles): # first vector @@ -116,6 +126,10 @@ def CrossProduct(Vector1,Vector2): # function to compute the Cross product of tw VectorB[1] = PBC(y[katom[i]] - y[jatom[i]],yprd) VectorB[2] = PBC(z[katom[i]] - z[jatom[i]],zprd) normVectorB = math.sqrt(InnerProduct(VectorB,VectorB)) + + VectorD[0] = VectorB[0]/normVectorB + VectorD[1] = VectorB[1]/normVectorB + VectorD[2] = VectorB[2]/normVectorB # third vector VectorC[0] = PBC(x[latom[i]] - x[katom[i]],xprd) @@ -124,7 +138,7 @@ def CrossProduct(Vector1,Vector2): # function to compute the Cross product of tw # compute the signed dihedral angle HelpVector = CrossProduct(CrossProduct(VectorA,VectorB),CrossProduct(VectorB,VectorC)) - Argument1 = InnerProduct(HelpVector,VectorB/normVectorB) + Argument1 = InnerProduct(HelpVector,VectorD) Argument2 = InnerProduct(CrossProduct(VectorA,VectorB),CrossProduct(VectorB,VectorC)) theta = 180.0*math.atan2(Argument1,Argument2)/math.pi From 48ad5db7e6bc3148580a1be55add42187b3714d1 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Wed, 9 May 2018 11:42:02 +0200 Subject: [PATCH 11/23] Update dihedral_distribute.py --- scripts/dihedral_distribute.py | 29 +++++++++++++++-------------- 1 file changed, 15 insertions(+), 14 deletions(-) diff --git a/scripts/dihedral_distribute.py b/scripts/dihedral_distribute.py index ff84fe6..1ef25de 100644 --- a/scripts/dihedral_distribute.py +++ b/scripts/dihedral_distribute.py @@ -2,7 +2,8 @@ # Script: dihedral_distribute.py # Purpose: binned signed dihedral distributions by dihedral type -# Syntax: dihedral_distribute.py datafile nbin theta_min theta_max outfile files ... +# Syntax: dihedral_distribute.py dihedral/improper datafile nbin theta_min theta_max outfile files ... +# dihedral/improper = either dihedral or improper keyword to compute the respective distribution # datafile = lammps data file # nbin = # of bins per dihedral type # theta_min = min expected angle @@ -49,16 +50,7 @@ def CrossProduct(Vector1,Vector2): # function to compute the Cross product of tw if len(argv) < 8: raise StandardError, \ "Syntax: dihedral_distribute.py datafile nbin theta_min theta_max outfile files ..." - -# get the angles from the data file -if argv[1] == "dihedrals": - angle = dt.get("Dihedrals") -elif argv[1] == "impropers": - angle = dt.get("Impropers") -else: - raise StandardError, "The second keyword is neither 'dihedrals' nor 'impropers' " - sys.exit() - + dt = data(argv[2]) nbins = int(argv[3]) theta_min = float(argv[4]) @@ -67,7 +59,13 @@ def CrossProduct(Vector1,Vector2): # function to compute the Cross product of tw files = ' '.join(argv[7:]) # get the angles from the data file -angle = dt.get("Dihedrals") +if argv[1] == "dihedrals": + angle = dt.get("Dihedrals") +elif argv[1] == "impropers": + angle = dt.get("Impropers") +else: + raise StandardError, "The second keyword is neither 'dihedrals' nor 'impropers' " + sys.exit() nangles = len(angle) atype = nangles * [0] @@ -153,8 +151,11 @@ def CrossProduct(Vector1,Vector2): # function to compute the Cross product of tw print time, print -print "Printing normalized dihedral distributions to",outfile - +if argv[1] == "dihedrals": + print "Printing normalized dihedral distributions to",outfile +else: + print "Printing normalized improper distributions to",outfile + fp = open(outfile,"w") theta_range = theta_max - theta_min for i in xrange(nbins): From 73086bd9ccf8b5fe27294758aeadfcf18f5de429 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Wed, 9 May 2018 11:49:09 +0200 Subject: [PATCH 12/23] Update dihedral_distribute.py --- scripts/dihedral_distribute.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/scripts/dihedral_distribute.py b/scripts/dihedral_distribute.py index 1ef25de..fce94ab 100644 --- a/scripts/dihedral_distribute.py +++ b/scripts/dihedral_distribute.py @@ -3,7 +3,7 @@ # Script: dihedral_distribute.py # Purpose: binned signed dihedral distributions by dihedral type # Syntax: dihedral_distribute.py dihedral/improper datafile nbin theta_min theta_max outfile files ... -# dihedral/improper = either dihedral or improper keyword to compute the respective distribution +# dihedrals/impropers = either dihedral or improper keyword to compute the respective distribution # datafile = lammps data file # nbin = # of bins per dihedral type # theta_min = min expected angle @@ -49,7 +49,7 @@ def CrossProduct(Vector1,Vector2): # function to compute the Cross product of tw if len(argv) < 8: raise StandardError, \ - "Syntax: dihedral_distribute.py datafile nbin theta_min theta_max outfile files ..." + "Syntax: dihedral_distribute.py dihedrals/impropers datafile nbin theta_min theta_max outfile files ..." dt = data(argv[2]) nbins = int(argv[3]) From 9c4e64451418f7879db72632e79e0867ae9d371e Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Wed, 9 May 2018 15:00:47 +0200 Subject: [PATCH 13/23] Update bond_distribute.py replace the open and close statements of the output file with a "with" construct --- scripts/bond_distribute.py | 21 ++++++++++----------- 1 file changed, 10 insertions(+), 11 deletions(-) diff --git a/scripts/bond_distribute.py b/scripts/bond_distribute.py index a4c78ab..17cbc16 100644 --- a/scripts/bond_distribute.py +++ b/scripts/bond_distribute.py @@ -110,14 +110,13 @@ print print "Printing bond distance normalized distribution to",outfile -fp = open(outfile,"w") -rrange = rmax - rmin -for i in xrange(nbins): - print >>fp, rmin + rrange*float(i)/float(nbins), - for j in xrange(ntypes): - if (ncount[j] > 0): - print >>fp, float(bin[i][j])/float(ncount[j])/rrange, - else: - print >>fp, 0.0, - print >>fp -fp.close() +with open(outfile,"w") as fp: + rrange = rmax - rmin + for i in xrange(nbins): + print >>fp, rmin + rrange*float(i)/float(nbins), + for j in xrange(ntypes): + if (ncount[j] > 0): + print >>fp, float(bin[i][j])/float(ncount[j])/rrange, + else: + print >>fp, 0.0, + print >>fp From 324367dc9bae317dcd647ad6bfb4a75496c1fcc3 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Wed, 9 May 2018 15:03:27 +0200 Subject: [PATCH 14/23] Update bond_distribute.py removing the box variable since it was not really used in the code --- scripts/bond_distribute.py | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/scripts/bond_distribute.py b/scripts/bond_distribute.py index 17cbc16..72eaa89 100644 --- a/scripts/bond_distribute.py +++ b/scripts/bond_distribute.py @@ -62,12 +62,9 @@ time = d.next() if time == -1: break - box = (d.snaps[-1].xlo,d.snaps[-1].ylo,d.snaps[-1].zlo, - d.snaps[-1].xhi,d.snaps[-1].yhi,d.snaps[-1].zhi) - - xprd = box[3] - box[0] - yprd = box[4] - box[1] - zprd = box[5] - box[2] + xprd = d.snaps[-1].xhi - d.snaps[-1].xlo + yprd = d.snaps[-1].yhi - d.snaps[-1].ylo + zprd = d.snaps[-1].zhi - d.snaps[-1].zlo d.unscale() d.sort() From 71332497a94416627f963661114228168e966d29 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Wed, 9 May 2018 15:11:43 +0200 Subject: [PATCH 15/23] Update angle_distribute.py remove the open and close statements of the output file and replace it with an "with" construct modified one printing statement to use the format option --- scripts/angle_distribute.py | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/scripts/angle_distribute.py b/scripts/angle_distribute.py index e1988af..81bfc1b 100644 --- a/scripts/angle_distribute.py +++ b/scripts/angle_distribute.py @@ -142,16 +142,16 @@ print time, print -print "Printing normalized angle distributions to",outfile +print("Printing normalized angle distributions to {}".format(outfile)) -fp = open(outfile,"w") -theta_range = theta_max - theta_min -for i in xrange(nbins): - print >>fp, theta_min + theta_range*float(i)/float(nbins), - for j in xrange(ntypes): - if (ncount[j] > 0): - print >>fp, float(bin[i][j])/float(ncount[j])/theta_range, - else: - print >>fp, 0.0, - print >>fp -fp.close() +with open(outfile,"w") as fp: + theta_range = theta_max - theta_min + for i in xrange(nbins): + print >>fp, theta_min + theta_range*float(i)/float(nbins), + for j in xrange(ntypes): + if (ncount[j] > 0): + print >>fp, float(bin[i][j])/float(ncount[j])/theta_range, + else: + print >>fp, 0.0, + print >>fp + From 68667d1af00b86dceb150b86dc86ee0ef3123e93 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Wed, 9 May 2018 21:09:37 +0200 Subject: [PATCH 16/23] Update bond_distribute.py --- scripts/bond_distribute.py | 147 ++++++++++++++++++------------------- 1 file changed, 72 insertions(+), 75 deletions(-) diff --git a/scripts/bond_distribute.py b/scripts/bond_distribute.py index 72eaa89..5ef7c26 100644 --- a/scripts/bond_distribute.py +++ b/scripts/bond_distribute.py @@ -14,100 +14,94 @@ # enable script to run from Python directly w/out Pizza.py -import sys +from argparse import ArgumentParser +import numpy from dump import dump from data import data from math import sqrt -if not globals().has_key("argv"): argv = sys.argv # main script +def compute_bond_distribution(): + parser = ArgumentParser(description='A python script to bond distribution') -if len(argv) < 7: - raise StandardError, \ - "Syntax: bond_distribute.py datafile nbin rmin rmax outfile files ..." + parser.add_argument("input_data_file", help="The name of the lammps input data file") + parser.add_argument("nbins", type=int, help="# of bins per bond type") + parser.add_argument("rmin", type=float, help="min expected bond length") + parser.add_argument("rmax", type=float, help="max expected bond length") + parser.add_argument("output_file", help="The name of the file to write stats") + parser.add_argument("dump_files", nargs='+', help="series of dump files") -dt = data(argv[1]) -nbins = int(argv[2]) -rmin = float(argv[3]) -rmax = float(argv[4]) -outfile = argv[5] -files = ' '.join(argv[6:]) + args = parser.parse_args() -# get the bonds from the data file + dt = data(args.input_data_file) + nbins = args.nbins + rmin = args.rmin + rmax = args.rmax + files = ' '.join(args.dump_files[:]) -bond = dt.get("Bonds") -nbonds = len(bond) -btype = nbonds * [0] -iatom = nbonds * [0] -jatom = nbonds * [0] -for i in xrange(nbonds): - btype[i] = int(bond[i][1] - 1) - iatom[i] = int(bond[i][2] - 1) - jatom[i] = int(bond[i][3] - 1) + # get the bonds from the data file -ntypes = 0 -for i in xrange(nbonds): ntypes = max(bond[i][1],ntypes) -ntypes = int(ntypes) -ncount = ntypes * [0] -bin = nbins * [0] -for i in xrange(nbins): - bin[i] = ntypes * [0] + bond = dt.get("Bonds") + nbonds = len(bond) + btype = numpy.asarray(bond[:][1], dtype=int) - 1 + iatom = numpy.asarray(bond[:][2], dtype=int) - 1 + jatom = numpy.asarray(bond[:][3], dtype=int) - 1 -# read snapshots one-at-a-time + ntypes = 0 + for i in xrange(nbonds): ntypes = max(bond[i][1],ntypes) + ntypes = int(ntypes) + ncount = numpy.zeros(ntypes, dtype=int) + bin = nbins * [0] + for i in xrange(nbins): + bin[i] = ntypes * [0] -d = dump(files,0) -d.map(1,"id",2,"type",3,"x",4,"y",5,"z") + # read snapshots one-at-a-time + d = dump(files,0) + d.map(1,"id",2,"type",3,"x",4,"y",5,"z") -while 1: - time = d.next() - if time == -1: break + while True: + time = d.next() + if time == -1: + break - xprd = d.snaps[-1].xhi - d.snaps[-1].xlo - yprd = d.snaps[-1].yhi - d.snaps[-1].ylo - zprd = d.snaps[-1].zhi - d.snaps[-1].zlo + xprd = d.snaps[-1].xhi - d.snaps[-1].xlo + yprd = d.snaps[-1].yhi - d.snaps[-1].ylo + zprd = d.snaps[-1].zhi - d.snaps[-1].zlo - d.unscale() - d.sort() - x,y,z = d.vecs(time,"x","y","z") - - for i in xrange(nbonds): + d.unscale() + d.sort() + x,y,z = d.vecs(time,"x","y","z") - delx = x[jatom[i]] - x[iatom[i]] - dely = y[jatom[i]] - y[iatom[i]] - delz = z[jatom[i]] - z[iatom[i]] - - if abs(delx) > 0.5*xprd: - if delx < 0.0: - delx += xprd - else: - delx -= xprd - if abs(dely) > 0.5*yprd: - if dely < 0.0: - dely += yprd - else: - dely -= yprd - if abs(delz) > 0.5*zprd: - if delz < 0.0: - delz += zprd - else: - delz -= zprd - - r = sqrt(delx*delx + dely*dely + delz*delz) + x = numpy.asarray(x) + y = numpy.asarray(y) + z = numpy.asarray(z) + + delx = x[jatom[:]] - x[iatom[:]] + dely = y[jatom[:]] - y[iatom[:]] + delz = z[jatom[:]] - z[iatom[:]] + + delx -= xprd*numpy.rint((x[jatom[:]] - x[iatom[:]])/xprd) + dely -= yprd*numpy.rint((y[jatom[:]] - y[iatom[:]])/yprd) + delz -= zprd*numpy.rint((z[jatom[:]] - z[iatom[:]])/zprd) + + r = numpy.sqrt(delx*delx + dely*dely + delz*delz) + + for i in xrange(nbonds): - ibin = int(nbins*(r - rmin)/(rmax - rmin) + 0.5) - if ((ibin >= 0) and (ibin <= nbins-1)): - bin[ibin][btype[i]] += nbins - ncount[btype[i]] += 1 - else: - print "Warning: bond distance outside specified range" - print "Bond type:", btype[i]+1 - print "Bond number:", i - print time, + ibin = int(nbins*(r[i] - rmin)/(rmax - rmin) + 0.5) + if ((ibin >= 0) and (ibin <= nbins-1)): + bin[ibin][btype[i]] += nbins + ncount[btype[i]] += 1 + else: + print "Warning: bond distance outside specified range" + print "Bond type:", btype[i]+1 + print "Bond number:", i + print time, -print -print "Printing bond distance normalized distribution to",outfile + print + print("Printing bond distance normalized distribution to {}".format(args.output_file)) -with open(outfile,"w") as fp: + with open(args.output_file,"w") as fp: rrange = rmax - rmin for i in xrange(nbins): print >>fp, rmin + rrange*float(i)/float(nbins), @@ -117,3 +111,6 @@ else: print >>fp, 0.0, print >>fp + +if __name__ == "__main__": + compute_bond_distribution() From 410919eec84a40ba775a1e31b358c67aa5ac1314 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Wed, 9 May 2018 21:59:01 +0200 Subject: [PATCH 17/23] Update bond_distribute.py large modifications to the code using numpy for the calculations - the goal is to speed up the computation --- scripts/bond_distribute.py | 77 +++++++++++++++++--------------------- 1 file changed, 35 insertions(+), 42 deletions(-) diff --git a/scripts/bond_distribute.py b/scripts/bond_distribute.py index 5ef7c26..435f42f 100644 --- a/scripts/bond_distribute.py +++ b/scripts/bond_distribute.py @@ -1,28 +1,30 @@ #!/usr/bin/python +""" +Script: bond_distribute.py +Purpose: binned bond length distributions by bond type +Syntax: bond_distribute.py datafile nbin rmin rmax outfile files ... + datafile = lammps data file + nbin = # of bins per bond type + rmin = min expected bond length + rmax = max expected bond length + outfile = file to write stats to + files = series of dump files +Example: bond_distribute.py pore.data 1000 1.5 1.85 bonds.out pore.dump.1 +Author: Paul Crozier (Sandia) -# Script: bond_distribute.py -# Purpose: binned bond length distributions by bond type -# Syntax: bond_distribute.py datafile nbin rmin rmax outfile files ... -# datafile = lammps data file -# nbin = # of bins per bond type -# rmin = min expected bond length -# rmax = max expected bond length -# outfile = file to write stats to -# files = series of dump files -# Example: bond_distribute.py pore.data 1000 1.5 1.85 bonds.out pore.dump.1 -# Author: Paul Crozier (Sandia) - -# enable script to run from Python directly w/out Pizza.py +enable script to run from Python directly w/out Pizza.py +""" from argparse import ArgumentParser import numpy from dump import dump from data import data -from math import sqrt # main script def compute_bond_distribution(): - parser = ArgumentParser(description='A python script to bond distribution') + """The function doing the actual calculation of the bond distribution""" + + parser = ArgumentParser(description='A python script to compute bond distribution using pizza.py') parser.add_argument("input_data_file", help="The name of the lammps input data file") parser.add_argument("nbins", type=int, help="# of bins per bond type") @@ -34,35 +36,33 @@ def compute_bond_distribution(): args = parser.parse_args() dt = data(args.input_data_file) - nbins = args.nbins - rmin = args.rmin - rmax = args.rmax files = ' '.join(args.dump_files[:]) # get the bonds from the data file - bond = dt.get("Bonds") nbonds = len(bond) btype = numpy.asarray(bond[:][1], dtype=int) - 1 iatom = numpy.asarray(bond[:][2], dtype=int) - 1 jatom = numpy.asarray(bond[:][3], dtype=int) - 1 - ntypes = 0 - for i in xrange(nbonds): ntypes = max(bond[i][1],ntypes) - ntypes = int(ntypes) + ntypes = int(numpy.max(btype)) + 1 + bin = numpy.zeros((args.nbins, ntypes), dtype=int) + ncount = numpy.zeros(ntypes, dtype=int) - bin = nbins * [0] - for i in xrange(nbins): - bin[i] = ntypes * [0] - + for itype in range(0, ntypes): + ncount[itype] = numpy.sum(btype==itype) + # read snapshots one-at-a-time d = dump(files,0) d.map(1,"id",2,"type",3,"x",4,"y",5,"z") + nconfs = 0 while True: time = d.next() if time == -1: break + + nconfs += 1 xprd = d.snaps[-1].xhi - d.snaps[-1].xlo yprd = d.snaps[-1].yhi - d.snaps[-1].ylo @@ -86,28 +86,21 @@ def compute_bond_distribution(): r = numpy.sqrt(delx*delx + dely*dely + delz*delz) - for i in xrange(nbonds): + for itype in range(0, ntypes): + xx, hist_edge = numpy.histogram(r[btype == itype],bins=args.nbins, range=(args.rmin, args.rmax)) + bin[:][itype] += xx - ibin = int(nbins*(r[i] - rmin)/(rmax - rmin) + 0.5) - if ((ibin >= 0) and (ibin <= nbins-1)): - bin[ibin][btype[i]] += nbins - ncount[btype[i]] += 1 - else: - print "Warning: bond distance outside specified range" - print "Bond type:", btype[i]+1 - print "Bond number:", i - print time, + print("{} ".format(time)) - print print("Printing bond distance normalized distribution to {}".format(args.output_file)) with open(args.output_file,"w") as fp: - rrange = rmax - rmin - for i in xrange(nbins): - print >>fp, rmin + rrange*float(i)/float(nbins), - for j in xrange(ntypes): + rrange = args.rmax - args.rmin + for i in range(0, args.nbins): + print >>fp, args.rmin + rrange*float(i)/float(args.nbins), + for j in range(0, ntypes): if (ncount[j] > 0): - print >>fp, float(bin[i][j])/float(ncount[j])/rrange, + print >>fp, float(bin[i][j])/float(ncount[j]*nconfs), else: print >>fp, 0.0, print >>fp From 239329a58ca14dc4f7cb172e68e2fd4f82ca119b Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Wed, 9 May 2018 22:03:41 +0200 Subject: [PATCH 18/23] Update bond_distribute.py --- scripts/bond_distribute.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/bond_distribute.py b/scripts/bond_distribute.py index 435f42f..9ddc095 100644 --- a/scripts/bond_distribute.py +++ b/scripts/bond_distribute.py @@ -100,7 +100,7 @@ def compute_bond_distribution(): print >>fp, args.rmin + rrange*float(i)/float(args.nbins), for j in range(0, ntypes): if (ncount[j] > 0): - print >>fp, float(bin[i][j])/float(ncount[j]*nconfs), + print >>fp, float(bin[i][j])/float(ncount[j]*nconfs)/rrange, else: print >>fp, 0.0, print >>fp From 889600ddb6a1e8e2a6220701d2473d46542d12c1 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Wed, 9 May 2018 22:11:07 +0200 Subject: [PATCH 19/23] Update bond_distribute.py --- scripts/bond_distribute.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/scripts/bond_distribute.py b/scripts/bond_distribute.py index 9ddc095..2db2cf4 100644 --- a/scripts/bond_distribute.py +++ b/scripts/bond_distribute.py @@ -94,16 +94,16 @@ def compute_bond_distribution(): print("Printing bond distance normalized distribution to {}".format(args.output_file)) - with open(args.output_file,"w") as fp: + with open(args.output_file,"w") as output_file: rrange = args.rmax - args.rmin for i in range(0, args.nbins): - print >>fp, args.rmin + rrange*float(i)/float(args.nbins), + output_file.write("{} ".format(args.rmin + rrange*float(i)/float(args.nbins))) for j in range(0, ntypes): if (ncount[j] > 0): - print >>fp, float(bin[i][j])/float(ncount[j]*nconfs)/rrange, + output_file.write("{} ".format(float(bin[i][j])/float(ncount[j]*nconfs)/rrange)) else: - print >>fp, 0.0, - print >>fp + output_file.write("{} ".format(0.0)) + output_file.write("\n") if __name__ == "__main__": compute_bond_distribution() From a203cb156951e9790ddcca75d621c172cfa342b0 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Wed, 9 May 2018 22:55:16 +0200 Subject: [PATCH 20/23] Update bond_distribute.py --- scripts/bond_distribute.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/scripts/bond_distribute.py b/scripts/bond_distribute.py index 2db2cf4..3cb7246 100644 --- a/scripts/bond_distribute.py +++ b/scripts/bond_distribute.py @@ -88,6 +88,9 @@ def compute_bond_distribution(): for itype in range(0, ntypes): xx, hist_edge = numpy.histogram(r[btype == itype],bins=args.nbins, range=(args.rmin, args.rmax)) + if numpy.sum(xx) != ncount[itype]: + print("Warning: bond distance outside specified range ") + print("Bond type: {}".format(itype+1)) bin[:][itype] += xx print("{} ".format(time)) From 4a56200dcbfc047a395282dde45b2eaa3b6ced4f Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Wed, 9 May 2018 22:58:02 +0200 Subject: [PATCH 21/23] Update bond_distribute.py --- scripts/bond_distribute.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/bond_distribute.py b/scripts/bond_distribute.py index 3cb7246..219cbb6 100644 --- a/scripts/bond_distribute.py +++ b/scripts/bond_distribute.py @@ -40,7 +40,7 @@ def compute_bond_distribution(): # get the bonds from the data file bond = dt.get("Bonds") - nbonds = len(bond) + btype = numpy.asarray(bond[:][1], dtype=int) - 1 iatom = numpy.asarray(bond[:][2], dtype=int) - 1 jatom = numpy.asarray(bond[:][3], dtype=int) - 1 From 7e6f11f3e39e51a139435f7f631d9dad6207bad8 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Sat, 12 May 2018 23:13:13 +0200 Subject: [PATCH 22/23] Update angle_distribute.py rewritting angle_distribute using numpy and in more pythonic style than c modifications are similar to those in bond_distribute.py script --- scripts/angle_distribute.py | 227 ++++++++++++++++-------------------- 1 file changed, 100 insertions(+), 127 deletions(-) diff --git a/scripts/angle_distribute.py b/scripts/angle_distribute.py index 81bfc1b..bde7031 100644 --- a/scripts/angle_distribute.py +++ b/scripts/angle_distribute.py @@ -1,157 +1,130 @@ #!/usr/bin/python +""" +Script: angle_distribute.py +Purpose: binned angle distributions by angle type +Syntax: angle_distribute.py datafile nbin theta_min theta_max outfile files ... + datafile = lammps data file + nbin = # of bins per angle type + theta_min = min expected angle + theta_max = max expected angle length + outfile = file to write stats to + files = series of dump files +Example: angle_distribute.py pore.data 1000 110. 120. angles.out pore.dump.1 +Author: Paul Crozier (Sandia) -# Script: angle_distribute.py -# Purpose: binned angle distributions by angle type -# Syntax: angle_distribute.py datafile nbin theta_min theta_max outfile files ... -# datafile = lammps data file -# nbin = # of bins per angle type -# theta_min = min expected angle -# theta_max = max expected angle length -# outfile = file to write stats to -# files = series of dump files -# Example: angle_distribute.py pore.data 1000 110. 120. angles.out pore.dump.1 -# Author: Paul Crozier (Sandia) +enable script to run from Python directly w/out Pizza.py +""" -# enable script to run from Python directly w/out Pizza.py - -import sys +from argparse import ArgumentParser +import numpy from dump import dump from data import data -from math import sqrt,acos,atan -if not globals().has_key("argv"): argv = sys.argv # main script +def compute_angle_distribution(): + """The function doing the actual calculation of the bond distribution""" + + parser = ArgumentParser(description='A python script to compute bond distribution using pizza.py') -if len(argv) < 7: - raise StandardError, \ - "Syntax: angle_distribute.py datafile nbin theta_min theta_max outfile files ..." - -dt = data(argv[1]) -nbins = int(argv[2]) -theta_min = float(argv[3]) -theta_max = float(argv[4]) -outfile = argv[5] -files = ' '.join(argv[6:]) - -# get the angles from the data file + parser.add_argument("input_data_file", help="The name of the lammps input data file") + parser.add_argument("nbins", type=int, help="# of bins per bond type") + parser.add_argument("theta_min", type=float, help="min expected angle") + parser.add_argument("theta_max", type=float, help="max expected angle") + parser.add_argument("output_file", help="The name of the file to write stats") + parser.add_argument("dump_files", nargs='+', help="series of dump files") -angle = dt.get("Angles") -nangles = len(angle) -atype = nangles * [0] -iatom = nangles * [0] -jatom = nangles * [0] -katom = nangles * [0] -for i in xrange(nangles): - atype[i] = int(angle[i][1] - 1) - iatom[i] = int(angle[i][2] - 1) - jatom[i] = int(angle[i][3] - 1) - katom[i] = int(angle[i][4] - 1) - -ntypes = 0 -for i in xrange(nangles): ntypes = max(angle[i][1],ntypes) -ntypes = int(ntypes) -ncount = ntypes * [0] -bin = nbins * [0] -for i in xrange(nbins): - bin[i] = ntypes * [0] + args = parser.parse_args() + + dt = data(args.input_data_file) + files = ' '.join(args.dump_files[:]) -# read snapshots one-at-a-time + # get the angles from the data file + angle = dt.get("Angles") + nangles = len(angle) + atype = numpy.asarray(angle[:][1], dtype=int) - 1 + iatom = numpy.asarray(angle[:][2], dtype=int) - 1 + jatom = numpy.asarray(angle[:][3], dtype=int) - 1 + katom = numpy.asarray(angle[:][4], dtype=int) - 1 + + ntypes = int(numpy.max(btype)) + 1 + bin = numpy.zeros((args.nbins, ntypes), dtype=int) + + ncount = numpy.zeros(ntypes, dtype=int) + for itype in range(0, ntypes): + ncount[itype] = numpy.sum(atype==itype) -d = dump(files,0) -d.map(1,"id",2,"type",3,"x",4,"y",5,"z") + # read snapshots one-at-a-time -PI = 4.0*atan(1.0) + d = dump(files,0) + d.map(1,"id",2,"type",3,"x",4,"y",5,"z") -while 1: + while True: time = d.next() - if time == -1: break - - box = (d.snaps[-1].xlo,d.snaps[-1].ylo,d.snaps[-1].zlo, - d.snaps[-1].xhi,d.snaps[-1].yhi,d.snaps[-1].zhi) - - xprd = box[3] - box[0] - yprd = box[4] - box[1] - zprd = box[5] - box[2] + if time == -1: + break + + xprd = d.snaps[-1].xhi - d.snaps[-1].xlo + yprd = d.snaps[-1].yhi - d.snaps[-1].ylo + zprd = d.snaps[-1].zhi - d.snaps[-1].zlo d.unscale() d.sort() x,y,z = d.vecs(time,"x","y","z") - - for i in xrange(nangles): - delx1 = x[iatom[i]] - x[jatom[i]] - dely1 = y[iatom[i]] - y[jatom[i]] - delz1 = z[iatom[i]] - z[jatom[i]] - - if abs(delx1) > 0.5*xprd: - if delx1 < 0.0: - delx1 += xprd - else: - delx1 -= xprd - if abs(dely1) > 0.5*yprd: - if dely1 < 0.0: - dely1 += yprd - else: - dely1 -= yprd - if abs(delz1) > 0.5*zprd: - if delz1 < 0.0: - delz1 += zprd - else: - delz1 -= zprd - - r1 = sqrt(delx1*delx1 + dely1*dely1 + delz1*delz1) + x = numpy.asarray(x) + y = numpy.asarray(y) + z = numpy.asarray(z) - delx2 = x[katom[i]] - x[jatom[i]] - dely2 = y[katom[i]] - y[jatom[i]] - delz2 = z[katom[i]] - z[jatom[i]] + delx = x[iatom[:]] - x[jatom[:]] + dely = y[iatom[:]] - y[jatom[:]] + delz = z[iatom[:]] - z[jatom[:]] + + delx -= xprd*numpy.rint((x[iatom[:]] - x[jatom[:]])/xprd) + dely -= yprd*numpy.rint((y[iatom[:]] - y[jatom[:]])/yprd) + delz -= zprd*numpy.rint((z[iatom[:]] - z[jatom[:]])/zprd) - if abs(delx2) > 0.5*xprd: - if delx2 < 0.0: - delx2 += xprd - else: - delx2 -= xprd - if abs(dely2) > 0.5*yprd: - if dely2 < 0.0: - dely2 += yprd - else: - dely2 -= yprd - if abs(delz2) > 0.5*zprd: - if delz2 < 0.0: - delz2 += zprd - else: - delz2 -= zprd + r1 = numpy.sqrt(delx*delx + dely*dely + delz*delz) + + delx2 = x[katom[:]] - x[jatom[:]] + dely2 = y[katom[:]] - y[jatom[:]] + delz2 = z[katom[:]] - z[jatom[:]] + + delx2 -= xprd*numpy.rint((x[katom[:]] - x[jatom[:]])/xprd) + dely2 -= yprd*numpy.rint((y[katom[:]] - y[jatom[:]])/yprd) + delz2 -= zprd*numpy.rint((z[katom[:]] - z[jatom[:]])/zprd) - r2 = sqrt(delx2*delx2 + dely2*dely2 + delz2*delz2) + r2 = numpy.sqrt(delx2*delx2 + dely2*dely2 + delz2*delz2) - c = delx1*delx2 + dely1*dely2 + delz1*delz2 - c /= r1*r2 + c = delx1*delx2 + dely1*dely2 + delz1*delz2 + c /= r1*r2 - if (c > 1.0): c = 1.0 - if (c < -1.0): c = -1.0 + if (c > 1.0): c = 1.0 + if (c < -1.0): c = -1.0 - theta = 180.0*acos(c)/PI - - ibin = int(nbins*(theta - theta_min)/(theta_max - theta_min) + 0.5) - if ((ibin >= 0) and (ibin <= nbins-1)): - bin[ibin][atype[i]] += nbins - ncount[atype[i]] += 1 - else: - print "Warning: angle outside specified range" - print "angle type:", atype[i]+1 - print "angle number:", i - print time, + theta = 180.0*numpy.arccos(c)/numpy.pi + + for itype in range(0, ntypes): + xx, hist_edge = numpy.histogram(theta[atype == itype],bins=args.nbins, range=(args.theta_min, args.theta_max)) + if numpy.sum(xx) != ncount[itype]: + print("Warning: angle distance outside specified range ") + print("Angle type: {}".format(itype+1)) + bin[:][itype] += xx + + print("{} ".format(time)) -print -print("Printing normalized angle distributions to {}".format(outfile)) + print + print("Printing normalized angle distributions to {}".format(args.output_file)) -with open(outfile,"w") as fp: - theta_range = theta_max - theta_min - for i in xrange(nbins): - print >>fp, theta_min + theta_range*float(i)/float(nbins), - for j in xrange(ntypes): + with open(args.output_file,"w") as output_file: + theta_range = args.theta_max - args.theta_min + for i in range(0, args.nbins): + output_file.write("{} ".format(theta_min + theta_range*float(i)/float(nbins))) + for j in range(0, ntypes): if (ncount[j] > 0): - print >>fp, float(bin[i][j])/float(ncount[j])/theta_range, + output_file.write("{} ".format(float(bin[i][j])/float(ncount[j]*nconfs)/theta_range)) else: - print >>fp, 0.0, - print >>fp + output_file.write("{} ".format(0.0)) + output_file.write("\n") +if __name__ == "__main__": + compute_angle_distribution() From fed1b3b9c574672e26295ae1744cd8195dc82722 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Sat, 12 May 2018 23:16:18 +0200 Subject: [PATCH 23/23] Update dihedral_distribute.py --- scripts/dihedral_distribute.py | 35 +++++++++++++++++----------------- 1 file changed, 18 insertions(+), 17 deletions(-) diff --git a/scripts/dihedral_distribute.py b/scripts/dihedral_distribute.py index fce94ab..7646cc3 100644 --- a/scripts/dihedral_distribute.py +++ b/scripts/dihedral_distribute.py @@ -1,19 +1,20 @@ #!/usr/bin/python - -# Script: dihedral_distribute.py -# Purpose: binned signed dihedral distributions by dihedral type -# Syntax: dihedral_distribute.py dihedral/improper datafile nbin theta_min theta_max outfile files ... -# dihedrals/impropers = either dihedral or improper keyword to compute the respective distribution -# datafile = lammps data file -# nbin = # of bins per dihedral type -# theta_min = min expected angle -# theta_max = max expected angle length -# outfile = file to write stats to -# files = series of dump files -# Example: dihedral_distribute.py pore.data 1000 110. 120. angles.out pore.dump.1 -# Author: Evangelos Voyiatzis (TU Darmstadt) - -# enable script to run from Python directly w/out Pizza.py +""" +Script: dihedral_distribute.py +Purpose: binned signed dihedral distributions by dihedral type +Syntax: dihedral_distribute.py dihedral/improper datafile nbin theta_min theta_max outfile files ... + dihedrals/impropers = either dihedral or improper keyword to compute the respective distribution + datafile = lammps data file + nbin = # of bins per dihedral type + theta_min = min expected angle + theta_max = max expected angle length + outfile = file to write stats to + files = series of dump files +Example: dihedral_distribute.py pore.data 1000 110. 120. angles.out pore.dump.1 +Author: Evangelos Voyiatzis (TU Darmstadt) + +enable script to run from Python directly w/out Pizza.py +""" import sys import math @@ -152,9 +153,9 @@ def CrossProduct(Vector1,Vector2): # function to compute the Cross product of tw print if argv[1] == "dihedrals": - print "Printing normalized dihedral distributions to",outfile + print("Printing normalized dihedral distributions to {}".format(outfile) else: - print "Printing normalized improper distributions to",outfile + print("Printing normalized improper distributions to {}".format(outfile) fp = open(outfile,"w") theta_range = theta_max - theta_min