Skip to content

GitLab

Unverified Commit a5e6a0c6 authored by Sudarsan Surendralal's avatar Sudarsan Surendralal Committed by GitHub
Browse files

Delete ex_xx_validation.ipynb

parent c989319a
Hide whitespace changes
Inline Side-by-side
day_3/ex_xx_validation.ipynb deleted 100644 → 0
View file @ c989319a
%% Cell type:markdown id:neural-wisdom tags:
# [**Workflows for atomistic simulations**](http://potentials.rub.de/)
%% Cell type:markdown id:chemical-stocks tags:
## **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 fitting
The aim of this exercise is to make you familiar with:
* Potential validation techniques
%% Cell type:code id:greek-heading tags:
``` python
import numpy as np
%matplotlib inline
import matplotlib.pylab as plt
import pandas as pd
```
%% Cell type:code id:middle-things tags:
``` python
from pyiron import Project
```
%% Cell type:code id:steady-genius tags:
``` python
pr = Project("validation")
#pr.remove_jobs(recursive=True)
```
%% Cell type:code id:contained-maryland tags:
``` python
# The list of potentials to iterate over
potential_list = ['2012--Mendelev-M-I--Cu--LAMMPS--ipr1', '2004--Zhou-X-W--Cu-Ag-Au--LAMMPS--ipr2', '1985--Foiles-S-M--Ni-Cu--LAMMPS--ipr1']
```
%% Cell type:code id:demonstrated-checkout tags:
``` python
# Do Murnaghan, ElasticMatrix job, vac formation energy, binding energy, surface energies, comparison with dataset forces, energies
```
%% Cell type:code id:incorporated-majority tags:
``` python
def clean_project_name(name):
return name.replace("-", "_")
```
%% Cell type:code id:ecological-narrow tags:
``` python
%%time
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 = pot
job_ref.calc_minimize()
murn_job = job_ref.create_job(pr_pot.job_type.Murnaghan, "murn_job")
murn_job.run()
murn_job.plot()
```
%% Cell type:code id:right-enhancement tags:
``` python
murn_job["output/equilibrium_energy"]
```
%% Cell type:code id:matched-string tags:
``` python
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"])
```
%% Cell type:code id:fifteen-history tags:
``` python
%%time
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
```
%% Cell type:code id:stainless-chance tags:
``` python
%%time
for pot in data_murn.potential.to_list():
group_name = clean_project_name(pot)
pr_pot = pr.create_group(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 = pot
job_ref.calc_minimize()
elastic_job = job_ref.create_job(pr_pot.job_type.ElasticMatrixJob, "elastic_job")
elastic_job.run()
```
%% Cell type:code id:signal-island tags:
``` python
plt.matshow(elastic_job["output/elasticmatrix"]["C"]);
```
%% Cell type:code id:photographic-balance tags:
``` python
elastic_job["output/elasticmatrix"]["C"].flatten()
```
%% Cell type:code id:vocal-topic tags:
``` python
def get_elastic_constants(job_path):
return job_path["output/elasticmatrix"]["C"].flatten()
def get_only_elastic(job_table):
return job_table.hamilton == "ElasticMatrixJob"
```
%% Cell type:code id:patent-peripheral tags:
``` python
%%time
table = pr.create_table("table_elastic", delete_existing_job=True)
table.db_filter_function = get_only_elastic
table.add["potential"] = get_potential
table.add["C_list"] = get_elastic_constants
table.run()
data_elastic = table.get_dataframe()
data_elastic
```
%% Cell type:code id:sunset-weather tags:
``` python
# Maybe too complex and intimidating for new python users
elastic_matrix = data_elastic.C_list.apply(pd.Series, index=["c{}{}".format(i, j) for i in range(1, 7) for j in range(1, 7)])
data_elastic_matrix = pd.concat([elastic_matrix, data_elastic], axis=1)
data_elastic_matrix
```
%% Cell type:code id:subsequent-contractor tags:
``` python
data_elastic_matrix_short = data_elastic_matrix[['potential', 'c11', 'c12', 'c13', 'c33', 'c44', 'c66']]
data_elastic_matrix_short
```
%% Cell type:code id:written-sense tags:
``` python
%%time
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 = pot
job_lammps.calc_minimize()
job_lammps.run()
```
%% Cell type:code id:hungry-bermuda tags:
``` python
pr_pot
```
%% Cell type:code id:agricultural-resort tags:
``` python
surface.plot3d()
```
%% Cell type:code id:standard-clinton tags:
``` python
len(surface) / 8
```
%% Cell type:code id:polar-surgeon tags:
``` python
pr_pot.job_table().job.str.contains("fcc")
```
%% Cell type:code id:pressing-indie tags:
``` python
def is_a_surface(job_table):
return (job_table.hamilton == "Lammps") & (job_table.job.str.contains("fcc"))
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]))
```
%% Cell type:code id:civilian-maryland tags:
``` python
get_surface_type(pr_pot.inspect('surf_fcc100'))
```
%% Cell type:code id:embedded-visitor tags:
``` python
%%time
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
```
%% Cell type:code id:changed-smith tags:
``` python
```
%% Cell type:code id:surrounded-channel tags:
``` python
data_merged = pd.merge(data_surf, data_murn, on="potential")
data_merged
```
%% Cell type:code id:under-sodium tags:
``` python
data_merged["surface_energy"] = data_merged.energy_tot - (data_merged.eq_energy * data_merged.Number_of_atoms)
```
%% Cell type:code id:increasing-newsletter tags:
``` python
data_merged["surface_energy_in_mJ_per_sq_m"] = data_merged.surface_energy / data_merged.surface_area / 2 * 16.0219 * 1e3
```
%% Cell type:code id:wooden-response tags:
``` python
data_merged
```
%% Cell type:markdown id:sound-transmission tags:
## Finite temperature thermodynamics (Harmonic approximation)
%% Cell type:code id:disturbed-engagement tags:
``` python
%%time
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 = pot
job_ref.calc_static()
phonopy_job = job_ref.create_job(pr_pot.job_type.PhonopyJob, "phonopy_job")
phonopy_job.run()
```
%% Cell type:code id:developmental-helen tags:
``` python
pr
```
%% Cell type:code id:abstract-combination tags:
``` python
fig, ax_list = plt.subplots(ncols=len(potential_list), nrows=3, sharey="row", sharex="row")
fig.set_figwidth(20)
fig.set_figheight(12)
for i, pot in enumerate(potential_list):
group_name = clean_project_name(pot)
ax = ax_list[0][i]
ax.set_title(pot)
phonopy_job = pr[group_name+"/phonopy_job"]
thermo = phonopy_job.get_thermal_properties(t_min=0, t_max=800)
ax = ax_list[0][i]
ax.plot(thermo.temperatures, thermo.free_energies)
ax.set_xlabel("Temperatures [K]")
ax.set_ylabel("Free energies [eV]")
ax = ax_list[1][i]
ax.plot(thermo.temperatures, thermo.entropy)
ax.set_xlabel("Temperatures [K]")
ax.set_ylabel("Entropy [eV/K]")
ax = ax_list[2][i]
ax.plot(thermo.temperatures, thermo.cv)
ax.set_xlabel("Temperatures [K]")
ax.set_ylabel("Heat capacity (C$_\mathrm{V}$) [eV/K]")
fig.subplots_adjust(wspace=0.05, hspace=0.3);
```
%% Cell type:code id:friendly-physics tags:
``` python
```
%% Cell type:markdown id:trying-royal tags:
## **Validating against datasets**
%% Cell type:code id:consistent-leone tags:
``` python
pr_import = Project("../datasets")
if len(pr_import.job_table()) == 0:
pr_import.unpack("Cu_training_archive")
```
%% Cell type:code id:novel-lambda tags:
``` python
container = pr_import['Cu_database/df1_A1_A2_A3_EV_elast_phon']
```
%% Cell type:code id:homeless-circuit tags:
``` python
training_dataset = container.to_pandas()
```
%% Cell type:code id:smoking-controversy tags:
``` python
training_dataset
```
%% Cell type:code id:subsequent-bacteria tags:
``` python
from pyiron import ase_to_pyiron
```
%% Cell type:code id:perfect-raleigh tags:
``` python
structure_list = training_dataset.atoms.apply(ase_to_pyiron).to_list()
energy_list = training_dataset.energy.to_list()
force_list = training_dataset.forces.to_list()
num_atoms_list = training_dataset.number_of_atoms.to_list()
energy_per_atom_list = np.array(energy_list) / np.array(num_atoms_list)
```
%% Cell type:code id:intermediate-novelty tags:
``` python
%%time
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 = 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 = 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()
# correct for energy of isolated atom
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
```
%% Cell type:code id:premier-plain tags:
``` python
fig, ax_list = plt.subplots(ncols=len(potential_list), nrows=4, sharey="row", sharex="row")
fig.set_figwidth(20)
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_x_orig = np.hstack([f[:, 0] for f in force_list[::stride]])
force_x_pred = np.hstack([f[:, 0] for f in force_pred_list])
ax.plot(force_x_orig, force_x_pred, "x")
ax.plot(force_x_orig, force_x_orig)
ax.set_xlabel("Force DFT [eV/$\mathrm{\AA}$]")
ax.set_ylabel("Force pred [eV/$\mathrm{\AA}$]")
ax.set_title(pot + " (Force-x)")
ax = ax_list[2][i]
force_y_orig = np.hstack([f[:, 1] for f in force_list[::stride]])
force_y_pred = np.hstack([f[:, 1] for f in force_pred_list])
ax.plot(force_y_orig, force_y_pred, "x")
ax.plot(force_y_orig, force_y_orig)
ax.set_xlabel("Force DFT [eV/$\mathrm{\AA}$]")
ax.set_ylabel("Force pred [eV/$\mathrm{\AA}$]")
ax.set_title(pot + " (Force-y)")
ax = ax_list[3][i]
force_z_orig = np.hstack([f[:, 2] for f in force_list[::stride]])
force_z_pred = np.hstack([f[:, 2] for f in force_pred_list])
ax.plot(force_z_orig, force_z_pred, "x")
ax.plot(force_z_orig, force_z_orig)
ax.set_xlabel("Force DFT [eV/$\mathrm{\AA}$]")
ax.set_ylabel("Force pred [eV/$\mathrm{\AA}$]")
ax.set_title(pot + " (Force-z)")
fig.subplots_adjust(wspace=0.05, hspace=0.6)
# ax.set_ylim(x_lim)
```
%% Cell type:code id:significant-blank tags:
``` python
```
%% Cell type:code id:imported-digit tags:
``` python
```
%% Cell type:code id:given-colleague tags:
``` python
```
%% Cell type:code id:spread-width tags:
``` python
```
%% Cell type:code id:grave-kingdom tags:
``` python
```
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment