DFT issues

[Shahxad-Akram]

Hi, I'm trying to run a DFT calculation to get a Potential at different distances but around 15-20A Pyscf is showing very weird behavior. I ran the same calculation with Turbomole and it's perfectly fine.

This was Benzene and He interaction and then I tried He He interaction and it was the same, particularly for the larger basis sets e.g. def2-QZVPPD

Here's my code:

import numpy as np
from src.constants import *

from gpu4pyscf import dft
from pyscf import gto

res = []

def dft_calc():

    import dftd4.pyscf as disp

    results = {}
    for functional in ["SCAN"]:
        for basis_set in ["def2-TZVPPD"]:
            point_en = []
            for distance in pes_points()[3:]:
                benzhe = gto.Mole(
                    atom=f"""He     0.000000    0.000000    0.000000
                He     0.000000    0.000000    {distance:.5f}""",
                    basis=basis_set,
                    symmetry=True,
                    verbose=1,
                )

                dft_ins = dft.rks.RKS(benzhe)
                dft_ins.xc = functional
                dft_ins.grids.level = 5
                dft_ins.kernel()

                d4 = disp.DFTD4Dispersion(benzhe, xc=functional)
                correc = d4.kernel()[0]
                res.append((dft_ins.e_tot+ correc, distance))

    print(res)

dft_calc()

[wxj6000]

I tried to use the same setting as you did, however the behavior was not reproduced. Can you please provide a runnable example? The case of He dimer should be enough.

[sunqm]

I cannot reproduce the issue. I guess it is not a problem of gpu4pyscf. It's probably just a common issue for PES jobs.

There are two tricks may be helpful for PES scanning. You can use the .as_scanner() method to generate a scanner object, then feed the geometry into the object for the energy. The scanner object will use the initial guess from the previous point. It helps to produce smooth curve.

dft_scanner = dft.RKS(mol, xc=xxx).as_scanner()
for distance in [...]:
    mol = gto.M(atom=f'...{distance}')
    e_dft = dft_scanner(mol)

Another trick is to scan in the reversed direction. This sometimes helps the PES stay on the adiabatic state, than jumping to another state at conical intersections.

[Shahxad-Akram]

I tried to use the same setting as you did, however the behavior was not reproduced. Can you please provide a runnable example? The case of He dimer should be enough.

I'm sorry, please run this code:

import numpy as np

from gpu4pyscf import dft
from pyscf import gto
import matplotlib.pyplot as plt

res = []

def dft_calc():

    for functional in ["b3lyp"]:
        for basis_set in ["def2-svp"]:
            point_en = []
            for distance in np.linspace(4.5, 20, 50):
                benzhe = gto.Mole(
                    atom=f"""He     0.000000    0.000000    0.000000
                He     0.000000    0.000000    {distance:.5f}""",
                    basis=basis_set,
                    symmetry=True,
                    verbose=1,
                )

                dft_ins = dft.rks.RKS(benzhe)
                dft_ins.xc = functional
                dft_ins.grids.level = 5
                dft_ins.kernel()

                res.append((dft_ins.e_tot, distance))

dft_calc()

res = np.array(res)

plt.plot(res[:, 1], res[:, 0])

It's returning the following,

[wxj6000]

I can confirm the issue now. The same issue can be reproduced in PySCF on CPU too. The fluctuation (around 1e-6) can not be reduced by tightening the criteria, removing pruning, or adding grids. Do you have any insights? @sunqm

[sunqm]

The fluctuation is caused by the discontinuity in the weights of DFT grids wrt geometry. If you use a very small step size in scan, it probably will produce a curve with a more smooth but periodic structure. If you use a very dense radial grids, such as 300 or 500, the curve can be more smooth. And ultimately, it will be a flatten line.

[wxj6000]

Thanks @sunqm. Those fluctuation is indeed decreasing when I use denser radial grids. The decreasing trend is just slow when we use the regular DFT grids. Here is the curve when (1000,590) grids is used. @shah-xad