Source code for gmso.formats.lammpsdata

"""Read and write LAMMPS data files."""
from __future__ import division

import datetime
import warnings

import numpy as np
import unyt as u
from sympy import sympify
from unyt.array import allclose_units

from gmso.core.angle import Angle
from gmso.core.angle_type import AngleType
from gmso.core.atom import Atom
from gmso.core.atom_type import AtomType
from gmso.core.bond import Bond
from gmso.core.bond_type import BondType
from gmso.core.box import Box
from gmso.core.element import element_by_mass
from gmso.core.topology import Topology
from gmso.lib.potential_templates import PotentialTemplateLibrary
from gmso.utils.conversions import (
    convert_opls_to_ryckaert,
    convert_ryckaert_to_opls,
)


[docs]def write_lammpsdata(topology, filename, atom_style="full"): """Output a LAMMPS data file. Outputs a LAMMPS data file in the 'full' atom style format. Assumes use of 'real' units. See http://lammps.sandia.gov/doc/atom_style.html for more information on atom styles. Parameters ---------- Topology : `Topology` A Topology Object filename : str Path of the output file atom_style : str, optional, default='full' Defines the style of atoms to be saved in a LAMMPS data file. The following atom styles are currently supported: 'full', 'atomic', 'charge', 'molecular' see http://lammps.sandia.gov/doc/atom_style.html for more information on atom styles. Notes ----- See http://lammps.sandia.gov/doc/2001/data_format.html for a full description of the LAMMPS data format. This is a work in progress, as only atoms, masses, and atom_type information can be written out. Some of this function has been adopted from `mdtraj`'s support of the LAMMPSTRJ trajectory format. See https://github.com/mdtraj/mdtraj/blob/master/mdtraj/formats/lammpstrj.py for details. """ if atom_style not in ["atomic", "charge", "molecular", "full"]: raise ValueError( 'Atom style "{}" is invalid or is not currently supported'.format( atom_style ) ) # TODO: Support various unit styles box = topology.box with open(filename, "w") as data: data.write( "{} written by topology at {}\n\n".format( topology.name if topology.name is not None else "", str(datetime.datetime.now()), ) ) data.write("{:d} atoms\n".format(topology.n_sites)) if atom_style in ["full", "molecular"]: if topology.n_bonds != 0: data.write("{:d} bonds\n".format(topology.n_bonds)) else: data.write("0 bonds\n") if topology.n_angles != 0: data.write("{:d} angles\n".format(topology.n_angles)) else: data.write("0 angles\n") if topology.n_dihedrals != 0: data.write("{:d} dihedrals\n\n".format(topology.n_dihedrals)) else: data.write("0 dihedrals\n\n") data.write("\n{:d} atom types\n".format(len(topology.atom_types))) data.write("{:d} bond types\n".format(len(topology.bond_types))) data.write("{:d} angle types\n".format(len(topology.angle_types))) data.write("{:d} dihedral types\n".format(len(topology.dihedral_types))) data.write("\n") # Box data if allclose_units( box.angles, u.unyt_array([90, 90, 90], "degree"), rtol=1e-5, atol=1e-8, ): warnings.warn("Orthorhombic box detected") box.lengths.convert_to_units(u.angstrom) for i, dim in enumerate(["x", "y", "z"]): data.write( "{0:.6f} {1:.6f} {2}lo {2}hi\n".format( 0, box.lengths.value[i], dim ) ) else: warnings.warn("Non-orthorhombic box detected") box.lengths.convert_to_units(u.angstrom) box.angles.convert_to_units(u.radian) vectors = box.get_vectors() a, b, c = box.lengths alpha, beta, gamma = box.angles lx = a xy = b * np.cos(gamma) xz = c * np.cos(beta) ly = np.sqrt(b ** 2 - xy ** 2) yz = (b * c * np.cos(alpha) - xy * xz) / ly lz = np.sqrt(c ** 2 - xz ** 2 - yz ** 2) xhi = vectors[0][0] yhi = vectors[1][1] zhi = vectors[2][2] xy = vectors[1][0] xz = vectors[2][0] yz = vectors[2][1] xlo = u.unyt_array(0, xy.units) ylo = u.unyt_array(0, xy.units) zlo = u.unyt_array(0, xy.units) xlo_bound = xlo + u.unyt_array( np.min([0.0, xy, xz, xy + xz]), xy.units ) xhi_bound = xhi + u.unyt_array( np.max([0.0, xy, xz, xy + xz]), xy.units ) ylo_bound = ylo + u.unyt_array(np.min([0.0, yz]), xy.units) yhi_bound = yhi + u.unyt_array(np.max([0.0, yz]), xy.units) zlo_bound = zlo zhi_bound = zhi data.write( "{0:.6f} {1:.6f} xlo xhi\n".format( xlo_bound.value, xhi_bound.value ) ) data.write( "{0:.6f} {1:.6f} ylo yhi\n".format( ylo_bound.value, yhi_bound.value ) ) data.write( "{0:.6f} {1:.6f} zlo zhi\n".format( zlo_bound.value, zhi_bound.value ) ) data.write( "{0:.6f} {1:.6f} {2:.6f} xy xz yz\n".format( xy.value, xz.value, yz.value ) ) # TODO: Get a dictionary of indices and atom types if topology.is_typed(): # Write out mass data data.write("\nMasses\n\n") for atom_type in topology.atom_types: data.write( "{:d}\t{:.6f}\t# {}\n".format( topology.atom_types.index(atom_type) + 1, atom_type.mass.in_units(u.g / u.mol).value, atom_type.name, ) ) # TODO: Modified cross-interactions # Pair coefficients data.write("\nPair Coeffs # lj\n\n") for idx, param in enumerate(topology.atom_types): data.write( "{}\t{:.5f}\t{:.5f}\n".format( idx + 1, param.parameters["epsilon"] .in_units(u.Unit("kcal/mol")) .value, param.parameters["sigma"].in_units(u.angstrom).value, ) ) if topology.bonds: data.write("\nBond Coeffs\n\n") for idx, bond_type in enumerate(topology.bond_types): data.write( "{}\t{:.5f}\t{:.5f}\n".format( idx + 1, bond_type.parameters["k"] .in_units(u.Unit("kcal/mol/angstrom**2")) .value / 2, bond_type.parameters["r_eq"] .in_units(u.Unit("angstrom")) .value, ) ) if topology.angles: data.write("\nAngle Coeffs\n\n") for idx, angle_type in enumerate(topology.angle_types): data.write( "{}\t{:.5f}\t{:.5f}\n".format( idx + 1, angle_type.parameters["k"] .in_units(u.Unit("kcal/mol/radian**2")) .value / 2, angle_type.parameters["theta_eq"] .in_units(u.Unit("degree")) .value, ) ) # TODO: Write out multiple dihedral styles if topology.dihedrals: data.write("\nDihedral Coeffs\n\n") for idx, dihedral_type in enumerate(topology.dihedral_types): rbtorsion = PotentialTemplateLibrary()[ "RyckaertBellemansTorsionPotential" ] if ( dihedral_type.expression == sympify(rbtorsion.expression) or dihedral_type.name == rbtorsion.name ): dihedral_type = convert_ryckaert_to_opls(dihedral_type) data.write( "{}\t{:.5f}\t{:5f}\t{:5f}\t{:.5f}\n".format( idx + 1, dihedral_type.parameters["k1"] .in_units(u.Unit("kcal/mol")) .value, dihedral_type.parameters["k2"] .in_units(u.Unit("kcal/mol")) .value, dihedral_type.parameters["k3"] .in_units(u.Unit("kcal/mol")) .value, dihedral_type.parameters["k4"] .in_units(u.Unit("kcal/mol")) .value, ) ) # Atom data data.write("\nAtoms\n\n") if atom_style == "atomic": atom_line = "{index:d}\t{type_index:d}\t{x:.6f}\t{y:.6f}\t{z:.6f}\n" elif atom_style == "charge": atom_line = "{index:d}\t{type_index:d}\t{charge:.6f}\t{x:.6f}\t{y:.6f}\t{z:.6f}\n" elif atom_style == "molecular": atom_line = "{index:d}\t{zero:d}\t{type_index:d}\t{x:.6f}\t{y:.6f}\t{z:.6f}\n" elif atom_style == "full": atom_line = "{index:d}\t{zero:d}\t{type_index:d}\t{charge:.6f}\t{x:.6f}\t{y:.6f}\t{z:.6f}\n" for i, site in enumerate(topology.sites): data.write( atom_line.format( index=topology.sites.index(site) + 1, type_index=topology.atom_types.index(site.atom_type) + 1, zero=0, charge=site.charge.to(u.elementary_charge).value, x=site.position[0].in_units(u.angstrom).value, y=site.position[1].in_units(u.angstrom).value, z=site.position[2].in_units(u.angstrom).value, ) ) if topology.bonds: data.write("\nBonds\n\n") for i, bond in enumerate(topology.bonds): data.write( "{:d}\t{:d}\t{:d}\t{:d}\n".format( i + 1, topology.bond_types.index(bond.connection_type) + 1, topology.sites.index(bond.connection_members[0]) + 1, topology.sites.index(bond.connection_members[1]) + 1, ) ) if topology.angles: data.write("\nAngles\n\n") for i, angle in enumerate(topology.angles): data.write( "{:d}\t{:d}\t{:d}\t{:d}\t{:d}\n".format( i + 1, topology.angle_types.index(angle.connection_type) + 1, topology.sites.index(angle.connection_members[0]) + 1, topology.sites.index(angle.connection_members[1]) + 1, topology.sites.index(angle.connection_members[2]) + 1, ) ) if topology.dihedrals: data.write("\nDihedrals\n\n") for i, dihedral in enumerate(topology.dihedrals): data.write( "{:d}\t{:d}\t{:d}\t{:d}\t{:d}\t{:d}\n".format( i + 1, topology.dihedral_types.index(dihedral.connection_type) + 1, topology.sites.index(dihedral.connection_members[0]) + 1, topology.sites.index(dihedral.connection_members[1]) + 1, topology.sites.index(dihedral.connection_members[2]) + 1, topology.sites.index(dihedral.connection_members[3]) + 1, ) )
def read_lammpsdata( filename, atom_style="full", unit_style="real", potential="lj" ): """Read in a lammps data file as a GMSO topology. Parameters ---------- filename : str LAMMPS data file atom_style : str, optional, default='full' Inferred atom style defined by LAMMPS potential: str, optional, default='lj' Potential type defined in data file Returns ------- top : GMSO Topology A GMSO Topology object Notes ----- See http://lammps.sandia.gov/doc/2001/data_format.html for a full description of the LAMMPS data format. This is a work in progress, as only several atom styles, potential styles, and unit styes are currently supported. Currently supporting the following atom styles: 'full' Currently supporting the following unit styles: 'real' Currently supporting the following potential styles: 'lj' Proper dihedrals can be read in but is currently not tested. Currently not supporting improper dihedrals. """ # TODO: Add argument to ask if user wants to infer bond type top = Topology() # Validate 'atom_style' if atom_style not in ["full"]: raise ValueError( 'Atom Style "{}" is invalid or is not currently supported'.format( atom_style ) ) # Validate 'unit_style' if unit_style not in ["real"]: raiseValueError( 'Unit Style "{}" is invalid or is not currently supported'.format( unit_style ) ) # Parse box information _get_box_coordinates(filename, unit_style, top) # Parse atom type information top, type_list = _get_ff_information(filename, unit_style, top) # Parse atom information _get_atoms(filename, top, unit_style, type_list) # Parse connection (bonds, angles, dihedrals) information # TODO: Add more atom styles if atom_style in ["full"]: _get_connection(filename, top, unit_style, connection_type="bond") _get_connection(filename, top, unit_style, connection_type="angle") top.update_topology() return top def get_units(unit_style): """Get units for specific LAMMPS unit style.""" # Need separate angle units for harmonic force constant and angle unit_style_dict = { "real": { "mass": u.g, "distance": u.angstrom, "energy": u.kcal / u.mol, "angle_k": u.radian, "angle": u.degree, "charge": u.elementary_charge, } } return unit_style_dict[unit_style] def _get_connection(filename, topology, unit_style, connection_type): """Parse connection types.""" with open(filename, "r") as lammps_file: types = False for i, line in enumerate(lammps_file): if connection_type in line.split(): n_connection_types = int(line.split()[0]) types = True if connection_type.capitalize() in line.split(): break if types == False: return topology connection_type_lines = open(filename, "r").readlines()[ i + 2 : i + n_connection_types + 2 ] connection_type_list = list() for line in connection_type_lines: if connection_type == "bond": c_type = BondType(name=line.split()[0]) # Multiply 'k' by 2 since LAMMPS includes 1/2 in the term c_type.parameters["k"] = ( float(line.split()[1]) * u.Unit( get_units(unit_style)["energy"] / get_units(unit_style)["distance"] ** 2 ) * 2 ) c_type.parameters["r_eq"] = float(line.split()[2]) * ( get_units(unit_style)["distance"] ** 2 ) elif connection_type == "angle": c_type = AngleType(name=line.split()[0]) # Multiply 'k' by 2 since LAMMPS includes 1/2 in the term c_type.parameters["k"] = ( float(line.split()[1]) * u.Unit( get_units(unit_style)["energy"] / get_units(unit_style)["angle_k"] ** 2 ) * 2 ) c_type.parameters["theta_eq"] = float(line.split()[2]) * u.Unit( get_units(unit_style)["angle"] ) connection_type_list.append(c_type) with open(filename, "r") as lammps_file: for i, line in enumerate(lammps_file): if connection_type + "s" in line.split(): n_connections = int(line.split()[0]) if connection_type.capitalize() + "s" in line.split(): break connection_lines = open(filename, "r").readlines()[ i + 2 : i + n_connections + 2 ] # Determine number of sites to generate if connection_type == "bond": n_sites = 2 elif connection_type == "angle": n_sites = 3 else: n_sites = 4 for i, line in enumerate(connection_lines): site_list = list() for j in range(n_sites): site = topology.sites[int(line.split()[j + 2]) - 1] site_list.append(site) if connection_type == "bond": connection = Bond( connection_members=site_list, bond_type=connection_type_list[int(line.split()[1]) - 1], ) elif connection_type == "angle": connection = Angle( connection_members=site_list, angle_type=connection_type_list[int(line.split()[1]) - 1], ) topology.add_connection(connection) return topology def _get_atoms(filename, topology, unit_style, type_list): """Parse the atom information in the LAMMPS data file.""" with open(filename, "r") as lammps_file: for i, line in enumerate(lammps_file): if "atoms" in line.split(): n_atoms = int(line.split()[0]) if "Atoms" in line.split(): break atom_lines = open(filename, "r").readlines()[i + 2 : i + n_atoms + 2] for line in atom_lines: atom_line = line.split() atom_type = atom_line[2] charge = u.unyt_quantity( float(atom_line[3]), get_units(unit_style)["charge"] ) coord = u.angstrom * u.unyt_array( [float(atom_line[4]), float(atom_line[5]), float(atom_line[6])] ) site = Atom( charge=charge, position=coord, atom_type=type_list[int(atom_type) - 1], ) element = element_by_mass(site.atom_type.mass.value) site.name = element.name site.element = element topology.add_site(site) topology.update_sites() return topology def _get_box_coordinates(filename, unit_style, topology): """Parse box information.""" with open(filename, "r") as lammps_file: for line in lammps_file: if "xlo" in line.split(): break x_line = line.split() y_line = lammps_file.readline().split() z_line = lammps_file.readline().split() x = float(x_line[1]) - float(x_line[0]) y = float(y_line[1]) - float(y_line[0]) z = float(z_line[1]) - float(z_line[0]) # Check if box is triclinic tilts = lammps_file.readline().split() if "xy" in tilts: xy = float(tilts[0]) xz = float(tilts[1]) yz = float(tilts[2]) xhi = float(x_line[1]) - np.max([0.0, xy, xz, xy + xz]) xlo = float(x_line[0]) - np.min([0.0, xy, xz, xy + xz]) yhi = float(y_line[1]) - np.max([0.0, yz]) ylo = float(y_line[0]) - np.min([0.0, yz]) zhi = float(z_line[1]) zlo = float(z_line[0]) lx = xhi - xlo ly = yhi - ylo lz = zhi - zlo c = np.sqrt(lz ** 2 + xz ** 2 + yz ** 2) b = np.sqrt(ly ** 2 + xy ** 2) a = lx alpha = np.arccos((yz * ly + xy * xz) / (b * c)) beta = np.arccos(xz / c) gamma = np.arccos(xy / b) # Box Information lengths = u.unyt_array([a, b, c], get_units(unit_style)["distance"]) angles = u.unyt_array([alpha, beta, gamma], u.radian) angles.to(get_units(unit_style)["angle"]) topology.box = Box(lengths, angles) else: # Box Information lengths = u.unyt_array([x, y, z], get_units(unit_style)["distance"]) topology.box = Box(lengths) return topology def _get_ff_information(filename, unit_style, topology): """Parse atom-type information.""" with open(filename, "r") as lammps_file: types = False for i, line in enumerate(lammps_file): if "atom" in line: n_atomtypes = int(line.split()[0]) types = True elif "Masses" in line: break if types == False: return topology mass_lines = open(filename, "r").readlines()[i + 2 : i + n_atomtypes + 2] type_list = list() for line in mass_lines: atom_type = AtomType( name=line.split()[0], mass=float(line.split()[1]) * get_units(unit_style)["mass"], ) type_list.append(atom_type) with open(filename, "r") as lammps_file: for i, line in enumerate(lammps_file): if "Pair" in line: break # Need to figure out if we're going have mixing rules printed out # Currently only reading in LJ params pair_lines = open(filename, "r").readlines()[i + 2 : i + n_atomtypes + 2] for i, pair in enumerate(pair_lines): if len(pair.split()) == 3: type_list[i].parameters["sigma"] = ( float(pair.split()[2]) * get_units(unit_style)["distance"] ) type_list[i].parameters["epsilon"] = ( float(pair.split()[1]) * get_units(unit_style)["energy"] ) elif len(pair.split()) == 4: warnings.warn("Currently not reading in mixing rules") return topology, type_list