API

Pure-Python implementation

pymbd.atomic_polarizabilities(coords, species, volumes, beta, lattice=None, kind='TS')

Calculate static atomic polarizabilities.

Parameters:

• coords (array-like) – (a.u.) atom coordinates in rows

• species (array-like) – atom types (elements)

• volumes (array-like) – ratios of Hirshfeld volumes in molecule and vacuum

• beta (float) – MBD damping parameter \(\beta\) (see J. Chem. Phys. 140, 18A508 (2014))

• lattice (array-like) – (a.u.) lattice vectors in rows

• kind (str) – one of ‘TS’, ‘BG’ or ‘TSsurf’

Returns the isotropic per-atom screened static polarizabilities of shape (N,) (a.u.).

pymbd.from_volumes(species, volumes, kind='TS')

Scale free-atom vdW parameters by Hirshfeld volume ratios.

Parameters:

• species (array-like) – atom types (elements)

• volumes (array-like) – ratios of Hirshfeld volumes in molecule and vacuum

• kind (str) – one of ‘TS’, ‘BG’ or ‘TSsurf’

Returns partitioned atomic polarizabilities, \(C_6\) coefficients, and \(R_\mathrm{vdw}\) coefficients (a.u.).

pymbd.mbd_energy(coords, alpha_0, C6, R_vdw, beta, lattice=None, k_grid=None, nfreq=15)

Calculate an MBD energy.

Use with Hirshfeld-partitioned atomic polarizabilities, C6 coefficients and van der Waal radii. Otherwise use mbd_energy_species().

Parameters:

• coords (array-like) – (a.u.) atomic coordinates in rows

• alpha_0 (array-like) – (a.u.) atomic polarizabilities

• C6 (array-like) – (a.u.) atomic \(C_6\) coefficients

• R_vdw (array-like) – (a.u.) atomic vdW radii

• beta (float) – MBD damping parameter \(\beta\)

• lattice (array-like) – (a.u.) lattice vectors in rows

• k_grid (array-like) – number of \(k\)-points along reciprocal axes

• nfreq (int) – number of grid points for frequency quadrature

pymbd.mbd_energy_species(coords, species, volume_ratios, beta, **kwargs)

Calculate an MBD energy from atom types and Hirshfeld-volume ratios.

Parameters:

• coords (array-like) – (a.u.) atom coordinates in rows

• species (array-like) – atom types (elements)

• volume_ratios (array-like) – ratios of Hirshfeld volumes in molecule and vacuum

• beta (float) – MBD damping parameter \(\beta\) (use :math:’beta = 0.83’ for PBE xc-functional)

• kwargs – see mbd_energy()

pymbd.molecular_polarizability(coords, species, volumes, beta, lattice=None, kind='TS')

Calculate the static molecular polarizability tensor.

Parameters:

• coords (array-like) – (a.u.) atom coordinates in rows

• species (array-like) – atom types (elements)

• volumes (array-like) – ratios of Hirshfeld volumes in molecule and vacuum

• beta (float) – MBD damping parameter \(\beta\) (see J. Chem. Phys. 140, 18A508 (2014))

• lattice (array-like) – (a.u.) lattice vectors in rows

• kind (str) – one of ‘TS’, ‘BG’ or ‘TSsurf’

Returns the \(3\times3\) static polarizability tensor (a.u.).

pymbd.screening(coords, alpha_0, C6, R_vdw, beta, lattice=None, nfreq=15)

Screen atomic polarizabilities.

Parameters:

• coords (array-like) – (a.u.) atom coordinates in rows

• alpha_0 (array-like) – (a.u.) atomic polarizabilities

• C6 (array-like) – (a.u.) atomic \(C_6\) coefficients

• R_vdw (array-like) – (a.u.) atomic vdW radii

• beta (float) – MBD damping parameter \(\beta\)

• lattice (array-like) – (a.u.) lattice vectors in rows

• nfreq (int) – number of grid points for frequency quadrature

Returns static polarizabilities, \(C_6\) coefficients, and \(R_\mathrm{vdw}\) coefficients (a.u.).

pymbd.screening_matrix(coords, alpha_0, R_vdw, beta, lattice=None)

Build the screened non-local polarizability matrix.

Returns the \((3N, 3N)\) matrix \(\bar{\boldsymbol\alpha}=(\boldsymbol\alpha_0^{-1}+\mathbf T_\mathrm{GG})^{-1}\), whose \(3\times3\) blocks \(\bar{\boldsymbol\alpha}_{ij}\) are the dipole-screened polarizabilities coupling atoms \(i\) and \(j\). This is the shared core of screening(), atomic_polarizabilities(), and molecular_polarizability().

Parameters:

• coords (array-like) – (a.u.) atom coordinates in rows

• alpha_0 (array-like) – (a.u.) atomic polarizabilities

• R_vdw (array-like) – (a.u.) atomic vdW radii

• beta (float) – MBD damping parameter \(\beta\)

• lattice (array-like) – (a.u.) lattice vectors in rows

pymbd.pymbd.ang = 1.8897259885789233

(a.u.) angstrom

Can be imported directly as pymbd.ang.

PySCF bridge

Bridge between PySCF and pyMBD via Hirshfeld volume rescaling.

This optional module turns a converged PySCF Kohn–Sham object into an MBD (or TS) dispersion energy. It computes per-atom Hirshfeld volume ratios

\[v_A = \frac{\int w_A(\mathbf r)\,|\mathbf r-\mathbf R_A|^3\, \rho(\mathbf r)\,d\mathbf r} {\int |\mathbf r-\mathbf R_A|^3\, \rho_A^\mathrm{free}(\mathbf r)\,d\mathbf r}, \qquad w_A = \frac{\rho_A^\mathrm{free}}{\sum_B \rho_B^\mathrm{free}},\]

directly from the molecular density and a promolecule of spherically averaged free-atom densities, and rescales the free-atom vdW parameters (\(\alpha\propto v\), \(C_6\propto v^2\), \(R_\mathrm{vdw}\propto v^{1/3}\)) before calling pymbd.fortran.MBDGeom.

The free-atom reference is a spherically averaged Kohn–Sham atom (AtomSphAverageRKS) computed with the same functional and the same basis as the molecule. Sharing the basis is what makes the ratio a pure measure of the chemical environment: the reference then cancels the basis-set incompleteness that is also present in the molecular density, so well-separated atoms come out at exactly 1 in any basis. For a one-electron atom (hydrogen) the reference is the diffuse, self-consistent DFT density, matching the FHI-aims/libMBD convention. MBD is only meaningful on top of a (semi)local/hybrid KS functional, so a Hartree–Fock mean-field is rejected.

These choices are validated point by point against independent all-electron FHI-aims MBD energies over S66x8; devtools/pyscf_s66x8_figshare_check.py reproduces that comparison.

Requires PySCF (pip install pymbd[pyscf]):

from pyscf import dft, gto

from pymbd.pyscf import mbd_energy

mol = gto.M(atom='He 0 0 0; He 0 0 3', basis='def2-tzvp')
mf = dft.RKS(mol, xc='PBE').run()
ene = mbd_energy(mf)  # beta is taken from the functional

class pymbd.pyscf.AtomSphAverageRKS(*args: Any, **kwargs: Any)

Spherically averaged neutral free-atom solver with a Kohn–Sham potential.

Composed the way PySCF builds RKS = KohnShamDFT + RHF, but with the spherical-average eig/get_occ of pyscf.scf.atom_hf. AtomSphAverageRHF, so the free-atom density is spherical (angular averaging of the Fock matrix plus fractional occupation of open shells) and consistent with the molecular functional.

pymbd.pyscf.hirshfeld_volume_ratios(mf)

Per-atom Hirshfeld volume ratios from a converged PySCF Kohn–Sham object.

Parameters:

mf – converged PySCF KS object (restricted or unrestricted)

Returns one ratio per atom carrying a nucleus, in order. Ghost centers have no free-atom reference and are skipped, so adding counterpoise ghosts leaves the ratios (and hence the MBD energy) almost unchanged – they then differ only through the improvement the extra basis functions make to the molecular density, not through the partitioning itself.

Requires an all-electron KS-DFT mean-field; Hartree–Fock (no exchange-correlation functional) and ECPs are rejected, as is a per-element basis, since the free-atom reference is built in the molecular basis. The functional, basis and integration grid are all taken from mf.

The ratios are much less sensitive to the integration grid than the energy is, because the same grid integrates both the in-molecule and the free-atom volume and the quadrature error largely cancels in their ratio. The cancellation is not unconditional, though, and how far it reaches depends on the elements present: hydrogen carries the coarsest atomic grid, so hydrogen-rich systems need a finer one than the cancellation alone would suggest. Choosing the grid for the density is enough in practice; PySCF’s coarsest level is not.

The free-atom reference is deliberately computed in the molecular basis. Any other choice mixes basis-set incompleteness into the ratio, so that a well-separated dimer no longer comes out at exactly 1 and the error carries straight into \(C_6\). A fixed reference basis is correct only at the basis-set limit.

pymbd.pyscf.mbd_energy(mf, beta=None, **kwargs)

MBD dispersion energy (a.u.) for a converged PySCF Kohn–Sham object.

Parameters:

• mf – converged PySCF KS object

• beta (float) – MBD range-separation parameter; if None it is looked up from the molecular functional (currently PBE, PBE0, HSE)

• kwargs – passed on to mbd_energy(), whose defaults (MBD@rsSCS with Fermi dipole damping) are used as they stand

Free-atom vdW parameters come from libMBD’s built-in table, rescaled by the Hirshfeld volume ratios with from_volumes().

Fortran bindings

class pymbd.fortran.MBDGeom(coords, lattice=None, k_grid=None, custom_k_pts=None, n_freq=None, do_rpa=False, get_spectrum=False, get_rpa_orders=False, rpa_rescale_eigs=False, max_atoms_per_block=None, ewald_cutoff_scaling=(1.0, 1.0))

Represents an initialized libMBD geom_t object.

Parameters:

• coords (array-like) – (a.u.) atomic coordinates as rows

• lattice (array-like) – (a.u.) lattice vectors as rows

• k_grid (array-like) – number of \(k\)-points per reciprocal lattice vector

• custom_k_pts (array-like) – (a.u.) custom \(k\)-points as rows

• n_freq (int) – number of quadrature points for frequency integration

• do_rpa (bool) – whether to calculate MBD energy via frequency integration

• get_spectrum (bool) – whether to return eigenvalues and eigenvectors

• get_rpa_orders (bool) – whether to return RPA order decomposition

• rpa_rescale_eigs (bool) – whether to rescale RPA eigenvalues

property coords

(a.u.) Atom coordinates in rows.

has_lattice()

Whether structure is a crystal.

int_density(pts, q, m, w_t, modes)

Evaluate the density of the interacting QDOs.

Returns the charge density of the quantum Drude oscillators (QDOs) with the dipole coupling switched on. Unlike nonint_density(), each per-atom contribution is a general (anisotropic) Gaussian obtained by integrating out the remaining oscillators from the coupled ground-state wavefunction (see QDO charge densities for the derivation and a plotting example). Like the non-interacting density, it integrates to \(\sum_A q_A\).

The coupled frequencies w_t and eigenmodes modes are the MBD normal modes, obtained from mbd_energy() with a geometry constructed with get_spectrum=True (w_t are the square roots of the returned eigenvalues).

Parameters:

• pts (array-like) – (a.u.) points at which to evaluate the density, shape (n_pts, 3)

• q (array-like) – (a.u.) oscillator charges

• m (array-like) – (a.u.) oscillator masses

• w_t (array-like) – (a.u.) coupled MBD frequencies \(\tilde\omega\), shape (3 n_atoms,)

• modes (array-like) – MBD eigenmodes, shape (3 n_atoms, 3 n_atoms)

Returns:

(a.u.) density evaluated at pts, shape (n_pts,)

property lattice

(a.u.) Lattice vectors in rows.

mbd_energy(alpha_0, C6, R_vdw=None, beta=0.0, a=6.0, sigma=None, damping='fermi,dip', variant='rsscs', force=False, vdw_params_grad=False, intermediates=False)

Calculate an MBD energy.

Parameters:

• alpha_0 (array-like) – (a.u.) atomic polarizabilities

• C6 (array-like) – (a.u.) atomic \(C_6\) coefficients

• R_vdw (array-like) – (a.u.) atomic vdW radii

• sigma (array-like) – (a.u.) oscillator widths

• beta (float) – MBD damping parameter \(\beta\)

• a (float) – MBD damping parameter \(a\)

• str (variant) – type of damping

• str – one of ‘plain’, ‘scs’, ‘rsscs’

• bool (vdw_params_grad) – if True, calculate energy gradients

• bool – if True, calculate the energy gradients with respect to alpha_0, C6 and R_vdw

mbd_energy_species(species, volume_ratios, beta, **kwargs)

Calculate an MBD energy from atom types and Hirshfed-volume ratios.

Parameters:

• species (array-like) – atom types (elements)

• volume_ratios (array-like) – ratios of Hirshfeld volumes in molecule and vacuum

• beta (float) – MBD damping parameter \(\beta\)

• kwargs – see mbd_energy()

nonint_density(pts, q, m, w)

Evaluate the density of the non-interacting QDOs.

Returns the charge density of the quantum Drude oscillators (QDOs) in the absence of the dipole coupling, that is a sum of independent spherical Gaussians, one per atom,

\[n(\mathbf r)=\sum_A q_A\left(\frac{m_A\omega_A}\pi\right)^\frac32 \exp\big(-m_A\omega_A|\mathbf r-\mathbf R_A|^2\big)\]

which integrates to \(\sum_A q_A\). The oscillator parameters are related to the vdW parameters by \(\omega_A=\tfrac43 C_{6,A}/\alpha_{0,A}^2\) and \(\alpha_{0,A}=q_A^2/(m_A\omega_A^2)\). The MBD convention is \(m_A=1\), which gives \(q_A=\omega_A\sqrt{\alpha_{0,A}}\). See QDO charge densities for a plotting example.

Parameters:

• pts (array-like) – (a.u.) points at which to evaluate the density, shape (n_pts, 3)

• q (array-like) – (a.u.) oscillator charges

• m (array-like) – (a.u.) oscillator masses

• w (array-like) – (a.u.) oscillator frequencies \(\omega\)

Returns:

(a.u.) density evaluated at pts, shape (n_pts,)

print_timing()

Print timing from libMBD.

ts_energy(alpha_0, C6, R_vdw, sR, d=20.0, damping='fermi', force=False, vdw_params_grad=False)

Calculate a TS energy.

Parameters:

• alpha_0 (array-like) – (a.u.) atomic polarizabilities

• C6 (array-like) – (a.u.) atomic \(C_6\) coefficients

• R_vdw (array-like) – (a.u.) atomic vdW radii

• sR (float) – TS damping parameter \(s_R\)

• d (float) – TS damping parameter \(d\)

• str (damping) – type of damping

• bool (vdw_params_grad) – if True, calculate energy gradients

• bool – if True, calculate the energy gradients with respect to alpha_0, C6 and R_vdw

ts_energy_species(species, volume_ratios, beta, **kwargs)

Calculate a TS energy from atom types and Hirshfed-volume ratios.

Parameters:

• species (array-like) – atom types (elements)

• volume_ratios (array-like) – ratios of Hirshfeld volumes in molecule and vacuum

• sR (float) – TS damping parameter \(s_R\)

• kwargs – see ts_energy()

pymbd.fortran.with_mpi

Whether libMBD was compiled with MPI

pymbd.fortran.with_scalapack

Whether libMBD was compiled with Scalapack