Opls cleanup

examples/opls/atoms_to_lammps.py

@@ -0,0 +1,122 @@
+#!/usr/bin/env python3
+#
+# Copyright 2023 Andreas Klemenz (Fraunhofer IWM)
+#
+# matscipy - Materials science with Python at the atomic-scale
+# https://github.com/libAtoms/matscipy
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 2 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+
+
+# This script demonstrates how to create an atomic configuration in a
+# python script and generate input files for LAMMPS from it.
+# Alternatively, atomic configurations can be read from an extended xyz
+# file. See example 'extxyz_to_lammps.py'.
+
+
+import numpy as np
+import ase
+import matscipy.opls
+import matscipy.io.opls
+
+
+# Simple example: ethane molecules
+# 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
+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 = [10., 10., 10.]
+    )
+a1.translate([0., 3., 0.])
+
+# Alternative: 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 = [10., 10., 10.]
+    )
+a2 = matscipy.opls.OPLSStructure(a2)
+
+
+a = matscipy.opls.OPLSStructure(cell=[10., 10., 10.])
+a.extend(a1)
+a.extend(a2)
+a.center()
+
+# Specify atomic types. Notice the difference between type and element.
+# In this example we are using 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, all types of pair, angle and
+# dihedral interactions must be specified manually. Usually this means
+# searching the literature for suitable parameters, which can be a
+# tedious task. If it is known which interactions are present in a
+# system, this can be much easier. 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.
+cutoffs = matscipy.io.opls.read_cutoffs('cutoffs.in')
+a.set_cutoffs(cutoffs)
+
+bond_types, _ = a.get_bonds()
+print('Pairwise interactions:')
+print(bond_types)
+
+angle_types, _ = a.get_angles()
+print('\nAngular interactions:')
+print(angle_types)
+
+dih_types, _ = a.get_dihedrals()
+print('\nDihedral interactions:')
+print(dih_types)
+
+
+
+# Once the parameters of all interactions are known, they can be written
+# to a file. This 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)
+
+a.get_bonds(bond_data)
+a.get_angles(angle_data)
+a.get_dihedrals(dih_data)
+
+# 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. Modify this file for more
+# complex simulations.
+matscipy.io.opls.write_lammps('example', a)
+

examples/opls/cutoffs.in

@@ -0,0 +1,7 @@
+# Cutoffs 
+C1-C1 1.85
+C1-H1 1.15
+C1-H2 1.15
+H1-H1 0.0
+H2-H2 0.0
+H1-H2 0.0

examples/opls/ethane.extxyz

@@ -0,0 +1,18 @@
+16
+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
+C        5.5       6.5       5.0       1        C1
+H        3.5       6.5       5.0       1        H1
+H        4.5       7.5       5.0       1        H1
+H        4.5       5.5       5.0       1        H1
+H        5.5       6.5       6.0       1        H2
+H        5.5       6.5       4.0       1        H2
+H        6.5       6.5       5.0       1        H2
+C        4.5       3.5       5.0       2        C1
+C        5.5       3.5       5.0       2        C1
+H        3.5       3.5       5.0       2        H1
+H        4.5       4.5       5.0       2        H1
+H        4.5       2.5       5.0       2        H1
+H        5.5       3.5       6.0       2        H2
+H        5.5       3.5       4.0       2        H2
+H        6.5       3.5       5.0       2        H2

examples/opls/extxyz_to_lammps.py

@@ -0,0 +1,81 @@
+#!/usr/bin/env python3
+#
+# Copyright 2023 Andreas Klemenz (Fraunhofer IWM)
+#
+# matscipy - Materials science with Python at the atomic-scale
+# https://github.com/libAtoms/matscipy
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 2 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+
+
+# This script demonstrates how to read an atomic configuration from an
+# extended xyz file and create input files for LAMMPS from it.
+# Alternatively, atomic configurations can be created directly in a
+# python script. See example 'atoms_to_lammps.py'.
+
+
+import matscipy.io.opls
+
+# Read the atomic configuration from an Extended XYZ file with labeled
+# atoms. The file should contain the following columns: element (1 or 2
+# characters), x(float), y(float), z (float), molecule id (int), name
+# (1 or 2 characters). See the file "ethane.extxyz" for an example.
+# A full description of the extended xyz format can be found for example
+# in the ASE documentation.
+a = matscipy.io.opls.read_extended_xyz('ethane.extxyz')
+
+# To perform a non-reactive simulation, all types of pair, angle and
+# dihedral interactions must be specified manually. Usually this means
+# searching the literature for suitable parameters, which can be a
+# tedious task. If it is known which interactions are present in a
+# system, this can be much easier. 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.
+cutoffs = matscipy.io.opls.read_cutoffs('cutoffs.in')
+a.set_cutoffs(cutoffs)
+
+bond_types, _ = a.get_bonds()
+print('Pairwise interactions:')
+print(bond_types)
+
+angle_types, _ = a.get_angles()
+print('\nAngular interactions:')
+print(angle_types)
+
+dih_types, _ = a.get_dihedrals()
+print('\nDihedral interactions:')
+print(dih_types)
+
+
+
+# Once the parameters of all interactions are known, they can be written
+# to a file. This 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)
+
+a.get_bonds(bond_data)
+a.get_angles(angle_data)
+a.get_dihedrals(dih_data)
+
+# 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. Modify this file for more
+# complex simulations.
+matscipy.io.opls.write_lammps('example', a)
+

examples/opls/parameters.in

@@ -0,0 +1,33 @@
+# the blocks are separated by empty lines
+# comments are allowed
+#
+# All potential parameters in this file are for illustrative
+# purposes only! Do not use them for real simulations!
+#
+# Lennard-Jones and Charges
+# Element Epsilon(eV) Sigma(A) Charge(e)
+C1 0.0028 3.50 -0.06
+H1 0.0013 2.50  0.06
+H2 0.0013 2.42  0.12
+
+# Bonds (style harmonic)
+C1-C1 13.4 1.5
+C1-H1 14.3 1.1
+C1-H2  9.0 1.2
+
+# Angles (style harmonic)
+C1-C1-H1 1.5 110.0
+C1-C1-H2 1.1 115.0
+H1-C1-H1 0.7 110.0
+H2-C1-H2 0.7 110.0
+
+# Dihedrals 
+H1-C1-C1-H2 0.0 0.0 0.016 0.0
+
+# Cutoffs 
+C1-C1 1.85
+C1-H1 1.15
+C1-H2 1.15
+H1-H1 0.0
+H2-H2 0.0
+H1-H2 0.0