# coding: utf-8
# Copyright (c) Pymatgen Development Team.
# Distributed under the terms of the MIT License.
"""
This module implements an interface to the VAMPIRE code for atomistic
simulations of magnetic materials.
This module depends on a compiled vampire executable available in the path.
Please download at https://vampire.york.ac.uk/download/ and
follow the instructions to compile the executable.
If you use this module, please cite the following:
"Atomistic spin model simulations of magnetic nanomaterials."
R. F. L. Evans, W. J. Fan, P. Chureemart, T. A. Ostler, M. O. A. Ellis
and R. W. Chantrell. J. Phys.: Condens. Matter 26, 103202 (2014)
"""
import subprocess
import logging
import pandas as pd
from monty.dev import requires
from monty.os.path import which
from monty.json import MSONable
from pymatgen.analysis.magnetism.heisenberg import HeisenbergMapper
__author__ = "ncfrey"
__version__ = "0.1"
__maintainer__ = "Nathan C. Frey"
__email__ = "ncfrey@lbl.gov"
__status__ = "Development"
__date__ = "June 2019"
VAMPEXE = which("vampire-serial")
[docs]class VampireCaller:
"""
Run Vampire on a material with magnetic ordering and exchange parameter information to compute the critical
temperature with classical Monte Carlo.
"""
@requires(
VAMPEXE,
"VampireCaller requires vampire-serial to be in the path."
"Please follow the instructions at https://vampire.york.ac.uk/download/.",
)
def __init__(
self,
ordered_structures=None,
energies=None,
mc_box_size=4.0,
equil_timesteps=2000,
mc_timesteps=4000,
save_inputs=False,
hm=None,
avg=True,
user_input_settings=None,
):
"""
user_input_settings is a dictionary that can contain:
* start_t (int): Start MC sim at this temp, defaults to 0 K.
* end_t (int): End MC sim at this temp, defaults to 1500 K.
* temp_increment (int): Temp step size, defaults to 25 K.
Args:
ordered_structures (list): Structure objects with magmoms.
energies (list): Energies of each relaxed magnetic structure.
mc_box_size (float): x=y=z dimensions (nm) of MC simulation box
equil_timesteps (int): number of MC steps for equilibrating
mc_timesteps (int): number of MC steps for averaging
save_inputs (bool): if True, save scratch dir of vampire input files
hm (HeisenbergModel): object already fit to low energy
magnetic orderings.
avg (bool): If True, simply use <J> exchange parameter estimate.
If False, attempt to use NN, NNN, etc. interactions.
user_input_settings (dict): optional commands for VAMPIRE Monte Carlo
Parameters:
sgraph (StructureGraph): Ground state graph.
unique_site_ids (dict): Maps each site to its unique identifier
nn_interacations (dict): {i: j} pairs of NN interactions
between unique sites.
ex_params (dict): Exchange parameter values (meV/atom)
mft_t (float): Mean field theory estimate of critical T
mat_name (str): Formula unit label for input files
mat_id_dict (dict): Maps sites to material id # for vampire
indexing.
TODO:
* Create input files in a temp folder that gets cleaned up after run terminates
"""
self.mc_box_size = mc_box_size
self.equil_timesteps = equil_timesteps
self.mc_timesteps = mc_timesteps
self.save_inputs = save_inputs
self.avg = avg
if not user_input_settings: # set to empty dict
self.user_input_settings = {}
else:
self.user_input_settings = user_input_settings
# Get exchange parameters and set instance variables
if not hm:
hmapper = HeisenbergMapper(
ordered_structures, energies, cutoff=3.0, tol=0.02
)
hm = hmapper.get_heisenberg_model()
# Attributes from HeisenbergModel
self.hm = hm
self.structure = hm.structures[0] # ground state
self.sgraph = hm.sgraphs[0] # ground state graph
self.unique_site_ids = hm.unique_site_ids
self.nn_interactions = hm.nn_interactions
self.dists = hm.dists
self.tol = hm.tol
self.ex_params = hm.ex_params
self.javg = hm.javg
# Full structure name before reducing to only magnetic ions
self.mat_name = hm.formula
# Switch to scratch dir which automatically cleans up vampire inputs files unless user specifies to save them
# with ScratchDir('/scratch', copy_from_current_on_enter=self.save_inputs,
# copy_to_current_on_exit=self.save_inputs) as temp_dir:
# os.chdir(temp_dir)
# Create input files
self._create_mat()
self._create_input()
self._create_ucf()
# Call Vampire
process = subprocess.Popen(
["vampire-serial"], stdout=subprocess.PIPE, stderr=subprocess.PIPE
)
stdout, stderr = process.communicate()
stdout = stdout.decode()
if stderr:
vanhelsing = stderr.decode()
if len(vanhelsing) > 27: # Suppress blank warning msg
logging.warning(vanhelsing)
if process.returncode != 0:
raise RuntimeError(
"Vampire exited with return code {}.".format(process.returncode)
)
self._stdout = stdout
self._stderr = stderr
# Process output
nmats = max(self.mat_id_dict.values())
parsed_out, critical_temp = VampireCaller.parse_stdout("output", nmats)
self.output = VampireOutput(parsed_out, nmats, critical_temp)
def _create_mat(self):
structure = self.structure
mat_name = self.mat_name
magmoms = structure.site_properties["magmom"]
# Maps sites to material id for vampire inputs
mat_id_dict = {}
nmats = 0
for key in self.unique_site_ids:
spin_up, spin_down = False, False
nmats += 1 # at least 1 mat for each unique site
# Check which spin sublattices exist for this site id
for site in key:
m = magmoms[site]
if m > 0:
spin_up = True
if m < 0:
spin_down = True
# Assign material id for each site
for site in key:
m = magmoms[site]
if spin_up and not spin_down:
mat_id_dict[site] = nmats
if spin_down and not spin_up:
mat_id_dict[site] = nmats
if spin_up and spin_down:
# Check if spin up or down shows up first
m0 = magmoms[key[0]]
if m > 0 and m0 > 0:
mat_id_dict[site] = nmats
if m < 0 and m0 < 0:
mat_id_dict[site] = nmats
if m > 0 and m0 < 0:
mat_id_dict[site] = nmats + 1
if m < 0 and m0 > 0:
mat_id_dict[site] = nmats + 1
# Increment index if two sublattices
if spin_up and spin_down:
nmats += 1
mat_file = ["material:num-materials=%d" % (nmats)]
for key in self.unique_site_ids:
i = self.unique_site_ids[key] # unique site id
for site in key:
mat_id = mat_id_dict[site]
# Only positive magmoms allowed
m_magnitude = abs(magmoms[site])
if magmoms[site] > 0:
spin = 1
if magmoms[site] < 0:
spin = -1
atom = structure[i].species.reduced_formula
mat_file += ["material[%d]:material-element=%s" % (mat_id, atom)]
mat_file += [
"material[%d]:damping-constant=1.0" % (mat_id),
"material[%d]:uniaxial-anisotropy-constant=1.0e-24"
% (mat_id), # xx - do we need this?
"material[%d]:atomic-spin-moment=%.2f !muB" % (mat_id, m_magnitude),
"material[%d]:initial-spin-direction=0,0,%d" % (mat_id, spin),
]
mat_file = "\n".join(mat_file)
mat_file_name = mat_name + ".mat"
self.mat_id_dict = mat_id_dict
with open(mat_file_name, "w") as f:
f.write(mat_file)
def _create_input(self):
structure = self.structure
mcbs = self.mc_box_size
equil_timesteps = self.equil_timesteps
mc_timesteps = self.mc_timesteps
mat_name = self.mat_name
input_script = ["material:unit-cell-file=%s.ucf" % (mat_name)]
input_script += ["material:file=%s.mat" % (mat_name)]
# Specify periodic boundary conditions
input_script += [
"create:periodic-boundaries-x",
"create:periodic-boundaries-y",
"create:periodic-boundaries-z",
]
# Unit cell size in Angstrom
abc = structure.lattice.abc
ucx, ucy, ucz = abc[0], abc[1], abc[2]
input_script += ["dimensions:unit-cell-size-x = %.10f !A" % (ucx)]
input_script += ["dimensions:unit-cell-size-y = %.10f !A" % (ucy)]
input_script += ["dimensions:unit-cell-size-z = %.10f !A" % (ucz)]
# System size in nm
input_script += [
"dimensions:system-size-x = %.1f !nm" % (mcbs),
"dimensions:system-size-y = %.1f !nm" % (mcbs),
"dimensions:system-size-z = %.1f !nm" % (mcbs),
]
# Critical temperature Monte Carlo calculation
input_script += [
"sim:integrator = monte-carlo",
"sim:program = curie-temperature",
]
# Default Monte Carlo params
input_script += [
"sim:equilibration-time-steps = %d" % (equil_timesteps),
"sim:loop-time-steps = %d" % (mc_timesteps),
"sim:time-steps-increment = 1",
]
# Set temperature range and step size of simulation
if "start_t" in self.user_input_settings:
start_t = self.user_input_settings["start_t"]
else:
start_t = 0
if "end_t" in self.user_input_settings:
end_t = self.user_input_settings["end_t"]
else:
end_t = 1500
if "temp_increment" in self.user_input_settings:
temp_increment = self.user_input_settings["temp_increment"]
else:
temp_increment = 25
input_script += [
"sim:minimum-temperature = %d" % (start_t),
"sim:maximum-temperature = %d" % (end_t),
"sim:temperature-increment = %d" % (temp_increment),
]
# Output to save
input_script += [
"output:temperature",
"output:mean-magnetisation-length",
"output:material-mean-magnetisation-length",
"output:mean-susceptibility",
]
input_script = "\n".join(input_script)
with open("input", "w") as f:
f.write(input_script)
def _create_ucf(self):
structure = self.structure
mat_name = self.mat_name
abc = structure.lattice.abc
ucx, ucy, ucz = abc[0], abc[1], abc[2]
ucf = ["# Unit cell size:"]
ucf += ["%.10f %.10f %.10f" % (ucx, ucy, ucz)]
ucf += ["# Unit cell lattice vectors:"]
a1 = list(structure.lattice.matrix[0])
ucf += ["%.10f %.10f %.10f" % (a1[0], a1[1], a1[2])]
a2 = list(structure.lattice.matrix[1])
ucf += ["%.10f %.10f %.10f" % (a2[0], a2[1], a2[2])]
a3 = list(structure.lattice.matrix[2])
ucf += ["%.10f %.10f %.10f" % (a3[0], a3[1], a3[2])]
nmats = max(self.mat_id_dict.values())
ucf += ["# Atoms num_materials; id cx cy cz mat cat hcat"]
ucf += ["%d %d" % (len(structure), nmats)]
# Fractional coordinates of atoms
for site, r in enumerate(structure.frac_coords):
# Back to 0 indexing for some reason...
mat_id = self.mat_id_dict[site] - 1
ucf += ["%d %.10f %.10f %.10f %d 0 0" % (site, r[0], r[1], r[2], mat_id)]
# J_ij exchange interaction matrix
sgraph = self.sgraph
ninter = 0
for i, node in enumerate(sgraph.graph.nodes):
ninter += sgraph.get_coordination_of_site(i)
ucf += ["# Interactions"]
ucf += ["%d isotropic" % (ninter)]
iid = 0 # counts number of interaction
for i, node in enumerate(sgraph.graph.nodes):
connections = sgraph.get_connected_sites(i)
for c in connections:
jimage = c[1] # relative integer coordinates of atom j
dx = jimage[0]
dy = jimage[1]
dz = jimage[2]
j = c[2] # index of neighbor
dist = round(c[-1], 2)
# Look up J_ij between the sites
if self.avg is True: # Just use <J> estimate
j_exc = self.hm.javg
else:
j_exc = self.hm._get_j_exc(i, j, dist)
# Convert J_ij from meV to Joules
j_exc *= 1.6021766e-22
j_exc = str(j_exc) # otherwise this rounds to 0
ucf += ["%d %d %d %d %d %d %s" % (iid, i, j, dx, dy, dz, j_exc)]
iid += 1
ucf = "\n".join(ucf)
ucf_file_name = mat_name + ".ucf"
with open(ucf_file_name, "w") as f:
f.write(ucf)
[docs] @staticmethod
def parse_stdout(vamp_stdout, nmats):
"""Parse stdout from Vampire.
Args:
vamp_stdout (txt file): Vampire 'output' file.
nmats (int): Num of materials in Vampire sim.
Returns:
parsed_out (DataFrame): MSONable vampire output.
critical_temp (float): Calculated critical temp.
"""
names = (
["T", "m_total"]
+ ["m_" + str(i) for i in range(1, nmats + 1)]
+ ["X_x", "X_y", "X_z", "X_m", "nan"]
)
# Parsing vampire MC output
df = pd.read_csv(vamp_stdout, sep="\t", skiprows=9, header=None, names=names)
df.drop("nan", axis=1, inplace=True)
parsed_out = df.to_json()
# Max of susceptibility <-> critical temp
critical_temp = df.iloc[df.X_m.idxmax()]["T"]
return parsed_out, critical_temp
[docs]class VampireOutput(MSONable):
"""
This class processes results from a Vampire Monte Carlo simulation
and returns the critical temperature.
"""
def __init__(self, parsed_out=None, nmats=None, critical_temp=None):
"""
Args:
parsed_out (json): json rep of parsed stdout DataFrame.
nmats (int): Number of distinct materials (1 for each specie and up/down spin).
critical_temp (float): Monte Carlo Tc result.
"""
self.parsed_out = parsed_out
self.nmats = nmats
self.critical_temp = critical_temp