Setups for non-reactive simulations with LAMMPS
This notebook shows how to create an atomic configuration in a python script and how to generate input files for LAMMPS from it. As a simple example, we will set up ethane molecules and demonstrate the basic functionality of the matscipy.opls module.
matscipy.opls.OPLSStructure is a subclass of ase.Atoms. OPLSStructure objects can therefore be constructed and manipulated in the same way as ase.Atoms objects. The full documentation can be found in the ASE documentation.
import numpy as np
import ase
import ase.visualize
import matscipy.opls
a1 = matscipy.opls.OPLSStructure(
'C2H6',
positions = [[1., 0., 0.],
[2., 0., 0.],
[0., 0., 0.],
[1., 1., 0.],
[1., -1., 0.],
[2., 0., 1.],
[2., 0., -1.],
[3., 0., 0.]],
cell = [5., 5., 5.]
)
ase.visualize.view(a1, viewer='ngl')
/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}"
Alternatively, we can construct an ase.Atoms object and convert it to a matscipy.opls.OPLSStructure object:
a2 = ase.Atoms(
'C2H6',
positions = [[1., 0., 0.],
[2., 0., 0.],
[0., 0., 0.],
[1., 1., 0.],
[1., -1., 0.],
[2., 0., 1.],
[2., 0., -1.],
[3., 0., 0.]],
cell = [5., 5., 5.]
)
a2 = matscipy.opls.OPLSStructure(a2)
We can combine the two structures as we can for ase.Atoms objects.
a = matscipy.opls.OPLSStructure(cell=[10., 10., 10.])
a1.translate([0., 4., 0.])
a.extend(a1)
a.extend(a2)
a.center()
ase.visualize.view(a, viewer='ngl')
Next, we specify atomic types. Note the difference between type and element. In this example we use two different types of hydrogen atoms.
a.set_types(['C1', 'C1', 'H1', 'H1', 'H1', 'H2', 'H2', 'H2',
'C1', 'C1', 'H1', 'H1', 'H1', 'H2', 'H2', 'H2'])
To perform a non-reactive simulation, manual specification of all pairs, angles, and dihedrals is required. Typically this involves searching the literature for appropriate parameters. This can be a tedious task. However, if the interactions present in the system are already known, this process can be greatly simplified. Lists of all existing interactions can be generated based on the distance of the atoms from each other. The maximum distances up to which two atoms are considered to interact can be read from a file which will typically look like this:
# Cutoffs
C1-C1 1.85
C1-H1 1.15
...
Here, we read these cutoffs from the file cutoffs.in:
import matscipy.io.opls
cutoffs = matscipy.io.opls.read_cutoffs('cutoffs.in')
a.set_cutoffs(cutoffs)
bond_types, _ = a.get_bonds()
print(bond_types)
['C1-C1' 'C1-H1' 'C1-H2']
The same procedure applies to angles and dihedrals:
angle_types, _ = a.get_angles()
print(angle_types)
['C1-C1-H1', 'H1-C1-H1', 'C1-C1-H2', 'H2-C1-H2']
dih_types, _ = a.get_dihedrals()
print(dih_types)
['H1-C1-C1-H2']
Once all interaction parameters are known, they can be written to a file that looks like this:
# Element
C1 0.0028 3.50 -0.06
H1 0.0013 2.50 0.06
...
# Bonds
C1-C1 13.4 1.5
C1-H1 14.3 1.1
...
# Angles
C1-C1-H1 1.5 110.0
C1-C1-H2 1.1 115.0
...
# Dihedrals
H1-C1-C1-H2 0.0 0.0 0.016 0.0
...
# Cutoffs
C1-C1 1.85
...
Such a file can be used to generate the lists of all interactions and to create input files for LAMMPS:
cutoffs, atom_data, bond_data, angle_data, dih_data = matscipy.io.opls.read_parameter_file('parameters.in')
a.set_cutoffs(cutoffs)
a.set_atom_data(atom_data)
After reading in all parameters, we can construct the lists of bonds, angles and dihedrals:
a.get_bonds(bond_data)
(array(['C1-C1', 'C1-H1', 'C1-H2'], dtype=object),
array([[ 0, 1, 0],
[ 0, 9, 8],
[ 1, 2, 0],
[ 1, 3, 0],
[ 1, 4, 0],
[ 1, 10, 8],
[ 1, 11, 8],
[ 1, 12, 8],
[ 2, 5, 1],
[ 2, 6, 1],
[ 2, 7, 1],
[ 2, 13, 9],
[ 2, 14, 9],
[ 2, 15, 9]]))
a.get_angles(angle_data)
(['C1-C1-H1', 'H1-C1-H1', 'C1-C1-H2', 'H2-C1-H2'],
[[0, 2, 0, 1],
[1, 2, 0, 4],
[0, 1, 0, 4],
[1, 2, 0, 3],
[0, 1, 0, 3],
[1, 4, 0, 3],
[2, 0, 1, 6],
[2, 0, 1, 7],
[3, 6, 1, 7],
[2, 0, 1, 5],
[3, 6, 1, 5],
[3, 7, 1, 5],
[0, 10, 8, 9],
[1, 10, 8, 12],
[0, 9, 8, 12],
[1, 10, 8, 11],
[0, 9, 8, 11],
[1, 12, 8, 11],
[2, 8, 9, 14],
[2, 8, 9, 15],
[3, 14, 9, 15],
[2, 8, 9, 13],
[3, 14, 9, 13],
[3, 15, 9, 13]])
a.get_dihedrals(dih_data)
(['H1-C1-C1-H2'],
[[0, 2, 0, 1, 6],
[0, 2, 0, 1, 7],
[0, 2, 0, 1, 5],
[0, 4, 0, 1, 6],
[0, 4, 0, 1, 7],
[0, 4, 0, 1, 5],
[0, 3, 0, 1, 6],
[0, 3, 0, 1, 7],
[0, 3, 0, 1, 5],
[0, 10, 8, 9, 14],
[0, 10, 8, 9, 15],
[0, 10, 8, 9, 13],
[0, 12, 8, 9, 14],
[0, 12, 8, 9, 15],
[0, 12, 8, 9, 13],
[0, 11, 8, 9, 14],
[0, 11, 8, 9, 15],
[0, 11, 8, 9, 13]])
We have to call all these methods manually. The constructed lists are automatically stored in the matscipy.OPLSStructure-object, so there is no need to capture the return values at this point. Note that the construction of these lists can be very time-consuming for complex systems, especially the construction of the dihedrals list.
Now we can write the atomic structure, the potential definitions and a sample input script for LAMMPS (3 files in total). The input script contains a simple relaxation of the atomic position. You can use this file as a starting point to write scripts for more complex simulations.
matscipy.io.opls.write_lammps('example', a)
Alternatively, atomic configurations can be read from any type of file ASE can read. Extended xyz-files are particularly useful for this purpose since they are human-readable and can contain information about the atom type. A typical extended XYZ file for non-reactive simulations will look like this
42
Lattice="10.0 0.0 0.0 0.0 10.0 0.0 0.0 0.0 10.0" Properties=species:S:1:pos:R:3:molid:I:1:type:S:1 pbc="F F F"
C 4.5 6.5 5.0 1 C1
H 4.5 5.5 5.0 1 H1
H 5.5 3.5 4.0 2 H2
...
The file should include the following columns: element (1 or 2 characters), x(float), y(float), z (float), molecule id (int), type (1 or 2 characters). A full description of the extended XYZ format can be found in the ASE documentation.
b = matscipy.io.opls.read_extended_xyz('ethane.extxyz')
ase.visualize.view(b, viewer='ngl')