Sim controller

A central executable script to handle all other modules based on input settings files.

ase-test/FeCrAl.py

-import os
-import sys
-import time
-from math import *
-import cProfile
-
-import numpy as np
-
-import ase.lattice.cubic
-import ase.io
-import ase.units
-
-from sajuhak import KokkoPotential, SajuhakDynamics
-
-
-def format_time(seconds):
-    """
-    Takes a time in seconds as a float and formats it into a string with separated units.
-
-    :param seconds: time in seconds
-    :return: time as a formatted string
-    """
-    milliseconds = int((seconds - floor(seconds)) * 1000)
-    minutes, seconds = divmod(int(seconds), 60)
-    hours, minutes = divmod(int(minutes), 60)
-    if minutes != 0:
-        if hours != 0:
-            return "{h}h {m}m {s}s {ms}ms".format(h=hours, m=minutes, s=seconds, ms=milliseconds)
-        else:
-            return "{m}m {s}s {ms}ms".format(m=minutes, s=seconds, ms=milliseconds)
-    else:
-        return "{s}s {ms}ms".format(s=seconds, ms=milliseconds)
-
-
-start_time = time.time()
-print("Creating lattice...", end="")
-sys.stdout.flush()
-
-# file output
-sim_filename = 'FeCrAl-dyn.{i}.xyz'
-time_string = time.strftime("%Y%m%d-%H'%M'%S")
-out_folder = 'sim_output ' + time_string
-
-# ------------------------------
-# SIMULATION PARAMETERS
-# ------------------------------
-bulk_scale = 10
-x_scale = 1
-y_scale = 1
-z_scale = 1
-Cr_concentration = 0.10
-Al_concentration = 0.05
-O_concentration = 0.1
-O_density = 2.7e25  # atoms per m^3
-empty_space = 1/2
-n_frames = 100
-n_exchanges = 10000
-kT = 300
-
-# ------------------------------
-# CREATE LATTICE
-# ------------------------------
-bulk_size = [bulk_scale*x_scale, bulk_scale*y_scale, bulk_scale*z_scale]
-bulk = ase.lattice.cubic.BodyCenteredCubic('Fe', size=bulk_size, pbc=(True, True, False))
-
-# project the cell vectors onto x, y, and z axes to get the bounding box
-cell_vectors = bulk.get_cell()
-cell_dimensions = [[], [], []]
-for cell_dir in [(1, 0, 0), (0, 1, 0), (0, 0, 1)]:
-    for j in [0, 1, 2]:
-        cell_dimensions[j].append(np.dot(cell_vectors[j], cell_dir))
-for i, dim_list in enumerate(cell_dimensions):
-    cell_dimensions[i] = max(dim_list)
-
-# calculate bulk volume as a*(bxc)
-bulk_volume = np.dot(cell_vectors[0], np.cross(cell_vectors[1], cell_vectors[2]))
-volume_per_atom = bulk_volume/len(bulk)
-
-# initialize calculator with the nearest neighbor distance
-all_dist = bulk.get_all_distances()
-nbor_dist = np.min(all_dist[np.nonzero(all_dist)])
-nbor_cutoffs = [nbor_dist] * len(bulk)
-dyn = SajuhakDynamics(bulk, nbor_cutoffs, 'potentials.txt')
-bulk.set_calculator(KokkoPotential(nbor_cutoffs, 'potentials.txt'))
-
-# initialize atom modifications
-symbols = dyn.atoms.get_chemical_symbols()
-positions = dyn.atoms.get_positions()
-max_miller = bulk_scale * z_scale * (1 - empty_space)
-
-# replace atoms with empty spaces above a certain z position
-for i, position in enumerate(positions):
-    if position[2] > dyn.atoms.millerbasis[2][2]*max_miller:
-        symbols[i] = 'X'
-
-# replace the remaining atoms with dopants until desired concentration is reached
-# first store the indexes of the remaining atoms to exclude empty space
-Fe_indexes = [index for index, symbol in enumerate(symbols) if symbol == 'Fe']
-changed = True
-while changed:
-    changed = False
-    if symbols.count('Cr') / len(Fe_indexes) < Cr_concentration:
-        index = np.random.randint(len(Fe_indexes))
-        symbols[index] = 'Cr'
-        changed = True
-    if symbols.count('Al') / len(Fe_indexes) < Al_concentration:
-        index = np.random.choice(Fe_indexes)
-        symbols[index] = 'Al'
-        changed = True
-
-# add oxygen to the empty space until desired concentration is reached
-X_indexes = [index for index, symbol in enumerate(symbols) if symbol == 'X']
-O_atom_concentration = O_density / 1.0e30 * volume_per_atom
-while True:
-    if symbols.count('O') / len(X_indexes) < O_concentration:
-        index = np.random.choice(X_indexes)
-        symbols[index] = 'O'
-        continue
-    break
-
-dyn.atoms.set_chemical_symbols(symbols)
-
-# finished modifying the lattice
-end_time = time.time()
-particles_time = format_time(end_time-start_time)
-print("{particles} particles in {time}".format(particles=len(bulk), time=particles_time))
-print("")
-
-# ------------------------------
-# SIMULATION
-# ------------------------------
-print("starting simulation")
-n_steps = int(n_exchanges / n_frames)
-n_frames_width = len(str(n_frames - 1))
-if not os.path.exists(out_folder):
-    os.mkdir(out_folder)
-sim_path = os.path.join(out_folder, sim_filename.format(i='*' * n_frames_width))
-print("writing to files \"" + sim_path + "\"")
-# ase.io.write(sim_path, bulk)
-# profiler = cProfile.Profile()
-# profiler.enable()
-start_time = time.time()
-start_atoms = dyn.atoms.copy()
-print_time = time.time()
-for i in range(n_frames):
-    i = str(i).zfill(n_frames_width)
-    dyn.run(kT, n_steps)
-    sim_path = os.path.join(out_folder, sim_filename.format(i=i))
-    ase.io.write(sim_path, bulk)
-    if time.time() - print_time >= 0.5:
-        print("\rwrote frame {i}/{n}".format(i=i, n=str(n_frames)), end="")
-        print_time = time.time()
-print("\rwrote frame {}/{}".format(str(n_frames), str(n_frames)), flush=True)
-print("simulation done")
-print("")
-end_time = time.time()
-sim_time = format_time(end_time - start_time)
-print("Simulation took " + sim_time)
-# profiler.print_stats(sort=2)
-
-# ------------------------------
-# SIMULATION REPORT
-# ------------------------------
-n_atoms = len(bulk)
-n_Fe = bulk.get_chemical_symbols().count('Fe')
-n_Cr = bulk.get_chemical_symbols().count('Cr')
-n_Al = bulk.get_chemical_symbols().count('Al')
-sum_bulk = n_Fe + n_Cr + n_Al
-n_O = bulk.get_chemical_symbols().count('O')
-n_X = bulk.get_chemical_symbols().count('X')
-sum_vacuum = n_O + n_X
-
-report_filename = "sim_report.txt"
-report_path = os.path.join(out_folder, report_filename)
-
-report_file = open(report_path, 'w', encoding='utf-8')
-report_string = """Atoms: {atoms}
-Cell dimensions: {cell_x:.2f} Å x {cell_y:.2f} Å x {cell_z:.2f} Å
-    Bulk:
-        Fe: {n_Fe} ({perc_Fe:.2f}%)
-        Cr: {n_Cr} ({perc_Cr:.2f}%)
-        Al: {n_Al} ({perc_Al:.2f}%)
-    Vacuum:
-        O: {n_O} ({perc_O:.2f}%)
-        Empty space: {n_X} ({perc_X:.2f}%)
-Particles created in {particles_time}.
-
-Metropolis steps: {n_exchanges}
-kT: {kT}
-Frames: {n_frames}.
-Simulation took {sim_time}.
-""".format(atoms=n_atoms,
-           cell_x=cell_dimensions[0],
-           cell_y=cell_dimensions[1],
-           cell_z=cell_dimensions[2],
-
-           n_Fe=n_Fe, perc_Fe=n_Fe/sum_bulk*100,
-           n_Cr=n_Cr, perc_Cr=n_Cr/sum_bulk*100,
-           n_Al=n_Al, perc_Al=n_Al/sum_bulk*100,
-
-           n_O=n_O, perc_O=n_O/sum_vacuum*100,
-           n_X=n_X, perc_X=n_X/sum_vacuum*100,
-
-           particles_time=particles_time,
-
-           n_exchanges=n_exchanges,
-           kT=kT,
-           n_frames=n_frames,
-           sim_time=sim_time)
-report_file.write(report_string)
-report_file.close()