Source code for gmso.formats.lammpsdata

from __future__ import division

import warnings
import numpy as np
import unyt as u
import datetime

from unyt.array import allclose_units
from gmso.core.site import Site
from gmso.core.atom_type import AtomType
from gmso.core.bond_type import BondType
from gmso.core.angle_type import AngleType
from gmso.core.bond import Bond
from gmso.core.angle import Angle
from gmso.core.topology import Topology
from gmso.core.box import Box
from gmso.core.element import element_by_mass


[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') for idx, dihedral_type in enumerate(topology.dihedral_types): if dihedral_type.name == 'RyckaertBellemansTorsionPotential': 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']/2, dihedral_type.parameters['k2']/2, dihedral_type.parameters['k3']/2, dihedral_type.parameters['k4']/2 )) # 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): """General function to 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, connection_type=connection_type_list[int(line.split()[1])-1], ) elif connection_type == 'angle': connection = Angle( connection_members=site_list, connection_type=connection_type_list[int(line.split()[1])-1], ) topology.add_connection(connection) return topology def _get_atoms(filename, topology, unit_style, type_list): """Function to 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 = Site( 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): """Function to 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): """Function to 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