diff --git a/examples/Au18/structopt.in.json b/examples/Au18/structopt.in.json index 606cd20..0d587e2 100644 --- a/examples/Au18/structopt.in.json +++ b/examples/Au18/structopt.in.json @@ -57,10 +57,11 @@ }, "mutations": { "options": [ - "lattice_alteration", - "lattice_alteration_group", - "rotation_geo", - "rotation" + "move_atoms", + "swap_positions", + "swap_species", + "rotate_atoms", + "rotate_cluster" ] }, "fingerprinters": { diff --git a/examples/Au18/structopt_load_xyz.in.json b/examples/Au18/structopt_load_xyz.in.json new file mode 100644 index 0000000..5391fb6 --- /dev/null +++ b/examples/Au18/structopt_load_xyz.in.json @@ -0,0 +1,96 @@ +{ + "globals": + { + "USE_MPI4PY": true + }, + "generators": { + "structure_type": "cluster", + "initializers": [ + { + "number_of_individuals": 5, + "data": { + "filenames": [ + "/home/jmaldonis/StructOpt_modular/examples/Au18/structure1.xyz", + "/home/jmaldonis/StructOpt_modular/examples/Au18/structure2.xyz", + "/home/jmaldonis/StructOpt_modular/examples/Au18/structure3.xyz", + "/home/jmaldonis/StructOpt_modular/examples/Au18/structure4.xyz", + "/home/jmaldonis/StructOpt_modular/examples/Au18/structure5.xyz" + ] + } + } + ] + }, + "fitnesses": + { + "modules": [ + "LAMMPS", + "FEMSIM" + ], + "weights": [ + 1.0, + 1.0 + ], + "FEMSIM": { + "parameter_filename": "/home/maldonis/USIT/StructOpt/StructOpt/fitness/FEMSIM/femsim-hrmc/parameters/femsim.in", + "vk_data_filename": "/home/maldonis/USIT/StructOpt/StructOpt/fitness/FEMSIM/femsim-hrmc/data/fem_exp_Zr50Cu35Al15_t3.txt", + "xsize": 28.2842712474619, + "ysize": 28.2842712474619, + "zsize": 28.2842712474619, + "Q": 0.0305, + "nphi": 1, + "npsi": 40, + "ntheta": 20, + "thickness_scaling_factor": 28.6378 + }, + "LAMMPS":{ + "keep_files": true, + "min_style": "cg\nmin_modify line quadratic", + "minimize": "1e-8 1e-8 5000 10000", + "pair_style": "eam", + "parallel": true, + "potential_file": "/home/jmaldonis/StructOpt/examples/Au18/Au_u3.eam", + "thermo_steps": 1000 + } + }, + "relaxations": + { + "modules": [ + "LAMMPS" + ], + "LAMMPS":{ + "keep_files": true, + "min_style": "cg\nmin_modify line quadratic", + "minimize": "1e-8 1e-8 5000 10000", + "pair_style": "eam", + "parallel": true, + "potential_file": "/home/jmaldonis/StructOpt/examples/Au18/Au_u3.eam", + "thermo_steps": 1000 + } + }, + "mutations": { + "options": [ + "move_atoms", + "swap_positions", + "swap_species", + "rotate_atoms", + "rotate_cluster" + ] + }, + "fingerprinters": { + "options": [] + }, + + "crossovers": { + "options": ["do_nothing"] + }, + "predators": { + "options": ["do_nothing"] + }, + "selections": { + "options": ["do_nothing"] + }, + + "fileio": {}, + "postprocessing": {}, + "tools": {} +} diff --git a/examples/Au18/structure1.xyz b/examples/Au18/structure1.xyz new file mode 100644 index 0000000..99ea026 --- /dev/null +++ b/examples/Au18/structure1.xyz @@ -0,0 +1,6 @@ +4 +10.0 10.0 10.0 +Zr 1.0 1.0 1.0 +Zr 2.0 2.0 2.0 +Zr 3.0 3.0 3.0 +Zr 4.0 4.0 4.0 diff --git a/examples/Au18/structure2.xyz b/examples/Au18/structure2.xyz new file mode 100644 index 0000000..99ea026 --- /dev/null +++ b/examples/Au18/structure2.xyz @@ -0,0 +1,6 @@ +4 +10.0 10.0 10.0 +Zr 1.0 1.0 1.0 +Zr 2.0 2.0 2.0 +Zr 3.0 3.0 3.0 +Zr 4.0 4.0 4.0 diff --git a/examples/Au18/structure3.xyz b/examples/Au18/structure3.xyz new file mode 100644 index 0000000..99ea026 --- /dev/null +++ b/examples/Au18/structure3.xyz @@ -0,0 +1,6 @@ +4 +10.0 10.0 10.0 +Zr 1.0 1.0 1.0 +Zr 2.0 2.0 2.0 +Zr 3.0 3.0 3.0 +Zr 4.0 4.0 4.0 diff --git a/examples/Au18/structure4.xyz b/examples/Au18/structure4.xyz new file mode 100644 index 0000000..99ea026 --- /dev/null +++ b/examples/Au18/structure4.xyz @@ -0,0 +1,6 @@ +4 +10.0 10.0 10.0 +Zr 1.0 1.0 1.0 +Zr 2.0 2.0 2.0 +Zr 3.0 3.0 3.0 +Zr 4.0 4.0 4.0 diff --git a/examples/Au18/structure5.xyz b/examples/Au18/structure5.xyz new file mode 100644 index 0000000..99ea026 --- /dev/null +++ b/examples/Au18/structure5.xyz @@ -0,0 +1,6 @@ +4 +10.0 10.0 10.0 +Zr 1.0 1.0 1.0 +Zr 2.0 2.0 2.0 +Zr 3.0 3.0 3.0 +Zr 4.0 4.0 4.0 diff --git a/structopt/__init__.py b/structopt/__init__.py index 33225c0..436a33b 100644 --- a/structopt/__init__.py +++ b/structopt/__init__.py @@ -1,8 +1,10 @@ import time import logging +import os from . import fileio -from optimizer import Optimizer +from .optimizer import Optimizer +from . import common def setup(parameter_file): parameters = fileio.parameters.read(parameter_file) @@ -40,6 +42,9 @@ def setup(parameter_file): if logging_level <= logging.DEBUG: debug_logger = fileio.logger_utils.initialize_logger(filename='{}-by-rank-debug.log'.format(parameters.globals.loggername), name="by-rank-debug", level=logging_level) + if not os.path.exists(parameters.globals.output_filename): + os.mkdir(parameters.globals.output_filename) + # Set defaults for parameters and globally scope them parameters = fileio.parameters.set_default(parameters) globals()["parameters"] = parameters diff --git a/structopt/clean.sh b/structopt/clean.sh new file mode 100755 index 0000000..6a2a4c0 --- /dev/null +++ b/structopt/clean.sh @@ -0,0 +1,57 @@ +rm cluster/*.pyc +rm cluster/individual/*.pyc +rm cluster/population/*.pyc +rm common/*.pyc +rm common/individual/*.pyc +rm common/individual/fingerprinters/*.pyc +rm common/individual/fitnesses/*.pyc +rm common/individual/generators/*.pyc +rm common/individual/mutations/*.pyc +rm common/individual/relaxations/*.pyc +rm common/population/*.pyc +rm common/population/crossovers/*.pyc +rm common/population/fitnesses/*.pyc +rm common/population/predators/*.pyc +rm common/population/relaxations/*.pyc +rm common/population/selections/*.pyc +rm crystal/*.pyc +rm crystal/individual/*.pyc +rm crystal/population/*.pyc +rm defect/*.pyc +rm defect/individual/*.pyc +rm defect/population/*.pyc +rm fileio/*.pyc +rm postprocessing/*.pyc +rm surface/*.pyc +rm surface/individual/*.pyc +rm surface/population/*.pyc +rm tools/*.pyc + +rm -rf cluster/__pycache__ +rm -rf cluster/individual/pycache__ +rm -rf cluster/population/pycache__ +rm -rf common/pycache__ +rm -rf common/individual/pycache__ +rm -rf common/individual/fingerprinters/pycache__ +rm -rf common/individual/fitnesses/pycache__ +rm -rf common/individual/generators/pycache__ +rm -rf common/individual/mutations/pycache__ +rm -rf common/individual/relaxations/pycache__ +rm -rf common/population/pycache__ +rm -rf common/population/crossovers/pycache__ +rm -rf common/population/fitnesses/pycache__ +rm -rf common/population/predators/pycache__ +rm -rf common/population/relaxations/pycache__ +rm -rf common/population/selections/pycache__ +rm -rf crystal/pycache__ +rm -rf crystal/individual/pycache__ +rm -rf crystal/population/pycache__ +rm -rf defect/pycache__ +rm -rf defect/individual/pycache__ +rm -rf defect/population/pycache__ +rm -rf fileio/pycache__ +rm -rf postprocessing/pycache__ +rm -rf surface/pycache__ +rm -rf surface/individual/pycache__ +rm -rf surface/population/pycache__ +rm -rf tools/pycache__ diff --git a/structopt/cluster/__init__.py b/structopt/cluster/__init__.py index 5edab68..8eb57bd 100644 --- a/structopt/cluster/__init__.py +++ b/structopt/cluster/__init__.py @@ -2,5 +2,7 @@ class Cluster(Individual): - """A stucture containing a cluster of atoms.""" + """ A stucture containing a non-periodic cluster of atoms. Relaxation algorithms such as LAMMPS and VASP require + a periodic structure, so when run though these types of algorithms, the cluster is embedded in a larger box. + """ diff --git a/structopt/common/__init__.py b/structopt/common/__init__.py index e69de29..3a3141d 100644 --- a/structopt/common/__init__.py +++ b/structopt/common/__init__.py @@ -0,0 +1 @@ +from . import individual, population diff --git a/structopt/common/individual/__init__.py b/structopt/common/individual/__init__.py index 1360feb..9283261 100644 --- a/structopt/common/individual/__init__.py +++ b/structopt/common/individual/__init__.py @@ -2,14 +2,19 @@ import random import ase from importlib import import_module +import numpy as np import structopt +from . import relaxations, fitnesses, mutations, fingerprinters, mutations class Individual(ase.Atoms): """An abstract base class for a structure.""" - def __init__(self, **kwargs): + def __init__(self, index=None, **kwargs): + self._kwargs = kwargs # Store the parameters necessary for initializing for making a copy of self + self.index = index + cls_name = self.__class__.__name__.lower() # Load in the appropriate functionality fingerprinters = import_module('structopt.{cls_name}.individual.fingerprinters'.format(cls_name=cls_name)) @@ -23,22 +28,40 @@ def __init__(self, **kwargs): self.relaxations = relaxations.Relaxations() # Initialize the ase.Atoms structure + super().__init__() generators = import_module('structopt.{cls_name}.individual.generators'.format(cls_name=cls_name)) generators.generate(self, **kwargs) + def mutate(self): self.mutations.select_mutation() return self.mutations.mutate(self) + def relax(self): return self.relaxations.relax(self) + def fitness(self): return self.fitnesses.fitness(self) + def fingerprint(self): return self.fingerprinters.fingerprint(self) - def __repr__(self): - return "<{cls} object at {loc}>".format(cls=self.__class__.__name__, loc=hex(id(self))) + + def get_atom_indices_within_distance_of_atom(self, atom_index, distance): + dists = self.get_distances(atom_index, slice(None, None, None)) + return np.where(dists < distance) + + def get_nearest_atom_indices(self, atom_index, count): + dists = self.get_distances(atom_index, slice(None, None, None))[0] + return np.argsort(dists)[:count] + + def copy(self): + new = self.__class__() + new.index = self.index # TODO + new.extend(self) + return new + diff --git a/structopt/common/individual/generators/__init__.py b/structopt/common/individual/generators/__init__.py index a043696..20f11c8 100644 --- a/structopt/common/individual/generators/__init__.py +++ b/structopt/common/individual/generators/__init__.py @@ -1,15 +1,41 @@ import functools +import numpy as np import structopt from structopt.common.individual import Individual +from .read_xyz import read_xyz - -def generate(individual, *args, **kwargs): +def generate(individual, **kwargs): """ Uses the relevant parameters from structopt to intialize the input Individual by modifying it in-place. Args: individual (Individual): an Individual that is uninitialized - *args, **kwargs: arguments for either ase.Atoms or a different generator function + **kwargs: keyword arguments for either ase.Atoms or a different generator function """ + if 'filenames' in kwargs: + filename = kwargs['filenames'][individual.index] + atoms = read_xyz(filename) + individual.extend(atoms) + #individual.set_atomic_numbers(atoms.get_atomic_numbers()) + #individual.set_charges(atoms.get_charges()) + #individual.set_chemical_symbols(atoms.get_chemical_symbols()) + #individual.set_initial_magnetic_moments(atoms.get_initial_magnetic_moments()) + #individual.set_masses(atoms.get_masses()) + #individual.set_momenta(atoms.get_momenta()) + #individual.set_positions(atoms.get_positions()) + #individual.set_scaled_positions(atoms.get_scaled_positions()) + #individual.set_tags(atoms.get_tags()) + #individual.set_velocities(atoms.get_velocities()) + + individual.set_pbc(True) + with open(filename) as of: + of.readline() # number of atoms + sizes = of.readline() # comment == box size + sizes = sizes.split()[:3] + sizes = [float(x) for x in sizes] + cell = np.identity(3) + np.fill_diagonal(cell, sizes) + individual.set_cell(cell) + return None diff --git a/structopt/common/individual/generators/read_xyz.py b/structopt/common/individual/generators/read_xyz.py new file mode 100644 index 0000000..24a50b7 --- /dev/null +++ b/structopt/common/individual/generators/read_xyz.py @@ -0,0 +1,4 @@ +import structopt.fileio + +def read_xyz(filename): + return structopt.fileio.read_xyz(filename) diff --git a/structopt/common/individual/mutations/__init__.py b/structopt/common/individual/mutations/__init__.py index 33eecca..e16324c 100644 --- a/structopt/common/individual/mutations/__init__.py +++ b/structopt/common/individual/mutations/__init__.py @@ -2,13 +2,11 @@ import random import structopt -from .add_atoms import add_atoms -from .remove_atoms import remove_atoms -from .remove_surface_atoms import remove_surface_atoms -from .lattice_alteration import lattice_alteration -from .lattice_alteration_group import lattice_alteration_group -from .rotation import rotation -from .rotation_geo import rotation_geo +from .swap_positions import swap_positions +from .swap_species import swap_species +from .move_atoms import move_atoms +from .rotate_atoms import rotate_atoms +from .rotate_cluster import rotate_cluster class Mutations(object): @@ -28,37 +26,27 @@ def mutate(self, individual): return self.selected_mutation(individual) @staticmethod - @functools.wraps(add_atoms) - def add_atoms(individual): - return add_atoms(individual) + @functools.wraps(swap_positions) + def swap_positions(individual): + return swap_positions(individual) @staticmethod - @functools.wraps(remove_atoms) - def remove_atoms(individual): - return remove_atoms(individual) + @functools.wraps(swap_species) + def swap_species(individual): + return swap_species(individual) @staticmethod - @functools.wraps(remove_surface_atoms) - def remove_surface_atoms(individual): - return remove_surface_atoms(individual) + @functools.wraps(move_atoms) + def move_atoms(individual): + return move_atoms(individual) @staticmethod - @functools.wraps(lattice_alteration) - def lattice_alteration(individual): - return lattice_alteration(individual) + @functools.wraps(rotate_atoms) + def rotate_atoms(individual): + return rotate_atoms(individual) @staticmethod - @functools.wraps(lattice_alteration_group) - def lattice_alteration_group(individual): - return lattice_alteration_group(individual) - - @staticmethod - @functools.wraps(rotation) - def rotation(individual): - return rotation(individual) - - @staticmethod - @functools.wraps(rotation_geo) - def rotation_geo(individual): - return rotation_geo(individual) + @functools.wraps(rotate_cluster) + def rotate_cluster(individual): + return rotate_cluster(individual) diff --git a/structopt/common/individual/mutations/add_atoms.py b/structopt/common/individual/mutations/add_atoms.py deleted file mode 100644 index cbdb2ec..0000000 --- a/structopt/common/individual/mutations/add_atoms.py +++ /dev/null @@ -1,2 +0,0 @@ -def add_atoms(individual): - return individual diff --git a/structopt/common/individual/mutations/lattice_alteration.py b/structopt/common/individual/mutations/lattice_alteration.py deleted file mode 100644 index e324e5b..0000000 --- a/structopt/common/individual/mutations/lattice_alteration.py +++ /dev/null @@ -1,2 +0,0 @@ -def lattice_alteration(individual): - return individual diff --git a/structopt/common/individual/mutations/lattice_alteration_group.py b/structopt/common/individual/mutations/lattice_alteration_group.py deleted file mode 100644 index 64ee03a..0000000 --- a/structopt/common/individual/mutations/lattice_alteration_group.py +++ /dev/null @@ -1,2 +0,0 @@ -def lattice_alteration_group(individual): - return individual diff --git a/structopt/common/individual/mutations/move_atoms.py b/structopt/common/individual/mutations/move_atoms.py new file mode 100644 index 0000000..e571d87 --- /dev/null +++ b/structopt/common/individual/mutations/move_atoms.py @@ -0,0 +1,30 @@ +import random +import numpy as np + + +def move_atoms(individual, max_natoms=0.20): + """Randomly moves atoms within the individual (in place). + + Args: + individual (Individual): an individual + max_natoms (float or int): if float, the maximum number of atoms that will be moved is max_natoms*len(individual) + if int, the maximum number of atoms that will be moved is max_natoms + default: 0.20 + """ + if len(individual): + cell_max = np.maximum.reduce(individual.get_cell()) + cell_min = np.minimum.reduce(individual.get_cell()) + + if isinstance(max_natoms, float): + max_natoms = int(len(individual)*max_natoms) + max_natoms = max(max_natoms, 1) + natoms_to_move = random.randint(1, max_natoms) + atom_indices = list(range(len(individual))) + random.shuffle(atom_indices) # Using random.shuffle on the indices guarantees no duplicates + atom_indices = atom_indices[:natoms_to_move] + for index in atom_indices: + atom = individual[index] + atom.x, atom.y, atom.z = (random.uniform(cell_min[0], cell_max[0]), + random.uniform(cell_min[1], cell_max[1]), + random.uniform(cell_min[2], cell_max[2])) + return None diff --git a/structopt/common/individual/mutations/remove_atoms.py b/structopt/common/individual/mutations/remove_atoms.py deleted file mode 100644 index 944d07b..0000000 --- a/structopt/common/individual/mutations/remove_atoms.py +++ /dev/null @@ -1,2 +0,0 @@ -def remove_atoms(individual): - return individual diff --git a/structopt/common/individual/mutations/remove_surface_atoms.py b/structopt/common/individual/mutations/remove_surface_atoms.py deleted file mode 100644 index bc6afbd..0000000 --- a/structopt/common/individual/mutations/remove_surface_atoms.py +++ /dev/null @@ -1,2 +0,0 @@ -def remove_surface_atoms(individual): - return individual diff --git a/structopt/common/individual/mutations/rotate_atoms.py b/structopt/common/individual/mutations/rotate_atoms.py new file mode 100644 index 0000000..cb91ffc --- /dev/null +++ b/structopt/common/individual/mutations/rotate_atoms.py @@ -0,0 +1,30 @@ +import random +import numpy as np + + +def rotate_atoms(individual, max_natoms=0.20): + """Randomly rotates a number of random atoms within the individual (in place). + + Args: + individual (Individual): an individual + max_natoms (float or int): if float, the maximum number of atoms that will be rotated is max_natoms*len(individual) + if int, the maximum number of atoms that will be rotated is max_natoms + default: 0.20 + """ + if len(individual): + if isinstance(max_natoms, float): + max_natoms = int(len(individual)*max_natoms) + max_natoms = max(max_natoms, 1) + natoms_to_rotate = random.randint(1, max_natoms) + atom_indices = list(range(len(individual))) + random.shuffle(atom_indices) # Using random.shuffle on the indices guarantees no duplicates + atom_indices = atom_indices[:natoms_to_rotate] + atoms = individual[atom_indices] # This creates a copy of the atoms I think + del individual[atom_indices] + + axis = random.choice(['x', '-x', 'y', '-y', 'z', '-z']) + angle = random.uniform(30, 180) + atoms.rotate(axis, a=angle, center='COM', rotate_cell=False) + individual.extend(atoms) + return None + diff --git a/structopt/common/individual/mutations/rotate_cluster.py b/structopt/common/individual/mutations/rotate_cluster.py new file mode 100644 index 0000000..87355a1 --- /dev/null +++ b/structopt/common/individual/mutations/rotate_cluster.py @@ -0,0 +1,39 @@ +import random +import numpy as np +from ase import Atom + + +def rotate_cluster(individual, max_natoms=0.20): + """Randomly rotates a random cluster of atoms within the individual (in place). + + Args: + individual (Individual): an individual + max_natoms (float or int): if float, the maximum number of atoms that will be rotated is max_natoms*len(individual) + if int, the maximum number of atoms that will be rotated is max_natoms + default: 0.20 + """ + if len(individual): + cell_max = np.maximum.reduce(individual.get_cell()) + cell_min = np.minimum.reduce(individual.get_cell()) + + if isinstance(max_natoms, float): + max_natoms = int(len(individual)*max_natoms) + max_natoms = max(max_natoms, 1) + natoms_to_rotate = random.randint(1, max_natoms) + + point = (random.uniform(cell_min, cell_max), + random.uniform(cell_min, cell_max), + random.uniform(cell_min, cell_max)) + atom = Atom('Si', point) + nearest_indices = individual.get_nearest_atom_indices(atom_index=atom.index, count=natoms_to_rotate) + + + atoms = individual[nearest_indices] # This creates a copy of the atoms I think + del individual[nearest_indices] + + axis = random.choice(['x', '-x', 'y', '-y', 'z', '-z']) + angle = random.uniform(30, 180) + atoms.rotate(axis, a=angle, center='COM', rotate_cell=False) + individual.extend(atoms) + return None + diff --git a/structopt/common/individual/mutations/rotation.py b/structopt/common/individual/mutations/rotation.py deleted file mode 100644 index 401b7c5..0000000 --- a/structopt/common/individual/mutations/rotation.py +++ /dev/null @@ -1,2 +0,0 @@ -def rotation(individual): - return individual diff --git a/structopt/common/individual/mutations/rotation_geo.py b/structopt/common/individual/mutations/rotation_geo.py deleted file mode 100644 index be67064..0000000 --- a/structopt/common/individual/mutations/rotation_geo.py +++ /dev/null @@ -1,2 +0,0 @@ -def rotation_geo(individual): - return individual diff --git a/structopt/common/individual/mutations/swap_positions.py b/structopt/common/individual/mutations/swap_positions.py new file mode 100644 index 0000000..2fe0fe1 --- /dev/null +++ b/structopt/common/individual/mutations/swap_positions.py @@ -0,0 +1,38 @@ +import random +import numpy as np + + +def swap_positions(individual, max_natoms=0.20): + """Randomly swaps the positions atoms within the individual (in place). + + Args: + individual (Individual): an individual + max_natoms (float or int): if float, the maximum number of atoms whose positions will be swapped is max_natoms*len(individual) + if int, the maximum number of atoms whose positions will be swapped is max_natoms + if the number of atoms to be swapped is (or evaluates to) an odd integer, it is rounded down to an even integer + max_natoms corresponds to the maximum number of atoms whose positions will change + default: 0.20 + """ + if len(individual): + cell_max = np.maximum.reduce(individual.get_cell()) + cell_min = np.minimum.reduce(individual.get_cell()) + + if isinstance(max_natoms, float): + max_natoms = int(len(individual)*max_natoms) + max_natoms = max(max_natoms, 2) + natoms_to_move = random.randint(2, max_natoms) + if natoms_to_move % 2: # Make even + natoms_to_move -= 1 + atom_indices = list(range(len(individual))) + random.shuffle(atom_indices) # Using random.shuffle on the indices guarantees no duplicates + atom_indices = atom_indices[:natoms_to_move] + # Convert e.g. [1, 2, 3, 4, 5, 6, 7, 8] to [(1, 2), (3, 4), (5, 6), (7, 8)] + pairs = [(atom_indices[2*i], atom_indices[2*i+1]) for i in range(max_natoms//2)] + for index1, index2 in pairs: + atom1 = individual[index1] + atom2 = individual[index2] + atom1.x, atom2.x = atom2.x, atom1.x + atom1.y, atom2.y = atom2.y, atom1.y + atom1.z, atom2.z = atom2.z, atom1.z + return None + diff --git a/structopt/common/individual/mutations/swap_species.py b/structopt/common/individual/mutations/swap_species.py new file mode 100644 index 0000000..b248560 --- /dev/null +++ b/structopt/common/individual/mutations/swap_species.py @@ -0,0 +1,36 @@ +import random +import numpy as np + + +def swap_species(individual, max_natoms=0.20): + """Randomly swaps the species of atoms within the individual (in place). + + Args: + individual (Individual): an individual + max_natoms (float or int): if float, the maximum number of atoms that will be swapped is max_natoms*len(individual) + if int, the maximum number of atoms that will be swapped is max_natoms + if the number of atoms to be swapped is (or evaluates to) an odd integer, it is rounded down to an even integer + max_natoms corresponds to the maximum number of atoms whose species will change + default: 0.20 + """ + if len(individual): + cell_max = np.maximum.reduce(individual.get_cell()) + cell_min = np.minimum.reduce(individual.get_cell()) + + if isinstance(max_natoms, float): + max_natoms = int(len(individual)*max_natoms) + max_natoms = max(max_natoms, 2) + natoms_to_move = random.randint(2, max_natoms) + if natoms_to_move % 2: # Make even + natoms_to_move -= 1 + atom_indices = list(range(len(individual))) + random.shuffle(atom_indices) # Using random.shuffle on the indices guarantees no duplicates + atom_indices = atom_indices[:natoms_to_move] + # Convert e.g. [1, 2, 3, 4, 5, 6, 7, 8] to [(1, 2), (3, 4), (5, 6), (7, 8)] + pairs = [(atom_indices[2*i], atom_indices[2*i+1]) for i in range(max_natoms//2)] + for index1, index2 in pairs: + atom1 = individual[index1] + atom2 = individual[index2] + atom1.symbol, atom2.symbol = atom2.symbol, atom1.symbol + return None + diff --git a/structopt/common/individual/relaxations/LAMMPS.py b/structopt/common/individual/relaxations/LAMMPS.py index 9a5b9ce..17b694c 100644 --- a/structopt/common/individual/relaxations/LAMMPS.py +++ b/structopt/common/individual/relaxations/LAMMPS.py @@ -1,11 +1,123 @@ +import logging +import os +import shutil + +import structopt +import structopt.tools.lammps + + class LAMMPS(object): + def __init__(self): + self.parameters = structopt.parameters.relaxations.LAMMPS + + def get_command(self, individual): - return 'who' + raise NotImplementedError - def get_energy(self, individual): - return 0.0 def relax(self, individual): - #subprocess.call(command, shell=True) + logger = logging.getLogger('by-rank') + structure = individual #structure = compose_structure(individual) # TODO + calculator = self.setup_lammps(self.parameters, relax=True) + structure.set_calculator(calculator) + structure.set_pbc(True) + + # Perform Energy Minimization + try: + cwd = os.getcwd() + # Run LAMMPS + result = structure.calc.calculate(structure) + new_structure = result['atoms'] # TODO This is not a full Individual object, it's just an ASE atoms object I think; need to modify the input structure I think? + new_structure.set_pbc(True) + pea = result['pea'] + energy = result['thermo'][-1]['pe'] + pressure = 0 # should be modified if enthalpy_fit TODO I don't know what this means (bc of divide by zero error maybe?) + volume = new_structure.get_volume() + logger.info('Finished relaxation of individual{0} @ rank {1}: energy = {2}'.format(individual.index, structopt.parameters.globals.rank, energy)) + except Exception as error: + logger.critical('Error in energy evaluation: {0}'.format(error), exc_info=True) + # Copy files to TroubledLammps directory + path = os.path.join(cwd,'TroubledLammps') + if not os.path.exists(path): + os.mkdir(path) + calc = structure.get_calculator() + shutil.copyfile(calc.lammps_traj, os.path.join(path, os.path.basename(calc.lammps_traj))) + shutil.copyfile(calc.lammps_in, os.path.join(path, os.path.basename(calc.lammps_in))) + shutil.copyfile(calc.lammps_log, os.path.join(path, os.path.basename(calc.lammps_log))) + shutil.copyfile(calc.lammps_data, os.path.join(path, os.path.basename(calc.lammps_data))) + raise RuntimeError from error + + #individual, buli = decompose_structure(new_structure, individual) # TODO + buli = None + individual = structure + + individual.buli = buli + individual.energy = energy + individual.pressure = pressure + individual.volume = volume + + structure.calc.clean() + + #if relax: return individual + #else: + # return energy + + + def setup_lammps(self, parameters, relax): + logger = logging.getLogger('by-rank') + # TODO I only copied one of the pair_style options + if parameters["pair_style"] == 'eam/alloy': + parcoff = '* * {0}'.format(parameters["pot_file"]) + for one in atomlist: + parcoff += ' {0}'.format(one[0]) + pair_coeff = [parcoff] + lammps_parameters = {'pair_style': parameters["pair_style"], + 'pair_coeff': pair_coeff} + files = [parameters["potential_file"]] + elif parameters["pair_style"] == 'eam': + pair_coeff = ['* * {0}'.format(parameters["potential_file"])] + lammps_parameters = {'pair_style': parameters["pair_style"], + 'pair_coeff': pair_coeff} + files = [parameters["potential_file"]] + + if parameters["minimize"] != None: + try: + lammps_parameters['mass'][len(lammps_parameters['mass'])-1] += '\nmin_style {0}'.format(parameters["min_style"]) + except KeyError: + lammps_parameters['pair_coeff'][0] += '\nmin_style {0}'.format(parameters["min_style"]) + lammps_parameters['minimize'] = parameters["minimize"] + + if not relax: + lammps_parameters['minimize'] = "1e-8 1e-8 0 0" # Set maximim steps to 0 + lammps_parameters['thermosteps'] = parameters["thermo_steps"] + + # Set up kwargs for LAMMPS calculator + kwargs = { + 'parameters': lammps_parameters + } + + if files != None: + kwargs['files'] = files + + if parameters["keep_files"]: + # Set up directory for saving files + #path = os.path.join(os.getcwd(), '{0}-rank0'.format(structopt.parameters.globals.output_filename)) + path = os.path.join(os.getcwd(), structopt.parameters.globals.output_filename) + rank = structopt.parameters.globals.rank + logger.debug('Setting up directory for keeping LAMMPS files') + if structopt.parameters.globals.USE_MPI4PY: + if not os.path.exists(os.path.join(path, 'LAMMPSFiles')): + os.mkdir(os.path.join(path, 'LAMMPSFiles')) + logger.info('Making directory: {0}'.format(os.path.join(path, 'LAMMPSFiles'))) + + # Update kwargs + kwargs['keep_tmp_files'] = True + if structopt.parameters.globals.USE_MPI4PY: + tmp_dir = os.path.join(os.path.join(path, 'LAMMPSFiles'), 'rank-{0}'.format(rank)) + else: + tmp_dir = os.path.join(path, 'LAMMPSFiles') + kwargs['tmp_dir'] = tmp_dir + + return structopt.tools.lammps.LAMMPS(**kwargs) diff --git a/structopt/common/individual/relaxations/__init__.py b/structopt/common/individual/relaxations/__init__.py index d87c226..8559637 100644 --- a/structopt/common/individual/relaxations/__init__.py +++ b/structopt/common/individual/relaxations/__init__.py @@ -1,4 +1,5 @@ import structopt +from . import LAMMPS class Relaxations(object): diff --git a/structopt/common/population/__init__.py b/structopt/common/population/__init__.py index 8e9ab20..2abd627 100644 --- a/structopt/common/population/__init__.py +++ b/structopt/common/population/__init__.py @@ -32,7 +32,7 @@ def __init__(self): # Generate/load initial structures for structure_information in structopt.parameters.generators.initializers: for i in range(structure_information.number_of_individuals): - structure = Structure(**structure_information.data) + structure = Structure(index=i, **structure_information.data) self.append(structure) diff --git a/structopt/common/population/fitnesses/FEMSIM.py b/structopt/common/population/fitnesses/FEMSIM.py index 7362f5f..ac5c03e 100644 --- a/structopt/common/population/fitnesses/FEMSIM.py +++ b/structopt/common/population/fitnesses/FEMSIM.py @@ -17,7 +17,7 @@ def fitness(population): # Run the parallelized `mpiexec` command commands = ['-np 1 {command}'.format(command=command) for command in commands] # TODO: correctly allocate cores command = 'mpiexec {}'.format(' : '.join(commands)) - subprocess.call(command, shell=True) + subprocess.call(command, shell=True, stdout=subprocess.DEVNULL) # Collect the results for each chisq and return them for i, individual in enumerate(population): diff --git a/structopt/common/population/fitnesses/LAMMPS.py b/structopt/common/population/fitnesses/LAMMPS.py index f17b72e..f2fa423 100644 --- a/structopt/common/population/fitnesses/LAMMPS.py +++ b/structopt/common/population/fitnesses/LAMMPS.py @@ -15,13 +15,13 @@ def fitness(population): for i, individual in enumerate(population): if structopt.parameters.globals.USE_MPI4PY and structopt.parameters.globals.rank == i: command = individual.fitnesses.LAMMPS.get_command(individual) - subprocess.call(command, shell=True) + subprocess.call(command, shell=True, stdout=subprocess.DEVNULL) energy = individual.fitnesses.LAMMPS.get_energy(individual) logger.info('Individual {0} for LAMPPS evaluation had energy {1}'.format(i, energy)) energies = MPI.COMM_WORLD.gather(energy, root=0) else: command = individual.fitnesses.LAMMPS.get_command(individual) - subprocess.call(command, shell=True) + subprocess.call(command, shell=True, stdout=subprocess.DEVNULL) energy = individual.fitnesses.LAMMPS.get_energy(individual) logger.info('Individual {0} for LAMPPS evaluation had energy {1}'.format(i, energy)) energies.append(energy) diff --git a/structopt/common/population/relaxations/LAMMPS.py b/structopt/common/population/relaxations/LAMMPS.py index 2a5547c..0a3fef7 100644 --- a/structopt/common/population/relaxations/LAMMPS.py +++ b/structopt/common/population/relaxations/LAMMPS.py @@ -1,13 +1,11 @@ import subprocess -from mpi4py import MPI import structopt def relax(population): + # TODO update this so that it correctly distributes which core relaxes which individual for i, individual in enumerate(population): if structopt.parameters.globals.rank == i: - print(individual.relaxations.LAMMPS) individual.relaxations.LAMMPS.relax(individual) - return None diff --git a/structopt/crystal/__init__.py b/structopt/crystal/__init__.py index d1e5036..2e2f255 100644 --- a/structopt/crystal/__init__.py +++ b/structopt/crystal/__init__.py @@ -2,5 +2,5 @@ class Crystal(Individual): - """A crystal stucture.""" + """A stucture with periodic boundary conditions.""" diff --git a/structopt/fileio/__init__.py b/structopt/fileio/__init__.py index 276cd57..620b12b 100644 --- a/structopt/fileio/__init__.py +++ b/structopt/fileio/__init__.py @@ -1 +1,2 @@ from . import parameters, logger_utils +from .read_xyz import read_xyz diff --git a/structopt/fileio/parameters.py b/structopt/fileio/parameters.py index b2ea4a0..a59579d 100644 --- a/structopt/fileio/parameters.py +++ b/structopt/fileio/parameters.py @@ -33,7 +33,6 @@ def set_default(parameters): parameters.globals.setdefault('seed', None) - print(parameters.globals) if 'fitnesses' not in parameters or not parameters['fitnesses']: raise ValueError('Fitnesses must be specified in the parameter file.') diff --git a/structopt/fileio/read_xyz.py b/structopt/fileio/read_xyz.py new file mode 100644 index 0000000..48fe1ce --- /dev/null +++ b/structopt/fileio/read_xyz.py @@ -0,0 +1,6 @@ +import ase.io + + +def read_xyz(filename, index=None, format=None, **kwargs): + return ase.io.read(filename, index, format, **kwargs) + diff --git a/structopt/tools/__init__.py b/structopt/tools/__init__.py index e69de29..8a2e7b8 100644 --- a/structopt/tools/__init__.py +++ b/structopt/tools/__init__.py @@ -0,0 +1 @@ +from . import lammps diff --git a/structopt/tools/lammps.py b/structopt/tools/lammps.py new file mode 100644 index 0000000..b15fdf3 --- /dev/null +++ b/structopt/tools/lammps.py @@ -0,0 +1,822 @@ +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 . + + +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 + +__all__ = ['LAMMPS', 'write_lammps_data'] + +# "End mark" used to indicate that the calculation is done +CALCULATION_END_MARK = '__end_of_ase_invoked_calculation__' + +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 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 = [] + + self.pea = [] + + 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 -") + print(" translation to proper LAMMPS input cell might fail") + cell = self.atoms.get_cell() + self.prism = prism(cell) + self.run() + + return {'thermo': self.thermo_content, 'atoms': self.atoms, 'pea': self.pea} + + 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): + """Method which explicitely 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 = '{}{:>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() + + # Save the names of the files + self.lammps_trj = lammps_trj_fd.name + self.lammps_in = lammps_in + self.lammps_log = lammps_log + self.lammps_data = lammps_data + + # 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 = special_tee(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 = special_tee(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=True) # Changed from ASE version: Added `set_atoms=True` + 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 == None): + lammps_data = 'data.' + self.label + write_lammps_data(lammps_data, self.atoms, self.specorder, + force_skew=self.always_triclinic, prismobj=self.prism) + + def write_lammps_in(self, lammps_in=None, lammps_trj=None, lammps_data=None): + """Method which writes a LAMMPS in file with run parameters and settings.""" + + if isinstance(lammps_in, str): + f = paropen(lammps_in, 'w') + 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 "{}"\n' + 'variable data_file string "{}"\n'.format(lammps_trj, lammps_data).encode('utf-8')) + + parameters = self.parameters + pbc = self.atoms.get_pbc() + f.write('units metal \n'.encode('utf-8')) + if ('boundary' in parameters): + f.write('boundary {} \n'.format(parameters['boundary']).encode('utf-8')) + else: + f.write('boundary {} {} {} \n'.format(*('sp'[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('{} {} \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(['# {.16} {.16} {.16}\n'.format(*x) for x in self.atoms.get_cell()]).encode('utf-8')) + + p = self.prism + f.write('lattice sc 1.0\n'.encode('utf-8')) + xhi, yhi, zhi, xy, xz, yz = p.get_lammps_prism_str() + if self.always_triclinic or p.is_skewed(): + f.write('region asecell prism 0.0 {} 0.0 {} 0.0 {} '.format(xhi, yhi, zhi).encode('utf-8')) + f.write('{} {} {} side in units box\n'.format(xy, xz, yz).encode('utf-8')) + else: + f.write('region asecell block 0.0 {} 0.0 {} 0.0 {} 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 {} 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: {.16} {.16} {.16}\n'.format(*tuple(pos)).encode('utf-8')) + f.write('create_atoms {} single {} {} {} units box\n'.format(*((species_i[s],)+p.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 {}\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 {} \n'.format(pair_style).encode('utf-8')) + for pair_coeff in parameters['pair_coeff']: + f.write('pair_coeff {} \n'.format(pair_coeff).encode('utf-8')) + if 'mass' in parameters: + for mass in parameters['mass']: + f.write('mass {} \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 'thermosteps' in parameters: + f.write('thermo_style custom {}\n' + 'thermo_modify flush yes\n' + 'thermo {}\n'.format(' '.join(self._custom_thermo_args), repr(parameters['thermosteps'])).encode('utf-8')) + else: + f.write('thermo_style custom {}\n' + 'thermo_modify flush yes\n' + 'thermo 1\n'.format(' '.join(self._custom_thermo_args)).encode('utf-8')) + + f.write('\n### run\n'.encode('utf-8')) + if 'lammps_command' in parameters: + f.write('fix fix_nve all nve\n'.encode('utf-8')) + f.write('minimize {}\n'.format(parameters['minimize']).encode('utf-8')) + f.write('unfix fix_nve \n'.encode('utf-8')) + + f.write('{} \n'.format(parameters['lammps_command']).encode('utf-8')) + f.write('minimize {}\n'.format(parameters['minimize']).encode('utf-8')) + else: + f.write('fix fix_nve all nve\n'.encode('utf-8')) + + if 'minimize' in parameters: + f.write('minimize {}\n'.format(parameters['minimize']).encode('utf-8')) + if 'run' in parameters: + f.write('run {}\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(('thermo_style custom {}\n' + # 'thermo_modify flush yes\n' + # 'thermo 1\n').format(' '.join(self._custom_thermo_args)).encode('utf-8')) + f.write('compute pea all pe/atom \n'.encode('utf-8')) + #f.write('dump dump_all all custom 1 {} id type x y z vx vy vz fx fy fz\n'.format(lammps_trj).encode('utf-8')) # This was up higher + #f.write('dump dump_all all custom 2 {} id type x y z vx vy vz fx fy fz\n'.format(lammps_trj).encode('utf-8')) + f.write('dump dump_all all custom 2 {} id type x y z c_pea \n'.format(lammps_trj).encode('utf-8')) + #f.write('dump dump_all all custom 2 {} id type x y z \n'.format(lammps_trj.encode('utf-8')) + f.write('run 1\n'.encode('utf-8')) + + f.write('print "{}"\n'.format(CALCULATION_END_MARK).encode('utf-8')) + f.write('log /dev/stdout\n'.encode('utf-8')) # Force LAMMPS to flush log + f.write('undump dump_all\n'.encode('utf-8')) # Force LAMMPS to flush trj + + if close_in_file: + f.flush() + f.close() + + def read_lammps_log(self, lammps_log=None, PotEng_first=False): + """Method which reads a LAMMPS output log file.""" + + if (lammps_log == None): + lammps_log = self.label + '.log' + + if isinstance(lammps_log, str): + f = paropen(lammps_log, 'w') + 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 == None): + lammps_trj = self.label + '.lammpstrj' + + f = paropen(lammps_trj, 'r') + nlines = 0 + while True: + line = f.readline() + + if not line: + break + nlines += 1 + + #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 = [] ; pea = [] #; 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'] ] ) + pea.append( [ float(fields[atom_attributes[x]]) for x in ['c_pea'] ] ) + #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'] ] ) + f.close() + + if nlines==0: + raise RuntimeError('Failed to retreive any output from trajectory file: {0}\n'.format(lammps_trj)) + + # 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]] + + # assume that LAMMPS does not reorder atoms internally + cell_atoms = np.array(cell) + type_atoms = np.array(type) + + if self.atoms: + #cell_atoms = self.atoms.get_cell() # Changed from ASE version: commented + + # 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.array( [np.dot(np.array(r), rotation_lammps2ase) for r in positions] ) + #velocities_atoms = np.array( [np.dot(np.array(v), rotation_lammps2ase) for v in velocities] ) + #forces_atoms = np.array( [np.dot(np.array(f), rotation_lammps2ase) for f in forces] ) + + if (set_atoms): + # assume periodic boundary conditions here (like also below in write_lammps) + self.atoms = Atoms(type_atoms, positions=positions_atoms, cell=cell_atoms) + self.pea = pea + + #self.forces = forces_atoms + + +class special_tee: + """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: + 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 precission 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), return a tuple of str" + # Two-stage fold, first into box, then into semi-open interval + # (within the given precission). + 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 pos_to_lammps_str(self, position): + "Rotate an ase-cell position to the lammps cell orientation, return tuple of strs" + return tuple([self.f2s(x) for x in np.dot(position, self.R)]) + + def pos_to_lammps_fold_str(self, position): + "Rotate and fold an ase-cell position into the lammps cell, return tuple of strs" + 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): + """Method which writes atomic structure data to a LAMMPS data file.""" + if isinstance(fileobj, str): + f = paropen(fileobj, 'w') + 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('{} (written by ASE) \n\n'.format(f.name).encode('utf-8')) + + symbols = atoms.get_chemical_symbols() + n_atoms = len(symbols) + f.write('{} \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 predictively 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('{} 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 {} xlo xhi\n'.format(xhi).encode('utf-8')) + f.write('0.0 {} ylo yhi\n'.format(yhi).encode('utf-8')) + f.write('0.0 {} zlo zhi\n'.format(zhi).encode('utf-8')) + + if force_skew or p.is_skewed(): + f.write('{} {} {} 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(map(p.pos_to_lammps_str, + atoms.get_positions())): + s = species.index(symbols[i]) + 1 + f.write('{:>6} {:>3} {} {} {}\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()): + f.write('{:>6} {} {} {}\n'.format(*(i+1,)+tuple(v)).encode('utf-8')) + + 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 = {}'.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 : {} , total energy : {}'.format(a, bulk.get_potential_energy())) + + calc.clean() +