There are 5 phases in Cu-Mg that will be described with the following models
| Phase | Model |
|---|---|
| $\textrm{liquid}$ | $\textrm{(Cu, Mg)}$ |
| $\textrm{fcc A1}$ | $\textrm{(Cu, Mg)}$ |
| $\textrm{hcp A3}$ | $\textrm{(Cu, Mg)}$ |
| $\textrm{Cu}\textrm{Mg}_2$ | $\textrm{(Cu)}\textrm{(Mg)}_2$ |
| $\textrm{Laves C15}$ | $\textrm{(Cu, Mg)}_2 \textrm{(Cu, Mg)}$ |
Four of these phases are solid phases that can be investigated with DFT.
The others are modeled as mixing phases, and thus the mixing properties must be calculated.
All compounds, solution endmembers and mixing structures will have their thermodynamic properties calculated via the Debye model.
Currently, the contributions to the Debye model do not include the electronic contribution to the heat capacity. An electronic DOS can be calculated later and this contribution added.
Before you start anything, you'll have to install prlworkflows.
Right now there is no release version on PyPI or Anaconda, so you need to get clone and install it that way.
Since this code simply adds structures to your LaunchPad, you can run this porition of prlworkflows and atomate on any computer, so long as you set up the databases as in the atomate installation instructions.
Either way, I suggest using pip to install and doing the following in your terminal:
cd ~/work/atomate/codesgit clone https://github.com/phasesresearchlab/prlworkflowscd prlworkflowspip install -e .Finally, you need to do a little configuration on the system(s) you will be adding and running workflows from:
- prlworkflows to the list of ADD_USER_PACKAGES in $MY_ATOMATE_INSTALL/config/FW_CONFIG.yamlChange the force convergence criteria in the drone analysis. If you fail to do this, relaxations that have some residual force will defuse the next Fireworks. If you don't have atomate installed as editiable from git, then you need to
a. Go to your virtual environment directory cd ~/work/atomate/atomate_env
b. vi lib/python2.7/site-packages/atomate/vasp/drones.py
c. Find the VaspDrone.set_analysis method and change the default arguments from
def set_analysis(d, max_force_threshold=0.5, volume_change_threshold=0.2):
to
def set_analysis(d, max_force_threshold=10, volume_change_threshold=0.2):
Now we can get to running the code. First, set up the LaunchPad using your own credentials:
In [1]:
from fireworks import LaunchPad
# lpad = LaunchPad.auto_load()
lpad = LaunchPad.from_file('/Users/brandon/.fireworks/my_launchpad.yaml')
In [3]:
from pymatgen import MPRester
with MPRester() as mpr: # provide an API key if you don't have on set up in your .pmgrc.yaml
structure = mpr.get_structure_by_material_id('mp-134')
In [ ]:
from pymatgen import Structure
structure = Structure.from_file('POSCAR')
Until the SQS database is complete and queryable, these will have to be created 'by hand' from ATAT lattice.in type files. To get all of these,
atat/data/sqsdb, there will be a list of structures and levelslat_in_to_sqs creating an AbstractSQS
In [2]:
from prlworkflows.sqs_db import lat_in_to_sqs
my_struct_best_sqs_path = '/Users/brandon/Downloads/atat/data/sqsdb/MGCU2_C15/sqsdb_lev=0_b=1_c=1/bestsqs.out'
with open(my_struct_best_sqs_path) as fp:
lat_in_txt = fp.read()
abstract_sqs = lat_in_to_sqs(lat_in_txt)
structure = abstract_sqs.get_concrete_sqs([['Cu'], ['Cu']])
The robust optimization will do a constrained optimization of your structures.
There are typically three kinds of degrees of freedom we can control:
If we fix any one of these degrees of freedom, we are doing a constrained relaxation, otherwise if all degrees of freedom are considered in a relaxation, it is called an unconstrained relaxation.
Most stable structures that are already close to their minimum energy configuration are typically easy to perform unconstrained relaxations on. All of the structures in the Materials Project are this way.
We would like to find the minimum energy configuration of our endmembers, SQS, and dilute mixing structures. However, many SQS, some endmembers, and some dilute mixing structures are unstable, and the minimum energy structures found in a unconstrtained optimization are not representative of the real structure and energy.
Thus we have developed a series of constrained optimizations to find the lowest energy structure while maintaining the integrity of the structure. It is called the robust optimization workflow, which performs the following steps:
The minimum energy structure we will chose is the one that has gone the furthest in the series, but still maintained the symmetry within a tolerance. The volume relaxation always will, but the ionic positions and cell shape + ionic positions relaxations may result in broken structures.
For now these are all run one after another and we have to find the structure that progressed the furthest 'by hand' in code. In the future we will modify the workflow to check that the structure has not broken after step 2 and 3 to find the constrained structure of minimal energy and preserved symmetry.
Note: for now we only run VASP one time in step 2 and 3. It is possible that we would like to run a series of these relaxations as in step 1, to make sure that each relaxation type is fully converged w.r.t. itself, but there is currently no machinery for that. It would have to be implemented either in custodian like the volume relaxation, or it would need to be implemented by adding a detour Firework to do another relaxation if some criteria isn't met.
In [3]:
from prlworkflows.prl_workflows import get_wf_robust_optimization
metadata = {
'phase': 'LAVES_C15',
'sublattice_model': {
# if you are using an SQS object, you can use these.
# otherwise you should set them by hand so you can query on them later
'configuration': structure.sublattice_configuration, # [['Cu'], ['Cu']]
'occupancies': structure.sublattice_occupancies, # [[1.0], [1.0]]
'site_ratios': structure.sublattice_site_ratios # [1, 2]
},
'output': 'HM_FORM' # To help filter things later. This is formation data, rather than mixing, because it's an endmembmer
}
wf = get_wf_robust_optimization(structure, metadata=metadata, db_file='>>db_file<<', vasp_cmd='>>vasp_cmd<<')
In [4]:
lpad.add_wf(wf)
Out[4]:
In [5]:
my_struct_best_sqs_path = '/Users/brandon/Downloads/atat/data/sqsdb/HCP_A3/sqsdb_lev=1_c=0.5,0.5/bestsqs.out'
with open(my_struct_best_sqs_path) as fp:
lat_in_txt = fp.read()
abstract_sqs = lat_in_to_sqs(lat_in_txt)
structure_mixing = abstract_sqs.get_concrete_sqs([['Cu', 'Mg']])
In [6]:
from prlworkflows.prl_workflows import get_wf_robust_optimization
metadata = {
'phase': 'HCP_A3',
'sublattice_model': {
# if you are using an SQS object, you can use these.
# otherwise you should set them by hand so you can query on them later
'configuration': structure_mixing.sublattice_configuration, # [['Cu', 'Mg']]
'occupancies': structure_mixing.sublattice_occupancies, # [[0.5, 0.5]]
'site_ratios': structure_mixing.sublattice_site_ratios # [2]
},
'output': 'HM_MIX' # To help filter things later
}
wf = get_wf_robust_optimization(structure_mixing, metadata=metadata, db_file='>>db_file<<', vasp_cmd='>>vasp_cmd<<')
In [7]:
lpad.add_wf(wf)
Out[7]:
Now that we've completed the robust calculations, we want to find the minimal energy strucure for each robust optimization.
The general steps are:
db.json)
In [13]:
from atomate.vasp.database import VaspCalcDb
# create the atomate db from your db.json
PATH_TO_MY_DB_JSON = '/Users/brandon/.fireworks/db.json'
atomate_db = VaspCalcDb.from_db_file(PATH_TO_MY_DB_JSON)
In [73]:
# the unique query to find a structure
structure_query = {'metadata.phase': 'HCP_A3',
'metadata.sublattice_model.configuration': [['Cu']]}
t = atomate_db.db.tasks.find(structure_query)
print('Calculations found : {}'.format(t.count())) # should have found 3 if this calculation is complete
In [123]:
# the oneline dict expressions require Python3
vol_task = atomate_db.db.tasks.find_one({**structure_query, 'orig_inputs.incar.ISIF': 7})
vol_relax_structure = Structure.from_dict(vol_task['output']['structure'])
ions_task = atomate_db.db.tasks.find_one({**structure_query, 'orig_inputs.incar.ISIF': 2})
ions_structure = Structure.from_dict(vol_task['output']['structure'])
shape_ions_task = atomate_db.db.tasks.find_one({**structure_query, 'orig_inputs.incar.ISIF': 4})
shape_ions_structure = Structure.from_dict(vol_task['output']['structure'])
relaxed_structs = [vol_relax_structure, ions_structure, shape_ions_structure]
In [131]:
import numpy as np
def check_ionic_positions(s1, s2, tolerance=0.1):
# this will totally fail if the order of the ions in the structure changed
norm = np.linalg.norm(s1.distance_matrix - s2.distance_matrix)
return norm < tolerance
def check_cell_shape(s1, s2, tolerance=0.1):
m1 = s1.lattice.matrix
m2 = s2.lattice.matrix
norm = np.linalg.norm(m1 - m2)
return norm < tolerance
def check_pure_symmetry(s1, s2, tolerance=0.1):
symm1 = make_pure_struct(s1).get_space_group_info(symprec=tolerance)
symm2 = make_pure_struct(s2).get_space_group_info(symprec=tolerance)
return symm1[1] == symm2[1]
def make_pure_struct(struct):
el0 = struct.types_of_specie[0].name
replacement_dict = {el.name: el0 for el in struct.types_of_specie[1:]}
struct_copy = struct.copy()
struct_copy.replace_species(replacement_dict)
return struct_copy
def find_stable_from_robust(vol_struct, ions_struct, shape_ions_struct):
if not check_ionic_positions(vol_struct, ions_struct):
print('Volume (step 1) is most stable')
return vol_struct
ions = check_ionic_positions(ions_structure, shape_ions_struct)
shape = check_cell_shape(ions_structure, shape_ions_struct)
symmetry = check_pure_symmetry(ions_structure, shape_ions_struct)
if not (ions and shape and symmetry):
print('Ionic position (step 2) is most stable')
return ions_struct
else:
print('Shape and ions (step 3) is most stable')
return shape_ions_struct
In [132]:
stable_struct = find_stable_from_robust(*relaxed_structs)
With the minimum energy constrained structures, we want to find
The reason for the finite temperature data is clear, but the unconstrained formation energy may seem less obvious. For unstable structures we often cannot reasonably reach the true minimum energy and structure (the so called lattice stability, see van de Walle et al. Nat. Commun. 6, 7559 (2015)). This minimum energy lies somewhere between the minimum constrained energy and the unconstrained energy. It's not uncommon for the constrained and unconstrained energies to differ by 10 kJ/mol-atom. Knowing unconstrained energy allows us to have some insight into how precise our constrained energy is compared to the true lattice stability.
In [ ]:
from prlworkflows.prl_workflows import wf_gibbs_free_energy
config_dict = {
'OPTIMIZE': False,
'T_MIN': 5,
'T_MAX': 2000,
'T_STEP' 5,
'POISSON': 0.32
}
wf = wf_gibbs_free_energy(stable_struct, config_dict)
In [ ]: