from __future__ import print_function
# lammps.py (2011/03/29)
# An ASE calculator for the LAMMPS classical MD code available from
# http://lammps.sandia.gov/
# The environment variable LAMMPS_COMMAND must be defined to point to the
# LAMMPS binary.
#
# Copyright (C) 2009 - 2011 Joerg Meyer, joerg.meyer@ch.tum.de
#
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
# License as published by the Free Software Foundation; either
# version 2.1 of the License, or (at your option) any later version.
#
# This library 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
# Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public
# License along with this file; if not, write to the Free Software
# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301
# USA or see <http://www.gnu.org/licenses/>.
import os
import shutil
import shlex
from subprocess import Popen, PIPE
from threading import Thread
from re import compile as re_compile, IGNORECASE
from tempfile import mkdtemp, NamedTemporaryFile, mktemp as uns_mktemp
import numpy as np
import decimal as dec
from ase import Atoms
from ase.parallel import paropen
from ase.units import GPa, Ang, fs
from ase.utils import basestring
__all__ = ['LAMMPS', 'write_lammps_data']
# "End mark" used to indicate that the calculation is done
CALCULATION_END_MARK = '__end_of_ase_invoked_calculation__'
[docs]class LAMMPS:
def __init__(self, label='lammps', tmp_dir=None, parameters={},
specorder=None, files=[], always_triclinic=False,
keep_alive=True, keep_tmp_files=False,
no_data_file=False):
"""The LAMMPS calculators object
files: list
Short explanation XXX
parameters: dict
Short explanation XXX
specorder: list
Short explanation XXX
keep_tmp_files: bool
Retain any temporary files created. Mostly useful for debugging.
tmp_dir: str
path/dirname (default None -> create automatically).
Explicitly control where the calculator object should create
its files. Using this option implies 'keep_tmp_files'
no_data_file: bool
Controls whether an explicit data file will be used for feeding
atom coordinates into lammps. Enable it to lessen the pressure on
the (tmp) file system. THIS OPTION MIGHT BE UNRELIABLE FOR CERTAIN
CORNER CASES (however, if it fails, you will notice...).
keep_alive: bool
When using LAMMPS as a spawned subprocess, keep the subprocess
alive (but idling when unused) along with the calculator object.
always_triclinic: bool
Force use of a triclinic cell in LAMMPS, even if the cell is
a perfect parallelepiped.
"""
self.label = label
self.parameters = parameters
self.specorder = specorder
self.files = files
self.always_triclinic = always_triclinic
self.calls = 0
self.forces = None
self.keep_alive = keep_alive
self.keep_tmp_files = keep_tmp_files
self.no_data_file = no_data_file
# if True writes velocities from atoms.get_velocities() to LAMMPS input
self.write_velocities = False
# file object, if is not None the trajectory will be saved in it
self.trajectory_out = None
# period of system snapshot saving (in MD steps)
self.dump_period = 1
if tmp_dir is not None:
# If tmp_dir is pointing somewhere, don't remove stuff!
self.keep_tmp_files = True
self._lmp_handle = None # To handle the lmp process
# read_log depends on that the first (three) thermo_style custom args
# can be capitilized and matched against the log output. I.e.
# don't use e.g. 'ke' or 'cpu' which are labeled KinEng and CPU.
self._custom_thermo_args = ['step', 'temp', 'press', 'cpu',
'pxx', 'pyy', 'pzz', 'pxy', 'pxz', 'pyz',
'ke', 'pe', 'etotal',
'vol', 'lx', 'ly', 'lz', 'atoms']
self._custom_thermo_mark = ' '.join([x.capitalize() for x in
self._custom_thermo_args[0:3]])
# Match something which can be converted to a float
f_re = r'([+-]?(?:(?:\d+(?:\.\d*)?|\.\d+)(?:e[+-]?\d+)?|nan|inf))'
n = len(self._custom_thermo_args)
# Create a re matching exactly N white space separated floatish things
self._custom_thermo_re = re_compile(
r'^\s*' + r'\s+'.join([f_re] * n) + r'\s*$', flags=IGNORECASE)
# thermo_content contains data "written by" thermo_style.
# It is a list of dictionaries, each dict (one for each line
# printed by thermo_style) contains a mapping between each
# custom_thermo_args-argument and the corresponding
# value as printed by lammps. thermo_content will be
# re-populated by the read_log method.
self.thermo_content = []
if tmp_dir is None:
self.tmp_dir = mkdtemp(prefix='LAMMPS-')
else:
self.tmp_dir = os.path.realpath(tmp_dir)
if not os.path.isdir(self.tmp_dir):
os.mkdir(self.tmp_dir, 0o755)
for f in files:
shutil.copy(f, os.path.join(self.tmp_dir, os.path.basename(f)))
def clean(self, force=False):
self._lmp_end()
if not self.keep_tmp_files:
shutil.rmtree(self.tmp_dir)
def get_potential_energy(self, atoms):
self.update(atoms)
return self.thermo_content[-1]['pe']
def get_forces(self, atoms):
self.update(atoms)
return self.forces.copy()
def get_stress(self, atoms):
self.update(atoms)
tc = self.thermo_content[-1]
# 1 bar (used by lammps for metal units) = 1e-4 GPa
return np.array([tc[i] for i in ('pxx', 'pyy', 'pzz', 'pyz', 'pxz',
'pxy')]) * (-1e-4 * GPa)
def update(self, atoms):
if not hasattr(self, 'atoms') or self.atoms != atoms:
self.calculate(atoms)
def calculate(self, atoms):
self.atoms = atoms.copy()
pbc = self.atoms.get_pbc()
if all(pbc):
cell = self.atoms.get_cell()
elif not any(pbc):
# large enough cell for non-periodic calculation -
# LAMMPS shrink-wraps automatically via input command
# "periodic s s s"
# below
cell = 2 * np.max(np.abs(self.atoms.get_positions())) * np.eye(3)
else:
print("WARNING: semi-periodic ASE cell detected - translation")
print(" to proper LAMMPS input cell might fail")
cell = self.atoms.get_cell()
self.prism = Prism(cell)
self.run()
def _lmp_alive(self):
# Return True if this calculator is currently handling a running
# lammps process
return self._lmp_handle and not isinstance(
self._lmp_handle.poll(), int)
def _lmp_end(self):
# Close lammps input and wait for lammps to end. Return process
# return value
if self._lmp_alive():
self._lmp_handle.stdin.close()
return self._lmp_handle.wait()
def run(self, set_atoms=False):
"""Method which explicitly runs LAMMPS."""
self.calls += 1
# set LAMMPS command from environment variable
if 'LAMMPS_COMMAND' in os.environ:
lammps_cmd_line = shlex.split(os.environ['LAMMPS_COMMAND'])
if len(lammps_cmd_line) == 0:
self.clean()
raise RuntimeError('The LAMMPS_COMMAND environment variable '
'must not be empty')
# want always an absolute path to LAMMPS binary when calling from
# self.dir
lammps_cmd_line[0] = os.path.abspath(lammps_cmd_line[0])
else:
self.clean()
raise RuntimeError(
'Please set LAMMPS_COMMAND environment variable')
if 'LAMMPS_OPTIONS' in os.environ:
lammps_options = shlex.split(os.environ['LAMMPS_OPTIONS'])
else:
lammps_options = shlex.split('-echo log -screen none')
# change into subdirectory for LAMMPS calculations
cwd = os.getcwd()
os.chdir(self.tmp_dir)
# setup file names for LAMMPS calculation
label = '{0}{1:>06}'.format(self.label, self.calls)
lammps_in = uns_mktemp(prefix='in_' + label, dir=self.tmp_dir)
lammps_log = uns_mktemp(prefix='log_' + label, dir=self.tmp_dir)
lammps_trj_fd = NamedTemporaryFile(
prefix='trj_' + label, dir=self.tmp_dir,
delete=(not self.keep_tmp_files))
lammps_trj = lammps_trj_fd.name
if self.no_data_file:
lammps_data = None
else:
lammps_data_fd = NamedTemporaryFile(
prefix='data_' + label, dir=self.tmp_dir,
delete=(not self.keep_tmp_files))
self.write_lammps_data(lammps_data=lammps_data_fd)
lammps_data = lammps_data_fd.name
lammps_data_fd.flush()
# see to it that LAMMPS is started
if not self._lmp_alive():
# Attempt to (re)start lammps
self._lmp_handle = Popen(
lammps_cmd_line + lammps_options + ['-log', '/dev/stdout'],
stdin=PIPE, stdout=PIPE)
lmp_handle = self._lmp_handle
# Create thread reading lammps stdout (for reference, if requested,
# also create lammps_log, although it is never used)
if self.keep_tmp_files:
lammps_log_fd = open(lammps_log, 'wb')
fd = SpecialTee(lmp_handle.stdout, lammps_log_fd)
else:
fd = lmp_handle.stdout
thr_read_log = Thread(target=self.read_lammps_log, args=(fd,))
thr_read_log.start()
# write LAMMPS input (for reference, also create the file lammps_in,
# although it is never used)
if self.keep_tmp_files:
lammps_in_fd = open(lammps_in, 'wb')
fd = SpecialTee(lmp_handle.stdin, lammps_in_fd)
else:
fd = lmp_handle.stdin
self.write_lammps_in(lammps_in=fd, lammps_trj=lammps_trj,
lammps_data=lammps_data)
if self.keep_tmp_files:
lammps_in_fd.close()
# Wait for log output to be read (i.e., for LAMMPS to finish)
# and close the log file if there is one
thr_read_log.join()
if self.keep_tmp_files:
lammps_log_fd.close()
if not self.keep_alive:
self._lmp_end()
exitcode = lmp_handle.poll()
if exitcode and exitcode != 0:
cwd = os.getcwd()
raise RuntimeError('LAMMPS exited in {} with exit code: {}.'
''.format(cwd, exitcode))
# A few sanity checks
if len(self.thermo_content) == 0:
raise RuntimeError('Failed to retrieve any thermo_style-output')
if int(self.thermo_content[-1]['atoms']) != len(self.atoms):
# This obviously shouldn't happen, but if prism.fold_...() fails,
# it could
raise RuntimeError('Atoms have gone missing')
self.read_lammps_trj(lammps_trj=lammps_trj, set_atoms=set_atoms)
lammps_trj_fd.close()
if not self.no_data_file:
lammps_data_fd.close()
os.chdir(cwd)
def write_lammps_data(self, lammps_data=None):
"""Method which writes a LAMMPS data file with atomic structure."""
if lammps_data is None:
lammps_data = 'data.' + self.label
write_lammps_data(
lammps_data, self.atoms, self.specorder,
force_skew=self.always_triclinic, prismobj=self.prism,
velocities=self.write_velocities)
def write_lammps_in(self, lammps_in=None, lammps_trj=None,
lammps_data=None):
"""Write a LAMMPS in_ file with run parameters and settings."""
if isinstance(lammps_in, basestring):
f = paropen(lammps_in, 'wb')
close_in_file = True
else:
# Expect lammps_in to be a file-like object
f = lammps_in
close_in_file = False
if self.keep_tmp_files:
f.write('# (written by ASE)\n'.encode('utf-8'))
# Write variables
f.write(('clear\n'
'variable dump_file string "{0}"\n'
'variable data_file string "{1}"\n'
).format(lammps_trj, lammps_data).encode('utf-8'))
parameters = self.parameters
if 'package' in parameters:
f.write(('\n'.join(['package {0}'.format(p)
for p in parameters['package']]) +
'\n').encode('utf-8'))
pbc = self.atoms.get_pbc()
f.write('units metal \n'.encode('utf-8'))
if 'boundary' in parameters:
f.write('boundary {0} \n'.format(
parameters['boundary']).encode('utf-8'))
else:
f.write('boundary {0} {1} {2} \n'.format(
*tuple('sp'[int(x)] for x in pbc)).encode('utf-8'))
f.write('atom_modify sort 0 0.0 \n'.encode('utf-8'))
for key in ('neighbor', 'newton'):
if key in parameters:
f.write('{0} {1} \n'.format(
key, parameters[key]).encode('utf-8'))
f.write('\n'.encode('utf-8'))
# If self.no_lammps_data,
# write the simulation box and the atoms
if self.no_data_file:
if self.keep_tmp_files:
f.write('## Original ase cell\n'.encode('utf-8'))
f.write(''.join(['# {0:.16} {1:.16} {2:.16}\n'.format(*x)
for x in self.atoms.get_cell()]
).encode('utf-8'))
f.write('lattice sc 1.0\n'.encode('utf-8'))
xhi, yhi, zhi, xy, xz, yz = self.prism.get_lammps_prism_str()
if self.always_triclinic or self.prism.is_skewed():
f.write('region asecell prism 0.0 {0} 0.0 {1} 0.0 {2} '
''.format(xhi, yhi, zhi).encode('utf-8'))
f.write('{0} {1} {2} side in units box\n'
''.format(xy, xz, yz).encode('utf-8'))
else:
f.write(('region asecell block 0.0 {0} 0.0 {1} 0.0 {2} '
'side in units box\n').format(
xhi, yhi, zhi).encode('utf-8'))
symbols = self.atoms.get_chemical_symbols()
if self.specorder is None:
# By default, atom types in alphabetic order
species = sorted(set(symbols))
else:
# By request, specific atom type ordering
species = self.specorder
n_atom_types = len(species)
species_i = dict([(s, i + 1) for i, s in enumerate(species)])
f.write('create_box {0} asecell\n'.format(
n_atom_types).encode('utf-8'))
for s, pos in zip(symbols, self.atoms.get_positions()):
if self.keep_tmp_files:
f.write('# atom pos in ase cell: {0:.16} {1:.16} {2:.16}'
'\n'.format(*tuple(pos)).encode('utf-8'))
f.write('create_atoms {0} single {1} {2} {3} units box\n'
''.format(
*((species_i[s],) +
self.prism.pos_to_lammps_fold_str(pos))
).encode('utf-8'))
# if NOT self.no_lammps_data, then simply refer to the data-file
else:
f.write('read_data {0}\n'.format(lammps_data).encode('utf-8'))
# Write interaction stuff
f.write('\n### interactions \n'.encode('utf-8'))
if ('pair_style' in parameters) and ('pair_coeff' in parameters):
pair_style = parameters['pair_style']
f.write('pair_style {0} \n'.format(pair_style).encode('utf-8'))
for pair_coeff in parameters['pair_coeff']:
f.write('pair_coeff {0} \n'
''.format(pair_coeff).encode('utf-8'))
if 'mass' in parameters:
for mass in parameters['mass']:
f.write('mass {0} \n'.format(mass).encode('utf-8'))
else:
# simple default parameters
# that should always make the LAMMPS calculation run
f.write('pair_style lj/cut 2.5 \n'
'pair_coeff * * 1 1 \n'
'mass * 1.0 \n'.encode('utf-8'))
if 'group' in parameters:
f.write(('\n'.join(['group {0}'.format(p)
for p in parameters['group']]) +
'\n').encode('utf-8'))
f.write(
'\n### run\n'
'fix fix_nve all nve\n'.encode('utf-8'))
if 'fix' in parameters:
f.write(('\n'.join(['fix {0}'.format(p)
for p in parameters['fix']]) +
'\n').encode('utf-8'))
f.write(
'dump dump_all all custom {1} {0} id type x y z vx vy vz '
'fx fy fz\n'
''.format(lammps_trj, self.dump_period).encode('utf-8'))
f.write('thermo_style custom {0}\n'
'thermo_modify flush yes\n'
'thermo 1\n'.format(
' '.join(self._custom_thermo_args)).encode('utf-8'))
if 'timestep' in parameters:
f.write('timestep {0}\n'.format(
parameters['timestep']).encode('utf-8'))
if 'minimize' in parameters:
f.write('minimize {0}\n'.format(
parameters['minimize']).encode('utf-8'))
if 'run' in parameters:
f.write('run {0}\n'.format(parameters['run']).encode('utf-8'))
if not (('minimize' in parameters) or ('run' in parameters)):
f.write('run 0\n'.encode('utf-8'))
f.write('print "{0}" \n'.format(CALCULATION_END_MARK).encode('utf-8'))
# Force LAMMPS to flush log
f.write('log /dev/stdout\n'.encode('utf-8'))
f.flush()
if close_in_file:
f.close()
def read_lammps_log(self, lammps_log=None, PotEng_first=False):
"""Method which reads a LAMMPS output log file."""
if lammps_log is None:
lammps_log = self.label + '.log'
if isinstance(lammps_log, basestring):
f = paropen(lammps_log, 'wb')
close_log_file = True
else:
# Expect lammps_in to be a file-like object
f = lammps_log
close_log_file = False
thermo_content = []
line = f.readline().decode('utf-8')
while line and line.strip() != CALCULATION_END_MARK:
# get thermo output
if line.startswith(self._custom_thermo_mark):
m = True
while m:
line = f.readline().decode('utf-8')
m = self._custom_thermo_re.match(line)
if m:
# create a dictionary between each of the
# thermo_style args and it's corresponding value
thermo_content.append(
dict(zip(self._custom_thermo_args,
map(float, m.groups()))))
else:
line = f.readline().decode('utf-8')
if close_log_file:
f.close()
self.thermo_content = thermo_content
def read_lammps_trj(self, lammps_trj=None, set_atoms=False):
"""Method which reads a LAMMPS dump file."""
if lammps_trj is None:
lammps_trj = self.label + '.lammpstrj'
f = paropen(lammps_trj, 'r')
while True:
line = f.readline()
if not line:
break
# TODO: extend to proper dealing with multiple steps in one
# trajectory file
if 'ITEM: TIMESTEP' in line:
n_atoms = 0
lo = []
hi = []
tilt = []
id = []
type = []
positions = []
velocities = []
forces = []
if 'ITEM: NUMBER OF ATOMS' in line:
line = f.readline()
n_atoms = int(line.split()[0])
if 'ITEM: BOX BOUNDS' in line:
# save labels behind "ITEM: BOX BOUNDS" in triclinic case
# (>=lammps-7Jul09)
tilt_items = line.split()[3:]
for i in range(3):
line = f.readline()
fields = line.split()
lo.append(float(fields[0]))
hi.append(float(fields[1]))
if len(fields) >= 3:
tilt.append(float(fields[2]))
if 'ITEM: ATOMS' in line:
# (reliably) identify values by labels behind "ITEM: ATOMS"
# - requires >=lammps-7Jul09
# create corresponding index dictionary before iterating over
# atoms to (hopefully) speed up lookups...
atom_attributes = {}
for (i, x) in enumerate(line.split()[2:]):
atom_attributes[x] = i
for n in range(n_atoms):
line = f.readline()
fields = line.split()
id.append(int(fields[atom_attributes['id']]))
type.append(int(fields[atom_attributes['type']]))
positions.append([float(fields[atom_attributes[x]])
for x in ['x', 'y', 'z']])
velocities.append([float(fields[atom_attributes[x]])
for x in ['vx', 'vy', 'vz']])
forces.append([float(fields[atom_attributes[x]])
for x in ['fx', 'fy', 'fz']])
# Re-order items according to their 'id' since running in
# parallel can give arbitrary ordering.
type = [x for _,x in sorted(zip(id, type))]
positions = [x for _,x in sorted(zip(id, positions))]
velocities = [x for _,x in sorted(zip(id, velocities))]
forces = [x for _,x in sorted(zip(id, forces))]
# determine cell tilt (triclinic case!)
if len(tilt) >= 3:
# for >=lammps-7Jul09 use labels behind "ITEM: BOX BOUNDS"
# to assign tilt (vector) elements ...
if len(tilt_items) >= 3:
xy = tilt[tilt_items.index('xy')]
xz = tilt[tilt_items.index('xz')]
yz = tilt[tilt_items.index('yz')]
# ... otherwise assume default order in 3rd column
# (if the latter was present)
else:
xy = tilt[0]
xz = tilt[1]
yz = tilt[2]
else:
xy = xz = yz = 0
xhilo = (hi[0] - lo[0]) - xy - xz
yhilo = (hi[1] - lo[1]) - yz
zhilo = (hi[2] - lo[2])
# The simulation box bounds are included in each snapshot and
# if the box is triclinic (non-orthogonal), then the tilt
# factors are also printed; see the region prism command for
# a description of tilt factors.
# For triclinic boxes the box bounds themselves (first 2
# quantities on each line) are a true "bounding box" around
# the simulation domain, which means they include the effect of
# any tilt.
# [ http://lammps.sandia.gov/doc/dump.html , lammps-7Jul09 ]
#
# This *should* extract the lattice vectors that LAMMPS uses
# from the true "bounding box" printed in the dump file.
# It might fail in some cases (negative tilts?!) due to the
# MIN / MAX construction of these box corners:
#
# void Domain::set_global_box()
# [...]
# if (triclinic) {
# [...]
# boxlo_bound[0] = MIN(boxlo[0],boxlo[0]+xy);
# boxlo_bound[0] = MIN(boxlo_bound[0],boxlo_bound[0]+xz);
# boxlo_bound[1] = MIN(boxlo[1],boxlo[1]+yz);
# boxlo_bound[2] = boxlo[2];
# # boxhi_bound[0] = MAX(boxhi[0],boxhi[0]+xy);
# boxhi_bound[0] = MAX(boxhi_bound[0],boxhi_bound[0]+xz);
# boxhi_bound[1] = MAX(boxhi[1],boxhi[1]+yz);
# boxhi_bound[2] = boxhi[2];
# }
# [ lammps-7Jul09/src/domain.cpp ]
#
cell = [[xhilo, 0, 0], [xy, yhilo, 0], [xz, yz, zhilo]]
# These have been put into the correct order
cell_atoms = np.array(cell)
type_atoms = np.array(type)
if self.atoms:
cell_atoms = self.atoms.get_cell()
# BEWARE: reconstructing the rotation from the LAMMPS
# output trajectory file fails in case of shrink
# wrapping for a non-periodic direction
# -> hence rather obtain rotation from prism object
# used to generate the LAMMPS input
# rotation_lammps2ase = np.dot(
# np.linalg.inv(np.array(cell)), cell_atoms)
rotation_lammps2ase = np.linalg.inv(self.prism.R)
type_atoms = self.atoms.get_atomic_numbers()
positions_atoms = np.dot(positions, rotation_lammps2ase)
velocities_atoms = np.dot(velocities, rotation_lammps2ase)
forces_atoms = np.dot(forces, rotation_lammps2ase)
if set_atoms:
# assume periodic boundary conditions here (as in
# write_lammps)
self.atoms = Atoms(type_atoms, positions=positions_atoms,
cell=cell_atoms)
self.atoms.set_velocities(velocities_atoms
* (Ang/(fs*1000.)))
self.forces = forces_atoms
if self.trajectory_out is not None:
tmp_atoms = Atoms(type_atoms, positions=positions_atoms,
cell=cell_atoms)
tmp_atoms.set_velocities(velocities_atoms)
self.trajectory_out.write(tmp_atoms)
f.close()
class SpecialTee(object):
"""A special purpose, with limited applicability, tee-like thing.
A subset of stuff read from, or written to, orig_fd,
is also written to out_fd.
It is used by the lammps calculator for creating file-logs of stuff
read from, or written to, stdin and stdout, respectively.
"""
def __init__(self, orig_fd, out_fd):
self._orig_fd = orig_fd
self._out_fd = out_fd
self.name = orig_fd.name
def write(self, data):
self._orig_fd.write(data)
self._out_fd.write(data)
self.flush()
def read(self, *args, **kwargs):
data = self._orig_fd.read(*args, **kwargs)
self._out_fd.write(data)
return data
def readline(self, *args, **kwargs):
data = self._orig_fd.readline(*args, **kwargs)
self._out_fd.write(data)
return data
def readlines(self, *args, **kwargs):
data = self._orig_fd.readlines(*args, **kwargs)
self._out_fd.write(''.join(data))
return data
def flush(self):
self._orig_fd.flush()
self._out_fd.flush()
class Prism(object):
def __init__(self, cell, pbc=(True, True, True), digits=10):
"""Create a lammps-style triclinic prism object from a cell
The main purpose of the prism-object is to create suitable
string representations of prism limits and atom positions
within the prism.
When creating the object, the digits parameter (default set to 10)
specify the precision to use.
lammps is picky about stuff being within semi-open intervals,
e.g. for atom positions (when using create_atom in the in-file),
x must be within [xlo, xhi).
"""
a, b, c = cell
an, bn, cn = [np.linalg.norm(v) for v in cell]
alpha = np.arccos(np.dot(b, c) / (bn * cn))
beta = np.arccos(np.dot(a, c) / (an * cn))
gamma = np.arccos(np.dot(a, b) / (an * bn))
xhi = an
xyp = np.cos(gamma) * bn
yhi = np.sin(gamma) * bn
xzp = np.cos(beta) * cn
yzp = (bn * cn * np.cos(alpha) - xyp * xzp) / yhi
zhi = np.sqrt(cn**2 - xzp**2 - yzp**2)
# Set precision
self.car_prec = dec.Decimal('10.0') ** \
int(np.floor(np.log10(max((xhi, yhi, zhi)))) - digits)
self.dir_prec = dec.Decimal('10.0') ** (-digits)
self.acc = float(self.car_prec)
self.eps = np.finfo(xhi).eps
# For rotating positions from ase to lammps
Apre = np.array(((xhi, 0, 0),
(xyp, yhi, 0),
(xzp, yzp, zhi)))
self.R = np.dot(np.linalg.inv(cell), Apre)
# Actual lammps cell may be different from what is used to create R
def fold(vec, pvec, i):
p = pvec[i]
x = vec[i] + 0.5 * p
n = (np.mod(x, p) - x) / p
return [float(self.f2qdec(a)) for a in (vec + n * pvec)]
Apre[1, :] = fold(Apre[1, :], Apre[0, :], 0)
Apre[2, :] = fold(Apre[2, :], Apre[1, :], 1)
Apre[2, :] = fold(Apre[2, :], Apre[0, :], 0)
self.A = Apre
self.Ainv = np.linalg.inv(self.A)
if self.is_skewed() and \
(not (pbc[0] and pbc[1] and pbc[2])):
raise RuntimeError('Skewed lammps cells MUST have '
'PBC == True in all directions!')
def f2qdec(self, f):
return dec.Decimal(repr(f)).quantize(self.car_prec, dec.ROUND_DOWN)
def f2qs(self, f):
return str(self.f2qdec(f))
def f2s(self, f):
return str(dec.Decimal(repr(f)).quantize(self.car_prec,
dec.ROUND_HALF_EVEN))
def dir2car(self, v):
"""Direct to cartesian coordinates"""
return np.dot(v, self.A)
def car2dir(self, v):
"""Cartesian to direct coordinates"""
return np.dot(v, self.Ainv)
def fold_to_str(self, v):
"""Fold a position into the lammps cell (semi open)
Returns tuple of str.
"""
# Two-stage fold, first into box, then into semi-open interval
# (within the given precision).
d = [x % (1 - self.dir_prec) for x in
map(dec.Decimal,
map(repr, np.mod(self.car2dir(v) + self.eps, 1.0)))]
return tuple([self.f2qs(x) for x in
self.dir2car(list(map(float, d)))])
def get_lammps_prism(self):
A = self.A
return A[0, 0], A[1, 1], A[2, 2], A[1, 0], A[2, 0], A[2, 1]
def get_lammps_prism_str(self):
"""Return a tuple of strings"""
p = self.get_lammps_prism()
return tuple([self.f2s(x) for x in p])
def positions_to_lammps_strs(self, positions):
"""Rotate an ase-cell position to the lammps cell orientation
Returns tuple of str.
"""
rot_positions = np.dot(positions, self.R)
return [tuple([self.f2s(x) for x in position])
for position in rot_positions]
def pos_to_lammps_fold_str(self, position):
"""Rotate and fold an ase-cell position into the lammps cell
Returns tuple of str.
"""
return self.fold_to_str(np.dot(position, self.R))
def is_skewed(self):
acc = self.acc
prism = self.get_lammps_prism()
axy, axz, ayz = [np.abs(x) for x in prism[3:]]
return (axy >= acc) or (axz >= acc) or (ayz >= acc)
def write_lammps_data(fileobj, atoms, specorder=None, force_skew=False,
prismobj=None, velocities=False):
"""Write atomic structure data to a LAMMPS data_ file."""
if isinstance(fileobj, basestring):
f = paropen(fileobj, 'wb')
close_file = True
else:
# Presume fileobj acts like a fileobj
f = fileobj
close_file = False
if isinstance(atoms, list):
if len(atoms) > 1:
raise ValueError(
'Can only write one configuration to a lammps data file!')
atoms = atoms[0]
f.write('{0} (written by ASE) \n\n'.format(f.name).encode('utf-8'))
symbols = atoms.get_chemical_symbols()
n_atoms = len(symbols)
f.write('{0} \t atoms \n'.format(n_atoms).encode('utf-8'))
if specorder is None:
# This way it is assured that LAMMPS atom types are always
# assigned predictably according to the alphabetic order
species = sorted(set(symbols))
else:
# To index elements in the LAMMPS data file
# (indices must correspond to order in the potential file)
species = specorder
n_atom_types = len(species)
f.write('{0} atom types\n'.format(n_atom_types).encode('utf-8'))
if prismobj is None:
p = Prism(atoms.get_cell())
else:
p = prismobj
xhi, yhi, zhi, xy, xz, yz = p.get_lammps_prism_str()
f.write('0.0 {0} xlo xhi\n'.format(xhi).encode('utf-8'))
f.write('0.0 {0} ylo yhi\n'.format(yhi).encode('utf-8'))
f.write('0.0 {0} zlo zhi\n'.format(zhi).encode('utf-8'))
if force_skew or p.is_skewed():
f.write('{0} {1} {2} xy xz yz\n'.format(xy, xz, yz).encode('utf-8'))
f.write('\n\n'.encode('utf-8'))
f.write('Atoms \n\n'.encode('utf-8'))
for i, r in enumerate(p.positions_to_lammps_strs(atoms.get_positions())):
s = species.index(symbols[i]) + 1
f.write('{0:>6} {1:>3} {2} {3} {4}\n'.format(
*(i + 1, s) + tuple(r)).encode('utf-8'))
if velocities and atoms.get_velocities() is not None:
f.write('\n\nVelocities \n\n'.encode('utf-8'))
for i, v in enumerate(atoms.get_velocities() / (Ang/(fs*1000.))):
f.write('{0:>6} {1} {2} {3}\n'.format(
*(i + 1,) + tuple(v)).encode('utf-8'))
f.flush()
if close_file:
f.close()
if __name__ == '__main__':
pair_style = 'eam'
Pd_eam_file = 'Pd_u3.eam'
pair_coeff = ['* * ' + Pd_eam_file]
parameters = {'pair_style': pair_style, 'pair_coeff': pair_coeff}
files = [Pd_eam_file]
calc = LAMMPS(parameters=parameters, files=files)
a0 = 3.93
b0 = a0 / 2.0
if True:
bulk = Atoms(
['Pd'] * 4,
positions=[(0, 0, 0), (b0, b0, 0), (b0, 0, b0), (0, b0, b0)],
cell=[a0] * 3, pbc=True)
# test get_forces
print('forces for a = {0}'.format(a0))
print(calc.get_forces(bulk))
# single points for various lattice constants
bulk.set_calculator(calc)
for n in range(-5, 5, 1):
a = a0 * (1 + n / 100.0)
bulk.set_cell([a] * 3)
print('a : {0} , total energy : {1}'.format(
a, bulk.get_potential_energy()))
calc.clean()