From 5a7c27adda5bd331fc194e1cf2eef0e741521430 Mon Sep 17 00:00:00 2001 From: Jan Janssen Date: Wed, 1 May 2024 10:35:31 -0600 Subject: [PATCH 1/9] Use Jupyterbook rather than Sphinx --- .ci_support/environment-docs.yml | 1 + .readthedocs.yml | 10 +- docs/Makefile | 20 --- docs/_config.yml | 14 ++ docs/_static/pyiron_logo.ico | Bin 21238 -> 0 bytes docs/_toc.yml | 7 + docs/{source => }/installation.md | 0 docs/make.bat | 35 ---- docs/{source => }/materialproperties.md | 10 +- docs/{_static => pictures}/pyiron-logo.png | Bin docs/{source => }/simulationcodes.md | 0 docs/source/conf.py | 41 ----- docs/source/index.rst | 180 --------------------- docs/{source => }/workflows.md | 6 +- 14 files changed, 39 insertions(+), 285 deletions(-) delete mode 100644 docs/Makefile create mode 100644 docs/_config.yml delete mode 100644 docs/_static/pyiron_logo.ico create mode 100644 docs/_toc.yml rename docs/{source => }/installation.md (100%) delete mode 100644 docs/make.bat rename docs/{source => }/materialproperties.md (99%) rename docs/{_static => pictures}/pyiron-logo.png (100%) rename docs/{source => }/simulationcodes.md (100%) delete mode 100644 docs/source/conf.py delete mode 100644 docs/source/index.rst rename docs/{source => }/workflows.md (99%) diff --git a/.ci_support/environment-docs.yml b/.ci_support/environment-docs.yml index f86cb8c6..1dcc4dc4 100644 --- a/.ci_support/environment-docs.yml +++ b/.ci_support/environment-docs.yml @@ -18,3 +18,4 @@ dependencies: - pandas =2.2.2 - pylammpsmpi =0.2.18 - jinja2 =3.1.3 +- jupyter-book =1.0.0 diff --git a/.readthedocs.yml b/.readthedocs.yml index c2d96ee1..f96f30ba 100644 --- a/.readthedocs.yml +++ b/.readthedocs.yml @@ -9,10 +9,16 @@ build: os: "ubuntu-22.04" tools: python: "mambaforge-22.9" + jobs: + pre_build: + # Generate the Sphinx configuration for this Jupyter Book so it builds. + - "cp README.md docs" + - "jupyter-book config sphinx docs/" # Build documentation in the docs/ directory with Sphinx sphinx: - configuration: docs/source/conf.py + builder: html + fail_on_warning: true # Optionally build your docs in additional formats such as PDF and ePub formats: [] @@ -26,3 +32,5 @@ python: install: - method: pip path: . + extra_requirements: + - sphinx diff --git a/docs/Makefile b/docs/Makefile deleted file mode 100644 index d0c3cbf1..00000000 --- a/docs/Makefile +++ /dev/null @@ -1,20 +0,0 @@ -# Minimal makefile for Sphinx documentation -# - -# You can set these variables from the command line, and also -# from the environment for the first two. -SPHINXOPTS ?= -SPHINXBUILD ?= sphinx-build -SOURCEDIR = source -BUILDDIR = build - -# Put it first so that "make" without argument is like "make help". -help: - @$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) - -.PHONY: help Makefile - -# Catch-all target: route all unknown targets to Sphinx using the new -# "make mode" option. $(O) is meant as a shortcut for $(SPHINXOPTS). -%: Makefile - @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) diff --git a/docs/_config.yml b/docs/_config.yml new file mode 100644 index 00000000..6cff8b60 --- /dev/null +++ b/docs/_config.yml @@ -0,0 +1,14 @@ +title: atomistics +author: Jan Janssen +logo: images/pyiron-logo.png + +execute: + execute_notebooks : off + +repository: + url : https://github.com/pyiron/atomistics + path_to_book : "" + +launch_buttons: + notebook_interface : jupyterlab + binderhub_url : https://mybinder.org \ No newline at end of file diff --git a/docs/_static/pyiron_logo.ico b/docs/_static/pyiron_logo.ico deleted file mode 100644 index f5270a0c8742724905f1bd7947fe7b8c7003f5a4..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 21238 zcmeI44{TLcp2uH-MMGF5hPcKMXjnpsF{baey{&|;4U4QYEFvOQMP-pe>^MXbu>6r$ z2aQ1l#${Q;ILIn8NYHVdu5DOG5gA2g5E&QI09GW3$gs5_md>MRKi~WAd2`?Gdwuue ztFvx=$*;fj=XZX8&iS4H=Qa{4h;)k-7Dm8@k&S&Lks*;tq_Q#*A37ou*-hF|(jbYD z?9zck{l-P3u^Xb%qEN(LKs2tXXfg8Zx`1^kmN$6v(V;KUONC5hSvR1Euw~NT=qN22 z_h>EwyJg&R4Q@G?PWdWh%V)x#{U0Oy1K7oPI=vQh(|8pXk~ccL3^a2_VS)JvzH z%T59VBp*D|G*EvH4rQIJ@?}4L;23i61zsh-3K&lK#4oCg$HVUgd@ch^{!Y{A)L_>Y z)c;Z9vw>&G-vaLD$rs&X;4pMq+aEa&Fbb#vW&k~bF@WRQH!>m58I8;d@Lfo}7jQZ8 zMaZsJJVEP7-?2a)X?f|>@R1+O#|BGAzJAnsH*8!)`!AHuo$zK^){}O*I2*0y z1D$+n4g4qAG(RR!^0xs7zPq8j!k@3ObiA>+c=BWS-~Y~YD_6Gt=GkY>hE=QLOJ>jB zbz@1%?cxzELs}0&zH|mKz>`558z-+9up0TI6a6g7wE5#>vB~dl-)`F5+f(3=@4jnx zz43;*skC&er4@$H5bW3iAMtAifw z`70ahN;KxVuztbyDf0cc2g7ePHmdE01$pKYpHLiKc;+`6Hxu%oo%@gG+_`h1bU54E zYO1cfYH=>|g0w%mCi;|>EqvgYrDomZkC;OrePoV)`l9}*jm4Qg_}1!X2+5e^Xeld=6Ar>#U-Zkp%SxtPL=u7@#86a%)RZl4Y|k* z(%w=YeJuW;f2seVF8=s8FFR)a6A!1zzw1XoS{>w>Db9GVIcYWadn9~!! zsj}-AO*Cz7ZKkEU*(_hM;3)Dh_REELZ=ed=mw;9Pxy~k)O+rK29&~c+fbnZ2JT@`@ ztOc)}@(m{c*GWI6rBxG(RUK2=D9ZaF@ zY#6FHmt*76Km5GZeDN#Cs4h;R&wu|ri;%40X_p`~yK6((m=l_jvZCPAY(D=%Y5O2CfDA1ZgRM8M3u5wDx+~ z^&5BXZp%PZ3B*%r9pX+3ou$u)u}$mX9jW{hdYt~Xvy8TJ*Uo-E(AK4)SHJA-rwbFt zcDoJ_+P;-}MRW4rARn*H=2RV=%}Ni_hKW<|O6pJhCaZ(it(PFXH*iTh9?;4cwt8jO z(0*ChOX>L7H006moptbQA{=TyhOQ?t05~tZt?V$^@H^!}^T9ViIrYqLto7=*zG@ygFg~Tv)$N77YfYrUBp_ni*&bR@+VV(2>jjZW)1!sG zE9<88seHY8;6tw))4j>GcWVc8GJRtwv)Uwe>4UfKz8f6D{8Gdt=e=@pp5qON4>lfkIynMp3?qqTAthA z>+eVWIw_3p8VA)c4*TMI%`@BZ)ExdGeUw!zLq((D-M-W`$G zpR^+$?ThhO9No{z*LTyWZNGKx#}4`W6|x&}$4~GYm%h>2PS`zTVl3;4jjIxNIlIB` zxIn&Vf7t=T?im4iT%vs0al*;*Lv?t*KWy=j|MGzo*uD#B_4>tL(lizdb^(TYQs?VB z6B8T;)F#uV1K~vfciYlEhtQp#dp1TrT}&H(AbQ|5W$<}NwkUoQNcsg~+yBM;IPeZ> znh!<0?tJYHmctJ1fA#>fZ45Sl0M?kVJBamwyWaS979ct1p_kf3bJ!+)o3?f) z65a!-O?97#{(rgmLifJFlcd{uP-CdhDu(!VBY!qPy_G!-G!o`a$32hI9OC+ppFUq< z?AAG8FWTSr=Un)rXNpfM4)mf(`2dpXRHjfaWT#O-=*un$gl^ z!_aqwkH&ku=Mg>*j0bYol`=XTo&{)bIB}EXn5m;inQ7NvYpR9~6Ymy4bH`ldraeE- zrnB-XKMdFp`1RI&+~G_<8$a}r@1G3vBEFWqzS+tv-zd%ptlC6(YG2oq3ao1%}4$=ckP$&lBYLMmxa+cNJF_n z@O{VQldMzXTK6T_w83&V4&6P~*I#cAA2^V_zYN}|2jhni9yBTg|5=t#XBfI+=$O*x zK7Ch_pYt7yUmoQzg#QX~(tbD0c?ft3vO4TlRGz#(ylH8XJwY3TaozoIT)o-_g?0Zv zsl2=ux$e5yua6Cr7wqSL{d7k=8VKz}mTx8uFYR^H^#`pDJIv{%jjFlhj;yw;F806w zep2`Sd-k|xmd~HBw0qL&LweeNPP+Eun}L*lS~?!*lLr4mw6}ZaAz3~@(_UaU{Wdwb zkw5Nd_uP}{)MrJ*#{*0 zLs5~&*p)yHpf>0Mj92?-vse8+9Xz#UiF{f7dVuvNZ!+O~fz9}to!d1=%o0sk+M0h|f*y8gL%+_arN8w&L9>)^h9X5N%3>K9KTw@dqZn7TuU!Bu6^ z)+5iRzSm*@ZFWA9us@ss=bHcZsZ-|1*Ig%Hs0^cHSJI%v3f{-){G|hQmHn!(-#+)8 z#)geuX6`C-e6R9>>pQn!exUxLyP$ua zR%~i-ikU?dV&=gcV`ky_m|1W`%*+Sojf6J#(`YZ_OB^5(?`V2mKR>=bl9k4?Zi58Ojv#MO^5x1weGQfJ)gK~U8%XP zVz~R>I8$GkGkLPwKkFGpzU{XBklmAYZYys-yQB`5+!iwpn>Jku=YdlNb z*>+q5HOCR`u+(X8+`SA3-P(YR>0OxJ#sso?PIsz zRoMH+&TP|N{LpOmP5$7^La}v9I#*a zQC09hA(=jTlB<2-WWefU^M|9OTDSbSNs4PAPX{#FO@WwwPPXqoI8uo&Bl@_{#sL;scuh6ZKzIBwu8~ zR5m}f?TJ3WNIX?XszGCw<_AHic`5~~Fv981XK>jG(-r57{+eWQxL&`Ag?<fZNy8>c^Qib@Pt) zeJ%9;FuX;hy3*V~u{7F3y$(&s+omDE2Q(wlsSvR5g+j@K{!z<&De>mokt4hs(G2XAKE;f2M zNSWP$WBGK3(L-+ocZJd0N3)f39ekFwLDa(;A6>?}gf@9M)0uwC!xs6=RA4o*1=t51 z!9R|HwPx5!J8lHeCr#r^*18p6_l;EIjdVsm)7WwZ-7~)BAa6bprVr`8S+f5~pYQsI z{OQy|7+pOYr;Re5ek;eAzR_Mjt-`8X$@7COYwto%E#xtM!=%x3*VRUQ3m-QBI9KP?a+rXj5 zZus-2H2t^E{<>$(_yF3J^MribnvNIaP9JpG;MXHv7`>9dPTOmrIuRUdevw?tUjpAw z_7e{OtEepNn3B$hZ<~f(^#k?OOmlX|^ z%fa`PfA3YFEd|%zSzgq?yD=+1UtL~tK~gM%_q>S3ZZv+iP-|^`plKRtnrquHFm)H0 zZMCMmiF7x0@N170w8O7yQKTRazvhZaLGwV0^obNSQMi2oaq-v|iTtf_U4(-CXcOUf zF#l>HU-+E(!ml~O|I=-?uDpV}+R8}eqV|HiiUE;GbG)E#ILoA_nu5AZ74KWH?ZS#k zBwhqw#G-aeLDS!0m)Ta(l(4gVb9dPvxu8|Ln#e`15m)Dd@kp&P@yOq`N3ItBdk^1b Z6mN&9j(pYAHjv@JM0_gznf|2w{ukjmYlQ#+ diff --git a/docs/_toc.yml b/docs/_toc.yml new file mode 100644 index 00000000..d314203b --- /dev/null +++ b/docs/_toc.yml @@ -0,0 +1,7 @@ +format: jb-book +root: README +chapters: +- file: installation.md +- file: simulationcodes.md +- file: workflows.md +- file: materialproperties.md \ No newline at end of file diff --git a/docs/source/installation.md b/docs/installation.md similarity index 100% rename from docs/source/installation.md rename to docs/installation.md diff --git a/docs/make.bat b/docs/make.bat deleted file mode 100644 index dc1312ab..00000000 --- a/docs/make.bat +++ /dev/null @@ -1,35 +0,0 @@ -@ECHO OFF - -pushd %~dp0 - -REM Command file for Sphinx documentation - -if "%SPHINXBUILD%" == "" ( - set SPHINXBUILD=sphinx-build -) -set SOURCEDIR=source -set BUILDDIR=build - -%SPHINXBUILD% >NUL 2>NUL -if errorlevel 9009 ( - echo. - echo.The 'sphinx-build' command was not found. Make sure you have Sphinx - echo.installed, then set the SPHINXBUILD environment variable to point - echo.to the full path of the 'sphinx-build' executable. Alternatively you - echo.may add the Sphinx directory to PATH. - echo. - echo.If you don't have Sphinx installed, grab it from - echo.https://www.sphinx-doc.org/ - exit /b 1 -) - -if "%1" == "" goto help - -%SPHINXBUILD% -M %1 %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O% -goto end - -:help -%SPHINXBUILD% -M help %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O% - -:end -popd diff --git a/docs/source/materialproperties.md b/docs/materialproperties.md similarity index 99% rename from docs/source/materialproperties.md rename to docs/materialproperties.md index f41feeea..68f05180 100644 --- a/docs/source/materialproperties.md +++ b/docs/materialproperties.md @@ -528,7 +528,7 @@ plt.ylabel("Temperature (K)") ``` The result is visualized in the following graph: -![Compare Thermal Expansion](../pictures/thermalexpansion.png) +![Compare Thermal Expansion](pictures/thermalexpansion.png) Both the [Moruzzi, V. L. et al.](https://link.aps.org/doi/10.1103/PhysRevB.37.790) and the quantum mechanical version of the quasi-harmonic approach start at a larger equilibrium volume as they include the zero point vibrations, resulting in @@ -724,7 +724,7 @@ plt.xlabel("Temperature (K)") plt.ylabel("Helmholtz Free Energy (eV)") plt.legend() ``` -![Helmholtz Free Energy](../pictures/free_energy_hqh.png) +![Helmholtz Free Energy](pictures/free_energy_hqh.png) ``` plt.title("Entropy") @@ -734,7 +734,7 @@ plt.xlabel("Temperature (K)") plt.ylabel("Entropy (J/K/mol)") plt.legend() ``` -![Entropy](../pictures/entropy_hqh.png) +![Entropy](pictures/entropy_hqh.png) ``` plt.title("Heat Capacity") @@ -745,7 +745,7 @@ plt.xlabel("Temperature (K)") plt.ylabel("Heat Capacity (J/K/mol)") plt.legend() ``` -![Heat Capacity](../pictures/heat_capacity_hqh.png) +![Heat Capacity](pictures/heat_capacity_hqh.png) ### Thermo Dynamic Integration To include the anharmonic contribution to the free energy thermo dynamic integration is used to integrate the internal @@ -845,7 +845,7 @@ for lattice_constant, temperature in zip(lattice_constant_lst, temperature_lst): ``` ### Summary The Helmholtz free energy for all three approximations is compared in the following: -![Heat Capacity](../pictures/thermo_int.png) +![Heat Capacity](pictures/thermo_int.png) ## Phase Diagram One of the goals of the `atomistics` package is to be able to calculate phase diagrams with ab-initio precision. Coming diff --git a/docs/_static/pyiron-logo.png b/docs/pictures/pyiron-logo.png similarity index 100% rename from docs/_static/pyiron-logo.png rename to docs/pictures/pyiron-logo.png diff --git a/docs/source/simulationcodes.md b/docs/simulationcodes.md similarity index 100% rename from docs/source/simulationcodes.md rename to docs/simulationcodes.md diff --git a/docs/source/conf.py b/docs/source/conf.py deleted file mode 100644 index 13da9531..00000000 --- a/docs/source/conf.py +++ /dev/null @@ -1,41 +0,0 @@ -# Configuration file for the Sphinx documentation builder. -# -# For the full list of built-in configuration values, see the documentation: -# https://www.sphinx-doc.org/en/master/usage/configuration.html - -# -- Project information ----------------------------------------------------- -# https://www.sphinx-doc.org/en/master/usage/configuration.html#project-information - -project = 'atomistics' -copyright = '2023, Jan Janssen' -author = 'Jan Janssen' - -# -- General configuration --------------------------------------------------- -# https://www.sphinx-doc.org/en/master/usage/configuration.html#general-configuration - -extensions = ["myst_parser", 'sphinx.ext.autodoc', 'sphinx.ext.napoleon'] - -templates_path = ['_templates'] -exclude_patterns = [] - - - -# -- Options for HTML output ------------------------------------------------- -# https://www.sphinx-doc.org/en/master/usage/configuration.html#options-for-html-output - -try: - import sphinx_rtd_theme - html_theme = 'sphinx_rtd_theme' - html_logo = "../_static/pyiron-logo.png" - html_favicon = "../_static/pyiron_logo.ico" -except ImportError: - html_theme = 'alabaster' - -html_static_path = ['_static'] - - -# -- Generate API documentation ---------------------------------------------- -# https://www.sphinx-doc.org/en/master/man/sphinx-apidoc.html - -from sphinx.ext.apidoc import main -main(['-e', '-o', 'apidoc', '../../atomistics/', '--force']) diff --git a/docs/source/index.rst b/docs/source/index.rst deleted file mode 100644 index 9401e2cc..00000000 --- a/docs/source/index.rst +++ /dev/null @@ -1,180 +0,0 @@ -==================================================================== -atomistics - Interfaces for atomistic simulation codes and workflows -==================================================================== - -:Author: Jan Janssen -:Contact: janssen@mpie.de - -The :code:`atomistics` package consists of two primary components. On the one hand it provides interfaces to atomistic -simulation codes - named :code:`calculators`. The supported simulation codes in alphabetical order are: - -* `Abinit `_ - Plane wave density functional theory -* `EMT `_ - Effective medium theory potential -* `GPAW `_ - Density functional theory Python code based on the projector-augmented wave method -* `LAMMPS `_ - Molecular Dynamics -* `Quantum Espresso `_ - Integrated suite of Open-Source computer codes for electronic-structure calculations -* `Siesta `_ - Electronic structure calculations and ab initio molecular dynamics - -For majority of these simulation codes the :code:`atomistics` package use the `Atomic Simulation Environment `_ -to interface the underlying C/ C++ and Fortran Codes with the Python programming language. Still this approach limits -the functionality of the simulation code to calculating the energy and forces, so by adding custom interfaces the -:code:`atomistics` package can support built-in features of the simulation code like structure optimization and molecular -dynamics. - -On the other hand the :code:`atomistics` package also provides :code:`workflows` to calculate material properties on the -atomistic scales, these include: - -* `Equation of State `_ - to calculate equilibrium properties like the equilibrium energy, equilibrium volume, equilibrium bulk modulus and its pressure derivative. -* `Elastic Matrix `_ - to calculate the elastic constants and elastic moduli. -* `Harmonic and Quasi-harmonic Approximation `_ - to calculate the density of states, vibrational free energy and thermal expansion based on the finite displacements method implemented in `phonopy `_. -* `Molecular Dynamics `_ - to calculate finite temperature properties like thermal expansion including the anharmonic contributions. - -All these :code:`workflows` can be coupled with all the simulation codes implemented in the :code:`atomistics` package. -In contrast to the `Atomic Simulation Environment `_ which provides similar functionality -the focus of the :code:`atomistics` package is not to reimplement existing functionality but rather simplify the process -of coupling existing simulation codes with existing workflows. Here the `phonopy `_ -workflow is a great example to enable the calculation of thermodynamic properties with the harmonic and quasi-harmonic -approximation. - -Example -------- -Use the Equation of State to calculate the equilibrium properties like the equilibrium volume, equilibrium energy, -equilibrium bulk modulus and its derivative using the `GPAW `_ simulation code:: - - from ase.build import bulk - from atomistics.calculators import evaluate_with_ase - from atomistics.workflows import EnergyVolumeCurveWorkflow - from gpaw import GPAW, PW - - workflow = EnergyVolumeCurveWorkflow( - structure=bulk("Al", a=4.05, cubic=True), - num_points=11, - fit_type='polynomial', - fit_order=3, - vol_range=0.05, - axes=['x', 'y', 'z'], - strains=None, - ) - task_dict = workflow.generate_structures() - print(task_dict) - >>> {'calc_energy': OrderedDict([ - >>> (0.95, Atoms(symbols='Al4', pbc=True, cell=[3.9813426685908118, 3.9813426685908118, 3.9813426685908118])), - >>> (0.96, Atoms(symbols='Al4', pbc=True, cell=[3.9952635604153612, 3.9952635604153612, 3.9952635604153612])), - >>> (0.97, Atoms(symbols='Al4', pbc=True, cell=[4.009088111958974, 4.009088111958974, 4.009088111958974])), - >>> (0.98, Atoms(symbols='Al4', pbc=True, cell=[4.022817972936038, 4.022817972936038, 4.022817972936038])), - >>> (0.99, Atoms(symbols='Al4', pbc=True, cell=[4.036454748321015, 4.036454748321015, 4.036454748321015])), - >>> (1.0, Atoms(symbols='Al4', pbc=True, cell=[4.05, 4.05, 4.05])), - >>> (1.01, Atoms(symbols='Al4', pbc=True, cell=[4.063455248345461, 4.063455248345461, 4.063455248345461])), - >>> (1.02, Atoms(symbols='Al4', pbc=True, cell=[4.076821973718458, 4.076821973718458, 4.076821973718458])), - >>> (1.03, Atoms(symbols='Al4', pbc=True, cell=[4.0901016179023415, 4.0901016179023415, 4.0901016179023415])), - >>> (1.04, Atoms(symbols='Al4', pbc=True, cell=[4.1032955854717175, 4.1032955854717175, 4.1032955854717175])), - >>> (1.05, Atoms(symbols='Al4', pbc=True, cell=[4.1164052451001565, 4.1164052451001565, 4.1164052451001565])) - >>> ])} - -In the first step the :code:`EnergyVolumeCurveWorkflow` object is initialized including all the parameters to generate -the strained structures and afterwards fit the resulting energy volume curve. This allows the user to see all relevant -parameters at one place. After the initialization the function :code:`generate_structures()` is called without any -additional parameters. This function returns the task dictionary :code:`task_dict` which includes the tasks which should -be executed by the calculator. In this case the task is to calculate the energy :code:`calc_energy` of the eleven -generated structures. Each structure is labeled by the ratio of compression or elongation. In the second step the -:code:`task_dict` is evaluate with the `GPAW `_ simulation code using the -:code:`evaluate_with_ase()` function:: - - result_dict = evaluate_with_ase( - task_dict=task_dict, - ase_calculator=GPAW( - xc="PBE", - mode=PW(300), - kpts=(3, 3, 3) - ) - ) - print(result_dict) - >>> {'energy': { - >>> 0.95: -14.895378072824752, - >>> 0.96: -14.910819737657118, - >>> 0.97: -14.922307241122466, - >>> 0.98: -14.930392279321056, - >>> 0.99: -14.935048569964911, - >>> 1.0: -14.936666396364169, - >>> 1.01: -14.935212782128556, - >>> 1.02: -14.931045138839849, - >>> 1.03: -14.924165445706581, - >>> 1.04: -14.914703574005678, - >>> 1.05: -14.902774559134226 - >>> }} - -In analogy to the :code:`task_dict` which defines the tasks to be executed by the simulation code the :code:`result_dict` -summarizes the results of the calculations. In this case the energies calculated for the specific strains. By ordering -both the :code:`task_dict` and the :code:`result_dict` with the same labels, the :code:`EnergyVolumeCurveWorkflow` object -is able to match the calculation results to the corresponding structure. Finally, in the third step the :code:`analyse_structures()` -function takes the :code:`result_dict` as an input and fits the Equation of State with the fitting parameters defined in -the first step:: - - fit_dict = workflow.analyse_structures(output_dict=result_dict) - print(fit_dict) - >>> {'poly_fit': array([-9.30297838e-05, 2.19434659e-02, -1.68388816e+00, 2.73605421e+01]), - >>> 'fit_type': 'polynomial', - >>> 'fit_order': 3, - >>> 'volume_eq': 66.44252286131888, - >>> 'energy_eq': -14.93670322204575, - >>> 'bulkmodul_eq': 72.38919826304497, - >>> 'b_prime_eq': 4.45383655040775, - >>> 'least_square_error': 4.432974529908853e-09, - >>> 'volume': [63.10861874999998, 63.77291999999998, ..., 69.75163125000002], - >>> 'energy': [-14.895378072824752, -14.910819737657118, ..., -14.902774559134226] - >>> } - -As a result the equilibrium parameters are returned plus the parameters of the polynomial and the set of volumes and -energies which were fitted to achieve these results. The important step here is that while the interface between the -first and the second as well as between the second and the third step is clearly defined independent of the specific -workflow, the initial parameters for the workflow to initialize the :code:`EnergyVolumeCurveWorkflow` object as well as -the final output of the :code:`fit_dict` are workflow specific. - -Disclaimer ----------- -While we try to develop a stable and reliable software library, the development remains a opensource project under the -BSD 3-Clause License without any warranties:: - - BSD 3-Clause License - - Copyright (c) 2023, Jan Janssen - All rights reserved. - - Redistribution and use in source and binary forms, with or without - modification, are permitted provided that the following conditions are met: - - * Redistributions of source code must retain the above copyright notice, this - list of conditions and the following disclaimer. - - * Redistributions in binary form must reproduce the above copyright notice, - this list of conditions and the following disclaimer in the documentation - and/or other materials provided with the distribution. - - * Neither the name of the copyright holder nor the names of its - contributors may be used to endorse or promote products derived from - this software without specific prior written permission. - - THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" - AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE - DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE - FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR - SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER - CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, - OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE - OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. - - -Documentation -------------- - -.. toctree:: - :maxdepth: 2 - - installation - simulationcodes - workflows - materialproperties - -* :ref:`modindex` \ No newline at end of file diff --git a/docs/source/workflows.md b/docs/workflows.md similarity index 99% rename from docs/source/workflows.md rename to docs/workflows.md index 985af6e3..e3126103 100644 --- a/docs/source/workflows.md +++ b/docs/workflows.md @@ -426,7 +426,7 @@ finite temperature phonon spectrum with the 0K phonon spectrum calulated with th ``` calculation.plot_renormalized_phonon_dispersion_bands() ``` -![finite temperature band_structure](../pictures/lammps_md_phonons.png) +![finite temperature band_structure](pictures/lammps_md_phonons.png) ### Langevin Thermostat In addition to the molecular dynamics implemented in the LAMMPS simulation code, the `atomistics` package also provides @@ -570,13 +570,13 @@ This band structure can also be visualised using the built-in plotting function: ``` workflow.plot_band_structure() ``` -![band_structure](../pictures/phonon_bands_al.png) +![band_structure](pictures/phonon_bands_al.png) Just like the desnsity of states which can be plotted using: ``` workflow.plot_dos() ``` -![density of states](../pictures/phonon_dos_al.png) +![density of states](pictures/phonon_dos_al.png) ### Quasi-harmonic Approximation To include the volume expansion with finite temperature the `atomistics` package implements the `QuasiHarmonicWorkflow`: From fbdc5ab01f852083a67e9284de7e4250f9733465 Mon Sep 17 00:00:00 2001 From: Jan Janssen Date: Wed, 1 May 2024 10:41:28 -0600 Subject: [PATCH 2/9] Update .readthedocs.yml --- .readthedocs.yml | 1 - 1 file changed, 1 deletion(-) diff --git a/.readthedocs.yml b/.readthedocs.yml index f96f30ba..ffe36490 100644 --- a/.readthedocs.yml +++ b/.readthedocs.yml @@ -18,7 +18,6 @@ build: # Build documentation in the docs/ directory with Sphinx sphinx: builder: html - fail_on_warning: true # Optionally build your docs in additional formats such as PDF and ePub formats: [] From 3121dc4529e71914511ff434c31b58d6c1d0686c Mon Sep 17 00:00:00 2001 From: Jan Janssen Date: Wed, 1 May 2024 11:01:31 -0600 Subject: [PATCH 3/9] Use jupyter notebooks rather than markdown --- docs/_toc.yml | 11 +- docs/materialproperties.md | 850 +------------------------------ docs/phasediagram.md | 3 + docs/simulationcodes.md | 189 ------- docs/workflows.md | 690 ------------------------- notebooks/simulation_codes.ipynb | 1 + 6 files changed, 13 insertions(+), 1731 deletions(-) create mode 100644 docs/phasediagram.md delete mode 100644 docs/simulationcodes.md delete mode 100644 docs/workflows.md create mode 100644 notebooks/simulation_codes.ipynb diff --git a/docs/_toc.yml b/docs/_toc.yml index d314203b..66724c71 100644 --- a/docs/_toc.yml +++ b/docs/_toc.yml @@ -2,6 +2,11 @@ format: jb-book root: README chapters: - file: installation.md -- file: simulationcodes.md -- file: workflows.md -- file: materialproperties.md \ No newline at end of file +- file: ../notebooks/simulation_codes.ipynb +- file: ../notebooks/lammps_workflows.ipynb +- file: materialproperties.md + sections: + - file: ../notebooks/bulk_modulus_with_gpaw.ipynb + - file: ../thermal_expansion_with_lammps.ipynb + - file: ../free_energy_calculation.ipynb + - file: phasediagram.md \ No newline at end of file diff --git a/docs/materialproperties.md b/docs/materialproperties.md index 68f05180..0ba7476d 100644 --- a/docs/materialproperties.md +++ b/docs/materialproperties.md @@ -1,852 +1,4 @@ # Materials Properties Demonstrate the calculation of common material properties with the `atomistics` package. The examples use different simulation codes, still the examples are not simulation code specific. It is one of the core features of the `atomistics` -package that all simulation workflow to calculate a specific material property can be executed with all simulation codes. - -## Elastic Properties -Calculate the bulk modulus for Aluminium using the [GPAW](https://wiki.fysik.dtu.dk/gpaw/) DFT code: - -### Equation of State -One way to calculate the bulk modulus is using the Equation of State to calculate the equilibrium properties: -``` -from ase.build import bulk -from atomistics.calculators import evaluate_with_ase -from atomistics.workflows import EnergyVolumeCurveWorkflow -from gpaw import GPAW, PW - -workflow = EnergyVolumeCurveWorkflow( - structure=bulk("Al", a=4.05, cubic=True), - num_points=11, - fit_type='polynomial', - fit_order=3, - vol_range=0.05, - axes=['x', 'y', 'z'], - strains=None, -) -task_dict = workflow.generate_structures() -print(task_dict) ->>> {'calc_energy': OrderedDict([ ->>> (0.95, Atoms(symbols='Al4', pbc=True, cell=[3.9813426685908118, 3.9813426685908118, 3.9813426685908118])), ->>> (0.96, Atoms(symbols='Al4', pbc=True, cell=[3.9952635604153612, 3.9952635604153612, 3.9952635604153612])), ->>> (0.97, Atoms(symbols='Al4', pbc=True, cell=[4.009088111958974, 4.009088111958974, 4.009088111958974])), ->>> (0.98, Atoms(symbols='Al4', pbc=True, cell=[4.022817972936038, 4.022817972936038, 4.022817972936038])), ->>> (0.99, Atoms(symbols='Al4', pbc=True, cell=[4.036454748321015, 4.036454748321015, 4.036454748321015])), ->>> (1.0, Atoms(symbols='Al4', pbc=True, cell=[4.05, 4.05, 4.05])), ->>> (1.01, Atoms(symbols='Al4', pbc=True, cell=[4.063455248345461, 4.063455248345461, 4.063455248345461])), ->>> (1.02, Atoms(symbols='Al4', pbc=True, cell=[4.076821973718458, 4.076821973718458, 4.076821973718458])), ->>> (1.03, Atoms(symbols='Al4', pbc=True, cell=[4.0901016179023415, 4.0901016179023415, 4.0901016179023415])), ->>> (1.04, Atoms(symbols='Al4', pbc=True, cell=[4.1032955854717175, 4.1032955854717175, 4.1032955854717175])), ->>> (1.05, Atoms(symbols='Al4', pbc=True, cell=[4.1164052451001565, 4.1164052451001565, 4.1164052451001565])) ->>> ])} -``` -In the first step the `EnergyVolumeCurveWorkflow` object is initialized including all the parameters to generate -the strained structures and afterwards fit the resulting energy volume curve. This allows the user to see all relevant -parameters at one place. After the initialization the function `generate_structures()` is called without any -additional parameters. This function returns the task dictionary `task_dict` which includes the tasks which should -be executed by the calculator. In this case the task is to calculate the energy `calc_energy` of the eleven generated -structures. Each structure is labeled by the ratio of compression or elongation. In the second step the `task_dict` -is evaluated with the [GPAW](https://wiki.fysik.dtu.dk/gpaw/) simulation code using the `evaluate_with_ase()` function: -``` -result_dict = evaluate_with_ase( - task_dict=task_dict, - ase_calculator=GPAW( - xc="PBE", - mode=PW(300), - kpts=(3, 3, 3) - ) -) -print(result_dict) ->>> {'energy': { ->>> 0.95: -14.895378072824752, ->>> 0.96: -14.910819737657118, ->>> 0.97: -14.922307241122466, ->>> 0.98: -14.930392279321056, ->>> 0.99: -14.935048569964911, ->>> 1.0: -14.936666396364169, ->>> 1.01: -14.935212782128556, ->>> 1.02: -14.931045138839849, ->>> 1.03: -14.924165445706581, ->>> 1.04: -14.914703574005678, ->>> 1.05: -14.902774559134226 ->>> }} -``` -In analogy to the `task_dict` which defines the tasks to be executed by the simulation code the `result_dict` summarizes -the results of the calculations. In this case the energies calculated for the specific strains. By ordering both the -`task_dict` and the `result_dict` with the same labels, the `EnergyVolumeCurveWorkflow` object is able to match the -calculation results to the corresponding structure. Finally, in the third step the `analyse_structures()` function takes -the `result_dict` as an input and fits the Equation of State with the fitting parameters defined in the first step: -``` -fit_dict = workflow.analyse_structures(output_dict=result_dict) -print(fit_dict) ->>> {'poly_fit': array([-9.30297838e-05, 2.19434659e-02, -1.68388816e+00, 2.73605421e+01]), ->>> 'fit_type': 'polynomial', ->>> 'fit_order': 3, ->>> 'volume_eq': 66.44252286131888, ->>> 'energy_eq': -14.93670322204575, ->>> 'bulkmodul_eq': 72.38919826304497, ->>> 'b_prime_eq': 4.45383655040775, ->>> 'least_square_error': 4.432974529908853e-09, ->>> 'volume': [63.10861874999998, 63.77291999999998, ..., 69.75163125000002], ->>> 'energy': [-14.895378072824752, -14.910819737657118, ..., -14.902774559134226] ->>> } -``` -The bulk modulus for Aluminium is calculated using the [GPAW](https://wiki.fysik.dtu.dk/gpaw/) simulation code by fitting -the Equation of State with a third order polynomial over a volume range of +/-5% to be 72.3GPa. - -### Elastic Matrix -An alternative approach to calculate the bulk modulus is based on the relation `B = (1/3) (C11 + 2 C12 )`. The bulk -modulus can be calculated based on the sum of the first elastic constant `C11` and twice the second elastic constant `C12` -divided by there. -``` -from ase.build import bulk -from atomistics.calculators import evaluate_with_ase -from atomistics.workflows import ElasticMatrixWorkflow -from gpaw import GPAW, PW - -workflow = ElasticMatrixWorkflow( - structure=bulk("Al", a=4.05, cubic=True), - num_of_point=5, - eps_range=0.05, - sqrt_eta=True, - fit_order=2 -) -task_dict = workflow.generate_structures() -print(task_dict) ->>> {'calc_energy': OrderedDict([ ->>> ('s_e_0', Atoms(symbols='Al4', pbc=True, cell=[4.05, 4.05, 4.05])), ->>> ('s_01_e_m0_05000', Atoms(symbols='Al4', pbc=True, cell=[3.8421673571095107, 3.8421673571095107, 3.8421673571095107])), ->>> ('s_01_e_m0_02500', Atoms(symbols='Al4', pbc=True, cell=[3.94745170964797, 3.94745170964797, 3.94745170964797])), ->>> ('s_01_e_0_02500', Atoms(symbols='Al4', pbc=True, cell=[4.150015060213919, 4.150015060213919, 4.150015060213919])), ->>> ('s_01_e_0_05000', Atoms(symbols='Al4', pbc=True, cell=[4.247675835085893, 4.247675835085893, 4.247675835085893])), ->>> ('s_08_e_m0_05000', Atoms(symbols='Al4', pbc=True, cell=[3.8421673571095107, 3.8421673571095107, 4.05])), ->>> ('s_08_e_m0_02500', Atoms(symbols='Al4', pbc=True, cell=[3.94745170964797, 3.94745170964797, 4.05])), ->>> ('s_08_e_0_02500', Atoms(symbols='Al4', pbc=True, cell=[4.150015060213919, 4.150015060213919, 4.05])), ->>> ('s_08_e_0_05000', Atoms(symbols='Al4', pbc=True, cell=[4.247675835085893, 4.247675835085893, 4.05])), ->>> ('s_23_e_m0_05000', Atoms(symbols='Al4', pbc=True, cell=[ ->>> [4.039260597921188, -0.2084152371679185, -0.2084152371679185], ->>> [-0.2084152371679185, 4.039260597921188, -0.2084152371679185], ->>> [-0.2084152371679185, -0.2084152371679185, 4.039260597921188] ->>> ])), ->>> ('s_23_e_m0_02500', Atoms(symbols='Al4', pbc=True, cell=[ ->>> [4.047399159178924, -0.1026159010347065, -0.1026159010347065], ->>> [-0.1026159010347065, 4.047399159178924, -0.1026159010347065], ->>> [-0.1026159010347065, -0.1026159010347065, 4.047399159178924] ->>> ])), ->>> ('s_23_e_0_02500', Atoms(symbols='Al4', pbc=True, cell=[ ->>> [4.047526418127057, 0.1000747084794181, 0.1000747084794181], ->>> [0.1000747084794181, 4.047526418127057, 0.1000747084794181], ->>> [0.1000747084794181, 0.1000747084794181, 4.047526418127057] ->>> ])), ->>> ('s_23_e_0_05000', Atoms(symbols='Al4', pbc=True, cell=[ ->>> [4.0402958099962145, 0.19812845289162093, 0.19812845289162093], ->>> [0.19812845289162093, 4.0402958099962145, 0.19812845289162093], ->>> [0.19812845289162093, 0.19812845289162093, 4.0402958099962145] ->>> ])) ->>> ])} -``` -In analogy to the example with the `EnergyVolumeCurveWorkflow` above, the `ElasticMatrixWorkflow` is initialized with all -the parameters required to generate the atomistic structures and afterwards fit the resulting energies. By calling the -`generate_structures()` function the task dictionary `task_dict` is generated. The task dictionary specifies that the -energy should be calculated for a total of thirteen structures with different displacements. In the second step the -structures are again evaluated with the [GPAW](https://wiki.fysik.dtu.dk/gpaw/) simulation code: -``` -result_dict = evaluate_with_ase( - task_dict=task_dict, - ase_calculator=GPAW( - xc="PBE", - mode=PW(300), - kpts=(3, 3, 3) - ) -) -print(result_dict) ->>> {'energy': { ->>> 's_e_0': -14.936666396364958, ->>> 's_01_e_m0_05000': -14.509157650668122, ->>> 's_01_e_m0_02500': -14.841982287144095, ->>> 's_01_e_0_02500': -14.861851384196036, ->>> 's_01_e_0_05000': -14.667794842771894, ->>> 's_08_e_m0_05000': -14.761984597147846, ->>> 's_08_e_m0_02500': -14.915410385310373, ->>> 's_08_e_0_02500': -14.906256779097374, ->>> 's_08_e_0_05000': -14.792358225782438, ->>> 's_23_e_m0_05000': -14.276020694686991, ->>> 's_23_e_m0_02500': -14.82856618064028, ->>> 's_23_e_0_02500': -14.919070452898067, ->>> 's_23_e_0_05000': -14.61301941504334 ->>> }} -``` -The atomistic structures are evaluated with the `evaluate_with_ase()` function, which returns the `result_dict`. This -`result_dict` in analogy to the `task_dict` contains the same keys as well as the energies calculated with the -[GPAW](https://wiki.fysik.dtu.dk/gpaw/) simulation code. Finally, the `result_dict` is provided as an input to the -`analyse_structures()` function to calculate the corresponding elastic constants: -``` -elastic_dict = workflow.analyse_structures(output_dict=result_dict) -print(elastic_dict) ->>> OrderedDict([ ->>> ('SGN', 225), ->>> ('v0', 66.43012500000002), ->>> ('LC', 'CI'), ->>> ('Lag_strain_list', ['01', '08', '23']), ->>> ('epss', array([-0.05 , -0.025, 0. , 0.025, 0.05 ])), ->>> ('e0', -14.936666396364958), ->>> ('strain_energy', [ ->>> [ ->>> (-0.05, -14.509157650668122), ->>> (-0.025, -14.841982287144095), ->>> (0.0, -14.936666396364958), ->>> (0.02500000000000001, -14.861851384196036), ->>> (0.05, -14.667794842771894) ->>> ], ->>> [ ->>> (-0.05, -14.761984597147846), ->>> (-0.025, -14.915410385310373), ->>> (0.0, -14.936666396364958), ->>> (0.02500000000000001, -14.906256779097374), ->>> (0.05, -14.792358225782438) ->>> ], ->>> [ ->>> (-0.05, -14.276020694686991), ->>> (-0.025, -14.82856618064028), ->>> (0.0, -14.936666396364958), ->>> (0.02500000000000001, -14.919070452898067), ->>> (0.05, -14.61301941504334) ->>> ] ->>> ]), ->>> ('C', array([ ->>> [98.43569795, 63.17412931, 63.17412931, 0. , 0. , 0. ], ->>> [63.17412931, 98.43569795, 63.17412931, 0. , 0. , 0. ], ->>> [63.17412931, 63.17412931, 98.43569795, 0. , 0. , 0. ], ->>> [ 0. , 0. , 0. , 84.66136128, 0. , 0. ], ->>> [ 0. , 0. , 0. , 0. , 84.66136128, 0. ], ->>> [ 0. , 0. , 0. , 0. , 0. , 84.66136128] ->>> ])), ->>> ('A2', array([2.10448666, 1.0086892 , 3.17048793])), ->>> ('BV', 74.92798552228444), ->>> ('GV', 57.8491304939761), ->>> ('EV', 138.02584019743034), ->>> ('nuV', 0.19298111327535986), ->>> ('S', array([ ->>> [ 0.02038923, -0.00797026, -0.00797026, 0. , 0. , 0. ], ->>> [-0.00797026, 0.02038923, -0.00797026, 0. , 0. , 0. ], ->>> [-0.00797026, -0.00797026, 0.02038923, 0. , 0. , 0. ], ->>> [ 0. , 0. , 0. , 0.01181176, 0. , 0. ], ->>> [ 0. , 0. , 0. , 0. , 0.01181176, 0. ], ->>> [ 0. , 0. , 0. , 0. , 0. , 0.01181176] ->>> ])), ->>> ('BR', 74.9279855222844), ->>> ('GR', 33.5856196420454), ->>> ('ER', 87.65941305083547), ->>> ('nuR', 0.30501408020913495), ->>> ('BH', 74.92798552228442), ->>> ('GH', 45.71737506801075), ->>> ('EH', 113.97207240565497), ->>> ('nuH', 0.246485304942663), ->>> ('AVR', 26.536421673199147), ->>> ('C_eigval', EigResult( ->>> eigenvalues=array([ 35.26156864, 224.78395657, 35.26156864, 84.66136128, 84.66136128, 84.66136128]), ->>> eigenvectors=array([ ->>> [-0.81649658, 0.57735027, -0.15564171, 0. , 0. , 0. ], ->>> [ 0.40824829, 0.57735027, -0.61632016, 0. , 0. , 0. ], ->>> [ 0.40824829, 0.57735027, 0.77196187, 0. , 0. , 0. ], ->>> [ 0. , 0. , 0. , 1. , 0. , 0. ], ->>> [ 0. , 0. , 0. , 0. , 1. , 0. ], ->>> [ 0. , 0. , 0. , 0. , 0. , 1. ] ->>> ]) ->>> )) ->>> ]) -``` -The bulk modulus calculated from the elastic constants `C11` and `C12` based on a strain of +/- 5% is calculated with -the [GPAW](https://wiki.fysik.dtu.dk/gpaw/) simulation code to be 74.9GPa. This differs from the bulk modulus calculated -from the Equation of State above by 2.6GPa. In comparison to the experimental bulk modulus for Aluminium which is -[reported to be 76GPa](https://periodictable.com/Elements/013/data.html) the calculation based on the elastic constants -seem to be more precise, still this is more likely related to error cancellation. In general elastic properties calculated -from density functional theory are expected to have errors of about 5-10% unless carefully converged. - -## Thermal Expansion -Calculate the thermal expansion for a Morse Pair potential using the [LAMMPS](https://www.lammps.org/) molecular dynamics -simulation code. In the following three methods to calculate the thermal expansion are introduced and compared for a -Morse Pair Potential for Aluminium. - -As a first step the potential is defined for the [LAMMPS](https://www.lammps.org/) molecular dynamics simulation code -by specifying the `pair_style` and `pair_coeff` commands for the [Morse Pair Potential](https://docs.lammps.org/pair_morse.html) -as well as the Aluminium bulk structure: -``` -from ase.build import bulk -import pandas - -potential_dataframe = pandas.DataFrame({ - "Config": [[ - "pair_style morse/smooth/linear 9.0", - "pair_coeff * * 0.5 1.8 2.95" - ]], - "Filename": [[]], - "Model": ["Morse"], - "Name": ["Morse"], - "Species": [["Al"]], -}) - -structure = bulk("Al", a=4.05, cubic=True) -``` -The `pandas.DataFrame` based format to specify interatomic potentials is the same `pylammpsmpi` uses to interface with -the [NIST database for interatomic potentials](https://www.ctcms.nist.gov/potentials). In comparison to just providing -the `pair_style` and `pair_coeff` commands, this extended format enables referencing specific files for the interatomic -potentials `"Filename": [[]],` as well as the atomic species `"Species": [["Al"]],` to enable consistency checks if the -interatomic potential implements all the interactions to simulate a given atomic structure. - -Finally, the last step of the preparation before starting the actual calculation is optimizing the interatomic structure. -While for the Morse potential used in this example this is not necessary, it is essential for extending this example to -other interactomic potentials. For the structure optimization the `optimize_positions_and_volume()` function is imported -and applied on the `ase.atoms.Atoms` bulk structure for Aluminium: -``` -from atomistics.workflows import optimize_positions_and_volume - -task_dict = optimize_positions_and_volume(structure=structure) -task_dict ->>> {'optimize_positions_and_volume': Atoms(symbols='Al4', pbc=True, cell=[4.05, 4.05, 4.05])} -``` -It returns a `task_dict` with a single task, the optimization of the positions and the volume of the Aluminium structure. -This task is executed with the [LAMMPS](https://www.lammps.org/) molecular dynamics simulation code using the -`evaluate_with_lammps()` function: -``` -from atomistics.calculators import evaluate_with_lammps - -result_dict = evaluate_with_lammps( - task_dict=task_dict, - potential_dataframe=potential_dataframe, -) -structure_opt = result_dict["structure_with_optimized_positions_and_volume"] -``` -The `result_dict` just contains a single element, the `ase.atoms.Atoms` structure object with optimized positions and -volume. After this step the preparation is completed and the three different approximations can be compared in the following. - -### Equation of State -The first approximation to calculate the thermal expansion is based on the Equation of State derived by [Moruzzi, V. L. et al.](https://link.aps.org/doi/10.1103/PhysRevB.37.790). -So in analogy to the previous example of calculating the elastic properties from the Equation of State, the `EnergyVolumeCurveWorkflow` -is initialized with the default parameters: -``` -from atomistics.workflows import EnergyVolumeCurveWorkflow - -workflow_ev = EnergyVolumeCurveWorkflow( - structure=structure_opt, - num_points=11, - fit_type='polynomial', - fit_order=3, - vol_range=0.05, - axes=['x', 'y', 'z'], - strains=None, -) -structure_dict = workflow_ev.generate_structures() -print(structure_dict) ->>> {'calc_energy': OrderedDict([ ->>> (0.95, Atoms(symbols='Al4', pbc=True, cell=[3.9813426685908118, 3.9813426685908118, 3.9813426685908118])), ->>> (0.96, Atoms(symbols='Al4', pbc=True, cell=[3.9952635604153612, 3.9952635604153612, 3.9952635604153612])), ->>> (0.97, Atoms(symbols='Al4', pbc=True, cell=[4.009088111958974, 4.009088111958974, 4.009088111958974])), ->>> (0.98, Atoms(symbols='Al4', pbc=True, cell=[4.022817972936038, 4.022817972936038, 4.022817972936038])), ->>> (0.99, Atoms(symbols='Al4', pbc=True, cell=[4.036454748321015, 4.036454748321015, 4.036454748321015])), ->>> (1.0, Atoms(symbols='Al4', pbc=True, cell=[4.05, 4.05, 4.05])), ->>> (1.01, Atoms(symbols='Al4', pbc=True, cell=[4.063455248345461, 4.063455248345461, 4.063455248345461])), ->>> (1.02, Atoms(symbols='Al4', pbc=True, cell=[4.076821973718458, 4.076821973718458, 4.076821973718458])), ->>> (1.03, Atoms(symbols='Al4', pbc=True, cell=[4.0901016179023415, 4.0901016179023415, 4.0901016179023415])), ->>> (1.04, Atoms(symbols='Al4', pbc=True, cell=[4.1032955854717175, 4.1032955854717175, 4.1032955854717175])), ->>> (1.05,Atoms(symbols='Al4', pbc=True, cell=[4.1164052451001565, 4.1164052451001565, 4.1164052451001565])) ->>> ])} -``` -After the initialization the `generate_structures()` function is called to generate the atomistic structures which are -then in the second step evaluated with the [LAMMPS](https://www.lammps.org/) molecular dynamics simulation code to derive -the equilibrium properties: -``` -result_dict = evaluate_with_lammps( - task_dict=structure_dict, - potential_dataframe=potential_dataframe -) -print(result_dict): ->>> {'energy': { ->>> 0.95: -14.619170288727801, ->>> 0.96: -14.664457483479836, ->>> 0.97: -14.697945635153152, ->>> 0.98: -14.720448033206749, ->>> 0.99: -14.732723972540498, ->>> 1.0: -14.73548275794779, ->>> 1.01: -14.729389420395107, ->>> 1.02: -14.715066161138207, ->>> 1.03: -14.693095226824505, ->>> 1.04: -14.664021603093682, ->>> 1.05: -14.628355511604163 ->>> }} -``` -While in the previous example the fit of the energy volume curve was used directly, here the output of the fit, in -particular the derived equilibrium properties are the input for the Debye model as defined by [Moruzzi, V. L. et al.](https://link.aps.org/doi/10.1103/PhysRevB.37.790): -``` -import numpy as np - -workflow_ev.analyse_structures(output_dict=result_dict) -thermal_properties_dict = workflow_ev.get_thermal_properties( - temperatures=np.arange(1, 1500, 50), - output_keys=["temperatures", "volumes"], -) -temperatures_ev, volume_ev = thermal_properties_dict["temperatures"], thermal_properties_dict["volumes"] -``` -The output of the Debye model provides the change of the temperature specific optimal volume `volume_ev` which can be -plotted over the temperature `temperatures_ev` to determine the thermal expansion. - -### Quasi-Harmonic Approximation -While the [Moruzzi, V. L. et al.](https://link.aps.org/doi/10.1103/PhysRevB.37.790) approach based on the Einstein crystal -is limited to a single frequency, the quasi-harmonic model includes the volume dependent free energy. Inside the -`atomistics` package the harmonic and quasi-harmonic model are implemented based on an interface to the [Phonopy](https://phonopy.github.io/phonopy/) -framework. Still the user interface is still structured in the same three steps of (1) generating structures, (2) evaluating -these structures and (3) fitting the corresponding model. Starting with the initialization of the `QuasiHarmonicWorkflow` -which combines the `PhonopyWorkflow` with the `EnergyVolumeCurveWorkflow`: -``` -from atomistics.workflows import QuasiHarmonicWorkflow -from phonopy.units import VaspToTHz - -workflow_qh = QuasiHarmonicWorkflow( - structure=structure_opt, - num_points=11, - vol_range=0.05, - interaction_range=10, - factor=VaspToTHz, - displacement=0.01, - dos_mesh=20, - primitive_matrix=None, - number_of_snapshots=None, -) -structure_dict = workflow_qh.generate_structures() -print(structure_dict) ->>> { ->>> 'calc_energy': OrderedDict([ ->>> (0.95, Atoms(symbols='Al4', pbc=True, cell=[3.9813426685908118, 3.9813426685908118, 3.9813426685908118])), ->>> (0.96, Atoms(symbols='Al4', pbc=True, cell=[3.9952635604153612, 3.9952635604153612, 3.9952635604153612])), ->>> (0.97, Atoms(symbols='Al4', pbc=True, cell=[4.009088111958974, 4.009088111958974, 4.009088111958974])), ->>> (0.98, Atoms(symbols='Al4', pbc=True, cell=[4.022817972936038, 4.022817972936038, 4.022817972936038])), ->>> (0.99, Atoms(symbols='Al4', pbc=True, cell=[4.036454748321015, 4.036454748321015, 4.036454748321015])), ->>> (1.0, Atoms(symbols='Al4', pbc=True, cell=[4.05, 4.05, 4.05])), ->>> (1.01, Atoms(symbols='Al4', pbc=True, cell=[4.063455248345461, 4.063455248345461, 4.063455248345461])), ->>> (1.02, Atoms(symbols='Al4', pbc=True, cell=[4.076821973718458, 4.076821973718458, 4.076821973718458])), ->>> (1.03, Atoms(symbols='Al4', pbc=True, cell=[4.0901016179023415, 4.0901016179023415, 4.0901016179023415])), ->>> (1.04, Atoms(symbols='Al4', pbc=True, cell=[4.1032955854717175, 4.1032955854717175, 4.1032955854717175])), ->>> (1.05,Atoms(symbols='Al4', pbc=True, cell=[4.1164052451001565, 4.1164052451001565, 4.1164052451001565])) ->>> ]), ->>> 'calc_forces': { ->>> (0.95, 0): Atoms(symbols='Al108', pbc=True, cell=[11.944028005772434, 11.944028005772434, 11.944028005772434]), ->>> (0.96, 0): Atoms(symbols='Al108', pbc=True, cell=[11.985790681246083, 11.985790681246083, 11.985790681246083]), ->>> (0.97, 0): Atoms(symbols='Al108', pbc=True, cell=[12.027264335876922, 12.027264335876922, 12.027264335876922]), ->>> (0.98, 0): Atoms(symbols='Al108', pbc=True, cell=[12.068453918808114, 12.068453918808114, 12.068453918808114]), ->>> (0.99, 0): Atoms(symbols='Al108', pbc=True, cell=[12.109364244963045, 12.109364244963045, 12.109364244963045]), ->>> (1.0, 0): Atoms(symbols='Al108', pbc=True, cell=[12.149999999999999, 12.149999999999999, 12.149999999999999]), ->>> (1.01, 0): Atoms(symbols='Al108', pbc=True, cell=[12.190365745036383, 12.190365745036383, 12.190365745036383]), ->>> (1.02, 0): Atoms(symbols='Al108', pbc=True, cell=[12.230465921155373, 12.230465921155373, 12.230465921155373]), ->>> (1.03, 0): Atoms(symbols='Al108', pbc=True, cell=[12.270304853707025, 12.270304853707025, 12.270304853707025]), ->>> (1.04, 0): Atoms(symbols='Al108', pbc=True, cell=[12.309886756415153, 12.309886756415153, 12.309886756415153]), ->>> (1.05, 0): Atoms(symbols='Al108', pbc=True, cell=[12.349215735300469, 12.349215735300469, 12.349215735300469]) ->>> } ->>> } -``` -In contrast to the previous workflows which only used the `calc_energy` function of the simulation codes the `PhonopyWorkflow` -and correspondingly also the `QuasiHarmonicWorkflow` require the calculation of the forces `calc_forces` in addition to -the calculation of the energy. Still the general steps of the workflow remain the same: -``` -result_dict = evaluate_with_lammps( - task_dict=structure_dict, - potential_dataframe=potential_dataframe, -) -``` -The `structure_dict` is evaluated with the [LAMMPS](https://www.lammps.org/) molecular dynamics simulation code to -calculate the corresponding energies and forces. The output is not plotted here as the forces for the 108 atom cells -result in 3x108 outputs per cell. Still the structure of the `result_dict` again follows the labels of the -`structure_dict` as explained before. Finally, in the third step the individual free energy curves at the different -temperatures are fitted to determine the equilibrium volume at the given temperature using the `analyse_structures()` -and `get_thermal_properties()` functions: -``` -workflow_qh.analyse_structures(output_dict=result_dict) -thermal_properties_dict_qm = workflow_qh.get_thermal_properties( - temperatures=np.arange(1, 1500, 50), - output_keys=["temperatures", "volumes"], - quantum_mechanical=True -) -temperatures_qh_qm, volume_qh_qm = thermal_properties_dict_qm["temperatures"], thermal_properties_dict_qm["volumes"] -``` -The optimal volume at the different `temperatures` is stored in the `volume_qh_qm` in analogy to the previous section. -Here the extension `_qm` indicates that the quantum-mechanical harmonic oszillator is used. -``` -thermal_properties_dict_cl = workflow_qh.get_thermal_properties( - temperatures=np.arange(1, 1500, 50), - output_keys=["temperatures", "volumes"], - quantum_mechanical=False, -) -temperatures_qh_cl, volume_qh_cl = thermal_properties_dict_cl["temperatures"], thermal_properties_dict_cl["volumes"] -``` -For the classical harmonic oszillator the resulting volumes are stored as `volume_qh_cl`. - -### Molecular Dynamics -Finally, the third and most commonly used method to determine the volume expansion is using a molecular dynamics -calculation. While the `atomistics` package already includes a `LangevinWorkflow` at this point we use the [Nose-Hoover -thermostat implemented in LAMMPS](https://docs.lammps.org/fix_nh.html) directly via the LAMMPS calculator interface. -``` -from atomistics.calculators import calc_molecular_dynamics_thermal_expansion_with_lammps - -structure_md = structure_opt.repeat(11) -result_dict = calc_molecular_dynamics_thermal_expansion_with_lammps( - structure=structure_md, # atomistic structure - potential_dataframe=potential_dataframe, # interatomic potential defined as pandas.DataFrame - Tstart=15, # temperature to for initial velocity distribution - Tstop=1500, # final temperature - Tstep=5, # temperature step - Tdamp=0.1, # temperature damping of the thermostat - run=100, # number of MD steps for each temperature - thermo=100, # print out from the thermostat - timestep=0.001, # time step for molecular dynamics - Pstart=0.0, # initial pressure - Pstop=0.0, # final pressure - Pdamp=1.0, # barostat damping - seed=4928459, # random seed - dist="gaussian", # Gaussian velocity distribution -) -temperature_md_lst, volume_md_lst = result_dict["temperatures"], result_dict["volumes"] -``` -The `calc_molecular_dynamics_thermal_expansion_with_lammps()` function defines a loop over a vector of temperatures in -5K steps. For each step 100 molecular dynamics steps are executed before the temperature is again increased by 5K. For -~280 steps with the Morse Pair Potential this takes approximately 5 minutes on a single core. These simulations can be -further accelerated by adding the `cores` parameter. The increase in computational cost is on the one hand related to -the large number of force and energy calls and on the other hand to the size of the atomistic structure, as these -simulations are typically executed with >5000 atoms rather than the 4 or 108 atoms in the other approximations. The -volume for the individual temperatures is stored in the `volume_md_lst` list. - -### Summary -To visually compare the thermal expansion predicted by the three different approximations, the [matplotlib](https://matplotlib.org) -is used to plot the volume over the temperature: -``` -import matplotlib.pyplot as plt -plt.plot(np.array(volume_md_lst)/len(structure_md) * len(structure_opt), temperature_md_lst, label="Molecular Dynamics", color="C2") -plt.plot(volume_qh_qm, temperatures_qh_qm, label="Quasi-Harmonic (qm)", color="C3") -plt.plot(volume_qh_cl, temperatures_qh_cl, label="Quasi-Harmonic (classic)", color="C0") -plt.plot(volume_ev, temperatures_ev, label="Moruzzi Model", color="C1") -plt.axvline(structure_opt.get_volume(), linestyle="--", color="red") -plt.legend() -plt.xlabel("Volume ($\AA^3$)") -plt.ylabel("Temperature (K)") -``` -The result is visualized in the following graph: - -![Compare Thermal Expansion](pictures/thermalexpansion.png) - -Both the [Moruzzi, V. L. et al.](https://link.aps.org/doi/10.1103/PhysRevB.37.790) and the quantum mechanical version of -the quasi-harmonic approach start at a larger equilibrium volume as they include the zero point vibrations, resulting in -an over-prediction of the volume expansion with increasing temperature. The equilibrium volume is indicated by the -dashed red line. Finally, the quasi-harmonic approach with the classical harmonic oscillator agrees very well with the -thermal expansion calculated from molecular dynamics for this example of using the Morse Pair Potential. - -## Helmholtz Free Energy -The thermodynamic phase stability can be evaluated by comparing the Helmholtz free energy of different phases. So being -able to calculate the Helmholtz free energy is an essential step towards calculating the full phase diagram. In the -following there approximations are introduced to calculate the Helmholtz free energy, starting with the harmonic -approximation, followed by the quasi-harmonic approximation which also includes the volume expansion and finally -thermodynamic integration is used to quantify the an-harmonic contributions. This addiabative approach to calculate the -Helmholtz free energy is typically used in combination with ab-initio simulation methods like density functional theory, -still here the approach is demonstrated with the [LAMMPS](https://lammps.org/) simulation code and the Morse potential. - -Starting with the import of the required python modules: -``` -from ase.build import bulk -import numpy as np -from atomistics.workflows import optimize_positions_and_volume, LangevinWorkflow, PhonopyWorkflow, QuasiHarmonicWorkflow -from atomistics.calculators import evaluate_with_lammps_library, evaluate_with_lammps, get_potential_by_name, evaluate_with_hessian -from pylammpsmpi import LammpsASELibrary -from phonopy.units import VaspToTHz -from tqdm import tqdm -import matplotlib.pyplot as plt -import pandas -import scipy.constants -``` - -Followed by the selection of the reference element and the definition of the Morse interatomic potential: -``` -structure_bulk = bulk("Al", cubic=True) -df_pot_selected = pandas.DataFrame({ - "Config": [[ - "pair_style morse/smooth/linear 9.0", - "pair_coeff * * 0.5 1.8 2.95" - ]], - "Filename": [[]], - "Model": ["Morse"], - "Name": ["Morse"], - "Species": [["Al"]], -}) -``` - -As a prerequisite the volume and positions of the atomistic structure are optimized to calculate the Helmholtz of the -equilibrium structure: -``` -task_dict = optimize_positions_and_volume(structure=structure_bulk) -result_dict = evaluate_with_lammps( - task_dict=task_dict, - potential_dataframe=df_pot_selected, -) -structure_opt = result_dict["structure_with_optimized_positions_and_volume"] -``` -Finally, the size of the atomistic structure is increased by repeating the structure in all three directions three times: -``` -structure = structure_opt.repeat([3, 3, 3]) -``` -The increased structure is required for all three approximations, for the harmonic and quasi-harmonic approximation it -is required to evaluate the phonons up to an interaction range of 10 Angstrom and for the thermodynamic integration the -larger supercell provides improved thermodynamic stability when calculating molecular dynamics trajectories at finite -temperatures. - -### Harmonic Approximation -The harmonic approximation is calculated using the finite displacement method using the [phonopy](https://phonopy.github.io/phonopy/) -package. In the `atomistics` package the [phonopy](https://phonopy.github.io/phonopy/) workflow is implemented in the -`PhonopyWorkflow` object. Following the typical three step process of generating the structures `generate_structures()` -evaluating them with the [LAMMPS](https://lammps.org/) simulation code using the `evaluate_with_lammps()` function and -analysing the results using the `analyse_structures()` function the thermodynamic properties can be calculated using the -`get_thermal_properties()` function. -``` -workflow_fix = PhonopyWorkflow( - structure=structure_opt, - interaction_range=10, - factor=VaspToTHz, - displacement=0.01, - dos_mesh=20, - primitive_matrix=None, - number_of_snapshots=None, -) -structure_dict = workflow_fix.generate_structures() -result_dict = evaluate_with_lammps( - task_dict=structure_dict, - potential_dataframe=df_pot_selected, -) -workflow_fix.analyse_structures(output_dict=result_dict) -term_base_dict = workflow_fix.get_thermal_properties( - t_min=1, - t_max=1500, - t_step=250, - temperatures=None, - cutoff_frequency=None, - pretend_real=False, - band_indices=None, - is_projection=False, -) -print(term_base_dict) ->>> { ->>> 'temperatures': array([1.000e+00, 2.510e+02, 5.010e+02, 7.510e+02, 1.001e+03, 1.251e+03, 1.501e+03]), ->>> 'free_energy': array([ 0.3054683 , 0.26911103, 0.08807443, -0.2090548 , -0.58871765, -1.03150448, -1.52522672]), ->>> 'entropy': array([5.12294911e-14, 4.12565469e+01, 9.50730111e+01, 1.32185181e+02, 1.59664975e+02, 1.81349904e+02, 1.99221616e+02]), ->>> 'heat_capacity': array([ nan, 63.97199066, 88.27044589, 94.37971636, 96.67969872, 97.77513519, 98.37869945]) ->>> } -``` -The dictionary returned by the `get_thermal_properties()` function contains the temperature dependence of the following -thermodynamic properties: Helmholtz free energy `free_energy`, Entropy `entropy` and the heat capcity at constant volume -`heat_capacity`. In addition, the temperature is also included in the result dictionary to simplify the plotting of the -results. The results are compared to the other approximations in the following. - -### Quasi-Harmonic Approximation -The quasi-harmonic approximation extends the harmonix approximation by including the volume expansion. In the -`atomistics` package this is implemented with the `QuasiHarmonicWorkflow` which internally is a combination of the -`EnergyVolumeCurveWorkflow` and the `PhonopyWorkflow` introduced above. Consequently, the input parameters for the -`QuasiHarmonicWorkflow` are a superset of the inputs of these workflows: -``` -workflow_qh = QuasiHarmonicWorkflow( - structure=structure_opt, - num_points=11, - vol_range=0.05, - interaction_range=10, - factor=VaspToTHz, - displacement=0.01, - dos_mesh=20, - primitive_matrix=None, - number_of_snapshots=None, -) -structure_dict = workflow_qh.generate_structures() -result_dict = evaluate_with_lammps( - task_dict=structure_dict, - potential_dataframe=df_pot_selected, -) -workflow_qh.analyse_structures(output_dict=result_dict) -term_qh_dict = workflow_qh.get_thermal_properties( - t_min=1, - t_max=1500, - t_step=250, - temperatures=None, - cutoff_frequency=None, - pretend_real=False, - band_indices=None, - is_projection=False, - quantum_mechanical=True, -) -print(term_qh_dict) ->>> { ->>> 'free_energy': array([ 0.30135185, 0.26231475, 0.07072364, -0.24471851, -0.65110999, -1.12998338, -1.67042388]), ->>> 'entropy': array([3.31166059e-07, 4.25802580e+01, 9.77768688e+01, 1.36322159e+02, 1.65329534e+02, 1.88653045e+02, 2.08300252e+02]), ->>> 'heat_capacity': array([ nan, 64.86754373, 88.84137042, 94.80149704, 97.01828281, 98.06138269, 98.62976854]), ->>> 'volumes': array([66.78514133, 66.90026195, 67.24614141, 67.66763429, 68.13272595, 68.63838791, 69.19002682]), ->>> 'temperatures': array([1.000e+00, 2.510e+02, 5.010e+02, 7.510e+02, 1.001e+03, 1.251e+03, 1.501e+03]) ->>> } -``` -In analogy to the `get_thermal_properties()` function of the `PhonopyWorkflow` the `get_thermal_properties()` function -of the `QuasiHarmonicWorkflow` returns a dictionary with the temperature dependence of the thermodynamic properties. The -volume at finite temperature is included as an additional output. - -Furthermore, the `QuasiHarmonicWorkflow` adds the ability to calculate the volume expansion with the classical harmonic -oscillator rather than the quantum mechanical oscillator which is used by default: -``` -term_qh_cl_dict = workflow_qh.get_thermal_properties( - t_min=1, - t_max=1500, - t_step=250, - temperatures=None, - cutoff_frequency=None, - pretend_real=False, - band_indices=None, - is_projection=False, - quantum_mechanical=False, -) -print(term_qh_cl_dict) ->>> { ->>> 'free_energy': array([ 0.00653974, 0.20391846, 0.04183779, -0.26289749, -0.66376056, -1.13921853]), ->>> 'volumes': array([66.29913676, 66.70960432, 67.13998119, 67.59365054, 68.07505673, 68.5902383 ]), ->>> 'temperatures': array([ 1, 251, 501, 751, 1001, 1251]) ->>> } -``` -In this case the dictionary of thermal properties calculated by the `get_thermal_properties()` function only contains -the Helmholtz free energy `free_energy`, the volume at finite temperature `volumes` and the corresponding temperature -`termpatures`. - -For the comparison of the harmonic approximation and the quasi-harmonic approximation the thermodynamic properties are -visualized in the following. Given the coarse sampling of the temperature dependence these graphs look a bit rough. A -better choice for the temperature step would be 50K rather than 250K. Here the larger temperature step is chosen to -minimize the computational cost. -``` -plt.title("Helmholtz Free Energy") -plt.plot(term_base_dict['temperatures'], term_base_dict['free_energy'], label="Harmonic") -plt.plot(term_qh_dict['temperatures'], term_qh_dict['free_energy'], label="Quasi-Harmonic (qm)") -plt.plot(term_qh_cl_dict['temperatures'], term_qh_cl_dict['free_energy'], label="Quasi-Harmonic (cl)") -plt.xlabel("Temperature (K)") -plt.ylabel("Helmholtz Free Energy (eV)") -plt.legend() -``` -![Helmholtz Free Energy](pictures/free_energy_hqh.png) - -``` -plt.title("Entropy") -plt.plot(term_base_dict['temperatures'], term_base_dict['entropy'], label="harmonic") -plt.plot(term_qh_dict['temperatures'], term_qh_dict['entropy'], label="quasi-harmonic (qm)") -plt.xlabel("Temperature (K)") -plt.ylabel("Entropy (J/K/mol)") -plt.legend() -``` -![Entropy](pictures/entropy_hqh.png) - -``` -plt.title("Heat Capacity") -plt.plot(term_base_dict['temperatures'], term_base_dict['heat_capacity'], label="harmonic") -plt.plot(term_qh_dict['temperatures'], term_qh_dict['heat_capacity'], label="quasi-harmonic") -plt.axhline(3 * scipy.constants.Boltzmann * scipy.constants.Avogadro * len(structure_opt), color="black", linestyle="--", label="$3Nk_B$") -plt.xlabel("Temperature (K)") -plt.ylabel("Heat Capacity (J/K/mol)") -plt.legend() -``` -![Heat Capacity](pictures/heat_capacity_hqh.png) - -### Thermo Dynamic Integration -To include the anharmonic contribution to the free energy thermo dynamic integration is used to integrate the internal -energy between the quasi-harmonic reference and the full anharmonic molecular dynamics calculation. For the choice of -computational efficiency molecular dynamics trajectories with 3000 steps are using and a set of five lambda values -between 0.0 and 1.0. Again, increasing the number of lambda values from 5 to 11 can improve the approximation. -``` -steps = 3000 -steps_lst = list(range(steps)) -lambda_lst = np.linspace(0.0, 1.0, 5) -``` -From the finite temperature volume of the `QuasiHarmonicWorkflow` above the finite temperature lattice constant is -calculated: -``` -lattice_constant_lst = np.array(term_qh_dict['volumes']) ** (1/3) -temperature_lst = term_qh_dict['temperatures'] -``` -The thermodynamic integration workflow consists of two loops. The first loop iterates over the different temperatures -and the corresponding finite temperature lattice constants while the second inner loop iterates over the different -lambda parameters. In the outter loop the `PhonopyWorkflow` workflow is used to calculate the finite temperature force -constants which can be accessed from the `PhonopyWorkflow` using the `get_hesse_matrix()` function. In addition, the -`evaluate_with_lammps()` function is used in the outter loop as well to calculate the equilibrium energy at finite -volume. Finally, in the inner loop the `LangevinWorkflow` is used to calculate the molecular dynamics trajectory with -both the forces and the energy being mixed in dependence of the lambda parameter from the Hessian calculation -`evaluate_with_hessian()` and the [LAMMPS](https://lammps.org/) calculation `evaluate_with_lammps_library()`. Here the -[LAMMPS](https://lammps.org/) library interface is used to evaluate multiple structures one after another. Finally, the -lambda dependence is fitted and the integral from 0 to 1 is calculated: -``` -free_energy_lst, eng_lambda_dependence_lst = [], [] -for lattice_constant, temperature in zip(lattice_constant_lst, temperature_lst): - structure = bulk("Al", a=lattice_constant, cubic=True).repeat([3, 3, 3]) - equilibrium_lammps = evaluate_with_lammps(task_dict={"calc_energy": structure}, potential_dataframe=df_pot_selected)['energy'] - workflow_fix = PhonopyWorkflow( - structure=structure, - interaction_range=10, - factor=VaspToTHz, - displacement=0.01, - dos_mesh=20, - primitive_matrix=None, - number_of_snapshots=None, - ) - structure_dict = workflow_fix.generate_structures() - result_dict = evaluate_with_lammps( - task_dict=structure_dict, - potential_dataframe=df_pot_selected, - ) - workflow_fix.analyse_structures(output_dict=result_dict) - energy_pot_all_lst, energy_mean_lst, energy_kin_all_lst = [], [], [] - for lambda_parameter in lambda_lst: - thermo_eng_pot_lst, thermo_eng_kin_lst = [], [] - workflow_md_thermo = LangevinWorkflow( - structure=structure, - temperature=temperature, - overheat_fraction=2.0, - damping_timescale=100.0, - time_step=1, - ) - lmp = LammpsASELibrary( - working_directory=None, - cores=1, - comm=None, - logger=None, - log_file=None, - library=None, - diable_log_file=True, - ) - for i in tqdm(steps_lst): - task_dict = workflow_md_thermo.generate_structures() - hessian_dict = evaluate_with_hessian( - task_dict=task_dict, - structure_equilibrium=structure, - force_constants=workflow_fix.get_hesse_matrix(), - bulk_modulus=0, - shear_modulus=0, - ) - lammps_dict = evaluate_with_lammps_library( - task_dict=task_dict, - potential_dataframe=df_pot_selected, - lmp=lmp, - ) - result_dict = { - "forces": {0: (1-lambda_parameter) * hessian_dict["forces"][0] + lambda_parameter * lammps_dict["forces"][0]}, - "energy": {0: (1-lambda_parameter) * hessian_dict["energy"][0] + lambda_parameter * (lammps_dict["energy"][0] - equilibrium_lammps)}, - } - eng_pot, eng_kin = workflow_md_thermo.analyse_structures(output_dict=result_dict) - thermo_eng_pot_lst.append(eng_pot) - thermo_eng_kin_lst.append(eng_kin) - lmp.close() - thermo_energy = np.array(thermo_eng_pot_lst) + np.array(thermo_eng_kin_lst) - energy_mean_lst.append(np.mean(thermo_energy[1000:])) - energy_pot_all_lst.append(thermo_eng_pot_lst) - energy_kin_all_lst.append(thermo_eng_kin_lst) - eng_lambda_dependence_lst.append(np.array(energy_mean_lst) / len(structure) * 1000) - fit = np.poly1d(np.polyfit(lambda_lst, np.array(energy_mean_lst) / len(structure) * 1000, 3)) - integral = np.polyint(fit) - free_energy_lst.append(integral(1.0) - integral(0.0)) -``` -### Summary -The Helmholtz free energy for all three approximations is compared in the following: -![Heat Capacity](pictures/thermo_int.png) - -## Phase Diagram -One of the goals of the `atomistics` package is to be able to calculate phase diagrams with ab-initio precision. Coming -soon. +package that all simulation workflow to calculate a specific material property can be executed with all simulation codes. \ No newline at end of file diff --git a/docs/phasediagram.md b/docs/phasediagram.md new file mode 100644 index 00000000..20622eef --- /dev/null +++ b/docs/phasediagram.md @@ -0,0 +1,3 @@ +# Phase Diagram +One of the goals of the `atomistics` package is to be able to calculate phase diagrams with ab-initio precision. Coming +soon. \ No newline at end of file diff --git a/docs/simulationcodes.md b/docs/simulationcodes.md deleted file mode 100644 index b86e1e49..00000000 --- a/docs/simulationcodes.md +++ /dev/null @@ -1,189 +0,0 @@ -# Simulation Codes - -## ASE -At the current stage the majority of simulation codes are interfaced using the [Atomic Simulation Environment (ASE)](https://wiki.fysik.dtu.dk/ase/). -The limitation of the ASE based interfaces is that the simulation codes are only used to calculate energies, forces and -stresses, while more complex computations like structure optimization or molecular dynamics are implemented in python. - -### Abinit -[Abinit](https://www.abinit.org) - Plane wave density functional theory: -``` -from ase.calculators.abinit import Abinit -from ase.units import Ry -from atomistics.calculators import evaluate_with_ase - -result_dict = evaluate_with_ase( - task_dict={}, - ase_calculator=Abinit( - label='abinit_evcurve', - nbands=32, - ecut=10 * Ry, - kpts=(3, 3, 3), - toldfe=1.0e-2, - v8_legacy_format=False, - ) -) -``` -The full documentation of the corresponding interface is available on the [Atomic Simulation Environment](https://wiki.fysik.dtu.dk/ase/ase/calculators/abinit.html) -website. - -### EMT -[EMT](https://wiki.fysik.dtu.dk/ase/ase/calculators/emt.html) - Effective medium theory: -``` -from ase.calculators.emt import EMT -from atomistics.calculators import evaluate_with_ase - -result_dict = evaluate_with_ase( - task_dict={}, - ase_calculator=EMT() -) -``` -The full documentation of the corresponding interface is available on the [Atomic Simulation Environment](https://wiki.fysik.dtu.dk/ase/ase/calculators/emt.html) -website. - -### GPAW -[GPAW](https://wiki.fysik.dtu.dk/gpaw/) - Density functional theory Python code based on the projector-augmented wave -method: -``` -from gpaw import GPAW, PW -from atomistics.calculators import evaluate_with_ase - -result_dict = evaluate_with_ase( - task_dict={}, - ase_calculator=GPAW( - xc="PBE", - mode=PW(300), - kpts=(3, 3, 3) - ) -) -``` -The full documentation of the corresponding interface is available on the [GPAW](https://wiki.fysik.dtu.dk/gpaw/) -website. - -### Quantum Espresso -[Quantum Espresso](https://www.quantum-espresso.org) - Integrated suite of Open-Source computer codes for -electronic-structure calculations: -``` -from ase.calculators.espresso import Espresso -from atomistics.calculators import evaluate_with_ase - -result_dict = evaluate_with_ase( - task_dict={}, - ase_calculator=Espresso( - pseudopotentials={"Al": "Al.pbe-n-kjpaw_psl.1.0.0.UPF"}, - tstress=True, - tprnfor=True, - kpts=(3, 3, 3), - ) -) -``` -The full documentation of the corresponding interface is available on the [Atomic Simulation Environment](https://wiki.fysik.dtu.dk/ase/ase/calculators/espresso.html) -website. - -### Siesta -[Siesta](https://siesta-project.org) - Electronic structure calculations and ab initio molecular dynamics: -``` -from ase.calculators.siesta import Siesta -from ase.units import Ry -from atomistics.calculators import evaluate_with_ase - -result_dict = evaluate_with_ase( - task_dict={}, - ase_calculator=Siesta( - label="siesta", - xc="PBE", - mesh_cutoff=200 * Ry, - energy_shift=0.01 * Ry, - basis_set="DZ", - kpts=(5, 5, 5), - fdf_arguments={"DM.MixingWeight": 0.1, "MaxSCFIterations": 100}, - pseudo_path=os.path.abspath("tests/static/siesta"), - pseudo_qualifier="", - ) -) -``` -The full documentation of the corresponding interface is available on the [Atomic Simulation Environment](https://wiki.fysik.dtu.dk/ase/ase/calculators/siesta.html) -website. - -## LAMMPS -[LAMMPS](https://www.lammps.org) - Molecular Dynamics: -``` -from ase.build import bulk -from atomistics.calculators import evaluate_with_lammps, get_potential_by_name - -structure = bulk("Al", cubic=True) -potential_dataframe = get_potential_by_name( - potential_name='1999--Mishin-Y--Al--LAMMPS--ipr1' -) - -result_dict = evaluate_with_lammps( - task_dict={}, - potential_dataframe=potential_dataframe, -) -``` -The [LAMMPS](https://www.lammps.org) interface is based on the [pylammpsmpi](https://github.com/pyiron/pylammpsmpi) -package which couples a [LAMMPS](https://www.lammps.org) instance which is parallelized via the Message Passing Interface -(MPI) with a serial python process or jupyter notebook. The challenging part about molecular dynamics simulation is -identifying a suitable interatomic potential. - -To address this challenge the `atomistics` package is leveraging the [NIST database of interatomic potentials](https://www.ctcms.nist.gov/potentials). -It is recommended to install this database `iprpy-data` via the `conda` package manager, then the `resource_path` is -automatically set to `${CONDA_PREFIX}/share/iprpy`. Alternatively, the `resource_path` can be specified manually as an -optional parameter of the `get_potential_by_name()` function. - -In addition, the `get_potential_dataframe(structure)` function which takes an `ase.atoms.Atoms` object as input can be -used to query the [NIST database of interatomic potentials](https://www.ctcms.nist.gov/potentials) for potentials, which -include the interatomic interactions required to simulate the atomic structure defined by the `ase.atoms.Atoms` object. -It returns a `pandas.DataFrame` with all the available potentials and the `resource_path` can again be specified as -optional parameter. - -Finally, another option to specify the interatomic potential for a LAMMPS simulation is by defining the `potential_dataframe` -directly: -``` -potential_dataframe = pandas.DataFrame({ - "Config": [[ - "pair_style morse/smooth/linear 9.0", - "pair_coeff * * 0.5 1.8 2.95" - ]], - "Filename": [[]], - "Model": ["Morse"], - "Name": ["Morse"], - "Species": [["Al"]], -}) -``` - -## Quantum Espresso -[Quantum Espresso](https://www.quantum-espresso.org) - Integrated suite of Open-Source computer codes for -electronic-structure calculations: -``` -from atomistics.calculators import evaluate_with_qe - -result_dict = evaluate_with_qe( - task_dict={}, - calculation_name="espresso", - working_directory=".", - kpts=(3, 3, 3), - pseudopotentials={ - "Al": "Al.pbe-n-kjpaw_psl.1.0.0.UPF" - }, - tstress=True, - tprnfor=True, - ecutwfc=40.0, # kinetic energy cutoff (Ry) for wavefunctions - conv_thr=1e-06, # Convergence threshold for selfconsistency - diagonalization='david', - electron_maxstep=100, # maximum number of iterations in a scf step. - nstep=200, # number of molecular-dynamics or structural optimization steps performed in this run. - etot_conv_thr=1e-4, # Convergence threshold on total energy (a.u) for ionic minimization - forc_conv_thr=1e-3, # Convergence threshold on forces (a.u) for ionic minimization - smearing='gaussian', # ordinary Gaussian spreading (Default) -) -``` -This secondary interface for [Quantum Espresso](https://www.quantum-espresso.org) is based on the input writer from the -[Atomic Simulation Environment (ASE)](https://wiki.fysik.dtu.dk/ase/) and the output is parsed using the [pwtools](https://elcorto.github.io/pwtools/). -The executable can be set using the `ASE_ESPRESSO_COMMAND` environment variable: -``` -export ASE_ESPRESSO_COMMAND="pw.x -in PREFIX.pwi > PREFIX.pwo" -``` -The full list of possible keyword arguments is available in the [Quantum Espresso Documentation](https://www.quantum-espresso.org/Doc/INPUT_PW.html). -Finally, the [Standard solid-state pseudopotentials (SSSP)](https://www.materialscloud.org/discover/sssp/table/efficiency) -for quantum espresso are distributed via the materials cloud. \ No newline at end of file diff --git a/docs/workflows.md b/docs/workflows.md deleted file mode 100644 index e3126103..00000000 --- a/docs/workflows.md +++ /dev/null @@ -1,690 +0,0 @@ -# Workflows -To demonstrate the workflows implemented in the `atomistics` package, the [LAMMPS](https://www.lammps.org/) molecular -dynamics simulation code is used in the following demonstrations. Still the same `workflows` can also be used with other -simulation codes: -``` -from atomistics.calculators import evaluate_with_lammps, get_potential_by_name - -potential_dataframe = get_potential_by_name( - potential_name='1999--Mishin-Y--Al--LAMMPS--ipr1' -) -result_dict = evaluate_with_lammps( - task_dict={}, - potential_dataframe=potential_dataframe, -) -``` -The interatomic potential for Aluminium from Mishin named `1999--Mishin-Y--Al--LAMMPS--ipr1` is used in the evaluation -with [LAMMPS](https://www.lammps.org/) `evaluate_with_lammps()`. - -## Elastic Matrix -The elastic constants and elastic moduli can be calculated using the `ElasticMatrixWorkflow`: -``` -from ase.build import bulk -from atomistics.calculators import evaluate_with_lammps, get_potential_by_name -from atomistics.workflows import ElasticMatrixWorkflow - -potential_dataframe = get_potential_by_name( - potential_name='1999--Mishin-Y--Al--LAMMPS--ipr1' -) -workflow = ElasticMatrixWorkflow( - structure=bulk("Al", cubic=True), - num_of_point=5, - eps_range=0.005, - sqrt_eta=True, - fit_order=2, -) -task_dict = workflow.generate_structures() -result_dict = evaluate_with_lammps( - task_dict=task_dict, - potential_dataframe=potential_dataframe, -) -fit_dict = workflow.analyse_structures(output_dict=result_dict) -print(fit_dict) -``` -The `ElasticMatrixWorkflow` takes an `ase.atoms.Atoms` object as `structure` input as well as the number of points -`num_of_point` for each compression direction. Depending on the symmetry of the input `structure` the number of -calculations required to calculate the elastic matrix changes. The compression and elongation range is defined by the -`eps_range` parameter. Furthermore, `sqrt_eta` and `fit_order` describe how the change in energy over compression and -elongation is fitted to calculate the resulting pressure. - -## Energy Volume Curve -The `EnergyVolumeCurveWorkflow` can be used to calculate the equilibrium properties: equilibrium volume, equilibrium -energy, equilibrium bulk modulus and the pressure derivative of the equilibrium bulk modulus. -``` -from ase.build import bulk -from atomistics.calculators import evaluate_with_lammps, get_potential_by_name -from atomistics.workflows import EnergyVolumeCurveWorkflow - -potential_dataframe = get_potential_by_name( - potential_name='1999--Mishin-Y--Al--LAMMPS--ipr1' -) -workflow = EnergyVolumeCurveWorkflow( - structure=bulk("Al", cubic=True), - num_points=11, - fit_type="polynomial", - fit_order=3, - vol_range=0.05, - axes=("x", "y", "z"), - strains=None, -) -task_dict = workflow.generate_structures() -result_dict = evaluate_with_lammps( - task_dict=task_dict, - potential_dataframe=potential_dataframe, -) -fit_dict = workflow.analyse_structures(output_dict=result_dict) -print(fit_dict) -``` -The input parameters for the `EnergyVolumeCurveWorkflow` in addition to the `ase.atoms.Atoms` object defined -as `structure` are: - -* `num_points` the number of strains to calculate energies and volumes. -* `fit_type` the type of the fit which should be used to calculate the equilibrium properties. This can either be a - `polynomial` fit or a specific equation of state like the Birch equation (`birch`), the Birch-Murnaghan equation - (`birchmurnaghan`) the Murnaghan equation (`murnaghan`), the Pourier Tarnatola eqaution (`pouriertarantola`) or the - Vinet equation (`vinet`). -* `fit_order` for the `polynomial` fit type the order of the polynomial can be set, for the other fit types this - parameter is ignored. -* `vol_range` specifies the amount of compression and elongation to be applied relative to the absolute volume. -* `axes` specifies the axes which are compressed, typically a uniform compression is applied. -* `strains` specifies the strains directly rather than deriving them from the range of volume compression `vol_range`. - -Beyond calculating the equilibrium properties the `EnergyVolumeCurveWorkflow` can also be used to calculate the thermal -properties using the [Moruzzi, V. L. et al.](https://link.aps.org/doi/10.1103/PhysRevB.37.790) model: -``` -tp_dict = workflow.get_thermal_properties( - t_min=1, - t_max=1500, - t_step=50, - temperatures=None, - constant_volume=False, -) -print(tp_dict) -``` -Or alternatively directly calculate the thermal expansion: -``` -thermal_properties_dict = workflow.get_thermal_properties( - t_min=1, - t_max=1500, - t_step=50, - temperatures=None, - constant_volume=False, - output_keys=["temperatures", "volumes"], -) -temperatures, volumes = thermal_properties_dict["temperatures"], thermal_properties_dict["volumes"] -``` -The [Moruzzi, V. L. et al.](https://link.aps.org/doi/10.1103/PhysRevB.37.790) model is a quantum mechanical approximation, so the equilibrium volume at 0K is not -the same as the equilibrium volume calculated by fitting the equation of state. - -## Molecular Dynamics -Just like the structure optimization also the molecular dynamics calculation can either be implemented inside the -simulation code or in the `atomistics` package. The latter has the advantage that it is the same implementation for all -different simulation codes, while the prior has the advantage that it is usually faster and computationally more -efficient. - -### Implemented in the Simulation Code -The [LAMMPS](https://lammps.org/) simulation code implements a wide range of different simulation workflows, this -includes molecular dynamics. In the `atomistics` package these can be directly accessed via the python interface. - -#### Langevin Thermostat -The Langevin thermostat is currently the only thermostat which is available as both a stand-alone python interface and -an integrated interface inside the [LAMMPS](https://lammps.org/) simulation code. The latter is introduced here: -``` -from ase.build import bulk -from atomistics.calculators import ( - calc_molecular_dynamics_langevin_with_lammps, - get_potential_by_name, -) - -potential_dataframe = get_potential_by_name( - potential_name='1999--Mishin-Y--Al--LAMMPS--ipr1' -) -result_dict = calc_molecular_dynamics_langevin_with_lammps( - structure=bulk("Al", cubic=True).repeat([10, 10, 10]), - potential_dataframe=potential_dataframe, - Tstart=100, - Tstop=100, - Tdamp=0.1, - run=100, - thermo=10, - timestep=0.001, - seed=4928459, - dist="gaussian", - output_keys=("positions", "cell", "forces", "temperature", "energy_pot", "energy_tot", "pressure", "velocities"), -) -``` -In addition to the typical LAMMPS input parameters like the atomistic structure `structure` as `ase.atoms.Atoms` object -and the `pandas.DataFrame` for the interatomic potential `potential_dataframe` are: - -* `Tstart` start temperature -* `Tstop` end temperature -* `Tdamp` temperature damping parameter -* `run` number of molecular dynamics steps to be executed during one temperature step -* `thermo` refresh rate for the thermo dynamic properties, this should typically be the same as the number of molecular - dynamics steps. -* `timestep` time step - typically 1fs defined as `0.001`. -* `seed` random seed for the molecular dynamics -* `dist` initial velocity distribution -* `lmp` Lammps library instance as `pylammpsmpi.LammpsASELibrary` object -* `output` the output quantities which are extracted from the molecular dynamics simulation - -#### Nose Hoover Thermostat -Canonical ensemble (nvt) - volume and temperature constraints molecular dynamics: -``` -from ase.build import bulk -from atomistics.calculators import ( - calc_molecular_dynamics_nvt_with_lammps, - get_potential_by_name, -) - -potential_dataframe = get_potential_by_name( - potential_name='1999--Mishin-Y--Al--LAMMPS--ipr1' -) -result_dict = calc_molecular_dynamics_nvt_with_lammps( - structure=bulk("Al", cubic=True).repeat([10, 10, 10]), - potential_dataframe=potential_dataframe, - Tstart=100, - Tstop=100, - Tdamp=0.1, - run=100, - thermo=10, - timestep=0.001, - seed=4928459, - dist="gaussian", - output_keys=("positions", "cell", "forces", "temperature", "energy_pot", "energy_tot", "pressure", "velocities"), -) -``` -In addition to the typical LAMMPS input parameters like the atomistic structure `structure` as `ase.atoms.Atoms` object -and the `pandas.DataFrame` for the interatomic potential `potential_dataframe` are: - -* `Tstart` start temperature -* `Tstop` end temperature -* `Tdamp` temperature damping parameter -* `run` number of molecular dynamics steps to be executed during one temperature step -* `thermo` refresh rate for the thermo dynamic properties, this should typically be the same as the number of molecular - dynamics steps. -* `timestep` time step - typically 1fs defined as `0.001`. -* `seed` random seed for the molecular dynamics -* `dist` initial velocity distribution -* `lmp` Lammps library instance as `pylammpsmpi.LammpsASELibrary` object -* `output` the output quantities which are extracted from the molecular dynamics simulation - -Isothermal-isobaric ensemble (npt) - pressure and temperature constraints molecular dynamics: -``` -from ase.build import bulk -from atomistics.calculators import ( - calc_molecular_dynamics_npt_with_lammps, - get_potential_by_name, -) - -potential_dataframe = get_potential_by_name( - potential_name='1999--Mishin-Y--Al--LAMMPS--ipr1' -) -result_dict = calc_molecular_dynamics_npt_with_lammps( - structure=bulk("Al", cubic=True).repeat([10, 10, 10]), - potential_dataframe=potential_dataframe, - Tstart=100, - Tstop=100, - Tdamp=0.1, - run=100, - thermo=100, - timestep=0.001, - Pstart=0.0, - Pstop=0.0, - Pdamp=1.0, - seed=4928459, - dist="gaussian", - output_keys=("positions", "cell", "forces", "temperature", "energy_pot", "energy_tot", "pressure", "velocities"), -) -``` -The input parameters for the isothermal-isobaric ensemble (npt) are the same as for the canonical ensemble (nvt) plus: - -* `Pstart` start pressure -* `Pstop` end pressure -* `Pdamp` pressure damping parameter - -Isenthalpic ensemble (nph) - pressure and helmholtz-energy constraints molecular dynamics: -``` -from ase.build import bulk -from atomistics.calculators import ( - calc_molecular_dynamics_nph_with_lammps, - get_potential_by_name, -) - -potential_dataframe = get_potential_by_name( - potential_name='1999--Mishin-Y--Al--LAMMPS--ipr1' -) -result_dict = calc_molecular_dynamics_nph_with_lammps( - structure=bulk("Al", cubic=True).repeat([10, 10, 10]), - potential_dataframe=potential_dataframe, - run=100, - thermo=100, - timestep=0.001, - Tstart=100, - Pstart=0.0, - Pstop=0.0, - Pdamp=1.0, - seed=4928459, - dist="gaussian", - output_keys=("positions", "cell", "forces", "temperature", "energy_pot", "energy_tot", "pressure", "velocities"), -) -``` - -#### Thermal Expansion -One example of a molecular dynamics calculation with the LAMMPS simulation code is the calculation of the thermal -expansion: -``` -from ase.build import bulk -from atomistics.calculators import ( - calc_molecular_dynamics_thermal_expansion_with_lammps, - get_potential_by_name, -) - -potential_dataframe = get_potential_by_name( - potential_name='1999--Mishin-Y--Al--LAMMPS--ipr1' -) -temperatures_md, volumes_md = calc_molecular_dynamics_thermal_expansion_with_lammps( - structure=bulk("Al", cubic=True).repeat([10, 10, 10]), - potential_dataframe=potential_dataframe, - Tstart=100, - Tstop=1000, - Tstep=100, - Tdamp=0.1, - run=100, - thermo=100, - timestep=0.001, - Pstart=0.0, - Pstop=0.0, - Pdamp=1.0, - seed=4928459, - dist="gaussian", - lmp=None, -) -``` -In addition to the typical LAMMPS input parameters like the atomistic structure `structure` as `ase.atoms.Atoms` object -and the `pandas.DataFrame` for the interatomic potential `potential_dataframe` are: - -* `Tstart` start temperature -* `Tstop` end temperature -* `Tstep` temperature step -* `Tdamp` temperature damping parameter -* `run` number of molecular dynamics steps to be executed during one temperature step -* `thermo` refresh rate for the thermo dynamic properties, this should typically be the same as the number of molecular - dynamics steps. -* `timestep` time step - typically 1fs defined as `0.001`. -* `Pstart` start pressure -* `Pstop` end pressure -* `Pdamp` pressure damping parameter -* `seed` random seed for the molecular dynamics -* `dist` initial velocity distribution -* `lmp` Lammps library instance as `pylammpsmpi.LammpsASELibrary` object - -These input parameters are based on the LAMMPS fix `nvt/npt`, you can read more about the specific implementation on the -[LAMMPS website](https://docs.lammps.org/fix_nh.html). - -#### Phonons from Molecular Dynamics -The softening of the phonon modes is calculated for Silicon using the [Tersoff interatomic potential](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.38.9902) -which is available via the [NIST potentials repository](https://www.ctcms.nist.gov/potentials/entry/1988--Tersoff-J--Si-c/). -Silicon is chosen based on its diamond crystal lattice which requires less calculation than the face centered cubic (fcc) -crystal of Aluminium. The simulation workflow consists of three distinct steps: - -* Starting with the optimization of the equilibrium structure. -* Followed by the calculation of the 0K phonon spectrum. -* Finally, the finite temperature phonon spectrum is calculated using molecular dynamics. - -The finite temperature phonon spectrum is calculated using the [DynaPhoPy](https://abelcarreras.github.io/DynaPhoPy/) -package, which is integrated inside the `atomistics` package. As a prerequisite the dependencies, imported and the bulk -silicon diamond structure is created and the Tersoff interatomic potential is loaded: -``` -from ase.build import bulk -from atomistics.calculators import ( - calc_molecular_dynamics_phonons_with_lammps, - evaluate_with_lammps, -) -from atomistics.workflows import optimize_positions_and_volume, PhonopyWorkflow -from dynaphopy import Quasiparticle -import pandas -from phonopy.units import VaspToTHz -import spglib - -structure_bulk = bulk("Si", cubic=True) -potential_dataframe = get_potential_by_name( - potential_name='1988--Tersoff-J--Si-c--LAMMPS--ipr1' -) -``` - -The first step is optimizing the Silicon diamond structure to match the lattice specifications implemented in the Tersoff -interatomic potential: -``` -task_dict = optimize_positions_and_volume(structure=structure_bulk) -result_dict = evaluate_with_lammps( - task_dict=task_dict, - potential_dataframe=potential_dataframe, -) -structure_ase = result_dict["structure_with_optimized_positions_and_volume"] -``` - -As a second step the 0K phonons are calculated using the `PhonopyWorkflow` which is explained in more detail below in -the section on [Phonons](https://atomistics.readthedocs.io/en/latest/workflows.html#phonons). -``` -cell = (structure_ase.cell.array, structure_ase.get_scaled_positions(), structure_ase.numbers) -primitive_matrix = spglib.standardize_cell(cell=cell, to_primitive=True)[0] / structure_ase.get_volume() ** (1/3) -workflow = PhonopyWorkflow( - structure=structure_ase, - interaction_range=10, - factor=VaspToTHz, - displacement=0.01, - dos_mesh=20, - primitive_matrix=primitive_matrix, - number_of_snapshots=None, -) -task_dict = workflow.generate_structures() -result_dict = evaluate_with_lammps( - task_dict=task_dict, - potential_dataframe=potential_dataframe, -) -workflow.analyse_structures(output_dict=result_dict) -``` - -The calcualtion of the finite temperature phonons starts by computing the molecular dynamics trajectory using the -`calc_molecular_dynamics_phonons_with_lammps()` function. This function is internally linked to [DynaPhoPy](https://abelcarreras.github.io/DynaPhoPy/) -to return an `dynaphopy.dynamics.Dynamics` object: -``` -trajectory = calc_molecular_dynamics_phonons_with_lammps( - structure_ase=structure_ase, - potential_dataframe=potential_dataframe, - force_constants=workflow.phonopy.get_force_constants(), - phonopy_unitcell=workflow.phonopy.get_unitcell(), - phonopy_primitive_matrix=workflow.phonopy.get_primitive_matrix(), - phonopy_supercell_matrix=workflow.phonopy.get_supercell_matrix(), - total_time=2, # ps - time_step=0.001, # ps - relaxation_time=5, # ps - silent=True, - supercell=[2, 2, 2], - memmap=False, - velocity_only=True, - temperature=100, -) -``` -When a total of 2 picoseconds is selected to compute the finite temperature phonons with a timestep of 1 femto second -then this results in a total of 2000 molecular dynamics steps. While more molecular dynamics steps result in more precise -predictions they also require more computational resources. - -The postprocessing is executed using the [DynaPhoPy](https://abelcarreras.github.io/DynaPhoPy/) package: -``` -calculation = Quasiparticle(trajectory) -calculation.select_power_spectra_algorithm(2) # select FFT algorithm -calculation.get_renormalized_phonon_dispersion_bands() -renormalized_force_constants = calculation.get_renormalized_force_constants().get_array() -renormalized_force_constants -``` -It calculates the re-normalized force constants which can then be used to calculate the finite temperature properties. - -In addition the [DynaPhoPy](https://abelcarreras.github.io/DynaPhoPy/) package can be used to directly compare the -finite temperature phonon spectrum with the 0K phonon spectrum calulated with the finite displacement method: -``` -calculation.plot_renormalized_phonon_dispersion_bands() -``` -![finite temperature band_structure](pictures/lammps_md_phonons.png) - -### Langevin Thermostat -In addition to the molecular dynamics implemented in the LAMMPS simulation code, the `atomistics` package also provides -the `LangevinWorkflow` which implements molecular dynamics independent of the specific simulation code. -``` -from ase.build import bulk -from atomistics.calculators import evaluate_with_lammps_library, get_potential_by_name -from atomistics.workflows import LangevinWorkflow -from pylammpsmpi import LammpsASELibrary - -steps = 300 -potential_dataframe = get_potential_by_name( - potential_name='1999--Mishin-Y--Al--LAMMPS--ipr1' -) -workflow = LangevinWorkflow( - structure=bulk("Al", cubic=True).repeat([2, 2, 2]), - temperature=1000.0, - overheat_fraction=2.0, - damping_timescale=100.0, - time_step=1, -) -lmp = LammpsASELibrary( - working_directory=None, - cores=1, - comm=None, - logger=None, - log_file=None, - library=None, - diable_log_file=True, -) -eng_pot_lst, eng_kin_lst = [], [] -for i in range(steps): - task_dict = workflow.generate_structures() - result_dict = evaluate_with_lammps_library( - task_dict=task_dict, - potential_dataframe=potential_dataframe, - lmp=lmp, - ) - eng_pot, eng_kin = workflow.analyse_structures(output_dict=result_dict) - eng_pot_lst.append(eng_pot) - eng_kin_lst.append(eng_kin) -lmp.close() -``` -The advantage of this implementation is that the user can directly interact with the simulation between the individual -molecular dynamics simulation steps. This provides a lot of flexibility to prototype new simulation methods. The input -parameters of the `LangevinWorkflow` are: - -* `structure` the `ase.atoms.Atoms` object which is used as initial structure for the molecular dynamics calculation -* `temperature` the temperature of the molecular dynamics calculation given in Kelvin -* `overheat_fraction` the over heating fraction of the Langevin thermostat -* `damping_timescale` the damping timescale of the Langevin thermostat -* `time_step` the time steps of the Langevin thermostat - -## Harmonic Approximation -The harmonic approximation is implemented in two variations, once with constant volume and once including the volume -expansion at finite temperature also known as quasi-harmonic approximation. Both of these are based on the [phonopy](https://phonopy.github.io/phonopy/) -package. - -### Phonons -To calculate the phonons at a fixed volume the `PhonopyWorkflow` is used: -``` -from ase.build import bulk -from atomistics.calculators import evaluate_with_lammps, get_potential_by_name -from atomistics.workflows import PhonopyWorkflow -from phonopy.units import VaspToTHz - -potential_dataframe = get_potential_by_name( - potential_name='1999--Mishin-Y--Al--LAMMPS--ipr1' -) -workflow = PhonopyWorkflow( - structure=bulk("Al", cubic=True), - interaction_range=10, - factor=VaspToTHz, - displacement=0.01, - dos_mesh=20, - primitive_matrix=None, - number_of_snapshots=None, -) -task_dict = workflow.generate_structures() -result_dict = evaluate_with_lammps( - task_dict=task_dict, - potential_dataframe=potential_dataframe, -) -phonopy_dict = workflow.analyse_structures(output_dict=result_dict) -``` -The `PhonopyWorkflow` takes the following inputs: - -* `structure` the `ase.atoms.Atoms` object to calculate the phonon spectrum -* `interaction_range` the cutoff radius to consider for identifying the interaction between the atoms -* `factor` conversion factor, typically just `phonopy.units.VaspToTHz` -* `displacement` displacement to calculate the forces -* `dos_mesh` mesh for the density of states -* `primitive_matrix` primitive matrix -* `number_of_snapshots` number of snapshots to calculate - -In addition to the phonon properties, the `PhonopyWorkflow` also enables the calculation of thermal properties: -``` -tp_dict = workflow.get_thermal_properties( - t_min=1, - t_max=1500, - t_step=50, - temperatures=None, - cutoff_frequency=None, - pretend_real=False, - band_indices=None, - is_projection=False, -) -print(tp_dict) -``` -The calculation of the thermal properties takes additional inputs: - -* `t_min` minimum temperature -* `t_max` maximum temperature -* `t_step` temperature step -* `temperatures` alternative to `t_min`, `t_max` and `t_step` the array of temperatures can be defined directly -* `cutoff_frequency` cutoff frequency to exclude the contributions of frequencies below a certain cut off -* `pretend_real` use the absolute values of the phonon frequencies -* `band_indices` select bands based on their indices -* `is_projection` multiplies the squared eigenvectors - not recommended - -Furthermore, also the dynamical matrix can be directly calculated with the `PhonopyWorkflow`: -``` -mat = workflow.get_dynamical_matrix() -``` - -Or alternatively the hesse matrix: -``` -mat = workflow.get_hesse_matrix() -``` - -Finally, also the function to calculate the band structure is directly available on the `PhonopyWorkflow`: -``` -band_structure = workflow.get_band_structure( - npoints=101, - with_eigenvectors=False, - with_group_velocities=False -) -``` - -This band structure can also be visualised using the built-in plotting function: -``` -workflow.plot_band_structure() -``` -![band_structure](pictures/phonon_bands_al.png) - -Just like the desnsity of states which can be plotted using: -``` -workflow.plot_dos() -``` -![density of states](pictures/phonon_dos_al.png) - -### Quasi-harmonic Approximation -To include the volume expansion with finite temperature the `atomistics` package implements the `QuasiHarmonicWorkflow`: -``` -from ase.build import bulk -from atomistics.calculators import evaluate_with_lammps, get_potential_by_name -from atomistics.workflows import QuasiHarmonicWorkflow - -potential_dataframe = get_potential_by_name( - potential_name='1999--Mishin-Y--Al--LAMMPS--ipr1' -) -workflow = QuasiHarmonicWorkflow( - structure=bulk("Al", cubic=True), - num_points=11, - vol_range=0.05, - interaction_range=10, - factor=VaspToTHz, - displacement=0.01, - dos_mesh=20, - primitive_matrix=None, - number_of_snapshots=None, -) -task_dict = workflow.generate_structures() -result_dict = evaluate_with_lammps( - task_dict=task_dict, - potential_dataframe=potential_dataframe, -) -fit_dict = workflow.analyse_structures(output_dict=result_dict) -``` -The `QuasiHarmonicWorkflow` is a combination of the `EnergyVolumeCurveWorkflow` and the `PhonopyWorkflow`. Consequently, -the inputs are a superset of the inputs of these two workflows. - -Based on the `QuasiHarmonicWorkflow` the thermal properties can be calculated: -``` -tp_dict = workflow.get_thermal_properties( - t_min=1, - t_max=1500, - t_step=50, - temperatures=None, - cutoff_frequency=None, - pretend_real=False, - band_indices=None, - is_projection=False, - quantum_mechanical=True, -) -``` -This requires the same inputs as the calculation of the thermal properties `get_thermal_properties()` with the -`PhonopyWorkflow`. The additional parameter `quantum_mechanical` specifies whether the classical harmonic oscillator or -the quantum mechanical harmonic oscillator is used to calculate the free energy. - -And finally also the thermal expansion can be calculated: -``` -tp_dict = workflow.get_thermal_properties( - t_min=1, - t_max=1500, - t_step=50, - temperatures=None, - cutoff_frequency=None, - pretend_real=False, - band_indices=None, - is_projection=False, - quantum_mechanical=True, - output_keys=["temperatures", "volumes"], -) -temperatures, volumes = tp_dict["temperatures"], tp_dict["volumes"] -``` - -## Structure Optimization -In analogy to the molecular dynamics calculation also the structure optimization could in principle be defined inside -the simulation code or on the python level. Still currently the `atomistics` package only supports the structure -optimization defined inside the simulation codes. - -### Volume and Positions -To optimize both the volume of the supercell as well as the positions inside the supercell the `atomistics` package -implements the `optimize_positions_and_volume()` workflow: -``` -from ase.build import bulk -from atomistics.calculators import evaluate_with_lammps, get_potential_by_name -from atomistics.workflows import optimize_positions_and_volume - -structure = bulk("Al", a=4.0, cubic=True) -potential_dataframe = get_potential_by_name( - potential_name='1999--Mishin-Y--Al--LAMMPS--ipr1' -) -result_dict = evaluate_with_lammps( - task_dict=optimize_positions_and_volume(structure=structure), - potential_dataframe=potential_dataframe, -) -structure_opt = result_dict["structure_with_optimized_positions_and_volume"] -``` -The result is the optimized atomistic structure as part of the result dictionary. - -### Positions -The optimization of the positions inside the supercell without the optimization of the supercell volume is possible with -the `optimize_positions()` workflow: -``` -from ase.build import bulk -from atomistics.calculators import evaluate_with_lammps, get_potential_by_name -from atomistics.workflows import optimize_positions - -structure = bulk("Al", a=4.0, cubic=True) -potential_dataframe = get_potential_by_name( - potential_name='1999--Mishin-Y--Al--LAMMPS--ipr1' -) -result_dict = evaluate_with_lammps( - task_dict=optimize_positions(structure=structure), - potential_dataframe=potential_dataframe, -) -structure_opt = result_dict["structure_with_optimized_positions"] -``` -The result is the optimized atomistic structure as part of the result dictionary. \ No newline at end of file diff --git a/notebooks/simulation_codes.ipynb b/notebooks/simulation_codes.ipynb new file mode 100644 index 00000000..93b86337 --- /dev/null +++ b/notebooks/simulation_codes.ipynb @@ -0,0 +1 @@ +{"metadata":{"language_info":{"name":"python","version":"3.10.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kernelspec":{"name":"python3","display_name":"Python 3 (ipykernel)","language":"python"}},"nbformat_minor":5,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Simulation Codes","metadata":{},"id":"dcc94f22-3be1-4cd9-9f1c-e61f7e03a1f5"},{"cell_type":"markdown","source":"## ASE\nAt the current stage the majority of simulation codes are interfaced using the [Atomic Simulation Environment (ASE)](https://wiki.fysik.dtu.dk/ase/).\nThe limitation of the ASE based interfaces is that the simulation codes are only used to calculate energies, forces and\nstresses, while more complex computations like structure optimization or molecular dynamics are implemented in python.","metadata":{},"id":"5ffec124-6252-4965-9267-cc771c79570f"},{"cell_type":"markdown","source":"### Abinit\n[Abinit](https://www.abinit.org) - Plane wave density functional theory:","metadata":{},"id":"ece0e223-3b49-41d7-99ab-e0518ef01c2b"},{"cell_type":"code","source":"from ase.calculators.abinit import Abinit\nfrom ase.units import Ry\nfrom atomistics.calculators import evaluate_with_ase\n\nresult_dict = evaluate_with_ase(\n task_dict={},\n ase_calculator=Abinit(\n label='abinit_evcurve',\n nbands=32,\n ecut=10 * Ry,\n kpts=(3, 3, 3),\n toldfe=1.0e-2,\n v8_legacy_format=False,\n )\n)","metadata":{"trusted":true,"tags":[]},"execution_count":1,"outputs":[],"id":"59904794-1beb-43f5-b6e2-f79404ca8b46"},{"cell_type":"markdown","source":"The full documentation of the corresponding interface is available on the [Atomic Simulation Environment](https://wiki.fysik.dtu.dk/ase/ase/calculators/abinit.html)\nwebsite. ","metadata":{},"id":"53687acd-cc66-470c-b88b-c4a2a7f4182c"},{"cell_type":"markdown","source":"### EMT\n[EMT](https://wiki.fysik.dtu.dk/ase/ase/calculators/emt.html) - Effective medium theory: ","metadata":{},"id":"a765a024-fc04-4add-b6c2-eef026e4dbd2"},{"cell_type":"code","source":"from ase.calculators.emt import EMT\nfrom atomistics.calculators import evaluate_with_ase\n\nresult_dict = evaluate_with_ase(\n task_dict={}, \n ase_calculator=EMT()\n)","metadata":{"trusted":true,"tags":[]},"execution_count":2,"outputs":[],"id":"a7467756-1366-42c2-b595-e8c7893437ab"},{"cell_type":"markdown","source":"The full documentation of the corresponding interface is available on the [Atomic Simulation Environment](https://wiki.fysik.dtu.dk/ase/ase/calculators/emt.html)\nwebsite. ","metadata":{},"id":"e635174f-be67-462f-a3aa-1d3103c31d97"},{"cell_type":"markdown","source":"### GPAW\n[GPAW](https://wiki.fysik.dtu.dk/gpaw/) - Density functional theory Python code based on the projector-augmented wave \nmethod:","metadata":{},"id":"af758cea-e4ad-475e-8fcc-2afbf497334d"},{"cell_type":"code","source":"from gpaw import GPAW, PW\nfrom atomistics.calculators import evaluate_with_ase\n\nresult_dict = evaluate_with_ase(\n task_dict={}, \n ase_calculator=GPAW(\n xc=\"PBE\",\n mode=PW(300),\n kpts=(3, 3, 3)\n )\n)","metadata":{"trusted":true,"tags":[]},"execution_count":3,"outputs":[{"name":"stdout","text":"\n ___ ___ ___ _ _ _ \n | | |_ | | | | \n | | | | | . | | | | \n |__ | _|___|_____| 24.1.0\n |___|_| \n\nUser: jovyan@jupyter-pyiron-2datomistics-2dha1r3gln\nDate: Wed May 1 16:53:34 2024\nArch: x86_64\nPid: 648\nCWD: /home/jovyan/notebooks\nPython: 3.10.12\ngpaw: /srv/conda/envs/notebook/lib/python3.10/site-packages/gpaw\n_gpaw: /srv/conda/envs/notebook/lib/python3.10/site-packages/\n _gpaw.cpython-310-x86_64-linux-gnu.so\nase: /srv/conda/envs/notebook/lib/python3.10/site-packages/ase (version 3.22.1)\nnumpy: /srv/conda/envs/notebook/lib/python3.10/site-packages/numpy (version 1.26.4)\nscipy: /srv/conda/envs/notebook/lib/python3.10/site-packages/scipy (version 1.13.0)\nlibxc: 6.2.2\nunits: Angstrom and eV\ncores: 1\nOpenMP: True\nOMP_NUM_THREADS: 1\n\nInput parameters:\n kpts: [3 3 3]\n mode: {ecut: 300.0,\n name: pw}\n xc: PBE\n\nMemory usage: 143.04 MiB\nDate: Wed May 1 16:53:34 2024\n","output_type":"stream"}],"id":"10a942f6-b4c0-4710-856f-e736a643ce17"},{"cell_type":"markdown","source":"The full documentation of the corresponding interface is available on the [GPAW](https://wiki.fysik.dtu.dk/gpaw/)\nwebsite.","metadata":{},"id":"a38a03dd-630c-4bff-a256-d0380da29db1"},{"cell_type":"markdown","source":"### Quantum Espresso \n[Quantum Espresso](https://www.quantum-espresso.org) - Integrated suite of Open-Source computer codes for \nelectronic-structure calculations:","metadata":{},"id":"f8a20ee4-153d-4c0b-b89f-a3ac64cfbcbc"},{"cell_type":"code","source":"from ase.calculators.espresso import Espresso\nfrom atomistics.calculators import evaluate_with_ase\n\nresult_dict = evaluate_with_ase(\n task_dict={}, \n ase_calculator=Espresso(\n pseudopotentials={\"Al\": \"Al.pbe-n-kjpaw_psl.1.0.0.UPF\"},\n tstress=True,\n tprnfor=True,\n kpts=(3, 3, 3),\n )\n)","metadata":{"trusted":true,"tags":[]},"execution_count":4,"outputs":[],"id":"40b9384e-02a0-41d8-adff-ad4dd6316a21"},{"cell_type":"markdown","source":"The full documentation of the corresponding interface is available on the [Atomic Simulation Environment](https://wiki.fysik.dtu.dk/ase/ase/calculators/espresso.html)\nwebsite. ","metadata":{},"id":"c7f74b95-e161-4ee6-8a9b-58c11ba3ca47"},{"cell_type":"markdown","source":"### Siesta\n[Siesta](https://siesta-project.org) - Electronic structure calculations and ab initio molecular dynamics:","metadata":{},"id":"cdc1c88c-62f8-4af1-a135-9266ba1ab1ee"},{"cell_type":"code","source":"import os\nfrom ase.calculators.siesta import Siesta\nfrom ase.units import Ry\nfrom atomistics.calculators import evaluate_with_ase\n\nresult_dict = evaluate_with_ase(\n task_dict={}, \n ase_calculator=Siesta(\n label=\"siesta\",\n xc=\"PBE\",\n mesh_cutoff=200 * Ry,\n energy_shift=0.01 * Ry,\n basis_set=\"DZ\",\n kpts=(5, 5, 5),\n fdf_arguments={\"DM.MixingWeight\": 0.1, \"MaxSCFIterations\": 100},\n pseudo_path=os.path.abspath(\"tests/static/siesta\"),\n pseudo_qualifier=\"\",\n )\n)","metadata":{"trusted":true,"tags":[]},"execution_count":5,"outputs":[],"id":"5fa73a10-043f-4b46-a162-3f1da7399cba"},{"cell_type":"markdown","source":"The full documentation of the corresponding interface is available on the [Atomic Simulation Environment](https://wiki.fysik.dtu.dk/ase/ase/calculators/siesta.html)\nwebsite.","metadata":{},"id":"ed3bba6d-8fd9-4d26-8f0d-437fb5559397"},{"cell_type":"markdown","source":"## LAMMPS\n[LAMMPS](https://www.lammps.org) - Molecular Dynamics:","metadata":{},"id":"0c09855c-80ca-4fb9-9f66-458b8d145e59"},{"cell_type":"code","source":"from ase.build import bulk\nfrom atomistics.calculators import evaluate_with_lammps, get_potential_by_name\n\nstructure = bulk(\"Al\", cubic=True)\npotential_dataframe = get_potential_by_name(\n potential_name='1999--Mishin-Y--Al--LAMMPS--ipr1',\n resource_path=\"../tests/static/lammps\"\n)\n\nresult_dict = evaluate_with_lammps(\n task_dict={},\n potential_dataframe=potential_dataframe,\n)","metadata":{"trusted":true,"tags":[]},"execution_count":6,"outputs":[{"name":"stderr","text":"/srv/conda/envs/notebook/lib/python3.10/site-packages/atomistics/calculators/lammps/potential.py:299: SettingWithCopyWarning: \nA value is trying to be set on a copy of a slice from a DataFrame.\nTry using .loc[row_indexer,col_indexer] = value instead\n\nSee the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n df_pot[\"Config\"] = config_lst\n","output_type":"stream"}],"id":"a95e1c9c-1ad1-4765-aa13-790db6971dab"},{"cell_type":"markdown","source":"The [LAMMPS](https://www.lammps.org) interface is based on the [pylammpsmpi](https://github.com/pyiron/pylammpsmpi)\npackage which couples a [LAMMPS](https://www.lammps.org) instance which is parallelized via the Message Passing Interface\n(MPI) with a serial python process or jupyter notebook. The challenging part about molecular dynamics simulation is \nidentifying a suitable interatomic potential. \n\nTo address this challenge the `atomistics` package is leveraging the [NIST database of interatomic potentials](https://www.ctcms.nist.gov/potentials). \nIt is recommended to install this database `iprpy-data` via the `conda` package manager, then the `resource_path` is\nautomatically set to `${CONDA_PREFIX}/share/iprpy`. Alternatively, the `resource_path` can be specified manually as an\noptional parameter of the `get_potential_by_name()` function.\n\nIn addition, the `get_potential_dataframe(structure)` function which takes an `ase.atoms.Atoms` object as input can be\nused to query the [NIST database of interatomic potentials](https://www.ctcms.nist.gov/potentials) for potentials, which\ninclude the interatomic interactions required to simulate the atomic structure defined by the `ase.atoms.Atoms` object. \nIt returns a `pandas.DataFrame` with all the available potentials and the `resource_path` can again be specified as \noptional parameter.\n\nFinally, another option to specify the interatomic potential for a LAMMPS simulation is by defining the `potential_dataframe`\ndirectly: ","metadata":{},"id":"ed5732c4-4221-4081-83d8-2c1df793d8b7"},{"cell_type":"code","source":"import pandas \n\npotential_dataframe = pandas.DataFrame({\n \"Config\": [[\n \"pair_style morse/smooth/linear 9.0\",\n \"pair_coeff * * 0.5 1.8 2.95\"\n ]],\n \"Filename\": [[]],\n \"Model\": [\"Morse\"],\n \"Name\": [\"Morse\"],\n \"Species\": [[\"Al\"]],\n})","metadata":{"trusted":true,"tags":[]},"execution_count":7,"outputs":[],"id":"4ab54444-0d8b-42dd-8a97-c0743598a7bd"},{"cell_type":"markdown","source":"## Quantum Espresso\n[Quantum Espresso](https://www.quantum-espresso.org) - Integrated suite of Open-Source computer codes for \nelectronic-structure calculations:","metadata":{},"id":"0fbec738-da33-44bf-ab6f-b222f1d56186"},{"cell_type":"markdown","source":"from atomistics.calculators import evaluate_with_qe\n\nresult_dict = evaluate_with_qe(\n task_dict={},\n calculation_name=\"espresso\",\n working_directory=\".\",\n kpts=(3, 3, 3),\n pseudopotentials={\n \"Al\": \"Al.pbe-n-kjpaw_psl.1.0.0.UPF\"\n },\n tstress=True,\n tprnfor=True,\n ecutwfc=40.0, # kinetic energy cutoff (Ry) for wavefunctions\n conv_thr=1e-06, # Convergence threshold for selfconsistency\n diagonalization='david', \n electron_maxstep=100, # maximum number of iterations in a scf step. \n nstep=200, # number of molecular-dynamics or structural optimization steps performed in this run.\n etot_conv_thr=1e-4, # Convergence threshold on total energy (a.u) for ionic minimization\n forc_conv_thr=1e-3, # Convergence threshold on forces (a.u) for ionic minimization\n smearing='gaussian', # ordinary Gaussian spreading (Default)\n)","metadata":{"tags":[]},"id":"e7a8b26d-5547-4e69-b224-d65f38060c0a"},{"cell_type":"markdown","source":"This secondary interface for [Quantum Espresso](https://www.quantum-espresso.org) is based on the input writer from the\n[Atomic Simulation Environment (ASE)](https://wiki.fysik.dtu.dk/ase/) and the output is parsed using the [pwtools](https://elcorto.github.io/pwtools/).\nThe executable can be set using the `ASE_ESPRESSO_COMMAND` environment variable:","metadata":{},"id":"727f824e-97a2-42ae-94f2-b8080c023e24"},{"cell_type":"raw","source":"export ASE_ESPRESSO_COMMAND=\"pw.x -in PREFIX.pwi > PREFIX.pwo\"","metadata":{},"id":"a0207a31-fa83-4729-9b4d-6a025af78980"},{"cell_type":"markdown","source":"The full list of possible keyword arguments is available in the [Quantum Espresso Documentation](https://www.quantum-espresso.org/Doc/INPUT_PW.html).\nFinally, the [Standard solid-state pseudopotentials (SSSP)](https://www.materialscloud.org/discover/sssp/table/efficiency) \nfor quantum espresso are distributed via the materials cloud.","metadata":{"tags":[]},"id":"54ee1c7b-a829-49b9-8746-d57496b339e6"}]} \ No newline at end of file From 5ab613d3fe1947882db720792dddd49ce0314388 Mon Sep 17 00:00:00 2001 From: Jan Janssen Date: Wed, 1 May 2024 11:11:36 -0600 Subject: [PATCH 4/9] copy notebooks --- .readthedocs.yml | 1 + docs/_toc.yml | 10 +++++----- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/.readthedocs.yml b/.readthedocs.yml index ffe36490..e08dac99 100644 --- a/.readthedocs.yml +++ b/.readthedocs.yml @@ -13,6 +13,7 @@ build: pre_build: # Generate the Sphinx configuration for this Jupyter Book so it builds. - "cp README.md docs" + - "cp notebooks/*.ipynb docs" - "jupyter-book config sphinx docs/" # Build documentation in the docs/ directory with Sphinx diff --git a/docs/_toc.yml b/docs/_toc.yml index 66724c71..208a1dec 100644 --- a/docs/_toc.yml +++ b/docs/_toc.yml @@ -2,11 +2,11 @@ format: jb-book root: README chapters: - file: installation.md -- file: ../notebooks/simulation_codes.ipynb -- file: ../notebooks/lammps_workflows.ipynb +- file: simulation_codes.ipynb +- file: lammps_workflows.ipynb - file: materialproperties.md sections: - - file: ../notebooks/bulk_modulus_with_gpaw.ipynb - - file: ../thermal_expansion_with_lammps.ipynb - - file: ../free_energy_calculation.ipynb + - file: bulk_modulus_with_gpaw.ipynb + - file: thermal_expansion_with_lammps.ipynb + - file: free_energy_calculation.ipynb - file: phasediagram.md \ No newline at end of file From fd0b031343580a440cc97319c0a5bdd9dc9ef866 Mon Sep 17 00:00:00 2001 From: Jan Janssen Date: Wed, 1 May 2024 11:46:21 -0600 Subject: [PATCH 5/9] update mybinder environment --- README.md | 2 +- binder/environment.yml | 5 +- notebooks/simulation_codes.ipynb | 460 ++++++++++++++++++++++++++++++- 3 files changed, 463 insertions(+), 4 deletions(-) diff --git a/README.md b/README.md index 2b824ac2..e87e638c 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -# atomistics - Interfaces for atomistic simulation codes and workflows +# atomistics [![Unittest](https://github.com/pyiron/atomistics/actions/workflows/unittests.yml/badge.svg)](https://github.com/pyiron/atomistics/actions/workflows/unittests.yml) [![Coverage Status](https://coveralls.io/repos/github/pyiron/atomistics/badge.svg?branch=main)](https://coveralls.io/github/pyiron/atomistics?branch=main) [![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/pyiron/atomistics/HEAD) diff --git a/binder/environment.yml b/binder/environment.yml index 7716f4fc..0b2b0088 100644 --- a/binder/environment.yml +++ b/binder/environment.yml @@ -2,7 +2,6 @@ channels: - conda-forge dependencies: - ase =3.22.1 -- coverage - numpy =1.26.4 - scipy =1.13.0 - spglib =2.4.0 @@ -10,7 +9,9 @@ dependencies: - structuretoolkit =0.0.22 - seekpath =2.1.0 - gpaw =24.1.0 -- lammps =2023.11.21 +- lammps =2024.02.07 - pandas =2.2.2 - pylammpsmpi =0.2.18 - jinja2 =3.1.3 +- qe =7.2 +- pwtools =1.2.3 diff --git a/notebooks/simulation_codes.ipynb b/notebooks/simulation_codes.ipynb index 93b86337..064a7c43 100644 --- a/notebooks/simulation_codes.ipynb +++ b/notebooks/simulation_codes.ipynb @@ -1 +1,459 @@ -{"metadata":{"language_info":{"name":"python","version":"3.10.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kernelspec":{"name":"python3","display_name":"Python 3 (ipykernel)","language":"python"}},"nbformat_minor":5,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Simulation Codes","metadata":{},"id":"dcc94f22-3be1-4cd9-9f1c-e61f7e03a1f5"},{"cell_type":"markdown","source":"## ASE\nAt the current stage the majority of simulation codes are interfaced using the [Atomic Simulation Environment (ASE)](https://wiki.fysik.dtu.dk/ase/).\nThe limitation of the ASE based interfaces is that the simulation codes are only used to calculate energies, forces and\nstresses, while more complex computations like structure optimization or molecular dynamics are implemented in python.","metadata":{},"id":"5ffec124-6252-4965-9267-cc771c79570f"},{"cell_type":"markdown","source":"### Abinit\n[Abinit](https://www.abinit.org) - Plane wave density functional theory:","metadata":{},"id":"ece0e223-3b49-41d7-99ab-e0518ef01c2b"},{"cell_type":"code","source":"from ase.calculators.abinit import Abinit\nfrom ase.units import Ry\nfrom atomistics.calculators import evaluate_with_ase\n\nresult_dict = evaluate_with_ase(\n task_dict={},\n ase_calculator=Abinit(\n label='abinit_evcurve',\n nbands=32,\n ecut=10 * Ry,\n kpts=(3, 3, 3),\n toldfe=1.0e-2,\n v8_legacy_format=False,\n )\n)","metadata":{"trusted":true,"tags":[]},"execution_count":1,"outputs":[],"id":"59904794-1beb-43f5-b6e2-f79404ca8b46"},{"cell_type":"markdown","source":"The full documentation of the corresponding interface is available on the [Atomic Simulation Environment](https://wiki.fysik.dtu.dk/ase/ase/calculators/abinit.html)\nwebsite. ","metadata":{},"id":"53687acd-cc66-470c-b88b-c4a2a7f4182c"},{"cell_type":"markdown","source":"### EMT\n[EMT](https://wiki.fysik.dtu.dk/ase/ase/calculators/emt.html) - Effective medium theory: ","metadata":{},"id":"a765a024-fc04-4add-b6c2-eef026e4dbd2"},{"cell_type":"code","source":"from ase.calculators.emt import EMT\nfrom atomistics.calculators import evaluate_with_ase\n\nresult_dict = evaluate_with_ase(\n task_dict={}, \n ase_calculator=EMT()\n)","metadata":{"trusted":true,"tags":[]},"execution_count":2,"outputs":[],"id":"a7467756-1366-42c2-b595-e8c7893437ab"},{"cell_type":"markdown","source":"The full documentation of the corresponding interface is available on the [Atomic Simulation Environment](https://wiki.fysik.dtu.dk/ase/ase/calculators/emt.html)\nwebsite. ","metadata":{},"id":"e635174f-be67-462f-a3aa-1d3103c31d97"},{"cell_type":"markdown","source":"### GPAW\n[GPAW](https://wiki.fysik.dtu.dk/gpaw/) - Density functional theory Python code based on the projector-augmented wave \nmethod:","metadata":{},"id":"af758cea-e4ad-475e-8fcc-2afbf497334d"},{"cell_type":"code","source":"from gpaw import GPAW, PW\nfrom atomistics.calculators import evaluate_with_ase\n\nresult_dict = evaluate_with_ase(\n task_dict={}, \n ase_calculator=GPAW(\n xc=\"PBE\",\n mode=PW(300),\n kpts=(3, 3, 3)\n )\n)","metadata":{"trusted":true,"tags":[]},"execution_count":3,"outputs":[{"name":"stdout","text":"\n ___ ___ ___ _ _ _ \n | | |_ | | | | \n | | | | | . | | | | \n |__ | _|___|_____| 24.1.0\n |___|_| \n\nUser: jovyan@jupyter-pyiron-2datomistics-2dha1r3gln\nDate: Wed May 1 16:53:34 2024\nArch: x86_64\nPid: 648\nCWD: /home/jovyan/notebooks\nPython: 3.10.12\ngpaw: /srv/conda/envs/notebook/lib/python3.10/site-packages/gpaw\n_gpaw: /srv/conda/envs/notebook/lib/python3.10/site-packages/\n _gpaw.cpython-310-x86_64-linux-gnu.so\nase: /srv/conda/envs/notebook/lib/python3.10/site-packages/ase (version 3.22.1)\nnumpy: /srv/conda/envs/notebook/lib/python3.10/site-packages/numpy (version 1.26.4)\nscipy: /srv/conda/envs/notebook/lib/python3.10/site-packages/scipy (version 1.13.0)\nlibxc: 6.2.2\nunits: Angstrom and eV\ncores: 1\nOpenMP: True\nOMP_NUM_THREADS: 1\n\nInput parameters:\n kpts: [3 3 3]\n mode: {ecut: 300.0,\n name: pw}\n xc: PBE\n\nMemory usage: 143.04 MiB\nDate: Wed May 1 16:53:34 2024\n","output_type":"stream"}],"id":"10a942f6-b4c0-4710-856f-e736a643ce17"},{"cell_type":"markdown","source":"The full documentation of the corresponding interface is available on the [GPAW](https://wiki.fysik.dtu.dk/gpaw/)\nwebsite.","metadata":{},"id":"a38a03dd-630c-4bff-a256-d0380da29db1"},{"cell_type":"markdown","source":"### Quantum Espresso \n[Quantum Espresso](https://www.quantum-espresso.org) - Integrated suite of Open-Source computer codes for \nelectronic-structure calculations:","metadata":{},"id":"f8a20ee4-153d-4c0b-b89f-a3ac64cfbcbc"},{"cell_type":"code","source":"from ase.calculators.espresso import Espresso\nfrom atomistics.calculators import evaluate_with_ase\n\nresult_dict = evaluate_with_ase(\n task_dict={}, \n ase_calculator=Espresso(\n pseudopotentials={\"Al\": \"Al.pbe-n-kjpaw_psl.1.0.0.UPF\"},\n tstress=True,\n tprnfor=True,\n kpts=(3, 3, 3),\n )\n)","metadata":{"trusted":true,"tags":[]},"execution_count":4,"outputs":[],"id":"40b9384e-02a0-41d8-adff-ad4dd6316a21"},{"cell_type":"markdown","source":"The full documentation of the corresponding interface is available on the [Atomic Simulation Environment](https://wiki.fysik.dtu.dk/ase/ase/calculators/espresso.html)\nwebsite. ","metadata":{},"id":"c7f74b95-e161-4ee6-8a9b-58c11ba3ca47"},{"cell_type":"markdown","source":"### Siesta\n[Siesta](https://siesta-project.org) - Electronic structure calculations and ab initio molecular dynamics:","metadata":{},"id":"cdc1c88c-62f8-4af1-a135-9266ba1ab1ee"},{"cell_type":"code","source":"import os\nfrom ase.calculators.siesta import Siesta\nfrom ase.units import Ry\nfrom atomistics.calculators import evaluate_with_ase\n\nresult_dict = evaluate_with_ase(\n task_dict={}, \n ase_calculator=Siesta(\n label=\"siesta\",\n xc=\"PBE\",\n mesh_cutoff=200 * Ry,\n energy_shift=0.01 * Ry,\n basis_set=\"DZ\",\n kpts=(5, 5, 5),\n fdf_arguments={\"DM.MixingWeight\": 0.1, \"MaxSCFIterations\": 100},\n pseudo_path=os.path.abspath(\"tests/static/siesta\"),\n pseudo_qualifier=\"\",\n )\n)","metadata":{"trusted":true,"tags":[]},"execution_count":5,"outputs":[],"id":"5fa73a10-043f-4b46-a162-3f1da7399cba"},{"cell_type":"markdown","source":"The full documentation of the corresponding interface is available on the [Atomic Simulation Environment](https://wiki.fysik.dtu.dk/ase/ase/calculators/siesta.html)\nwebsite.","metadata":{},"id":"ed3bba6d-8fd9-4d26-8f0d-437fb5559397"},{"cell_type":"markdown","source":"## LAMMPS\n[LAMMPS](https://www.lammps.org) - Molecular Dynamics:","metadata":{},"id":"0c09855c-80ca-4fb9-9f66-458b8d145e59"},{"cell_type":"code","source":"from ase.build import bulk\nfrom atomistics.calculators import evaluate_with_lammps, get_potential_by_name\n\nstructure = bulk(\"Al\", cubic=True)\npotential_dataframe = get_potential_by_name(\n potential_name='1999--Mishin-Y--Al--LAMMPS--ipr1',\n resource_path=\"../tests/static/lammps\"\n)\n\nresult_dict = evaluate_with_lammps(\n task_dict={},\n potential_dataframe=potential_dataframe,\n)","metadata":{"trusted":true,"tags":[]},"execution_count":6,"outputs":[{"name":"stderr","text":"/srv/conda/envs/notebook/lib/python3.10/site-packages/atomistics/calculators/lammps/potential.py:299: SettingWithCopyWarning: \nA value is trying to be set on a copy of a slice from a DataFrame.\nTry using .loc[row_indexer,col_indexer] = value instead\n\nSee the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n df_pot[\"Config\"] = config_lst\n","output_type":"stream"}],"id":"a95e1c9c-1ad1-4765-aa13-790db6971dab"},{"cell_type":"markdown","source":"The [LAMMPS](https://www.lammps.org) interface is based on the [pylammpsmpi](https://github.com/pyiron/pylammpsmpi)\npackage which couples a [LAMMPS](https://www.lammps.org) instance which is parallelized via the Message Passing Interface\n(MPI) with a serial python process or jupyter notebook. The challenging part about molecular dynamics simulation is \nidentifying a suitable interatomic potential. \n\nTo address this challenge the `atomistics` package is leveraging the [NIST database of interatomic potentials](https://www.ctcms.nist.gov/potentials). \nIt is recommended to install this database `iprpy-data` via the `conda` package manager, then the `resource_path` is\nautomatically set to `${CONDA_PREFIX}/share/iprpy`. Alternatively, the `resource_path` can be specified manually as an\noptional parameter of the `get_potential_by_name()` function.\n\nIn addition, the `get_potential_dataframe(structure)` function which takes an `ase.atoms.Atoms` object as input can be\nused to query the [NIST database of interatomic potentials](https://www.ctcms.nist.gov/potentials) for potentials, which\ninclude the interatomic interactions required to simulate the atomic structure defined by the `ase.atoms.Atoms` object. \nIt returns a `pandas.DataFrame` with all the available potentials and the `resource_path` can again be specified as \noptional parameter.\n\nFinally, another option to specify the interatomic potential for a LAMMPS simulation is by defining the `potential_dataframe`\ndirectly: ","metadata":{},"id":"ed5732c4-4221-4081-83d8-2c1df793d8b7"},{"cell_type":"code","source":"import pandas \n\npotential_dataframe = pandas.DataFrame({\n \"Config\": [[\n \"pair_style morse/smooth/linear 9.0\",\n \"pair_coeff * * 0.5 1.8 2.95\"\n ]],\n \"Filename\": [[]],\n \"Model\": [\"Morse\"],\n \"Name\": [\"Morse\"],\n \"Species\": [[\"Al\"]],\n})","metadata":{"trusted":true,"tags":[]},"execution_count":7,"outputs":[],"id":"4ab54444-0d8b-42dd-8a97-c0743598a7bd"},{"cell_type":"markdown","source":"## Quantum Espresso\n[Quantum Espresso](https://www.quantum-espresso.org) - Integrated suite of Open-Source computer codes for \nelectronic-structure calculations:","metadata":{},"id":"0fbec738-da33-44bf-ab6f-b222f1d56186"},{"cell_type":"markdown","source":"from atomistics.calculators import evaluate_with_qe\n\nresult_dict = evaluate_with_qe(\n task_dict={},\n calculation_name=\"espresso\",\n working_directory=\".\",\n kpts=(3, 3, 3),\n pseudopotentials={\n \"Al\": \"Al.pbe-n-kjpaw_psl.1.0.0.UPF\"\n },\n tstress=True,\n tprnfor=True,\n ecutwfc=40.0, # kinetic energy cutoff (Ry) for wavefunctions\n conv_thr=1e-06, # Convergence threshold for selfconsistency\n diagonalization='david', \n electron_maxstep=100, # maximum number of iterations in a scf step. \n nstep=200, # number of molecular-dynamics or structural optimization steps performed in this run.\n etot_conv_thr=1e-4, # Convergence threshold on total energy (a.u) for ionic minimization\n forc_conv_thr=1e-3, # Convergence threshold on forces (a.u) for ionic minimization\n smearing='gaussian', # ordinary Gaussian spreading (Default)\n)","metadata":{"tags":[]},"id":"e7a8b26d-5547-4e69-b224-d65f38060c0a"},{"cell_type":"markdown","source":"This secondary interface for [Quantum Espresso](https://www.quantum-espresso.org) is based on the input writer from the\n[Atomic Simulation Environment (ASE)](https://wiki.fysik.dtu.dk/ase/) and the output is parsed using the [pwtools](https://elcorto.github.io/pwtools/).\nThe executable can be set using the `ASE_ESPRESSO_COMMAND` environment variable:","metadata":{},"id":"727f824e-97a2-42ae-94f2-b8080c023e24"},{"cell_type":"raw","source":"export ASE_ESPRESSO_COMMAND=\"pw.x -in PREFIX.pwi > PREFIX.pwo\"","metadata":{},"id":"a0207a31-fa83-4729-9b4d-6a025af78980"},{"cell_type":"markdown","source":"The full list of possible keyword arguments is available in the [Quantum Espresso Documentation](https://www.quantum-espresso.org/Doc/INPUT_PW.html).\nFinally, the [Standard solid-state pseudopotentials (SSSP)](https://www.materialscloud.org/discover/sssp/table/efficiency) \nfor quantum espresso are distributed via the materials cloud.","metadata":{"tags":[]},"id":"54ee1c7b-a829-49b9-8746-d57496b339e6"}]} \ No newline at end of file +{ + "cells": [ + { + "cell_type": "markdown", + "id": "dcc94f22-3be1-4cd9-9f1c-e61f7e03a1f5", + "metadata": {}, + "source": [ + "# Simulation Codes" + ] + }, + { + "cell_type": "markdown", + "id": "5ffec124-6252-4965-9267-cc771c79570f", + "metadata": {}, + "source": [ + "## ASE\n", + "At the current stage the majority of simulation codes are interfaced using the [Atomic Simulation Environment (ASE)](https://wiki.fysik.dtu.dk/ase/).\n", + "The limitation of the ASE based interfaces is that the simulation codes are only used to calculate energies, forces and\n", + "stresses, while more complex computations like structure optimization or molecular dynamics are implemented in python." + ] + }, + { + "cell_type": "markdown", + "id": "ece0e223-3b49-41d7-99ab-e0518ef01c2b", + "metadata": {}, + "source": [ + "### Abinit\n", + "[Abinit](https://www.abinit.org) - Plane wave density functional theory:" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "59904794-1beb-43f5-b6e2-f79404ca8b46", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "from ase.calculators.abinit import Abinit\n", + "from ase.units import Ry\n", + "from atomistics.calculators import evaluate_with_ase\n", + "\n", + "result_dict = evaluate_with_ase(\n", + " task_dict={},\n", + " ase_calculator=Abinit(\n", + " label='abinit_evcurve',\n", + " nbands=32,\n", + " ecut=10 * Ry,\n", + " kpts=(3, 3, 3),\n", + " toldfe=1.0e-2,\n", + " v8_legacy_format=False,\n", + " )\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "53687acd-cc66-470c-b88b-c4a2a7f4182c", + "metadata": {}, + "source": [ + "The full documentation of the corresponding interface is available on the [Atomic Simulation Environment](https://wiki.fysik.dtu.dk/ase/ase/calculators/abinit.html)\n", + "website. " + ] + }, + { + "cell_type": "markdown", + "id": "a765a024-fc04-4add-b6c2-eef026e4dbd2", + "metadata": {}, + "source": [ + "### EMT\n", + "[EMT](https://wiki.fysik.dtu.dk/ase/ase/calculators/emt.html) - Effective medium theory: " + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "a7467756-1366-42c2-b595-e8c7893437ab", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "from ase.calculators.emt import EMT\n", + "from atomistics.calculators import evaluate_with_ase\n", + "\n", + "result_dict = evaluate_with_ase(\n", + " task_dict={}, \n", + " ase_calculator=EMT()\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "e635174f-be67-462f-a3aa-1d3103c31d97", + "metadata": {}, + "source": [ + "The full documentation of the corresponding interface is available on the [Atomic Simulation Environment](https://wiki.fysik.dtu.dk/ase/ase/calculators/emt.html)\n", + "website. " + ] + }, + { + "cell_type": "markdown", + "id": "af758cea-e4ad-475e-8fcc-2afbf497334d", + "metadata": {}, + "source": [ + "### GPAW\n", + "[GPAW](https://wiki.fysik.dtu.dk/gpaw/) - Density functional theory Python code based on the projector-augmented wave \n", + "method:" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "10a942f6-b4c0-4710-856f-e736a643ce17", + "metadata": { + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + " ___ ___ ___ _ _ _ \n", + " | | |_ | | | | \n", + " | | | | | . | | | | \n", + " |__ | _|___|_____| 24.1.0\n", + " |___|_| \n", + "\n", + "User: jovyan@jupyter-pyiron-2datomistics-2dha1r3gln\n", + "Date: Wed May 1 16:53:34 2024\n", + "Arch: x86_64\n", + "Pid: 648\n", + "CWD: /home/jovyan/notebooks\n", + "Python: 3.10.12\n", + "gpaw: /srv/conda/envs/notebook/lib/python3.10/site-packages/gpaw\n", + "_gpaw: /srv/conda/envs/notebook/lib/python3.10/site-packages/\n", + " _gpaw.cpython-310-x86_64-linux-gnu.so\n", + "ase: /srv/conda/envs/notebook/lib/python3.10/site-packages/ase (version 3.22.1)\n", + "numpy: /srv/conda/envs/notebook/lib/python3.10/site-packages/numpy (version 1.26.4)\n", + "scipy: /srv/conda/envs/notebook/lib/python3.10/site-packages/scipy (version 1.13.0)\n", + "libxc: 6.2.2\n", + "units: Angstrom and eV\n", + "cores: 1\n", + "OpenMP: True\n", + "OMP_NUM_THREADS: 1\n", + "\n", + "Input parameters:\n", + " kpts: [3 3 3]\n", + " mode: {ecut: 300.0,\n", + " name: pw}\n", + " xc: PBE\n", + "\n", + "Memory usage: 143.04 MiB\n", + "Date: Wed May 1 16:53:34 2024\n" + ] + } + ], + "source": [ + "from gpaw import GPAW, PW\n", + "from atomistics.calculators import evaluate_with_ase\n", + "\n", + "result_dict = evaluate_with_ase(\n", + " task_dict={}, \n", + " ase_calculator=GPAW(\n", + " xc=\"PBE\",\n", + " mode=PW(300),\n", + " kpts=(3, 3, 3)\n", + " )\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "a38a03dd-630c-4bff-a256-d0380da29db1", + "metadata": {}, + "source": [ + "The full documentation of the corresponding interface is available on the [GPAW](https://wiki.fysik.dtu.dk/gpaw/)\n", + "website." + ] + }, + { + "cell_type": "markdown", + "id": "f8a20ee4-153d-4c0b-b89f-a3ac64cfbcbc", + "metadata": {}, + "source": [ + "### Quantum Espresso \n", + "[Quantum Espresso](https://www.quantum-espresso.org) - Integrated suite of Open-Source computer codes for \n", + "electronic-structure calculations:" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "40b9384e-02a0-41d8-adff-ad4dd6316a21", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "from ase.calculators.espresso import Espresso\n", + "from atomistics.calculators import evaluate_with_ase\n", + "\n", + "result_dict = evaluate_with_ase(\n", + " task_dict={}, \n", + " ase_calculator=Espresso(\n", + " pseudopotentials={\"Al\": \"Al.pbe-n-kjpaw_psl.1.0.0.UPF\"},\n", + " tstress=True,\n", + " tprnfor=True,\n", + " kpts=(3, 3, 3),\n", + " )\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "c7f74b95-e161-4ee6-8a9b-58c11ba3ca47", + "metadata": {}, + "source": [ + "The full documentation of the corresponding interface is available on the [Atomic Simulation Environment](https://wiki.fysik.dtu.dk/ase/ase/calculators/espresso.html)\n", + "website. " + ] + }, + { + "cell_type": "markdown", + "id": "cdc1c88c-62f8-4af1-a135-9266ba1ab1ee", + "metadata": {}, + "source": [ + "### Siesta\n", + "[Siesta](https://siesta-project.org) - Electronic structure calculations and ab initio molecular dynamics:" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "5fa73a10-043f-4b46-a162-3f1da7399cba", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "import os\n", + "from ase.calculators.siesta import Siesta\n", + "from ase.units import Ry\n", + "from atomistics.calculators import evaluate_with_ase\n", + "\n", + "result_dict = evaluate_with_ase(\n", + " task_dict={}, \n", + " ase_calculator=Siesta(\n", + " label=\"siesta\",\n", + " xc=\"PBE\",\n", + " mesh_cutoff=200 * Ry,\n", + " energy_shift=0.01 * Ry,\n", + " basis_set=\"DZ\",\n", + " kpts=(5, 5, 5),\n", + " fdf_arguments={\"DM.MixingWeight\": 0.1, \"MaxSCFIterations\": 100},\n", + " pseudo_path=os.path.abspath(\"tests/static/siesta\"),\n", + " pseudo_qualifier=\"\",\n", + " )\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "ed3bba6d-8fd9-4d26-8f0d-437fb5559397", + "metadata": {}, + "source": [ + "The full documentation of the corresponding interface is available on the [Atomic Simulation Environment](https://wiki.fysik.dtu.dk/ase/ase/calculators/siesta.html)\n", + "website." + ] + }, + { + "cell_type": "markdown", + "id": "0c09855c-80ca-4fb9-9f66-458b8d145e59", + "metadata": {}, + "source": [ + "## LAMMPS\n", + "[LAMMPS](https://www.lammps.org) - Molecular Dynamics:" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "a95e1c9c-1ad1-4765-aa13-790db6971dab", + "metadata": { + "tags": [] + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/srv/conda/envs/notebook/lib/python3.10/site-packages/atomistics/calculators/lammps/potential.py:299: SettingWithCopyWarning: \n", + "A value is trying to be set on a copy of a slice from a DataFrame.\n", + "Try using .loc[row_indexer,col_indexer] = value instead\n", + "\n", + "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n", + " df_pot[\"Config\"] = config_lst\n" + ] + } + ], + "source": [ + "from ase.build import bulk\n", + "from atomistics.calculators import evaluate_with_lammps, get_potential_by_name\n", + "\n", + "structure = bulk(\"Al\", cubic=True)\n", + "potential_dataframe = get_potential_by_name(\n", + " potential_name='1999--Mishin-Y--Al--LAMMPS--ipr1',\n", + " resource_path=\"../tests/static/lammps\"\n", + ")\n", + "\n", + "result_dict = evaluate_with_lammps(\n", + " task_dict={},\n", + " potential_dataframe=potential_dataframe,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "ed5732c4-4221-4081-83d8-2c1df793d8b7", + "metadata": {}, + "source": [ + "The [LAMMPS](https://www.lammps.org) interface is based on the [pylammpsmpi](https://github.com/pyiron/pylammpsmpi)\n", + "package which couples a [LAMMPS](https://www.lammps.org) instance which is parallelized via the Message Passing Interface\n", + "(MPI) with a serial python process or jupyter notebook. The challenging part about molecular dynamics simulation is \n", + "identifying a suitable interatomic potential. \n", + "\n", + "To address this challenge the `atomistics` package is leveraging the [NIST database of interatomic potentials](https://www.ctcms.nist.gov/potentials). \n", + "It is recommended to install this database `iprpy-data` via the `conda` package manager, then the `resource_path` is\n", + "automatically set to `${CONDA_PREFIX}/share/iprpy`. Alternatively, the `resource_path` can be specified manually as an\n", + "optional parameter of the `get_potential_by_name()` function.\n", + "\n", + "In addition, the `get_potential_dataframe(structure)` function which takes an `ase.atoms.Atoms` object as input can be\n", + "used to query the [NIST database of interatomic potentials](https://www.ctcms.nist.gov/potentials) for potentials, which\n", + "include the interatomic interactions required to simulate the atomic structure defined by the `ase.atoms.Atoms` object. \n", + "It returns a `pandas.DataFrame` with all the available potentials and the `resource_path` can again be specified as \n", + "optional parameter.\n", + "\n", + "Finally, another option to specify the interatomic potential for a LAMMPS simulation is by defining the `potential_dataframe`\n", + "directly: " + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "4ab54444-0d8b-42dd-8a97-c0743598a7bd", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "import pandas \n", + "\n", + "potential_dataframe = pandas.DataFrame({\n", + " \"Config\": [[\n", + " \"pair_style morse/smooth/linear 9.0\",\n", + " \"pair_coeff * * 0.5 1.8 2.95\"\n", + " ]],\n", + " \"Filename\": [[]],\n", + " \"Model\": [\"Morse\"],\n", + " \"Name\": [\"Morse\"],\n", + " \"Species\": [[\"Al\"]],\n", + "})" + ] + }, + { + "cell_type": "markdown", + "id": "0fbec738-da33-44bf-ab6f-b222f1d56186", + "metadata": {}, + "source": [ + "## Quantum Espresso\n", + "[Quantum Espresso](https://www.quantum-espresso.org) - Integrated suite of Open-Source computer codes for \n", + "electronic-structure calculations:" + ] + }, + { + "cell_type": "raw", + "id": "e0bea89e-9e2b-47cf-ba3b-e3b86fdcb823", + "metadata": { + "tags": [] + }, + "source": [ + "from atomistics.calculators import evaluate_with_qe\n", + "\n", + "result_dict = evaluate_with_qe(\n", + " task_dict={},\n", + " calculation_name=\"espresso\",\n", + " working_directory=\".\",\n", + " kpts=(3, 3, 3),\n", + " pseudopotentials={\n", + " \"Al\": \"Al.pbe-n-kjpaw_psl.1.0.0.UPF\"\n", + " },\n", + " tstress=True,\n", + " tprnfor=True,\n", + " ecutwfc=40.0, # kinetic energy cutoff (Ry) for wavefunctions\n", + " conv_thr=1e-06, # Convergence threshold for selfconsistency\n", + " diagonalization='david', \n", + " electron_maxstep=100, # maximum number of iterations in a scf step. \n", + " nstep=200, # number of molecular-dynamics or structural optimization steps performed in this run.\n", + " etot_conv_thr=1e-4, # Convergence threshold on total energy (a.u) for ionic minimization\n", + " forc_conv_thr=1e-3, # Convergence threshold on forces (a.u) for ionic minimization\n", + " smearing='gaussian', # ordinary Gaussian spreading (Default)\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "727f824e-97a2-42ae-94f2-b8080c023e24", + "metadata": {}, + "source": [ + "This secondary interface for [Quantum Espresso](https://www.quantum-espresso.org) is based on the input writer from the\n", + "[Atomic Simulation Environment (ASE)](https://wiki.fysik.dtu.dk/ase/) and the output is parsed using the [pwtools](https://elcorto.github.io/pwtools/).\n", + "The executable can be set using the `ASE_ESPRESSO_COMMAND` environment variable:" + ] + }, + { + "cell_type": "raw", + "id": "a0207a31-fa83-4729-9b4d-6a025af78980", + "metadata": {}, + "source": [ + "export ASE_ESPRESSO_COMMAND=\"pw.x -in PREFIX.pwi > PREFIX.pwo\"" + ] + }, + { + "cell_type": "markdown", + "id": "54ee1c7b-a829-49b9-8746-d57496b339e6", + "metadata": { + "tags": [] + }, + "source": [ + "The full list of possible keyword arguments is available in the [Quantum Espresso Documentation](https://www.quantum-espresso.org/Doc/INPUT_PW.html).\n", + "Finally, the [Standard solid-state pseudopotentials (SSSP)](https://www.materialscloud.org/discover/sssp/table/efficiency) \n", + "for quantum espresso are distributed via the materials cloud." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.8" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} From 708f4e6d97ff4dd006bb40a392db36629393b41c Mon Sep 17 00:00:00 2001 From: Jan Janssen Date: Wed, 1 May 2024 14:13:47 -0600 Subject: [PATCH 6/9] Update environment.yml --- binder/environment.yml | 2 -- 1 file changed, 2 deletions(-) diff --git a/binder/environment.yml b/binder/environment.yml index 0b2b0088..4e40b3e0 100644 --- a/binder/environment.yml +++ b/binder/environment.yml @@ -13,5 +13,3 @@ dependencies: - pandas =2.2.2 - pylammpsmpi =0.2.18 - jinja2 =3.1.3 -- qe =7.2 -- pwtools =1.2.3 From f9680ceb214dee7c53ead3d46ffad16d9c0c3475 Mon Sep 17 00:00:00 2001 From: Jan Janssen Date: Wed, 1 May 2024 14:19:22 -0600 Subject: [PATCH 7/9] Update build_notebooks.sh --- .ci_support/build_notebooks.sh | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.ci_support/build_notebooks.sh b/.ci_support/build_notebooks.sh index 81a89d84..f8d36f2f 100755 --- a/.ci_support/build_notebooks.sh +++ b/.ci_support/build_notebooks.sh @@ -1,10 +1,10 @@ # execute notebooks i=0; -for notebook in $(ls notebooks/*.ipynb); do +for notebook in $(ls *.ipynb); do papermill ${notebook} ${notebook%.*}-out.${notebook##*.} -k "python3" || i=$((i+1)); done; # push error to next level if [ $i -gt 0 ]; then exit 1; -fi; \ No newline at end of file +fi; From 2042c5d4aaa00babc89774c15c7a118a507a427b Mon Sep 17 00:00:00 2001 From: Jan Janssen Date: Wed, 1 May 2024 14:19:45 -0600 Subject: [PATCH 8/9] Update notebooks.yml --- .github/workflows/notebooks.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/notebooks.yml b/.github/workflows/notebooks.yml index 643f39dc..e153219b 100644 --- a/.github/workflows/notebooks.yml +++ b/.github/workflows/notebooks.yml @@ -31,4 +31,5 @@ jobs: run: | pip install versioneer[toml]==0.29 pip install . --no-deps --no-build-isolation + cd notebooks ./.ci_support/build_notebooks.sh From 55fbefde10b6993fbf579fcf84cb02cd91cac0e4 Mon Sep 17 00:00:00 2001 From: Jan Janssen Date: Wed, 1 May 2024 14:24:53 -0600 Subject: [PATCH 9/9] Update notebooks.yml --- .github/workflows/notebooks.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/notebooks.yml b/.github/workflows/notebooks.yml index e153219b..18e490db 100644 --- a/.github/workflows/notebooks.yml +++ b/.github/workflows/notebooks.yml @@ -32,4 +32,4 @@ jobs: pip install versioneer[toml]==0.29 pip install . --no-deps --no-build-isolation cd notebooks - ./.ci_support/build_notebooks.sh + ../.ci_support/build_notebooks.sh