how to use keywords with pyscf

[small-bamboo]

I found in the doc this is no example for a keyword with pyscf, I wang to use the pbe0(D3)/x2c-TZVPall with sfx2c1e methond in pyscf to optimize structure, how can I create a yaml file?

[eljost]

Dear small-bamboo,

can you please provide me an example how to do this type of calculation in PySCF via a python script? Then I can see, if this can be translated to an pysisyphus-input, or how I should modify the code, to make this calculation possible.

Personally I don't really use PySCF, beside in the test suite of pysisyphus.

All the best
Johannes

[small-bamboo]

Dear eljost:
here is a simple example

from pyscf import gto, scf, cc,lib,dft,dftd3
import basis_set_exchange
from pyscf.lib import logger
import numpy as np
from pyscf.qsdopt.qsd_optimizer import QSD

mol = gto.M(
atom = '''

C -2.01237621 0.87642747 0.00764201
H -2.03228125 0.27575325 -0.91175712
H -2.93571568 1.45458485 0.00419744
H -2.05515335 0.25181909 0.91035923
H -0.60208942 1.59377024 0.00655957

''',

basis = {
'C': gto.load(basis_set_exchange.api.get_basis('x2c-TZVPall', elements='C', fmt='nwchem'), 'C'),
'H': gto.load(basis_set_exchange.api.get_basis('x2c-TZVPall', elements='H', fmt='nwchem'), 'H'),
},
charge = 0,
spin = 0,
)

mol.verbose = 4

#mf = scf.HF(mol).sfx2c1e().newton()

mf= dft.UKS(mol)
mf.xc = 'pbe0'
mf = mf.sfx2c1e()
mf.conv_tol = 1e-6
mf.max_cycle = 200
mf.kernel()
mf = stable_opt_internal(mf)
optimizer = QSD(mf, stationary_point="TS")
optimizer.kernel()

[eljost]

Dear small-bamboo,

I implemented some changes on the wigner branch, that should make your calculation possible, but not via YAML input. Only via python.

You have to subclass the PySCF-calculator and override prepare_mol and prepare_mf. See below:

import basis_set_exchange

# import dftd3.pyscf as d3
from pyscf import gto

from pysisyphus.calculators.PySCF import PySCF
from pysisyphus.helpers import geom_loader


class ModPySCF(PySCF):
    def prepare_mol(self, atoms, coords, build=True):
        # Don't build yet, as we want to modify the basis.
        mol = super().prepare_mol(atoms, coords, build=False)
        mol.basis = {
            "C": gto.load(
                basis_set_exchange.api.get_basis(
                    "x2c-TZVPall", elements="C", fmt="nwchem"
                ),
                "C",
            ),
            "H": gto.load(
                basis_set_exchange.api.get_basis(
                    "x2c-TZVPall", elements="H", fmt="nwchem"
                ),
                "H",
            ),
        }
        # Build after customization.
        mol.build(parse_arg=False)
        return mol

    def prepare_mf(self, mf):
        mf = mf.sfx2c1e()
        # mf = d3.energy(mf)  # This does not seem to work with forces/gradients
        return mf


xyz_str = """5

C -2.01237621 0.87642747 0.00764201
H -2.03228125 0.27575325 -0.91175712
H -2.93571568 1.45458485 0.00419744
H -2.05515335 0.25181909 0.91035923
H -0.60208942 1.59377024 0.00655957"""
geom = geom_loader(xyz_str)
calc = ModPySCF(xc="pbe0", basis="dummy", unrestricted=True, verbose=4)
mol = calc.prepare_mol(geom.atoms, geom.cart_coords)
geom.set_calculator(calc)
print("energy", geom.energy)
print("forces", geom.forces)

DFTD3-Gradients don't seem to work though...

Let me know if this helps, or even if it doesn't ;) It is unlikely that I'll provide a simpler way to do this, as I don't really use the PySCF calculator by myself.

All the best
Johannes