Skip to content

model

thmd.model

Modules:

D1tube

Functions:

  • lattice_CNT

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

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

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

coord_rotation

Classes:

  • CoordTransform

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

Functions:

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"
  1. Bower, Allan F. Applied Mechanics of Solids. CRC Press, 2009. page 711
  2. https://link.aps.org/doi/10.1103/PhysRevB.92.180102
  3. 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 = 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
  1. The polar(theta) angle defined from from Z-axis down (zero at the North pole to 180° at the South pole)
  2. adapted from this

  3. 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 = [-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:

Attributes:

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
  1. Due to mbuild cannot be installed with python 3.10, so import this package in functions to avoid checking in thmd

Functions:

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