<jupyter_start><jupyter_text>[**Workflows for atomistic simulations**](http://potentials.rub.de/)   **Day 3 - Validation of various potentials** **Exercise 1: Validation of generated potentials using pyiron based workflows**Before the excercise, you should:* Have run the notebooks from day 1 and day2* Be familiar with working with pyiron and the basics of potential fittingThe aim of this exercise is to make you familiar with:* Potential validation workflows  **Importing the modules and creating the project**<jupyter_code>import numpy as np
%matplotlib inline
import matplotlib.pylab as plt
import pandas as pd
from pyiron import Project
pr = Project("validation")
# Potentials from fit the previous day
potential_list  = ['df1_cut5_pyace', 'EAM', 'RuNNer-Cu']

#  Typical EAM and well fit potentials
potential_list  = ['2012--Mendelev-M-I--Cu--LAMMPS--ipr1', '2004--Zhou-X-W--Cu-Ag-Au--LAMMPS--ipr2', 
                   'Cu-ace', 'Cu-runner-df4', 'Cu-atomicrex-df1-107-25'] #, 'Cu-runner-df1', 'Cu-atomicrex-df1-0-13']
# Function to return appropriate potentials generated yesterday

def get_pyiron_potential(potential):
    if potential in ["df1_cut5_pyace"]:
        return pr["../../day_2/03-ace/pacemaker_fit/df1_cut5_pyace"].get_lammps_potential()
    
    if potential in ["EAM"]:
        return pr["../../day_2/01-atomicrex/WorkshopPotential/PotentialDF1"].lammps_potential
    
    elif potential in ["RuNNer-Cu"]:
        return pr["../../day_2/02-runner/runner_fit/fit"].lammps_potential
    # to returns potentials already in DB
    else:
        return potential<jupyter_output><empty_output><jupyter_text>**Step 1: Get the equilibrium Cu lattice constant determined by the potentials**First we compute energy-volume curves to get the equilibrium lattice constants determined by each potential<jupyter_code># Necessary to remove wierd characters in the potential name
def clean_project_name(name):
    return name.replace("-", "_")
for pot in potential_list:
    group_name = clean_project_name(pot)
    pr_pot = pr.create_group(pot)
    job_ref = pr_pot.create_job(pr_pot.job_type.Lammps, "ref_job")
    job_ref.structure = pr_pot.create_ase_bulk("Cu")
    job_ref.potential = get_pyiron_potential(pot)
    job_ref.calc_minimize()
    murn_job = job_ref.create_job(pr_pot.job_type.Murnaghan, "murn_job")
    murn_job.run()
    print("Potential: ", pot)
    murn_job.plot()
murn_job["output/equilibrium_energy"]
# Define functions to get data

# Only work with Murnaghan jobs
def get_only_murn(job_table):
    return job_table.hamilton == "Murnaghan"

def get_eq_vol(job_path):
    return job_path["output/equilibrium_volume"]

def get_eq_lp(job_path):
    return np.linalg.norm(job_path["output/structure/cell/cell"][0]) * np.sqrt(2)

def get_eq_bm(job_path):
    return job_path["output/equilibrium_bulk_modulus"]

def get_potential(job_path):
    return job_path["ref_job/input/potential/Name"]

def get_eq_energy(job_path):
    return job_path["output/equilibrium_energy"]

def get_n_atoms(job_path):
    return len(job_path["output/structure/positions"])
# Compile data using pyiron tables
table = pr.create_table("table_murn", delete_existing_job=True)
table.db_filter_function = get_only_murn
table.add["potential"] = get_potential
table.add["a"] = get_eq_lp
table.add["eq_vol"] = get_eq_vol
table.add["eq_bm"] = get_eq_bm
table.add["eq_energy"] = get_eq_energy
table.add["n_atoms"] = get_n_atoms
table.run()
data_murn = table.get_dataframe()
data_murn<jupyter_output>100%|██████████| 5/5 [00:00<00:00, 997.41it/s]
100%|██████████| 5/5 [00:00<00:00, 43.30it/s]<jupyter_text>**Compute elastic constants**Using the equilibrium lattice constant, compute the elastic constant matrix using the automated workflow<jupyter_code>for pot in data_murn.potential.to_list():
    group_name = clean_project_name(pot)
    pr_pot = pr.create_group(pot)
    print(pot)
    job_id = int(data_murn[data_murn.potential==pot].job_id)
    job_ref = pr_pot.create_job(pr_pot.job_type.Lammps, "ref_job")
    job_ref.structure = pr_pot.inspect(job_id)["output/structure"].to_object()
    job_ref.potential = get_pyiron_potential(pot)
    job_ref.calc_minimize()
    elastic_job = job_ref.create_job(pr_pot.job_type.ElasticMatrixJob, "elastic_job")
    elastic_job.input["eps_range"] = 0.05
    elastic_job.run()
# Inspecting the elastic matrix
plt.matshow(elastic_job["output/elasticmatrix"]["C"]);
# Define functions to compute them in a table

# Only operate on ElasticMatrix jobs
def filter_elastic(job_table):
    return (job_table.hamilton == "ElasticMatrixJob") & (job_table.status == "finished")

# Get corresponding lattice constants
def get_c11(job_path):
    return job_path["output/elasticmatrix"]["C"][0, 0]

def get_c12(job_path):
    return job_path["output/elasticmatrix"]["C"][0, 1]

def get_c44(job_path):
    return job_path["output/elasticmatrix"]["C"][3, 3]
table = pr.create_table("table_elastic", delete_existing_job=True)
table.db_filter_function = filter_elastic
table.add["potential"] = get_potential
table.add["C11"] = get_c11
table.add["C12"] = get_c12
table.add["C44"] = get_c44

table.run()
data_elastic = table.get_dataframe()
data_elastic<jupyter_output>100%|██████████| 5/5 [00:00<00:00, 993.02it/s]
 20%|██        | 1/5 [00:00<00:00,  7.41it/s]<jupyter_text>**Compare computed elastic constants with that from DFT calculations**Comparing the computing elastic constant values with DFT benchmarks from the http://atomistictools.org/ database<jupyter_code># Querying the database http://atomistictools.org/ from using the API interface

from structdbrest import StructDBLightRester
rest = StructDBLightRester(token="<PASSWORD>")

# Querying calculations done using the FHI-AIMS (https://aimsclub.fhi-berlin.mpg.de/) DFT code
fhi_calc = rest.query_calculator_types("FHI%aims%")[0]

# Querying for elastic properties
dft_elast_prop = rest.query_properties(rest.PropertyTypes.ELASTIC_MATRIX, composition="Cu-%", strukturbericht="A1", 
                      calculator_name=fhi_calc.NAME)[0]

# Querying for phonon properties
dft_phon_prop = rest.query_properties(rest.PropertyTypes.PHONONS, composition="Cu-%", strukturbericht="A1", 
                      calculator_name=fhi_calc.NAME)[0]
C_dft = dft_elast_prop.value["C"]
print("DFT C11={:.1f} GPa".format(C_dft[0][0]))
print("DFT C12={:.1f} GPa".format(C_dft[0][1]))
print("DFT C44={:.1f} GPa".format(C_dft[3][3]))<jupyter_output>DFT C11=176.9 GPa
DFT C12=131.7 GPa
DFT C44=82.5 GPa<jupyter_text>**Calculation of surface energies**We can also use the potentials to calculate the surface energies of the Cu (111), (110), and (100) surfaces<jupyter_code>surface_type_list = ["fcc111", "fcc110", "fcc100"]

for i, pot in enumerate(data_murn.potential.to_list()):
    group_name = clean_project_name(pot)
    pr_pot = pr.create_group(group_name)
    a = data_murn.a.to_list()[i]
    for surface_type in surface_type_list:
        surface = pr.create_surface("Cu", surface_type=surface_type, size=(8, 8, 8), a=a, orthogonal=True, vacuum=12)
        job_lammps = pr_pot.create_job(pr_pot.job_type.Lammps, "surf_{}".format(surface_type))
        job_lammps.structure = surface
        job_lammps.potential = get_pyiron_potential(pot)
        job_lammps.calc_minimize()
        job_lammps.run()
# Function to filter only surface calculations
def is_a_surface(job_table):
    return (job_table.hamilton == "Lammps") & (job_table.job.str.contains("fcc")) & (job_table.status == "finished")

def get_potential_lammps_job(job_path):
    return job_path["input/potential/Name"]

def get_surface_type(job_path):
    surf_list =  ["fcc111", "fcc110", "fcc100"]
    conditions = [val in job_path.job_name for val in surf_list]
    return surf_list[np.where(conditions)[0].tolist()[0]]

def get_area(job_path):
    cell = job_path["output/structure/cell/cell"]
    return np.linalg.norm(np.cross(cell[0], cell[1]))
table = pr.create_table("table_surface", delete_existing_job=True)
table.db_filter_function = is_a_surface
table.add["potential"] = get_potential_lammps_job
table.add["surface_type"] = get_surface_type
table.add["surface_area"] = get_area
table.add.get_total_number_of_atoms

table.add.get_energy_tot
table.run()
data_surf = table.get_dataframe()
data_surf<jupyter_output>100%|██████████| 15/15 [00:00<00:00, 1154.69it/s]
 60%|██████    | 9/15 [00:00<00:00, 80.53it/s]<jupyter_text>Compute surface energies using data from the bulk fcc crystal<jupyter_code>eV_A2_to_mJ_m2_factor = 16.0219 * 1e3
data_merged = pd.merge(data_surf, data_murn, on="potential")
data_merged["surface_energy"] = data_merged.energy_tot - (data_merged.eq_energy * data_merged.Number_of_atoms)
data_merged["surface_energy_in_mJ_per_sq_m"] = data_merged.surface_energy / data_merged.surface_area / 2 * eV_A2_to_mJ_m2_factor
data_merged[["potential", "surface_type", "surface_energy_in_mJ_per_sq_m"]]<jupyter_output><empty_output><jupyter_text>Compare to DFT reference data queried from atomistictools.org<jupyter_code>surface_props = rest.query_properties(rest.PropertyTypes.SURFACE_ENERGY, 
                                      composition="Cu-%", 
                                      strukturbericht="A1",
                                      calculator_name=fhi_calc.NAME
                     )
surf_dft_data=[]
for surf_prop in surface_props:
    surf_dft_data.append([surf_prop.VALUE["surface_orientation"], surf_prop.VALUE["number_of_atoms"], surf_prop.VALUE["surface_energy"]])
surf_dft_df=pd.DataFrame(surf_dft_data, columns=["surface_type","number_of_atoms","surface_energy(eV/A^2)"])
surf_dft_df.sort_values(["surface_type", "number_of_atoms"], ascending=[False, False], inplace=True)
surf_dft_df.drop_duplicates(["surface_type"], keep="first", inplace=True)
surf_dft_df["surface_type"] = "fcc" + surf_dft_df["surface_type"]
surf_dft_df["surface_energy_in_mJ_per_sq_m"] = surf_dft_df["surface_energy(eV/A^2)"]*eV_A2_to_mJ_m2_factor
surf_dft_df<jupyter_output><empty_output><jupyter_text>**Finite temperature thermodynamics (Harmonic approximation)**We can also compute free energies and other finite temperature thermodynamic properties within the Harmonic approximation using phonopy (http://phonopy.github.io)<jupyter_code>for pot in data_murn.potential.to_list():
    group_name = clean_project_name(pot)
    pr_pot = pr.create_group(group_name)
    job_id = int(data_murn[data_murn.potential==pot].job_id)
    job_ref = pr_pot.create_job(pr_pot.job_type.Lammps, "ref_job")
    job_ref.structure = pr_pot.inspect(job_id)["output/structure"].to_object()
    job_ref.potential = get_pyiron_potential(pot)
    job_ref.calc_static()
    phonopy_job = job_ref.create_job(pr_pot.job_type.PhonopyJob, "phonopy_job")
    phonopy_job.run()<jupyter_output>The job phonopy_job was saved and received the ID: 238
The job ref_job_0 was saved and received the ID: 239
The job phonopy_job was saved and received the ID: 240
The job ref_job_0 was saved and received the ID: 241
The job phonopy_job was saved and received the ID: 242
The job ref_job_0 was saved and received the ID: 243
The job phonopy_job was saved and received the ID: 244
The job ref_job_0 was saved and received the ID: 245
The job phonopy_job was saved and received the ID: 246
The job ref_job_0 was saved and received the ID: 247<jupyter_text>Plotting the thermodynamic properties using phonopy<jupyter_code>fig, ax_list = plt.subplots(ncols=3, nrows=1, sharex="row")
fig.set_figwidth(20)
#fig.set_figheight(20)


for i, pot in enumerate(potential_list):
    group_name = clean_project_name(pot)
    phonopy_job = pr[group_name+"/phonopy_job"]
    thermo = phonopy_job.get_thermal_properties(t_min=0, t_max=800)
    
    ax = ax_list[0]
    ax.plot(thermo.temperatures, thermo.free_energies, label=pot)
    ax.set_xlabel("Temperatures [K]")
    ax.set_ylabel("Free energies [eV]")
    
    ax = ax_list[1]
    ax.plot(thermo.temperatures, thermo.entropy, label=pot)
    ax.set_xlabel("Temperatures [K]")
    ax.set_ylabel("Entropy [eV/K]")
    
    ax = ax_list[2]
    ax.plot(thermo.temperatures, thermo.cv, label=pot)
    ax.set_xlabel("Temperatures [K]")
    ax.set_ylabel("Heat capacity (C$_\mathrm{V}$) [eV/K]")
    
ax_list[0].legend()
ax_list[1].legend()
ax_list[2].legend()

fig.subplots_adjust(wspace=0.2);<jupyter_output><empty_output><jupyter_text>We can also compute the phonon density of states and compare with that from benchmarked DFT calculations<jupyter_code>fig, ax = plt.subplots(ncols=1, nrows=1, sharex="row")

for i, pot in enumerate(potential_list):
    group_name = clean_project_name(pot)
    phonopy_job = pr[group_name+"/phonopy_job"]
    thermo = phonopy_job.get_thermal_properties(t_min=0, t_max=800)
    phonopy_job.plot_dos(label=pot, ax=ax)
ax.plot(dft_phon_prop.value["dos_energies"],  dft_phon_prop.value["dos_total"], label="DFT")
ax.legend();<jupyter_output><empty_output><jupyter_text>**Validating against original dataset**Finally we see how well the potentials reproduce the energies and forces from the dataset they are fit on<jupyter_code>pr_import = Project("../datasets")
if len(pr_import.job_table()) == 0:
    pr_import.unpack("Cu_training_archive")
pr_import.job_table()
container = pr_import.load("df1_A1_A2_A3_EV_elast_phon")
container.to_pandas()<jupyter_output><empty_output><jupyter_text>Now we iterate over the structures in the dataset and compute the energies and forces for comparison<jupyter_code>structure_list, energy_list, force_list, num_atoms_list =  container.to_list()

energy_per_atom_list = np.array(energy_list) / np.array(num_atoms_list)

energy_pred_dict = dict()
force_pred_dict = dict()

for pot in potential_list:
    group_name = clean_project_name(pot)
    pr_pot = pr.create_group(pot)
    energy_pred_list = list()
    force_pred_list = list()
    
    stride = 10
    for i, struct in enumerate(structure_list[::stride]):
        job = pr_pot.create_job(pr.job_type.Lammps, "lammps_struct_{}".format(i))
        job.potential = get_pyiron_potential(pot)
        job.structure = struct
        job.calc_static()
        job.run()
        energy_pred_list.append(job["output/generic/energy_tot"][-1] / len(struct))
        force_pred_list.append(job["output/generic/forces"][-1])
    job_box = pr_pot.create_job(pr.job_type.Lammps, "lammps_box")
    job_box.potential = get_pyiron_potential(pot)
    job_box.structure = pr_pot.create_atoms("Cu", scaled_positions=[[0.5, 0.5, 0.5]], cell=np.eye(3)*10, pbc=True)
    job_box.calc_static()
    job_box.run()
    if "runner" in pot.lower():
        energy_pred_list = np.array(energy_pred_list)
        
    # correct for energy of isolated atom (no need for NN potentials)
    else:
        energy_pred_list = np.array(energy_pred_list) - job_box["output/generic/energy_tot"][-1]
        
    energy_pred_dict[pot] = energy_pred_list
    force_pred_dict[pot] = force_pred_list
fig, ax_list = plt.subplots(ncols=len(potential_list), nrows=2, sharey="row", sharex="row")
fig.set_figwidth(18)
fig.set_figheight(12)

for i, (pot, energy_pred) in enumerate(energy_pred_dict.items()):
    
    ax = ax_list[0][i]
    ax.plot(energy_per_atom_list[::stride], energy_pred, "x")
    ax.plot(energy_per_atom_list[::stride], energy_per_atom_list[::stride])
    ax.set_title(pot + " (Energies)")
    ax.set_xlabel("Energy DFT [eV]")
    ax.set_ylabel("Energy pred [eV]")
    
    ax = ax_list[1][i]
    force_orig = np.hstack([f.flatten() for f in force_list[::stride]])
    force_pred = np.hstack([f.flatten() for f in force_pred_dict[pot]])
    
    ax.plot(force_orig, force_pred, "x")
    ax.plot(force_orig, force_orig)
    ax.set_xlabel("Force DFT [eV/$\mathrm{\AA}$]")
    ax.set_ylabel("Force pred [eV/$\mathrm{\AA}$]")
    ax.set_title(pot + " (Forces)")
    
fig.subplots_adjust(wspace=0.1, hspace=0.2)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12708372947962224, 'fraction_numerical': 0.022691259428281676, 'mean_word_length': 5.177760375880971, 'pattern_counts': {'":': 0, '<': 51, '<?xml version=': 0, '>': 45, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '564388', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 6057, 'n_tokens_neox': 5658, 'n_words': 1075, 'path': 'day_3/ex_01_validation.ipynb', 'repo_name': 'pyiron/pyiron_potentialf', 'size': '15779'}