irpjast/jastlib.py

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()