API¶

`thmd` ¶

Python package for building models, pre-processing and post-processing Molecular Dynamics simulations.

Developed and maintained by C.Thang Nguyen

Modules:

calc –
io –

This module contains classes and functions to read/write data in various formats.
latex –
model –
plot –

This module provides functions and class to produce publication-quality figures. The main idea is to create a handy class which setup all necessary parameters for plotting, and then just call its methods to produce the publication-quality figures.
post –
qflow –
recipe –
util –

Attributes:

__description__ –
__long_description__ –
__author__ –
mySTR –

`description = 'Python package'` `module-attribute` ¶

`__long_description__ = 'Python package for pre- and post-processing data from MD simulations.'` `module-attribute` ¶

`author = 'thangckt'` `module-attribute` ¶

mySTR = '\n```shell\n thmd_package\n │ README.md\n │ LICENSE.md\n │ setup.py\n │\n │ thmd\n │ init.py\n │ data.py\n │\n │ io\n │ │ init.py\n │ │__ define_script.py\n │ │ LmpFrame.py\n │ │ ...\n │\n │ model\n │ │ init.py\n │ │__ box_orientation.py\n │ │ crystal3D.py\n │ │ ...\n │\n │ colvar\n │ │ init.py\n │ │__ cv_fccCUBIC.py\n │ │__ cv_localCRYSTALLINITY.py\n │ │ ...\n │\n │ free_energy\n │ │ init.py\n │ │ Helmholtz_excess_UF.py\n │ │__ replica_logPD_intergration.py\n │ │ ...\n │\n │ utils\n │ │ init.py\n │ │ coord_rotation.py\n │ │__ unit_convert.py\n │ │__ compute_distance.py\n │ │ fitting.py\n │ │ ...\n │\n```\n' `module-attribute` ¶

`__version` ¶

Attributes:

TYPE_CHECKING –
VERSION_TUPLE –
version (str) –
__version__ (str) –
__version_tuple__ (VERSION_TUPLE) –
version_tuple (VERSION_TUPLE) –

`TYPE_CHECKING = False` `module-attribute` ¶

`VERSION_TUPLE = Tuple[Union[int, str], ...]` `module-attribute` ¶

`version: str = '0.1.dev533+g053f72e.d20250107'` `module-attribute` ¶

`__version_tuple__: VERSION_TUPLE = (0, 1, 'dev533', 'g053f72e.d20250107')` `module-attribute` ¶

`version_tuple: VERSION_TUPLE = (0, 1, 'dev533', 'g053f72e.d20250107')` `module-attribute` ¶

`calc` ¶

Modules:

colvar –

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

`colvar` ¶

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

Modules:

cv_CoordNum –
cv_PairANGLE –
cv_Steinhardt –
cv_fccCUBIC –
cv_localCRYSTALLINITY –
cv_slip_atom –
find_neighbor –
sph_harmonics –
switch_function –
voronoi –

Classes:

RATIONAL –

Create an Object of SWITCHING FUNCTION
HEAVISIDE –
CUBIC –

Create an Object of SWITCHING FUNCTION
SMAP –

`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 –
d0 –
n –
m –

`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 –

`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 –
dmax –

`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:

eval –
findDmax –
findDmin –
find_Dmin_Dmax –

find Dmin and Dmax of function based on given tolerance
estimate_ab –

Estimate a and b parameters of SMAP function based on given R0, Dmax and Dmin.
estimate_ab_old_manual –

Estimate a and b parameters of SMAP function based on given R0, Dmax and Dmin.

Attributes:

a –
b –
d0 –
r0 –
tol –

`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

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 –

Function to Calculate Order Parameter with multi_vectors K.
compute_Gvectors_FCC –

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

`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

Visualizing the real forms of the spherical harmonics
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 –

Create an Object of SWITCHING FUNCTION
HEAVISIDE –
CUBIC –

Create an Object of SWITCHING FUNCTION
SMAP –

`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 –
d0 –
n –
m –

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 –

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 –
dmax –

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:

eval –
findDmax –
findDmin –
find_Dmin_Dmax –

find Dmin and Dmax of function based on given tolerance
estimate_ab –

Estimate a and b parameters of SMAP function based on given R0, Dmax and Dmin.
estimate_ab_old_manual –

Estimate a and b parameters of SMAP function based on given R0, Dmax and Dmin.

Attributes:

a –
b –
d0 –
r0 –
tol –

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:

Voro3D –

Voro ++ library

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 –

compute atomic-volume of each atom in Bulk models
fAtomicVol_Plate_gen –

compute atomic-volume of each atom in Plate models
fAtomicVol_Bulk –

compute atomic-volume of each atom in Bulk models
fAtomicVol_Plate –

compute atomic-volume of each atom in Plate models

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 2^nd and 3^rd 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 2^nd and 3^rd 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

`free_energy` ¶

Modules:

TImethod –
logmfd –

`TImethod` ¶

Functions:

AS_integration –

Compute Free energy difference in TI method - Lambda integration along Adiabatic Switching path
RS_integration –

Compute Free energy difference in TI method - Lambda integration along Reversible Scaling path
Helmholtz_excess_UF –

this func. is to compute the excess Helmholtz freeEnergy as Eq.(25) by R.Paula Leite 2016.

Attributes:

kB –

`kB = sc.value('Boltzmann constant in eV/K')` `module-attribute` ¶

`AS_integration(forwFile, backFile)` ¶

Compute Free energy difference in TI method - Lambda integration along Adiabatic Switching path

Parameters:

forwFile (str) –

file name of forward integration (must contain 2 first columns are: dE, lambda)
backFile (str) –

file name of backward integration (must contain 2 first columns are: dE, lambda)

Returns:

W ( float ) –

irreversible TI work
Q ( float ) –

dissipation

`RS_integration(forwFile, backFile, T0=0, F0=0, kB=kB)` ¶

Compute Free energy difference in TI method - Lambda integration along Reversible Scaling path

Parameters:

forwFile (str) –

file name of forward integration (must contain 2 first columns are: dE, lambda)
backFile (str) –

file name of backward integration (must contain 2 first columns are: dE, lambda)
T0 (float, default: 0 ) –

reference temperature. Default is 0.
F0 (float, default: 0 ) –

free energy at reference temperature. Default is 0.
kB (float, default: kB ) –

Boltzmann constant (chose unit to consist with F0). Default is kB in unit eV/K.

Returns:

T ( array ) –

temperature
Ft ( array ) –

free energy as a function of temperature.
W ( array ) –

cumulative work

`Helmholtz_excess_UF(p, x)` ¶

this func. is to compute the excess Helmholtz freeEnergy as Eq.(25) by R.Paula Leite 2016. J.Chem.Phys.145,no.19:194101. https://doi.org/10.1063/1.4967775. p : # UFM p-parameter x = b*rho : # adimensional variable, involved UFM-sigma parameter, rho in unit [1/A^3]

the excess_Helmholtz free enegy in eV : (beta*Fexc)/N

the pressure : beta*b*P

`logmfd` ¶

Classes:

LogMFD_FCCUBIC –

Create an Object (class) of Potential, contain some pre-setup information for Energy barrier for LogMFD calculation.

Functions:

replica_logPD_integration –

The function to compute LogPD-based MeanForce
replica_MD_average –

compute Replica_MD_Average from output of MD.
exp_normalize –
read_df –

define lamda(x)
replica_SteerMD –

compute Average Work from output of SteerMD.
find_basin_logmfd_profile –

Args:

`LogMFD_FCCUBIC(Element, Potential_Name, modelType='Bulk', zDirect='001')` ¶

Create an Object (class) of Potential, contain some pre-setup information for Energy barrier for LogMFD calculation. Energy barrier of System does not depend on CV (same value for energy CV: meanCV, nAtomLiquid,...) * Attributes: Element : Cutoff : force cutoff of potential Potential_Name : * Energy barrier Coeff: an estimation of per-atom energy barrier for melting, for a system of N atoms 𝐹=𝑓∗𝑁^(⅔) , then perAtom barrier 𝑓=𝐹/𝑁^(⅔) This method return f(T) * histo_point_coeff: the CV-value at which histogram of solid and liquid meets, coeff of function f(T) = a0 + a1T

Methods: Latt_Const : compute lattice constant at a specific temperature

Element : 'Al', 'Cu',... Potential_Name : 'Cu' : 'Mishin-2001'; 'Foiles-1986';... 'Al' : 'Mishin-1999'; 'Sheng-2011';... modelType : 'Bulk' or 'Plate' zDirect : '001' or '110' or '111' or thickness : thickness of plate

Methods:

meltingBarrier –

for a system of N atoms 𝐹=𝑓∗𝑁^(⅔) , then perAtom barrier 𝑓=𝐹/𝑁^(⅔). This method return f(T)
histo_point –

Value of CV at intersection point of histogram as function y = a + bx

Attributes:

Element –
Potential_Name –
melt_barrier_coeff –
histo_point_coeff –

Element = Element instance-attribute ¶

Potential_Name = Potential_Name instance-attribute ¶

melt_barrier_coeff = melt_barrier_coeff[melt_barrier_key] instance-attribute ¶

histo_point_coeff = histo_point_coeff instance-attribute ¶

meltingBarrier(Temp) ¶

for a system of N atoms 𝐹=𝑓∗𝑁^(⅔) , then perAtom barrier 𝑓=𝐹/𝑁^(⅔). This method return f(T)

histo_point(Temp) ¶

Value of CV at intersection point of histogram as function y = a + bx

`replica_logPD_integration(logmfd_files: list[str], replica_files: list[str], beta: float = 1.0)` ¶

The function to compute LogPD-based MeanForce

Parameters:

logmfd_files (list[str]) –

list of logmfd.out files
replica_files (list[str]) –

list of replica.out files
beta (float, default: 1.0 ) –

kB is Boltzmann constant (can be set to 1.0, regardless of kB unit). Defaults to beta = 1.0/(TEMP*kB)

Returns:

file ( obj ) –

contains logPD-based MeanForce

Examples:

free_energy.replica_logPD_intergration(logmfd_files, replica_files)

Requisites:

1. Run logMFD simulations to produce `replica_*/logmfd.out` and `replica_*/replica.out`

```
<logmfd.out>
1:iter_md, 2:Flog(t), …, 6: X(t), 7: V(t), 8: Force(t)
1   F(1), …, X(1), V(1), Force(0)
2   F(3), …, X(2), V(2), Force(1)
```

```
<replica.out>
iter_md, work, weight, cv
1  work(1)   weight(1)  cv(0)
2  work(2)   weight(2)  cv(1)
```

Notes

About the printed values in <replica.out> and <logmfd.out> as emails replied by Tetsuya Morishita. (check thangckt email)

Specify type of function cumulative_trapezoid:np.float64 to be used in numba

cumulative_trapezoid:np.float64 = cumulative_trapezoid,   # this to specify `type` for function `cumulative_trapezoid`
logMFD:pd.DataFrame = read_data.logMFD,
matrix_auto:pd.DataFrame = read_data.matrix_auto):

Refs

[1].https://pubs.acs.org/doi/10.1021/acs.jctc.7b00252 Free Energy Reconstruction from Logarithmic Mean-Force Dynamics Using Multiple Nonequilibrium TrajectoriesFree [2] Exp-normalize trick: https://timvieira.github.io/blog/post/2014/02/11/exp-normalize-trick/

2023 Apr

CV-value in replica.out and logmfd.out are the same, so don't need to shift value in replica.out as below:

Force = [elem['CV1_force'].shift(-1) for elem in logmfd ]

`replica_MD_average(MD_out_files: list[str])` ¶

compute Replica_MD_Average from output of MD.

Parameters:

MD_out_files (list[str]) –

list of "MDout_replica.txt" files

Returns:

–

logPD file: contains logPD-based MeanForce

Requisites

Replica_* files from separate MD simulations

Thang, Jul2020 (update: Sep 2021)

`exp_normalize(x)` ¶

`read_df(file, engine='LAMMPS')` ¶

define lamda(x)

`replica_SteerMD(SteerMD_files, beta=1.0, engine='Lammps')` ¶

compute Average Work from output of SteerMD.

Parameters:

SteerMD_files –

|list| of "SteerMD.txt" files
beta (= 1.0/(TEMP*kB, default: 1.0 ) –

kB is Boltzmann constant (can be set to 1.0, regardless of kB unit)

Returns:

–

aveSteerMD file: contains logPD-based MeanForce

Requisites: 1. Replica_* files from separate MD simulations

Refs

[1]. https://github.com/sandeshkalantre/jarzynski/blob/master/code/Simulations%20on%20Harmonic%20Oscillator%20Model.ipynb [2]. https://www.plumed.org/doc-v2.6/user-doc/html/belfast-5.html#belfast-5-work [3]. Exp-normalize trick: https://timvieira.github.io/blog/post/2014/02/11/exp-normalize-trick/

`find_basin_logmfd_profile(ldf: list[pd.DataFrame], mfd_interval: list[float], hill_region: list = [0.2, 0.6]) -> list[pd.DataFrame]` ¶

Parameters:

ldf (list[DataFrame]) –

list of pd.DataFrames of LogMFD.out files
mfd_interval (list) –

interval of MFDstep in PLUMED input files
hill_region (1x2 list, default: [0.2, 0.6] ) –

the region of CV that contains saddle point. The Left/Right basin is located at left/right of this region. This region should exclude the bassins.

Returns:

df_min_right ( DataFrame ) –

DataFrame of right basin minima
df_min_left ( DataFrame ) –

DataFrame of left basin minima
df_max ( DataFrame ) –

DataFrame of maxima energy
df_diff ( DataFrame ) –

DataFrame of energy difference between minima and maxima

`io` ¶

This module contains classes and functions to read/write data in various formats.

Modules:

read_block –
read_data –

This module contains functions to read numeric data from various formats of TEXT files.
script –

This module contains functions to read/write some specific scripts.
traj –

This module contains classes and functions to process molecular dynamics trajectories.

`read_block` ¶

Classes:

LmpLogFile –

Create an Object of LOG file.
LmpRDF –

class to read Radial Distribution Fuction (RDF) file from Lammps compute
LmpAveChunk –

class to read Radial Distribution Fuction (RDF) file from Lammps compute
PlumHistogram –

Create an Object of DUMP file

Functions:

average_df_list –

compute average of list of DataFrame

`LmpLogFile(logfile=None)` ¶

Create an Object of LOG file.

Notes

run 0 without data

Udpate

2024-10-10: use Polars DataFrame instead of Pandas Series

Parameters:

logfile (str, default: None ) –

file_name of LOG file

Returns:

Obj ( LmpLogFile ) –

LmpLogFile object

Methods:

read_log –

Read LAMMPS logfile

`read_log(logfile)` ¶

Read LAMMPS logfile Args: logfile (str): input LOG file

Returns:

Obj ( LmpLogFile ) –

LmpLogFile object

`LmpRDF(file_name: str)` ¶

class to read Radial Distribution Fuction (RDF) file from Lammps compute

Attributes:

file_name (str) –

file name
frame (DataFrame) –

3d pandas Frame (multi-row-index DataFrame)

Methods:

ReadRDF –

read RDF file
AverageRDF –

the Average RDF

initiate object

Parameters:

file_name (str) –

file_name

Returns:

Obj ( LmpAveChunk ) –

LmpAveChunk object

Examples:

RDF = LmpAveChunk('rdf.txt')

Methods:

read_RDF –

Args:
compute_AveRDF –

compute average of RDF over all frames

`read_RDF(file_name: str)` ¶

Parameters:

file_name (str) –

input RDF file

Returns:

Obj ( LmpRDF ) –

LmpRDF object

`compute_AveRDF()` ¶

compute average of RDF over all frames

Returns:

df ( DataFrame ) –

Average of RDF

`LmpAveChunk(file_name)` ¶

class to read Radial Distribution Fuction (RDF) file from Lammps compute

Attributes:

file_name (str) –

file name.
frame (DataFrame) –

3d pandas Frame (multi-row-index DataFrame).

Methods:

ReadRDF –

read RDF file
AverageRDF –

the Average RDF

initiate object

Parameters:

file_name (str) –

file_name

Returns:

Obj ( LmpAveChunk ) –

LmpAveChunk object

Examples:

RDF = LmpAveChunk('LmpAveChunk.txt')

Methods:

read_AveChunk –

Args:
compute_AveChunk –

compute average of RDF over all frames

`read_AveChunk(file_name)` ¶

Parameters:

file_name (str) –

input RDF file

Returns:

Obj ( LmpAveChunk ) –

LmpAveChunk object

`compute_AveChunk()` ¶

compute average of RDF over all frames

Returns:

df ( DataFrame ) –

Average of RDF

`PlumHistogram(file_name)` ¶

Create an Object of DUMP file

Methods:

read_histogram –

read Histogram file
average_histogram –

the Average Histogram
AreaHisto –

Area under pdf curve
find_tail –

find limit of histogram
find_center –

find center of histogram

Examples:

from thmd.io.read_block import PlumHistogram
RDF = PlumHistogram(file_name='myRDF.txt')

initiate object

Parameters:

file_name (str) –

file_name

Returns:

Obj ( PlumHistogram ) –

PlumHistogram object

Methods:

read_histogram –

Args:
compute_average_histogram –

compute average of histogram over all frames
areaHisto –
fit_std_gaussian –

Fit the average-histogarm to Standard Gaussian function
find_tail –

Find tail of distribution function
find_center –

Find tail of distribution function

`read_histogram(file_name: str)` ¶

Parameters:

file_name (str) –

input HISTOGRAM file

Returns:

Obj ( PlumHistogram ) –

update PlumHistogram object

`compute_average_histogram()` ¶

compute average of histogram over all frames

Returns:

df ( DataFrame ) –

DataFrame of avergave histogram

`areaHisto()` ¶

`fit_std_gaussian()` ¶

Fit the average-histogarm to Standard Gaussian function

Returns:

amp, miu, sigma) (tuple –

parameters of Gaussian function

`find_tail(tol=0.0001, gridSize=1e-06)` ¶

Find tail of distribution function

Parameters:

tol (float, default: 0.0001 ) –

tolerance
gridSize (float, default: 1e-06 ) –

size of grid

Returns:

left_tail ( float ) –

limit on the left side
right_tail ( float ) –

limit on the right side

`find_center(gridSize=1e-06)` ¶

Find tail of distribution function

Parameters:

gridSize (float, default: 1e-06 ) –

size of grid

Returns:

Xcenter ( float ) –

center of the distribution function

`average_df_list(list_df: list[pl.DataFrame]) -> pl.DataFrame` ¶

compute average of list of DataFrame

Parameters:

list_df (list) –

list of DataFrame

Returns:

df ( DataFrame ) –

Average of list of DataFrame

`read_data` ¶

This module contains functions to read numeric data from various formats of TEXT files.

Functions:

matrix_lost –

Function to read data in matrix form, in which number of values in each line are NOT equal (missing values)
matrix –

Function to read Data that is as a regular matrix.
logMFD –

Function to read data from LogMFD calculation.
lammps_var –

Function to extract variable values from LAMMPS input file.
plumed_var –

Function to extract variable values from PLUMED input file.
list_matrix_in_dir –

read data from all *.txt files in current and sub-folders.

`matrix_lost(file_name: str, header_line: int = None, column_names: list[str] = None, comment: str = '#', sep: str = ' ', read_note: bool = False) -> pl.DataFrame` ¶

Function to read data in matrix form, in which number of values in each line are NOT equal (missing values) This cannot be read by Numpy, polars,...

The names of columns are extracted from header_line or set by column_names. If both column_names and header_line are not available, the default column's name is: 0 1 2...

Parameters:

file_name (str) –

the text file.
header_line (int, default: None ) –

the lines to extract column-names. Defaults to None.
column_names (list, default: None ) –

Names of columns to extract. Defaults to None.
comment (str, default: '#' ) –

comment-line mark. Defaults to "#".
sep (str, default: ' ' ) –

separator. Defaults to " ".
read_note (bool, default: False ) –

read 'note' column (any text beyond comment mark). Defaults to False.

Returns:

df ( DataFrame ) –

polars DataFrame

Notes

To return 2 lists from list comprehension, it is better (may faster) running 2 separated list comprehensions.
.strip() function removes trailing and leading space in string.

`matrix(file_name: str, header_line: int = None, column_names: list[str] = None, usecols: tuple[int] = None) -> pl.DataFrame` ¶

Function to read Data that is as a regular matrix. The names of columns are exatract based on column_names or header_line. If both column_names and header_line are not available, the default column's name is: 0 1 2...

Parameters:

file_name (str) –

the text file.
header_line (int, default: None ) –

the line to extract column-names. Defaults to None.
column_names (list[str], default: None ) –

Names of columns to extract. Defaults to None.
usecols (tuple[int], default: None ) –

only extract some columns. Defaults to None.

Returns:

df ( DataFrame ) –

polars DataFrame

`logMFD(file_name, dim=1) -> pl.DataFrame` ¶

Function to read data from LogMFD calculation.

Parameters:

file_name (str) –

the logmfd.out file.
dim (int, default: 1 ) –

dimension of LogMFD calulation. Defaults to 1.

Raises:

Exception –

description

Returns:

df ( DataFrame ) –

polars DataFrame

`lammps_var(file_name, var_names=None)` ¶

Function to extract variable values from LAMMPS input file.

Parameters:

file_name (str) –

the text file in LAMMPS input format.
var_names (list, default: None ) –

list of varibalbes to be extracted. Default to None. mean extract all variables.

Returns:

df ( DataFrame ) –

polars DataFrame contains variable in Lammps file

`plumed_var(file_name, var_name, block_name=None)` ¶

Function to extract variable values from PLUMED input file.

Parameters:

file_name (str) –

the text file in LAMMPS input format.
var_name (str) –

list of keyworks in PLUMED, ex: INTERVAL,...
block_name (str, default: None ) –

block command in Plumed, ex: METAD, LOGMFD. Defaults to None.

Returns:

value ( float ) –

value of plumed_var.

Refs

Include negative decimal numbers in regular expression

`list_matrix_in_dir(search_key='deform_', file_ext='.txt', read_note=False, recursive=True)` ¶

read data from all *.txt files in current and sub-folders.

Parameters:

search_key (str, default: 'deform_' ) –

a string to search file_name.
file_ext (str, default: '.txt' ) –

file extension. Default to '.txt'
read_note (bool, default: False ) –

read 'note' column in pl.DataFrame. Default to False.
recursive (bool, default: True ) –

search in sub-folders. Default to True.

Returns:

ldf ( list ) –

list of DataFrames.
files ( list ) –

list of filenames.

`script` ¶

This module contains functions to read/write some specific scripts.

Functions:

read_lines –

Function to read lines in a script that match some KEY_WORDs.
read_plumed_block –

Function to read block_command in PLUMED script.
write_lines –

Funtion to write a list of strings into file.
write_list –

Funtion to write a list of strings into file.

`read_lines(file_name: str, keywords: list = [])` ¶

Function to read lines in a script that match some KEY_WORDs.

Parameters:

file_name (str) –

a text file of any format.
keywords (list, default: [] ) –

list-of-Keywords to extract a line, ex: METAD, LOGMFD. Default to [], mean read all lines.

Returns:

lines ( list ) –

a list of lines.
rest_lines ( list ) –

a list of rest lines.

`read_plumed_block(file_name: str, block_name: str = ' ') -> list` ¶

Function to read block_command in PLUMED script.

Parameters:

file_name (str) –

a text file of PLUMED format.
block_name (str, default: ' ' ) –

block command in PLUMED, ex: METAD, LOGMFD

Returns:

lines ( list ) –

block_of_commandlines

`write_lines(filename: str, lines: list)` ¶

Funtion to write a list of strings into file.

Parameters:

filename (str) –

file name.
lines (list) –

list of strings.

`write_list(lines: list, filename: str)` ¶

Funtion to write a list of strings into file.

Parameters:

filename (str) –

file name.
lines (list) –

list of strings.

`traj` ¶

This module contains classes and functions to process molecular dynamics trajectories.

Modules:

md_traj –
mod_calc_global_properties –
mod_calc_local_properties –
mod_change_box –
mod_change_data –
parser_gromacs –
parser_pdb_xyz –

`md_traj` ¶

Classes:

Frame –

Create an Object for a single-FRAME of trajectories from MD simulation.
Traj –

Create an Object for a multi-FRAMEs of trajectories from MD simulation.

`Frame(dump_file: str = None, data_file: str = None, atom_style: str = 'auto', pdb_file: str = None, xyz_file: str = None, from_df: pd.DataFrame = None, box: np.array = None, box_angle: np.array = None)` ¶

Create an Object for a single-FRAME of trajectories from MD simulation.

This class create a data-object (single configuration) for the analysis of computing data from LAMMPS. The file formats implemented in this class

This class implemented several ways to create Frame object

create an empty data object
create_DATA object with input data
read from DUMP file
read from DATA file
read frome PDB file

Attributes:

filename (str) –

name of input file
timestep (int) –

the timestep of configuration
box (array) –

3x2 array, the box size
box_angle (array) –

1x3 array, the box angle
atom (DataFrame) –

DataFrame of per-atom values
prop_key (list) –

column-names of properties
mass (DataFrame) –

DataFrame of per-type masses
fmt (str) –

default format for float numbers, don't use %g because it will lost precision

Examples:

from thmd.traj  import Frame

da = Frame()                        # empty object
da = Frame(from_df=df)              # oject with input data
da = Frame(dump_file='test.cfg')    # from DUMP file
da = Frame(data_file='mydata.dat')  # from DATA file
da = Frame(pdb_file='test.pdb')     # from PDB file

Refs: [1]. Use chain mutator calls

initilize the Frame object

Parameters:

dump_file (str, default: None ) –

filename of DUMP file.
data_file (str, default: None ) –

filename of DATA file.
pdb_file (str, default: None ) –

filename of PBD file.
xyz_file (str, default: None ) –

filename of XYZ file.
from_df (DataFrame, default: None ) –

create FRAME from data.
atom_style (str, default: 'auto' ) –

atom_style of system. Only need when data_file is used. Possible values: 'atomic', 'molecular', 'charge', 'full', 'auto'
box (np.array list, default: None ) –

Define simulation box. Only need when from_df is used.
box_angle (np.array list, default: None ) –

Define angle of simulation box. Only need when from_df is used.

Returns:

Obj ( Frame ) –

object of trajectories

Notes

Use mutator, so do not use self.* when define value

Methods:

create_DATA –

The method to create new FRAME object with input data.
copy –

The method to make an indepedent copy of Frame Obj. Then, the change values of the fields of the new object, the old object should not be affected by that.
check_exist –

The method to check whether something is existed in system or not.

Attributes:

filename (str) –
timestep (int) –
decimals (int) –
box –
box_angle –
prop_key (list[str]) –
atom –
mass –
pair_coeff –
bond_coeff –
angle_coeff –
dihedral_coeff –
improper_coeff –
bondBond_coeff –
bondAngle_coeff –
middleBondTorsion_coeff –
endBondTorsion_coeff –
angleTorsion_coeff –
angleAngleTorsion_coeff –
bondBond13_coeff –
angleAngle_coeff –
bond –
angle –
dihedral –
improper –

filename: str = 'md_frame' instance-attribute ¶

timestep: int = None instance-attribute ¶

decimals: int = 7 instance-attribute ¶

box = np.asarray([[0, 1], [0, 1], [0, 1]], dtype=float) instance-attribute ¶

box_angle = np.asarray([0, 0, 0], dtype=float) instance-attribute ¶

prop_key: list[str] = None instance-attribute ¶

_num = {'n_atoms': 0, 'n_atom_types': 0} instance-attribute ¶

_style = {'atom_style': '', 'pair_style': '', 'bond_style': '', 'angle_style': '', 'dihedral_style': '', 'improper_style': ''}

instance-attribute ¶

atom = None instance-attribute ¶

mass = None instance-attribute ¶

pair_coeff = None instance-attribute ¶

bond_coeff = None instance-attribute ¶

angle_coeff = None instance-attribute ¶

dihedral_coeff = None instance-attribute ¶

improper_coeff = None instance-attribute ¶

bondBond_coeff = None instance-attribute ¶

bondAngle_coeff = None instance-attribute ¶

middleBondTorsion_coeff = None instance-attribute ¶

endBondTorsion_coeff = None instance-attribute ¶

angleTorsion_coeff = None instance-attribute ¶

angleAngleTorsion_coeff = None instance-attribute ¶

bondBond13_coeff = None instance-attribute ¶

angleAngle_coeff = None instance-attribute ¶

bond = None instance-attribute ¶

angle = None instance-attribute ¶

dihedral = None instance-attribute ¶

improper = None instance-attribute ¶

create_DATA(df: pd.DataFrame, box: np.array = None, box_angle: np.array = None) ¶

The method to create new FRAME object with input data.

Parameters:

DataFrame (DataFrame) –

of input data
box (array, default: None ) –

3x2 array, option to input boxSize.
box_angle (array, default: None ) –

1x3 array, option to input box_angle.

Returns:

Obj ( Frame ) –

update Frame

Examples:

da = Frame()
da.create_DATA(DataFrame=df)
# or
da = Frame(from_df=df)

copy() ¶

The method to make an indepedent copy of Frame Obj. Then, the change values of the fields of the new object, the old object should not be affected by that.

Returns:

Obj ( Frame ) –

new Frame Obj.

Examples:

da1 = da.copy()

Refs: [1]. "shallow copying" vs "deep copying": https://stackoverflow.com/questions/3975376/understanding-dict-copy-shallow-or-deep/3975388#3975388

check_exist(atom_types=None, mass_types=None) ¶

The method to check whether something is existed in system or not.

Parameters:

atom_types (list, default: None ) –

list-of-int of atom-types. Default to None.
mass_types (list, default: None ) –

list-of-int of atom-types. Default to None.

Returns:

mgs ( str ) –

raise Message if error.

Examples:

da.isExist(atom_types=[2,3])

Notes

set() also return unique values.

`Traj(**kwargs)` ¶

Create an Object for a multi-FRAMEs of trajectories from MD simulation. - read frome XYZ file

initilize the TrajFrame object

Notes

Use mutator, so do not use self.* when define value

Methods:

readXYZ –

The method create Multi-FRAME object by reading XYZ file.

Attributes:

decimals (int) –

decimals: int = 6 instance-attribute ¶

readXYZ(filename) ¶

The method create Multi-FRAME object by reading XYZ file.

Parameters:

filename (str) –

name of input file

Returns:

Obj ( TrajFrame ) –

update FRAME

Examples:

da = io.TrajFrame(pdb_file='mydata.pdb')

`mod_calc_global_properties` ¶

Functions:

compute_mass –

The method to compute mass of selected atom_types.
compute_wt_percent –

The method to compute weight percentage of some atom_types.

`compute_mass(self, atom_types=[])` ¶

The method to compute mass of selected atom_types.

Parameters:

atom_types (list, default: [] ) –

atom-types to compute masses. Defaults to [], mean all-types.

Returns:

m ( float ) –

total mass of selected atoms.

Examples:

da.compute_mass(atom_types=[2,3])

`compute_wt_percent(self, atom_types)` ¶

The method to compute weight percentage of some atom_types.

Parameters:

atom_types (list) –

atom-types compute percentage of weight.

Returns:

wt ( float ) –

weight percentage of chosen atoms.

Examples:

da.compute_wt_percent(atom_types=[2,3])

`mod_calc_local_properties` ¶

Functions:

compute_colvar_sph_harm –

Compute per-atom vector of spherical harmonics

`compute_colvar_sph_harm(self, l, kind='real', normalization='4pi', bound_cond=[1, 1, 1], cutoff_neighbor=None, k_neighbor=None, k_cutoff=20, keep_ref=False, calc_coord=False)` ¶

Compute per-atom vector of spherical harmonics

Parameters:

l (int) ) –

degree of Spherical Harmonic
form ( (str) –

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

`mod_change_box` ¶

Functions:

combine_frame –

The method to combine 2 Lammps Frames.
unwrap_coord_DATA –

The method to upwrap coords in DATA file.
wrap_coords_DUMP –

The method to flip coords over the center.
flip_coords –

The method to flip coords over the center.
replicate –

The method to flip coords over the center.
scale_box –

The method to change size of simulation box.

`combine_frame(self, TrajFrame, merge_type=False, alignment='comXYZ', shift_XYZ=None, separate_XYZ=None, merge_box=True, use_box='box1')` ¶

The method to combine 2 Lammps Frames.

Parameters:

TrajFrame (TrajFrame Obj) –

an Object of TrajFrame
merge_type (bool, default: False ) –

merge the same type in 2 TrajFrame.
alignment (str, default: 'comXYZ' ) –

choose how to align 2 frame. Defaults to 'comXYZ'. + 'comXYZ': align based on COM + 'minXYZ': align based on left corner + 'maxXYZ': align based on right corner
shift_XYZ (list, default: None ) –

shift a distance from COM aligment. Defaults to [0,0,0].
separate_XYZ (list, default: None ) –

Separate 2 frame with a specific value. Defaults to [0,0,0].
merge_box (bool, default: True ) –

choose to merge box or not. Defaults to True.
use_box (str, default: 'box1' ) –

be used as the box size if merge_box=False. Defaults to 'box1'.

Returns:

Obj ( TrajFrame ) –

Update TrajFrame da1

Examples:

da1.combine_frame(da2)

Todo

combine box_angle

Refs:

[1]. Deep copy: https://stackoverflow.com/questions/3975376/understanding-dict-copy-shallow-or-deep/3975388#3975388

`unwrap_coord_DATA(self, imgFlag=['x', 'y', 'z'], atom_types=[])` ¶

The method to upwrap coords in DATA file.

Parameters:

imgFlag (list, default: ['x', 'y', 'z'] ) –

image Flags in data file. Defaults to ['x','y','z'].
atom_types (list, default: [] ) –

just unwrap some atom-types. Defaults to [], mean unwrap all-types.

Returns:

Obj ( TrajFrame ) –

update FRAME

Notes

cannot unwrap_coord_data if imgFlags are not available.

`wrap_coords_DUMP(self, dim=[1, 1, 1], inplace=True)` ¶

The method to flip coords over the center.

Parameters:

dim (list, default: [1, 1, 1] ) –

choose the dimenstion to take flip. Defaults to [1,1,1].
inplace (bool, default: True ) –

if True it will replace (or addd x y z columns). If False it will return xyz array.

Returns:

Obj ( TrajFrame ) –

update FRAME

`flip_coords(self, dim=[1, 1, 1])` ¶

The method to flip coords over the center.

Parameters:

dim (list, default: [1, 1, 1] ) –

choose the dimenstion to take flip. Defaults to [1,1,1].

Returns:

Obj ( TrajFrame ) –

update FRAME

TODOs

Remove pandas Warning.

`replicate(self, dim=[1, 1, 1])` ¶

The method to flip coords over the center.

Parameters:

dim (list, default: [1, 1, 1] ) –

choose the dimenstion to take flip. Defaults to [1,1,1].

Returns:

Obj ( TrajFrame ) –

update FRAME

`scale_box(self, scale=None, final=None, delta=None, remap=True)` ¶

The method to change size of simulation box.

Parameters:

scale (list, default: None ) –

to set scale ratio on each dimension of the box. scale = [0.7, 0.7, None] : if one dimension is set "None" its length is not changed.
final (list, default: None ) –

to set final length on each dimension of the box.
delta (list, default: None ) –

to set amount of change on each dimension of the box.
remap (bool, default: True ) –

remap atom coordinate. Default to True.

Returns:

Obj ( TrajFrame ) –

update FRAME

Examples:

da.scale_box(scale=[0.7, 0.7, None])

`mod_change_data` ¶

Functions:

add_column –

The method to add new columns to da.atom.
delete_column –

The method to delete columns from da.atom.
set_mass –

The method to set masses of atoms in system. Before use it, need to define element_dict with 2 keys: 'type', 'atom_symbol'
change_atom_type –

The method to change types of atoms in system.
merge_atom_type –

The method to merge types of atoms in system.

`add_column(self, data, newColumn=None, replace=False)` ¶

The method to add new columns to da.atom.

Parameters:

data (pd.DataFrame pd.Series list) –

Nxm data of new columns
newColumn (list, default: None ) –

1xN list contains names of columns. Default to None, mean it will take columnNames from DataFrame
replace (bool, default: False ) –

replace column if existed.

Returns:

Obj ( TrajFrame ) –

Update da.atom

Examples:

da.add_column(df, myColumn=['col1','col2'], replace=True)

`delete_column(self, delColumns)` ¶

The method to delete columns from da.atom.

Parameters:

delColumns (list) –

1xN list contains names of columns to be deleted.

Returns:

Obj ( TrajFrame ) –

Update da.atom

Examples:

da.delete_column(delColumns=['col1','col2'])

`set_mass(self, element_dict)` ¶

The method to set masses of atoms in system. Before use it, need to define element_dict with 2 keys: 'type', 'atom_symbol' element_dict={'type': list_values, 'atom_symbol':list_values}

Parameters:

element_dict (dict) –

a dict to define atom-types and atom-symbols.

Returns:

Obj ( TrajFrame ) –

Update da.atom

Examples:

da.set_mass(element_dict={'type':[1,2,3], 'atom_symbol':['C','H','N']})

`change_atom_type(self, old_type, new_type, save_old_type=True)` ¶

The method to change types of atoms in system.

Parameters:

old_type (list) –

a list of old-types.
new_type (int) –

one new-type.
save_old_type (bool, default: True ) –

to back up old types. Default to True.

Returns:

Obj ( TrajFrame ) –

update FRAME

Examples:

da.chage_atom_type([1,2,3], 2)

`merge_atom_type(self, old_type, save_old_type=True)` ¶

The method to merge types of atoms in system.

Parameters:

old_type (list) –

a list of old-types.
save_old_type (bool, default: True ) –

to back up old types. Default to True.

Returns:

Obj ( TrajFrame ) –

update FRAME

Examples:

da.chage_atom_type([1,2,3], 2)

`parser_gromacs` ¶

`parser_pdb_xyz` ¶

Functions:

read_PDB –

The method to create FRAME object by reading PDB file.
read_XYZ –

The method to create FRAME object by reading XYZ file.
write_XYZ –

The method to write XYZ file.
write_PDB –

The method to write PDB file

`read_PDB(self, filename)` ¶

The method to create FRAME object by reading PDB file.

Parameters:

filename (str) –

name of input file

Returns:

Obj ( TrajFrame ) –

update FRAME record_name (str): atom_symbol (str): same as column 'type' in DUMP format residue_name (str): residue_id (int): chain (str): occupancy (float): beta (float):

Examples:

da = TrajFrame(pdb_file='mydata.pdb')

`read_XYZ(self, filename)` ¶

The method to create FRAME object by reading XYZ file.

Parameters:

filename (str) –

name of input file

Returns:

Obj ( TrajFrame ) –

update FRAME

Examples:

da = TrajFrame(pdb_file='mydata.pdb')

`write_XYZ(self, filename, column=['X', 'xu', 'yu', 'zu'], fmt=None)` ¶

The method to write XYZ file.

Parameters:

filename (str) –

name of input file
column (list, default: ['X', 'xu', 'yu', 'zu'] ) –

list-of-str contains columns to be written. Defaults to ['X','xu','yu','zu']
fmt (str, default: None ) –

string format for output values. Defaults to None, mean use self.fmt

Returns:

file ( obj ) –

the XYZ file

Examples:

da.write_XYZ('test.xyz')

`write_PDB(self, filename, writeBox=False)` ¶

The method to write PDB file

Parameters:

filename (str) –

name of input file
writeBox (bool, default: False ) –

write box or not.

Returns:

file ( obj ) –

the PDB file

Examples:

da.write_PDB('test.pdb')

`latex` ¶

Functions:

get_citekeys –

Extract all citekeys from a list of .tex files.
select_entries –

Select bibliography entries from a Bibtex .bib file
write_bibtex –

Write bibliography data to a Bibtex .bib file
minimize_bibtex –

Minimize the .bib file so that it contains only entries that are cited in the .tex files.
check_bibtex –

Check for missing fields in the bibtex entries. Based on the type of the entry, there are mandatory fields. This function is to check if any of mandatory fields is missing.
tex2doc –

Convert a Latex source .tex file to a Microsoft Word .docx file.
tex2pdf –

Convert a Latex source .tex file to a Microsoft Word .docx file.
replace_acronyms –

Check and replace if acronyms are used in a .tex file.
count_cited_journals –

Count the number of times each journal is cited in the .bib file.

`_DATA_PATH = Path(file).parents[0].resolve() / '_data'` `module-attribute` ¶

`_check_pybtex() -> None` ¶

Check if the pybtex library is installed or not. If not, raise an error.

`_check_pandoc() -> None` ¶

check if pandoc is installed or not. If not, raise an error.

`get_citekeys(tex_files: Union[str, list]) -> list` ¶

Extract all citekeys from a list of .tex files.

Parameters:

tex_files (Union[str, list]) –

paths to .tex files

Returns:

citekeys ( list ) –

list of citekeys

`select_entries(citekeys: list, bib_file: str) -> BibliographyData` ¶

Select bibliography entries from a Bibtex .bib file

Parameters:

citekeys (list) –

list of citekeys
bib_file (str) –

path to a .bib file

Returns:

BibliographyData ( BibliographyData ) –

selected bibliography entries. See pybtex.database.BibliographyData

`write_bibtex(bib_data: BibliographyData, out_file: str = 'reference.bib') -> None` ¶

Write bibliography data to a Bibtex .bib file

Parameters:

bib_data (BibliographyData) –

bibliography entries. See pybtex.database.BibliographyData
out_file (str, default: 'reference.bib' ) –

path to the output .bib file

`minimize_bibtex(tex_files: Union[str, list], bib_file: str, out_file: str = 'reference_min.bib') -> None` ¶

Minimize the .bib file so that it contains only entries that are cited in the .tex files.

Parameters:

tex_files (Union[str, list]) –

path to .tex files
bib_file (str) –

path to the original .bib file
out_file (str, default: 'reference_min.bib' ) –

path to the new .bib file

`_check_missing_fields(entry: object, fields: list) -> list` ¶

`check_bibtex(bib_file: str, verbose: bool = True, logfile: bool = False, lower_key=False) -> list` ¶

Check for missing fields in the bibtex entries. Based on the type of the entry, there are mandatory fields. This function is to check if any of mandatory fields is missing.

Parameters:

bib_file (str) –

path to the original .bib file.
verbose (bool, default: True ) –

print output to the console.
logfile (bool, default: False ) –

write output to a log file.
lower_key (bool, default: False ) –

convert the keys in bib_file to lower case.

Examples:

from thml.latex import check_bibtex
check_bibtex('reference.bib', verbose=True)

`tex2doc(tex_file: str, out_file: str = 'output.docx', bib_file: str = None, **kwargs)` ¶

Convert a Latex source .tex file to a Microsoft Word .docx file. Need to install pandoc and panflute first.

condac install -c conda-forge pandoc
pip install panflute

Parameters:

tex_file (str) –

path to the .tex file
out_file (str, default: 'output.docx' ) –

path to the output .docx file
bib_file (str, default: None ) –

path to the .bib file

Other Parameters:

cite_style (str = 'elsevier_vancouver.csl') –

path to the citation style .cls' file. Can download from [here](https://github.com/citation-style-language/styles) or inZotero_folder/styles`
reference_doc (str) –

path to the reference .docx file
resource_path (str = './figure') –

path to the folder containing the figures and other resources.
verbose (bool = False) –

print the output of the command
pandoc-crossref (bool = True) –

use the pandoc-crossref extension for equation numbering

`tex2pdf(tex_file: str, out_file: str = 'output.pdf', bib_file: str = None, **kwargs)` ¶

Convert a Latex source .tex file to a Microsoft Word .docx file.

Parameters:

tex_file (str) –

path to the .tex file
out_file (str, default: 'output.pdf' ) –

path to the output .docx file
bib_file (str, default: None ) –

path to the .bib file

Other Parameters:

cite_style (str = 'elsevier_vancouver.csl') –

path to the citation style `.cls' file. Can download from here
resource_path (str = './figure') –

path to the folder containing the figures and other resources.
verbose (bool = False) –

print the output of the command

`replace_acronyms(tex_file: str, acronyms: dict) -> list` ¶

Check and replace if acronyms are used in a .tex file.

Not implemented yet.

`count_cited_journals(bib_file: str) -> dict` ¶

Count the number of times each journal is cited in the .bib file.

Parameters:

bib_file (str) –

path to the .bib file

Returns:

dict ( dict ) –

a dictionary containing the journal names and the number of times they are cited.

`model` ¶

Modules:

D1tube –
D2haxagonal –
D3crystal –
box_tool –
combining_LJ_interface –
coord_rotation –
forcefield_info –

This module contains some data for various ForceField. Data obtained from simulation
polymer_mbuild –

This module contains classes and functions to build models of atomic polymers
polymer_pysimm –

This module contains classes to build models of atomic polymers

`D1tube` ¶

Functions:

lattice_CNT –

Calculates the 3D Cartesian coordinates of atoms of of (n,m) CNT.

`_UnitCell_CNT(m, n, bond_CC=1.421, basis_atom='AB', diameter=None)` ¶

Calculates the 3D Cartesian coordinates of atoms of 1 units cell of (n,m) CNT, which which n >= m >= 0

thangckt, Aug 2022

Parameters:

n,m (int) –

Chiral indices n>=m>=0
bond_CC (float, default: 1.421 ) –

Length of C-C bonds
basis_atom (str, default: 'AB' ) –

'AB'/'A'/'B': to create semi-Graphene lattice for adoptting Hydrogen atoms - 'AB': full Garaphene-like crystal - 'A': semi Graphene-like with atom at A-position - 'B': semi Graphene-like with atom at B-position

Returns:

df ( DataFrame ) –

Nx3 array, contain positions of atoms of 1 units cell of (n,m)graphene sheet.
unitbox ( array ) –

size of unit box
param ( dict ) –

dict contains characteristic parameters of graphene lattice. 'Chiral_len' (float): length of Translational unit vector, corresponding to length on Y direction 'Translate_len' (float): length of Chiral vector, corresponding to length on X direction 'Chiral_ang' (float): Chiral angle of Graphene sheet in Degree. 'Chiral_vector' (array_1x2): Chiral vector 'Translate_vector' (array_1x2): translations vector

Notes

n=1, m=0 : is the unit cell for Zigzag
n=1, m=1 : is the unit cell for Armchair
n>m      : unit cell is automatically computed

`lattice_CNT(m, n, bond_CC=1.421, aspect=1, basis_atom='AB')` ¶

Calculates the 3D Cartesian coordinates of atoms of of (n,m) CNT.

thangckt, Aug 2022

Parameters:

n,m (int) –

Chiral indices n>=m>=0
bond_CC (float, default: 1.421 ) –

Length of C-C bonds
aspect (int, default: 1 ) –

The nanotube aspect ratio L/D.
basis_atom (str, default: 'AB' ) –

'AB'/'A'/'B': to create semi-Graphene lattice for adoptting Hydrogen atoms - 'AB': full Garaphene-like crystal - 'A': semi Graphene-like with atom at A-position - 'B': semi Graphene-like with atom at B-position

Returns:

df ( DataFrame ) –

Nx3 array, contain positions of atoms of 1 units cell of (n,m)graphene sheet.
box ( array ) –

Simulation box
param ( dict ) –

dict contains characteristic parameters of graphene lattice. 'Chiral_len' (float): length of Translational unit vector, corresponding to length on Y direction 'Translate_len' (float): length of Chiral vector, corresponding to length on X direction 'Chiral_ang' (float): Chiral angle of Graphene sheet in Degree. 'Chiral_vector' (array): Chiral vector 'Translate_vector' (array): translations vector

Notes

n=1, m=0 : is the unit cell for Zigzag
n=1, m=1 : is the unit cell for Armchair
n>m      : unit cell is automatically computed

`D2haxagonal` ¶

Functions:

lattice_Graphene –

Calculates the 3D Cartesian coordinates of atoms of (n,m)graphene sheet/ Graphite

`_UnitCell_2Dhoneycomb(m, n, bond_CC=1.421, basis_atom='AB')` ¶

Calculates the 3D Cartesian coordinates of atoms of 1 units cell of (n,m)graphene sheet, which which n >= m >= 0

thangckt, Nov 2019 (update 2022)

Parameters:

n,m ( (int) –

Chiral indices n>=m>=0
bond_CC (float, default: 1.421 ) –

Length of C-C bonds
basis_atom (str, default: 'AB' ) –

'AB'/'A'/'B': to create semi-Graphene lattice for adoptting Hydrogen atoms - 'AB': full Garaphene-like crystal - 'A': semi Graphene-like with atom at A-position - 'B': semi Graphene-like with atom at B-position

Returns:

–

df (DataFrame) : Nx3 array, contain positions of atoms of 1 units cell of (n,m)graphene sheet.
–

unitbox (array) : size of unit box
–

param (dict) : dict contains characteristic parameters of graphene lattice. 'Chiral_len' (float) : length of Translational unit vector, corresponding to length on Y direction 'Translate_len' (float) : length of Chiral vector, corresponding to length on X direction 'Chiral_ang' (float) : Chiral angle of Graphene sheet in Degree. 'Chiral_vector' (array_1x2) : Chiral vector 'Translate_vector' (array_1x2): translations vector

Notes

n=1, m=0 : is the unit cell for Zigzag
n=1, m=1 : is the unit cell for Armchair
n>m      : unit cell is automatically computed

Refs

Dresselhaus et al. “Physics of Carbon Nanotubes.”, 1995, doi:10.1016/0008-6223(95)00017-8.
Antonsen, and Thomas Garm Pedersen. “Characterisation and Modelling of Carbon Nanotubes,” 2013.

`_UnitCell_Graphene(m, n, bond_CC=1.421, basis_atom='AB')` ¶

`lattice_Graphene(m, n, bond_CC=1.421, sheet_size=[1, 1], sheet_number=1, layer_bond=3.35, basis_atom='AB')` ¶

Calculates the 3D Cartesian coordinates of atoms of (n,m)graphene sheet/ Graphite

thangckt, Nov 2019 (update 2022)

Parameters:

n (int) –

Chiral indices n>=m>=0
m (int) –

Chiral indices n>=m>=0
bond_CC (float, default: 1.421 ) –

Length of C-C bonds
sheet_size (list, default: [1, 1] ) –

[Xsize, Ysize], size of graphene sheet
sheet_number (int, default: 1 ) –

number of sheets
layer_bond (float, default: 3.35 ) –

Length of plane-plane bonds
basis_atom (str, default: 'AB' ) –

'AB'/'A'/'B': to create semi-Graphene lattice for adoptting Hydrogen atoms - 'AB': full Garaphene-like crystal - 'A': semi Graphene-like with atom at A-position - 'B': semi Graphene-like with atom at B-position

Returns:

df ( DataFrame ) –

Nx3 array, contain positions of atoms of 1 units cell of (n,m)graphene sheet.
box ( array ) –

Simulation box
param ( dict ) –

dict contains characteristic parameters of graphene lattice. 'Chiral_len' (float): length of Translational unit vector, corresponding to length on Y direction 'Translate_len' (float): length of Chiral vector, corresponding to length on X direction 'Chiral_ang' (float): Chiral angle of Graphene sheet in Degree. 'Chiral_vector' (array): Chiral vector 'Translate_vector' (array): translations vector

Notes

n=1, m=0 : is the unit cell for Zigzag
n=1, m=1 : is the unit cell for Armchair
n>m      : unit cell is automatically computed

`D3crystal` ¶

Functions:

lattice_orthoRHOMBIC –

Function to create atomic coordinates for crystal model
lattice_CUBIC –

Shortcut to create CUBIC crystal, as subclass of ortthoRHOMBIC

`_UnitCell_orthoRHOMBIC(crystal_type, lattice_constant=[1, 1, 1])` ¶

a DICT, contain 1 conventional cell of a Simple crystal UnitCell(FCC, BCC,...), in pricipal axes [100][010] [001] - V2O5 crystal based on: https://materialsproject.org/materials/mvc-11944/ but change order of lattice constants: c,a,b - This order is the same as: 10.1016/j.triboint.2020.106750

Parameters:

crystal_type (str) –

'V2O5', 'FCC', 'BCC'
lattice_constant (list, default: [1, 1, 1] ) –

lattice constant [a,b,c] corresponding to [x,y,z]

Returns:

R ( array ) –

Nx3 array, contain positions of atoms in conventional unit cell.

`lattice_orthoRHOMBIC(crystal_type, lattice_constant, orient=[[1, 0, 0], [0, 1, 0], [0, 0, 1]], size=[1, 1, 1], bound_cond=[1, 1, 1], tol_on_bound=0.1)` ¶

Function to create atomic coordinates for crystal model

Parameters:

crystal_type (str) –

'V2O5', 'FCC', 'BCC'
lattice_constant (list) –

lattice constant [a,b,c] corresponding to [x,y,z]
orient (list, default: [[1, 0, 0], [0, 1, 0], [0, 0, 1]] ) –

3x3 array, contain direction vectors define crystal orientation, ex: ([[1,0,0], [0,1,0], [0,0,1]])
size (list, default: [1, 1, 1] ) –

[Nx Ny Nz] 1x3 array, size of model, Nx is X-size in lattice constant unit
bound_cond (list, default: [1, 1, 1] ) –

1x3 array contain convention for boundary conditions: 1 is peridic; 0 is not

Returns:

points ( array ) –

Nx3 array contain positions of atoms.
box ( array ) –

3x2 array contain size of box contain lattice ([[xlo, xhi], [ylo, yhi], [zlo, zhi]])
unit_box ( array ) –

3x2 array contain size of unit cell

`lattice_CUBIC(crystal_type, lattice_constant, orient=[[1, 0, 0], [0, 1, 0], [0, 0, 1]], size=[1, 1, 1], bound_cond=[1, 1, 1], tol_on_bound=0.1)` ¶

Shortcut to create CUBIC crystal, as subclass of ortthoRHOMBIC

Parameters:

crystal_type (str) –

'FCC', 'BCC'
lattice_constant (float) –

lattice constant a
orient (list, default: [[1, 0, 0], [0, 1, 0], [0, 0, 1]] ) –

3x3 array, contain direction vectors define crystal orientation, ex: ([[1,0,0], [0,1,0], [0,0,1]])
size (list, default: [1, 1, 1] ) –

[Nx Ny Nz] 1x3 array, size of model, Nx is X-size in lattice constant unit
bound_cond (list, default: [1, 1, 1] ) –

1x3 array contain convention for boundary conditions: 1 is peridic; 0 is not

Returns:

points ( array ) –

Nx3 array contain positions of atoms.
box ( array ) –

3x2 array contain size of box contain lattice ([[xlo, xhi], [ylo, yhi], [zlo, zhi]])
unit_box ( array ) –

3x2 array contain size of unit cell

`box_tool` ¶

Functions:

add_periodic_image –

Function to add "Periodic Images" of atoms at Periodic Boundaries (with a specific cutoff distance)
wrap_coord_PBC –

Function to wrap atom positions at Periodic Boundaries
shell_fcc –

Compute nearest-neighbor shells for FCC crystal
box_generate –

Generate orientation and dimensions of simulation box.
box_orient –

Generate orirentations of simulation box.

`add_periodic_image(points: pd.DataFrame | np.ndarray | list, box: np.ndarray, bound_cond: list = [1, 1, 1], cutoff: float = 6.5) -> pd.DataFrame` ¶

Function to add "Periodic Images" of atoms at Periodic Boundaries (with a specific cutoff distance) By Thang, June 2019 (update 2022)

Parameters:

points (2d-list np.array pd.DataFrame) –

Mx3 Matrix contain positions of atoms before Wrapping
box (3d-list array) –

3x2 Matrix contain simulation box bounds
bound_cond (list, default: [1, 1, 1] ) –

1x3 list contains convention of Peridic bounds(ex: bound_cond = [1 1 1])
cutoff (float, default: 6.5 ) –

Cutoff distance

Returns:

df ( DataFrame ) –

contains original atoms and image atoms with remark colum df['image'].

Examples:

```py

df = add_periodic_image(P, box, bound_cond=[1 1 0], cutoff=5)
    ```

`wrap_coord_PBC(points: pd.DataFrame | np.ndarray | list, box: np.ndarray, bound_cond: list = [1, 1, 1]) -> pd.DataFrame` ¶

Function to wrap atom positions at Periodic Boundaries By Thang, June 2019 (update 2022)

Parameters:

points (2d-list np.array pd.DataFrame) –

Mx3 Matrix contain positions of atoms before Wrapping
box (3d-list array) –

3x2 Matrix contain simulation box bounds
bound_cond (list, default: [1, 1, 1] ) –

1x3 list contains convention of Peridic bounds(ex: bound_cond = [1 1 1])

Returns:

df ( DataFrame ) –

contains atom positions.

Examples: py df = wrap_coord_PBC(P, box, bound_cond=[1 1 0], cutoff=5)

`shell_fcc(a)` ¶

Compute nearest-neighbor shells for FCC crystal

Parameters:

a (float) –

lattice constant

Returns:

shell ( list ) –

5 nearest-neighbor shells

`box_generate(box_size: list = [1, 1, 1], zDirect: str = '001', xDirect: str = None) -> Dict[str, List]` ¶

Generate orientation and dimensions of simulation box.

Parameters:

box_size (list, default: [1, 1, 1] ) –

dimension of box on each side as in [100] direction
zDirect (str, default: '001' ) –

specify the direction of z-side. Defaults to '001', mean nothing is happen.
xDirect (str, default: None ) –

specify the direction of z-side. Defaults to None.

Returns:

box ( dict ) –

dictionary contain orientation and box_size - 'orient' (list[list]): list-of-vectors of 3 directional vectors. - 'box_size' (list[list]): dimension of box on each side.

Examples:

box = box_generate(box_size=[1, 1, 1], zDirect='001')
orient = box['orient']
box_size = box['box_size']

`box_orient(zDirect: str = '001', xDirect: str = None) -> List[List]` ¶

Generate orirentations of simulation box. Args: zDirect (str): specify the direction of z-side. Defaults to '001', mean nothing is happen. xDirect (str): specify the direction of z-side. Defaults to None.

Returns:

orient ( list[list] ) –

list-of-vectors of 3 directional vectors.

`combining_LJ_interface` ¶

Functions:

pair_LJ –

compute parameters (epsilon & sigma) of LJ potential at interface

`pair_LJ(dict_group1, dict_group2, unit_style, combining_rule='geometric', pair_style='lj/cut')` ¶

compute parameters (epsilon & sigma) of LJ potential at interface Note that in LAMMPS: 'Lorentz_Berthelot'='arithmetic' https://tinyurl.com/yzpwg2hs

Parameters:

dict_group1, (dict_group2) –

Dicts contain sig & eps of each element of 2 surfaces. Must contain keys: 'atom_name', 'type', 'sigma', 'epsilon'
unit_style –

'real' or 'metal'

combining_rule='arithmetic' (also 'Lorentz_Berthelot')
            + 'arithmetic'/'Lorentz_Berthelot'
            + 'geometric'
            + 'sixthpower'
        pair_style='lj/cut': pair_style of Lammps  lj/cut/coul/long
        external_interaction: require

Return

list-of-string: contain pair_coeffs for LAMMPS

Notes

energy unit is kcal/mol, but in OPLSaa of Foyer is kJ/mol. types in 2 dict must either completely different or indentical

Examples: PMMA/h_BN interface

dict_group1 = {'element':['CT','CT','CT','CT','HC','HC','C_2','O_2','OS','CT','HC'],
        'atom_name':['opls_135','opls_136','opls_137','opls_139','opls_140','opls_282','opls_465','opls_466','opls_467','opls_468','opls_469'],
        'type':[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11],
        'sigma': [3.5, 3.5, 3.5, 3.5, 2.5, 2.42, 3.75, 2.96, 3.0, 3.5, 2.42],
        'epsilon':[0.066, 0.066, 0.066, 0.066, 0.03, 0.015, 0.105, 0.21, 0.17, 0.066, 0.015]}
dict_group2 = {'element':['B','N'],
                'atom_name':['B','N'],
                'type':[12,13],
                'sigma': [3.453, 3.365],
                'epsilon':[0.094988, 0.1448671]}
combining_LJ(dict_group1, dict_group2, combining_rule='Lorentz_Berthelot', pair_style='lj/cut/coul/long')

`_lorentz(eps1, eps2, sig1, sig2)` ¶

`_geometric(eps1, eps2, sig1, sig2)` ¶

`_sixthpower(eps1, eps2, sig1, sig2)` ¶

`_compute_LJ_param(combining_rule, eps1, eps2, sig1, sig2)` ¶

`coord_rotation` ¶

Classes:

CoordTransform –

We can express a rotation using direction-cosines-matrix (DCM) or Euler-angles (phi,theta,psi)

Functions:

rot1axis –

Rotate array of points about 1 axis
check_right_hand –

check right_hand_rule orthogonal of 3 vectors
guess_right_hand –

give 2 vectors, then guess the third vector that satisfy right_hand_rule
cartesian2spherical –

Convert cartesian coordinates to Spherical coordinates

`CoordTransform(old_orient=[[1, 0, 0], [0, 1, 0], [0, 0, 1]], new_orient=[[1, 0, 0], [0, 1, 0], [0, 0, 1]])` ¶

We can express a rotation using direction-cosines-matrix (DCM) or Euler-angles (phi,theta,psi)

Parameters:

old_orient (array / list, default: [[1, 0, 0], [0, 1, 0], [0, 0, 1]] ) –

3x3 array/list, contains 3 mutully orthotropic unit vectors of the OLD basis
new_orient (array / list, default: [[1, 0, 0], [0, 1, 0], [0, 0, 1]] ) –

3x3 array/list, contains 3 mutully orthotropic unit vectors of the NEW basis (all input vector will be normalized to unit vectors)

Returns:

obj ( obj ) –

3x3 array, the rotation matrix or matrix of direction cosines

Examples:

oldAxis = array([[1, 0, 0], [0, 1, 0], [0, 0, 1]])
newAxis = array([[1, -1, 0], [1, 1, -2], [1, 1, 1]])
BT = thmd.CoordTransform(old_orient=oldAxis, new_orient=newAxis)

"Refs"

Bower, Allan F. Applied Mechanics of Solids. CRC Press, 2009. page 711
https://link.aps.org/doi/10.1103/PhysRevB.92.180102
https://en.wikipedia.org/wiki/Euler_angles

Methods:

direction_cosine_matrix –

Calculate direction-cosines-matrix (DCM) between 2 coordinates systems
rotation_matrix –

Calculate Rotation-matrix (R) as transpose of DCM
EA2ROT –

Calculate Rotation_Matrix Euler Angles (EA) between 2 coordinates systems (ZXZ proper Euler angles)
euler_angle –

Calculate Euler Angles (EA) between 2 coordinates systems (intrinsic ZXZ proper Euler angles)
euler_angle_PSpincal –

Calculate Euler Angles (EA) between 2 coordinates systems (proper Eulerian angles)
rotate_3d –

Rotate a set of points (or set of vectors) from a OLD-coords to NEW-coords

Attributes:

old_orient –
new_orient –
DCM –
ROT –
EA –

`old_orient = np.asarray(old_orient)` `instance-attribute` ¶

`new_orient = np.asarray(new_orient)` `instance-attribute` ¶

`DCM = None` `instance-attribute` ¶

`ROT = None` `instance-attribute` ¶

`EA = None` `instance-attribute` ¶

`direction_cosine_matrix()` ¶

Calculate direction-cosines-matrix (DCM) between 2 coordinates systems

Returns:

Q –

3x3 array, the rotation matrix or matrix of direction cosines

Examples:

BT = thmd.CoordTransform(old_orient=oldAxis, new_orient=newAxis)
Q = BT.direction_cosine_matrix()

`rotation_matrix()` ¶

Calculate Rotation-matrix (R) as transpose of DCM By Cao Thang, Apr 2019, Update: May2020

`EA2ROT(euler_angle, unit='rad')` ¶

Calculate Rotation_Matrix Euler Angles (EA) between 2 coordinates systems (ZXZ proper Euler angles) This is just for testing, since we dont know whether input angles yield orthogonal axis or not

Parameters:

euler_angle (list) –

1x3 array/list (phi,theta,psi) in Rad or Deg
unit (str, default: 'rad' ) –

'rad' or 'deg' (default is rad)

Returns:

Q ( array ) –

3x3 array, the rotation matrix or matrix of direction cosines

Examples:

BT = thmd.CoordTransform()
DCM = BT.EulerAngle(euler_angle=[90,], unit='deg')

Notes

don't use arctan2()

`euler_angle(unit='rad')` ¶

Calculate Euler Angles (EA) between 2 coordinates systems (intrinsic ZXZ proper Euler angles) https://en.wikipedia.org/wiki/Euler_angles#Definition_by_intrinsic_rotations https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.transform.Rotation.as_euler.html#r72d546869407-1

Parameters:

unit='rad' –

'rad' or 'deg' (default is rad)

Returns:

Angle –

1x3 array (phi,theta,psi) in Rad (apply intrinsic ZXZ proper Euler)

Examples:

BT = thmd.CoordTransform(old_orient=oldAxis, new_orient=newAxis) phi,theta,psi = BT.EulerAngle(unit='deg')

Notes

don't use arctan2()
Rotation Matrix is as to tranpose of DCM, use Rotation Matrix to compute EA
To avoid devide by zero, we use 1e-64 instead of 0.

`euler_angle_PSpincal(euler_order='zxz', unit='rad', tol=1e-07)` ¶

Calculate Euler Angles (EA) between 2 coordinates systems (proper Eulerian angles)

Parameters:

unit –

'rad', 'deg' (default is rad)
euler_order='zxz' –

rotation order, lowercase ["zyx","zxy","yxz","xzy","xyz","yzx","zyz","zxz","yxy","yzy","xyx","xzx"]

Returns:

Angle –

1x3 array (phi,theta,psi) in Rad (apply extrinsic ZXZ proper Euler)

Notes

this module may define psi as phi, and vice versa. So becareful
should not use PSpincalc, since it produce unknown value?

"Refs"

[1] navpy not use 'ZXZ': https://navpy.readthedocs.io/en/latest/code_docs/coordinate_transformations.html [2] use this https://pypi.org/project/PSpincalc/ [3] https://github.com/tuxcell/PSpincalc/blob/master/PSpincalc/PSpincalc.py Ex: https://github.com/tuxcell/PSpincalc/blob/master/examples/examplesPSpincalc.ipynb

`rotate_3d(points)` ¶

Rotate a set of points (or set of vectors) from a OLD-coords to NEW-coords

Parameters:

points –

Nx3 array, contain coords in OLD coordinates systems

Returns:

points –

Nx3 array, contain coords in NEW coordinates systems

Examples:

BT = thmd.CoordTransform(old_orient=oldAxis, new_orient=newAxis)
newP = BT.rotate_3d(P)

`rot1axis(P, theta, axis='X')` ¶

Rotate array of points about 1 axis

Parameters:

P ( ) –

Nx3 array, contain input poits
theta ( ) –

the rotation angle in Degree
axis ( , default: 'X' ) –

Rotation axis

Returns:

outP –

Nx3 array, contain points after rotation

`check_right_hand(list_3vec=[[1, 0, 0], [0, 1, 0], [0, 0, 1]])` ¶

check right_hand_rule orthogonal of 3 vectors

`guess_right_hand(list_2vec=[[1, 0, 0], [0, 1, 0]])` ¶

give 2 vectors, then guess the third vector that satisfy right_hand_rule

`cartesian2spherical(xyz)` ¶

Convert cartesian coordinates to Spherical coordinates

Parameters:

xyz (array) –

Mx3 array contain Cartesian coordinates (X, Y, Z)

Returns:

shpCoord ( array ) –

Mx3 array contain Spherical coordinates (R, theta, phi). Also (radial distance, polar angle, azimuthal(longitude) angle)

Notes

The polar(theta) angle defined from from Z-axis down (zero at the North pole to 180° at the South pole)
adapted from this
There are many conventions that angles can be
- In geography, angles are in latitude/longitude or elevation/azimuthal form, polar angle is called latitude, measuze from XY-plane (ranges from -90° at the south pole to 90° at the north pole, with 0° at the Equator)
- In mathematics and physics, polar angle measured from Z-axis (zero at the North pole to 180° at the South pole)
- In scipy, polar angle is defined as Colatitude, which is a non-negative quantity, ranging from zero at the North pole to 180° at the South pole (same as commonly used in mathematics and physics)

`forcefield_info` ¶

This module contains some data for various ForceField. Data obtained from simulation

Classes:

EAM –

Create an Object (class) of Potential, contain some pre-setup information
ReaxFF –

Create an Object (class) of Potential, contain some pre-setup information

`EAM(atom_symbol, force_field, model_type='BULK', zDirect='001', thickness=20)` ¶

Create an Object (class) of Potential, contain some pre-setup information

    - Cu: Mishin2001,  Mendelev2008, Foiles1986

Al: Laird2000, Mishin1999, Mendelev2008, LiuEA2004, Sheng2011
V: Olsson2009

Parameters:

atom_symbol (str) –

define element, e.g., 'Al', 'Cu',...
force_field (str) –

the name of potential 'Cu' : 'Mishin-2001'; 'Foiles-1986';... 'Al' : 'Mishin-1999'; 'Sheng-2011';...
model_type (str, default: 'BULK' ) –

type of model (BULK or PLATE)
zDirect (str, default: '001' ) –

define the crystal orient along the z-direction simulation box, e,g., '001'/ '110'/ '111'
thickness (int, default: 20 ) –

define thickness in case model_type=PALTE

Returns:

obj –

Attributes:

atom_symbol (str) –

element
force_field (str) –

forcefield name.
cutoff (float) –

return cutoff of forcefield.
thermal_coeff (list) –

return values thermal expansion coefficients of input Structure

Stored DATA (these data are compute from several simulations, or from papers)

Methods:

lattice_constant –

Compute lattice constant at a specific temperature T.
atomic_volume_FCC –

Compute atomic volume at a specific temperature T.
melt_barrier –

Compute free energy barrier/atom of melting. For a system of N atoms $F = f*N^(2/3)$ , then barrier/atom $f = F/N^(2/3)$. This method return f(T)
estimate_FCCUBIC_liqsol –

Estimate value of FCCUBIC parameter in bulk solid/liquid at a specific temperature T.

Attributes:

coeff_FCCUBIC_liqsol_coexist –
coeff_FCCUBIC_tail_solid –
coeff_FCCUBIC_tail_liquid –
atom_symbol –
force_field –
thermal_coeff –
cutoff –
melt_barrier_coeff –

`coeff_FCCUBIC_liqsol_coexist = [-0.000467305786, 0.88174401]` `instance-attribute` ¶

`coeff_FCCUBIC_tail_solid = [-0.0017082985, 1.80537817]` `instance-attribute` ¶

`coeff_FCCUBIC_tail_liquid = [-9.33394967e-06, 0.608209715]` `instance-attribute` ¶

`atom_symbol = atom_symbol` `instance-attribute` ¶

`force_field = force_field` `instance-attribute` ¶

`thermal_coeff = D[thermal_key]` `instance-attribute` ¶

`cutoff = Rcut` `instance-attribute` ¶

`melt_barrier_coeff = melt_barrier[melt_barrier_key]` `instance-attribute` ¶

`lattice_constant(temp)` ¶

Compute lattice constant at a specific temperature T.

Parameters:

temp (float) –

input temperature.

Returns:

a ( float ) –

lattice constant as input temperature.

`atomic_volume_FCC(temp)` ¶

Compute atomic volume at a specific temperature T.

Parameters:

temp (float) –

input temperature.

Returns:

V ( float ) –

atomic volume (volume/atom) as input temperature.

`melt_barrier(temp)` ¶

Compute free energy barrier/atom of melting. For a system of N atoms $F = f*N^(2/3)$ , then barrier/atom $f = F/N^(2/3)$. This method return f(T)

Parameters:

temp (float) –

input temperature.

Returns:

barrier ( float ) –

free energy barrier/atom of melting as input temperature.

`estimate_FCCUBIC_liqsol(temp)` ¶

Estimate value of FCCUBIC parameter in bulk solid/liquid at a specific temperature T.

`ReaxFF(atom_symbol, force_field, model_type='BULK', zDirect='001')` ¶

Create an Object (class) of Potential, contain some pre-setup information

Methods:

lattice_constant –

Attributes:

atom_symbol –
force_field –
thermal_coeff –
cutoff –
melt_barrier_coeff –

`atom_symbol = atom_symbol` `instance-attribute` ¶

`force_field = force_field` `instance-attribute` ¶

`thermal_coeff = D[thermal_key]` `instance-attribute` ¶

`cutoff = Rcut` `instance-attribute` ¶

`melt_barrier_coeff = melt_barrier[melt_barrier_key]` `instance-attribute` ¶

`lattice_constant(Temp)` ¶

`polymer_mbuild` ¶

This module contains classes and functions to build models of atomic polymers See this Python package: [1] mBuild: https://mbuild.mosdef.org/en/stable/

See the files: D:\code\code_simulate\polymer_c21_pickup_hBN_PMMA ef_using_mBuild_foyer.ipynb

NOTEs

Due to mbuild cannot be installed with python 3.10, so import this package in functions to avoid checking in thmd

Functions:

PMMA_chain –

build polyPMMA from monomers
PVC_chain –
packing_lammps –

Packing polymer chains into box, and write LAMMPS file

`PMMA_chain(chain_len)` ¶

build polyPMMA from monomers

Parameters:

chain_len (int) –

number of monomers in each polymer = degree of polymerization

Returns:

chain ( compound ) –

polymer chain

`PVC_chain(chain_len)` ¶

`packing_lammps(chain, chain_num, density=None, box_size=None, forcefield_name=None, forcefield_files=None, atom_style='full', unit_style='metal', combining_rule='geometric', file_name='polymer.dat')` ¶

Packing polymer chains into box, and write LAMMPS file Packing based on either density or box_size.

Parameters:

chain (compound) –

polymer chain
chain_num (int) –

number of chains to be packed.
density (float, default: None ) –

density, unit in kg/m3 (= 1e-3 g/cm3)
box_size (list, default: None ) –

box_size = [3,3,3]
forcefield_name (str, default: None ) –

should be 'oplsaa'.
forcefield_files (str, default: None ) –

path to the *.xml file.
atom_style (str, default: 'full' ) –

atom_style of LAMMPS.
unit_style (str, default: 'metal' ) –

can be 'metal'/'real'/'lj'

`polymer_pysimm` ¶

This module contains classes to build models of atomic polymers See this Python package: [1] pysimm: A python package for simulation of molecular systems, 10.1016/j.softx.2016.12.002 source code: https://github.com/polysimtools/pysimm

`plot` ¶

This module provides functions and class to produce publication-quality figures. The main idea is to create a handy class which setup all necessary parameters for plotting, and then just call its methods to produce the publication-quality figures.

The core class is the Plotter, which is built on top of matplotlib.pyplot. The Plotter, therefore, has all useful functions from matplotlib.pyplot.

Assume we have a data frame df:

import numpy as np
import polars as pl
x = np.arange(0, 1000)
df = pl.DataFrame({'x': x, 'y': x**2})

Then, use the Plotter().plt to plot the data as a normal matplotlib.pyplot object. Adapt all functions from matplotlib.pyplot.

from thmd.plot import Plotter
plt = Plotter(style='thang').plt
plt.plot(df['x'], df['y'])
plt.show()

Refs

Customizing Matplotlib with style sheets and rcParams
Figure size Elsevier
- Single column: W = 90 mm (~3.5 in). H = W*⅘ = 2.8
- Double column: W = 190 mm (~7.5 in). H = 6

Classes:

Plotter –

Class to produce publication-quality figures with matplotlib.

Functions:

custom_styles –

style for matplotlib

`_COLOR = ['black', 'red', 'blue', 'green', 'magenta', 'orange', 'lime', 'cyan', 'violet', 'purple', 'olive', 'gray', 'yellow', 'navy', 'saddlebrown', 'darkgreen', 'lawngreen', 'lightgreen', 'steelblue', 'darkcyan', 'plum', 'slateblue', 'indigo']` `module-attribute` ¶

`_MARKER = ['o', 's', 'D', 'p', '*', '^', 'd', 'H', 'X', '>', '<', '^']` `module-attribute` ¶

`_LINE = {'solid': 'solid', 'dotted': 'dotted', 'dashed': 'dashed', 'dashdot': 'dashdot', 'loosely dotted': (0, (1, 10)), 'densely dotted': (0, (1, 1)), 'loosely dashed': (0, (5, 10)), 'densely dashed': (0, (5, 1)), 'loosely dashdotted': (0, (3, 10, 1, 10)), 'dashdotted': (0, (3, 5, 1, 5)), 'densely dashdotted': (0, (3, 1, 1, 1)), 'dashdotdotted': (0, (3, 5, 1, 5, 1, 5)), 'loosely dashdotdotted': (0, (3, 10, 1, 10, 1, 10)), 'densely dashdotdotted': (0, (3, 1, 1, 1, 1, 1))}` `module-attribute` ¶

`Plotter()` ¶

Class to produce publication-quality figures with matplotlib.

Parameters:

style (str) –

the plotting style.

Examples:

from thmd.plot import Plotter
plt = Plotter().plt
plt.plot(df['x'], df['y'])

Methods:

get_plt –

Shortcut to the matplotlib.pyplot object.
set_minor_locator –

Set minor locator for axes
get_PdfPages –

shortcut to matplotlib.backends.backend_pdf.PdfPages.
avail_styles –

List all available styles.

Attributes:

COLOR –
MARKER –
LINE –
info –

`COLOR = _COLOR` `instance-attribute` ¶

`MARKER = _MARKER` `instance-attribute` ¶

`LINE = _LINE` `instance-attribute` ¶

`info` `property` ¶

`get_plt(style: str = 'thang')` ¶

Shortcut to the matplotlib.pyplot object.

`set_minor_locator(ax: object, xminor=1, yminor=1)` ¶

Set minor locator for axes

`get_PdfPages()` ¶

shortcut to matplotlib.backends.backend_pdf.PdfPages.

`avail_styles()` ¶

List all available styles.

`_set_style(style: str = 'thang')` ¶

`custom_styles()` ¶

style for matplotlib

`post` ¶

Modules:

aseAtoms –

`aseAtoms` ¶

`qflow` ¶

Modules:

help_function –

Define support functions for qflow package.
machine_unist –
qsub_script –

Define the template for qsub job submission.
resource_unist –

`help_function` ¶

Define support functions for qflow package.

Functions:

read_nodes –

Read node list from file.
format_ignore_string –

Format qsub_string for ignore nodes.
format_select_string –

Format qsub_string for selected nodes.

`read_nodes(file_name: str) -> list` ¶

Read node list from file. Support delimiter: '|', ',', ';', ':', ' ', ' '

Args:
    file_name (str): File name.

Returns:
    list: List of unique nodes.

`format_ignore_string(nodes: list) -> str` ¶

Format qsub_string for ignore nodes.

`format_select_string(nodes: list) -> str` ¶

Format qsub_string for selected nodes.

`machine_unist` ¶

Functions:

Tachyon –

Define Machine for running jobs on Tachyon cluster at UNIST.
Eagle –

Define Machine for running jobs on Tachyon cluster at UNIST.
Lion –

Define Machine for running jobs on Tachyon cluster at UNIST.
Leopard –

Define Machine for running jobs on Tachyon cluster at UNIST.
CanCentosGpu –

Define Machine for running jobs on CAN cluster at UNIST.

Attributes:

sge_template –

`sge_template = {'context_type': 'SSHContext', 'batch_type': 'SGE', 'local_root': './', 'remote_root': '/home1/p001cao/work/w24_tmp', 'remote_profile': {'hostname': '10.0.7.53', 'username': 'p001cao', 'password': 'p001cao', 'port': 2022, 'timeout': 20}}` `module-attribute` ¶

`Tachyon(remote_path: str = '/home1/p001cao/work/w24_alff')` ¶

Define Machine for running jobs on Tachyon cluster at UNIST.

Parameters:

remote_path (str, default: '/home1/p001cao/work/w24_alff' ) –

Remote working folder.

`Eagle(remote_path: str = '/uwork/p001cao/work/w24_alff')` ¶

Define Machine for running jobs on Tachyon cluster at UNIST.

Parameters:

remote_path (str, default: '/uwork/p001cao/work/w24_alff' ) –

Remote working folder.

`Lion(remote_path: str = '/uwork/p001cao/work/w24_alff')` ¶

Define Machine for running jobs on Tachyon cluster at UNIST.

Parameters:

remote_path (str, default: '/uwork/p001cao/work/w24_alff' ) –

Remote working folder.

`Leopard(remote_path: str = '/uwork/p001cao/work/w24_alff')` ¶

Define Machine for running jobs on Tachyon cluster at UNIST.

Parameters:

remote_path (str, default: '/uwork/p001cao/work/w24_alff' ) –

Remote working folder.

`CanCentosGpu(remote_path: str = '/home/tha/work/w24_alff')` ¶

Define Machine for running jobs on CAN cluster at UNIST.

Parameters:

remote_path (str, default: '/home/tha/work/w24_alff' ) –

Remote working folder.

`qsub_script` ¶

Define the template for qsub job submission.

Functions:

qsub_sge_job –

Return the example script for submiting SGE job.
winBAT_submit_dispatcher –

Return the windows batch script for submitting DP dispatcher.
winBAT_cmd –

Return the example script for running command in cmd.
winBAT_wsl –

Return the example script for running command in cmd.
python_dispatcher –

Return the example script for using dpdispatcher.

`_DATA_PATH = Path(file).parent / 'lib_shell'` `module-attribute` ¶

`qsub_sge_job() -> str` ¶

Return the example script for submiting SGE job.

`winBAT_submit_dispatcher(conda_env: str = None) -> str` ¶

Return the windows batch script for submitting DP dispatcher.

`winBAT_cmd(command: str, conda_env: str = 'py11') -> str` ¶

Return the example script for running command in cmd.

`winBAT_wsl(command: str, conda_env: str = 'py11') -> str` ¶

Return the example script for running command in cmd.

`python_dispatcher() -> str` ¶

Return the example script for using dpdispatcher.

`resource_unist` ¶

Functions:

Tachyon_lammps_ib –

Define Resources for running LAMMPS on Tachyon cluster at UNIST, only InfiniBand nodes.
Tachyon_lammps_tcp –

Define Resources for running LAMMPS on Tachyon cluster at UNIST, only use TCP.
Tachyon_gpaw_ib –

Define Resources for running GPAW on Tachyon cluster at UNIST.
Tachyon_gpaw_ucx –

Define Resources for running GPAW on Tachyon cluster at UNIST.
Tachyon_gpaw_tcp –
Tachyon_lammps_ase –
Tachyon_ase_ib –
Tachyon_ase_tcp –
Lion_lammps –

Define Resources for running LAMMPS on Eagle cluster at UNIST.
Leopard_lammps –
Eagle_lammps –

Define Resources for running LAMMPS on Eagle cluster at UNIST.
CanCentosGpu_lammps –

Define Resources for running LAMMPS on Tachyon cluster at UNIST.
CanCentosGpu_gpaw –

Define Resources for running LAMMPS on Tachyon cluster at UNIST.
CanCentosGpu_ase –

Attributes:

tachyon_template –

`_DATA_PATH = Path(file).parent / '_data_cluster_UNIST'` `module-attribute` ¶

`tachyon_template = {'group_size': 1, 'queue_name': 'ib.q', 'kwargs': {'pe_name': 'mpi_8', 'job_name': 'zDP'}, 'cpu_per_node': 8, 'custom_flags': ['#$ -l h_rt=168:00:00', '#$ -l mem_used=600M'], 'source_list': ['/etc/profile.d/modules.sh'], 'module_list': ['mpi/openmpi4.1.7-clang17-IB', 'lammps/llvmOMPI4-dev'], 'envs': {'OMP_NUM_THREADS': 1, 'OMPI_MCA_btl_openib_allow_ib': 1, 'OMPI_MCA_btl': '^tcp'}}` `module-attribute` ¶

`_error_nodes_ucx()` ¶

`_error_nodes_openib()` ¶

`Tachyon_lammps_ib(n_cpu: int = 16)` ¶

Define Resources for running LAMMPS on Tachyon cluster at UNIST, only InfiniBand nodes.

Parameters:

n_cpu (int, default: 16 ) –

Number of CPU to run job.

`Tachyon_lammps_tcp(n_cpu: int = 8)` ¶

Define Resources for running LAMMPS on Tachyon cluster at UNIST, only use TCP.

Parameters:

n_cpu (int, default: 8 ) –

Number of CPU to run job.

`Tachyon_gpaw_ib(n_cpu: int = 16)` ¶

Define Resources for running GPAW on Tachyon cluster at UNIST.

Parameters:

n_cpu (int, default: 16 ) –

Number of CPU to run job.

`Tachyon_gpaw_ucx(n_cpu: int = 16)` ¶

Define Resources for running GPAW on Tachyon cluster at UNIST.

Parameters:

n_cpu (int, default: 16 ) –

Number of CPU to run job.

`Tachyon_gpaw_tcp(n_cpu: int = 8)` ¶

`Tachyon_lammps_ase(n_cpu: int = 16)` ¶

`Tachyon_ase_ib(n_cpu: int = 16)` ¶

`Tachyon_ase_tcp(n_cpu: int = 8)` ¶

`Lion_lammps(n_cpu: int = 12, queue_name: str = 'lion-normal.q', pe_name: str = 'lion-normal', time: str = '168:00:00')` ¶

Define Resources for running LAMMPS on Eagle cluster at UNIST.

Parameters:

n_cpu (int, default: 12 ) –

Number of CPU to run job. 12x for Lion, 8x for Leopard, 10x for Eagle.
queue_name (str, default: 'lion-normal.q' ) –

SGE's queue name. Available queues: - Lion: lion-normal.q, lion-short.q, lion-long.q - Leopard: leopard-normal.q, leopard-short.q, leopard-long.q
pe_name (str, default: 'lion-normal' ) –

SGE's Parallel environment name: mpi, mpi_1~mpi_20 - Lion: lion-normal, lion-short, lion-long - Leopard: leopard-normal, leopard-short, leopard-long
time (str, default: '168:00:00' ) –

Time limit for the job. normal: '168:00:00', long: '540:00:00'

`Leopard_lammps(n_cpu: int = 8, queue_name: str = 'leopard-normal.q', pe_name: str = 'leopard-normal', time: str = '168:00:00')` ¶

`Eagle_lammps(n_cpu: int = 10, pe_name: str = 'mpi_10')` ¶

Define Resources for running LAMMPS on Eagle cluster at UNIST.

Parameters:

n_cpu (int, default: 10 ) –

Number of CPU to run job. 12x for Lion, 8x for Leopard, 10x for Eagle.
pe_name (str, default: 'mpi_10' ) –

SGE's Parallel environment name: mpi, mpi_1~mpi_20 - Eagle: mpi, mpi_1~mpi_20

`CanCentosGpu_lammps(n_cpu: int = 12, pe_name: str = 'mpi_12')` ¶

Define Resources for running LAMMPS on Tachyon cluster at UNIST.

Parameters:

pe_name (str, default: 'mpi_12' ) –

SGE's Parallel environment name.
n_cpu (int, default: 12 ) –

Number of CPU to run job.

`CanCentosGpu_gpaw(n_cpu: int = 12, pe_name: str = 'mpi_12')` ¶

Define Resources for running LAMMPS on Tachyon cluster at UNIST.

Parameters:

pe_name (str, default: 'mpi_12' ) –

SGE's Parallel environment name.
n_cpu (int, default: 12 ) –

Number of CPU to run job.

`CanCentosGpu_ase(n_cpu: int = 12, pe_name: str = 'mpi_12')` ¶

`recipe` ¶

Modules:

git –
gpaw –
lammps –
mace –

Some notes
mlff –
ovito –
plumed –
sevennet –

Some notes

`git` ¶

`gpaw` ¶

Modules:

arg_parser –

Some notes
cli_gpaw_check_PWcutoff –

Some notes
cli_gpaw_check_kpoints –

Some notes
cli_gpaw_check_kpoints_density –

Some notes
cli_gpaw_optimize –

Some notes
cli_gpaw_singlepoint –

Some notes

`arg_parser` ¶

Some notes

RFE

Python CLI variables: https://stackoverflow.com/questions/4033723/how-do-i-access-command-line-arguments
argparse: https://stackoverflow.com/questions/20063/whats-the-best-way-to-parse-command-line-arguments
Parse a list in argparse: https://stackoverflow.com/questions/15753701/how-can-i-pass-a-list-as-a-command-line-argument-with-argparse

Functions:

args_optimize –

`args_optimize()` ¶

`cli_gpaw_check_PWcutoff` ¶

Some notes - calc.new() will create a new calculator that inherits all parameters from the current calculator, except for txt and timer

RFE

convergence-checks: https://wiki.fysik.dtu.dk/gpaw/tutorialsexercises/structureoptimization/water/water.html#convergence-checks
Gpaw tools: https://github.com/lrgresearch/gpaw-tools/blob/main/optimizations/optimize_ecut.py

Attributes:

output_dir –
parallel_args –
parser –
args –
extxyz_file –
ecut –
atoms –
calc –
new_calc –
pe –
out_file –

`output_dir = Path('output')` `module-attribute` ¶

`parallel_args = {'sl_auto': True, 'use_elpa': True}` `module-attribute` ¶

`parser = argparse.ArgumentParser(description='Check convergence of PW energy cutoff')` `module-attribute` ¶

`args = parser.parse_args()` `module-attribute` ¶

`extxyz_file = args.poscar` `module-attribute` ¶

`ecut = args.ecutoff` `module-attribute` ¶

`atoms = read(extxyz_file, format='extxyz_file', index='-1')` `module-attribute` ¶

`calc = GPAW(mode=PW(550), xc='PBE', occupations=FermiDirac(0.01), kpts={'density': 15, 'gamma': False}, parallel=parallel_args)` `module-attribute` ¶

`new_calc = calc.new(mode=PW(ecut), txt=f'{output_dir}/calc_ecut_{ecut}.txt')` `module-attribute` ¶

`pe = atoms.get_potential_energy() / len(atoms)` `module-attribute` ¶

`out_file = f'{output_dir}/check_ecut_{ecut}.txt'` `module-attribute` ¶

`cli_gpaw_check_kpoints` ¶

Some notes

RFE

convergence-checks: https://wiki.fysik.dtu.dk/gpaw/tutorialsexercises/structureoptimization/water/water.html#convergence-checks
Gpaw tools: https://github.com/lrgresearch/gpaw-tools/blob/main/optimizations/optimize_kpoints.py
Parse a list in argparse: https://stackoverflow.com/questions/15753701/how-can-i-pass-a-list-as-a-command-line-argument-with-argparse

Attributes:

output_dir –
parallel_args –
parser –
args –
extxyz_file –
pbc –
atoms –
calc –
pe –
out_file –

`output_dir = Path('output')` `module-attribute` ¶

`parallel_args = {'sl_auto': True, 'use_elpa': True}` `module-attribute` ¶

`parser = argparse.ArgumentParser(description='Check convergence of PW energy cutoff')` `module-attribute` ¶

`args = parser.parse_args()` `module-attribute` ¶

`extxyz_file = args.poscar` `module-attribute` ¶

`pbc = [int(item) for item in args.pbc.split(' ')]` `module-attribute` ¶

`atoms = read(extxyz_file, format='extxyz_file', index='-1')` `module-attribute` ¶

`calc = GPAW(mode=PW(550), xc='PBE', occupations=FermiDirac(0.01), kpts={'size': (nkx, nky, nkz), 'gamma': True}, txt=f'{output_dir}/calc_kpoints_{nkx}x{nky}x{nkz}.txt', parallel=parallel_args)` `module-attribute` ¶

`pe = atoms.get_potential_energy() / len(atoms)` `module-attribute` ¶

`out_file = f'{output_dir}/check_kpoints_{nkx}x{nky}x{nkz}.txt'` `module-attribute` ¶

`cli_gpaw_check_kpoints_density` ¶

Some notes

RFE

convergence-checks: https://wiki.fysik.dtu.dk/gpaw/tutorialsexercises/structureoptimization/water/water.html#convergence-checks
Gpaw tools: https://github.com/lrgresearch/gpaw-tools/blob/main/optimizations/optimize_kpoints.py

Attributes:

output_dir –
parallel_args –
parser –
args –
extxyz_file –
kdensity –
pbc –
atoms –
calc –
pe –
out_file –

`output_dir = Path('output')` `module-attribute` ¶

`parallel_args = {'sl_auto': True, 'use_elpa': True}` `module-attribute` ¶

`parser = argparse.ArgumentParser(description='Check convergence of PW energy cutoff')` `module-attribute` ¶

`args = parser.parse_args()` `module-attribute` ¶

`extxyz_file = args.poscar` `module-attribute` ¶

`kdensity = args.kdensity` `module-attribute` ¶

`pbc = [int(item) for item in args.pbc.split(' ')]` `module-attribute` ¶

`atoms = read(extxyz_file, format='extxyz_file', index='-1')` `module-attribute` ¶

`calc = GPAW(mode=PW(550), xc='PBE', occupations=FermiDirac(0.01), kpts={'density': kdensity, 'gamma': True}, txt=f'{output_dir}/calc_kdensity_{kdensity}.txt', parallel=parallel_args)` `module-attribute` ¶

`pe = atoms.get_potential_energy() / len(atoms)` `module-attribute` ¶

`out_file = f'{output_dir}/check_kdensity_{kdensity}.txt'` `module-attribute` ¶

`cli_gpaw_optimize` ¶

Some notes - Must set txt='calc.txt' in GPAW calculator for backward files. - param_yaml must contain - a dict gpaw with GPAW parameters. - a dict optimize with ASE optimization parameters.

Attributes:

parser –
args –
configfile –
pdict –
extxyz_file –
atoms –
gpaw_arg –
params –
calc –
opt_args –
relax_dim –
pbc –
fmax –
nsteps –
atoms_filter –
opt –
pot_energy –
forces –
stress –
atoms_fake –
output_file –

`parser = argparse.ArgumentParser(description='Optimize structure using GPAW')` `module-attribute` ¶

`args = parser.parse_args()` `module-attribute` ¶

`configfile = args.param` `module-attribute` ¶

`pdict = load_setting_file(configfile)` `module-attribute` ¶

`extxyz_file = pdict['input_extxyz_path']` `module-attribute` ¶

`atoms = read(extxyz_file, format='extxyz', index='-1')` `module-attribute` ¶

`gpaw_arg = pdict.get('gpaw_arg', {})` `module-attribute` ¶

`params = {'mode': {'name': 'pw', 'ecut': 500}, 'xc': 'PBE', 'convergence': {'energy': 1e-06, 'density': 0.0001, 'eigenstates': 1e-08}, 'occupations': {'name': 'fermi-dirac', 'width': 0.01}, 'txt': 'calc_optimize.txt'}` `module-attribute` ¶

`calc = GPAW(**params)` `module-attribute` ¶

`opt_args = pdict.get('optimize', {})` `module-attribute` ¶

`relax_dim = opt_args.get('relax_dim', None)` `module-attribute` ¶

`pbc = atoms.get_pbc()` `module-attribute` ¶

`fmax = opt_args.get('fmax', 0.05)` `module-attribute` ¶

`nsteps = opt_args.get('nsteps', 10000)` `module-attribute` ¶

`atoms_filter = FrechetCellFilter(atoms, mask=relax_dim)` `module-attribute` ¶

`opt = BFGS(atoms_filter)` `module-attribute` ¶

`pot_energy = atoms.get_potential_energy()` `module-attribute` ¶

`forces = atoms.get_forces()` `module-attribute` ¶

`stress = atoms.get_stress(voigt=True)` `module-attribute` ¶

`atoms_fake = atoms.copy()` `module-attribute` ¶

`output_file = extxyz_file.replace('.extxyz', '_labeled.extxyz')` `module-attribute` ¶

`cli_gpaw_singlepoint` ¶

Some notes - Must set txt='calc.txt' in GPAW calculator for backward files. - param_yaml must contain - a dict gpaw with GPAW parameters.

Attributes:

parser –
args –
configfile –
pdict –
extxyz_file –
atoms –
gpaw_arg –
params –
calc –
pot_energy –
forces –
stress –
atoms_fake –
output_file –

`parser = argparse.ArgumentParser(description='Optimize structure using GPAW')` `module-attribute` ¶

`args = parser.parse_args()` `module-attribute` ¶

`configfile = args.param` `module-attribute` ¶

`pdict = load_setting_file(configfile)` `module-attribute` ¶

`extxyz_file = pdict['input_extxyz_path']` `module-attribute` ¶

`atoms = read(extxyz_file, format='extxyz', index='-1')` `module-attribute` ¶

`gpaw_arg = pdict.get('gpaw_arg', {})` `module-attribute` ¶

`params = {'mode': {'name': 'pw', 'ecut': 500}, 'xc': 'PBE', 'convergence': {'energy': 1e-06, 'density': 0.0001, 'eigenstates': 1e-08}, 'occupations': {'name': 'fermi-dirac', 'width': 0.01}, 'txt': 'calc_singlepoint.txt'}` `module-attribute` ¶

`calc = GPAW(**params)` `module-attribute` ¶

`pot_energy = atoms.get_potential_energy()` `module-attribute` ¶

`forces = atoms.get_forces()` `module-attribute` ¶

`stress = atoms.get_stress(voigt=True)` `module-attribute` ¶

`atoms_fake = atoms.copy()` `module-attribute` ¶

`output_file = extxyz_file.replace('.extxyz', '_labeled.extxyz')` `module-attribute` ¶

`lammps` ¶

`mace` ¶

Some notes

RFE

SevenNet repo: https://github.com/MDIL-SNU/SevenNet

Functions:

cli_mace_optimize –

Return filepath of the script for running optimization using SevenNet.

`_DATA_PATH = Path(file).parent / 'lib_mace'` `module-attribute` ¶

`cli_mace_optimize(copy_to: str = None) -> str` ¶

Return filepath of the script for running optimization using SevenNet.

`mlff` ¶

`ovito` ¶

Functions:

scale_RGB –

Function to convert RGB color code from scale 0-255 to scale 0-1.
mod_set_prop_atom_name –

Modifier to set atom names
mod_set_prop_atom_color_PMMAori –

Modifier to assign atom colors based on atom_names.
delete_pipelines –

delete all existed viewports, pipelines, and scene

`scale_RGB(RGB=(255, 255, 255))` ¶

Function to convert RGB color code from scale 0-255 to scale 0-1.

Parameters:

RGB (tuple, default: (255, 255, 255) ) –

RGB code in scale 0-255

Returns:

rgb ( tuple ) –

RGB code in scale 0-1

Examples:

rgb = scale_RGB((255,255,255)))

Quote

rgb-values-to-0-to-1-scale

`mod_set_prop_atom_name(frame, data)` ¶

Modifier to set atom names

Examples:

from thmd.visual.ovito_modifier import mod_set_prop_atom_name
from ovito.io import import_file

pipeline = import_file("test.cfg")
pipeline.add_to_scene()
## add mod
dict_name = {'type_id':[1, 2], 'atom_name':['C', 'H']}
pipeline.modifiers.append(mod_set_prop_atom_name)

Note

So far, can not a custom argument to modifier, see here. So we need to define a global variable before using this function
```
dict_name = {'type_id':(1, 2), 'atom_name':('C', 'H')}
```
Do not use 'return` in modifier
the underscore notation mean modifiable version of the quantity in ovito

Quote

Pass custom args to modifier
ovito.data.Property - type.id, type.name, type.color, type.radius

`mod_set_prop_atom_color_PMMAori(frame, data)` ¶

Modifier to assign atom colors based on atom_names.

Examples:

from thmd.visual.ovito_modifier import mod_set_prop_atom_color_PMMAori
from ovito.io import import_file

pipeline = import_file("test.cfg")
pipeline.add_to_scene()
## add mod
pipeline.modifiers.append(mod_set_prop_atom_color_PMMAori)

`delete_pipelines(viewports: list[object] = [], pipelines: list[object] = [], scene: object = None)` ¶

delete all existed viewports, pipelines, and scene

Parameters:

viewports (obj, default: [] ) –

list of ovito viewport objects
pipelines (obj, default: [] ) –

list of ovito pipeline objects
scene (obj, default: None ) –

ovito scene object

`plumed` ¶

Functions:

script_FCCUBIC –

PLUMED script to compute FCCUBIC
script_LOCAL_CRYSTALINITY –

PLUMED script to compute LOCAL_CRYSTALINITY
script_LOGMFD –

PLUMED script to compute LOGFMD

`script_FCCUBIC(a_fcc, zDirect, label='mcv', alpha=27, partialCompute=False, atoms='@mdatoms', atomsA=None, atomsB=None, options='')` ¶

PLUMED script to compute FCCUBIC

Parameters:

a_fcc (float) –

Lattice constant of FCC crystal
zDirect (str) –

specify the z-direction of crystal
label (str, default: 'mcv' ) –

label of PLUMED command
alpha (int, default: 27 ) –

ALPHA parameter to compute FCCUBIC colvar.
atoms (str, default: '@mdatoms' ) –

specify atom-ids in computed group.
partialCompute (bool, default: False ) –

compute for some atoms.
atomsA (str, default: None ) –

specify atom-ids in group A.
atomsB (str, default: None ) –

specify atom-ids in group B.
options (str, default: '' ) –

add options.

Returns:

list ( list ) –

list of strings.

`script_LOCAL_CRYSTALINITY(a_fcc, zDirect, label='mcv', vectors=[[1, 0, 0], [0, 1, 0], [0, 0, 1]], atoms='@mdatoms', options='')` ¶

PLUMED script to compute LOCAL_CRYSTALINITY

Parameters:

a_fcc (float) –

Lattice constant of FCC crystal
zDirect (str) –

specify the z-direction of crystal
label (str, default: 'mcv' ) –

label of PLUMED command
vectors (list, default: [[1, 0, 0], [0, 1, 0], [0, 0, 1]] ) –

2xN list of lists, to specify directions of reciprocal vectors.
atoms (str, default: '@mdatoms' ) –

specify atom-ids in computed group.
options (str, default: '' ) –

add options.

Returns:

list ( list ) –

list of strings.

`script_LOGMFD(ARG, FICT, FICT_MIN, FICT_MAX, TEMP, DELTA_T, INTERVAL, KAPPA, deltaF, deltaX, kB, label='mfd', FLOG=5000, MFDstat='VS')` ¶

PLUMED script to compute LOGFMD

Parameters:

ARG (str) –

the scalar input for this action
FICT (float) –

The initial values of the fictitious dynamical variables
FICT_MIN (float) –

Boundaries of CV_space
FICT_MAX (float) –

Boundaries of CV_space
TEMP (float) –

Temperature of the fictitious dynamical variables
DELTA_T (float) –

Time step for the fictitious dynamical variables (MFD step)
INTERVAL (int) –

Period of MD steps ( Nm) to update fictitious dynamical variables
KAPPA (int) –

Spring constant of the harmonic restraining potential for the fictitious dynamical variables
deltaF (float) –

Energy Barrier to estimate ALPHA (Alpha parameter for LogMFD)
deltaX (float) –

CV distance at each MFDstep, to estimate MFICT, VFICT (mass & velocity of fictitious dynamical variable)
kB (float) –

Boltzmann constant
label (str, default: 'mfd' ) –

label of PLUMED command
FLOG (float, default: 5000 ) –

The initial free energy value in the LogMFD, initial F(X)
MFDstat (str, default: 'VS' ) –

Type of thermostat for the fictitious dynamical variables. NVE, NVT, VS are available.

Returns:

list ( list ) –

list of strings.

`sevennet` ¶

Some notes

RFE

SevenNet repo: https://github.com/MDIL-SNU/SevenNet

Functions:

cli_7net_optimize –

Return filepath of the script for running optimization using SevenNet.

`_DATA_PATH = Path(file).parent / 'lib_7net'` `module-attribute` ¶

`cli_7net_optimize(copy_to: str = None) -> str` ¶

Return filepath of the script for running optimization using SevenNet.

`util` ¶

Modules:

check_installation –
compute_angle –
compute_distance –
compute_tensor –
detect_sign_change –
fit –
grid_box –
many_stuff –
row_operation –
string_index –
unit –

This module to convert unit of some physical properties

`check_installation` ¶

Functions:

check_install_ovito –
check_install_googlesearch –

`check_install_ovito()` ¶

`check_install_googlesearch()` ¶

`compute_angle` ¶

Functions:

angle_vector2vectors –

copmute angles between a vector with set of vectors

`angle_vector2vectors(fixVector, arrayVectors, unit='rad')` ¶

copmute angles between a vector with set of vectors

`compute_distance` ¶

Functions:

dist2_point2points –

Compute bond_len and postion_vetors from 1 point to a list of points
dist2_points2line –

Compute bond_len and postion_vetors from 1 point to a list of points
closest_points2line –

Find all points locate inside a checkin-distance "dist" from a line.
closest_points2multilines –

Find all points locate inside a checkin-distance "dist" from multilines.

`dist2_point2points(point, points)` ¶

Compute bond_len and postion_vetors from 1 point to a list of points

Parameters:

point (list array) –

coordinate of 1 point.
points (list array) –

2d-list of coordinates of points/point.

Returns:

df ( DataFrame ) –

pd.DataFrame constains distance and component of connecting vectors.

`dist2_points2line(points, line=[(0, 0, 0), (0, 0, 0)])` ¶

Compute bond_len and postion_vetors from 1 point to a list of points Ref: https://stackoverflow.com/questions/39840030/distance-between-point-and-a-line-from-two-points

Parameters:

points (list array DataFrame) –

list of coordinates of points/point.
line (list array, default: [(0, 0, 0), (0, 0, 0)] ) –

2d-array contains coordinates to define a line.

Returns:

d ( float list ) –

distances between points and a line.

`closest_points2line(points, line=[(0, 0, 0), (0, 0, 0)], distance=0, Xbound=None, Ybound=None, Zbound=None)` ¶

Find all points locate inside a checkin-distance "dist" from a line.

Parameters:

points (list array) –

list of coordinates of points/point.
line (list array, default: [(0, 0, 0), (0, 0, 0)] ) –

[[x1,y1,z1], [x1,y2,z2]]: 2d-list contains coordinates to define a line.
distance (float, default: 0 ) –

the checkin-distance.
Xbound (tuple, default: None ) –

define the boundaries for checking. Xbound='line': use the lengths of lines as bounds. Xbound=None: extend to INF. Xbound = (xlo, xhi)
Ybound (tuple, default: None ) –

define the boundaries for checking.
Zbound (tuple, default: None ) –

define the boundaries for checking.

Returns:

ds_idx ( Series ) –

Series of indices of points within the checkin-distance

`closest_points2multilines(points, multilines=[], distance=0, Xbound=None, Ybound=None, Zbound=None)` ¶

Find all points locate inside a checkin-distance "dist" from multilines. The Bound is set as the line-lengths.

Parameters:

points (list array) –

list of coordinates of points/point.
multilines (list, default: [] ) –

list of pair-points, each pair-point contains coordinates of 2 points to define a line used in 'closest_points2line'.
distance (float, default: 0 ) –

the checkin-distance.

Returns:

ds_idx ( Series ) –

Series of indices of points within the checkin-distance

`compute_tensor` ¶

Functions:

ke_tensor –

Compute Kinetic Energy tensors, and Temp
stress_tensor –

Compute local pressure/stress

`ke_tensor(vel, mass, kb)` ¶

Compute Kinetic Energy tensors, and Temp Args: vel (array): Nx3 array of per-atom velocity mass (array): Nx3 array of atomic mass inUNIT (str): ['angstrom','ps','amu','eV'], outUNIT=['eV','K'] Returns: Kinetic energy tensor, Kinetic scalar, Temperature scalar

`stress_tensor(per_atom_stress_tensor, atomic_volume, unitFac=1)` ¶

Compute local pressure/stress Args: per_atom_stress_tensor : Nx6 array of the per-atom stress tensor atomVol : Nx1 vector of atomVol inUNIT=['bar','angstrom'], outUNIT=['bar'] → unitFac=1e-4 for ['GPa'] Returns: pressure scalar Stress tensor

`detect_sign_change` ¶

Functions:

detect_sign_change –

determine points where line y=y(x) change its sign

`detect_sign_change(y, x=[])` ¶

determine points where line y=y(x) change its sign

Parameters:

y –

Nx1 arrays, contains dependent variable y
x –

(Optinal) Nx1 arrays, contains independent variable x of line y(x)

Returns:

idx –

1d array of indices where sign changes

`fit` ¶

Modules:

curve_intersect –
fit_root –
user_lmfit –

Functions:

curve_intersect_shapely –

find the intersection points between 2 lines
polyfit –

Polynomial fitting for y = f(x) relation, which returns uncertainties of coefficients.
find_slope –

Compute slope of a linear relation y = A + B*x.
find_roots –

find roots for y = f(x) relation in polynomial form.
find_extrema –

find extrema points (der=0) for y = f(x) relation in polynomial form.
find_convergence –

y is a function of x, then find value of x that y converge with a tolerance < tol.
extrapolate –

`curve_intersect_shapely(line1, line2)` ¶

find the intersection points between 2 lines Args: line1 (array like): Nx2 arrays, contains data points of curve 1 line2 (array like): Nx2 arrays, contains data points of curve 2

Returns:

Points ( array like ) –

Nx2 arrays, contains data points of intersection points

`polyfit(x, y, deg=1, sigma_y=None, uncert=False, **kwargs)` ¶

Polynomial fitting for y = f(x) relation, which returns uncertainties of coefficients. The fitted polynomial(s) are in the form

p(x) = a0 + a1*x + a2*x^2 + ... + an*x^n

Parameters:

x,y (array) –

1D array of x and y data
deg (int, default: 1 ) –

degree of polynomial
sigma_y (array, default: None ) –

1D array of standard deviation of y data. Use in weighted least square fitting.
uncert (bool, default: False ) –

return uncertainties of coefficients. Defaults to False.
**kwargs –

additional arguments, adapt all args from np.polyfit

Returns:

pars ( array ) –

1D array of coefficients
uncert ( array ) –

1D array of standard deviation of coefficients

`find_slope(x, y)` ¶

Compute slope of a linear relation y = A + B*x.

Parameters:

x (list) –

a list/array of x value
y (list) –

a list/array of y value

Return

slope (float): the slope to linear relation

`find_roots(x, y, order=1)` ¶

find roots for y = f(x) relation in polynomial form.

Parameters:

x (list) –

a list/array of x value
y (list) –

a list/array of y value
order (int, default: 1 ) –

order of polynomial. Defaults to 1.

Returns:

roots ( list ) –

list of roots

`find_extrema(x, y, order=2, retun_fit=False)` ¶

find extrema points (der=0) for y = f(x) relation in polynomial form.

Parameters:

x (list) –

a list/array of x value
y (list) –

a list/array of y value
order (int, default: 2 ) –

order of polynomial. Defaults to 2.
retun_fit (bool, default: False ) –

return fitted params. Defaults to False.

Returns:

p ( list ) –

list of flex points.

`find_convergence(x, y, tol=1e-06, grid_size=None, convergence_side='right')` ¶

y is a function of x, then find value of x that y converge with a tolerance < tol.

Parameters:

x (list) –

a list/array of x values
y (list) –

a list/array of y values
tol (float, default: 1e-06 ) –

tolerance of the slope (dy/dx). Defaults to 1e-6.
grid_size (float, default: None ) –

set grid size of x if need to interpolate finer data. Defaults to None.
convergence_side (str, default: 'right' ) –

find the first x on the left or right side. Defaults to "right".

Returns:

–

best_x, best_y (float): found x and y value that y converge with a tolerance < tol.

Notes

Value df['dy'] = df['y'].diff() depends on grid size. Should use df['dy'] = df['y'].diff()/ df['x'].diff().

References

https://docs.scipy.org/doc/scipy/tutorial/interpolate.html

`extrapolate(x, y, left_side: list = None, right_side: list = None, grid_size=0.1) -> tuple[list, list]` ¶

`curve_intersect` ¶

Functions:

curve_intersect_shapely –

find the intersection points between 2 lines

`curve_intersect_shapely(line1, line2)` ¶

find the intersection points between 2 lines Args: line1 (array like): Nx2 arrays, contains data points of curve 1 line2 (array like): Nx2 arrays, contains data points of curve 2

Returns:

Points ( array like ) –

Nx2 arrays, contains data points of intersection points

`___curve_intersect_numpy(curve1, curve2, degree=3, bounds=None)` ¶

find the intersection points between 2 curves

Parameters:

curve1 (array like) –

Nx2 arrays in form (x,y), contains data points of curve 1
curve2 (array like) –

Nx2 arrays in form (x,y), contains data points of curve 2
degree (int, default: 3 ) –

degree of polynomial function to fit the curve
bounds (tuple, default: None ) –

(min, max) of x, to select points in bounds

Returns:

Points ( array like ) –

Nx2 arrays, contains data points of intersection points

Refs

[1] http://t.ly/2_a-K

`fit_root` ¶

Functions:

polyfit –

Polynomial fitting for y = f(x) relation, which returns uncertainties of coefficients.
find_slope –

Compute slope of a linear relation y = A + B*x.
find_roots –

find roots for y = f(x) relation in polynomial form.
find_extrema –

find extrema points (der=0) for y = f(x) relation in polynomial form.
find_convergence –

y is a function of x, then find value of x that y converge with a tolerance < tol.
extrapolate –

`polyfit(x, y, deg=1, sigma_y=None, uncert=False, **kwargs)` ¶

Polynomial fitting for y = f(x) relation, which returns uncertainties of coefficients. The fitted polynomial(s) are in the form

p(x) = a0 + a1*x + a2*x^2 + ... + an*x^n

Parameters:

x,y (array) –

1D array of x and y data
deg (int, default: 1 ) –

degree of polynomial
sigma_y (array, default: None ) –

1D array of standard deviation of y data. Use in weighted least square fitting.
uncert (bool, default: False ) –

return uncertainties of coefficients. Defaults to False.
**kwargs –

additional arguments, adapt all args from np.polyfit

Returns:

pars ( array ) –

1D array of coefficients
uncert ( array ) –

1D array of standard deviation of coefficients

`find_slope(x, y)` ¶

Compute slope of a linear relation y = A + B*x.

Parameters:

x (list) –

a list/array of x value
y (list) –

a list/array of y value

Return

slope (float): the slope to linear relation

`find_roots(x, y, order=1)` ¶

find roots for y = f(x) relation in polynomial form.

Parameters:

x (list) –

a list/array of x value
y (list) –

a list/array of y value
order (int, default: 1 ) –

order of polynomial. Defaults to 1.

Returns:

roots ( list ) –

list of roots

`find_extrema(x, y, order=2, retun_fit=False)` ¶

find extrema points (der=0) for y = f(x) relation in polynomial form.

Parameters:

x (list) –

a list/array of x value
y (list) –

a list/array of y value
order (int, default: 2 ) –

order of polynomial. Defaults to 2.
retun_fit (bool, default: False ) –

return fitted params. Defaults to False.

Returns:

p ( list ) –

list of flex points.

`find_convergence(x, y, tol=1e-06, grid_size=None, convergence_side='right')` ¶

y is a function of x, then find value of x that y converge with a tolerance < tol.

Parameters:

x (list) –

a list/array of x values
y (list) –

a list/array of y values
tol (float, default: 1e-06 ) –

tolerance of the slope (dy/dx). Defaults to 1e-6.
grid_size (float, default: None ) –

set grid size of x if need to interpolate finer data. Defaults to None.
convergence_side (str, default: 'right' ) –

find the first x on the left or right side. Defaults to "right".

Returns:

–

best_x, best_y (float): found x and y value that y converge with a tolerance < tol.

Notes

Value df['dy'] = df['y'].diff() depends on grid size. Should use df['dy'] = df['y'].diff()/ df['x'].diff().

References

https://docs.scipy.org/doc/scipy/tutorial/interpolate.html

`extrapolate(x, y, left_side: list = None, right_side: list = None, grid_size=0.1) -> tuple[list, list]` ¶

`___uncertainty_weighted_fit_linear(x, y, sigma_y=None)` ¶

Compute uncertainties in coefficients A and B of y = A + B*x from the weighted least square fitting.

Apply when the measured data y_i has different uncertainties sigma_y_i, the weights w_i = 1/(sigma_y_i^2).

Parameters:

x,y (array) –

1D array of x and y data
sigma_y (array, default: None ) –

1D array of standard deviation of y data. Use in weighted least square fitting.

Returns:

pars ( array ) –

1D array of coefficients A and B
uncertainties ( tuple ) –

1D array (sigma_intercept, sigma_slope) of standard deviation of A and B

References

Taylor_1997_An introduction to error analysis: the study of uncertainties in physical measurements, page 198.
https://numpy.org/doc/stable/reference/generated/numpy.polynomial.polynomial.polyfit.html

`user_lmfit` ¶

Classes:

UserLmfit –

The class contains set of objective function of fitting use by LMFIT package

`UserLmfit` ¶

The class contains set of objective function of fitting use by LMFIT package NOTEs: - defined function followed the convection of LMFIT: the first argument of the function is taken as the independent variable, held in independent_vars, and the rest of the functions positional arguments (and, in certain cases, keyword arguments – see below) are used for Parameter names. https://lmfit.github.io/lmfit-py/model.html - This Class defines curve-forms that are not vailable in LMFIT's built-in models

Attributes: swType : (default='RATIONAL') Type of witching function, r0, d0 : The r_0 parameter of the switching function
Methods: fFunc : compute & return value and derivation of sw function fDmax : estimate value of Dmax

Ex: func = thmd.CurveLib.Linear(x)

Methods:

Linear –

this func is available in LMFIT, just play as an example here
inverseTemperature –
ExpDecay –
sizeEffect –

system size-dependence on term N^(⅔
unNormalGaussian –

The unNormalize Gaussian function
NormalGaussian –

The Normalize Gaussian function
sum_2unNormalGaussian –

The sum of 2 Gaussian function
sum_3unNormalGaussian –

The sum of 3 Gaussian function
sum_4unNormalGaussian –

The sum of 4 Gaussian function
sum_5unNormalGaussian –

The sum of 5 Gaussian function
sum_2NormalGaussian –

The sum of 2 Gaussian function
sum_3NormalGaussian –

The sum of 3 Gaussian function
sum_4NormalGaussian –

The sum of 4 Gaussian function
sum_5NormalGaussian –

The sum of 5 Gaussian function
DoseResp –

Dose-response curve with variable Hill slope given by parameter 'p'.
BiDoseResp –

Biphasic Dose Response Function,
Carreau –

Carreau-Yasuda model to describe pseudoplastic flow with asymptotic viscosities at zero and infinite shear rates
Cross –

Cross model to describe pseudoplastic flow with asymptotic viscosities at zero and infinite shear rates
GammaCFD –

Gamma cumulative distribution function

Linear(x, a0, a1) ¶

this func is available in LMFIT, just play as an example here

inverseTemperature(x, a, b) ¶

ExpDecay(x, A, lambd) ¶

sizeEffect(x, a, b) ¶

system size-dependence on term N^(⅔

unNormalGaussian(x, amp, cen, sig) ¶

The unNormalize Gaussian function

NormalGaussian(x, amp, cen, sig) ¶

The Normalize Gaussian function

sum_2unNormalGaussian(x, amp1, amp2, cen1, cen2, sig1, sig2) ¶

The sum of 2 Gaussian function

sum_3unNormalGaussian(x, amp1, amp2, amp3, cen1, cen2, cen3, sig1, sig2, sig3) ¶

The sum of 3 Gaussian function

sum_4unNormalGaussian(x, amp1, amp2, amp3, amp4, cen1, cen2, cen3, cen4, sig1, sig2, sig3, sig4) ¶

The sum of 4 Gaussian function

sum_5unNormalGaussian(x, amp1, amp2, amp3, amp4, amp5, cen1, cen2, cen3, cen4, cen5, sig1, sig2, sig3, sig4, sig5) ¶

The sum of 5 Gaussian function

sum_2NormalGaussian(x, amp1, amp2, cen1, cen2, sig1, sig2) ¶

The sum of 2 Gaussian function

sum_3NormalGaussian(x, amp1, amp2, amp3, cen1, cen2, cen3, sig1, sig2, sig3) ¶

The sum of 3 Gaussian function

sum_4NormalGaussian(x, amp1, amp2, amp3, amp4, cen1, cen2, cen3, cen4, sig1, sig2, sig3, sig4) ¶

The sum of 4 Gaussian function

sum_5NormalGaussian(x, amp1, amp2, amp3, amp4, amp5, cen1, cen2, cen3, cen4, cen5, sig1, sig2, sig3, sig4, sig5) ¶

The sum of 5 Gaussian function

DoseResp(x, A1=-3.3, A2=-2.9, LOGx0=480000, p=1.2) ¶

Dose-response curve with variable Hill slope given by parameter 'p'. Origin's Category: Pharmacology * Params: Names=A1,A2,LOGx0,p Meanings=bottom asymptote,top asymptote, center, hill slope Initiate params: pars = mod.make_params(A1=-3.3, A2=-2.9, LOGx0=480000, p=1.2)

BiDoseResp(x, A1=-3.3, A2=-2.9, LOGx01=175000, LOGx02=480000, h1=0.1, h2=0.2, p=0.5) ¶

Biphasic Dose Response Function, Origin's Category: Pharmacology * Params: Names=A1, A2, LOGx01, LOGx02, h1, h2, p Meanings=Bottom, Top, 1^st EC50, 2^nd EC50, slope1, slope2, proportion Initiate params: pars = mod.make_params(A1=0, A2=100, LOGx01=-8, LOGx02=-4, h1=0.8, h2=1.2, p=0.5)

Carreau(x, A1=60, A2=3, t=3.0, a=2.2, n=0.3) ¶

Carreau-Yasuda model to describe pseudoplastic flow with asymptotic viscosities at zero and infinite shear rates Origin's Category: Rheology * Params: Names = A1,A2,t,a,n >0 (lower bound) Meanings = zero shear viscosity,infinite shear viscosity,time constant,transition control factor,power index Initiate params: pars = mod.make_params(A1=-3.3, A2=-2.9, t=2.0, a=2.2, n=0.2)

Cross(x, A1=0.1, A2=3, t=1000, m=0.9) ¶

Cross model to describe pseudoplastic flow with asymptotic viscosities at zero and infinite shear rates Origin's Category: Rheology * Params: Names = A1,A2,t,m >0 (lower bound) Meanings = zero shear viscosity,infinite shear viscosity,time constant,power index Initiate params: pars = mod.make_params(A1=0.1, A2=3, t=1000, m=0.9)

GammaCFD(x, y0, A1, a, b) ¶

Gamma cumulative distribution function Origin's Category: Statistics * Params: Names = y0,A1,a,b (A1,a,b >0) Meanings = Offset,Amplitude,Shape,Scale Initiate params: pars = mod.make_params(A1=0.1, A2=3, t=1000, m=0.9)

`grid_box` ¶

Functions:

grid_box_2d –

devide box into 2d grid, return list of atom-IDs in each slab and list of slab-centers
grid_box_1d –

devide box into 1d slabs, return list of atom-IDs in each slab and list of slab-centers

`grid_box_2d(points, box, plane='XY', mode='bin_number', grid_size=[20, 20])` ¶

devide box into 2d grid, return list of atom-IDs in each slab and list of slab-centers Args: P : Nx3 array contain positions of atoms box : simulation box mode : "bin_number" or "bin_size" mode_value : corresponding 'Number-of-bins' or 'size-of-bin' plane : on which plane the box will be gridded Returns: atomIDinCell : 1xBinNumber array of 1xM-vector, contain indices of atoms of each Cell cellCenter : 1xBinNumber array of scalar, is center of each slab

`grid_box_1d(points, box, axis='Z', mode='bin_number', grid_size=20)` ¶

devide box into 1d slabs, return list of atom-IDs in each slab and list of slab-centers Args: P : Nx3 array contain positions of atoms box : simulation box mode : "bin_number" or "bin_size" mode_value : corresponding 'Number-of-bins' or 'size-of-bin' axis : on which axis the box will be slabbed Returns: atomIDinCell : 1xBinNumber array of 1xM arrays, contain indices of atoms of each Slab, array of arrays geoCenter : 1xBinNumber array of scalar, is geometry center of each slab massCenter : 1xBinNumber array of scalar, is mass center of each slab

`many_stuff` ¶

Functions:

memory_usage –

return the memory usage in MB
natSorted –

https://stackoverflow.com/questions/4836710/is-there-a-built-in-function-for-string-natural-sort
split_list –

Should use np.array_split instead
find_nearest_value –

`memory_usage()` ¶

return the memory usage in MB

`natSorted(mylist)` ¶

https://stackoverflow.com/questions/4836710/is-there-a-built-in-function-for-string-natural-sort

`split_list(a, n)` ¶

Should use np.array_split instead Args: a (list): list to be splitted n (int): number of chunks Returns: generator: a generator of splitted list

`find_nearest_value(array, value)` ¶

`row_operation` ¶

Functions:

unique_row –

find match_indices & mismatch_indices of arr(find_rows) in arr(X), return indices of X
match_row –

find match_indices & mismatch_indices of arr(find_rows) in arr(X), return indices of X
asvoid –

Base on: https://stackoverflow.com/questions/38674027/find-the-row-indexes-of-several-values-in-a-numpy-array

`unique_row(X, tol_decimal=2)` ¶

find match_indices & mismatch_indices of arr(find_rows) in arr(X), return indices of X Args: X, find_rows : NxN numpy arrays tol_decimal : number of digits for round off input data

`match_row(X, find_rows, tol_decimal=2)` ¶

find match_indices & mismatch_indices of arr(find_rows) in arr(X), return indices of X Args: X, find_rows : NxN numpy arrays tol_decimal : number of digits for round off input data

`asvoid(arr)` ¶

Base on: https://stackoverflow.com/questions/38674027/find-the-row-indexes-of-several-values-in-a-numpy-array

`string_index` ¶

Functions:

string_index –

groupSURF index by consecutive-series

`string_index(idx_list)` ¶

groupSURF index by consecutive-series

`unit` ¶

This module to convert unit of some physical properties pressure

Consider to use this module: https://unyt.readthedocs.io/en/stable/usage.html

Functions:

pressure –

convert unit of pressure
force –

convert unit of force
energy –

convert unit of energy
constant –

list of constants

`pressure(key_word='all_key')` ¶

convert unit of pressure Pa: Pascal atm: standard atmosphere at: technical atmosphere

kgf/cm2 = kg/cm2 1 Pa = 1 N/m^2 1 kgf/cm2 = 1

Parameters:

key_word (str, default: 'all_key' ) –

a string to specify units to be converted.

Returns:

factor ( float ) –

multiply factor of conversion

Examples:

key_word='Pa_atm': convert from Pa (Pascal) to atm (Standard atmosphere)

`force(key_word='all_key')` ¶

convert unit of force N: Newton kgf = m.g: kilogram-force (weight: one kilogram of mass in a 9.80665 m/s2 gravitational field) lbf: pound-force p: pond

1 N = 1 J/m (Work = Force.distance) 1 kcal = 4184 J = 4184 N.m = 4184.10^10 N.Angstrom 69.4786 pN = 1 kcal/mol Angstrom. https://tinyurl.com/yb2gnlhc

Parameters:

key_word (str, default: 'all_key' ) –

a string to specify units to be converted.

Returns:

factor ( float ) –

multiply factor of conversion

`energy(key_word='all_key')` ¶

convert unit of energy J: Joule W.h: watt-hour cal: calorie (th) hp.h: horsepower hour eV: electron-volt

1 J = 1 N.m (Work = Force.distance) 1J = 1 W.s

Parameters:

key_word (str, default: 'all_key' ) –

a string to specify units to be converted.

Returns:

factor ( float ) –

multiply factor of conversion

Notes

## convert eV to kcal/mol
eV2J = 1/unit_convert.energy('J_eV')
J2Jmol = unit_convert.constant('1/mol')
kj2kcal = 1/unit_convert.energy('kcal/mol_kJ/mol')
eV2kcalmol = eV2J * J2Jmol * 1e-3 *kj2kcal

`constant(key_word='all_key')` ¶

list of constants Na = 6.02214076e23 (=1/mol): Avogadro number

Parameters:

Ex –

key_word='Pa_atm': convert from Pa (Pascal) to atm (Standard atmosphere)

Returns: factor: float, multiply factor of conversion

API¶

thmd ¶

__description__ = 'Python package' module-attribute ¶

__long_description__ = 'Python package for pre- and post-processing data from MD simulations.' module-attribute ¶

__author__ = 'thangckt' module-attribute ¶

__version ¶

TYPE_CHECKING = False module-attribute ¶

VERSION_TUPLE = Tuple[Union[int, str], ...] module-attribute ¶

version: str = '0.1.dev533+g053f72e.d20250107' module-attribute ¶

__version__: str = '0.1.dev533+g053f72e.d20250107' module-attribute ¶

__version_tuple__: VERSION_TUPLE = (0, 1, 'dev533', 'g053f72e.d20250107') module-attribute ¶

version_tuple: VERSION_TUPLE = (0, 1, 'dev533', 'g053f72e.d20250107') module-attribute ¶

calc ¶

colvar ¶

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

r0 = r0 instance-attribute ¶

d0 = d0 instance-attribute ¶

n = n instance-attribute ¶

m = m instance-attribute ¶

dmax(dmax_tol) ¶

eval(x, compute_der=False) ¶

HEAVISIDE(r0) ¶

dmax = r0 instance-attribute ¶

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

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

CUBIC(d0, dmax) ¶

d0 = d0 instance-attribute ¶

dmax = dmax instance-attribute ¶

eval(x, compute_der=False) ¶

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

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

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

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

cv_CoordNum ¶

coord_number(Points, **kwargs) ¶

cv_PairANGLE ¶

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

cv_Steinhardt ¶

Ql_Steinhardt(ql_i) ¶

Local_Ql_Steinhardt(ql_i, qlm_j, SW) ¶

cv_fccCUBIC ¶

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

cv_localCRYSTALLINITY ¶

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

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

cv_slip_atom ¶

slip_atom(points, d_ref, d_diff) ¶

find_neighbor ¶

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_neighbor_list(points, box, bound_cond=[1, 1, 1], cutoff_neighbor=None, k_neighbor=None, k_cutoff=20, keep_ref=False) ¶

sph_harmonics ¶

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

switch_function ¶

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

HEAVISIDE(r0) ¶

CUBIC(d0, dmax) ¶

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

voronoi ¶

Voro3D() ¶

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

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

free_energy ¶

TImethod ¶

kB = sc.value('Boltzmann constant in eV/K') module-attribute ¶

AS_integration(forwFile, backFile) ¶

RS_integration(forwFile, backFile, T0=0, F0=0, kB=kB) ¶

Helmholtz_excess_UF(p, x) ¶

logmfd ¶

LogMFD_FCCUBIC(Element, Potential_Name, modelType='Bulk', zDirect='001') ¶

replica_logPD_integration(logmfd_files: list[str], replica_files: list[str], beta: float = 1.0) ¶

replica_MD_average(MD_out_files: list[str]) ¶

exp_normalize(x) ¶

`thmd` ¶

`description = 'Python package'` `module-attribute` ¶

`__long_description__ = 'Python package for pre- and post-processing data from MD simulations.'` `module-attribute` ¶

`author = 'thangckt'` `module-attribute` ¶

`__version` ¶

`TYPE_CHECKING = False` `module-attribute` ¶

`VERSION_TUPLE = Tuple[Union[int, str], ...]` `module-attribute` ¶

`version: str = '0.1.dev533+g053f72e.d20250107'` `module-attribute` ¶

`version: str = '0.1.dev533+g053f72e.d20250107'` `module-attribute` ¶

`__version_tuple__: VERSION_TUPLE = (0, 1, 'dev533', 'g053f72e.d20250107')` `module-attribute` ¶

`version_tuple: VERSION_TUPLE = (0, 1, 'dev533', 'g053f72e.d20250107')` `module-attribute` ¶

`calc` ¶

`colvar` ¶

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

`r0 = r0` `instance-attribute` ¶

`d0 = d0` `instance-attribute` ¶

`n = n` `instance-attribute` ¶

`m = m` `instance-attribute` ¶

`dmax(dmax_tol)` ¶

`eval(x, compute_der=False)` ¶

`HEAVISIDE(r0)` ¶

`dmax = r0` `instance-attribute` ¶

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

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

`CUBIC(d0, dmax)` ¶

`d0 = d0` `instance-attribute` ¶

`dmax = dmax` `instance-attribute` ¶

`eval(x, compute_der=False)` ¶

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

`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)` ¶

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

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

`cv_CoordNum` ¶

`coord_number(Points, **kwargs)` ¶

`cv_PairANGLE` ¶

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

`cv_Steinhardt` ¶

`Ql_Steinhardt(ql_i)` ¶

`Local_Ql_Steinhardt(ql_i, qlm_j, SW)` ¶

`cv_fccCUBIC` ¶

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

`cv_localCRYSTALLINITY` ¶

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

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

`cv_slip_atom` ¶

`slip_atom(points, d_ref, d_diff)` ¶

`find_neighbor` ¶

`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_neighbor_list(points, box, bound_cond=[1, 1, 1], cutoff_neighbor=None, k_neighbor=None, k_cutoff=20, keep_ref=False)` ¶

`sph_harmonics` ¶

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

`switch_function` ¶

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

`HEAVISIDE(r0)` ¶

`CUBIC(d0, dmax)` ¶

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

`voronoi` ¶

`Voro3D()` ¶

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

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

`free_energy` ¶

`TImethod` ¶

`kB = sc.value('Boltzmann constant in eV/K')` `module-attribute` ¶

`AS_integration(forwFile, backFile)` ¶

`RS_integration(forwFile, backFile, T0=0, F0=0, kB=kB)` ¶

`Helmholtz_excess_UF(p, x)` ¶

`logmfd` ¶

`LogMFD_FCCUBIC(Element, Potential_Name, modelType='Bulk', zDirect='001')` ¶

`replica_logPD_integration(logmfd_files: list[str], replica_files: list[str], beta: float = 1.0)` ¶

`replica_MD_average(MD_out_files: list[str])` ¶

`exp_normalize(x)` ¶

`read_df(file, engine='LAMMPS')` ¶