diff --git a/examples/neuronet_test/analyse.py b/examples/neuronet_test/analyse.py index 9e416d8ec..ebe110951 100644 --- a/examples/neuronet_test/analyse.py +++ b/examples/neuronet_test/analyse.py @@ -1,201 +1,295 @@ #!/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, pot, shift_x=0, shift_y=0, shift_z=0): +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: - input = get_sceleton(del_x, pot, shift_x, shift_y, shift_z) + 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 ("For {}:".format(name)) + 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 = [("Ryo", del_x, '' ), - ("Me", del_x, 'c')] + cases = [("Old", del_x, '' ), + ("New", del_x, 'c')] for case in cases: analyse_case(*case) - print("\n\nShifting solid motion\n") + 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()