Skip to content
Snippets Groups Projects
Commit 512f5869 authored by Germain Clavier's avatar Germain Clavier
Browse files

Suppressed plot_table for git_forces?

parent eb85ff76
No related branches found
No related tags found
No related merge requests found
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
'''
Plot LAMMPS forces
'''
import argparse
import numpy as np
import os
import logging
import sys
import GCcolors
from matplotlib import pyplot as plt
def compute_energy(tp):
r = tp[0]
fo = tp[2]
e = np.zeros(r.shape)
for i, (ri, fi) in enumerate(zip(r, fo)):
if i == 0:
continue
dr = ri - r[i-1]
e[i] = e[i-1]-dr*fo[i-1]
e -= e[-1]
return e
def main():
parser = argparse.ArgumentParser(
description="Plots Lammps tabulated forces."
)
parser.add_argument(
"-v",
"--verbose",
dest="verbosity",
default=0,
action="count",
help="display more information",
)
parser.add_argument(
"-f",
"--file",
dest="infile",
default="",
help="File to read",
)
parser.add_argument(
"-x",
dest="xrange",
default="",
help="xrange : separated (for negative values: -x=-3:10)",
)
parser.add_argument(
"-y",
dest="yrange",
default="",
help="yrange : separated",
)
parser.add_argument(
"-e",
dest="extract",
action="store_true",
help="Extract the forces in separate files",
)
args = parser.parse_args()
##########
# Manage arguments
# -v/--verbose
try:
utilsscript.init_logging(args.verbosity)
except NameError:
logging.basicConfig(level=logging.WARNING)
logging.warning("utilsscript lib not found, using default logging at warning level.")
infile = args.infile
if not os.path.isfile(infile):
logging.error("Input file not found")
sys.exit(1)
toplot = []
with open(infile, 'r') as f:
lines = iter(f.readlines())
while True:
try:
r = []
force = []
ener = []
tok = []
while not tok:
tok = next(lines).partition('#')[0].rstrip()
logging.info("Found {} token".format(tok))
infos = next(lines).split()
npoints = int(infos[1])
dummy = next(lines)
for i in range(npoints):
line = next(lines).split()
r.append(float(line[1]))
ener.append(float(line[2]))
force.append(float(line[3]))
r = np.array(r)
ener = np.array(ener)
force = np.array(force)
toplot.append([r, ener, force, tok])
tok = []
dummy = next(lines)
except StopIteration:
break
for tp in toplot:
tp[1] = compute_energy(tp)
fig, axes = plt.subplots(1,2)
for tp in toplot:
axes[0].plot(tp[0], tp[1], label=tp[3], linewidth=3)
axes[1].plot(tp[0], tp[2], label=tp[3], linewidth=3)
if args.xrange:
xmin, xmax = list(map(float, args.xrange.split(":")))
axes[0].set_xlim(xmin, xmax)
axes[1].set_xlim(xmin, xmax)
if args.yrange:
ymin, ymax = list(map(float, args.yrange.split(":")))
axes[0].set_ylim(ymin, ymax)
axes[1].set_ylim(ymin, ymax)
# Setting axes 0
axes[0].set_title("Estimated energy", fontsize=30)
axes[0].legend(frameon=False, fontsize=30) # Fat font, no frame
axes[0].set_xlabel("r [A]", fontname="Hack", fontsize=30) # ylabel name, font is Hack, size 30pts
axes[0].set_ylabel("E [kcal/mol]", fontname="Hack", fontsize=30) # ylabel name, font is Hack, size 30pts
axes[0].tick_params(axis='both', which='major', labelsize=30) # Biggers ticks, bigger ticks label!
# Setting axes 1
axes[1].set_title("Tabulated force", fontsize=30)
axes[1].legend(frameon=False, fontsize=30) # Fat font, no frame
axes[1].set_xlabel("r [A]", fontname="Hack", fontsize=30) # ylabel name, font is Hack, size 30pts
axes[1].set_ylabel("F [kcal/mol/A]", fontname="Hack", fontsize=30) # ylabel name, font is Hack, size 30pts
axes[1].tick_params(axis='both', which='major', labelsize=30) # Biggers ticks, bigger ticks label!
plt.subplots_adjust(wspace=0.3)
plt.show()
if args.extract:
for tp in toplot:
outfile = "".join([tp[3], '.plot'])
logging.info("Writing file {}".format(outfile))
with open(outfile, 'w') as f:
f.write("# {} force extracted from {}\n".format(tp[3], infile))
f.write("# {:^20} {:^20} {:^20}\n".format('r', 'energy', 'force'))
for a, b, c in zip(tp[0], tp[1], tp[2]):
f.write("{:>18.16e} {:>18.16e} {:>18.16e}\n".format(a, b, c))
return
if __name__ == "__main__":
try:
main()
except KeyboardInterrupt:
raise SystemExit("User interruption.")
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment