Topology building for LAMMPS non-reactive MD simulations of amorphous carbon
This is an example of how to generate LAMMPS input files for simulations of water-lubricated amorphous carbon (a-C) interfaces with the non-reactive interatomic potential described in Reichenbach et al. starting from a set of atomic coordinates. The parameter file can be found in matscipy’s docs/topology/ directory. The example structure file aC_H2O.xyz (provided in the same directory) consists of two H/OH-terminated a-C slabs that are serparated by a couple of water molecules. Periodic boundary conditions are used in the interface plane.
from nglview import show_ase, ASEStructure
from ase.io import read
atoms = read('aC_H2O.xyz')
atoms.rotate(90, 'x')
view = show_ase(atoms)
component = view.add_component(ASEStructure(atoms), default_representation=False)
component.add_ball_and_stick(cylinderOnly=True, radiusType='covalent', radiusScale=0.6, aspectRatio=0.1)
component.add_ball_and_stick(radiusType='covalent', radiusScale=0.6)
view.control.zoom(0.8)
view
/usr/lib/python3/dist-packages/scipy/__init__.py:146: UserWarning: A NumPy version >=1.17.3 and <1.25.0 is required for this version of SciPy (detected version 1.26.4
warnings.warn(f"A NumPy version >={np_minversion} and <{np_maxversion}"
Assignment of the different atom types
The following script assigns an atom type for each atom based on its neighbours and generates a new structure file struct.extxyz that includes these atom types, the atomic coordinates and molecular ids, which can at a later stage faciliate group assignments in LAMMPS. The script only assigns atom types for which parameters are available in parameters_Reichenbach2023.in and gives an error otherwise. Please modify at your convenience if you use other parameter sets. Note that the assignment is based on the cutoffs specified when building the neighbour list: ('C', 'C'): 1.85, ('C', 'H'): 1.15, ('C', 'O'): 1.55, ('O', 'H'): 1.3. In most cases, it is desirable to choose these cutoffs in line with the cutoffs that are specified in the parameter file, which is used to generate the bonding topology later (here parameters_Reichenbach2023.in). The naming convention of the atom types follows Table 1 in Reichenbach et al.. Note that the script is rather strict and raises an error as soon as it encounters atoms that are in an unphysical or chemically not very stable environment such as C atoms with 5 neighbours or reactive sp-hybridised C on the surfaces. In these situations you may want to passivate the reactive atoms on the surfaces e.g. with H. Regarding C atoms with 5 C neighbours in the bulk, it may be acceptable to simply assign the CD atom type to them as long as you ensure that the bulk’s elastic response is not affected by this.
from ase.io import write, read
import numpy as np
from matscipy.neighbours import neighbour_list
from ase import Atoms
# read initial structure and initialize neighbour list
a = read('aC_H2O.xyz', index=0)
types = [0] * len(a)
# create neighbour list, note that the cutoffs here should be
# consistent with the ones in the parameter file.
i, j = neighbour_list('ij', a, {('C', 'C'): 1.85, ('C', 'H'): 1.15,
('C', 'O'): 1.55, ('O', 'H'): 1.3})
# define surface region to
# assign CA atom types to sp2-C surface atoms in the contact region
h_surface_region = 2.0 # arbitrarily defined depth of surface region
zmid_C = 65 # z coordinate separating upper and lower slab
mask_C = a.numbers == 6
mask_C_lower_slab = np.logical_and(mask_C, a.get_positions()[:, 2] < zmid_C)
mask_C_upper_slab = np.logical_and(mask_C, a.get_positions()[:, 2] > zmid_C)
zmax_C_lower_slab = np.max(a[mask_C_lower_slab].get_positions()[:, 2])
zmin_C_upper_slab = np.min(a[mask_C_upper_slab].get_positions()[:, 2])
CA_region_lower_slab = (zmax_C_lower_slab - h_surface_region,
zmax_C_lower_slab)
CA_region_upper_slab = (zmin_C_upper_slab,
zmin_C_upper_slab + h_surface_region)
def atom_in_CA_region(atom):
# returns True if the atom is in the surface region
zpos = atom.position[2]
return ((CA_region_lower_slab[0] < zpos <= CA_region_lower_slab[1]) or
(CA_region_upper_slab[0] <= zpos < CA_region_upper_slab[1]))
# assign C atom types
for k, atom in enumerate(a):
neighs = j[i == k] # all neighbor indices of atom k
if atom.symbol == 'C':
assert len(neighs) <= 4
# count neighbors by element
nC = 0
nO = 0
nH = 0
for neigh in neighs:
if a[neigh].symbol == 'C':
nC += 1
elif a[neigh].symbol == 'O':
nO += 1
elif a[neigh].symbol == 'H':
nH += 1
else:
raise NotImplementedError('unknown atom type.')
if len(neighs) == 4: # sp3-C atom types
assert nO <= 1 # only allow at most 1 OH group per C
if nC == 4:
types[k] = 'CD'
elif nC == 3:
if nO == 1:
types[k] = 'C2'
elif nH == 1:
types[k] = 'C1'
elif nC == 2:
assert nO + nH == 2
if nO == 1:
assert nH == 1
types[k] = 'C4'
elif nO == 0:
types[k] = 'C3'
elif nC == 1:
assert nO + nH == 3
if nO == 1:
assert nH == 2
types[k] = 'C6'
elif nO == 0:
types[k] = 'C5'
else:
raise NotImplementedError('unknown atom type.')
elif len(neighs) == 3: # sp2-C atom types
assert nC + nO + nH == 3
if nC == 3:
N_sp3neigh = 0 # count number of sp3 neighbors
for neigh in neighs:
neighs2 = j[i == neigh]
if len(neighs2) == 4:
N_sp3neigh += 1
if atom_in_CA_region(atom):
types[k] = 'CA'
else:
if N_sp3neigh == 0:
types[k] = 'CG' # bulk sp2 C with LJ interaction
else:
types[k] = 'CB' # bulk sp2 C, no LJ interaction
elif nC == 2:
if nO == 1:
types[k] = 'CY'
elif nH == 1:
types[k] = 'CZ'
else:
raise NotImplementedError('unknown atom type.')
elif nC == 1:
assert nH == 2
types[k] = 'CX'
# only allow bonds to CA, CB, CG, CY, CZ with 3 bonds
for neigh in neighs:
if a[neigh].symbol == 'C':
neighs2 = j[i == neigh]
if len(neighs2) != 3:
raise NotImplementedError('CX bonds to CA, CB, CG,'
' CY, CZ with 3 bonds')
elif len(neighs) == 2: # sp atoms
assert nC == 2
N_sp3neigh = 0
for neigh in neighs:
neighs2 = j[i == neigh]
if len(neighs2) == 4:
N_sp3neigh += 1
if atom_in_CA_region(atom):
print(k)
raise NotImplementedError('Sp-C atom in surface region')
else:
if N_sp3neigh == 0:
types[k] = 'CG' # bulk sp C with LJ interaction
else:
types[k] = 'CB' # bulk sp C, no LJ interaction
# assign all O atom types
for k, atom in enumerate(a):
neighs = j[i == k]
if atom.symbol == 'O':
assert len(neighs) == 2
nH = 0
for neigh in neighs:
sy = a[neigh].symbol
if sy == 'C':
if types[neigh] == 'CY':
types[k] = 'O2'
elif types[neigh] in ['C2', 'C4', 'C6']:
types[k] = 'O1'
elif sy == 'H':
nH += 1
if nH == 2:
types[k] = 'OW'
# assign all H atom types
for k, atom in enumerate(a):
neighs = j[i == k]
if atom.symbol == 'H':
assert len(neighs) == 1
for neigh in neighs:
if a[neigh].symbol == 'C':
if types[neigh] in ['CZ', 'CX']:
types[k] = 'H2'
elif types[neigh] in ['C1', 'C3', 'C4', 'C5', 'C6']:
types[k] = 'H1'
elif a[neigh].symbol == 'O':
if types[neigh] == 'O1':
types[k] = 'H4'
elif types[neigh] == 'O2':
types[k] = 'H5'
elif types[neigh] == 'OW':
types[k] = 'HW'
else:
raise NotImplementedError('unknown atom type.')
a.set_array('type', np.array(types))
# assign molecule ids
molid1 = np.ones(np.sum(a.get_positions()[:, 2] < 63)) # id=1 for lower slab
molid2 = np.ones(np.sum(a.get_positions()[:, 2] > 67)) + 1 # id=2 for upper slab
molid_H2O = np.arange(180) // 3 + 3 # H2O ids from 3 to 63
molid = np.append(molid1, molid2)
molid = np.append(molid, molid_H2O)
a.set_array('molid', molid.astype(int))
# save structure to .extxyz file
write('struct.extxyz', a)
# check for missed atoms
for k, t in enumerate(types):
if t == 0:
neighs = j[i == k]
print(k, a[k].symbol, [(neigh, types[neigh]) for neigh in neighs])
assert 0 not in types
Generation of the bonding topology and LAMMPS input files
In the next step, we use matscipy’s routines to generate the bonding topology and create LAMMPS input files based on the parameters specified in parameters_Reichenbach2023.in.
This script generates three files:
struct.lammps.atoms, which contains the atomic structure and the bonding topology in a LAMMPS-readable way.struct.lammps.opls, which contains the pair_styles and force field parameters in a LAMMPS-readable way.struct.lammps.in, which is an example LAMMPS input script how to run a geometry optimisation. Note that the script gives warnings of missing Dihedrals that are not part of the force field described in Reichenbach et al. (e.g., C-C-C-C Dihedrals).
%%capture
import matscipy
import matscipy.opls
import matscipy.io.opls
parameter_file = 'parameters_Reichenbach2023.in'
s = matscipy.io.opls.read_extended_xyz('struct.extxyz')
parameter_data = matscipy.io.opls.read_parameter_file(parameter_file)
cutoffs, atom_data, bond_data, angle_data, dihed_data = parameter_data
s.set_cutoffs(cutoffs)
s.set_atom_data(atom_data)
assert abs(sum(s.get_charges())) < 0.01 # ensure charge neutrality
s.get_bonds(bond_data)
s.get_angles(angle_data)
s.get_dihedrals(dihed_data)
matscipy.io.opls.write_lammps('struct.lammps', s);
To remove all files created during this tutorial, you can use the following code.
import os
for tmp_file in ['struct.extxyz', 'struct.lammps.atoms',
'struct.lammps.opls', 'struct.lammps.in']:
try:
os.remove(tmp_file)
except FileNotFoundError:
continue