Skip to content
GitLab
Unverified Commit 72e92b42 authored by Jan Janssen's avatar Jan Janssen Committed by GitHub
Browse files

Update pacemaker-fit-tutorial.ipynb

parent b5871528
Hide whitespace changes
Inline Side-by-side
day_2/ace/pacemaker-fit-tutorial.ipynb
View file @ 72e92b42
......@@ -32,7 +32,7 @@
"metadata": {},
"outputs": [],
"source": [
"data_pr = Project(\"import_database\")\n",
"data_pr = Project(\"../../datasets\")\n",
"if len(data_pr.job_table()) == 0:\n",
" data_pr.unpack(\"Cu_training_archive\")"
]
......
%% Cell type:code id: tags:
``` python
%pylab inline
```
%%%% Output: stream
Populating the interactive namespace from numpy and matplotlib
%% Cell type:code id: tags:
``` python
from pyiron import Project
```
%% Cell type:code id: tags:
``` python
data_pr = Project("import_database")
data_pr = Project("../../datasets")
if len(data_pr.job_table()) == 0:
data_pr.unpack("Cu_training_archive")
```
%% Cell type:code id: tags:
``` python
data_pr.job_table()
```
%%%% Output: execute_result
id status chemicalformula job \
0 21 finished None df1_A1_A2_A3_EV_elast_phon
1 22 finished None df3_10k
2 23 finished None df2_1k
subjob projectpath \
0 /df1_A1_A2_A3_EV_elast_phon /home/yury/pyiron/projects/
1 /df3_10k /home/yury/pyiron/projects/
2 /df2_1k /home/yury/pyiron/projects/
project timestart \
0 pyiron-2021/import_database/Cu_database/ 2021-02-08 10:33:52.341472
1 pyiron-2021/import_database/Cu_database/ 2021-02-08 10:33:53.993230
2 pyiron-2021/import_database/Cu_database/ 2021-02-08 10:33:54.435308
timestop totalcputime computer hamilton hamversion parentid \
0 None None zora@cmti001#1 GenericJob 0.4 None
1 None None zora@cmti001#1 GenericJob 0.4 None
2 None None zora@cmti001#1 GenericJob 0.4 None
masterid
0 None
1 None
2 None
%% Cell type:code id: tags:
``` python
data_job = data_pr.load('df1_A1_A2_A3_EV_elast_phon')
```
%% Cell type:markdown id: tags:
# Fitting project
%% Cell type:code id: tags:
``` python
from pyiron_gpl.pacemaker.pacemaker import PaceMakerJob
```
%% Cell type:code id: tags:
``` python
fit_pr = Project("pacemaker_fit")
```
%% Cell type:code id: tags:
``` python
#fit_pr.remove_jobs_silently()#
```
%% Cell type:code id: tags:
``` python
job = fit_pr.create_job(job_type=PaceMakerJob, job_name="df1_cut5_pyace" ,delete_existing_job=True)
```
%% Cell type:markdown id: tags:
## Fit
%% Cell type:code id: tags:
``` python
cutoff = 5.0
```
%% Cell type:code id: tags:
``` python
job.input["potential"]= {
"deltaSplineBins": 0.001,
"element": "Cu",
"fs_parameters": [1, 1, 1, 0.5],
"npot": "FinnisSinclairShiftedScaled",
"NameOfCutoffFunction": "cos",
"rankmax": 3,
"nradmax": [7,2,1],
"lmax": [0,2,1],
"ndensity": 2,
"rcut": cutoff,
"dcut": 0.01,
"radparameters": [5.25],
"radbase": "ChebExpCos",
}
```
%% Cell type:code id: tags:
``` python
job.input["fit"]= {
'optimizer': 'BFGS',
'maxiter': 150,
'loss': {
'kappa': 0.5,
'L1_coeffs': 5e-7, # L1-regularization
'L2_coeffs': 5e-7, # L2-regularization
'w1_coeffs': 1,
'w2_coeffs': 1,
#radial smoothness regularization
'w0_rad': 1e-4,
'w1_rad': 1e-4,
'w2_rad': 1e-4,
},
}
```
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
job.input["cutoff"] = cutoff
```
%% Cell type:code id: tags:
``` python
job.structure_data=data_job
```
%% Cell type:code id: tags:
``` python
job.executable
```
%%%% Output: execute_result
'/home/yury/pyiron/resources/pacemaker/bin/run_pacemakerjob_pyace.sh'
%% Cell type:code id: tags:
``` python
#job.executable.version="tf_cpu" # select the CPU-forces tensorflow version
```
%% Cell type:code id: tags:
``` python
job.run()
```
%%%% Output: stream
2021-02-28 01:16:07,353 - root - INFO - structure_data is TrainingContainer
2021-02-28 01:16:07,363 - root - INFO - Saving training structures dataframe into /home/yury/pyiron/projects/pyiron-2021/pacemaker_fit/df1_cut5_pyace_hdf5/df1_cut5_pyace/df_fit.pckl.gzip with pickle protocol = 4, compression = gzip
%%%% Output: stream
The job df1_cut5_pyace was saved and received the ID: 256
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
fit_pr.job_table(full_table=True)
```
%%%% Output: execute_result
id status chemicalformula job subjob \
0 62 finished None df1_cut5 /df1_cut5
1 256 finished None df1_cut5_pyace /df1_cut5_pyace
projectpath project \
0 /home/yury/pyiron/projects/ pyiron-2021/pacemaker_fit/
1 /home/yury/pyiron/projects/ pyiron-2021/pacemaker_fit/
timestart timestop totalcputime \
0 2021-02-26 18:00:20.326459 2021-02-26 18:02:25.068815 124.0
1 2021-02-28 01:16:07.713672 2021-02-28 01:20:43.280739 275.0
computer hamilton hamversion parentid masterid
0 pyiron@dell-inspiron#1 PaceMakerJob 0.1 None None
1 pyiron@dell-inspiron#1 PaceMakerJob 0.1 None None
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
plt.plot(job["output/log/loss"])
plt.yscale('log')
plt.xlabel("# iter")
plt.ylabel("Loss")
```
%%%% Output: execute_result
Text(0, 0.5, 'Loss')
%%%% Output: display_data
%% Cell type:code id: tags:
``` python
plt.plot(job["output/log/rmse_energy"])
plt.yscale('log')
plt.xlabel("# iter")
plt.ylabel("Energy RMSE, meV/atom")
```
%%%% Output: execute_result
Text(0, 0.5, 'Energy RMSE, meV/atom')
%%%% Output: display_data
%% Cell type:code id: tags:
``` python
plt.plot(job["output/log/rmse_forces"])
plt.yscale('log')
plt.xlabel("# iter")
plt.ylabel("Forces RMSE, meV/Ang/structure")
```
%%%% Output: execute_result
Text(0, 0.5, 'Forces RMSE, meV/Ang/structure')
%%%% Output: display_data
%% Cell type:markdown id: tags:
# Overview of the fitted potential internals
%% Cell type:code id: tags:
``` python
from pyace import *
```
%% Cell type:code id: tags:
``` python
final_potential = job.get_final_potential()
final_basis_set = ACEBBasisSet(final_potential)
```
%% Cell type:markdown id: tags:
For single-species potential there is only one **species block** for *Cu*:
%% Cell type:code id: tags:
``` python
len(final_potential.funcspecs_blocks)
```
%%%% Output: execute_result
1
%% Cell type:code id: tags:
``` python
Cu_block = final_potential.funcspecs_blocks[0]
```
%% Cell type:markdown id: tags:
Basic definitions and notations:
* Radial functions: $R_{nl}(r) = \sum_k c_{nlk} g_k(r)$
* Spherical harmonics: $ Y_{lm}(\hat{\pmb{r}}_{ji}) $
* Basis function: $\phi_{\mu_j \mu_i nlm}(\pmb{r}_{ji}) = R_{nl}^{\mu_j \mu_i}(r_{ji}) Y_{lm}(\hat{\pmb{r}}_{ji}) $
* Atomic base (A-functions): $ A_{i \mu n l m} = \sum_j \delta_{\mu \mu_j} \phi_{\mu_j\mu_i nlm}(\pmb{r}_{ji}) $
* Product of atomic base: $ \pmb{A}_{i\pmb{\mu n l m}} = \prod_{t = 1}^{\nu} A_{i \mu_t n_t l_t m_t} $
* Equivariant basis (B-functions): $ {B}_{i\pmb{\mu n l L}} = \sum_{\pmb{m}} \left(
\begin{array}{c}
\pmb{l m} \\
\pmb{L M}
\end{array}
L_R
\right) \pmb{A}_{i\pmb{\mu n l m}} $ ,
where $ \left(\begin{array}{c} \pmb{l m} \\ \pmb{L M}\end{array} L_R\right) $ is *generalized Clebsh-Gordan coefficients*
* Atomic property (densities) $ \rho_i^{(p)} = \sum_{\pmb{\mu n l L}} {c}^{(p)}_{\mu_i\pmb{\mu n l L}} {B}_{i\pmb{\mu n l L}} $
* Atomic energy: $ E_i = F(\rho_i^{(1)}, \dots,\rho_i^{(P)} ) $
%% Cell type:code id: tags:
``` python
```
%% Cell type:markdown id: tags:
radial coefficients $c_{nlk}$:
%% Cell type:code id: tags:
``` python
np.shape(Cu_block.radcoefficients)
```
%%%% Output: execute_result
(2, 3, 7)
%% Cell type:code id: tags:
``` python
```
%% Cell type:markdown id: tags:
Visualize the radial basis functions ($g_k$) and radial functions ($R_{nl}$) and their derivatives:
%% Cell type:code id: tags:
``` python
RadialFunctionsVisualization(final_basis_set).plot()
```
%%%% Output: display_data
%% Cell type:markdown id: tags:
Total number of basis functions
%% Cell type:code id: tags:
``` python
len(Cu_block.funcspecs)
```
%%%% Output: execute_result
18
%% Cell type:markdown id: tags:
List of B-basis functions $ {B}_{i\pmb{\mu n l L}}$: $\pmb{\mu} = $ `elements`, $\pmb{n} = $ `ns`, $\pmb{l} = $ `ls`, $\pmb{L} = $ `LS`) and corresponding coefficients $ {c}^{(p)}_{\mu_i\pmb{\mu n l L}} =$ `coeffs` for two densities
%% Cell type:code id: tags:
``` python
Cu_block.funcspecs
```
%%%% Output: execute_result
[BBasisFunctionSpecification(elements=[Cu,Cu], ns=[1], ls=[0], coeffs=[-0.947664,-0.20673]),
BBasisFunctionSpecification(elements=[Cu,Cu], ns=[2], ls=[0], coeffs=[0.127744,0.0378252]),
BBasisFunctionSpecification(elements=[Cu,Cu], ns=[3], ls=[0], coeffs=[0.395046,0.128133]),
BBasisFunctionSpecification(elements=[Cu,Cu], ns=[4], ls=[0], coeffs=[0.490821,0.17195]),
BBasisFunctionSpecification(elements=[Cu,Cu], ns=[5], ls=[0], coeffs=[0.192457,0.100866]),
BBasisFunctionSpecification(elements=[Cu,Cu], ns=[6], ls=[0], coeffs=[-0.381239,-0.0435031]),
BBasisFunctionSpecification(elements=[Cu,Cu], ns=[7], ls=[0], coeffs=[-0.717072,-0.124286]),
BBasisFunctionSpecification(elements=[Cu,Cu,Cu], ns=[1,1], ls=[0,0], coeffs=[0.176241,0.135483]),
BBasisFunctionSpecification(elements=[Cu,Cu,Cu], ns=[1,1], ls=[1,1], coeffs=[0.0335965,0.00542081]),
BBasisFunctionSpecification(elements=[Cu,Cu,Cu], ns=[1,1], ls=[2,2], coeffs=[0.361527,0.101027]),
BBasisFunctionSpecification(elements=[Cu,Cu,Cu], ns=[2,1], ls=[0,0], coeffs=[0.953892,0.247374]),
BBasisFunctionSpecification(elements=[Cu,Cu,Cu], ns=[2,1], ls=[1,1], coeffs=[-0.000366411,-0.000581538]),
BBasisFunctionSpecification(elements=[Cu,Cu,Cu], ns=[2,1], ls=[2,2], coeffs=[0.0654394,0.0146316]),
BBasisFunctionSpecification(elements=[Cu,Cu,Cu], ns=[2,2], ls=[0,0], coeffs=[0.142427,0.0350483]),
BBasisFunctionSpecification(elements=[Cu,Cu,Cu], ns=[2,2], ls=[1,1], coeffs=[0.00472481,6.57905e-05]),
BBasisFunctionSpecification(elements=[Cu,Cu,Cu], ns=[2,2], ls=[2,2], coeffs=[0.00853967,-0.000426002]),
BBasisFunctionSpecification(elements=[Cu,Cu,Cu,Cu], ns=[1,1,1], ls=[0,0,0], LS=[0], coeffs=[-0.0106058,-0.0227657]),
BBasisFunctionSpecification(elements=[Cu,Cu,Cu,Cu], ns=[1,1,1], ls=[1,1,0], LS=[0], coeffs=[0.309698,0.0822887])]
%% Cell type:markdown id: tags:
# Test fitted potential
%% Cell type:code id: tags:
``` python
test_pr = Project("test_ace_potential")
```
%% Cell type:code id: tags:
``` python
test_pr.remove_jobs_silently()
```
%% Cell type:code id: tags:
``` python
test_pr.job_table()
```
%%%% Output: execute_result
Empty DataFrame
Columns: []
Index: []
%% Cell type:code id: tags:
``` python
cu_ace_potential = job.get_lammps_potential()
```
%% Cell type:code id: tags:
``` python
cu_ace_potential
```
%%%% Output: execute_result
Config \
0 [pair_style pace\n, pair_coeff * * /home/yury/pyiron/projects/pyiron-2021/pacemaker_fit/df1_cut5_pyace_hdf5/df1_cut5_pyace/df1_cut5_pyace.ace Cu\n]
Filename Model Name Species
0 ACE df1_cut5_pyace [Cu]
%% Cell type:markdown id: tags:
## Optimization
%% Cell type:code id: tags:
``` python
lammps_job = test_pr.create.job.Lammps("opt_lammps", delete_existing_job=True)
```
%% Cell type:code id: tags:
``` python
# lammps_job.executable.version="2020.12.24_pace"
```
%% Cell type:code id: tags:
``` python
lammps_job.potential = cu_ace_potential
```
%% Cell type:code id: tags:
``` python
lammps_job.structure = test_pr.create.structure.ase_bulk("Cu","fcc",cubic=True)
```
%% Cell type:code id: tags:
``` python
lammps_job.calc_minimize(pressure=0.0)
```
%% Cell type:code id: tags:
``` python
lammps_job.run()
```
%%%% Output: stream
The job opt_lammps was saved and received the ID: 257
%% Cell type:code id: tags:
``` python
test_pr.job_table()
```
%%%% Output: execute_result
id status chemicalformula job subjob \
0 257 finished Cu4 opt_lammps /opt_lammps
projectpath project \
0 /home/yury/pyiron/projects/ pyiron-2021/test_ace_potential/
timestart timestop totalcputime \
0 2021-02-28 01:20:51.247184 2021-02-28 01:20:52.024520 0.0
computer hamilton hamversion parentid masterid
0 pyiron@dell-inspiron#1 Lammps 0.1 None None
%% Cell type:markdown id: tags:
## Elastic matrix
%% Cell type:code id: tags:
``` python
elmat_job = test_pr.create.job.ElasticMatrixJob("elmat", delete_existing_job=True)
```
%% Cell type:code id: tags:
``` python
elmat_job.input["eps_range"] = 0.05
```
%% Cell type:code id: tags:
``` python
ref_job = test_pr.create.job.Lammps("ref_job", delete_existing_job=True)
```
%% Cell type:code id: tags:
``` python
# ref_job.executable.version="2020.12.24_pace"
```
%% Cell type:code id: tags:
``` python
ref_job.potential = cu_ace_potential
```
%% Cell type:code id: tags:
``` python
ref_job.structure = lammps_job.get_structure()
```
%% Cell type:code id: tags:
``` python
elmat_job.ref_job = ref_job
```
%% Cell type:code id: tags:
``` python
elmat_job.run()
```
%%%% Output: stream
The job elmat was saved and received the ID: 258
The job s_e_0 was saved and received the ID: 259
The job s_01_e_m0_05000 was saved and received the ID: 260
The job s_01_e_m0_02500 was saved and received the ID: 261
The job s_01_e_0_02500 was saved and received the ID: 262
The job s_01_e_0_05000 was saved and received the ID: 263
The job s_08_e_m0_05000 was saved and received the ID: 264
The job s_08_e_m0_02500 was saved and received the ID: 265
The job s_08_e_0_02500 was saved and received the ID: 266
The job s_08_e_0_05000 was saved and received the ID: 267
The job s_23_e_m0_05000 was saved and received the ID: 268
The job s_23_e_m0_02500 was saved and received the ID: 269
The job s_23_e_0_02500 was saved and received the ID: 270
The job s_23_e_0_05000 was saved and received the ID: 271
%% Cell type:code id: tags:
``` python
elast_val=elmat_job["output/elasticmatrix"]
```
%% Cell type:code id: tags:
``` python
ses = np.array(elast_val['strain_energy'])
for se in ses:
plt.plot(se[:,0], se[:,1])
plt.xlabel("Strain")
plt.ylabel("E, eV")
```
%%%% Output: execute_result
Text(0, 0.5, 'E, eV')
%%%% Output: display_data
%% Cell type:code id: tags:
``` python
C=elast_val["C"]
```
%% Cell type:code id: tags:
``` python
print("C11 = {:.1f} GPa".format(C[0,0]))
print("C12 = {:.1f} GPa".format(C[0,1]))
print("C44 = {:.1f} GPa".format(C[3,3]))
```
%%%% Output: stream
C11 = 210.4 GPa
C12 = 120.1 GPa
C44 = 83.0 GPa
%% Cell type:markdown id: tags:
## Phonons
%% Cell type:code id: tags:
``` python
phon_job = test_pr.create.job.PhonopyJob("phon_job", delete_existing_job=True)
ref_job = test_pr.create.job.Lammps("ref_job", delete_existing_job=True)
# ref_job.executable.version="2020.12.24_pace"
ref_job.potential = cu_ace_potential
ref_job.structure = lammps_job.get_structure()
phon_job.input['dos_mesh']=75
phon_job.ref_job = ref_job
phon_job.run()
```
%%%% Output: stream
The job phon_job was saved and received the ID: 272
The job ref_job_0 was saved and received the ID: 273
%% Cell type:markdown id: tags:
# Comparison with some DFT reference data
%% Cell type:code id: tags:
``` python
from structdbrest import StructDBLightRester
```
%% Cell type:code id: tags:
``` python
rest = StructDBLightRester(token="workshop2021")
```
%% Cell type:code id: tags:
``` python
fhi_calc = rest.query_calculator_types("FHI%aims%")[0]
```
%%%% Output: stream
Querying...done
Response successful, size = 7.20 kB, time = 0.27 s
1 entries received
%% Cell type:code id: tags:
``` python
dft_elast_prop = rest.query_properties(rest.PropertyTypes.ELASTIC_MATRIX, composition="Cu-%", strukturbericht="A1",
calculator_name=fhi_calc.NAME)[0]
```
%%%% Output: stream
Querying...done
Response successful, size = 37.78 kB, time = 0.45 s
1 entries received
%% Cell type:code id: tags:
``` python
C_dft=dft_elast_prop.value["C"]
```
%% Cell type:code id: tags:
``` python
print("C11 = {:.1f} GPa / DFT C11={:.1f} GPa".format(C[0,0], C_dft[0][0]))
print("C12 = {:.1f} GPa / DFT C11={:.1f} GPa".format(C[0,1], C_dft[0][1]))
print("C44 = {:.1f} GPa / DFT C11={:.1f} GPa".format(C[3,3], C_dft[3][3]))
```
%%%% Output: stream
C11 = 210.4 GPa / DFT C11=176.9 GPa
C12 = 120.1 GPa / DFT C11=131.7 GPa
C44 = 83.0 GPa / DFT C11=82.5 GPa
%% Cell type:code id: tags:
``` python
dft_phon_prop = rest.query_properties(rest.PropertyTypes.PHONONS, composition="Cu-%", strukturbericht="A1",
calculator_name=fhi_calc.NAME)[0]
```
%%%% Output: stream
Querying...done
Response successful, size = 22.50 kB, time = 0.24 s
1 entries received
%% Cell type:code id: tags:
``` python
plt.plot(phon_job["output/dos_energies"], phon_job["output/dos_total"]/4, label="ACE")
plt.plot(dft_phon_prop.value["dos_energies"], dft_phon_prop.value["dos_total"], label="DFT")
plt.legend()
```
%%%% Output: execute_result
<matplotlib.legend.Legend at 0x7f3664444f10>
%%%% Output: display_data
%% Cell type:code id: tags:
``` python
```
......
Supports Markdown
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