PySCF calculator using CASSCF
EmotionalSupportBurrito opened this issue · 4 comments
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!
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()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!
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()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.