Quantum Chemistry Program Interfaces based on integrals

Some quantum chemistry programs allow calculation of the SCF or a correlated wavefunction directly from integrals, avoiding the basis-set specific creation of integrals typically required. Often this involves reading in the integrals from an FCIDUMP file.

Running quantum chemistry from integrals directly would allow some freedom in performing quantum chemistry with integrals from different sources.

ASH has some basic support for facilitating this. Currently supported programs are:

• MRCC : Arbitrary order CC, explicit correlation CC (F12), local natural orbital CC etc.

• pyscf : open-source general QM program: SCF, CC, MCSCF etc.

• ccpy : Specialized CC methods

• Block2: DMRG

• Dice: semi-stochastic heat-bath CI

Read/write AO integrals

Here we will assume that we have the 1-electron and 2-electron integrals available as Numpy arrays.

• The 1-electron and overlap integrals should be 2-dimensional Numpy arrays with dimensions (numbf,numbf) where numbf is number of basis functions.

• The 2-electron integrals can be stored with different permutation symmetry: 8-fold, 4-fold or no symmetry.

  • For 8-fold symmetry we will have a dimensional array of (numbf*(numbf+1)/2)*(numbf*(numbf+1)/2+1)/2 elements.

  • For 4-fold symmetry we will have a 2-dimensional array of (numbf*(numbf+1)/2, numbf*(numbf+1)/2) elements.

  • For no-symmetry we will have a 4-dimensional array of (numbf,numbf,numbf,numbf) elements.

In the examples below we will usually have the 2-electron integrals be a non-symmetric 4-dimensional array with dimensions (numbf,numbf,numbf,numbf).

We can read and write them like this:

from ash import *

numbf = 4

#Define some dummy integrals (here random)
AO_CO_1el = np.random.random((numbf,numbf))
AO_CO_2el = np.random.random((numbf,numbf,numbf,numbf))

# Write the integrals to Numpy binary format
one_el_integrals = np.save("AO_CO_1el.npy", AO_CO_1el)
two_el_integrals = np.save("AO_CO_2el.npy", AO_CO_2el)

# Read the integrals from Numpy binary format
one_el_integrals = np.load("AO_CO_1el.npy")
two_el_integrals = np.load("AO_CO_2el.npy")

Running SCF from integrals, get mf object (pySCF) and MO integrals

In order to run a correlated WF calculation we need integrals in the MO-basis which means that we need to first solve the MO problem, i.e. the SCF. Here we utilize the open-source pySCF program for this purpose. It is possible to create a pySCF mean-field object directly from AO integrals in pySCF, without ever defining the standard GTO basis set. The ASH function create_pyscf_mol_and_mf is convenient for this purpose.

To do this, we need to specify how many electrons we have, the nuclear repulsion energy of the system and read in the integrals as Numpy arrays. Here we read in the integrals from Numpy binary files.

from ash import *
from ash.interfaces.interface_pyscf import create_pyscf_mol_and_mf

#Defining the basic parameters of the system
nuc_repulsion_energy = 22.51817918808511 # Nuclear repulsion energy in Eh
# Or use function: nuc_nuc_repulsion(coords, charges)  (coordinates in Angstrom)
num_el=14

scf_type="RHF" # RHF, UHF, ROHF
mult = 1 # Spin multiplicity
#num_corr_el = 2 # Number of active electrons (i.e. non-frozen) in post-HF.
#num_corr_orbs = 22 # Number of active orbitals in post-HF
#numocc = num_corr_el/2 #Number of occupied orbitals

#AO integrals read as Numpy arrays from disk (Numpy binary format)
one_el_integrals = np.load("AO_CO_1el.npy")
two_el_integrals = np.load("AO_CO_2el.npy")
overlap = np.load("AO_CO_overlap.npy")

num_basis=one_el_integrals.shape[0] #Number of basis functions

#Create mol and mf objects from integrals
mol, mf = create_pyscf_mol_and_mf(numel=num_el, mult=mult, nuc_repulsion_energy=nuc_repulsion_energy,
    one_el_integrals=one_el_integrals, two_el_integrals=two_el_integrals, overlap=overlap )

#Run the mf object (i.e. solve the SCF) to get the MOs
mf.kernel()
# MO coefficients : mf.mo_coeff
# MO occupations: mf.mo_occ
# MO energies: mf.mo_energy

Once the meanfield object has been run and we have the MOs, we can perform the AO->MO integral transformation. This can also be accomplished by tools inside pySCF, namely the ao2mo module, see https://pyscf.org/contributor/ao2mo_developer.html https://pyscf.org/pyscf_api_docs/pyscf.ao2mo.html

from pyscf import ao2mo
# 2-el and 1-el integral transformation from AO to MO
twoel_MObas = ao2mo.kernel(two_el_integrals, mf.mo_coeff)
oneel_MObas = np.einsum("ap, ab, bq -> pq", mf.mo_coeff, one_el_integrals, mf.mo_coeff)

Write MO integrals to disk

Now that we have the integrals in the MO-basis it is convenient to write them to disk. A common standard is the FCIDUMP format defined in the article: Knowles, P. J.; Handy, N. C. A Determinant Based Full Configuration Interaction Program. Computer Physics Communications 1989, 54, 75–83. https://doi.org/10.1016/0010-4655(89)90033-7. ↩

ASH features the function ASH_write_integralfile for this purpose.

def ASH_write_integralfile(two_el_integrals=None, one_el_integrals=None, nuc_repulsion_energy=None, header_format="MRCC",
                            num_corr_el=None, filename=None, int_threshold=1e-16, scf_type="RHF", mult=None,
                            symmetry_option=0, orbsym=None):

For writing a standard FCIDUMP file we would do:

from ash import *
from ash.interfaces.interface_pyscf import create_pyscf_mol_and_mf
from pyscf import ao2mo

nuc_repulsion_energy = 22.51817918808511 # Nuclear repulsion energy in Eh
num_el=14
scf_type="RHF"
mult = 1

#AO integrals read as Numpy arrays from disk (Numpy binary format)
one_el_integrals = np.load("AO_CO_1el.npy")
two_el_integrals = np.load("AO_CO_2el.npy")
overlap = np.load("AO_CO_overlap.npy")

# Run SCF via pySCF
mol, mf = create_pyscf_mol_and_mf(numel=num_el, mult=mult, nuc_repulsion_energy=nuc_repulsion_energy,
        one_el_integrals=one_el_integrals, two_el_integrals=two_el_integrals, overlap=overlap )
mf.kernel()

# AO->MO integral transformation
twoel_MObas = ao2mo.kernel(two_el_integrals, mf.mo_coeff)
oneel_MObas = np.einsum("ap, ab, bq -> pq", mf.mo_coeff, one_el_integrals, mf.mo_coeff)

# Write FCIDUMP file
ASH_write_integralfile(two_el_integrals=twoel_MObas, one_el_integrals=oneel_MObas,
    nuc_repulsion_energy=nuc_repulsion_energy, header_format="FCIDUMP",
    num_corr_el=num_el, int_threshold=1e-16, scf_type=scf_type, mult=mult)

For MRCC calculations (see next) we want the FCIDUMP file to have an MRCC-specific header and have the filename be fort.55 and so we would run the function like this:

nuc_repulsion_energy = 22.51817918808511 # Nuclear repulsion energy in Eh
num_el=14
scf_type="RHF"
mult = 1
ASH_write_integralfile(two_el_integrals=twoel_MObas, one_el_integrals=oneel_MObas,
    nuc_repulsion_energy=nuc_repulsion_energy, header_format="MRCC",
    num_corr_el=num_corr_el, filename="fort.55", int_threshold=1e-16, scf_type=scf_type, mult=mult)

Running CC from integrals (MRCC)

To run a MRCC calculation directly from MO integrals we need 2 files: the integral-file called fort.55 and a special basic inputfile named fort.56. We can write the inputfile like below where we have specified the excitation level to be 4 (corresponding to CCSDTQ), scf_type to be RHF etc. To request the calculation of reduced density matrices we specify dens = 2. We also need to specify the occupations to use for the WF calculation and define how many electrons should be correlated etc.

num_corr_el = 2 # Number of active electrons (i.e. non-frozen) in post-HF.
numocc = num_corr_el/2 #Number of occupied orbitals
num_basis=one_el_integrals.shape[0] #Number of basis functions
occupations = [2.0 if c < numocc else 0.0 for i,c in enumerate(range(num_basis))]

MRCC_write_basic_inputfile(occupations=occupations, filename="fort.56", scf_type="RHF",
                           ex_level=4, nsing=1, ntrip=0, rest=0, CC_CI=1, dens=2, CS=1,
                           spatial=1, HF=1, ndoub=0, nacto=0, nactv=0, tol=9, maxex=0,
                           sacc=0, freq=0.0000, symm=0, conver=0, diag=0, dboc=0, mem=1024)

The inputfile is written to disk as fort.56 and looks like this:

2    1    0     0    1    -2     0     0     0    1    1     1      0    0      0    9      0     0 0.0     0 1024
ex.lev, nsing, ntrip, rest, CC/CI, dens, conver, symm, diag, CS, spatial, HF, ndoub, nacto, nactv, tol, maxex, sacc, freq, dboc, mem
2.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0

The file can be modified as desired. Once both fort.56 and fort.55 have been created we can run MRCC directly.

from ash import *

numcores = 1
run_mrcc("/path/to/mrccdir","mrcc.out", "OMP", numcores)

The density matrices are then available in the file CCDENSITIES

Running CC calculations from integrals (pySCF)

pySCF has support for various correlated wavefunctions, including the ability to get RDM1 and RDM2 from CCSD and CCSD(T) wavefunctions. If we already have a pySCF mean-field object created, we can use the ASH wrapper around pySCF to create a PySCFTheory object and request a CC calculation.

from ash import *

#Fragment
coordsstring="""C 0 0 0
O 0 0 1.128
"""
frag = Fragment(coordsstring=coordsstring, charge=0, mult=1)


# Create ASH PySCFTheory object via previously created mf object (see earlier)
pyscfobject = PySCFTheory(scf_type="RHF", mf_object=mf, CC=True, CCmethod="CCSD(T)", do_pop_analysis=False, symmetry=None)

Singlepoint(theory=pyscfobject, fragment=frag)

Once the pySCF CC calculation has been done we can request calculation of the RDM1 and RDM2 in AO or MO basis and save them to disk.

# RDM1 in AO basis
rdm1_MO = pyscfobject.ccobject.make_rdm1(ao_repr=False)
# Convert RDM1 from MO to AO basis
new_rdm_AO = DM_MO_to_AO(rdm1_MO, mf.mo_coeff)

# RDM2 in AO and MO basis
rdm2_MO = theory.ccobject.make_rdm2(ao_repr=False)
rdm2_AO = theory.ccobject.make_rdm2(ao_repr=True)

# Write the RDMs to disk as Numpy binary files or text
np.save("rdm1_MObasis", rdm1_MO)
np.save("rdm1_newAObasis", new_rdm_AO)
np.savetxt("rdm1_MObasis.txt", rdm1_MO)
np.savetxt("rdm1_newAObasis.txt", new_rdm_AO)

np.save("rdm2_MObasis", rdm2_MO)
np.savetxt("rdm2_MObasis.txt", rdm2_MO)
np.save("rdm2_AObasis", rdm2_AO)
np.savetxt("rdm2_AObasis.txt", rdm2_AO)

Running DMRG from integrals (Block2)

A DMRG calculation with the Block2 interface can be run from integrals in a few different ways:

1. Via PySCFTheory object. If the pySCF mean-field object has been created (see earlier) we can create an ASH pySCFTheory object which is conveniently passed over to BlockTheory.

from ash import *

#Fragment
coordsstring="""C 0 0 0
O 0 0 1.128
"""
frag = Fragment(coordsstring=coordsstring, charge=0, mult=1)


# Create ASH PySCFTheory object via previously created mf object (see earlier)
pyscfobject = PySCFTheory(scf_type="RHF", mf_object=mf, do_pop_analysis=False, symmetry=None)

#Block2 DMRG calculation via input pyscftheory object. Here we define a DMRG-CASCI(2,10) calculation from the input SCF orbitals
blocktheory = BlockTheory(pyscftheoryobject=pyscfobject, active_space=[2,10],  macroiter=0,
    numcores=1, memory=50000, tol=1e-8, initial_orbitals='HF', block_parallelization='OpenMP',
    maxM=1000, singlet_embedding=True, DMRG_DoRDM=True, DMRG_DoRDM2=True)

# Run
Singlepoint(theory=blocktheory, fragment=frag)

2. Via MO-basis integrals in FCIDUMP format.

If we have the 1- and 2-electron integrals available in the MO-basis we can write an FCIDUMP file to disk and start a DMRG calculation directly from the FCIDUMP file. We first need to get the integrals in the MO-basis. Here we first run an SCF via pySCF as before and then convert the integrals from AO to MO basis.

from ash import *


#Fragment
coordsstring="""C 0 0 0
O 0 0 1.128
"""
frag = Fragment(coordsstring=coordsstring, charge=0, mult=1)

# Get 1- and 2-integrals in AO basis
one_el_integrals = np.load("AO_CO_1el.npy")
two_el_integrals = np.load("AO_CO_2el.npy")
overlap = np.load("AO_CO_overlap.npy")

#Create mol and mf objects from integrals
mol, mf = create_pyscf_mol_and_mf(numel=14, mult=1,
    nuc_repulsion_energy=22.51817918808511,
    one_el_integrals=one_el_integrals, two_el_integrals=two_el_integrals,
    overlap=overlap )


# AO->MO basis transformation via pyscf
#2-el and 1-el integral transformation from MO to AO
from pyscf import ao2mo
twoel_MObas = ao2mo.kernel(two_el_integrals, mf.mo_coeff)
oneel_MObas = np.einsum("ap, ab, bq -> pq", mf.mo_coeff, one_el_integrals, mf.mo_coeff)

###########################################################
# Write MO-basis integrals to disk as FCIDUMP-style file
###########################################################
ASH_write_integralfile(two_el_integrals=twoel_MObas, one_el_integrals=oneel_MObas,
    nuc_repulsion_energy=22.51817918808511, header_format="FCIDUMP",
    num_corr_el=14, filename="system.fcidump", int_threshold=1e-16,
    scf_type='RHF', mult=1)


#Block2 DMRG calculation via input FCIDUMP file. Here we define a DMRG-CASCI(2,10) calculation from the input SCF orbitals
blocktheory = BlockTheory(fcidumpfile="system.fcidump", active_space=[2,10],  macroiter=0,
    numcores=1, memory=50000, tol=1e-8, initial_orbitals='HF', block_parallelization='OpenMP',
    maxM=1000, singlet_embedding=True, DMRG_DoRDM=True, DMRG_DoRDM2=True)

# Run
Singlepoint(theory=blocktheory, fragment=frag)

Once a DMRG calculation has been run with RDMs requested, the RDMs are accessible from the BlockTheory object.

#Grab RDMs in MO basis and do RDM1-MO to AO conversion
rdm1_MO = blocktheory.properties["rdm1_MO"]
print("RDM1-MO:", rdm1_MO)
from ash.functions.functions_elstructure import DM_MO_to_AO
# MO coefficients come from mf object
rdm1_AO = DM_MO_to_AO(rdm1_MO, mf.mo_coeff)
print("rdm1_AO:", rdm1_AO)

Running SHCI from integrals (Dice)

A semistochastic heatbath CI calculation with the Dice interface can be run from integrals almost identically to the DMRG-Block route above. Simple follow the examples above but define a DiceTheory like this:

#From pyscf object
dicetheory = DiceTheory(pyscftheoryobject=pyscftheoryobject, SHCI_active_space=[2,10],
    numcores=1, memory=50000, initial_orbitals='HF',
    SHCI_DoRDM=True, SHCI_DoRDM2=True)

# Or FCIDUMP-file:
dicetheory = DiceTheory(fcidumpfile="system.fcidump", SHCI_active_space=[2,10],
    numcores=1, memory=50000, initial_orbitals='HF',
    SHCI_DoRDM=True, SHCI_DoRDM2=True)

Once a SHCI calculation has been run with RDMs requested, the RDMs are accessible from the BlockTheory object.

#Grab RDMs in MO basis and do RDM1-MO to AO conversion
rdm1_MO = dicetheory.properties["rdm1_MO"]
print("RDM1-MO:", rdm1_MO)
from ash.functions.functions_elstructure import DM_MO_to_AO
# MO coefficients come from mf object
rdm1_AO = DM_MO_to_AO(rdm1_MO, mf.mo_coeff)
print("rdm1_AO:", rdm1_AO)

Running CC from integrals (ccpy)

ccpy is a coupled cluster program in Python (with important routines in Fortran) that includes a number of interesting specialized CC methods. A coupled cluster calculation with the interface to the ccpy program can be run from integrals like this:

From a PySCFTheory object:

First create a mf object as shown before, then create a PySCFTheory object.

#Fragment
coordsstring="""C 0 0 0
O 0 0 1.128
"""
frag = Fragment(coordsstring=coordsstring, charge=0, mult=1)

# Create ASH PySCFTheory object via previously created mf object (see earlier)
pyscfobject = PySCFTheory(scf_type="RHF", mf_object=mf, do_pop_analysis=False, symmetry=None)

# A CCSD ccpyTheory object using pyscfobj as input
ccpy_theory = ccpyTheory(method="CCSD", pyscftheoryobject=pyscfobj, frozencore=True,
            cc_tol=1e-10, numcores=1, cc_maxiter=300)

result = Singlepoint(theory=ccpy_theory, fragment=frag)

From an FCIDUMP file:

If an FCIDUMP file in MO-basis is already available (see earlier) we can start a ccpy calculation directly from it.

#Fragment
coordsstring="""C 0 0 0
O 0 0 1.128
"""
frag = Fragment(coordsstring=coordsstring, charge=0, mult=1)

# A CCSD ccpyTheory object using FCIDUMP-file (MO-baiss) as input
theory = ccpyTheory(method="CCSD", fcidumpfile="FCIDUMP-file", frozencore=True,
            cc_tol=1e-10, numcores=1, cc_maxiter=300)

result = Singlepoint(theory=ccpy_theory, fragment=frag)