When using LAMMPS with Matlantis, you need the matlantis-lammps integration package to use PFP as the force field (installation guide here).
In this article, we have summarized how to create inputs for performing calculations with PFP in LAMMPS, after installing matlantis-lammps. Please also refer to the "lammps.in" file and/or README file included with the integration package after installation.
------
Input Format for using PFP Potential
- Lines starting with # are comment lines.
-
"units" should be "metal" (line 3).
- Unit conversion is required as needed. Please refer to the LAMMPS official manual for details: https://docs.lammps.org/units.html
- For example, note that the time unit for "metal" is in picoseconds.
- The Tdamp value (the third number) for the thermostat requires adjustment because it is in time units.
- When using "fix nvt temp 300.0 300.0 100.0" with units set to "real," it corresponds to "fix nvt temp 300.0 300.0 0.1" when using "metal."
-
The official manual can be found here: https://docs.lammps.org/fix_nh.html.
-
"atom_style" should be "atomic" (line 5).
- "charge" can also be used, it is synonymous with "atomic" since no charge is assigned.
- "molecular" and "full" can also be used.
- For more details about atom_style, please refer to the official manual: https://docs.lammps.org/atom_style.html.
- "pair_style" should be "pfp_api [half-width space] v.N.0.0 (PFP version) [half-width space] CRYSTAL_U0 (PFP calculation mode)" (line 18).
e.g.
pair_style pfp_api v7.0.0 CRYSTAL
- For LightPFP, use "light_pfp [half-width space] examplesv1hienca (model ID)."
e.g.
pair_style light_pfp examplesv1hienca
- For LightPFP, use "light_pfp [half-width space] examplesv1hienca (model ID)."
-
Specify "mass" for the constituent atoms (line 15, 16), and correspondingly, define the "pair_coeff" for the "species" with the constituent atoms (line 28).
- For example, if you define mass 3 16, you need to match the mass with pair_coeff like pair_coeff * * species C O O.
- Typically, LAMMPS requires a "data" file in addition to the "in" file. However, when using PFP, the settings for bonds, angles, dihedrals, and impropers, which are typically listed in the data file for classical force fields, are not needed. Therefore, "mass" and "pair_style" are specified in the "in"file. Of course, it is also possible to set "mass" and "pair_style" in the "data" file and execute without any issues.
Here is lammps.in
Creating LAMMPS Input File with ASE
This guide introduces how to create LAMMPS input files for PFP (Pairwise Force Field) using ASE, starting from a CIF file. Here, we will first write the structural information to a LAMMPS data file and then create an in file containing the necessary information to use the PFP potential.
Required Data: CIF file (here, test.cif)
Output Data: LAMMPS input files (here, test_opt.in, test.data)
1. Library Imports and Function Definition
We'll import the necessary libraries and define a function that writes a CIF file as a LAMMPS data file. This function will also calculate the molecular weight as a reference value.
The arguments for writing the LAMMPS data are set as follows:
- format:"lammps-data"
- masses :True
- units: "metal"
- atom_style: "atomic"
- specorder: None
import sys, os
import numpy as np
from ase.data import atomic_numbers
from ase.io import read, write, Trajectory
from ase.io.lammpsdata import read_lammps_data, write_lammps_data
from pfcc_extras.visualize import view_ngl
def cif_to_lammps_ase(cif_file, lammps_file, specorder=None, atom_style='atomic', units='metal'):
atoms = read(cif_file)
atoms.write(lammps_file,
format='lammps-data',
masses=True,
units=units,
atom_style=atom_style,
specorder=specorder
)
mol_weight = np.sum(atoms.get_masses())
return mol_weight
2. Output File Name Specification and Data File Creation
Specify the output file name and execute the function defined in step 1.
lammps_data_file = "test"
lammps_in_file = "test_opt"
mol_weight = cif_to_lammps_ase("IS(1).cif", lammps_data_file+".data")
print(mol_weight)
We'll visualize the loaded data file to confirm its contents. In this visualization, mask_slab is used to hide the elements that form the slab structure (here, "Rh" as an example).
atoms = read(lammps_data_file+'.data', format='lammps-data', units='metal')
mask_slab = atoms.numbers==atomic_numbers['Rh']
display(view_ngl([atoms, atoms[~mask_slab]], ['ball+stick'], replace_structure=True))
3. Load CIF File and Obtain "pair_coeff"
Load the CIF file as an ASE atoms object and get the values for "pair_coeff" input, as shown below.
ase_atoms = read("test.cif")
species = ase_atoms.symbols.species()
specorder = sorted(list(species), key=lambda x:-atomic_numbers[x])
pair_coeff = ' '.join(specorder)
print(species)
print(pair_coeff)
4. Writing the LAMMPS Input File
Write the LAMMPS input file.
lmp_data = f"./input_ase.data"
lmp_min_in = f"""
units metal
boundary p p p
atom_style atomic
atom_modify map yes
box tilt large
pair_style pfp_api v7.0.0 CRYSTAL_U0_PLUS_D3
read_data {lmp_data}
pair_coeff * * species {pair_coeff}
thermo 100
thermo_style custom step temp pe ke density
neigh_modify every 1 delay 0 check yes
minimize 1e-4 1e-6 10000 1000000
write_data lmp_min.data pair ij
"""
with open(lammps_in_file+".in", "w") as fo:
fo.write(lmp_min_in)
5. Running LAMMPS
!lmp_serial -i {lammps_in_file}.in -log ./opt_log.lammps