<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Cu-Mg workflows
Goal
Step1:  Get your structures
Three ways main ways you'll use to get structures
1. From the Materials Project via the MPRester API
Step2:  2. From a POSCAR file
Step3:  3. From an SQS in prlworkflows
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,
Download ATAT from https
Step4:  Run the robust optimization
The robust optimization will do a constrained optimization of your structures.
There are typically three kinds of degrees of freedom we can control
Step5:  Mixing robust optimization
Step6:  Finding the minimum energy structure
Now that we've completed the robust calculations, we want to find the minimal energy strucure for each robust optimization.
The general steps are
Step7:  Get all the structures
Step8:  Get the most optimized constrained structure
This is essentially a collection of heuristics on matrix norms and symmetry.
They are not optimized yet, don't trust the tolerances!
Step9:  Gibbs energies and unconstrained formation energies
With the minimum energy constrained structures, we want to find
The finite temperature data (Gibbs free energies) for the constrained structure
The unconstrained formation energy
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.
Gibbs free energy
<ASSISTANT_TASK:>
Python Code:
from fireworks import LaunchPad
# lpad = LaunchPad.auto_load()
lpad = LaunchPad.from_file('/Users/brandon/.fireworks/my_launchpad.yaml')
Explanation: Cu-Mg workflows
Goal: fully describe the Cu-Mg system with DFT calculations
Phases
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.
Approach
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.
General imports and setup
Configuration
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/codes
git clone https://github.com/phasesresearchlab/prlworkflows
cd prlworkflows
pip install -e .
Finally, you need to do a little configuration on the system(s) you will be adding and running workflows from:
Add - prlworkflows to the list of ADD_USER_PACKAGES in $MY_ATOMATE_INSTALL/config/FW_CONFIG.yaml
Change 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 
   python
    def set_analysis(d, max_force_threshold=0.5, volume_change_threshold=0.2):
   to 
   python
    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:
End of explanation
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')
Explanation: Get your structures
Three ways main ways you'll use to get structures
1. From the Materials Project via the MPRester API
End of explanation
from pymatgen import Structure
structure = Structure.from_file('POSCAR')
Explanation: 2. From a POSCAR file
End of explanation
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']])
    
Explanation: 3. From an SQS in prlworkflows
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,
Download ATAT from https://www.brown.edu/Departments/Engineering/Labs/avdw/atat/
Unzip the archive
Go to the folder atat/data/sqsdb, there will be a list of structures and levels
Pick the phase/level/structure that you want, then get the file path to the bestsqs.out file.
Process the file with lat_in_to_sqs creating an AbstractSQS
Make the AbstractSQS concrete with a sublattice model
End of explanation
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<<')
lpad.add_wf(wf)
Explanation: Run the robust optimization
The robust optimization will do a constrained optimization of your structures.
There are typically three kinds of degrees of freedom we can control:
Unit cell volume
Unit cell shape (lattice vectors; length and angles of unit cell)
Ion positions
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:
Volume relaxations until the volume has converged
An ionic position relaxation
A shape and ionic position relaxation
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.
Formation robust optimization
End of explanation
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']])
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<<')
lpad.add_wf(wf)
Explanation: Mixing robust optimization
End of explanation
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)
# 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
Explanation: Finding the minimum energy structure
Now that we've completed the robust calculations, we want to find the minimal energy strucure for each robust optimization.
The general steps are: 
Load the database where calculation results are stored (db.json)
Query for the structure
Check if the ion positions changed significantly from relaxation step 2. If it did not change continue, otherwise the minimum energy structure is the volume relaxed structure from relaxation step 1.
Check if the cell shape and ion positions changed significantly in relaxation step 3. If they did not change, the final structure is our minimum energy structure. If either changed, then our minimum energy structure is the cell shape relaxed (relaxation step 2) structure.
End of explanation
# 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]
Explanation: Get all the structures
End of explanation
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
stable_struct = find_stable_from_robust(*relaxed_structs)
Explanation: Get the most optimized constrained structure
This is essentially a collection of heuristics on matrix norms and symmetry.
They are not optimized yet, don't trust the tolerances!
End of explanation
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)
Explanation: Gibbs energies and unconstrained formation energies
With the minimum energy constrained structures, we want to find
The finite temperature data (Gibbs free energies) for the constrained structure
The unconstrained formation energy
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.
Gibbs free energy
End of explanation
<END_TASK>