Page MenuHomec4science
Files F74292316

analyse.py
No OneTemporary
Actions

File Metadata

Created
Fri, Jul 26, 23:15

analyse.py
View Options

#!/usr/bin/env python3
import numpy as np
import subprocess as sub
import shlex
import matplotlib.pyplot as plt
sceleton = """
# 3d Lennard-Jones melt
units metal
atom_style atomic
atom_modify map array
timestep 0
region box block 0 12 0 12 0 12
boundary p p p
create_box 3 box
create_atoms 1 single {pos1} units box
create_atoms 2 single {pos2} units box
create_atoms 3 single {pos3} units box
create_atoms 3 single {pos4} units box
mass 1 1.0
mass 2 1.0
mass 3 1.0
pair_style {pot}neuronet
pair_coeff * * in.const.NN_short in.params.NN_short
neighbor 0.3 bin
compute 1 all pe
variable pppp equal c_1
variable f1x equal "fx[1]"
variable f1y equal "fy[1]"
variable f1z equal "fz[1]"
variable f2x equal "fx[2]"
variable f2y equal "fy[2]"
variable f2z equal "fz[2]"
variable r1x equal "x[1]"
variable r1y equal "y[1]"
variable r1z equal "z[1]"
variable r2x equal "x[2]"
variable r2y equal "y[2]"
variable r2z equal "z[2]"
fix 1 all nve
fix 2 all print 1 "POTENTIAL ENERGY DIRECT ${pppp}"
fix 3 all print 1 "FORCE VECTOR DIRECT ${{f1x}} ${{f1y}} ${{f1z}} ${{f2x}} ${{f2y}} ${{f2z}} "
fix 4 all print 1 "POSITION VECTOR DIRECT ${{r1x}} ${{r1y}} ${{r1z}} ${{r2x}} ${{r2y}} ${{r2z}} "
dump 1 all custom 1 dump.out id xs ys zs fx fy fz
run 2
"""
def get_sceleton(del_x = 0., pot='', shift_x=0, shift_y=0, shift_z=0):
box_siz = 12.
pos1 = "{} {} {}".format((0 + shift_x)%box_siz, (0 + shift_y )%box_siz, (del_x + shift_z)%box_siz)
pos2 = "{} {} {}".format((3 + shift_x)%box_siz, (3 + del_x + shift_y)%box_siz, ( 0 + shift_z)%box_siz)
pos3 = "{} {} {}".format((3 + shift_x)%box_siz, (0 + shift_y )%box_siz, ( 3 + shift_z)%box_siz)
pos4 = "{} {} {}".format((3 + shift_x)%box_siz, (3 + shift_y )%box_siz, ( 3 + shift_z)%box_siz)
return sceleton.format(pos1=pos1, pos2=pos2, pos3=pos3, pos4=pos4, pot=pot, pppp='{pppp}')
def get_sceleton_rot(rot_mat, pot):
box_siz = 12.
pos = np.array([[0, 3, 3, 3],
[0, 3, 0, 3],
[0, 0, 3, 3]], dtype=float)
pos = np.dot(rot_mat, pos)%box_siz
pos1 = "{} {} {}".format(*pos[:, 0].ravel().tolist())
pos2 = "{} {} {}".format(*pos[:, 1].ravel().tolist())
pos3 = "{} {} {}".format(*pos[:, 2].ravel().tolist())
pos4 = "{} {} {}".format(*pos[:, 3].ravel().tolist())
return sceleton.format(pos1=pos1, pos2=pos2, pos3=pos3, pos4=pos4, pot=pot, pppp='{pppp}')
def nextline(fh, do_print=False):
line = fh.__next__()
if do_print:
print(line)
return line
def get_force_del_pos(fname = "dump.out"):
do_print = False
with open(fname, 'r') as fh:
line = ''
while not line.startswith("ITEM: NUMBER OF ATOMS"):
line = nextline(fh, do_print)
nb_atoms = int(nextline(fh,do_print))
while not line.startswith("ITEM: ATOMS id xs ys zs fx fy fz"):
line = nextline(fh, do_print)
force = []
pos_init = []
for i in range(nb_atoms):
vals = [float(val) for val in nextline(fh, do_print).split()[-6:]]
pos_init += vals[:3]
force += vals[-3:]
pos_fin = []
line = ''
while not line.startswith("ITEM: ATOMS id xs ys zs fx fy fz"):
line = nextline(fh, do_print)
for i in range(nb_atoms):
vals = [float(val) for val in nextline(fh, do_print).split()[-6:]]
pos_fin += vals[:3]
force = np.array(force)
pos_init = 12*np.array(pos_init)
pos_fin = 12*np.array(pos_fin)
del_pos = pos_fin - pos_init
print("force = {}".format (force))
print("pos = {}".format (pos_init))
return force, pos_init
def get_energy (fname = "log.lammps"):
do_print = False
with open(fname, 'r') as fh:
line = ""
while not line.startswith("potential energy"):
line = nextline(fh, do_print)
epot_init = float(line.split()[-1])
print("E_pot = {}".format(epot_init))
return epot_init
def get_everything(output):
do_print = False
def nextliner(do_print):
for line in output.decode().split('\n'):
if do_print:
print(line)
yield line
gen = nextliner(do_print)
def nextline(do_print):
return gen.__next__()
line = ""
while not line.startswith("POTENTIAL ENERGY DIRECT"):
line = nextline(do_print)
e_pot = float(line.split()[-1])
force = np.array([float(item) for item in nextline(do_print).split()[3:]])
pos = np.array([float(item) for item in nextline(do_print).split()[3:]])
return e_pot, force, pos
def eval_derivative(del_pos, force, del_epot):
estim = np.dot(-force, del_pos)
print("Δx = {}".format(del_pos))
print("ΔE = {}, Δx.-f = {}, err = {}".format(
del_epot, estim, 1-estim/del_epot))
def run(del_x=0, pot='', shift_x=0, shift_y=0, shift_z=0, rot_mat=None, nb_p=None):
command = "../../src/lmp_serial"
if nb_p is not None:
command = "mpirun -np {} ".format(nb_p) + "../../src/lmp_mpi"
with sub.Popen(shlex.split(command),
stdin=sub.PIPE, stdout=sub.PIPE) as proc:
if rot_mat is None:
input = get_sceleton(del_x, pot, shift_x, shift_y, shift_z)
else:
input = get_sceleton_rot(rot_mat, pot)
#print(input)
out, err = proc.communicate(input=input.encode())
proc.wait()
#print(out.decode())
return out
def analyse_case(name, del_x, pot):
print ("For {}:".format(name))
output = run(0., pot)
e_i, f, x_i = get_everything(output)
print("x_i = {}".format(x_i))
output = run (del_x, pot)
e_f, dummy, x_f = get_everything(output)
print("x_f = {}".format(x_f))
print("E_i = {}, E_f = {}".format(e_i, e_f))
eval_derivative(x_f-x_i, f, e_f - e_i)
def shift_case(name, del_x, pot, dir, nb_step):
print ("\nFor {}:".format(name))
e_pots = np.zeros(nb_step, dtype=float)
errs = np.zeros_like(e_pots)
estims = np.zeros_like(e_pots)
for i in range(nb_step):
shifts = tuple((d*i for d in dir))
output = run(0, pot, *shifts)
e_pots[i], f, x_i = get_everything(output)
output = run(del_x, pot, *shifts)
e_f, dummy, x_f = get_everything(output)
estims[i] = np.dot(-f, x_f-x_i)
errs[i] = 1-estims[i]/(e_f - e_pots[i])
del_epots = e_pots-e_pots[0]
print("E_pot0 = {}".format(e_pots[0]))
print("ΔE_pot = {}".format(del_epots))
print("estimated = {}".format(estims))
print("force_err = {}".format(errs))
if(any(abs(errs) > 1e-6)) :
print("FAILURE! Force calculation wrong for {}".format(name))
return del_epots
def rotate_case(name, pot, rot_mat):
print ("\nFor {}:".format(name))
r_mats = (np.eye(3), rot_mat)
f = []
e = []
for r_mat in r_mats:
output = run(pot=pot, rot_mat=r_mat)
epot, force, dummy = get_everything(output)
e.append(epot)
force.shape = (-1, 3)
f.append(force.T)
err_E = 1-e[1]/e[0]
err_f = abs(f[0]-np.dot(rot_mat.T, f[1])).max()
print("err_E = {}".format(err_E))
print("err_f = {}".format(err_f))
def proc_case(name, pot):
print ("\nFor {}:".format(name))
nb_procs = (None, 8, 27)
f = []
e = []
for nb_p in nb_procs:
output = run(pot=pot, nb_p=nb_p)
epot, force, dummy = get_everything(output)
e.append(epot)
force.shape = (-1, 3)
f.append(force.T)
e = np.array(e)
err_E = 1-e/e[0]
err_f = np.array([abs(force - f[0]).max() for force in f])
print("err_E = {}".format(err_E))
print("err_f = {}".format(err_f))
def get_rot_mat(invert=False):
def normalize(vec):
vec /= np.sqrt((vec**2).sum())
factor = -1 if invert else 1
rot_mat = np.random.random((3, 3))
normalize(rot_mat[0,:])
rot_mat[1, :] = np.cross(rot_mat[0, :], rot_mat[1, :])
normalize(rot_mat[1, :])
rot_mat[2, :] = factor*np.cross(rot_mat[0, :], rot_mat[1, :])
#print("NORM = {}".format(np.sqrt((rot_mat[2,:]**2).sum())))
#print(rot_mat)
#print(np.linalg.det(rot_mat))
#print(np.dot(rot_mat, rot_mat.T))
return(rot_mat)
def main():
del_x = 1.e-8
cases = [("Old", del_x, '' ),
("New", del_x, 'c')]
for case in cases:
analyse_case(*case)
print("\nShifting solid motion")
dir = [3, 4, 6.48]
for case in cases:
shift_case(*case, dir, 4)
print("\nRotating solid motion")
rot_mat = get_rot_mat()
cases = [("Old", '' , rot_mat),
("New", 'c', rot_mat)]
for case in cases:
rotate_case(*case)
print("\nInverting (chiral) 'motion'")
rot_mat = get_rot_mat(invert = True)
cases = [("Old", '' , rot_mat),
("New", 'c', rot_mat)]
for case in cases:
rotate_case(*case)
print("\nIntroducing processor boundaries")
cases = [("Old", '' ),
("New", 'c')]
for case in cases[1:]:
proc_case(*case)
if __name__ == "__main__":
main()

Event Timeline

c4science · Help