import numpy as np
import matplotlib.pyplot as plt
#Customize default plotting style
%matplotlib inline
import seaborn as sns
sns.set_context('talk')
plt.rcParams["figure.figsize"] = (10, 8)
spglib
, phonopy
) using minimal ASE-compatible APIf90wrap
adds derived type support to f2py
¶f90wrap
package addresses this by generating an additional layer of wrappers, adding support for derived types, module data, efficient array access, Python 2.6+ and 3.xf90wrap
schematic overview¶More details: J. R. Kermode, f90wrap: an automated tool for constructing deep Python interfaces to modern Fortran codes. J. Phys. Condens. Matter (2020) doi:10.1088/1361-648X/ab82d2
bader
code¶f90wrap
to automatically generate a deep Python interface with very little manual workGeneration and compilation of wrappers:
f90wrap -v -k kind_map -I init.py -m bader \
kind_mod.f90 matrix_mod.f90 \
ions_mod.f90 options_mod.f90 charge_mod.f90 \
chgcar_mod.f90 cube_mod.f90 io_mod.f90 \
bader_mod.f90 voronoi_mod.f90 multipole_mod.f90
f2py-f90wrap -c -m _bader f90wrap_*.f90 -L. -lbader
To install it yourself, run the commands below, taking care to adjust PY_INSTALL_DIR
according to your local setup:
git clone https://gitlab.com/jameskermode/bader
cd bader
export PY_INSTALL_DIR=~/.local/lib/python3.8/site-packages
make -f makefile.osx_gfortran python
bader
code (contd.)¶Restart a gpaw
DFT calculation (or run if necessary) and retrieve the density:
import os
from ase.build import bulk
from gpaw import GPAW, restart
if not os.path.exists('si-vac.gpw'):
si = bulk('Si', cubic=True)
del si[0] # create a vacancy
gpaw = GPAW(h=0.15)
si.set_calculator(gpaw)
si.get_potential_energy()
gpaw.write('si-vac.gpw')
si, gpaw = restart('si-vac.gpw')
rho = gpaw.get_pseudo_density()
plt.plot(si.positions[:, 0], si.positions[:, 1], 'r.', ms=50)
plt.imshow(rho[:,:,0], extent=[0, si.cell[0,0], 0, si.cell[1,1]]);
import bader
bdr = bader.bader(si, rho)
print('ionchg:', bdr.ionchg)
ionchg: [4.2901335 4.43094328 4.28875689 4.43093625 4.32306902 4.43071537 4.3392679 ] CALCULATING BADER CHARGE DISTRIBUTION 0 10 25 50 75 100 PERCENT DONE: ********************** REFINING AUTOMATICALLY ITERATION: 2 CHECKED POINTS: 0 REASSIGNED POINTS: 0 RUN TIME: 0.05 SECONDS CALCULATING MINIMUM DISTANCES TO ATOMS 0 10 25 50 75 100 PERCENT DONE: ********************** RUN TIME: 0.01 SECONDS
## TAB-complete to introspect inside `bdr` Fortran type
dtype('float64')
# collect Bader volumes associated with atom #5
atom = 4
rho3 = np.zeros_like(rho)
for v in (bdr.nnion == atom+1).nonzero()[0]:
rho3[bdr.volnum == v+1] = rho[bdr.volnum == v+1]
# write a CUBE file to allow visualisation of density (FIXME: can this be avoided?)
from ase.io.cube import write_cube
with open('rho3.cube', 'w') as f:
write_cube(f, si, rho3)
import nglview
v = nglview.show_ase(si, gui=True)
v.add_representation('unitcell')
v.add_component('rho3.cube')
v.component_1.update_surface()
v
NGLWidget()
Tab(children=(Box(children=(Box(children=(Box(children=(Label(value='step'), IntSlider(value=1, min=-100)), la…
f90wrap
can now wrap large codes like Castep to provide deep access to internal dataCasPyTep
in 2016, adding MPI support and optimising performance of Nudged Elastic Band algorithm, coded in PythonCastep_Python
module with a high level API, also built using f90wrap
gfortran
and ifort
)pip install numpy
pip install f90wrap
pip install ase
In principle:
cd castep
make
make python
export PY_INSTALL_DIR=/usr/local/python3.7/site-packages # or somewhere else
make python-install
For the high level bindings only - see README.INSTALL
for more details:
make castep_python
pip install -e obj/[ARCH]
Set of source files currently wrapped is as follows, can be easily expanded:
Utility: constants.F90 algor.F90 comms.serial.F90
io.F90 trace.F90 license.F90 buildinfo.f90
Fundamental: parameters.F90 cell.F90 basis.F90
ion.F90 density.F90 wave.F90
Functional: model.F90 electronic.F90 firstd.f90 xc.f90
Already far too much to wrap by hand:
current_cell
, etc (NB: must have target
attribute)unit_cell
, model_state
), including all elements, arrays, etc. within themcurrent_cell%ionic_positions
exposed directly in Pythonio_abort()
raises RuntimeError
exceptionCasPyTep(atoms)
caspytep
for a test drive¶import numpy as np
import caspytep
## uncomment line below and press TAB for autocompletion
caspytep.
#caspytep.cell.unit_cell.
## append a ? to access documentation strings
caspytep.model.model_wave_read?
#caspytep.cell.cell_read?
This uses the ASE-compatible interface provided by the CasPyTep
class.
from ase.build import bulk
atoms = bulk('Si', cubic=True)
calc = caspytep.calculator.CasPyTep(atoms=atoms)
atoms.set_calculator(calc)
e = atoms.get_potential_energy()
f = atoms.get_forces()
print(f'energy: {e:.3f} eV')
print(f'forces: {f}')
energy: -1351.306 eV forces: [[-1.77531190e-06 7.79031704e-06 -7.78323921e-06] [ 5.28965998e-05 -1.26374348e-05 -2.42809706e-05] [-4.99546145e-06 1.06564617e-05 5.72414230e-05] [-1.50490900e-05 -1.54993264e-05 -1.93367246e-06] [ 9.94049519e-06 3.30692511e-05 -7.02829293e-06] [-1.56277562e-05 7.36665153e-06 -1.61929791e-05] [ 1.36706212e-05 -1.19693116e-05 3.73698167e-05] [-3.90600967e-05 -1.87766085e-05 -3.73920853e-05]]
Unlike with standard ASE or other scripting approaches, after running a calculation, we can now poke around in all the internal arrays:
#calc.model.eigenvalues
#calc.model.cell.ionic_positions
#calc.model.cell.ionic_positions[...,0]
#calc.model.wvfn.beta_phi
#calc.model.cell.ionic_velocities.T
from ase.units import Hartree
calc.parameters.cut_off_energy
190.4797021705707
from ase.units import Bohr
p = calc.model.cell.ionic_positions.copy()
p = p[:, :, 0] # first species only
p = calc.model.cell.real_lattice.T @ p
xi, yi, zi = p * Bohr
plt.scatter(xi, yi, s=200, c='r')
plt.axis([0, atoms.cell[0,0], 0, atoms.cell[1,1]]); plt.axis("scaled");
# overlay the charge density
plt.scatter(xi, yi, s=200, c='r')
den = calc.model.den.real_charge.copy()
basis = caspytep.basis.get_current_basis()
den3 = (den.reshape((basis.ngx, basis.ngy,
basis.ngz), order='F') /
basis.total_grid_points)
plt.imshow(den3[:, :, basis.ngz//2],
extent=[0, atoms.cell[0,0], 0, atoms.cell[1,1]]);
.check
files etc.).from ase.build import bulk
from ase.optimize import LBFGS
from caspytep.calculator import CasPyTep
atoms = bulk("Si", cubic=True)
calc = CasPyTep(atoms=atoms)
atoms.set_calculator(calc)
atoms.rattle(0.01)
a0 = atoms.copy()
opt = LBFGS(atoms)
opt.run(fmax=0.1)
print(atoms.get_potential_energy())
Step Time Energy fmax LBFGS: 0 15:45:59 -1351.292226 0.3984 LBFGS: 1 15:46:02 -1351.299107 0.2784 LBFGS: 2 15:46:05 -1351.306328 0.0425 -1351.3063283114898
Python interface makes it quick to try out new high-level algorithms or connect things in new ways, e.g. testing preconditioned geometry optimizer [Packwood2016]
from ase.optimize.precon import PreconLBFGS
from caspytep.calculator import CasPyTep
atoms = a0.copy() # restart from same randomised positions as above
atoms.set_calculator(CasPyTep(atoms=atoms))
opt = PreconLBFGS(atoms, precon='Exp')
opt.run(fmax=0.05)
PreconLBFGS: 0 15:46:42 -1351.292226 0.3984 PreconLBFGS: 1 15:46:51 -1351.306183 0.0593 PreconLBFGS: 2 15:46:54 -1351.306486 0.0182
True
This is an example of using the native CasPyTep interface directly rather than the ASE compatibility layer. We increase the plane wave cutoff energy in steps of 10% until energy changes by less than $10^{-3}$ Hartree. (This isn't necessarily the best way to do convergence testing...)
from ase.build import bulk
from caspytep.calculator import CasPyTep
calc = CasPyTep(atoms=bulk("Si")) # 2-atom Si system
energy_tol = 1e-4
current_params = caspytep.parameters.get_current_params()
current_params.cut_off_energy = 7.0
cutoffs = []
energy = []
while True:
caspytep.basis.basis_initialise(current_params.cut_off_energy)
current_params.fine_gmax = (current_params.fine_grid_scale *
np.sqrt(2.0*current_params.cut_off_energy))
caspytep.ion.ion_real_initialise()
model = caspytep.model.model_state()
model.converged = caspytep.electronic.electronic_minimisation(model)
current_params.cut_off_energy *= 1.1
print('cutoff %.2f energy %.5f' % (current_params.cut_off_energy,
model.total_energy))
cutoffs.append(current_params.cut_off_energy)
energy.append(model.total_energy)
if len(energy) > 2 and abs(energy[-1] - energy[-2]) < energy_tol:
print('converged at cutoff', cutoffs[-1])
break
cutoff 7.70 energy -12.42124 cutoff 8.47 energy -12.42233 cutoff 9.32 energy -12.42250 cutoff 10.25 energy -12.42281 cutoff 11.27 energy -12.42289 converged at cutoff 11.273570000000007
from ase.units import Hartree
ecut = np.array(cutoffs) * Hartree
ediff = np.array(energy) * Hartree
ediff -= ediff[-1]
plt.plot(ecut, abs(ediff) * 1e3 / caspytep.cell.current_cell.num_ions, 'o-')
plt.xlabel('Cutoff / eV')
plt.ylabel('Energy Error / meV')
plt.axhline(1, linestyle='--');
current_cell
to be updated in place without having to call cell_read()
:cell_read()
into cell_read()
and cell_initialise()
mpirun -np N python script.py
works if CasPyTep compiled with MPI librariesDownload the image and have a play with CasPyTep in the practical session.
Install Docker CE (free download) on your Linux/Mac/Windows machine (root access required), then:
gzcat caspytep.tar.gz | docker load
docker run -v ~:/root/host -p 8899:8899 caspytep
Then point your web brower at http://localhost:8899 and browse to noteboooks > demo.ipynb to open this notebook.
To open a bash
shell instead:
docker run -v ~:/root/host -it caspytep /bin/bash
$ python
>>> import caspytep
Primary:
Secondary:
Longer term benefits: