From e921bbeeaecc09c5e4925ed7c3826543e4b70d9f Mon Sep 17 00:00:00 2001 From: vijay gopal chilkuri Date: Wed, 10 Mar 2021 14:04:07 +0100 Subject: [PATCH] Added files for generating data. --- generateData.py | 60 +++++++++++ jastlib.py | 271 ++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 331 insertions(+) create mode 100644 generateData.py create mode 100644 jastlib.py diff --git a/generateData.py b/generateData.py new file mode 100644 index 0000000..304c155 --- /dev/null +++ b/generateData.py @@ -0,0 +1,60 @@ +#!/usr/bin/python + +import numpy as np +import sys, getopt + +from jastlib import getCoefList, get_sphere_distribution, generateElsAroundPoints + +def main(argv): + Natom = 2 + Ratio = 5 + inputfile = '' + outputfile = '' + try: + opts, args = getopt.getopt(argv,"ha:r:",["atom=","ratio="]) + except getopt.GetoptError: + print('test.py -a -r ') + sys.exit(2) + for opt, arg in opts: + if opt == '-h': + print('test.py -a -r ') + sys.exit() + elif opt in ("-a", "--atom"): + Natom = int(arg) + elif opt in ("-r", "--ratio"): + Ratio = int(arg) + + + Nelec = Ratio + Nord = 5 + + L1 = 20.0 + n = Natom # number of points + dmin = 0.1 # min dist between points + Ls = np.array([L1,L1,L1]) # lengths of the box + shift = -10.0 + kappa = 2.0 + filename_atom = "geometry.txt" + filename_coeffsa = "jast_coeffs.txt" + (coeffsa, coeffsb, coeffsc) = getCoefList(Nord,n) + coeffsall = np.concatenate((coeffsa,coeffsb,coeffsc)) + + atomList,_,_ = get_sphere_distribution(n, dmin, Ls, maxiter=1e4, allow_wall=True) + + L1 = 15.0 + n = Nelec # number of points + dmin = 0.1 # min dist between points + Ls = np.array([L1,L1,L1]) # lengths of the box + shift = -10.0 + kappa = 2.0 + filename_elec = "elec_coord.txt" + + rlist = generateElsAroundPoints(n,atomList,dmin) + + # Save file + np.savetxt(filename_elec,rlist) + np.savetxt(filename_atom,atomList) + np.savetxt(filename_coeffsa,coeffsall) + +if __name__ == "__main__": + main(sys.argv[1:]) diff --git a/jastlib.py b/jastlib.py new file mode 100644 index 0000000..afd5587 --- /dev/null +++ b/jastlib.py @@ -0,0 +1,271 @@ +import numpy as np +import matplotlib.pyplot as plt + +def get_sphere_distribution(n, dmin, Ls, maxiter=1e4, allow_wall=True): + """Get random points in a box with given dimensions and minimum separation. + + Parameters: + + - n: number of points + - dmin: minimum distance + - Ls: dimensions of box, shape (3,) array + - maxiter: maximum number of iterations. + - allow_wall: whether to allow points on wall; + (if False: points need to keep distance dmin/2 from the walls.) + + Return: + + - ps: array (n, 3) of point positions, + with 0 <= ps[:, i] < Ls[i] + - n_iter: number of iterations + - dratio: average nearest-neighbor distance, divided by dmin. + + Note: with a fill density (sphere volume divided by box volume) above about + 0.53, it takes very long. (Random close-packed spheres have a fill density + of 0.64). + + Author: Han-Kwang Nienhuys (2020) + Copying: BSD, GPL, LGPL, CC-BY, CC-BY-SA + See Stackoverflow: https://stackoverflow.com/a/62895898/6228891 + """ + Ls = np.array(Ls).reshape(3) + if not allow_wall: + Ls -= dmin + + # filling factor; 0.64 is for random close-packed spheres + # This is an estimate because close packing is complicated near the walls. + # It doesn't work well for small L/dmin ratios. + sphere_vol = np.pi/6*dmin**3 + box_vol = np.prod(Ls + 0.5*dmin) + fill_dens = n*sphere_vol/box_vol + if fill_dens > 0.64: + msg = f'Too many to fit in the volume, density {fill_dens:.3g}>0.64' + raise ValueError(msg) + + # initial try + ps = np.random.uniform(size=(n, 3)) * Ls + + # distance-squared matrix (diagonal is self-distance, don't count) + dsq = ((ps - ps.reshape(n, 1, 3))**2).sum(axis=2) + dsq[np.arange(n), np.arange(n)] = np.infty + + for iter_no in range(int(maxiter)): + # find points that have too close neighbors + close_counts = np.sum(dsq < dmin**2, axis=1) # shape (n,) + n_close = np.count_nonzero(close_counts) + if n_close == 0: + break + + # Move the one with the largest number of too-close neighbors + imv = np.argmax(close_counts) + + # new positions + newp = np.random.uniform(size=3)*Ls + ps[imv]= newp + + # update distance matrix + new_dsq_row = ((ps - newp.reshape(1, 3))**2).sum(axis=-1) + dsq[imv, :] = dsq[:, imv] = new_dsq_row + dsq[imv, imv] = np.inf + else: + raise RuntimeError(f'Failed after {iter_no+1} iterations.') + + if not allow_wall: + ps += dmin/2 + + dratio = (np.sqrt(dsq.min(axis=1))/dmin).mean() + return ps, iter_no+1, dratio + +def generateElsAroundPoints(n,LS,dmin): + """ + Parameters: + + - n: number of points + - LS: list of position of all atoms + - dmin: minimum intra block distance + - shift: inter block distance + + Return: + + - r: array (n, 3) of point positions, + + """ + + xs = None + for Ls in LS: + + # Get list of random points around Ls + distrib,a,b = get_sphere_distribution(n,dmin,Ls) + + if xs is None: + xs = distrib[:,0] + ys = distrib[:,1] + zs = distrib[:,2] + else: + xs = np.concatenate((xs,distrib[:,0])) + ys = np.concatenate((ys,distrib[:,1])) + zs = np.concatenate((zs,distrib[:,2])) + + return((np.array((xs,ys,zs))).T) + +def getCoefList(Nord,Natom): + assert(Nord < 11) + dict = { + 0 : lambda x,y:x-y-2, + 1 : lambda x,y:x-y, + 2 : lambda x,y:x-y, + 3 : lambda x,y:x-y, + 4 : lambda x,y:x-y, + 5 : lambda x,y:x-y, + 6 : lambda x,y:x-y, + 7 : lambda x,y:x-y, + 8 : lambda x,y:x-y, + 9 : lambda x,y:x-y, + 10 : lambda x,y:x-y, + 11 : lambda x,y:x-y, + } + count = 0 + for p in range(2,Nord+1): + for k in range(p-1,-1,-1): + lmax = dict[k](p,k) + for l in range(lmax,-1,-1): + if (p-k-l) & 1 is 0: + count += 1 + coeflista = np.random.rand(Nord+1,Natom) + coeflistb = np.random.rand(Nord+1) + coeflistc = np.random.rand(count,Natom) + return (coeflista.reshape((Nord+1)*Natom),coeflistb,coeflistc.reshape(count*Natom)) + #return (coeflista,coeflistb,coeflistc) + +def get_sphere_distribution(n, dmin, Ls, maxiter=1e4, allow_wall=True): + """Get random points in a box with given dimensions and minimum separation. + + Parameters: + + - n: number of points + - dmin: minimum distance + - Ls: dimensions of box, shape (3,) array + - maxiter: maximum number of iterations. + - allow_wall: whether to allow points on wall; + (if False: points need to keep distance dmin/2 from the walls.) + + Return: + + - ps: array (n, 3) of point positions, + with 0 <= ps[:, i] < Ls[i] + - n_iter: number of iterations + - dratio: average nearest-neighbor distance, divided by dmin. + + Note: with a fill density (sphere volume divided by box volume) above about + 0.53, it takes very long. (Random close-packed spheres have a fill density + of 0.64). + + Author: Han-Kwang Nienhuys (2020) + Copying: BSD, GPL, LGPL, CC-BY, CC-BY-SA + See Stackoverflow: https://stackoverflow.com/a/62895898/6228891 + """ + Ls = np.array(Ls).reshape(3) + if not allow_wall: + Ls -= dmin + + # filling factor; 0.64 is for random close-packed spheres + # This is an estimate because close packing is complicated near the walls. + # It doesn't work well for small L/dmin ratios. + sphere_vol = np.pi/6*dmin**3 + box_vol = np.prod(Ls + 0.5*dmin) + fill_dens = n*sphere_vol/box_vol + if fill_dens > 0.64: + msg = f'Too many to fit in the volume, density {fill_dens:.3g}>0.64' + raise ValueError(msg) + + # initial try + ps = np.random.uniform(size=(n, 3)) * Ls + + # distance-squared matrix (diagonal is self-distance, don't count) + dsq = ((ps - ps.reshape(n, 1, 3))**2).sum(axis=2) + dsq[np.arange(n), np.arange(n)] = np.infty + + for iter_no in range(int(maxiter)): + # find points that have too close neighbors + close_counts = np.sum(dsq < dmin**2, axis=1) # shape (n,) + n_close = np.count_nonzero(close_counts) + if n_close == 0: + break + + # Move the one with the largest number of too-close neighbors + imv = np.argmax(close_counts) + + # new positions + newp = np.random.uniform(size=3)*Ls + ps[imv]= newp + + # update distance matrix + new_dsq_row = ((ps - newp.reshape(1, 3))**2).sum(axis=-1) + dsq[imv, :] = dsq[:, imv] = new_dsq_row + dsq[imv, imv] = np.inf + else: + raise RuntimeError(f'Failed after {iter_no+1} iterations.') + + if not allow_wall: + ps += dmin/2 + + dratio = (np.sqrt(dsq.min(axis=1))/dmin).mean() + return ps, iter_no+1, dratio + + +def scalingee(r,kappa=1.0): + return (numpy.ones_like(r) - numpy.exp(-kappa*r))/kappa + +def scalingen(r,kappa=1.0): + return numpy.exp(-kappa*r) + +if False: + Nord = 5 + L1 = 2.0 + n = 2 # number of points + dmin = 0.1 # min dist between points + Ls = np.array([L1,L1,L1]) # lengths of the box + shift = -10.0 + kappa = 2.0 + filename_atom = str(n) + "_geometry.txt" + filename_elec = str(n) + filename_coeffs = str(n) + "_jast_coeffs.txt" + (coeffsa, coeffsb, coeffsc) = getCoefList(Nord,n) + coeffsall = np.concatenate((coeffsa,coeffsb,coeffsc)) + print(coeffsa.shape,coeffsb.shape,coeffsc.shape) + + atomList,_,_ = get_sphere_distribution(n, dmin, Ls, maxiter=1e4, allow_wall=True) + #print(atomList) + + L1 = 5.0 + n = 5 # number of points + dmin = 0.1 # min dist between points + Ls = np.array([L1,L1,L1]) # lengths of the box + shift = -10.0 + kappa = 2.0 + filename_elec = filename_elec + "_" + str(n) + "_elec_coord.txt" + + #rlist = generateBlockRandomPointsAtShftApart(n,L1,dmin,shift) + rlist = generateElsAroundPoints(n,atomList,dmin) + + + # Save file + np.savetxt(filename_elec,rlist) + np.savetxt(filename_atom,atomList) + np.savetxt(filename_coeffs,coeffsall) + + fig = plt.figure() + ax = fig.add_subplot(111, projection='3d') + + xs = rlist.T[0] + ys = rlist.T[1] + zs = rlist.T[2] + ax.scatter(xs, ys, zs, marker='o') + + plt.show() + + rijScaled = np.array([[(lambda xval, yval: np.linalg.norm(xval-yval))(xval,yval) for yval in rlist] for xval in rlist]) + + plt.imshow(rijScaled) + plt.colorbar() + plt.show()