diff --git a/scripts/angle_distribute.py b/scripts/angle_distribute.py index 1f03398..bde7031 100644 --- a/scripts/angle_distribute.py +++ b/scripts/angle_distribute.py @@ -1,156 +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 math import sqrt,acos,atan -if not globals().has_key("argv"): argv = sys.argv +from data import data # 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",outfile + print + print("Printing normalized angle distributions to {}".format(args.output_file)) -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(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): + output_file.write("{} ".format(float(bin[i][j])/float(ncount[j]*nconfs)/theta_range)) + else: + output_file.write("{} ".format(0.0)) + output_file.write("\n") + +if __name__ == "__main__": + compute_angle_distribution() diff --git a/scripts/bond_distribute.py b/scripts/bond_distribute.py index de2a600..219cbb6 100644 --- a/scripts/bond_distribute.py +++ b/scripts/bond_distribute.py @@ -1,122 +1,112 @@ #!/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 - -import sys +from argparse import ArgumentParser +import numpy from dump import dump -from math import sqrt -if not globals().has_key("argv"): argv = sys.argv +from data import data # main script +def compute_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') -if len(argv) < 7: - raise StandardError, \ - "Syntax: bond_distribute.py datafile nbin rmin rmax outfile files ..." - -dt = data(argv[1]) -nbins = int(argv[2]) -rmin = float(argv[3]) -rmax = float(argv[4]) -outfile = argv[5] -files = ' '.join(argv[6:]) + 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") -# get the bonds from the data file + args = parser.parse_args() -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) + dt = data(args.input_data_file) + files = ' '.join(args.dump_files[:]) -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] + # get the bonds from the data file + bond = dt.get("Bonds") + + 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 = 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(btype==itype) + + # read snapshots one-at-a-time + d = dump(files,0) + d.map(1,"id",2,"type",3,"x",4,"y",5,"z") -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 + nconfs = 0 + while True: + time = d.next() + if time == -1: + break + + nconfs += 1 - 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() - 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 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 - 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, + print("{} ".format(time)) -print -print "Printing bond distance normalized distribution to",outfile + print("Printing bond distance normalized distribution to {}".format(args.output_file)) -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(args.output_file,"w") as output_file: + rrange = args.rmax - args.rmin + for i in range(0, args.nbins): + output_file.write("{} ".format(args.rmin + rrange*float(i)/float(args.nbins))) + for j in range(0, ntypes): + if (ncount[j] > 0): + output_file.write("{} ".format(float(bin[i][j])/float(ncount[j]*nconfs)/rrange)) + else: + output_file.write("{} ".format(0.0)) + output_file.write("\n") + +if __name__ == "__main__": + compute_bond_distribution() diff --git a/scripts/dihedral_distribute.py b/scripts/dihedral_distribute.py new file mode 100644 index 0000000..7646cc3 --- /dev/null +++ b/scripts/dihedral_distribute.py @@ -0,0 +1,170 @@ +#!/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 +""" + +import sys +import math +from data import data +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*boxlength: + 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 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] +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 + VectorD = [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)) + + 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) + 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,VectorD) + 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 +if argv[1] == "dihedrals": + print("Printing normalized dihedral distributions to {}".format(outfile) +else: + print("Printing normalized improper 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()