# -*- coding: utf-8 -*-
"""
Created on Wed Dec 3 14:25:06 2014
@author: Matti Ropo
@author: Henrik Levämäki
"""
from __future__ import print_function
import time
import os
import sys
import numpy as np
import pyemto.common.common as common
[docs]class System:
"""The main class which provides the basis for the pyEMTO scripts.
Somewhere in the beginning of a pyEMTO script a new instance of
the system class should be created. All subsequent communication
with the newly created system should be through the class methods,
which are described below.
:param folder: Main folder where the input and output files will
be stored. Use of absolute paths is recommended
(Default value = current working directory)
:type folder: str
:param EMTOdir: Path to the folder of the EMTO installation.
This entry can and should be modified by the user
inside the System.__init__ function
(Default value = /home/user/EMTO5.8)
:type EMTOdir: str
:param xc: Choice for the xc-functional can be set here.
(Default value = PBE)
:type xc: str
:returns: None
:rtype: None
"""
def __init__(self, folder=None, EMTOdir=None, xc=None):
# Import necessary packages
from pyemto.latticeinputs.latticeinputs import Latticeinputs
from pyemto.emtoinputs.emtoinputs import Emtoinputs
# Check input arguments
if folder is None:
self.folder = os.getcwd() # Use current folder
else:
self.folder = folder
if EMTOdir is None:
self.EMTOdir = "/home/hpleva/EMTO5.8"
else:
self.EMTOdir = EMTOdir
# Initialize default parameters
self.ca_range_default = np.linspace(1.50, 1.70, 7)
self.elastic_constants_points = 6
self.elastic_constants_deltas = np.linspace(0.0, 0.05,
self.elastic_constants_points)
self.RyBohr3_to_GPa = 14710.5065722
self.kappaw_default = [0.0, -20.0]
self.hcpo_relax_points = 5
self.hcpm_relax_points = 5
if xc is None:
self.xc = 'PBE'
else:
self.xc = xc
# Create working folders
common.check_folders(self.folder + '/kgrn', self.folder + '/kgrn/tmp',
self.folder + '/kfcd', self.folder + '/fit')
# BMDL, KSTR, SHAPE, KGRN and KFCD class instances
self.lattice = Latticeinputs()
self.emto = Emtoinputs()
return
[docs] def bulk(self, jobname=None, lat=None, atoms=None, concs=None, splts=None, sws=None,
latname=None, latpath=None, emtopath=None, ibz=None, bmod=None, xc=None, ca=None,
**kwargs):
"""Initializes the basic parameters for bulk systems.
Basic information concerning the system,
such as the types of atoms and the crystal structure should be given to this function and
it should be called right after the class instance has been created.
:param jobname: Name of the system (Default value = None)
:type jobname:
:param lat: The type of lattice structure (Default value = None)
:type lat:
:param atoms: List of atoms in the system (Default value = None)
:type atoms:
:param concs: List of concentrations of the elements in the
'atoms' list. This information is only used in CPA
calculations (Default value = None)
:type concs:
:param splts: List of initial magnetic moments of the elements in the
'atoms' list (Default value = None)
:type splts:
:param sws: The Wigner-Seitz radius of the system (Default value = None)
:type sws: float
:param latname: The 'jobname' of the BMDL, KSTR and SHAPE output files. These
structure output files have to be located in the 'latpath'
directory and they have to be named jobname.extention
(Default value = None)
:type latname:
:param latpath: The absolute path to the folder where the 'bmdl', 'kstr' and 'shape'
folders are located, which in turn contain the output files of the
structure calculation (Default value = None)
:type latpath:
:param emtopath: The absolute path to the folder where the EMTO installation is
located (Default value = None)
:type emtopath:
:param ibz: The code number indicating the Bravais lattice that the crystal
structure of the system has. For a list of possible values, please consult the
EMTO manual (Default value = None)
:type ibz:
:param bmod: The bulk modulus can be inputed here and if it is given,
it will be used by the elastic modulus routines (Default value = None)
:type bmod:
:param xc: The choice of the xc-functional. If None, PBE will be used as default
(Default value = None)
:type xc:
:param ca: The c/a ratio of hcp structures can be inputed here and if it is given,
it will be used by the elastic modulus routines (Default value = None)
:type ca:
:param **kwargs: Arbitrary other KGRN and KFCD input parameters can be given here
as keyword arguments. They will be passed down to the
self.emto.set_values() function
:type **kwargs: str,int,float,list(str),list(int),list(float)
:returns: None
:rtype: None
"""
if lat is None:
sys.exit('System.bulk(): \'lat\' has to be given!')
else:
self.lat = lat
if latname is None:
self.latname = self.lat
else:
self.latname = latname
if latpath is None:
self.latpath = "./"
else:
self.latpath = latpath
if emtopath is None:
self.emtopath = self.folder
else:
self.emtopath = emtopath
if atoms is None:
sys.exit('System.bulk(): \'atoms\' has to be given!')
else:
self.atoms = atoms
if concs is None:
# Assume equal concentrations for each element
self.concs = np.zeros(len(atoms))
self.concs[:] = 1.0 / float(len(atoms))
else:
self.concs = concs
if splts is None:
self.splts = np.zeros(len(atoms))
else:
self.splts = np.asarray(splts)
if sws is None:
self.sws = 0.0
#sys.exit('System.bulk(): \'sws\' has to be given!')
else:
self.sws = sws
if jobname is None:
self.jobname, self.fulljobname = self.create_jobname()
else:
self.jobname = jobname
self.fulljobname = self.create_jobname(jobname)
if ibz is None:
self.ibz = common.lat_to_ibz(self.lat)
else:
self.ibz = ibz
# Knowledge of the c/a lattice parameter for hcp systems
if ca is not None:
self.ca = ca
else:
self.ca = None
# Knowledge of the xc-functional we want to use
if xc is None:
self.xc = 'PBE'
else:
self.xc = xc
# Knowledge of the value of the bulk modulus, which
# is mainly needed in the elastic constant functions
self.bmod = bmod
# hcp requires that we "double" the atoms array and
# create a non-trivial iqs-array because hcp has a
# non-trivial two-atom basis.
if self.lat == 'hcp':
self.atoms = np.array([self.atoms, self.atoms]).flatten()
if concs is None:
self.concs = np.zeros(len(self.atoms))
self.concs[:] = 2.0 / float(len(self.atoms))
else:
self.concs = np.array([self.concs, self.concs]).flatten()
self.iqs = np.zeros(len(self.atoms), dtype='int32')
self.iqs[:len(self.iqs) // 2] = 1
self.iqs[len(self.iqs) // 2:] = 2
if splts is None:
self.splts = np.zeros(len(self.atoms))
else:
self.splts = np.array([self.splts, self.splts]).flatten()
self.itas = np.arange(1, len(self.atoms) / 2 + 1, dtype='int32')
self.itas = np.array([self.itas, self.itas]).flatten()
self.emto.set_values(jobname=self.fulljobname, sws=self.sws, atoms=self.atoms,
iqs=self.iqs, itas=self.itas, concs=self.concs, splts=self.splts,
ibz=self.ibz, latname=self.latname, latpath=self.latpath,
emtopath=self.emtopath, EMTOdir=self.EMTOdir, **kwargs)
# Special settings for the B2 structure. CPA currently not supported!!!
elif self.lat == 'B2':
self.iqs = np.array([1,2],dtype='int32')
self.concs = np.array([1.0,1.0])
self.its = np.array([1,2],dtype='int32')
self.itas = np.array([1,1],dtype='int32')
self.emto.set_values(jobname=self.fulljobname, sws=self.sws, atoms=self.atoms,
iqs=self.iqs, its=self.its, itas=self.itas, concs=self.concs,
splts=self.splts, ibz=self.ibz, latname=self.latname, latpath=self.latpath,
emtopath=self.emtopath, EMTOdir=self.EMTOdir, **kwargs)
else:
self.emto.set_values(jobname=self.fulljobname, sws=self.sws, atoms=self.atoms,
concs=self.concs, splts=self.splts, ibz=self.ibz,
latname=self.latname, latpath=self.latpath, emtopath=self.emtopath,
EMTOdir=self.EMTOdir, **kwargs)
# Make sure that structure files also receive the new input options,
# such as slurm_options.
self.lattice.set_values(**kwargs)
return
[docs] def bulk_new(self, jobname=None, lat=None, atoms=None, concs=None, splts=None, sws=None,
latname=None, latpath=None, emtopath=None, ibz=None, bmod=None, xc=None, ca=None,
iqs=None, its=None, itas=None,
**kwargs):
"""Initializes the basic parameters for bulk systems.
!!!A NEW VERSION OF THE OLD "bulk" ROUTINE!!!
Basic information concerning the system,
such as the types of atoms and the crystal structure should be given to this function and
it should be called right after the class instance has been created.
:param jobname: Name of the system (Default value = None)
:type jobname:
:param lat: The type of lattice structure (Default value = None)
:type lat:
:param atoms: List of atoms in the system (Default value = None)
:type atoms:
:param concs: List of concentrations of the elements in the
'atoms' list. This information is only used in CPA
calculations (Default value = None)
:type concs:
:param splts: List of initial magnetic moments of the elements in the
'atoms' list (Default value = None)
:type splts:
:param sws: The Wigner-Seitz radius of the system (Default value = None)
:type sws: float
:param latname: The 'jobname' of the BMDL, KSTR and SHAPE output files. These
structure output files have to be located in the 'latpath'
directory and they have to be named jobname.extention
(Default value = None)
:type latname:
:param latpath: The absolute path to the folder where the 'bmdl', 'kstr' and 'shape'
folders are located, which in turn contain the output files of the
structure calculation (Default value = None)
:type latpath:
:param emtopath: The absolute path to the folder where the EMTO installation is
located (Default value = None)
:type emtopath:
:param ibz: The code number indicating the Bravais lattice that the crystal
structure of the system has. For a list of possible values, please consult the
EMTO manual (Default value = None)
:type ibz:
:param bmod: The bulk modulus can be inputed here and if it is given,
it will be used by the elastic modulus routines (Default value = None)
:type bmod:
:param xc: The choice of the xc-functional. If None, PBE will be used as default
(Default value = None)
:type xc:
:param ca: The c/a ratio of hcp structures can be inputed here and if it is given,
it will be used by the elastic modulus routines (Default value = None)
:type ca:
:param **kwargs: Arbitrary other KGRN and KFCD input parameters can be given here
as keyword arguments. They will be passed down to the
self.emto.set_values() function
:type **kwargs: str,int,float,list(str),list(int),list(float)
:returns: None
:rtype: None
"""
if lat is None:
sys.exit('System.bulk(): \'lat\' has to be given!')
else:
self.lat = lat
if latname is None:
self.latname = self.lat
else:
self.latname = latname
if latpath is None:
self.latpath = "./"
else:
self.latpath = latpath
if emtopath is None:
self.emtopath = self.folder
else:
self.emtopath = emtopath
if atoms is None:
sys.exit('System.init_bulk(): \'atoms\' has to be given!')
else:
self.atoms = atoms
if concs is None:
# Assume equal concentrations for each element
self.concs = np.zeros(len(atoms))
self.concs[:] = 1.0 / float(len(atoms))
else:
self.concs = concs
if splts is None:
self.splts = np.zeros(len(atoms))
else:
self.splts = np.asarray(splts)
if sws is None:
self.sws = 0.0
#sys.exit('System.bulk(): \'sws\' has to be given!')
else:
self.sws = sws
if jobname is None:
self.jobname, self.fulljobname = self.create_jobname()
else:
self.jobname = jobname
self.fulljobname = self.create_jobname(jobname)
if ibz is None:
self.ibz = common.lat_to_ibz(self.lat)
else:
self.ibz = ibz
# Knowledge of the c/a lattice parameter for hcp systems
if ca is not None:
self.ca = ca
else:
self.ca = None
# Knowledge of the xc-functional we want to use
if xc is None:
self.xc = 'PBE'
else:
self.xc = xc
# Knowledge of the value of the bulk modulus, which
# is mainly needed in the elastic constant functions
self.bmod = bmod
# Construct input parameter arrays.
self.iqs = iqs
self.its = its
self.itas = itas
self.emto.set_values(jobname=self.fulljobname, sws=self.sws, atoms=self.atoms,
iqs=self.iqs, its=self.its, itas=self.itas, concs=self.concs,
splts=self.splts, ibz=self.ibz, latname=self.latname, latpath=self.latpath,
emtopath=self.emtopath, EMTOdir=self.EMTOdir, **kwargs)
return
[docs] def lattice_constants_analyze(self, sws=None, ca=None,prn=True,debug=False,method='morse',return_error=False):
"""Analyzes the output files generated using the
lattice_constants_batch_generate function.
The results are printed on screen.
:param sws: List of WS-radii (Default value = None)
:type sws: list(float)
:param ca: List hpc c/a ratios (Default value = None)
:type ca: list(float)
:param prn: True if results should be printed on screen, False if not (Default value = True)
:type prn: boolean
:returns: Equilibrium WS-radius, c/a (only hcp), bulk modulus, energy,
R (only hcp) and cs (only hcp)
:rtype: float, float (only hcp), float, float, float (only hcp), float (only hcp)
"""
from pyemto.EOS.EOS import EOS
eos = EOS(name=self.jobname, xc=self.xc, method=method, units='bohr')
if prn:
print('')
print('*****lattice_constants_analyze*****')
print('')
if sws is None:
sys.exit('System.lattice_constants_analyze(): An array of' +
' WS-radii \'sws\' has to be given!')
else:
self.lc_analyze_sws_range = np.asarray(sws)
if ca is None:
self.lc_analyze_ca_range = self.ca_range_default
else:
self.lc_analyze_ca_range = np.asarray(ca)
if self.lat == 'bcc' or self.lat == 'fcc' or self.lat == 'trig' or self.lat == 'stric':
energies = []
swses = []
for j in range(len(self.lc_analyze_sws_range)):
self.sws = self.lc_analyze_sws_range[j]
job = self.create_jobname(self.jobname)
en = self.get_energy(job, folder=self.folder, func=self.xc)
if isinstance(en, type(None)):
print('System.lattice_constants_analyze(): Warning:' +
' No output energy found for {0}'.format(job))
else:
energies.append(en)
swses.append(self.lc_analyze_sws_range[j])
if prn:
self.print_sws_ens(
'lattice_constants_analyze(cubic)', swses, energies)
sws0, e0, B0, grun, R_squared = eos.fit(swses, energies)
# These functions create files on disk about the data to be fitted
# as well as the results of the fit.
# eos.prepareData()
#sws0,e0,B,grun = eos.fit2file()
if prn:
print('lattice_constants_analyze(cubic):')
print('sws0 = {0:13.6f}'.format(sws0))
print('B0 = {0:13.6f}'.format(B0))
print('E0 = {0:13.6f}'.format(e0))
print('')
if return_error:
return sws0, B0, e0, grun, R_squared
else:
return sws0, B0, e0, grun
if self.lat == 'hcp':
# Fit an n'th order polynomial to the energy vs. c/a data.
ca_fit_order = 2
#latnames = ['hcp_ca1', 'hcp_ca2', 'hcp_ca3', 'hcp_ca4', 'hcp_ca5', 'hcp_ca6',
# 'hcp_ca7', 'hcp_ca8', 'hcp_ca9', 'hcp_ca10', 'hcp_ca11', 'hcp_ca12']
caname = 'hcp_ca'
# For the energies a 2D-array [i,j], where i = c/a axis and j = sws
# axis
energies = np.zeros(
(len(self.lc_analyze_ca_range), len(self.lc_analyze_sws_range)))
energies0 = [] # c/a optimized energy for a given WS-radius
swses = [] # List of WS-radii
cas0 = [] # Energetically optimized c/a's for a given WS-radius
# First collect all the output energies into the 2D array.
for i in range(len(self.lc_analyze_ca_range)):
for j in range(len(self.lc_analyze_sws_range)):
self.sws = self.lc_analyze_sws_range[j]
job = self.create_jobname(self.jobname + "_" + caname +str(i+1))
en = self.get_energy(job, folder=self.folder, func=self.xc)
if isinstance(en, type(None)):
print('System.lattice_constants_analyze(): Warning:' +
' No output energy found for {0}'.format(job))
else:
energies[i, j] = en
if debug:
print('Energy matrix (y axis = c/a axis, x axis = sws axis):'+"\n")
formaatti = " "
for i in range(len(self.lc_analyze_sws_range)):
formaatti = formaatti + "{0}{1}{2}{3} ".format("{",i,":8.6f","}")
print(formaatti.format(*self.lc_analyze_sws_range))
formaatti = ""
for i in range(len(self.lc_analyze_sws_range)+1):
formaatti = formaatti + "{0}{1}{2}{3} ".format("{",i,":13.6f","}")
for i in range(len(energies[:,0])):
print(formaatti.format(self.lc_analyze_ca_range[i],*energies[i,:]))
# Now we can start processing the 2D energy array.
# There might be some calculations that didn't converge.
# For those the energy will be zero in the 2D array so
# we have to leave those points out.
for j in range(len(self.lc_analyze_sws_range)):
good_energies = []
good_cas = []
for i in range(len(self.lc_analyze_ca_range)):
if energies[i, j] < -1.0:
good_energies.append(energies[i, j])
good_cas.append(self.lc_analyze_ca_range[i])
if len(good_energies) >= 3:
ca0, en0 = eos.ca_fit(good_cas, good_energies, ca_fit_order,
debug=debug,title='sws{0}'.format(j+1))
cas0.append(ca0)
energies0.append(en0)
swses.append(self.lc_analyze_sws_range[j])
if prn:
self.print_sws_ens_hcp(
'lattice_constants_analyze(hcp)', swses, energies0, cas0)
print('#'*80)
print('# Ground state EOS fit:'+' '*56+'#')
print('#'*80)
sws0, e0, B0, grun, R_squared = eos.fit(swses, energies0)
print('*'*80)
print('*'*80)
print('*'*80+'\n')
# These functions create files on disk about the data to be fitted
# as well as the results of the fit.
# eos.prepareData()
#sws0,e0,B0,grun = eos.fit2file()
# Now that we have the ground state WS-radius sws0 we can use the cas0 array
# to compute the corresponding ground state ca0
ca0_vs_sws = np.polyfit(swses, cas0, 2)
ca0_vs_sws = np.poly1d(ca0_vs_sws)
c_over_a0 = ca0_vs_sws(sws0)
# In order to calculate R = (c33 - c11 - c12 + c13) / cs
# we use Eq. (6.56) p. 108 in Vitos' book.
dca0_vs_dsws = np.polyder(ca0_vs_sws)
dca0dsws0 = dca0_vs_dsws(sws0)
R0 = -sws0 / 3.0 / c_over_a0 * dca0dsws0
# In order to calculate cs = c11 + c12 + 2c33 - 4c13
# we use Eq. (6.68) p. 109 in Vitos' book.
# We need the total energies as a function of c/a
# where sws is fixed to sws0. These can be obtained
# by doing a Morse fit for each c/a and evaluating
# the fitting function at sws0.
energies_cs = []
# There might be some calculations that didn't converge.
# For those the energy will be zero in the 2D array so
# we have to leave those points out.
for i in range(len(self.lc_analyze_ca_range)):
good_energies = []
good_swses = []
for j in range(len(self.lc_analyze_sws_range)):
if energies[i, j] < -1.0:
good_energies.append(energies[i, j])
good_swses.append(self.lc_analyze_sws_range[j])
print('#'*80)
print('# c/a{0} EOS fit:'.format(i+1)+' '*64+'#')
print('#'*80)
# _tmp variables are just dummies, we only want to
# update the EOS parameters of the "eos" instance.
sws_tmp, e_tmp, B_tmp, grun_tmp, R_squared_tmp = eos.fit(
good_swses, good_energies)
print('*'*80)
print('*'*80)
print('*'*80+'\n')
#e_cs = eos.fit_eval(sws0)
#e_cs = eos.fit_eval(sws_tmp)
e_cs = e_tmp
energies_cs.append(e_cs)
e_vs_ca_at_sws0 = np.polyfit(
self.lc_analyze_ca_range, energies_cs, 2)
e_vs_ca_at_sws0 = np.poly1d(e_vs_ca_at_sws0)
d2e_vs_dca2_at_sws0 = np.polyder(e_vs_ca_at_sws0, 2)
d2edca02sws0 = d2e_vs_dca2_at_sws0(c_over_a0)
vol0 = 4.0 / 3.0 * np.pi * sws0**3
cs0 = 9.0 / 2.0 * c_over_a0**2 / vol0 * \
d2edca02sws0 * self.RyBohr3_to_GPa
if prn:
print('hcp_lattice_constants_analyze(hcp):')
print('sws0 = {0:13.6f}'.format(sws0))
print('c/a0 = {0:13.6f}'.format(c_over_a0))
print('B0 = {0:13.6f}'.format(B0))
print('E0 = {0:13.6f}'.format(e0))
print('R = {0:13.6f}'.format(R0))
print('cs = {0:13.6f}'.format(cs0))
print('')
if return_error:
return sws0, c_over_a0, B0, e0, R0, cs0, grun, R_squared
else:
return sws0, c_over_a0, B0, e0, R0, cs0, grun
[docs] def lattice_constants_batch_generate(self, sws=None, ca=None, auto_ca=False):
"""Generates input files and writes them to disk.
Based on the input *sws* and *ca* lists jobnames are created and
then corresponding input files are generated and written
on disk. List of jobnames are returned.
:param sws: List of WS-radii (Default value = None)
:type sws: list(float)
:param ca: List of hcp c/a ratios (Default value = None)
:type ca: list(float)
:returns: List of jobnames
:rtype: list(str)
"""
if sws is None:
sys.exit('System.lattice_constants_batch(): An array of' +
' WS-radii \'sws\' has to be given!')
else:
self.lc_batch_sws_range = np.asarray(sws)
if ca is None:
self.lc_batch_ca_range = self.ca_range_default
else:
self.lc_batch_ca_range = np.asarray(ca)
jobnames = []
if self.lat == 'bcc' or self.lat == 'fcc':
for j in range(len(self.lc_batch_sws_range)):
self.sws = self.lc_batch_sws_range[j]
job = self.create_jobname(self.jobname)
jobnames.append(job)
self.emto.set_values(sws=self.sws, jobname=job)
common.check_folders(
self.folder, self.folder + "/kgrn", self.folder + "/kgrn/tmp")
common.check_folders(self.folder + "/kfcd")
common.check_folders(self.folder + "/fit")
self.emto.kgrn.write_input_file(folder=self.folder)
self.emto.kfcd.write_input_file(folder=self.folder)
self.emto.batch.write_input_file(folder=self.folder)
elif self.lat == 'hcp' and auto_ca is True:
for i in range(len(self.lc_batch_ca_range)):
for j in range(len(self.lc_batch_sws_range)):
caname = "hcp_ca"+str(i+1)
self.sws = self.lc_batch_sws_range[j]
job = self.create_jobname(self.jobname + "_" + caname)
jobnames.append(job)
self.emto.set_values(sws=self.sws, jobname=job)
self.emto.set_values(latname=caname)
common.check_folders(
self.folder, self.folder + "/kgrn", self.folder + "/kgrn/tmp")
common.check_folders(self.folder + "/kfcd")
common.check_folders(self.folder + "/fit")
self.emto.kgrn.write_input_file(folder=self.folder)
self.emto.kfcd.write_input_file(folder=self.folder)
self.emto.batch.write_input_file(folder=self.folder)
else:
for j in range(len(self.lc_batch_sws_range)):
self.sws = self.lc_batch_sws_range[j]
job = self.create_jobname(self.jobname)
jobnames.append(job)
self.emto.set_values(sws=self.sws, jobname=job)
common.check_folders(
self.folder, self.folder + "/kgrn", self.folder + "/kgrn/tmp")
common.check_folders(self.folder + "/kfcd")
common.check_folders(self.folder + "/fit")
self.emto.kgrn.write_input_file(folder=self.folder)
self.emto.kfcd.write_input_file(folder=self.folder)
self.emto.batch.write_input_file(folder=self.folder)
return jobnames
[docs] def lattice_constants_batch_calculate(self, sws=None, ca=None):
"""Calculates the ground state WS-radius using the parallelism of
the batch system by submitting one job for each entry in the *sws* list.
This is a combination of the batch_generate and batch_analyze functions.
At the end results are printed on screen.
:param sws: List of WS-radii (Default value = None)
:type sws: list(float)
:param ca: hpc c/a ratio (Default value = None)
:type ca: float
:returns: WS-radius, c/a (only hcp), bulk modulus, energy,
R (only hcp), cs (only hcp)
:rtype: float, float (only hcp), float, float,
float (only hcp), float (only hcp)
"""
if sws is None:
sys.exit('System.lattice_constants_calculate(): An array of' +
' WS-radii \'sws\' has to be given!')
else:
self.lc_batch_sws_range = np.asarray(sws)
if ca is None:
self.lc_batch_ca_range = self.ca_range_default
else:
self.lc_batch_ca_range = np.asarray(ca)
# Create input files
jobnames = self.lattice_constants_batch_generate(
sws=self.lc_batch_sws_range, ca=self.lc_batch_ca_range)
# Submit calculations to the batch queue system
jobIDs = self.submit_jobs(jobnames, folder=self.folder)
# Wait until all the jobs have finished
self.wait_for_jobs(jobIDs)
# Now we can analyze the results
if self.lat == 'bcc' or self.lat == 'fcc':
sws0, B0, e0, grun = self.lattice_constants_analyze(
sws=sws, ca=self.lc_batch_ca_range)
return sws0, B0, e0, grun
elif self.lat == 'hcp':
sws0, c_over_a0, B0, e0, R0, cs0,grun = self.lattice_constants_analyze(
sws=sws, ca=self.lc_batch_ca_range)
return sws0, c_over_a0, B0, e0, R0, cs0, grun
[docs] def lattice_constants_serial_calculate(self, sws=None, stype="simple", rerun=False, skip=False,
delta=0.01, refine=True):
"""Calculates the equilibrium lattice constants on one CPU **WITHOUT** the batch system.
Also the eq. bulk modulus, c/a ratio and energy are returned.
:param sws: Initial guess for the eq. WS-radius (Default value = None)
:type sws: float
:param stype: Type of the energy minimisation algorithm (Default value = "simple")
:type stype: str
:param rerun: True if (Default value = False)
:type rerun: boolean
:param skip: True if (Default value = False)
:type skip: boolean
:param delta: Step size (Default value = 0.01)
:type delta: float
:param refine: True if an optimized WS-radius vs. energy curve
should be computed, False if not (Default value = True)
:type refine: boolean
:returns: WS-radius, c/a, bulk modulus, energy
:rtype: float, float, float, float
"""
if sws is None:
sys.exit('System.lattice_constants_serial_calculate():' +
' Starting point \'sws\' has to be given!')
else:
self.lc_initial_sws = sws
self.lc_stype = stype
self.lc_rerun = rerun
self.lc_skip = skip
# Find initial volume and bulk module and refine those results.
if self.lat == 'bcc' or self.lat == 'fcc':
c_over_a = 0.0 # sc,bcc and fcc structures only need a
print('Running self.finc_lc()')
v, b, energy = self.find_lc(delta=delta, xc=self.xc)
print('System.find_lc(): sws,Bmod,energy = ', v, b, energy)
if refine:
print('Running self.refine_lc()')
v, b, energy = self.refine_lc(v, delta=delta, xc=self.xc)
print('lattice_constants_serial_calculate: sws,c/a,Bmod,energy = ',
v, c_over_a, b, energy)
elif self.lat == 'hcp':
print('Running self.find_lc_hcp()')
v, b, c_over_a, energy = self.find_lc_hcp(delta=delta, xc=self.xc)
print('System.find_lc_hcp(): sws,c/a,Bmod,energy = ',
v, c_over_a, b, energy)
if refine:
print('Running self.refine_lc_hcp()')
v, b, c_over_a, energy = self.refine_lc_hcp(
v, c_over_a, delta=delta, xc=self.xc)
print('lattice_constants_serial_calculate: sws,c/a,Bmod,energy = ',
v, c_over_a, b, energy)
print()
return v, c_over_a, b, energy
[docs] def elastic_constants_analyze(
self, sws=None, bmod=None, ca=None, R=None, cs=None, relax=True, debug=False):
"""Analyzes the output files generated using the
elastic_constants_batch_generate function.
The results are printed on screen.
:param sws: WS-radius (Default value = None)
:type sws: float
:param bmod: Bulk modulus (Default value = None)
:type bmod: float
:param ca: hpc c/a ratio (Default value = None)
:type ca: float
:param R: Dimensionless quantity of hcp systems (Default value = None)
:type R: float
:param cs: Second order c/a-derivative of the energy (Default value = None)
:type cs: float
:returns: None
:rtype: None
"""
from pyemto.EOS.EOS import EOS
eos = EOS(name=self.jobname, xc=self.xc, method='morse', units='bohr')
# Mission critical parameters
if sws is None:
sys.exit(
'System.elastic_constants_analyze(): \'sws\' has not been specified!')
elif sws is not None:
self.sws = sws
if bmod is None:
sys.exit(
'System.elastic_constants_analyze(): \'sws\' has not been specified!')
elif bmod is not None:
self.bmod = bmod
if ca is None and self.lat == 'hcp':
sys.exit(
'System.elastic_constants_analyze(): \'ca\' (c/a) has not been specified!')
else:
self.ec_analyze_ca = ca
if R is None and self.lat == 'hcp':
sys.exit(
'System.elastic_constants_analyze(): \'R\' has not been specified!')
else:
self.ec_analyze_R = R
if cs is None and self.lat == 'hcp':
sys.exit(
'System.elastic_constants_analyze(): \'cs\' has not been specified!')
else:
self.ec_analyze_cs = cs
deltas = self.elastic_constants_deltas
if self.lat == 'bcc' or self.lat == 'fcc' or self.lat == 'B2':
# Orthorhombic distortion for c' first
if self.lat == 'bcc':
jobname_dist = [
'_bcco0', '_bcco1', '_bcco2', '_bcco3', '_bcco4', '_bcco5']
latname_dist = [
'bcco0', 'bcco1', 'bcco2', 'bcco3', 'bcco4', 'bcco5']
elif self.lat == 'fcc':
jobname_dist = [
'_fcco0', '_fcco1', '_fcco2', '_fcco3', '_fcco4', '_fcco5']
latname_dist = [
'fcco0', 'fcco1', 'fcco2', 'fcco3', 'fcco4', 'fcco5']
elif self.lat == 'B2':
jobname_dist = [
'_B2o0', '_B2o1', '_B2o2', '_B2o3', '_B2o4', '_B2o5']
latname_dist = [
'B2o0', 'B2o1', 'B2o2', 'B2o3', 'B2o4', 'B2o5']
en_cprime = []
good_deltas_cprime = []
for i in range(len(jobname_dist)):
already = False
job = self.create_jobname(self.jobname + jobname_dist[i])
en = self.get_energy(job, folder=self.folder, func=self.xc)
if isinstance(en, type(None)):
print('System.cubic_elastic_constants_analyze(): Warning:' +
' No output energy found for {0}'.format(job))
else:
en_cprime.append(en)
good_deltas_cprime.append(deltas[i])
# Exit if we don't have enough points to do the fit
if len(en_cprime) < 3:
sys.exit(
'System.elastic_constants_analyze(): Not enough energy points to fit c\'!')
# Convert into a numpy array
en_cprime = np.asarray(en_cprime)
good_deltas_cprime = np.asarray(good_deltas_cprime)
# Next the monoclinic distortion for c44
if self.lat == 'bcc':
jobname_dist = [
'_bccm0', '_bccm1', '_bccm2', '_bccm3', '_bccm4', '_bccm5']
latname_dist = [
'bccm0', 'bccm1', 'bccm2', 'bccm3', 'bccm4', 'bccm5']
elif self.lat == 'fcc':
jobname_dist = [
'_fccm0', '_fccm1', '_fccm2', '_fccm3', '_fccm4', '_fccm5']
latname_dist = [
'fccm0', 'fccm1', 'fccm2', 'fccm3', 'fccm4', 'fccm5']
elif self.lat == 'B2':
jobname_dist = [
'_B2m0', '_B2m1', '_B2m2', '_B2m3', '_B2m4', '_B2m5']
latname_dist = [
'B2m0', 'B2m1', 'B2m2', 'B2m3', 'B2m4', 'B2m5']
en_c44 = []
good_deltas_c44 = []
for i in range(len(jobname_dist)):
already = False
job = self.create_jobname(self.jobname + jobname_dist[i])
en = self.get_energy(job, folder=self.folder, func=self.xc)
if isinstance(en, type(None)):
print('System.elastic_constants_analyze(): Warning:' +
' No output energy found for {0}'.format(job))
else:
en_c44.append(en)
good_deltas_c44.append(deltas[i])
# Exit if we don't have enough points to do the fit
if len(en_c44) < 3:
sys.exit(
'System.elastic_constants_analyze(): Not enough energy points to fit c44!')
# Convert into a numpy array
en_c44 = np.asarray(en_c44)
good_deltas_c44 = np.asarray(good_deltas_c44)
# All calculations have been done, now it's time to fit the results
popt_cprime, cprime_rsq = eos.distortion_fit(
good_deltas_cprime, en_cprime,title='cprime')
popt_c44, c44_rsq = eos.distortion_fit(
good_deltas_c44, en_c44,title='c44')
volume = 4.0 / 3.0 * np.pi * self.sws**3
cprime = popt_cprime[0] / 2.0 / volume * self.RyBohr3_to_GPa
c44 = popt_c44[0] / 2.0 / volume * self.RyBohr3_to_GPa
c11 = self.bmod + 4.0 / 3.0 * cprime
c12 = self.bmod - 2.0 / 3.0 * cprime
# Polycrystalline elastic constants
#
# B = bulk modulus
# G = shear modulus
# E = Young modulus
# v = Poisson ratio
# Voigt average
BV = (c11 + 2 * c12) / 3.0
#BV = self.bmod
GV = (c11 - c12 + 3 * c44) / 5.0
EV = 9 * BV * GV / (3 * BV + GV)
vV = (3 * BV - 2 * GV) / (6 * BV + 2 * GV)
# Reuss average
BR = BV
#BR = self.bmod
GR = 5 * (c11 - c12) * c44 / (4 * c44 + 3 * (c11 - c12))
ER = 9 * BR * GR / (3 * BR + GR)
vR = (3 * BR - 2 * GR) / (6 * BR + 2 * GR)
# Hill average
BH = (BV + BR) / 2.0
#BH = self.bmod
GH = (GV + GR) / 2.0
EH = 9 * BH * GH / (3 * BH + GH)
vH = (3 * BH - 2 * GH) / (6 * BH + 2 * GH)
# Elastic anisotropy
AVR = (GV - GR) / (GV + GR)
print("")
print('***cubic_elastic_constants***')
print("")
print(self.jobname)
print("")
print('sws(bohr) = {0:7.3f}'.format(self.sws))
print('B(GPa) = {0:6.2f}'.format(self.bmod))
print('c11(GPa) = {0:6.2f}'.format(c11))
print('c12(GPa) = {0:6.2f}'.format(c12))
print('c\'(GPa) = {0:6.2f}'.format(cprime))
print('c44(GPa) = {0:6.2f}'.format(c44))
print('R-squared(c\') = {0:8.6f}'.format(cprime_rsq))
print('R-squared(c44) = {0:8.6f}'.format(c44_rsq))
print("")
print('Voigt average:')
print("")
print('BV(GPa) = {0:6.2f}'.format(BV))
print('GV(GPa) = {0:6.2f}'.format(GV))
print('EV(GPa) = {0:6.2f}'.format(EV))
print('vV(GPa) = {0:6.2f}'.format(vV))
print("")
print('Reuss average:')
print("")
print('BR(GPa) = {0:6.2f}'.format(BR))
print('GR(GPa) = {0:6.2f}'.format(GR))
print('ER(GPa) = {0:6.2f}'.format(ER))
print('vR(GPa) = {0:6.2f}'.format(vR))
print("")
print('Hill average:')
print("")
print('BH(GPa) = {0:6.2f}'.format(BH))
print('GH(GPa) = {0:6.2f}'.format(GH))
print('EH(GPa) = {0:6.2f}'.format(EH))
print('vH(GPa) = {0:6.2f}'.format(vH))
print("")
print('Elastic anisotropy:')
print("")
print('AVR(GPa) = {0:6.2f}'.format(AVR))
return
elif self.lat == 'hcp':
# Orthorhombic distortion for c66 first
if relax:
jobname_dist = []
latname_dist = []
for i in range(self.elastic_constants_points):
tmp_array1 = []
tmp_array2 = []
for j in range(self.hcpo_relax_points):
tmp_str1 = 'hcpo{0}_r{1}'.format(i,j)
tmp_str2 = '_hcpo{0}_r{1}'.format(i,j)
tmp_array1.append(tmp_str1)
tmp_array2.append(tmp_str2)
# We don't need to relax the first structure
# because it's undistorted
if i == 0 and j == 0:
break
latname_dist.append(tmp_array1)
jobname_dist.append(tmp_array2)
else:
jobname_dist = ['_hcpo0_ur',
'_hcpo1_ur',
'_hcpo2_ur',
'_hcpo3_ur',
'_hcpo4_ur',
'_hcpo5_ur']
latname_dist = ['hcpo0_ur',
'hcpo1_ur',
'hcpo2_ur',
'hcpo3_ur',
'hcpo4_ur',
'hcpo5_ur']
en_c66 = []
good_deltas_c66 = []
if relax:
# First we have to determine how relaxation affects
# the energy. The first structure corresponds to the
# undistorted lattice, which means we don't need to
# relax it:
job = self.create_jobname(self.jobname + jobname_dist[0][0])
en_orig = self.get_energy(job, folder=self.folder, func=self.xc)
if isinstance(en_orig, type(None)):
print('System.elastic_constants_analyze(): Warning:' +
' No output energy found for {0}'.format(job))
else:
en_c66.append(en_orig)
good_deltas_c66.append(deltas[0])
# For all the rest of the structures we compute the
# relaxed ground state energy by fitting a polynomial
# to the atomic pos. vs. energy data:
for i in range(1,self.elastic_constants_points):
ens_tmp = []
xind_tmp = []
for j in range(len(jobname_dist[i])):
already = False
job = self.create_jobname(self.jobname + jobname_dist[i][j])
en = self.get_energy(job, folder=self.folder, func=self.xc)
if isinstance(en, type(None)):
print('System.elastic_constants_analyze(): Warning:' +
' No output energy found for {0}'.format(job))
else:
ens_tmp.append(en)
xind_tmp.append(j)
# Now we should find the relaxed energy by fitting
# a polynomial to the energy vs. atom coordinates data.
# Leave out the point if not enough calculations have
# converged.
if len(ens_tmp) >= 3:
relax_xmin, relax_emin = eos.relax_fit(xind_tmp,ens_tmp,2,debug=debug,
title='c66: hcpo{0}'.format(i))
en_c66.append(relax_emin)
# TEST: Use unrelaxed energies
#en_c66.append(ens_tmp[0])
good_deltas_c66.append(deltas[i])
else:
for i in range(self.elastic_constants_points):
already = False
job = self.create_jobname(self.jobname + jobname_dist[i])
en = self.get_energy(job, folder=self.folder, func=self.xc)
if isinstance(en, type(None)):
print('System.elastic_constants_analyze(): Warning:' +
' No output energy found for {0}'.format(job))
else:
en_c66.append(en)
good_deltas_c66.append(deltas[i])
# Exit if we don't have enough points to do the fit
if len(en_c66) < 3:
sys.exit('System.elastic_constants_analyze():' +
' Not enough energy points to fit hcp c66!')
# Convert into a numpy array
en_c66 = np.asarray(en_c66)
good_deltas_c66 = np.asarray(good_deltas_c66)
# Next the monoclinic distortion for c44
# For the monoclinic distortion relaxation effects
# are negligible:
relax = False
if relax:
jobname_dist = []
latname_dist = []
for i in range(self.elastic_constants_points):
tmp_array1 = []
tmp_array2 = []
for j in range(self.hcpo_relax_points):
tmp_str1 = 'hcpm{0}_r{1}'.format(i,j)
tmp_str2 = '_hcpm{0}_r{1}'.format(i,j)
tmp_array1.append(tmp_str1)
tmp_array2.append(tmp_str2)
latname_dist.append(tmp_array1)
jobname_dist.append(tmp_array2)
else:
jobname_dist = ['_hcpm0_ur',
'_hcpm1_ur',
'_hcpm2_ur',
'_hcpm3_ur',
'_hcpm4_ur',
'_hcpm5_ur']
latname_dist = ['hcpm0_ur',
'hcpm1_ur',
'hcpm2_ur',
'hcpm3_ur',
'hcpm4_ur',
'hcpm5_ur']
en_c44 = []
good_deltas_c44 = []
if relax:
# First we have to determine how relaxation affects
# the energy.
for i in range(self.elastic_constants_points):
ens_tmp = []
xind_tmp = []
for j in range(self.hcpo_relax_points):
already = False
job = self.create_jobname(self.jobname + jobname_dist[i][j])
en = self.get_energy(job, folder=self.folder, func=self.xc)
if isinstance(en, type(None)):
print('System.elastic_constants_analyze(): Warning:' +
' No output energy found for {0}'.format(job))
else:
ens_tmp.append(en)
xind_tmp.append(j)
# Now we should find the relaxed energy by fitting
# a polynomial to the energy vs. atom coordinates data.
# Leave out the point if not enough calculations have
# converged.
if len(ens_tmp) >= 3:
relax_xmin, relax_emin = eos.relax_fit(xind_tmp,ens_tmp,3,debug=debug,
title='c44: hcpm{0}'.format(i))
en_c44.append(relax_emix)
good_deltas_c44.append(deltas[i])
else:
for i in range(self.elastic_constants_points):
already = False
job = self.create_jobname(self.jobname + jobname_dist[i])
en = self.get_energy(job, folder=self.folder, func=self.xc)
if isinstance(en, type(None)):
print('System.elastic_constants_analyze(): Warning:' +
' No output energy found for {0}'.format(job))
else:
en_c44.append(en)
good_deltas_c44.append(deltas[i])
# Exit if we don't have enough points to do the fit
if len(en_c44) < 3:
sys.exit('System.elastic_constants_analyze():' +
' Not enough energy points to fit hcp c44!')
# Convert into a numpy array
en_c44 = np.asarray(en_c44)
good_deltas_c44 = np.asarray(good_deltas_c44)
# All calculations have been done, now it's time to fit the results
popt_c66, c66_rsq = eos.distortion_fit(good_deltas_c66, en_c66, title='c66')
popt_c44, c44_rsq = eos.distortion_fit(good_deltas_c44, en_c44, title='c44')
volume = 4.0 / 3.0 * np.pi * self.sws**3
c66 = popt_c66[0] / 2.0 / volume * self.RyBohr3_to_GPa
c44 = popt_c44[0] / 2.0 / volume * self.RyBohr3_to_GPa
c11 = self.bmod + c66 + self.ec_analyze_cs * \
(2 * self.ec_analyze_R - 1)**2 / 18.0
c12 = self.bmod - c66 + self.ec_analyze_cs * \
(2 * self.ec_analyze_R - 1)**2 / 18.0
c13 = self.bmod + 1.0 / 9.0 * self.ec_analyze_cs * \
(2 * self.ec_analyze_R**2 + self.ec_analyze_R - 1)
c33 = self.bmod + 2.0 / 9.0 * \
self.ec_analyze_cs * (self.ec_analyze_R + 1)**2
c2 = c33 * (c11 + c12) - 2.0 * c13**2
# Polycrystalline elastic constants
#
# B = bulk modulus
# G = shear modulus
# E = Young modulus
# v = Poisson ratio
# Voigt average
BV = (2 * c11 + 2 * c12 + 4 * c13 + c33) / 9.0
GV = (12 * c44 + 12 * c66 + self.ec_analyze_cs) / 30.0
EV = 9 * BV * GV / (3 * BV + GV)
vV = (3 * BV - 2 * GV) / (6 * BV + 2 * GV)
# Reuss average
BR = self.bmod
GR = 5.0 / 2.0 * (c44 * c66 * c2) / \
((c44 + c66) * c2 + 3.0 * BV * c44 * c66)
ER = 9 * BR * GR / (3 * BR + GR)
vR = (3 * BR - 2 * GR) / (6 * BR + 2 * GR)
# Hill average
BH = (BV + BR) / 2.0
#BH = self.bmod
GH = (GV + GR) / 2.0
EH = 9 * BH * GH / (3 * BH + GH)
vH = (3 * BH - 2 * GH) / (6 * BH + 2 * GH)
# Elastic anisotropy
AVR = (GV - GR) / (GV + GR)
print("")
print('***hcp_elastic_constants***')
print("")
print(self.jobname)
print("")
print('sws(bohr) = {0:7.3f}'.format(self.sws))
print('B(GPa) = {0:6.2f}'.format(self.bmod))
print('c11(GPa) = {0:6.2f}'.format(c11))
print('c12(GPa) = {0:6.2f}'.format(c12))
print('c13(GPa) = {0:6.2f}'.format(c13))
print('c33(GPa) = {0:6.2f}'.format(c33))
print('c44(GPa) = {0:6.2f}'.format(c44))
print('c66(GPa) = {0:6.2f}'.format(c66))
print('R-squared(c44) = {0:8.6f}'.format(c44_rsq))
print('R-squared(c66) = {0:8.6f}'.format(c66_rsq))
print("")
print('Voigt average:')
print("")
print('BV(GPa) = {0:6.2f}'.format(BV))
print('GV(GPa) = {0:6.2f}'.format(GV))
print('EV(GPa) = {0:6.2f}'.format(EV))
print('vV(GPa) = {0:6.2f}'.format(vV))
print("")
print('Reuss average:')
print("")
print('BR(GPa) = {0:6.2f}'.format(BR))
print('GR(GPa) = {0:6.2f}'.format(GR))
print('ER(GPa) = {0:6.2f}'.format(ER))
print('vR(GPa) = {0:6.2f}'.format(vR))
print("")
print('Hill average:')
print("")
print('BH(GPa) = {0:6.2f}'.format(BH))
print('GH(GPa) = {0:6.2f}'.format(GH))
print('EH(GPa) = {0:6.2f}'.format(EH))
print('vH(GPa) = {0:6.2f}'.format(vH))
print("")
print('Elastic anisotropy:')
print("")
print('AVR(GPa) = {0:6.2f}'.format(AVR))
return
[docs] def elastic_constants_batch_generate(self, sws=None, ca=None, relax=True):
"""Generates all the necessary input files based on the class data.
:param sws: WS-radius (Default value = None)
:type sws: float
:param ca: hpc c/a ratio (Default value = None)
:type ca: float
:returns: List of jobnames
:rtype: list(str)
"""
# Mission critical parameters
if sws is None and self.sws is None:
sys.exit(
'System.elastic_constants_batch_generate(): \'sws\' has not been specified!')
elif sws is not None:
self.sws = sws
if ca is None and self.lat == 'hcp':
sys.exit('System.elastic_constants_batch_generate():' +
' \'ca\' (c/a) has not been specified!')
elif ca is not None:
self.ca = ca
jobnames = []
if self.lat == 'B2':
# Orthorhombic distortion input files for c' first
jobname_dist = ['_B2o0','_B2o1','_B2o2','_B2o3','_B2o4','_B2o5']
latname_dist = ['B2o0','B2o1','B2o2','B2o3','B2o4','B2o5']
self.emto.set_values(ibz=8,nkx=17,nky=17,nkz=17) # High-quality
#self.emto.set_values(ibz=10,nkx=15,nky=15,nkz=15) # Normal quality
for i in range(len(jobname_dist)):
job = self.create_jobname(self.jobname + jobname_dist[i])
jobnames.append(job)
self.emto.set_values(
sws=self.sws, jobname=job, latname=latname_dist[i])
common.check_folders(
self.folder, self.folder + "/kgrn", self.folder + "/kgrn/tmp")
common.check_folders(self.folder + "/kfcd")
common.check_folders(self.folder + "/fit")
self.emto.kgrn.write_input_file(folder=self.folder)
self.emto.kfcd.write_input_file(folder=self.folder)
self.emto.batch.write_input_file(folder=self.folder)
# Next produce the input files of monoclinic distortion for c44
jobname_dist = ['_B2m0','_B2m1','_B2m2','_B2m3','_B2m4','_B2m5']
latname_dist = ['B2m0','B2m1','B2m2','B2m3','B2m4','B2m5']
self.emto.set_values(ibz=9,nkx=17,nky=17,nkz=17) # High-quality
#self.emto.set_values(ibz=11,nkx=15,nky=15,nkz=21) # Normal quality
for i in range(len(jobname_dist)):
job = self.create_jobname(self.jobname + jobname_dist[i])
jobnames.append(job)
self.emto.set_values(
sws=self.sws, jobname=job, latname=latname_dist[i])
self.emto.kgrn.write_input_file(folder=self.folder)
self.emto.kfcd.write_input_file(folder=self.folder)
self.emto.batch.write_input_file(folder=self.folder)
return jobnames
if self.lat == 'bcc' or self.lat == 'fcc':
# Orthorhombic distortion input files for c' first
if self.lat == 'bcc':
jobname_dist = ['_bcco0','_bcco1','_bcco2','_bcco3','_bcco4','_bcco5']
latname_dist = ['bcco0','bcco1','bcco2','bcco3','bcco4','bcco5']
self.emto.set_values(ibz=10,nkx=31,nky=31,nkz=31) # High-quality
#self.emto.set_values(ibz=10,nkx=15,nky=15,nkz=15) # Normal quality
elif self.lat == 'fcc':
jobname_dist = ['_fcco0','_fcco1','_fcco2','_fcco3','_fcco4','_fcco5']
latname_dist = ['fcco0','fcco1','fcco2','fcco3','fcco4','fcco5']
#self.emto.set_values(ibz=11,nkx=41,nky=41,nkz=41) # Super-high-quality
self.emto.set_values(ibz=11,nkx=31,nky=31,nkz=31) # High-quality
#self.emto.set_values(ibz=11,nkx=17,nky=17,nkz=17) # Normal-quality
for i in range(len(jobname_dist)):
job = self.create_jobname(self.jobname + jobname_dist[i])
jobnames.append(job)
self.emto.set_values(
sws=self.sws, jobname=job, latname=latname_dist[i])
common.check_folders(
self.folder, self.folder + "/kgrn", self.folder + "/kgrn/tmp")
common.check_folders(self.folder + "/kfcd")
common.check_folders(self.folder + "/fit")
self.emto.kgrn.write_input_file(folder=self.folder)
self.emto.kfcd.write_input_file(folder=self.folder)
self.emto.batch.write_input_file(folder=self.folder)
# Next produce the input files of monoclinic distortion for c44
if self.lat == 'bcc':
jobname_dist = ['_bccm0','_bccm1','_bccm2','_bccm3','_bccm4','_bccm5']
latname_dist = ['bccm0','bccm1','bccm2','bccm3','bccm4','bccm5']
self.emto.set_values(ibz=11,nkx=27,nky=27,nkz=37) # High-quality
#self.emto.set_values(ibz=11,nkx=15,nky=15,nkz=21) # Normal quality
elif self.lat == 'fcc':
jobname_dist = ['_fccm0','_fccm1','_fccm2','_fccm3','_fccm4','_fccm5']
latname_dist = ['fccm0','fccm1','fccm2','fccm3','fccm4','fccm5']
#self.emto.set_values(ibz=10,nkx=37,nky=53,nkz=37) # Super-high-quality
self.emto.set_values(ibz=10,nkx=27,nky=37,nkz=27) # High-quality
#self.emto.set_values(ibz=10,nkx=15,nky=21,nkz=15) # Normal quality
for i in range(len(jobname_dist)):
job = self.create_jobname(self.jobname + jobname_dist[i])
jobnames.append(job)
self.emto.set_values(
sws=self.sws, jobname=job, latname=latname_dist[i])
self.emto.kgrn.write_input_file(folder=self.folder)
self.emto.kfcd.write_input_file(folder=self.folder)
self.emto.batch.write_input_file(folder=self.folder)
return jobnames
elif self.lat == 'hcp':
# Orthorhombic distortion input files for c66 first
# With hcp the structure depends on the c/a ratio. Therefore we also have
# to generate the corresponding structure files.
if relax:
jobname_dist = []
latname_dist = []
for i in range(self.elastic_constants_points):
tmp_array1 = []
tmp_array2 = []
for j in range(self.hcpo_relax_points):
tmp_str1 = 'hcpo{0}_r{1}'.format(i,j)
tmp_str2 = '_hcpo{0}_r{1}'.format(i,j)
tmp_array1.append(tmp_str1)
tmp_array2.append(tmp_str2)
# We don't need to relax the first structure
# because it's undistorted
if i == 0 and j == 0:
break
latname_dist.append(tmp_array1)
jobname_dist.append(tmp_array2)
else:
jobname_dist = ['_hcpo0_ur', '_hcpo1_ur',
'_hcpo2_ur', '_hcpo3_ur', '_hcpo4_ur', '_hcpo5_ur']
latname_dist = ['hcpo0_ur', 'hcpo1_ur',
'hcpo2_ur', 'hcpo3_ur', 'hcpo4_ur', 'hcpo5_ur']
# Check whether Two-center Taylor expansion is on/off
if self.emto.kgrn.expan == 'M':
kappaw = self.kappaw_default
self.lattice.set_values(kappaw=kappaw)
common.check_folders(self.folder)
common.check_folders(self.folder + '/bmdl')
common.check_folders(self.folder + '/kstr')
common.check_folders(self.folder + '/shape')
if relax:
for i in range(self.elastic_constants_points):
for j in range(len(jobname_dist[i])):
if i == 0:
do_i_relax = False
else:
do_i_relax = True
self.lattice.distortion(lat='hcp', dist='ortho', ca=self.ca, index=i,
deltas=self.elastic_constants_deltas,
relax=do_i_relax,relax_index=j)
self.lattice.set_values(jobname_lat=latname_dist[i][j],latpath=self.folder)
self.lattice.bmdl.write_input_file(folder=self.folder)
self.lattice.kstr.write_input_file(folder=self.folder)
self.lattice.shape.write_input_file(folder=self.folder)
self.lattice.batch.write_input_file(folder=self.folder)
else:
for i in range(self.elastic_constants_points):
self.lattice.distortion(lat='hcp', dist='ortho', ca=self.ca, index=i,
deltas=self.elastic_constants_deltas,
relax=False)
self.lattice.set_values(jobname_lat=latname_dist[i],latpath=self.folder)
self.lattice.bmdl.write_input_file(folder=self.folder)
self.lattice.kstr.write_input_file(folder=self.folder)
self.lattice.shape.write_input_file(folder=self.folder)
self.lattice.batch.write_input_file(folder=self.folder)
#self.emto.set_values(ibz=9, nkx=25, nky=25, nkz=21)
self.emto.set_values(ibz=9, nkx=31, nky=31, nkz=25)
if relax:
for i in range(self.elastic_constants_points):
for j in range(len(jobname_dist[i])):
job = self.create_jobname(self.jobname + jobname_dist[i][j])
jobnames.append(job)
self.emto.set_values(sws=self.sws, jobname=job, latname=latname_dist[i][j],
latpath=self.folder)
common.check_folders(self.folder, self.folder + "/kgrn")
common.check_folders(self.folder + "/kgrn/tmp")
common.check_folders(self.folder + "/kfcd")
common.check_folders(self.folder + "/fit")
self.emto.kgrn.write_input_file(folder=self.folder)
self.emto.kfcd.write_input_file(folder=self.folder)
self.emto.batch.write_input_file(folder=self.folder)
else:
for i in range(self.elastic_constants_points):
job = self.create_jobname(self.jobname + jobname_dist[i])
jobnames.append(job)
self.emto.set_values(sws=self.sws, jobname=job, latname=latname_dist[i],
latpath=self.folder)
common.check_folders(self.folder, self.folder + "/kgrn")
common.check_folders(self.folder + "/kgrn/tmp")
common.check_folders(self.folder + "/kfcd")
common.check_folders(self.folder + "/fit")
self.emto.kgrn.write_input_file(folder=self.folder)
self.emto.kfcd.write_input_file(folder=self.folder)
self.emto.batch.write_input_file(folder=self.folder)
# Monoclinic distortion input files for c44 next
# For the monoclinic distortion relaxation effects are
# negligibly small; therefore
relax = False
# With hcp the structure depends on the c/a ratio. Therefore we also have
# to generate the corresponding structure files.
if relax:
jobname_dist = []
latname_dist = []
for i in range(self.elastic_constants_points):
tmp_array1 = []
tmp_array2 = []
for j in range(self.hcpm_relax_points):
tmp_str1 = 'hcpm{0}_r{1}'.format(i,j)
tmp_str2 = '_hcpm{0}_r{1}'.format(i,j)
tmp_array1.append(tmp_str1)
tmp_array2.append(tmp_str2)
# We don't need to relax the first structure
# because it's undistorted
if i == 0 and j == 0:
break
latname_dist.append(tmp_array1)
jobname_dist.append(tmp_array2)
else:
jobname_dist = ['_hcpm0_ur', '_hcpm1_ur',
'_hcpm2_ur', '_hcpm3_ur', '_hcpm4_ur', '_hcpm5_ur']
latname_dist = ['hcpm0_ur', 'hcpm1_ur',
'hcpm2_ur', 'hcpm3_ur', 'hcpm4_ur', 'hcpm5_ur']
# Check whether Two-center Taylor expansion is on/off
if self.emto.kgrn.expan == 'M':
kappaw = self.kappaw_default
self.lattice.set_values(kappaw=kappaw)
common.check_folders(self.folder)
common.check_folders(self.folder + '/bmdl')
common.check_folders(self.folder + '/kstr')
common.check_folders(self.folder + '/shape')
if relax:
for i in range(self.elastic_constants_points):
for j in range(len(jobname_dist[i])):
if i == 0:
do_i_relax = False
else:
do_i_relax = True
self.lattice.distortion(lat='hcp', dist='mono', ca=self.ca, index=i,
deltas=self.elastic_constants_deltas,
relax=do_i_relax,relax_index=j)
self.lattice.set_values(jobname_lat=latname_dist[i][j],latpath=self.folder)
self.lattice.bmdl.write_input_file(folder=self.folder)
self.lattice.kstr.write_input_file(folder=self.folder)
self.lattice.shape.write_input_file(folder=self.folder)
self.lattice.batch.write_input_file(folder=self.folder)
else:
for i in range(self.elastic_constants_points):
self.lattice.distortion(lat='hcp', dist='mono', ca=self.ca, index=i,
deltas=self.elastic_constants_deltas,
relax=False)
self.lattice.set_values(jobname_lat=latname_dist[i],latpath=self.folder)
self.lattice.bmdl.write_input_file(folder=self.folder)
self.lattice.kstr.write_input_file(folder=self.folder)
self.lattice.shape.write_input_file(folder=self.folder)
self.lattice.batch.write_input_file(folder=self.folder)
# Two-atom basis
#self.emto.set_values(ibz=13, nkx=41, nky=41, nkz=41)
if relax:
for i in range(self.elastic_constants_points):
for j in range(len(jobname_dist[i])):
job = self.create_jobname(self.jobname + jobname_dist[i][j])
jobnames.append(job)
self.emto.set_values(sws=self.sws, jobname=job, latname=latname_dist[i][j],
latpath=self.folder)
common.check_folders(self.folder, self.folder + "/kgrn")
common.check_folders(self.folder + "/kgrn/tmp")
common.check_folders(self.folder + "/kfcd")
common.check_folders(self.folder + "/fit")
self.emto.kgrn.write_input_file(folder=self.folder)
self.emto.kfcd.write_input_file(folder=self.folder)
self.emto.batch.write_input_file(folder=self.folder)
else:
# Four-atom basis
#self.emto.set_values(ibz=12, nkx=18, nky=18, nkz=12)
self.emto.set_values(ibz=12, nkx=22, nky=22, nkz=14)
###################################################################
# Atconf related arrays need to be modified because we now have #
# a four atom basis. #
###################################################################
self.atoms = np.array([self.atoms, self.atoms]).flatten()
self.concs = np.array([self.concs, self.concs]).flatten()
self.iqs = np.zeros(len(self.atoms), dtype='int32')
len_div = len(self.iqs) // 4
for i in range(4):
self.iqs[i * len_div:(i + 1) * len_div] = i + 1
self.splts = np.array([self.splts, self.splts]).flatten()
self.itas = np.array([self.itas, self.itas]).flatten()
self.emto.set_values(atoms=self.atoms, iqs=self.iqs, itas=self.itas,
concs=self.concs, splts=self.splts)
####################################################################
for i in range(self.elastic_constants_points):
job = self.create_jobname(self.jobname + jobname_dist[i])
jobnames.append(job)
self.emto.set_values(sws=self.sws, jobname=job, latname=latname_dist[i],
latpath=self.folder)
common.check_folders(self.folder, self.folder + "/kgrn")
common.check_folders(self.folder + "/kgrn/tmp")
common.check_folders(self.folder + "/kfcd")
common.check_folders(self.folder + "/fit")
self.emto.kgrn.write_input_file(folder=self.folder)
self.emto.kfcd.write_input_file(folder=self.folder)
self.emto.batch.write_input_file(folder=self.folder)
# These following lines are for the four-atom basis with
# simple monoclinic lattice.
"""
for i in range(len(self.elastic_constants_deltas)):
self.lattice.distortion(lat='hcp', dist='mono', ca=self.ca, index=i,
deltas=self.elastic_constants_deltas)
self.lattice.set_values(
jobname_lat=latname_dist[i], latpath=self.folder)
self.lattice.bmdl.write_input_file(folder=self.folder)
self.lattice.kstr.write_input_file(folder=self.folder)
self.lattice.shape.write_input_file(folder=self.folder)
self.lattice.batch.write_input_file(folder=self.folder)
self.emto.set_values(ibz=12, nkx=30, nky=20, nkz=20)
################################################################
# Atconf related arrays need to be modified because we now have
# a four atom basis.
################################################################
self.atoms = np.array([self.atoms, self.atoms]).flatten()
self.concs = np.array([self.concs, self.concs]).flatten()
self.iqs = np.zeros(len(self.atoms), dtype='int32')
len_div = len(self.iqs) // 4
for i in range(4):
self.iqs[i * len_div:(i + 1) * len_div] = i + 1
self.splts = np.array([self.splts, self.splts]).flatten()
self.itas = np.array([self.itas, self.itas]).flatten()
self.emto.set_values(atoms=self.atoms, iqs=self.iqs, itas=self.itas,
concs=self.concs, splts=self.splts)
for i in range(len(jobname_dist)):
job = self.create_jobname(self.jobname + jobname_dist[i])
jobnames.append(job)
self.emto.set_values(sws=self.sws, jobname=job, latname=latname_dist[i],
latpath=self.folder)
common.check_folders(
self.folder, self.folder + "/kgrn", self.folder + "/kgrn/tmp")
common.check_folders(self.folder + "/kfcd")
common.check_folders(self.folder + "/fit")
self.emto.kgrn.write_input_file(folder=self.folder)
self.emto.kfcd.write_input_file(folder=self.folder)
self.emto.batch.write_input_file(folder=self.folder)
"""
return jobnames
[docs] def elastic_constants_batch_calculate(
self, sws=None, bmod=None, ca=None, R=None, cs=None):
"""Calculates the elastic constants of a system using the parallelism
of the batch system.
This is a combination of the batch_generate and batch_analyze functions.
:param sws: WS-radius (Default value = None)
:type sws: float
:param bmod: Bulk modulus (Default value = None)
:type bmod: float
:param ca: hpc c/a ratio (Default value = None)
:type ca: float
:param R: The dimensionless quantity of hcp systems (Default value = None)
:type R: float
:param cs: Second order c/a-derivative of the energy (Default value = None)
:type cs: float
:returns: None
:rtype: None
"""
import time
# Mission critical parameters
if sws is None and self.sws is None:
sys.exit(
'System.elastic_constants_batch_calculate(): \'sws\' has not been specified!')
elif sws is not None:
self.sws = sws
if bmod is None and self.bmod is None:
sys.exit(
'System.elastic_constants_batch_calculate(): \'bmod\' has not been specified!')
elif bmod is not None:
self.bmod = bmod
if ca is None and self.lat == 'hcp':
sys.exit('System.elastic_constants_batch_calculate():' +
' \'ca\' (c/a) has not been specified!')
else:
self.ec_batch_calculate_ca = ca
if R is None and self.lat == 'hcp':
sys.exit(
'System.elastic_constants_batch_calculate(): \'R\' has not been specified!')
else:
self.ec_batch_calculate_R = R
if cs is None and self.lat == 'hcp':
sys.exit(
'System.elastic_constants_batch_calculate(): \'cs\' has not been specified!')
else:
self.ec_batch_calculate_cs = cs
# Generate input files
if self.lat == 'bcc' or self.lat == 'fcc':
jobnames = self.elastic_constants_batch_generate(sws=self.sws)
# Submit calculation to the batch system
jobIDs = self.submit_jobs(jobnames, folder=self.folder)
# Wait until all the jobs have finished
self.wait_for_jobs(jobIDs)
# Now we can analyze the results
self.elastic_constants_analyze(sws=self.sws, bmod=self.bmod)
return
elif self.lat == 'hcp':
jobnames = self.elastic_constants_batch_generate(
sws=self.sws, ca=self.ec_batch_calculate_ca)
jobnames_lat = ['hcpo0_ca0',
'hcpo1_ca0',
'hcpo2_ca0',
'hcpo3_ca0',
'hcpo4_ca0',
'hcpo5_ca0',
'hcpm0_ca0',
'hcpm1_ca0',
'hcpm2_ca0',
'hcpm3_ca0',
'hcpm4_ca0',
'hcpm5_ca0']
# First submit and run the lattice calculations
jobIDs_lat = self.submit_jobs(jobnames_lat, folder=self.folder)
# Wait until all the jobs have finished
self.wait_for_jobs(jobIDs_lat)
# Structure calculations have finished, now submit
# KGRN and KFCD calculation.
jobIDs = self.submit_jobs(jobnames, folder=self.folder)
# Wait until all the jobs have finished
self.wait_for_jobs(jobIDs)
# Now we can analyze the results
self.elastic_constants_analyze(
sws=self.sws, bmod=self.bmod, ca=self.ec_batch_calculate_ca,
R=self.ec_batch_calculate_R, cs=self.ec_batch_calculate_cs)
return
[docs] def elastic_constants_serial_calculate(self, sws=None, bmod=None, ca=None, R=None, cs=None):
"""Calculates elastic constants on one CPU without using the batch system.
At the end the results are printed on screen.
:param sws: WS-radius (Default value = None)
:type sws: float
:param bmod: Bulk modulus (Default value = None)
:type bmod: float
:param ca: hpc c/a ratio (Default value = None)
:type ca: float
:returns: None
:rtype: None
"""
from pyemto.EOS.EOS import EOS
eos = EOS(name=self.jobname, xc=self.xc, method='morse', units='bohr')
# Mission critical parameters
if sws is None and self.sws is None:
sys.exit(
'System.elastic_constants_serial_calculate(): \'sws\' has not been specified!')
elif sws is not None:
self.sws = sws
if bmod is None and self.bmod is None:
sys.exit('System.elastic_constants_serial_calculate():' +
' \'bmod\' has not been specified!')
elif bmod is not None:
self.bmod = bmod
if ca is None and self.lat == 'hcp':
sys.exit('System.elastic_constants_serial_calculate():' +
' \'ca\' (c/a) has not been specified!')
elif ca is not None:
self.ca = ca
if R is None and self.lat == 'hcp':
sys.exit(
'System.elastic_constants_analyze(): \'R\' has not been specified!')
else:
self.ec_analyze_R = R
if cs is None and self.lat == 'hcp':
sys.exit(
'System.elastic_constants_analyze(): \'cs\' has not been specified!')
else:
self.ec_analyze_cs = cs
deltas = self.elastic_constants_deltas
if self.lat == 'bcc' or self.lat == 'fcc':
# Orthorhombic distortion for c' first
if self.lat == 'bcc':
jobname_dist = [
'_bcco0', '_bcco1', '_bcco2', '_bcco3', '_bcco4', '_bcco5']
latname_dist = [
'bcco0', 'bcco1', 'bcco2', 'bcco3', 'bcco4', 'bcco5']
self.emto.set_values(ibz=10, nkx=27, nky=27, nkz=27)
elif self.lat == 'fcc':
jobname_dist = [
'_fcco0', '_fcco1', '_fcco2', '_fcco3', '_fcco4', '_fcco5']
latname_dist = [
'fcco0', 'fcco1', 'fcco2', 'fcco3', 'fcco4', 'fcco5']
self.emto.set_values(ibz=11, nkx=27, nky=27, nkz=27)
en_cprime = []
for i in range(len(jobname_dist)):
already = False
job = self.create_jobname(self.jobname + jobname_dist[i])
self.emto.set_values(
sws=self.sws, jobname=job, latname=latname_dist[i])
# check if calculations are already done
already = self.check_conv(job, folder=self.folder)
if all(already):
conv = (True, True)
else:
if already[0] and not already[1]:
conv = self.runemto(
jobname=job, folder=self.folder, onlyKFCD=True)
else:
conv = self.runemto(jobname=job, folder=self.folder)
# KGRN has crashed, find out why
if conv[0] == False:
self.which_error(job, folder=self.folder)
quit()
en = self.get_energy(job, folder=self.folder, func=self.xc)
en_cprime.append(en)
# Convert energies into a numpy array
en_cprime = np.asarray(en_cprime)
# Next calculate the monoclinic distortion for c44
if self.lat == 'bcc':
jobname_dist = [
'_bccm0', '_bccm1', '_bccm2', '_bccm3', '_bccm4', '_bccm5']
latname_dist = [
'bccm0', 'bccm1', 'bccm2', 'bccm3', 'bccm4', 'bccm5']
self.emto.set_values(ibz=11, nkx=27, nky=27, nkz=37)
elif self.lat == 'fcc':
jobname_dist = [
'_fccm0', '_fccm1', '_fccm2', '_fccm3', '_fccm4', '_fccm5']
latname_dist = [
'fccm0', 'fccm1', 'fccm2', 'fccm3', 'fccm4', 'fccm5']
self.emto.set_values(ibz=10, nkx=27, nky=37, nkz=27)
en_c44 = []
for i in range(len(jobname_dist)):
already = False
job = self.create_jobname(self.jobname + jobname_dist[i])
self.emto.set_values(
sws=self.sws, jobname=job, latname=latname_dist[i])
# check if calculations are already done
already = self.check_conv(job, folder=self.folder)
if all(already):
conv = (True, True)
else:
if already[0] and not already[1]:
conv = self.runemto(
jobname=job, folder=self.folder, onlyKFCD=True)
else:
conv = self.runemto(jobname=job, folder=self.folder)
# KGRN has crashed, find out why
if conv[0] == False:
self.which_error(job, folder=self.folder)
quit()
en = self.get_energy(job, folder=self.folder, func=self.xc)
en_c44.append(en)
# Convert energies into a numpy array
en_c44 = np.asarray(en_c44)
# All calculations have been done, now it's time to fit the results
popt_cprime, cprime_rsq = eos.distortion_fit(deltas, en_cprime)
popt_c44, c44_rsq = eos.distortion_fit(deltas, en_c44)
volume = 4.0 / 3.0 * np.pi * self.sws**3
cprime = popt_cprime[0] / 2.0 / volume * self.RyBohr3_to_GPa
c44 = popt_c44[0] / 2.0 / volume * self.RyBohr3_to_GPa
c11 = self.bmod + 4.0 / 3.0 * cprime
c12 = self.bmod - 2.0 / 3.0 * cprime
# Polycrystalline elastic constants
#
# B = bulk modulus
# G = shear modulus
# E = Young modulus
# v = Poisson ratio
# Voigt average
BV = (c11 + 2 * c12) / 3.0
#BV = self.bmod
GV = (c11 - c12 + 3 * c44) / 5.0
EV = 9 * BV * GV / (3 * BV + GV)
vV = (3 * BV - 2 * GV) / (6 * BV + 2 * GV)
# Reuss average
BR = BV
#BR = self.bmod
GR = 5 * (c11 - c12) * c44 / (4 * c44 + 3 * (c11 - c12))
ER = 9 * BR * GR / (3 * BR + GR)
vR = (3 * BR - 2 * GR) / (6 * BR + 2 * GR)
# Hill average
BH = (BV + BR) / 2.0
#BH = self.bmod
GH = (GV + GR) / 2.0
EH = 9 * BH * GH / (3 * BH + GH)
vH = (3 * BH - 2 * GH) / (6 * BH + 2 * GH)
# Elastic anisotropy
AVR = (GV - GR) / (GV + GR)
print("")
print('***cubic_elastic_constants***')
print("")
print(self.jobname)
print("")
print('c11(GPa) = {0:6.2f}'.format(c11))
print('c12(GPa) = {0:6.2f}'.format(c12))
print(
'c44(GPa) = {0:6.2f}, R-squared = {1:8.6f}'.format(c44, c44_rsq))
print(
'c\' (GPa) = {0:6.2f}, R-squared = {1:8.6f}'.format(cprime, cprime_rsq))
print('B (GPa) = {0:6.2f}'.format(self.bmod))
print("")
print('Voigt average:')
print("")
print('BV(GPa) = {0:6.2f}'.format(BV))
print('GV(GPa) = {0:6.2f}'.format(GV))
print('EV(GPa) = {0:6.2f}'.format(EV))
print('vV(GPa) = {0:6.2f}'.format(vV))
print("")
print('Reuss average:')
print("")
print('BR(GPa) = {0:6.2f}'.format(BR))
print('GR(GPa) = {0:6.2f}'.format(GR))
print('ER(GPa) = {0:6.2f}'.format(ER))
print('vR(GPa) = {0:6.2f}'.format(vR))
print("")
print('Hill average:')
print("")
print('BH(GPa) = {0:6.2f}'.format(BH))
print('GH(GPa) = {0:6.2f}'.format(GH))
print('EH(GPa) = {0:6.2f}'.format(EH))
print('vH(GPa) = {0:6.2f}'.format(vH))
print("")
print('Elastic anisotropy:')
print("")
print('AVR(GPa) = {0:6.2f}'.format(AVR))
return
elif self.lat == 'hcp':
# Check whether Two-center Taylor expansion is on/off
if self.emto.kgrn.expan == 'M':
kappaw = self.kappaw_default
self.lattice.set_values(kappaw=kappaw)
# Check whether we need to create some folders
# for the structure output files.
common.check_folders(self.folder)
common.check_folders(self.folder + '/bmdl')
common.check_folders(self.folder + '/kstr')
common.check_folders(self.folder + '/shape')
# Orthorhombic distortion for c66 first
jobname_dist = ['_hcpo0_ca0', '_hcpo1_ca0',
'_hcpo2_ca0', '_hcpo3_ca0', '_hcpo4_ca0', '_hcpo5_ca0']
latname_dist = ['hcpo0_ca0', 'hcpo1_ca0',
'hcpo2_ca0', 'hcpo3_ca0', 'hcpo4_ca0', 'hcpo5_ca0']
self.emto.set_values(ibz=9, nkx=31, nky=19, nkz=19)
en_c66 = []
for i in range(len(jobname_dist)):
# With hcp the structure depends on the c/a ratio. Therefore we also have
# to generate the corresponding structure files.
self.lattice.distortion(lat='hcp', dist='ortho', ca=self.ca, index=i,
deltas=self.elastic_constants_deltas)
self.lattice.set_values(
jobname_lat=latname_dist[i], latpath=self.folder)
self.runlattice(jobname=latname_dist[i], folder=self.folder)
already = False
job = self.create_jobname(self.jobname + jobname_dist[i])
self.emto.set_values(
sws=self.sws, jobname=job, latname=latname_dist[i])
# check if calculations are already done
already = self.check_conv(job, folder=self.folder)
if all(already):
conv = (True, True)
else:
if already[0] and not already[1]:
conv = self.runemto(
jobname=job, folder=self.folder, onlyKFCD=True)
else:
conv = self.runemto(jobname=job, folder=self.folder)
# KGRN has crashed, find out why
if conv[0] == False:
self.which_error(job, folder=self.folder)
quit()
en = self.get_energy(job, folder=self.folder, func=self.xc)
en_c66.append(en)
# Convert energies into a numpy array
en_c66 = np.asarray(en_cprime)
# Monoclinic distortion for c44 next
jobname_dist = ['_hcpm0_ca0',
'_hcpm1_ca0',
'_hcpm2_ca0',
'_hcpm3_ca0',
'_hcpm4_ca0',
'_hcpm5_ca0']
latname_dist = ['hcpm0_ca0',
'hcpm1_ca0',
'hcpm2_ca0',
'hcpm3_ca0',
'hcpm4_ca0',
'hcpm5_ca0']
self.emto.set_values(ibz=12, nkx=30, nky=20, nkz=20)
en_c44 = []
for i in range(len(jobname_dist)):
# With hcp the structure depends on the c/a ratio. Therefore we also have
# to generate the corresponding structure files.
self.lattice.distortion(lat='hcp', dist='mono', ca=self.ca, index=i,
deltas=self.elastic_constants_deltas)
self.lattice.set_values(
jobname_lat=latname_dist[i], latpath=self.folder)
self.runlattice(jobname=latname_dist[i], folder=self.folder)
###############################################################
# Atconf related arrays need to be modified because we now have
# a four atom basis.
###############################################################
self.atoms = np.array([self.atoms, self.atoms]).flatten()
self.concs = np.array([self.concs, self.concs]).flatten()
self.iqs = np.zeros(len(self.atoms), dtype='int32')
len_div = len(self.iqs) // 4
for i in range(4):
self.iqs[i * len_div:(i + 1) * len_div] = i + 1
self.splts = np.array([self.splts, self.splts]).flatten()
self.itas = np.array([self.itas, self.itas]).flatten()
self.emto.set_values(atoms=self.atoms, iqs=self.iqs, itas=self.itas,
concs=self.concs, splts=self.splts)
already = False
job = self.create_jobname(self.jobname + jobname_dist[i])
self.emto.set_values(
sws=self.sws, jobname=job, latname=latname_dist[i])
# check if calculations are already done
already = self.check_conv(job, folder=self.folder)
if all(already):
conv = (True, True)
else:
if already[0] and not already[1]:
conv = self.runemto(
jobname=job, folder=self.folder, onlyKFCD=True)
else:
conv = self.runemto(jobname=job, folder=self.folder)
# KGRN has crashed, find out why
if conv[0] == False:
self.which_error(job, folder=self.folder)
quit()
en = self.get_energy(job, folder=self.folder, func=self.xc)
en_c44.append(en)
# Convert energies into a numpy array
en_c44 = np.asarray(en_c44)
# All calculations have been done, now it's time to fit the results
popt_c66, c66_rsq = eos.distortion_fit(deltas, en_c66)
popt_c44, c44_rsq = eos.distortion_fit(deltas, en_c44)
volume = 4.0 / 3.0 * np.pi * self.sws**3
c66 = popt_c66[0] / 2.0 / volume * self.RyBohr3_to_GPa
c44 = popt_c44[0] / 2.0 / volume * self.RyBohr3_to_GPa
c11 = self.bmod + c66 + self.ec_analyze_cs * \
(2 * self.ec_analyze_R - 1)**2 / 18.0
c12 = self.bmod - c66 + self.ec_analyze_cs * \
(2 * self.ec_analyze_R - 1)**2 / 18.0
c13 = self.bmod + 1.0 / 9.0 * self.ec_analyze_cs * (
2 * self.ec_analyze_R**2 + self.ec_analyze_R - 1)
c33 = self.bmod + 2.0 / 9.0 * \
self.ec_analyze_cs * (self.ec_analyze_R + 1)**2
c2 = c33 * (c11 + c12) - 2.0 * c13**2
# Polycrystalline elastic constants
#
# B = bulk modulus
# G = shear modulus
# E = Young modulus
# v = Poisson ratio
# Voigt average
BV = (2 * c11 + 2 * c12 + 4 * c13 + c33) / 9.0
GV = (12 * c44 + 12 * c66 + self.ec_analyze_cs) / 30.0
EV = 9 * BV * GV / (3 * BV + GV)
vV = (3 * BV - 2 * GV) / (6 * BV + 2 * GV)
# Reuss average
BR = self.bmod
GR = 5.0 / 2.0 * (c44 * c66 * c2) / \
((c44 + c66) * c2 + 3.0 * BV * c44 * c66)
ER = 9 * BR * GR / (3 * BR + GR)
vR = (3 * BR - 2 * GR) / (6 * BR + 2 * GR)
# Hill average
BH = (BV + BR) / 2.0
#BH = self.bmod
GH = (GV + GR) / 2.0
EH = 9 * BH * GH / (3 * BH + GH)
vH = (3 * BH - 2 * GH) / (6 * BH + 2 * GH)
# Elastic anisotropy
AVR = (GV - GR) / (GV + GR)
print("")
print('***hcp_elastic_constants***')
print("")
print(self.jobname)
print("")
print('c11(GPa) = {0:6.2f}'.format(c11))
print('c12(GPa) = {0:6.2f}'.format(c12))
print('c13(GPa) = {0:6.2f}'.format(c13))
print('c33(GPa) = {0:6.2f}'.format(c33))
print(
'c44(GPa) = {0:6.2f}, R-squared = {1:8.6f}'.format(c44, c44_rsq))
print(
'c66(GPa) = {0:6.2f}, R-squared = {1:8.6f}'.format(c66, c66_rsq))
print('B (GPa) = {0:6.2f}'.format(self.bmod))
print("")
print('Voigt average:')
print("")
print('BV(GPa) = {0:6.2f}'.format(BV))
print('GV(GPa) = {0:6.2f}'.format(GV))
print('EV(GPa) = {0:6.2f}'.format(EV))
print('vV(GPa) = {0:6.2f}'.format(vV))
print("")
print('Reuss average:')
print("")
print('BR(GPa) = {0:6.2f}'.format(BR))
print('GR(GPa) = {0:6.2f}'.format(GR))
print('ER(GPa) = {0:6.2f}'.format(ER))
print('vR(GPa) = {0:6.2f}'.format(vR))
print("")
print('Hill average:')
print("")
print('BH(GPa) = {0:6.2f}'.format(BH))
print('GH(GPa) = {0:6.2f}'.format(GH))
print('EH(GPa) = {0:6.2f}'.format(EH))
print('vH(GPa) = {0:6.2f}'.format(vH))
print("")
print('Elastic anisotropy:')
print("")
print('AVR(GPa) = {0:6.2f}'.format(AVR))
return
##########################################################################
# #
# Internal routines start here #
# #
##########################################################################
[docs] def find_lc(self, delta=0.01, prn=True, xc='PBE'):
"""Computes initial estimates for the ground state quantities for cubic systems.
:param delta: Step size for the volume vs. energy array (Default value = 0.01)
:type delta: float
:param prn: True if results should be printed, False if not (Default value = True)
:type prn: boolean
:param xc: Choice of the xc-functional (Default value = 'PBE')
:type xc: str
:returns: WS-radius, bulk modulus, c/a and energy
:rtype: float, float, float, float
"""
from pyemto.EOS.EOS import EOS
eos = EOS(name=self.jobname, xc=xc, method='morse', units='bohr')
if prn:
print('')
print('*****find_lc*****')
print('')
#enough = False
energies = []
swses = []
# Compute first two points
# Start at the initial sws
already = False
self.sws = self.lc_initial_sws
job = self.create_jobname(self.jobname)
self.emto.set_values(sws=self.sws, jobname=job)
# check if calculations are already done
already = self.check_conv(job, folder=self.folder)
if self.lc_skip or (all(already) and not self.lc_rerun):
conv = (True, True)
else:
if already[0] and not already[1]:
conv = self.runemto(
jobname=job, folder=self.folder, onlyKFCD=True)
else:
conv = self.runemto(jobname=job, folder=self.folder)
# KGRN has crashed, find out why
if conv[0] == False:
self.which_error(job, folder=self.folder)
quit()
en = self.get_energy(job, folder=self.folder, func=xc)
energies.append(en)
swses.append(self.sws)
# Compute second point at sws+delta
already = False
self.sws = self.lc_initial_sws + delta
job = self.create_jobname(self.jobname)
self.emto.set_values(sws=self.sws, jobname=job)
# check if calculations are already done
already = self.check_conv(job, folder=self.folder)
#print('already = ',already)
if self.lc_skip or (all(already) and not self.lc_rerun):
conv = (True, True)
else:
if already[0] and not already[1]:
conv = self.runemto(
jobname=job, folder=self.folder, onlyKFCD=True)
else:
conv = self.runemto(jobname=job, folder=self.folder)
#print('conv = ',conv)
# KGRN has crashed, find out why
if conv[0] == False:
self.which_error(job, folder=self.folder)
quit()
en = self.get_energy(job, folder=self.folder, func=xc)
energies.append(en)
swses.append(self.sws)
# The initial 2 points are now ready and
# we use them to predict where to calculate next point
next_sws, minfound = self.predict_next_sws(swses, energies, delta)
# 2 points is not enought to calculate equation of state
enough = False
# Loop until we have enough points to calculate initial sws
iteration = 0
while not enough:
iteration += 1
# if 25 is not enough one should check initial sws
if iteration > 25:
print(
"SWS loop did not converge in {0} iterations!".format(iteration))
quit()
# Calculate next point
self.sws = next_sws
job = self.create_jobname(self.jobname)
self.emto.set_values(sws=self.sws, jobname=job)
# First check if calculations are already done
already = False
already = self.check_conv(job, folder=self.folder)
# Use existing calculations if available.
if self.lc_skip or (all(already) and not self.lc_rerun):
conv = (True, True)
else:
if already[0] and not already[1]:
conv = self.runemto(
jobname=job, folder=self.folder, onlyKFCD=True)
else:
conv = self.runemto(jobname=job, folder=self.folder)
# KGRN has crashed, find out why
if conv[0] == False:
self.which_error(job, folder=self.folder)
quit()
en = self.get_energy(job, folder=self.folder, func=xc)
energies.append(en)
swses.append(self.sws)
# Check do we have minumun and predict next sws
next_sws, minfound = self.predict_next_sws(swses, energies, delta)
# Check if we have mimimum and enough points
# 9 points should be enough for EOS
if minfound and len(swses) > 9:
enough = True
#print('debug find_lc: ',swses)
if prn:
self.print_sws_ens('find_lc', swses, energies)
sws0, e0, B, grun = eos.fit(swses, energies)
# These functions create files on disk about the data to be fitted
# as well as the results of the fit.
# eos.prepareData()
#sws0,e0,B,grun = eos.fit2file()
return sws0, B, e0
[docs] def refine_lc(self, sws, delta=0.01, prn=True, xc='PBE'):
"""Calculates a more accurate equilibrium volume for cubic systems.
:param sws: WS-radius
:type sws: float
:param delta: Step size for the volume vs. energy array (Default value = 0.01)
:type delta: float
:param prn: True if results should be printed, False if not (Default value = True)
:type prn: boolean
:param xc: Choice of the xc-functional (Default value = 'PBE')
:type xc: str
:returns: Ws-radius, bulk modulus, c/a and energy
:rtype: float, float, float, float
"""
from pyemto.EOS.EOS import EOS
eos = EOS(name=self.jobname, xc=xc, method='morse', units='bohr')
if prn:
print('')
print('*****refine_lc*****')
print('')
# make sws ranges around given sws
points = []
for i in range(-3, 7):
points.append(i * delta)
energies = []
swses = []
for p in points:
already = False
next_sws = sws + p
# Make inputs and job name
self.sws = next_sws
job = self.create_jobname(self.jobname)
self.emto.set_values(sws=self.sws, jobname=job)
# check if calculations are already done
already = self.check_conv(job, folder=self.folder)
if self.lc_skip or (all(already) and not self.lc_rerun):
conv = (True, True)
else:
if already[0] and not already[1]:
conv = self.runemto(
jobname=job, folder=self.folder, onlyKFCD=True)
else:
conv = self.runemto(jobname=job, folder=self.folder)
# KGRN has crashed, find out why
if conv[0] == False:
self.which_error(job, folder=self.folder)
quit()
en = self.get_energy(job, folder=self.folder, func=xc)
energies.append(en)
swses.append(self.sws)
if prn:
self.print_sws_ens('refine_lc', swses, energies)
sws0, e0, B, grun = eos.fit(swses, energies)
# These functions create files on disk about the data to be fitted
# as well as the results of the fit.
# eos.prepareData()
#sws0,e0,B,grun = eos.fit2file()
return sws0, B, e0
[docs] def find_lc_hcp(self, delta=0.01, prn=True, xc='PBE'):
"""Computes initial estimates for the ground state quantities for hcp systems.
:param delta: Step size for the volume vs. energy array (Default value = 0.01)
:type delta: float
:param prn: True if results should be printed, False if not (Default value = True)
:type prn: boolean
:param xc: Choice of the xc-functional (Default value = 'PBE')
:type xc: str
:returns: WS-radius, bulk modulus, c/a and energy
:rtype: float, float, float, float
"""
from pyemto.EOS.EOS import EOS
eos_hcp = EOS(name=self.jobname, xc=xc, method='morse', units='bohr')
if prn:
print('')
print('*****find_lc_hcp*****')
print('')
#enough = False
# Before anything else happens we generate the
# structure output files (BMDL, KSTR and SHAPE) for a reasonable c over
# a mesh
ca_mesh = self.ca_range_default
ca_mesh_len = len(ca_mesh)
# sws-optimization is done using this particular c/a value,
ca_sws_ind = 2
# all the other c/a's will use these same volumes.
# Fit an n'th order polynomial to the energy vs. c/a data.
ca_fit_order = 2
ca_latpaths = []
ca_prefixes = []
# A 2D-array [i,j], where i = c/a axis and j = sws axis
energies = np.zeros((ca_mesh_len, 50))
energies0 = [] # c/a optimized energy for a given WS-radius
swses = [] # List of WS-radii
cas0 = [] # Energetically optimized c/a's for a given WS-radius
for i in range(0, ca_mesh_len):
ca_prefixes.append("/ca_{0:8.6f}".format(ca_mesh[i]))
ca_latpaths.append(self.latpath + ca_prefixes[i])
# Check whether structure files already exists
# and if not run the calculations.
for i in range(0, ca_mesh_len):
self.lattice.set_values(ca=ca_mesh[i])
# self.lattice.bmdl.write_input_file(folder=ca_latpaths[i])
str_exists = self.check_str(self.latname, folder=ca_latpaths[i])
if str_exists == False:
self.runlattice(jobname=self.latname, folder=ca_latpaths[i])
# Compute first two points
# Start at the initial sws
already = False
self.sws = self.lc_initial_sws
job = self.create_jobname(self.jobname)
self.emto.set_values(sws=self.sws, jobname=job)
# check if calculations are already done
for i in range(ca_mesh_len):
# We have to remember to update the lattice path every
# time we change c/a
self.emto.set_values(latpath=ca_latpaths[i])
hcp_subfolder = self.folder + "/{0}".format(ca_prefixes[i])
already = self.check_conv(job, folder=hcp_subfolder)
if self.lc_skip or (all(already) and not self.lc_rerun):
conv = (True, True)
else:
if already[0] and not already[1]:
conv = self.runemto(
jobname=job, folder=hcp_subfolder, onlyKFCD=True)
else:
conv = self.runemto(jobname=job, folder=hcp_subfolder)
# KGRN has crashed, find out why
if conv[0] == False:
self.which_error(job, folder=hcp_subfolder)
quit()
en = self.get_energy(job, folder=hcp_subfolder, func=xc)
energies[i, 0] = en
swses.append(self.sws)
# Calculate energetically optimized c/a and the corresponding energy
# print(energies[:,:15])
ca0, en0 = eos_hcp.ca_fit(ca_mesh, energies[:, 0], ca_fit_order)
cas0.append(ca0)
energies0.append(en0)
# Compute second point at initial sws + delta
already = False
self.sws = self.lc_initial_sws + delta
job = self.create_jobname(self.jobname)
self.emto.set_values(sws=self.sws, jobname=job)
# check if calculations are already done
for i in range(ca_mesh_len):
# We have to remember to update the lattice path every
# time we change c/a
self.emto.set_values(latpath=ca_latpaths[i])
hcp_subfolder = self.folder + "/{0}".format(ca_prefixes[i])
already = self.check_conv(job, folder=hcp_subfolder)
if self.lc_skip or (all(already) and not self.lc_rerun):
conv = (True, True)
else:
if already[0] and not already[1]:
conv = self.runemto(
jobname=job, folder=hcp_subfolder, onlyKFCD=True)
else:
conv = self.runemto(jobname=job, folder=hcp_subfolder)
# KGRN has crashed, find out why
if conv[0] == False:
self.which_error(job, folder=hcp_subfolder)
quit()
en = self.get_energy(job, folder=hcp_subfolder, func=xc)
energies[i, 1] = en
swses.append(self.sws)
# Calculate energetically optimized c/a and the corresponding energy
# print(energies[:,:15])
ca0, en0 = eos_hcp.ca_fit(ca_mesh, energies[:, 1], ca_fit_order)
cas0.append(ca0)
energies0.append(en0)
# The initial 2 points are now ready and
# we use them to predict where to calculate next point
next_sws, minfound = self.predict_next_sws(swses, energies0, delta)
# 2 points is not enought to calculate equation of state
enough = False
# Loop until we have enough points to calculate initial sws
iteration = 0
while not enough:
iteration += 1
# if 25 is not enough one should check initial sws
if iteration > 25:
print(
"SWS loop did not converge in {0} iterations!".format(iteration))
quit()
# Calculate next point
already = False
self.sws = next_sws
job = self.create_jobname(self.jobname)
self.emto.set_values(sws=self.sws, jobname=job)
# First check if calculations are already done
for i in range(ca_mesh_len):
# We have to remember to update the lattice path every
# time we change c/a
self.emto.set_values(latpath=ca_latpaths[i])
hcp_subfolder = self.folder + "/{0}".format(ca_prefixes[i])
already = self.check_conv(job, folder=hcp_subfolder)
if self.lc_skip or (all(already) and not self.lc_rerun):
conv = (True, True)
else:
if already[0] and not already[1]:
conv = self.runemto(
jobname=job, folder=hcp_subfolder, onlyKFCD=True)
else:
conv = self.runemto(jobname=job, folder=hcp_subfolder)
# KGRN has crashed, find out why
if conv[0] == False:
self.which_error(job, folder=hcp_subfolder)
quit()
en = self.get_energy(job, folder=hcp_subfolder, func=xc)
energies[i, iteration + 1] = en
swses.append(self.sws)
# Calculate energetically optimized c/a and the corresponding energy
# print(energies[:,:15])
ca0, en0 = eos_hcp.ca_fit(
ca_mesh, energies[:, iteration + 1], ca_fit_order)
cas0.append(ca0)
energies0.append(en0)
# Check do we have minumun and predict next sws
next_sws, minfound = self.predict_next_sws(swses, energies0, delta)
# Check if we have mimimum and enough points
# 9 points should be enough for EOS
if minfound and len(swses) > 9:
enough = True
#print('debug find_lc_hcp: ',swses)
if prn:
self.print_sws_ens_hcp('find_lc_hcp', swses, energies0, cas0)
sws0, e0, B0, grun = eos_hcp.fit(swses, energies0)
# These functions create files on disk about the data to be fitted
# as well as the results of the fit.
# eos.prepareData()
#sws0,e0,B0,grun = eos.fit2file()
# Now that we have the ground state WS-radius sws0 we can use the cas0 array
# to compute the corresponding ground state ca0
ca0_vs_sws = np.polyfit(swses, cas0, 2)
ca0_vs_sws = np.poly1d(ca0_vs_sws)
c_over_a0 = ca0_vs_sws(sws0)
return sws0, c_over_a0, B0, e0
[docs] def refine_lc_hcp(self, sws, ca0, delta=0.01, prn=True, xc='PBE'):
"""Calculates a more accurate equilibrium volume for hcp systems.
:param sws: WS-radius
:type sws: float
:param ca0: Previously computed eq. c/a ratio
:type ca0: float
:param delta: Step size for the volume vs. energy array (Default value = 0.01)
:type delta: float
:param prn: True if results should be printed, False if not (Default value = True)
:type prn: boolean
:param xc: Choice of the xc-functional (Default value = 'PBE')
:type xc: str
:returns: Ws-radius, bulk modulus, c/a and energy
:rtype: float, float, float, float
"""
from pyemto.EOS.EOS import EOS
eos_hcp = EOS(name=self.jobname, xc=xc, method='morse', units='bohr')
if prn:
print('')
print('*****refine_lc_hcp*****')
print('')
# First compute the structure with the optimized
# c/a (found by the find_lc_hcp() function).
self.lattice.set_values(ca=ca0)
ca_prefix = "/ca_{0:8.6f}".format(ca0)
ca_latpath = self.latpath + ca_prefix
str_exists = self.check_str(self.latname, folder=ca_latpath)
if str_exists == False:
self.runlattice(jobname=self.latname, folder=ca_latpath)
# make sws ranges around given sws
points = []
for i in range(-3, 7):
points.append(i * delta)
energies = []
swses = []
self.emto.set_values(latpath=ca_latpath)
hcp_subfolder = self.folder + "/{0}".format(ca_prefix)
for p in points:
already = False
next_sws = sws + p
# Make inputs and job name
self.sws = next_sws
job = self.create_jobname(self.jobname)
self.emto.set_values(sws=self.sws, jobname=job)
# check if calculations are already done
already = self.check_conv(job, folder=hcp_subfolder)
if self.lc_skip or (all(already) and not self.lc_rerun):
conv = (True, True)
else:
if already[0] and not already[1]:
conv = self.runemto(
jobname=job, folder=hcp_subfolder, onlyKFCD=True)
else:
conv = self.runemto(jobname=job, folder=hcp_subfolder)
# KGRN has crashed, find out why
if conv[0] == False:
self.which_error(job, folder=hcp_subfolder)
quit()
en = self.get_energy(job, folder=hcp_subfolder, func=xc)
energies.append(en)
swses.append(self.sws)
if prn:
self.print_sws_ens('refine_lc_hcp', swses, energies)
sws0, e0, B, grun = eos_hcp.fit(swses, energies)
c_over_a = ca0
# These functions create files on disk about the data to be fitted
# as well as the results of the fit.
# eos.prepareData()
#sws0,e0,B,grun = eos.fit2file()
return sws0, B, c_over_a, e0
[docs] def predict_next_sws(self, swses, en, maxdelta=0.05):
"""Predict next WS-radius based on a simple gradient descent algorithm.
:param swses: List of current WS-radii
:type swses: list(float)
:param en: List of current energies
:type en: list(float)
:param maxdelta: Maximum step size (Default value = 0.05)
:type maxdelta: float
:returns: Next WS-radii and True if energy minimum has been found,
False if not yet found
:rtype: float, boolean
"""
# Check if we do have minimum here and predict directon for next point
m = 0 # The first one is minimum at start
minsws = 10000.0
maxsws = 0.0
for i in range(len(swses)):
# Check if we have minimum volume
if swses[i] < minsws:
minsws = swses[i]
elif swses[i] > maxsws:
maxsws = swses[i]
if en[i] < en[m]:
m = i
wehaveminumum = False
# Possible cases
if swses[m] == min(swses): # No minimum decrease sws
# One should make delta depend on energy slope at some point.
delta = -maxdelta
newsws = swses[m] + delta
elif swses[m] == max(swses): # No minimum increase sws
delta = maxdelta
newsws = swses[m] + delta
else:
wehaveminumum = True
# Minimum is inside the sws range. Decide where to add new point
larger = [s for s in swses if swses[m] < s]
smaller = [s for s in swses if swses[m] > s]
if len(larger) > len(smaller):
# Decrease volume
sws = minsws
delta = -maxdelta
newsws = sws + delta
else:
# Increase volume
sws = maxsws
delta = maxdelta
newsws = sws + delta
return newsws, wehaveminumum
[docs] def print_sws_ens(self, string, swses, energies):
"""Prints the WS-radii and calculated energies of cubic systems
:param string: Header for the printout
:type string: str
:param swses: List of WS-radii
:type swses: list(float)
:param energies: List of energies
:type energies: list(float)
:returns: None
:rtype: None
"""
str_len = len(string)
print('' + '\n')
print('************')
print(str_len * '*')
print(string)
print(str_len * '*')
print('************')
print(' SWS Energy')
for i in range(len(swses)):
print('{0:8.6f} {1:12.6f}'.format(swses[i], energies[i]))
print('' + '\n')
return
[docs] def print_sws_ens_hcp(self, string, swses, energies, cas):
"""Prints the WS-radii, ca/a ratios and calculated energies of hcp systems
:param string: Header for the printout
:type string: str
:param swses: List of WS-radii
:type swses: list(float)
:param energies: List of energies
:type energies: list(float)
:param cas: List of c/a ratios
:type cas: list(float)
:returns: None
:rtype: None
"""
str_len = len(string)
print('' + '\n')
print(str_len * '*')
print(string)
print(str_len * '*')
print(' SWS Energy0 c/a0')
for i in range(len(swses)):
print('{0:8.6f} {1:12.6f} {2:8.6f}'.format(
swses[i], energies[i], cas[i]))
print('' + '\n')
return
[docs] def check_conv(self, jobname, folder="./"):
"""Checks the convergence of given KGRN and KFCD calculations by reading
their output files.
:param jobname: Name of the output files
:type jobname: str
:param folder: Name of the folder where the output files are located (Default value = "./")
:type folder: str
:returns: Convergence of KGRN (True/False), convergence of KFCD (True/False)
:rtype: boolean, boolean
"""
folderKGRN = folder + '/kgrn/'
# Check if we got convergence in KGRN
prntfile = jobname + ".prn"
convergedKGRN = False
fn = os.path.join(folderKGRN, prntfile)
if not os.path.isfile(fn):
pass
else:
pfile = open(fn, "r")
for line in pfile:
if "Converged in" in line:
convergedKGRN = True
break
pfile.close()
folderKFCD = folder + '/kfcd/'
# Check if we got convergence in KFCD
prntfile = jobname + ".prn"
convergedKFCD = False
fn = os.path.join(folderKFCD, prntfile)
if not os.path.isfile(fn):
pass
else:
pfile = open(fn, "r")
for line in pfile:
if "Finished at:" in line:
convergedKFCD = True
break
pfile.close()
return convergedKGRN, convergedKFCD
[docs] def runlattice(self, jobname=None, folder="./", EMTODIR=None):
"""Run BMDL, KSTR and SHAPE calculation **WITHOUT** using the batch system.
:param jobname: Name of the input files (Default value = None)
:type jobname: str
:param folder: Name of the folder where the input files are located (Default value = "./")
:type folder: str
:param EMTODIR: Path to the EMTO installation folder (Default value = None)
:type EMTODIR: str
:returns: None
:rtype: None
"""
if EMTODIR is None:
EMTODIR = self.EMTOdir
if jobname is None:
sys.exit("System.runlattice: \'jobname\' has to be given!")
# Make sure folders exist
common.check_folders(folder, folder + "/bmdl/")
# Create input file
self.lattice.bmdl.write_input_file(folder=folder)
# Run BMDL
job = os.path.join(folder, jobname)
command = "cd {0}; ".format(
folder) + EMTODIR + "/bmdl/source/bmdl < " + job + ".bmdl"
print("running: " + command)
starttime = time.time()
os.system(command)
endtime = time.time()
timetaken = endtime - starttime
timehours = int(timetaken // 3600)
timeminutes = int((timetaken - timehours * 3600) // 60)
timeseconds = (timetaken - timehours * 3600) - timeminutes * 60
print("Finished running: " + command + '. Time: {0:3}h {1:2}m {2:5.2f}s '
.format(timehours, timeminutes, timeseconds) + '\n')
# Make sure folders exist
common.check_folders(folder, folder + "/kstr/")
# Create input file
self.lattice.kstr.write_input_file(folder=folder)
# Run KSTR
job = os.path.join(folder, jobname)
command = "cd {0}; ".format(
folder) + EMTODIR + "/kstr/source/kstr < " + job + ".kstr"
print("running: " + command)
starttime = time.time()
os.system(command)
endtime = time.time()
timetaken = endtime - starttime
timehours = int(timetaken // 3600)
timeminutes = int((timetaken - timehours * 3600) // 60)
timeseconds = (timetaken - timehours * 3600) - timeminutes * 60
print("Finished running: " + command + '. Time: {0:3}h {1:2}m {2:5.2f}s '
.format(timehours, timeminutes, timeseconds) + '\n')
if self.lattice.kstr.twocenter == True:
command = "cd {0}; ".format(
folder) + EMTODIR + "/kstr/source/kstr < " + job + "2.kstr"
print("running: " + command)
starttime = time.time()
os.system(command)
endtime = time.time()
timetaken = endtime - starttime
timehours = int(timetaken // 3600)
timeminutes = int((timetaken - timehours * 3600) // 60)
timeseconds = (timetaken - timehours * 3600) - timeminutes * 60
print("Finished running: " + command + '. Time: {0:3}h {1:2}m {2:5.2f}s '
.format(timehours, timeminutes, timeseconds) + '\n')
# Make sure folders exist
common.check_folders(folder, folder + "/shape/")
# Create input file
self.lattice.shape.write_input_file(folder=folder)
# Run SHAPE
job = os.path.join(folder, jobname)
command = "cd {0}; ".format(
folder) + EMTODIR + "/shape/source/shape < " + job + ".shape"
print("running: " + command)
starttime = time.time()
os.system(command)
endtime = time.time()
timetaken = endtime - starttime
timehours = int(timetaken // 3600)
timeminutes = int((timetaken - timehours * 3600) // 60)
timeseconds = (timetaken - timehours * 3600) - timeminutes * 60
print("Finished running: " + command + '. Time: {0:3}h {1:2}m {2:5.2f}s '
.format(timehours, timeminutes, timeseconds) + '\n')
return
[docs] def runemto(self, jobname, folder="./", EMTODIR=None, onlyKFCD=False):
"""Run KGRN and KFCD **WITHOUT** using the batch system and check convergence.
:param jobname: Name of the input files
:type jobname: str
:param folder: Name of the folder where the input files are located (Default value = "./")
:type folder: str
:param EMTODIR: Path to the EMTO installation folder (Default value = None)
:type EMTODIR: str
:param onlyKFCD: True if only KFCD needs to calculated,
False if also KGRN (Default value = False)
:type onlyKFCD: boolean
:returns: True if the calculations converged, False if not
:rtype: boolean
"""
if jobname is None:
sys.exit("System.runemto: \'jobname\' has to be given!")
if EMTODIR is None:
EMTODIR = self.EMTOdir
if onlyKFCD:
# Make sure folders exist
common.check_folders(folder, folder + "/kfcd")
# Run KFCD
self.emto.kfcd.write_input_file(folder=folder)
job = os.path.join(folder, jobname)
command = "cd {0}; ".format(
folder) + EMTODIR + "/kfcd/source/kfcd_cpa < " + job + ".kfcd"
print("running: " + command)
starttime = time.time()
os.system(command)
endtime = time.time()
timetaken = endtime - starttime
timehours = int(timetaken // 3600)
timeminutes = int((timetaken - timehours * 3600) // 60)
timeseconds = (timetaken - timehours * 3600) - timeminutes * 60
print("Finished running: " + command + '. Time: {0:3}h {1:2}m {2:5.2f}s '
.format(timehours, timeminutes, timeseconds) + '\n')
# Check if we got convergence
converged = self.check_conv(jobname, folder=folder)
else:
# Make sure folders exist
common.check_folders(folder, folder + "/kgrn", folder + "/kgrn/tmp")
# Run KGRN
self.emto.kgrn.write_input_file(folder=folder)
job = os.path.join(folder, jobname)
command = "cd {0}; ".format(
folder) + EMTODIR + "/kgrn/source/kgrn_cpa < " + job + ".kgrn"
print("running: " + command)
starttime = time.time()
os.system(command)
endtime = time.time()
timetaken = endtime - starttime
timehours = int(timetaken // 3600)
timeminutes = int((timetaken - timehours * 3600) // 60)
timeseconds = (timetaken - timehours * 3600) - timeminutes * 60
print("Finished running: " + command + '. Time: {0:3}h {1:2}m {2:5.2f}s '
.format(timehours, timeminutes, timeseconds))
# Check if we got convergence
converged = self.check_conv(jobname, folder=folder)
if converged[0]:
# Make sure folders exist
common.check_folders(folder, folder + "/kfcd")
# Run KFCD
self.emto.kfcd.write_input_file(folder=folder)
job = os.path.join(folder, jobname)
command = "cd {0}; ".format(folder) + EMTODIR +\
"/kfcd/source/kfcd_cpa < " + job + ".kfcd"
print("running: " + command)
starttime = time.time()
os.system(command)
endtime = time.time()
timetaken = endtime - starttime
timehours = int(timetaken // 3600)
timeminutes = int((timetaken - timehours * 3600) // 60)
timeseconds = (timetaken - timehours * 3600) - timeminutes * 60
print("Finished running: " + command + '. Time: {0:3}h {1:2}m {2:5.2f}s '
.format(timehours, timeminutes, timeseconds) + '\n')
converged = self.check_conv(jobname, folder=folder)
return converged
[docs] def which_error(self, jobname, folder="./"):
"""Tries to determine the reason why a given KGRN calculation did not converge.
The reason for the crash will be printed on screen.
:param jobname: Name of the KGRN output file
:type jobname: str
:param folder: Name of the folder where the output file is located (Default value = "./")
:type folder: str
:returns: None
:rtype: None
"""
folder = folder + '/kgrn/'
errdict = {'EFXMOM': 'EFXMOM:** Fermi level not found',
'NOTCONV': 'Not converged'}
prntfile = jobname + ".prn"
fn = os.path.join(folder, prntfile)
if not os.path.isfile(fn):
print('Function which_error: file {0} does not exist!'.format(fn))
quit()
pfile = open(fn, "r")
done = False
# knownError = failure caused by a known reason in the errdict.
# Otherwise cause is a more exotic KGRN error or the run crashed.
knownError = False
for line in pfile:
for key in errdict:
if errdict[key] in line:
errorkey = key
errormesg = line
done = True
knownError = True
break
if done:
break
pfile.close()
if knownError == True:
print('Problem in KGRN for {0}: {1}'.format(fn, errormesg))
else:
print('Problem in KGRN for {0}: Program has crashed or stopped due to a' +
' less common KGRN error. Check the output file.'.format(fn))
return
[docs] def get_energy(self, jobname=None, func="PBE", folder=None):
"""Extracts total energy from the KFCD output file.
Different total energies given by different xc-functionals can
be selected using the *func* input parameter. Default value is
'PBE'.
:param jobname: Name of the KFCD output file
:type jobname: str
:param func: Name of the xc-functional (Default value = "PBE")
:type func: str
:param folder: Name of the folder where the output file is located (Default value = "./")
:type folder: str
:returns: Total energy if it is found, otherwise return None
:rtype: float or None
"""
if folder == None:
folder = self.folder
if jobname == None:
jobname = self.fulljobname
fn = os.path.join(folder, "kfcd/")
fn = os.path.join(fn, jobname + ".prn")
try:
kfcdfile = open(fn, "r")
energy = "TOT-" + func
energyFound = False
for line in kfcdfile:
if energy in line:
energyFound = True
break
# Total energy of the unit cell
# return float(line.split()[1])
# Energy per WS-radius: Better to use this to get Bmod correctly
if energyFound:
return float(line.split()[3])
else:
return
except IOError:
print('System.get_energy(): {0} does not exist!'.format(fn))
return
[docs] def get_moments(self, jobname=None, func="PBE", folder=None):
"""Extracts magnetic moments from the KFCD output file.
:param jobname: Name of the KFCD output file
:type jobname: str
:param func: Name of the xc-functional (Default value = "PBE")
:type func: str
:param folder: Name of the folder where the output file is located (Default value = "./")
:type folder: str
:returns: Total energy if it is found, otherwise return None
:rtype: float or None
"""
if folder == None:
folder = self.folder
if jobname == None:
jobname = self.fulljobname
fn = os.path.join(folder, "kfcd/")
fn = os.path.join(fn, jobname + ".prn")
readyTag = "KFCD: Finished at:"
momTag = "Magnetic moment for IQ ="
try:
kfcdfile = open(fn, "r")
lines = kfcdfile.readlines()
kfcdfile.close()
moms = []
ready = False
for line in lines:
if readyTag in line:
ready = True
if ready:
for line in lines:
#if enTag in line:
# linesplit = line.split()
# sws = float(linesplit[6])
if momTag in line:
linesplit = line.split()
moms.append(float(linesplit[7]))
return moms
else:
return None
except IOError:
print('System.get_moments(): {0} does not exist!'.format(fn))
return None
[docs] def get_fdos(self, jobname = None,folder=None):
""" Extract density of state (DOS) at fermi level from KGRN output
:param jobname: Name of the KGRN output file
:type jobname: str
:param folder: Name of the folder where the output file is located (Default value = "./")
:type folder: str
:returns: DOS at fermi level
:rtype: float
"""
if folder == None:
folder = self.folder
if jobname == None:
jobname = self.fulljobname
fn = os.path.join(folder, "kgrn/")
fn = os.path.join(fn, jobname + ".prn")
file = open(fn,'r')
lines = file.readlines()
file.close()
nxtsws_tag = " NXTSWS: IQ IT ITA MMT Type NZ ION ELN QTR SPLIT FIX CONC"
alat_tag = "Alat ="
dos_tag = " Dos(Ef) ="
mag_tag = " Magn. mom. ="
hop_tag = " Hopfield ="
for i in range(len(lines)):
if nxtsws_tag in lines[i]:
indMin = i+2
elif alat_tag in lines[i]:
indMax = i-2
break
#
concs = np.zeros(indMax + 1 - indMin)
doses = np.zeros(indMax + 1 - indMin)
its = np.zeros(indMax + 1 - indMin)
#
ind_tmp = 0
for i in range(indMin,indMax+1):
concs[ind_tmp] = float(lines[i].split()[-1])
its[ind_tmp] = int(lines[i].split()[1])
ind_tmp += 1
#
num_sites = np.max(its)
ind_tmp = len(doses) - 1
# Because KGRN output file is different for non- and magnetic calculations,
# we have to do some additional checks to make sure we are reading the right
# values.
for i in range(len(lines)-1,indMax,-1):
if dos_tag in lines[i]:
#print(lines[i])
if mag_tag in lines[i+1] or hop_tag in lines[i+1]:
#print(lines[i+1])
doses[ind_tmp] = float(lines[i].split()[-1])
ind_tmp -= 1
if ind_tmp == -1:
break
#
#for i in range(len(doses)):
# print(doses[i])
#
ry2ev = 13.605698066
dos_tot = 0.0
for i in range(len(concs)):
dos_tot += concs[i]*doses[i]
dos_tot /= num_sites
dos_tot /= ry2ev
return dos_tot
[docs] def create_jobname(self, jobname=None):
"""Creates jobnames based on system information.
Creates a jobname based on the optional input prefix *jobname*.
The format of the full jobname in this case will be:
jobname_3.000000, where 3.000000 is the sws.
If *jobname* prefix is not given, full jobname is generated automatically
based on system data self.sws, self.atoms and self.concs.
The format of the jobname is: au0.50pd0.50_3.000000, where 0.50 are the
atomic concentrations and 3.000000 is the sws.
:param jobname: Optional prefix for the full jobname (Default value = None)
:type jobname: str
:returns: Newly created full jobname
:rtype: str
"""
if jobname is None:
# Prefix not specified => create a default jobname based on
# the types of atoms and their concentrations.
jobname = ''
for i in range(len(self.atoms)):
if jobname == "":
pass
else:
jobname = jobname + "_"
jobname = jobname + \
self.atoms[i].lower() + "%4.2f" % (self.concs[i])
fulljobname = jobname + "_%8.6f" % (self.sws)
return jobname, fulljobname
else:
fulljobname = jobname + "_{0:8.6f}".format(self.sws)
return fulljobname
[docs] def check_str(self, jobname, folder="./"):
"""Checks if a KSTR calculation file exists and has converged.
Only KSTR file is checked because BMDL and SHAPE are fast to
rerun in any case.
:param jobname: Filename of the structure output file
:type jobname: str
:param folder: Folder where the output file is located (Default value = "./")
:type folder: str
:returns: True if KSTR has finished and False if not
:rtype: boolean
"""
folderKSTR = folder + '/kstr/'
# Check if we got convergence in KSTR
prntfile = jobname + ".prn"
convergedKSTR = False
fn = os.path.join(folderKSTR, prntfile)
if not os.path.isfile(fn):
pass
else:
pfile = open(fn, "r")
for line in pfile:
if "KSTR: Finished at :" in line:
convergedKSTR = True
break
pfile.close()
return convergedKSTR
[docs] def wait_for_jobs(self, jobsdict, restart_partition='general', sleeptime=60, restart_z=None,
restart_stragglers_after=0.75, kill_if_all_ssusp=False):
"""Loops checking status until no jobs are waiting or running / all are finished.
wait/run states:
======= =========== ==================================================================================================
**Key** **Meaning** **Description**
------- ----------- --------------------------------------------------------------------------------------------------
CF CONFIGURING Job has been allocated resources, but are waiting for them to become ready for use (e.g. booting).
CG COMPLETING Job is in the process of completing. Some processes on some nodes may still be active.
PD PENDING Job is awaiting resource allocation.
R RUNNING Job currently has an allocation.
RS RESIZING Job is about to change size.
S SUSPENDED Job has an allocation, but execution has been suspended.
======= =========== ==================================================================================================
done states:
======= =========== ==============================================================================================================
**Key** **Meaning** **Description**
------- ----------- --------------------------------------------------------------------------------------------------------------
CA CANCELLED Job was explicitly cancelled by the user or system administrator. The job may or may not have been initiated.
CD COMPLETED Job has terminated all processes on all nodes.
F FAILED Job terminated with non-zero exit code or other failure condition.
NF NODE_FAIL Job terminated due to failure of one or more allocated nodes.
PR PREEMPTED Job terminated due to preemption.
TO TIMEOUT Job terminated upon reaching its time limit.
======= =========== ==============================================================================================================
:param jobsdict:
:type jobsdict:
:param restart_partition: (Default value = 'general')
:type restart_partition:
:param sleeptime: (Default value = 60)
:type sleeptime:
:param restart_z: (Default value = None)
:type restart_z:
:param restart_stragglers_after: (Default value = 0.75)
:type restart_stragglers_after:
:param kill_if_all_ssusp: (Default value = False)
:type kill_if_all_ssusp:
:returns: None
:rtype: None
"""
run_status = ['CONFIGURING', 'COMPLETING',
'PENDING', 'RUNNING', 'RESIZING', 'SUSPENDED']
done_status = ['CANCELLED', 'COMPLETED',
'FAILED', 'NODE_FAIL', 'PREEMPTED', 'TIMEOUT']
import time
import datetime
jobs_amount = len(jobsdict)
print()
print('wait_for_jobs: Submitted {0} jobs'.format(jobs_amount))
status = self.get_status_counts(jobsdict)
t = time.time()
maxllen = 0
print('wait_for_jobs: Will be requesting job statuses' +
' every {0} seconds'.format(sleeptime) + "\n")
while any([k in run_status for k in status.keys()]):
time.sleep(sleeptime)
status = self.get_status_counts(jobsdict)
pctdone = sum([status.get(rs, 0)
for rs in done_status]) / float(sum(status.values()))
# CHECK SUSPENDED; RESTART STRAGGLERS, ETC
outl = '%s %s (%3d%% completion)' % (str(
datetime.timedelta(seconds=int(time.time() - t))), status.__repr__(), pctdone * 100)
# if len(outl) < maxllen:
# pad = maxllen - len(outl)
# outl += ' '*pad
# else:
# maxllen = len(outl)
print(outl)
# sys.stderr.write(outl)
# sys.stderr.flush()
print('completed {0} batch jobs in {1}'.format(
jobs_amount, str(datetime.timedelta(seconds=int(time.time() - t)))))
return
[docs] def get_status_counts(self, jobids=None):
"""Returns the counts of all jobs by status category.
:param jobids: (Default value = None)
:type jobids:
:returns:
:rtype:
"""
from collections import defaultdict
jobs_status = self.get_jobs_status(jobids)
status_counts = defaultdict(int)
for jd in jobs_status.values():
status_counts[jd['State'].split()[0]] += 1
return dict(status_counts)
[docs] def get_jobs_status(self, jobids=None, toplevel=True):
"""Returns status of the jobs indicated
(jobsdict or list of job ids) or all jobs if no jobids supplied.
Set toplevel=False for job step data.
:param jobids: List of job IDs (Default value = None)
:type jobids: dict or list
:param toplevel: (Default value = True)
:type toplevel: boolean
:returns: Job statuses
:rtype: list(str)
"""
import subprocess
if jobids is None:
sacct_return = subprocess.Popen(
'sacct -p -l', shell=True, stdout=subprocess.PIPE).stdout.readlines()
else:
if isinstance(jobids, dict):
qjobs = jobids.keys()
else:
qjobs = jobids
sacct_return = subprocess.Popen(
'sacct -j %s -p -l' % (
','.join(qjobs),), shell=True, stdout=subprocess.PIPE).stdout.readlines()
jobs_status = {}
for el in sacct_return[1:]:
d = dict(
zip(sacct_return[0].strip().split('|'), el.strip().split('|')))
if not '.' in d['JobID'] or not toplevel:
jobs_status[d['JobID']] = d
return jobs_status
[docs] def submit_jobs(self, jobnames, folder=None):
"""Takes a list of jobnames and submits the corresponding
batch scripts to the batch system.
:param jobnames: List of jobnames
:type jobnames: list(float)
:param folder: Folder where the batch scripts are located (Default value = None)
:type folder: str
:returns: List of job ID numbers
:rtype: list(str)
"""
import time
from pyemto.utilities.utils import run_emto
sleeptime = 10
if folder is None:
folder = self.folder
job_ids = []
for i in range(len(jobnames)):
job_ids.append(run_emto(jobnames[i], folder=self.folder))
# Flatten job_ids list and convert the integers into strings
job_ids = [item for sublist in job_ids for item in sublist]
for i in range(len(job_ids)):
job_ids[i] = str(job_ids[i])
# Give SLURM some time to register the jobs.
# If continued straight away the self.wait_for_jobs
# script will likely think all the jobs have finished
# since it cannot find them yet.
time.sleep(sleeptime)
return job_ids