PyReaxFF for LAMMPS

pyreaxff.py

from lammps import PyLammps
from ase.data import atomic_numbers,atomic_masses

class PyReaxFF(object):
    def __init__(self,datafile="data.reaxc",atomtype=["C","H","O"],ffieldfilename = "ffield.reax.cho"):
        self.datafile=datafile
        self.atomtype=atomtype
        self.ffieldfilename = "ffield.reax.cho"

    def getforce(self,atoms):
        lmp = PyLammps()

        lmp.units("real")
        lmp.atom_style("charge")
        lmp.boundary("s s s")
        self.writedata(atoms)
        lmp.read_data(self.datafile)

        lmp.pair_style("reax/c lmp_control")
        lmp.pair_coeff("* * ", self.ffieldfilename, ' '.join(self.atomtype))
        lmp.fix(1, "all", "qeq/reax", 1, 0.0, 10.0, 1.0e-6, "reax/c")
        lmp.run(0)

        f = np.array(list(lmp.lmp.gather_atoms("f", 1, 3)))
        lmp.close()
        return f

    def writedata(self,atoms):
        buff=[]
        buff.extend(["LAMMPS PyReaxFF",f" {len(atoms)} atoms"," 0 bonds"," 0 dihedrals"," 0 impropers",f" {len(self.atomtype)} atom types"," 0 bond types"," 0 angle types"," 0 dihedral types"," 0 improper types"," -100 100 xlo xhi"," -100 100 ylo yhi"," -100 100 zlo zhi",""," Masses",""])
        buff.extend([f" {i} {atomic_masses[atomic_numbers[t]]} # {t}" for i,t in enumerate(self.atomtype,1)])
        buff.extend([""," Atoms # charge",""])
        buff.extend([f"{i} {self.atomtype.index(atom.symbol)+1} 0.0 {atom.position[0]} {atom.position[1]} {atom.position[2]} # {atom.symbol}" for i,atom in enumerate(atoms,1)])
        buff.append("")
        with open(self.datafile,'w') as f:
            f.write('\n'.join(buff))