PySCF calculator using CASSCF

[EmotionalSupportBurrito]

Hi! First of all, thank you very much for writing this code. I am trying to use pysisyphus with PySCF and multi-reference calculations based on MP2 natural orbitals. An example code for the MP2/CASSCF calculation would be:

from pyscf import gto,scf,mcscf,mp
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"""
mol      = gto.M()
mol.atom = atom
mol.build()
mf       = scf.RHF(mol)
mf.kernel()
mymp = mp.UMP2(mf)
mymp.kernel()
noons, natorbs = mcscf.addons.make_natural_orbitals(mymp)
ncas, nelecas  = (2,2)
mycas          = mcscf.CASSCF(mf, ncas, nelecas)
mycas.kernel(natorbs)

After looking at the PySCF calculator source code, it seems that right now there is only support for HF, DFT and MP2.

However, I was wondering if one can "hack" the code a little bit as you suggested some time ago here. The code from the link did not run, but I went through the API documentation and changed a few lines, which seems to work.

import pysisyphus.Geometry as Geometry
from pysisyphus.calculators.PySCF import PySCF
from pysisyphus.xyzloader import parse_xyz_str
from pyscf import gto,scf

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"""

class ModPySCF(PySCF):
    def prepare_mol(self, atoms, coords, build=True):
        mol = gto.M()
        mol.atom = atom
        mol.build()
        return mol

    def prepare_mf(self, mf):
        mf = scf.RHF(mol)
        #PUT CASSCF AND MP2 here ?
        return mf

atoms,coords  = parse_xyz_str(xyz_str,False)
geom          = Geometry.Geometry(atoms,coords)
calc          = ModPySCF(basis="STO3G", verbose=4)
geom.set_calculator(calc)
print("energy", geom.energy)
print("forces", geom.forces)

I was wondering if one could modify this code to include the MP2 and CASSCF calculations. PySCF has a as_scanner method for casscf, and I am wondering if one could use that.

Thank you very much!

[eljost]

Maybe something along these lines will do the job

import numpy as np

from pyscf import gto, scf, mcscf, mp

from pysisyphus.calculators.Calculator import Calculator
from pysisyphus.helpers import geom_loader

def as_atom_str(atoms, coords, fmt=".8f"):
    per_atom = list()
    for atom, (x, y, z) in zip(atoms, coords.reshape(-1, 3)):
        per_atom.append(f"{atom} {x:{fmt}} {y:{fmt}} {y:{fmt}}")
    mol_str = "; ".join(per_atom)
    return mol_str

def as_mol(atoms, coords):
    atom_str = as_atom_str(atoms, coords)
    return gto.M(atom=atom_str, unit="Bohr")

class PySCFScanner(Calculator):
    def __init__(self, mf, **kwargs):
        super().__init__(**kwargs)

        self.mf = mf
        self.scanner = mf.nuc_grad_method().as_scanner()

    def get_energy(self, atoms, coords):
        return self.get_forces(atoms, coords)

    def get_forces(self, atoms, coords):
        mol = as_mol(atoms, coords)
        energy, grad = self.scanner(mol)
        return {
            "energy": energy,
            "forces": -np.ravel(grad),
        }

def get_mf(atoms, cart_coords):
    mol = gto.M(verbose=4, basis="sto3g", output="scanner.log", unit="Bohr")
    mol.atom = as_atom_str(geom.atoms, geom.cart_coords)
    mol.build(parse_arg=False)
    mf = scf.RHF(mol)
    mf.kernel()
    print("Did RHF")
    mymp = mp.UMP2(mf)
    mymp.kernel()
    print("Did UMP2")
    noons, natorbs = mcscf.addons.make_natural_orbitals(mymp)
    print("Created natural orbitals")
    ncas, nelecas = (2, 2)
    mf = mcscf.CASSCF(mf, ncas, nelecas)
    # Seed with active space guess
    mf.kernel(natorbs)
    print("Initialized active space")
    return mf

geom = geom_loader("lib:h2o.xyz")
mf = get_mf(geom.atoms, geom.cart_coords)
calc = PySCFScanner(mf)
geom.set_calculator(calc)
from pysisyphus.optimizers.RFOptimizer import RFOptimizer

opt = RFOptimizer(geom, dump=True)
opt.run()

[EmotionalSupportBurrito]

Thank you very much for the quick reply! I tried the code and it works. In case somebody wants to use this with a different basis set, I changed the as_mol() function to

as_mol(atoms, coords):
    atom_str = as_atom_str(atoms, coords)
    return gto.M(atom=atom_str, unit="Ang",basis="def2-svp")  # ADDED basis

I wanted to use this calculator for NEB calculations. As a test, I wanted to calculate the isomerization of cyclobutadiene:

import pysisyphus.Geometry as Geometry
from pysisyphus.xyzloader import parse_xyz_str
from pysisyphus.cos.NEB import NEB

cyclobut_start= """
8

C   0.6782562916    0.7728795927    -0.0002401881
C   -0.6784975438   0.7726666140    -0.0002477805
C   0.6784985949    -0.7726661155   0.0002454522
C   -0.6782544577   -0.7728791470   0.0002395929
H   1.4429418388    1.5346909038    -0.0004996409
H   -1.4434223475   1.5342377028    -0.0004777765
H   1.4434231895    -1.5342374436   0.0004693093
H   -1.4429397589   -1.5346907390   0.0004935107
 """

cyclobut_end = """
8

C   0.7726716169    0.6784987395    -0.0002142999
C   -0.7728832231   0.6782560013    -0.0002152374
C   0.7728846696    -0.6782555549   0.0002132490
C   -0.7726701762   -0.6784982832   0.0002124153
H   1.5342470945    1.4434201329    -0.0004548814
H   -1.5346989285   1.4429381636    -0.0004561837
H   1.5347003826    -1.4429376985   0.0004543025
H   -1.5342456528   -1.4434196791   0.0004534706 """

atoms_start,coords_start  = parse_xyz_str(cyclobut_start,False)
geom_start                = Geometry.Geometry(atoms_start,coords_start)

atoms_end,coords_end    = parse_xyz_str(cyclobut_end,False)
geom_end                = Geometry.Geometry(atoms_end,coords)

image_list    = [geom_start,geom_end]
neb           = NEB(image_list)

What I am not sure about is how to set a calculator for the NEB class or how to run it. If I understand the documentation correctly, the NEB class will just return forces. Do I have to parse this to a force-minimizer class or do I have to run the ChainOfStates class?

Thank you very much for your help!

[eljost]

Please see the attached example, using XTB.

from pysisyphus.calculators import XTB
from pysisyphus.cos.NEB import NEB
from pysisyphus.helpers import geom_loader
from pysisyphus.interpolate import interpolate_all
from pysisyphus.optimizers.LBFGS import LBFGS

"""
Steps for COS calculations:
    1. Load geometries
    2. Interpolate between them to get a set of images
    3. Set calculators on the images
    4. Create COS object, e.g., NEB
    5. Pass COS object to optimizer
    6. ...
    7. Profit!
"""

cyclobut_start = """
8

C   0.6782562916    0.7728795927    -0.0002401881
C   -0.6784975438   0.7726666140    -0.0002477805
C   0.6784985949    -0.7726661155   0.0002454522
C   -0.6782544577   -0.7728791470   0.0002395929
H   1.4429418388    1.5346909038    -0.0004996409
H   -1.4434223475   1.5342377028    -0.0004777765
H   1.4434231895    -1.5342374436   0.0004693093
H   -1.4429397589   -1.5346907390   0.0004935107
 """

cyclobut_end = """
8

C   0.7726716169    0.6784987395    -0.0002142999
C   -0.7728832231   0.6782560013    -0.0002152374
C   0.7728846696    -0.6782555549   0.0002132490
C   -0.7726701762   -0.6784982832   0.0002124153
H   1.5342470945    1.4434201329    -0.0004548814
H   -1.5346989285   1.4429381636    -0.0004561837
H   1.5347003826    -1.4429376985   0.0004543025
H   -1.5342456528   -1.4434196791   0.0004534706 """

geom_start = geom_loader(cyclobut_start)
geom_end = geom_loader(cyclobut_end)
between = 11  # Will result in between + 2 images

# If you want to interpolate in internal coordinates you have to recreate the images
# with Cartesian coordinates, as NEB does not support internal coordinates.
# images = interpolate_all((geom_start, geom_end), between=between, kind="redund")
# images = [image.copy(coord_type="cartesian") for image in images]

# Or directly interpolate in Cartesian coordinates
images = interpolate_all((geom_start, geom_end), between=between)

def get_calc_getter():
    calc_number = 0

    def get_calc():
        nonlocal calc_number

        # Create your custom PySCF calculator here, instead of XTB
        calc = XTB(calc_number=calc_number)

        calc_number += 1
        return calc

    return get_calc

get_calc = get_calc_getter()
for image in images:
    calc = get_calc()
    image.set_calculator(calc)

cos_kwargs = {
    # "climb": True,
}
cos = NEB(images, **cos_kwargs)
opt_kwargs = {
    "dump": True,
    # Set thresholds here etc.
}
opt = LBFGS(cos, **opt_kwargs)
opt.run()

[eljost]

as_mol(atoms, coords):
    atom_str = as_atom_str(atoms, coords)
    return gto.M(atom=atom_str, unit="Ang",basis="def2-svp")  # ADDED basis

Internally pysisyphus uses atomic units throught, i.e., Bohr for length. So I guess it would be a bad idea to write functions that take arguments in Angstrom.

Nontheless, if you have a working module that uses my example PySCFScanner I would be happy if you share it, as I somehow got only strange results when minimizing molecular geometries with it ...

EDIT The basis should not matter in the as_mol function, as it only creates a kind of dummy Mol to be passed to the scanner. I'm pretty sure the scanner will utilize the level of theory (method/basis set etc.) of the underlying, original Mol.