PBC: Nuclear gradient of total energy with respect to nuclear coordiantes

[computative]

I want to calculate forces aka gradients of total energy with respect to nuclear coordinates in periodic boundary conditions. Can I do this with pyscf?

It should be possible to do this (feynman-hellman forces) just given the density merely by applying feynman-hellmann to the hamiltonian (see Mike Finnis book on interatomic forces in condensed matter page 79 eqn 3.2)

1. feynman hellmann forces are scf method independent as you can see. Can I get this integral with pyscf, perhaps using the cell.intor-method somehow?
2. Can I get the gradient of total energy with respect to nuclear coordinates for rks in pbc any other way?

[susilehtola]

The Hellmann-Feynman theorem is not satisfied in atomic-orbital basis sets due to the geometry dependence of the basis set; you also get so-called Pulay terms.

It doesn't seem like nuclear gradients are very well documented in the examples; however, the code in examples/pbc/27-multigrid.py does contain the following code snippet

# Nuclear Gradients
from pyscf.pbc.grad import rks as rks_grad
grad = rks_grad.Gradients(mf)
g = grad.kernel()

This driver should also calculate the Pulay terms and yield the correct gradient.

IIRC geometry optimization in periodic calculations also requires the calculation of the stress tensor, since the cell parameters also need to be able to change. I did not find anything on this in the examples.

[computative]

Thanks. That module is broken.

import numpy 
import pyscf, copy
from pyscf.pbc import gto, dft
from pyscf import data
import numpy as np
from pyscf.pbc.dft import numint

cell = gto.Cell()

R = np.random.rand(2,3)
A = np.eye(3)

cell.a = "\n".join([f"{ str(ai).strip('[]') }" for ai in A])
cell.atom = "\n".join([ f"C  { str(r).strip('[]') }" for r in R])
cell.verbose = 1
cell.basis = 'gth-szv'
cell.build()

mf = dft.RKS(cell)
mf.conv_tol = 1e-5
mf.kernel()

from pyscf.pbc.grad import rks as rks_grad
grad = rks_grad.Gradients(mf)
g = grad.kernel()

"""
$ python grad.py 
Traceback (most recent call last):
  File "~/Desktop/grad.py", line 28, in <module>
    g = grad.kernel()
        ^^^^^^^^^^^^^
  File "~/Programs/anaconda3/lib/python3.11/site-packages/pyscf/grad/rhf.py", line 415, in kernel
    de = self.grad_elec(mo_energy, mo_coeff, mo_occ, atmlst)
         ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "~/Programs/anaconda3/lib/python3.11/site-packages/pyscf/pbc/grad/rhf.py", line 46, in grad_elec
    vhf = mf_grad.get_veff(mol, dm0, kpt)
          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "~/Programs/anaconda3/lib/python3.11/site-packages/pyscf/pbc/grad/rhf.py", line 164, in get_veff
    return get_veff(self, mol, dm, kpt)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "~/Programs/anaconda3/lib/python3.11/site-packages/pyscf/pbc/grad/rhf.py", line 131, in get_veff
    return -mydf.get_veff_ip1(dm, xc_code=xc_code, kpts=kpts)
            ^^^^^^^^^^^^^^^^^
AttributeError: 'FFTDF' object has no attribute 'get_veff_ip1'
"""

But am able to do it with the MultiGridFFTDF2 as in the 27-multigrid2.py. Are the developers aware of this issue?