Introductory PyBigDFT tour - Basic functionalities and first calculations
This material is adapted from work by Luigi Genovese, Laura Ratcliff and others
We start here by presenting a walkthrough overview of using the PyBigDFT library to set up and run BigDFT calculations.
The main documentation of the BigDFT code can be found in the User Page, and informations about the PyBigDFT API can be found here
For this tutorial, we have to ensure at least the installation of the bigdft client.
The basic installation can be performed with those commands:
!pip install pybigdft py3dmol pyfutile
Collecting pybigdft
Using cached pybigdft-1.0.7-py3-none-any.whl.metadata (10.0 kB)
Collecting py3dmol
Using cached py3Dmol-2.0.4-py2.py3-none-any.whl.metadata (2.1 kB)
Collecting pyfutile
Using cached pyfutile-1.9.3-py3-none-any.whl (56 kB)
Requirement already satisfied: numpy in /opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages (from pybigdft) (1.24.4)
Requirement already satisfied: pyyaml in /opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages (from pybigdft) (6.0.1)
Collecting scipy (from pybigdft)
Using cached scipy-1.10.1-cp38-cp38-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (34.5 MB)
Using cached pybigdft-1.0.7-py3-none-any.whl (361 kB)
Using cached py3Dmol-2.0.4-py2.py3-none-any.whl (12 kB)
Installing collected packages: py3dmol, scipy, pyfutile, pybigdft
Successfully installed py3dmol-2.0.4 pybigdft-1.0.7 pyfutile-1.9.3 scipy-1.10.1
System Generation
In this first overview, we will introduce the data structures available in PyBigDFT for creating and manipulating the systems we want to study (whether molecular or solid state).
1) Atom Class
Any system we want to study is going to be made up of atoms. Considering the features available in Python, the best way to store information about atoms is a dictonary (dict), which are serializable in human readable formats like json or yaml. For example, in BigDFT, we use yaml dictonaries.
from yaml import dump
datm = {"sym": "H", "r": [1, 0, 0], "units": "angstroem"}
print(dump(datm))
r:
- 1
- 0
- 0
sym: H
units: angstroem
Nonetheless, just manipulating a dict by itself is error prone, and you might want some extra helpful subroutines. For these reasons, we’ve wrapped up a dict in the Atom class.
from BigDFT.Atoms import Atom
atom = Atom(datm)
# the following also works
atom = Atom({"sym": "H", "r": [1, 0, 0], "units": "angstroem"})
atom = Atom({"H": [1, 0, 0], "units": "angstroem"})
print(dump(atom))
!!python/object:BigDFT.Atoms.Atom
store:
H:
- 1
- 0
- 0
units: angstroem
Some of in-built subroutines for Atom class manipulation are demonstrated below.
print(atom.sym)
print(atom.atomic_number)
print(atom.get_position("angstroem"))
print(atom.get_position("bohr"))
H
1
[1.0, 0.0, 0.0]
[1.8897261245650618, 0.0, 0.0]
With this approach, we nonetheless retain the flexibility of a dict.
atom["source"] = "tutorial"
print(atom["source"])
for k, v in atom.items():
print(k, v)
tutorial
H [1, 0, 0]
units angstroem
source tutorial
2) Fragment Class
We won’t do many calculations involving a single atom, instead we want to put together groups of atoms. In this case, we will use a list as our model data structure, with the wrapper class refered to as a Fragment.
atm1 = Atom({"sym": "O", "r": [2.3229430273, 1.3229430273, 1.7139430273], "units": "angstroem"})
atm2 = Atom({"sym": "H", "r": [2.3229430273, 2.0801430273, 1.1274430273], "units": "angstroem"})
atm3 = Atom({"sym": "H", "r": [2.3229430273, 0.5657430273000001, 1.1274430273], "units": "angstroem"})
from BigDFT.Fragments import Fragment
frag1 = Fragment([atm1, atm2, atm3])
print(len(frag1))
print(frag1.centroid)
3
[4.389726124565062, 2.5000000000000004, 2.4999999999999996]
It’s also possible to build up a fragment in a more step by step process.
frag1 = Fragment()
frag1.append(atm1)
frag1 += Fragment([atm2])
frag1.extend(Fragment([atm3]))
We added features to translate and rotate a fragment.
from copy import deepcopy
frag2 = deepcopy(frag1)
frag2.translate([10, 0, 0])
frag2.rotate(x=90, units="degrees")
print(dump(frag2))
!!python/object:BigDFT.Fragments.Fragment
atoms:
- !!python/object:BigDFT.Atoms.Atom
store:
r:
- 14.389726124565062
- 1.7611170852950608
- 2.499999999999999
sym: O
units: bohr
- !!python/object:BigDFT.Atoms.Atom
store:
r:
- 14.389726124565062
- 2.86944145735247
- 3.930900621520664
sym: H
units: bohr
- !!python/object:BigDFT.Atoms.Atom
store:
r:
- 14.389726124565062
- 2.86944145735247
- 1.0690993784793343
sym: H
units: bohr
conmat: null
frozen: null
q1: null
q2: null
3) System Class
Many electronic structure codes top off at the list of atoms level, but in PyBigDFT we go one step further. At the top, we have the System class which is again based on a dict.
from BigDFT.Systems import System
sys = System()
sys["WAT:0"] = frag1
sys["WAT:1"] = frag2
print(dump(sys))
!!python/object:BigDFT.Systems.System
cell: !!python/object:BigDFT.UnitCells.UnitCell
cell:
- - .inf
- 0
- 0
- - 0
- .inf
- 0
- - 0
- 0
- .inf
conmat: null
store:
WAT:0: !!python/object:BigDFT.Fragments.Fragment
atoms:
- !!python/object:BigDFT.Atoms.Atom
store:
r:
- 2.3229430273
- 1.3229430273
- 1.7139430273
sym: O
units: angstroem
- !!python/object:BigDFT.Atoms.Atom
store:
r:
- 2.3229430273
- 2.0801430273
- 1.1274430273
sym: H
units: angstroem
- !!python/object:BigDFT.Atoms.Atom
store:
r:
- 2.3229430273
- 0.5657430273000001
- 1.1274430273
sym: H
units: angstroem
conmat: null
frozen: null
q1: null
q2: null
WAT:1: !!python/object:BigDFT.Fragments.Fragment
atoms:
- !!python/object:BigDFT.Atoms.Atom
store:
r:
- 14.389726124565062
- 1.7611170852950608
- 2.499999999999999
sym: O
units: bohr
- !!python/object:BigDFT.Atoms.Atom
store:
r:
- 14.389726124565062
- 2.86944145735247
- 3.930900621520664
sym: H
units: bohr
- !!python/object:BigDFT.Atoms.Atom
store:
r:
- 14.389726124565062
- 2.86944145735247
- 1.0690993784793343
sym: H
units: bohr
conmat: null
frozen: null
q1: null
q2: null
In principle, any dictionary key is fine to use for our System class, but in practice we follow the convention of giving it a name and identifier separated by a colon. To summarize the hierarchy, let’s iterate over our System.
for fragid, frag in sys.items():
print(fragid)
for atm in frag:
print(dict(atm))
WAT:0
{'sym': 'O', 'r': [2.3229430273, 1.3229430273, 1.7139430273], 'units': 'angstroem'}
{'sym': 'H', 'r': [2.3229430273, 2.0801430273, 1.1274430273], 'units': 'angstroem'}
{'sym': 'H', 'r': [2.3229430273, 0.5657430273000001, 1.1274430273], 'units': 'angstroem'}
WAT:1
{'sym': 'O', 'r': [14.389726124565062, 1.7611170852950608, 2.499999999999999], 'units': 'bohr'}
{'sym': 'H', 'r': [14.389726124565062, 2.86944145735247, 3.930900621520664], 'units': 'bohr'}
{'sym': 'H', 'r': [14.389726124565062, 2.86944145735247, 1.0690993784793343], 'units': 'bohr'}
Now that we’ve reached the top level, let’s visualize the system we have built.
_ = sys.display()
You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol
The visualization module has identified that there are two separate fragments, and colored them accordingly. Of course if we merged our fragments, the visualization would look different.
sys2 = System()
sys2["COM:0"] = sum(sys.values())
_ = sys2.display()
You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol
Multiple Views
We can take advantage of deepcopying to construct multiple views of the same system. For example, we might want to have two separate views of the same set of atoms. In one view, we split the set into two molecules, and the other we have just one big fragment. This might be convenient if, for example, we want to be able to rotate the entire system as a group.
#Reminder of standard shallow copy
a = {"x": "1"}
b = {"x": "2"}
my_list = [a, b]
print(my_list)
a["x"] = 3
print(my_list)
[{'x': '1'}, {'x': '2'}]
[{'x': 3}, {'x': '2'}]
#Deep copy
sep = deepcopy(sys)
joint = System()
joint["COM:0"] = sum(sep.values())
joint["COM:0"].rotate(y=90, units="degrees")
_ = sep.display()
You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol
_ = joint.display()
You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol
Extended Systems
For building and manipulating extended systems, the extra ingredient required is a UnitCell.
from BigDFT.UnitCells import UnitCell
sys.cell = UnitCell([7.0, 7.0, 7.0], units="angstroem")
_ = sys.display()
You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol
We see that the minimum image convention has wrapped the red fragment around so that it is now on the left side. We can inspect this position values closer from the Atom class.
for fragid, frag in sys.items():
print(fragid)
for at in frag:
print(at.get_position("angstroem"), at.get_position("angstroem", cell=sys.cell))
WAT:0
[2.3229430273, 1.3229430273, 1.7139430273] [2.3229430273, 1.3229430273, 1.7139430273]
[2.3229430273, 2.0801430273, 1.1274430273] [2.3229430273, 2.0801430273, 1.1274430273]
[2.3229430273, 0.5657430273000001, 1.1274430273] [2.3229430273, 0.5657430273000001, 1.1274430273]
WAT:1
[7.6147151365, 0.9319430273000001, 1.3229430272999996] [0.6147151365000005, 0.9319430273000001, 1.3229430272999996]
[7.6147151365, 1.5184430273000002, 2.0801430272999997] [0.6147151365000005, 1.5184430273000002, 2.0801430272999997]
[7.6147151365, 1.5184430273000002, 0.5657430272999997] [0.6147151365000005, 1.5184430273000002, 0.5657430272999997]
We also can get accessed to fractional units this way.
for fragid, frag in sys.items():
print(fragid)
for at in frag:
print(at.get_position("reduced", cell=sys.cell))
WAT:0
[0.331849, 0.188992, 0.244849]
[0.331849, 0.297163, 0.161063]
[0.331849, 0.08082, 0.161063]
WAT:1
[0.087816, 0.133135, 0.188992]
[0.087816, 0.21692, 0.297163]
[0.087816, 0.21692, 0.08082]
The sys.get_posinp() method provides the information written in YAML markup format:
sys.get_posinp()
{'positions': [{'frag': ['WAT', '0'],
'O': [2.3229430273, 1.3229430273, 1.7139430273]},
{'frag': ['WAT', '0'], 'H': [2.3229430273, 2.0801430273, 1.1274430273]},
{'frag': ['WAT', '0'],
'H': [2.3229430273, 0.5657430273000001, 1.1274430273]},
{'frag': ['WAT', '1'],
'O': [0.6147151365000005, 0.9319430273000001, 1.3229430272999996]},
{'frag': ['WAT', '1'],
'H': [0.6147151365000005, 1.5184430273000002, 2.0801430272999997]},
{'frag': ['WAT', '1'],
'H': [0.6147151365000005, 1.5184430273000002, 0.5657430272999997]}],
'units': 'angstroem',
'cell': [7.0, 7.0, 7.0]}
Calculation examples
1 - Simple finite system
2 - Periodic system
We now start to perform some basic calculations with the code. Such calculation can be executed:
A. On the local workstation if you have the bigdft code installed.
B. On google colab with the Installation method “full-suite” (see the 0.Installation.ipynb) notebook of this project.
C. We can use the dry-run approach where we do not actually run the calculation but we employ previously existing logfiles of the same systems.
We here employ the solution C. In order to do that we have to compy the data from the bigdft-school database repo.
!wget https://raw.githubusercontent.com/BigDFT-group/bigdft-school/main/packaging/install.py 2> /dev/null
import install
install.data('data/Enccs_basics.tar.xz')
Executing: wget https://gitlab.com/luigigenovese/bigdft-school/-/raw/main/data/Enccs_basics.tar.xz -O Enccs_basics.tar.xz
--2023-12-07 22:24:12-- https://gitlab.com/luigigenovese/bigdft-school/-/raw/main/data/Enccs_basics.tar.xz
Resolving gitlab.com (gitlab.com)... 172.65.251.78, 2606:4700:90:0:f22e:fbec:5bed:a9b9
Connecting to gitlab.com (gitlab.com)|172.65.251.78|:443... connected.
HTTP request sent, awaiting response... 200 OK
Length: 22228 (22K) [application/octet-stream]
Saving to: ‘Enccs_basics.tar.xz’
0K .......... .......... . 100% 104M=0s
2023-12-07 22:24:13 (104 MB/s) - ‘Enccs_basics.tar.xz’ saved [22228/22228]
Executing: mkdir -p /home/runner/work/efficient-materials-modelling-on-hpc/efficient-materials-modelling-on-hpc/content/code/bigdft
Executing: tar xJf Enccs_basics.tar.xz -C .
1 - Let us create a simple system using the tools described so far and run a calculation with a wavelet grid spacing of $0.4$ bohr (good for the large majority of the systems) and the PBE functional. A logfile is generated from which we can access the computed properties. This logfile has built in properties and can be accessed like a dictionary.
# Create Three Atoms
at1 = Atom({"H": [0, 0, 0], "units": "bohr"})
at2 = Atom({"H": [0, 0, 1.4], "units": "bohr"})
at3 = Atom({"He": [10, 0, 0], "units": "bohr"})
# Construct a System from Two Fragments (H2, He)
sys = System()
sys["H2:1"] = Fragment([at1, at2])
sys["He:2"] = Fragment([at3])
# Iterate Over The System to inspect the atoms positions
for at in sys.get_atoms():
print(at.sym, at.get_position())
# or alternatively, two-level iteration (fragment first)
for fragid, frag in sys.items():
print(fragid)
for at in frag:
print(' ', at.sym, at.get_position())
H [0.0, 0.0, 0.0]
H [0.0, 0.0, 1.4]
He [10.0, 0.0, 0.0]
H2:1
H [0.0, 0.0, 0.0]
H [0.0, 0.0, 1.4]
He:2
He [10.0, 0.0, 0.0]
# this is what we have created
sys.display()
You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol
<BigDFT.Visualization.InlineVisualizer at 0x7fd6a0f72e50>
from BigDFT.Inputfiles import Inputfile
inp = Inputfile()
inp.set_hgrid(0.4)
inp.set_xc("PBE")
# we remove the calculation of the forces
inp["perf"] = {"calculate_forces": False}
from BigDFT.Calculators import SystemCalculator
calc = SystemCalculator(skip=True, verbose=False)
log = calc.run(sys=sys, input=inp, name="quick", run_dir="scratch")
#some basic info
print(log.energy)
print(list(log.evals))
print(log.log["Memory Consumption Report"]
["Memory occupation"])
-4.0546151576846
[BandArray([[-0.57600357, -0.38118857]])]
{'Peak Value (MB)': 230.69, 'for the array': 'f_i', 'in the routine': 'vxcpostprocessing', 'Memory Peak of process': 'unknown'}
We have “executed” (the calculation result is provided with the school material here) a very simple calculation with two Kohn-Sham orbitals. The description of the possible input parameters and how to trigger larger calculations (even on remote computers) will be provided in the forthcoming tutorials.
2 - We setup a BCC unit cell of iron and perform the calculation using a 2x2x2 k-point grid with a Monkhorst-Pack grid.
from BigDFT.UnitCells import UnitCell
# one single periodic atom
pat = Atom({"Fe": [0, 0, 0], "units": "angstroem"})
psys = System({"CEL:0": Fragment([pat])})
psys.cell = UnitCell([2.867, 2.867, 2.867], units="angstroem")
psys.display()
You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol
<BigDFT.Visualization.InlineVisualizer at 0x7fd6a0edaf10>
# very small k-point mesh, just for comparison
inp = Inputfile()
inp.set_hgrid(0.3)
inp.set_xc("LDA") # can be omitted as this is the default
inp["kpt"] = {"method": "mpgrid", "ngkpt": [2, 2, 2]}
log = calc.run(sys=psys, input=inp, name="psys", run_dir="scratch")
_ = log.get_dos().plot()
