Skip to content

Calc colvar

thmd.calc.colvar

This module contains classes and functions to compute Order parameters, Collective variables,...

Modules:

Classes:

RATIONAL(r0, d0=0.0, n=6, m=12, dmax_tol=0.001)

Create an Object of SWITCHING FUNCTION

Parameters:

  • r0 (float) –

    The r_0 parameter of the switching function

  • d0 (float, default: 0.0 ) –

    The d_0 parameter of the switching function

  • n ( (float, default: 6 ) –

    The n parameter of the switching function

  • m ( (float, default: 12 ) –

    The m parameter of the switching function

Returns:

  • Obj ( object ) –

    Object of the switching function

Notes

Dmin != D0 D0, R0 : are the parameter of switching function Dmin, Dmax : are the bounds at which the switching take affect

Examples:

sw = thmd.SwitchFunc.RATIONAL(r0=6.3, d0=0.0, n=10)

create some intiatial attributes...

Methods:

  • dmax –
  • eval –

    compute & return value and derivation of sw function

Attributes:

r0 = r0 instance-attribute

d0 = d0 instance-attribute

n = n instance-attribute

m = m instance-attribute

dmax(dmax_tol)

eval(x, compute_der=False)

compute & return value and derivation of sw function Input x can be a scalar or a 1d numpy ndarray

HEAVISIDE(r0)

create some intiatial attributes...

Methods:

  • eval –

    compute value of sw function

  • der –

    compute derivative of sw function

Attributes:

dmax = r0 instance-attribute

eval(x: Union[float, list, np.ndarray]) -> np.ndarray

compute value of sw function

der(x: float | list | np.ndarray) -> np.ndarray

compute derivative of sw function

CUBIC(d0, dmax)

Create an Object of SWITCHING FUNCTION

Parameters:

  • d0 (float) –

    The r_0 parameter of the switching function

  • dmax (float) –

    The d_0 parameter of the switching function

Returns:

  • Obj ( object ) –

    Object of the switching function

Examples:

        sw = thmd.SwitchFunc.CUBIC(d0=0.0, dmax=2)

create some intiatial attributes...

Methods:

  • eval –

    compute & return value and derivation of sw function

Attributes:

d0 = d0 instance-attribute

dmax = dmax instance-attribute

eval(x, compute_der=False)

compute & return value and derivation of sw function

Parameters:

  • x (float or ndarray) –

    input values

Returns:

  • f ( float or ndarray ) –

    value of the switching function

  • df ( (float or ndarray, optional) ) –

    value of the derivative of the switching function. Just return if compute_der=True

SMAP(r0, a=10, b=20, d0=0, tol=0.0001)

create some intiatial attributes...

Methods:

Attributes:

a = a instance-attribute

b = b instance-attribute

d0 = d0 instance-attribute

r0 = r0 instance-attribute

tol = tol instance-attribute

eval(x, compute_der=False)

findDmax(tol=None, upper_bound=100, gridSize=0.0001)

findDmin(tol=None, lower_bound=None, gridSize=0.0001)

find_Dmin_Dmax(tol=None, gridSize=0.0001, upper_bound=50)

find Dmin and Dmax of function based on given tolerance

estimate_ab(target_Dmin, target_Dmax, init_a=None, init_b=None, tol=0.001)

Estimate a and b parameters of SMAP function based on given R0, Dmax and Dmin. This version employs scipy.optimize.minimize

Parameters:

  • target_Dmin (float) –

    target Dmin

  • target_Dmax (float) –

    target Dmax

  • search_step (float) –

    incremental value of A/B for searching. Defaults to 0.2.

  • tol (float, default: 0.001 ) –

    tolerance of Dmin/Dmax. Defaults to 1e-3. Should be smaller than 1e-3 (larger value will cause problem)

  • A (int) –

    initial value of A. Defaults to 10.

  • B (int) –

    initial value of B. Defaults to 20.

Returns:

  • A ( float ) –

    found A

  • B ( float ) –

    found B

estimate_ab_old_manual(target_Dmin, target_Dmax, init_a=None, init_b=None, tol=0.001, search_step=0.1)

Estimate a and b parameters of SMAP function based on given R0, Dmax and Dmin. This version use manual search.

Parameters:

  • target_Dmin (float) –

    target Dmin

  • target_Dmax (float) –

    target Dmax

  • search_step (float, default: 0.1 ) –

    incremental value of A/B for searching. Defaults to 0.2.

  • tol (float, default: 0.001 ) –

    tolerance of Dmin/Dmax. Defaults to 1e-3. Should be smaller than 1e-3 (larger value will cause problem)

  • A (int) –

    initial value of A. Defaults to 10.

  • B (int) –

    initial value of B. Defaults to 20.

Returns:

  • A ( float ) –

    found A

  • B ( float ) –

    found B

cv_CoordNum

Functions:

  • coord_number –

    The Coordination is the size of input "Points", this function just weight it with a switching function

coord_number(Points, **kwargs)

The Coordination is the size of input "Points", this function just weight it with a switching function * Compulsory Inputs: ** optional Inputs: switchFunc=[1,1...,1] : Nx1 array, contain values of switching function s(Rj) (Rj is positions of atom j) Returns: coord : scalar, Order Parameter Examples: S = thmd.OrderPara.Coordination([1,0,0; 0,1,0], SW=sw) By Cao Thang, Aug 2020

cv_PairANGLE

Functions:

  • PairANGLE –

    Order Parameter based on pair functions of Angles in the first shell:

PairANGLE(Points, CENTER, SIGMA, **kwargs)

Order Parameter based on pair functions of Angles in the first shell:

Agrs

Points : Nx3 Matrix, contain bonding vectors between neighboring atoms j and ref atom i CENTER=[pi/3, pi/2, 2*pi/3, pi] : list, centers of Gaussians SIGMA =[0.03,0.04,0.04,0.03] : list, sigmas of Gaussians switchFunc=[1,1...,1] : Nx1 array, contain values of switching function s(Rj) (Rj is positions of atom j)

Returns:

  • gamma ( float ) –

    Order Parameter

Examples:

S = thmd.OrderPara.FCCcubic([1,0,0; 0,1,0], SW=sw)
Notes

Require to best chose Rcut for Switching function

Refs
  1. Gobbob et al., "Nucleation of Molecular Crystals Driven by Relative Information Entropy"

cv_Steinhardt

Functions:

  • Ql_Steinhardt –

    compute origincal Stainhardt of l-th order

  • Local_Ql_Steinhardt –

    compute Local Stainhardt of l-th order (modified Steinhardt as: 10.1021/acs.jctc.6b01073)

Ql_Steinhardt(ql_i)

compute origincal Stainhardt of l-th order Args: ql_i : a vector of (2l+1) complex components, qlm(i) vector of atom i Returns: Ql : scalar value of l-th order Stainhardt parameter

Local_Ql_Steinhardt(ql_i, qlm_j, SW)

compute Local Stainhardt of l-th order (modified Steinhardt as: 10.1021/acs.jctc.6b01073) Args: ql_i : 1x(2l+1) array, vector of (2l+1) complex components, qlm(i) vector of atom i qlm_j : Nx(2l+1) array, rows are vectors of (2l+1) complex components, qlm(j) of all neighbors j of atom i Returns: Local_Ql_i : scalar value of l-th order Stainhardt parameter of atom i * PreRequire: compute ql_i complex vector for all atoms before this function can be used

cv_fccCUBIC

Functions:

  • fccCUBIC –

    Function to Calculate FCC CUBIC parameters.

fccCUBIC(points, alpha=27, zDirect='001', switch_function=None)

Function to Calculate FCC CUBIC parameters.

By thangckt, Mar 2020

Parameters:

  • points (Nx3 np.array) –

    array contains bonding vectors between neighboring atoms j and ref atom i

  • alpha (int, default: 27 ) –

    coefficient of harmonic function. Default to 27.

  • zDirect (str, default: '001' ) –

    direction of Z-direction, available '001' '110' '111'. Default to '001'.

  • switch_function (list, default: None ) –

    list contain values of switching function s(Rj) (Rj is positions of atom j). Default to None.

Returns:

  • param ( float ) –

    Order Parameter.

Examples:

S = thmd.colvar.fccCUBIC([1,0,0; 0,1,0], alpha=3, zDirect='001', switch_function=sw)
Notes

Must choose suitable Rcut for Switching function.

cv_localCRYSTALLINITY

Functions:

localCRYSTALLINITY(points: np.array, g_vectors: list[list], switch_function=None)

Function to Calculate Order Parameter with multi_vectors K.

By thangckt, Apr 2019

Parameters:

  • points (Nx3 np.array) –

    array contains bonding vectors between neighboring atoms j and ref atom i g_vectors (tuple): 2d-tuple contains "directions_vectors" for g_vectors (ex: ((4*pi/a)(1,0,0), (4*pi/a)(0,1,0)). The actual g_vectors will be computed in function. Default to ((1,0,0)).

  • switch_function (object, default: None ) –

    switching function s(Rj) (Rj is positions of atom j). Default to None.

Returns:

  • aveLC ( float ) –

    is average Order Parameter , tage over on input g_factors, 0 <= LC <=1

  • LC ( list ) –

    list of real numbers, are Order Parameters corresponding to each g-vector 0 <= LC <=1

  • S ( list ) –

    (not computed) Kx1 vetor of complex numbers, are Static Structure Factors corresponding to each g-vector

Examples:

S = thmd.colvar.localCRYSTALLINITY([1,0,0; 0,1,0], switch_function=sw, zDirect='001')
Notes

If multi g-vectors is input, then OrderPara is take by averaging over all g-vectors.

compute_Gvectors_FCC(g_directions: list[list] = [[1, 0, 0]], lattice_constant: float = 1.0, zDirect: str = '001') -> list[list]

Function to convert reciprocal vectors G to be used in localCRYSTALLINITY.

Parameters:

  • g_directions (list, default: [[1, 0, 0]] ) –

    2d-list contains "directions_vectors" for g_vectors (ex: ((1,0,0), (0,1,0)). The actual g_vectors will be computed. Default to ((1,0,0)).

  • lattice_constant (float, default: 1.0 ) –

    lattice constant of crystal. Default to 1.

  • zDirect (str, default: '001' ) –

    direction of Z-direction, available '001' '110' '111'. Default to '001'.

Returns:

  • g_vectors ( list ) –

    2d-list contains g_vectors

cv_slip_atom

Functions:

  • slip_atom –

    Function to Calculate FCC CUBIC parameters.

slip_atom(points, d_ref, d_diff)

Function to Calculate FCC CUBIC parameters.

Parameters:

  • points (Nx3 np.array) –

    array contains bonding vectors between neighboring atoms j and ref atom i

  • d_ref (float) –

    reference distance.

  • d_diff (float) –

    difference distance to referred as dislocation.

Returns:

  • param ( float ) –

    Order Parameter.

Examples:

S = thmd.colvar.fccCUBIC([1,0,0; 0,1,0], alpha=3, zDirect='001', switch_function=sw)
Notes

Must choose suitable Rcut for Switching function.

find_neighbor

Functions:

  • find_neighbor_gen –

    find Nearest_Neighbors, return generator of Nearest_IDs, "Nearest relative-Position vetors from atom i"

  • find_neighbor_list –

    find Nearest_Neighbors, return list of Nearest_IDs, "Nearest relative-Position vetors from atom i"

find_neighbor_gen(points: Union[np.array, list[list]], box: Union[np.array, list[list]], bound_cond: list = [1, 1, 1], cutoff_neighbor=None, k_neighbor=None, k_cutoff=20, keep_ref=False) -> Generator[list, list[list], None]

find Nearest_Neighbors, return generator of Nearest_IDs, "Nearest relative-Position vetors from atom i" Ver 2: spatial.cKDTree

thangckt, Sep 2019. Update: Aug 2022 to use generator

Parameters:

  • points (array) –

    Nx3 array contain positions of atoms

  • box (array) –

    simulation box

  • bound_cond (list, default: [1, 1, 1] ) –

    boundary condition

  • cutoff_neighbor (float, default: None ) –

    find neighbors within a Cutoff.

  • k_neighbor (int, default: None ) –

    find k nearest neighbors

  • keep_ref (bool, default: False ) –

    include referal-atom in result

Returns:

  • obj ( generator ) –

    this output a GEN contains (Idx_neigh, Rij_vectors)

Examples:

GEN = colvar.find_neighbors_gen(P, box, bound_cond = [1, 1, 1], cutoff_neighbor=9, keep_ref=False)
old version 
access items in generator with:
        for Near_ID, Rij_vector in GEN:
                print (Near_ID, Rij_vector)

- Idx_neigh    : Nx1 list of Mx1-vectors, contain Image_IDs(id of the original atoms before make periodicity) of Nearest atoms
- Rij_vectors : Nx1 list of Mx3-Matrices, contain Nearest Rij relative-Position vetors from Ref.atom i (Nearest Positions)

find_neighbor_list(points, box, bound_cond=[1, 1, 1], cutoff_neighbor=None, k_neighbor=None, k_cutoff=20, keep_ref=False)

find Nearest_Neighbors, return list of Nearest_IDs, "Nearest relative-Position vetors from atom i" Ver 2: spatial.cKDTree

Parameters:

  • points (array) –

    Nx3 array contain positions of atoms

  • box (array) –

    simulation box

  • bound_cond (list, default: [1, 1, 1] ) –

    boundary condition

  • cutoff_neighbor (float, default: None ) –

    find neighbors within a Cutoff.

  • k_neighbor (int, default: None ) –

    find k nearest neighbors

  • keep_ref (bool, default: False ) –

    include referal-atom in result

Returns:

  • Idx_neigh ( array ) –

    Nx1 list of Mx1-vectors, contain Image_IDs(id of the original atoms before make periodicity) of Nearest atoms

  • Rij_vectors ( array ) –

    Nx1 list of Mx3-Matrices, contain Nearest Rij relative-Position vetors from Ref.atom i (Nearest Positions)

Examples:

Idx_neigh, Rij_vectors = colvar.find_neighbors_list(P, box, bound_cond = [1, 1, 1], cutoff_neighbor=9, keep_ref=False)
Notes

don't compute Rij_Bond to save memory Rij_Bonds (np.array): Nx1 list of scalars, contain Rij_bonds from Ref.atom to Nearest_atoms (Nearest-bonds)

sph_harmonics

Functions:

  • yl_i –

    Compute vector of Spherical Harmonics for a set of point (ylm vector have (2l+1) components)

yl_i(l, Rij, SW=None, kind='real', normalization='4pi', deg=False)

Compute vector of Spherical Harmonics for a set of point (ylm vector have (2l+1) components)

Parameters:

  • l (int) ) –

    degree of Spherical Harmonic

  • Rij (array - like) –

    Nx3 array contain Rij of nearest neighbors compute from atom i

  • SW ( (array-like, default: None ) –

    Nx1 values of switching function. Default to 'None'

  • kind ( (str, default: 'real' ) –

    kind of return result. Possible complex/real. Default to complex

  • normalization (str, default: '4pi' ) –

    '4pi', 'ortho', 'schmidt', or 'unnorm' for geodesy 4pi normalized, orthonormalized, Schmidt semi-normalized, or unnormalized spherical harmonic functions, respectively. Default to '4pi'

  • deg (bool, default: False ) –

    If True, theta and phi are expressed in degrees. Default to False

Returns:

  • yl ( array - like ) –

    vector of (2l+1) components

Notes

This functions used the function spharm_lm() from pyshtools

Refs
  1. Visualizing the real forms of the spherical harmonics
  2. In scipy.special.sph_harm function the azimuthal coordinate, theta, comes before the polar coordinate, phi; anh may return complex number only

switch_function

Classes:

RATIONAL(r0, d0=0.0, n=6, m=12, dmax_tol=0.001)

Create an Object of SWITCHING FUNCTION

Parameters:

  • r0 (float) –

    The r_0 parameter of the switching function

  • d0 (float, default: 0.0 ) –

    The d_0 parameter of the switching function

  • n ( (float, default: 6 ) –

    The n parameter of the switching function

  • m ( (float, default: 12 ) –

    The m parameter of the switching function

Returns:

  • Obj ( object ) –

    Object of the switching function

Notes

Dmin != D0 D0, R0 : are the parameter of switching function Dmin, Dmax : are the bounds at which the switching take affect

Examples:

sw = thmd.SwitchFunc.RATIONAL(r0=6.3, d0=0.0, n=10)

create some intiatial attributes...

Methods:

  • eval –

    compute & return value and derivation of sw function

  • dmax –

Attributes:

r0 = r0 instance-attribute
d0 = d0 instance-attribute
n = n instance-attribute
m = m instance-attribute
eval(x, compute_der=False)

compute & return value and derivation of sw function Input x can be a scalar or a 1d numpy ndarray

dmax(dmax_tol)

HEAVISIDE(r0)

create some intiatial attributes...

Methods:

  • eval –

    compute value of sw function

  • der –

    compute derivative of sw function

Attributes:

dmax = r0 instance-attribute
eval(x: Union[float, list, np.ndarray]) -> np.ndarray

compute value of sw function

der(x: float | list | np.ndarray) -> np.ndarray

compute derivative of sw function

CUBIC(d0, dmax)

Create an Object of SWITCHING FUNCTION

Parameters:

  • d0 (float) –

    The r_0 parameter of the switching function

  • dmax (float) –

    The d_0 parameter of the switching function

Returns:

  • Obj ( object ) –

    Object of the switching function

Examples:

        sw = thmd.SwitchFunc.CUBIC(d0=0.0, dmax=2)

create some intiatial attributes...

Methods:

  • eval –

    compute & return value and derivation of sw function

Attributes:

d0 = d0 instance-attribute
dmax = dmax instance-attribute
eval(x, compute_der=False)

compute & return value and derivation of sw function

Parameters:

  • x (float or ndarray) –

    input values

Returns:

  • f ( float or ndarray ) –

    value of the switching function

  • df ( (float or ndarray, optional) ) –

    value of the derivative of the switching function. Just return if compute_der=True

SMAP(r0, a=10, b=20, d0=0, tol=0.0001)

create some intiatial attributes...

Methods:

Attributes:

a = a instance-attribute
b = b instance-attribute
d0 = d0 instance-attribute
r0 = r0 instance-attribute
tol = tol instance-attribute
eval(x, compute_der=False)
findDmax(tol=None, upper_bound=100, gridSize=0.0001)
findDmin(tol=None, lower_bound=None, gridSize=0.0001)
find_Dmin_Dmax(tol=None, gridSize=0.0001, upper_bound=50)

find Dmin and Dmax of function based on given tolerance

estimate_ab(target_Dmin, target_Dmax, init_a=None, init_b=None, tol=0.001)

Estimate a and b parameters of SMAP function based on given R0, Dmax and Dmin. This version employs scipy.optimize.minimize

Parameters:

  • target_Dmin (float) –

    target Dmin

  • target_Dmax (float) –

    target Dmax

  • search_step (float) –

    incremental value of A/B for searching. Defaults to 0.2.

  • tol (float, default: 0.001 ) –

    tolerance of Dmin/Dmax. Defaults to 1e-3. Should be smaller than 1e-3 (larger value will cause problem)

  • A (int) –

    initial value of A. Defaults to 10.

  • B (int) –

    initial value of B. Defaults to 20.

Returns:

  • A ( float ) –

    found A

  • B ( float ) –

    found B

estimate_ab_old_manual(target_Dmin, target_Dmax, init_a=None, init_b=None, tol=0.001, search_step=0.1)

Estimate a and b parameters of SMAP function based on given R0, Dmax and Dmin. This version use manual search.

Parameters:

  • target_Dmin (float) –

    target Dmin

  • target_Dmax (float) –

    target Dmax

  • search_step (float, default: 0.1 ) –

    incremental value of A/B for searching. Defaults to 0.2.

  • tol (float, default: 0.001 ) –

    tolerance of Dmin/Dmax. Defaults to 1e-3. Should be smaller than 1e-3 (larger value will cause problem)

  • A (int) –

    initial value of A. Defaults to 10.

  • B (int) –

    initial value of B. Defaults to 20.

Returns:

  • A ( float ) –

    found A

  • B ( float ) –

    found B

voronoi

Classes:

Functions:

  • layer_extractor –

    Extract atoms of outermost layers, based on Voronoi analysis

  • surface_detect –

    Extract atoms on free surface, are atoms have Voronoi with max(faceArea) >= max_edge**2

Voro3D()

Voro ++ library

Methods:

fAtomicVol_Bulk_gen(P, box, coord_number=False)

compute atomic-volume of each atom in Bulk models Args: P : Nx3 Matrix contain positions of atoms box : simulation box coord_number=False : compute coordiation number Returns: atomicVol : 1xM array of atomic volume coord, cell_neighbor : 1xN list, Nxlist index of neighbors of each atom in original input points

Examples:

gen2 = thmd.fAtomicVol_Bulk_gen(P, box, coord_number=False)   # gen2 = (tomicVol_i, cell_neighbor_i)
gen3 = thmd.fAtomicVol_Bulk_gen(P, box, coord_number=True)    # gen3 = (tomicVol_i, cell_neighbor_i, coord_i)
fAtomicVol_Plate_gen(P, box, bound_cond=(1, 1, 0), max_edge=3.1, surf_deep=15, virtual_surface=True, offset=False, coord_number=False)

compute atomic-volume of each atom in Plate models update Ver2: copy layer 2nd and 3rd to cover free surface, offset distance = inter-layer distance. So, no need to input offset value Args: P : Nx3 Matrix contain positions of atoms box : simulation box

**Optional: bound_cond=(1, 1, 0) : tuple of boundary condtions max_edge=3.1 : value to to compute face-area = max_edge**2 surf_deep=15 : = max(peak) - min(valey), unit in Angstrom, use to reduce input data to save memory virtual_surface=True : apply virtual surface to cover real surface before compute Voronoi + offset=False: extract 2 layers near outmost layer, and use to cover surface, offset_distance is auto computed based on atomic arrangement + offset=true: extract 1 layer near outmost layer, and use to cover surface with offset_distance = offset value coord_number=False : compute coordiation number Returns: atomicVol : 1xM array of atomic volume coord, cell_neighbor : 1xN list, Nxlist index of neighbors of each atom in original input points

fAtomicVol_Bulk(P, box, coord_number=False)

compute atomic-volume of each atom in Bulk models Args: P : Nx3 Matrix contain positions of atoms box : simulation box **Optional: coord_number=False : compute coordiation number Returns: atomicVol : 1xM array of atomic volume coord, cell_neighbor : 1xN list, Nxlist index of neighbors of each atom in original input points

fAtomicVol_Plate(P, box, bound_cond=(1, 1, 0), max_edge=3.1, surf_deep=15, virtual_surface=True, offset=False, coord_number=False)

compute atomic-volume of each atom in Plate models update Ver2: copy layer 2nd and 3rd to cover free surface, offset distance = inter-layer distance. So, no need to input offset value Args: P : Nx3 Matrix contain positions of atoms box : simulation box

**Optional: bound_cond=(1, 1, 0) : tuple of boundary condtions max_edge=3.1 : value to to compute face-area = max_edge**2 surf_deep=15 : = max(peak) - min(valey), unit in Angstrom, use to reduce input data to save memory virtual_surface=True : apply virtual surface to cover real surface before compute Voronoi + offset=False: extract 2 layers near outmost layer, and use to cover surface, offset_distance is auto computed based on atomic arrangement + offset=true: extract 1 layer near outmost layer, and use to cover surface with offset_distance = offset value coord_number=False : compute coordiation number Returns: atomicVol : 1xM array of atomic volume coord, cell_neighbor : 1xN list, Nxlist index of neighbors of each atom in original input points

layer_extractor(P, bound_cond=(1, 1, 0), layer_num=1, max_edge=3.0, surf_deep=15, method='max_face_perimeter')

Extract atoms of outermost layers, based on Voronoi analysis Args: P : Nx3 Matrix contain positions of atoms bound_cond=(1, 1, 0) : tuple of boundary condtions layer_num : number of Layers need to extract (layer_num=0 will extract all layers) max_edge : value to to compute face-area = max_edge**2 surf_deep : = max(peak) - min(valey), unit in Angstrom, use to reduce input data to save memory method : 'max_face_perimeter' or 'max_face_area' Returns: hiLayerIndex, loLayerIndex: list of lists (1xM index of atoms in each layer) By Cao Thang, Jan 2020

surface_detect(P, bound_cond=(1, 1, 0), max_edge=3.0, surf_deep=15, method='max_face_perimeter')

Extract atoms on free surface, are atoms have Voronoi with max(faceArea) >= max_edge**2 Agrs: P : Nx3 Matrix contain positions of atoms bound_cond=(1, 1, 0) : tuple of boundary condtions max_edge : value to to compute face-area = max_edge**2 surf_deep : = max(peak) - min(valey), unit in Angstrom, use to reduce input data to save memory method : 'max_face_perimeter' or 'max_face_area'

Returns:

  • –

    hiSurIndex, loSurIndex: 1xM array, index of surface atoms in the original input points

Notes

experimental choose: max_edge=0.73*latticeConst only 1 pair of surface is detect each time