#!/usr/bin/env python
# Copyright 2019-2021 The PySCF Developers. All Rights Reserved.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#
# Author: Elvira R. Sayfutyarova
#
''' When using results of this code for publications, please cite the following paper:
"Constructing molecular pi-orbital active spaces for multireference calculations of conjugated systems"
E. R. Sayfutyarova and S. Hammes-Schiffer, J. Chem. Theory Comput., 15, 1679 (2019).
'''
import numpy as np
import numpy.linalg as la
from pyscf import gto
[docs]
def mdot(*args):
"""chained matrix product: mdot(A,B,C,..) = A*B*C*...
No attempt is made to optimize the contraction order."""
r = args[0]
for a in args[1:]:
r = np.dot(r,a)
return r
[docs]
def MakeShellsForElement(mol, Element):
"""make a list with MINAO basis set data for a given element in the mol object"""
PyScfShellInfo = gto.basis.load(mol.basis, Element)
AtShells = []
for ShellData in PyScfShellInfo:
# see: gto/basis/minao.py on format of data.
# data comes as a list of:
# [l, (...), (...), (...)]
# where each (...) specifies both exponents and contractions and has the same length.
l, ExponentsAndContractions = ShellData[0], np.asarray(ShellData[1:])
# primitive exponents
# Exp = ExponentsAndContractions[:,0]
# contraction coefficient matrix
CGTO = ExponentsAndContractions[:,1:]
nCGTO = CGTO.shape[1]
nLk = nCGTO*(2*l+1)
AtShells.append((l, nCGTO, nLk))
return AtShells
[docs]
def MakeShells(mol, Elements):
"""collect MINAO basis set data for all elements"""
Shells = []
for iAt in range(mol.natm):
ElementShells = MakeShellsForElement(mol, Elements[iAt])
Shells.append(ElementShells)
return Shells
[docs]
def Atoms_w_Coords(mol):
"""collect info about atoms' positions"""
AtomCoords = mol.atom_coords()
assert (AtomCoords.shape == (mol.natm, 3))
Elements = [mol.atom_pure_symbol(iAt) for iAt in range(mol.natm)]
Elements =np.asarray(Elements)
return Elements,AtomCoords
[docs]
def MakePiOS(mol,mf,PiAtomsList, nPiOcc=None,nPiVirt=None):
np.set_printoptions(precision=4,linewidth=10000,edgeitems=3,suppress=False)
print("================================CONSTRUCTING PI-ORBITALS===========================================")
# PiAtoms list contains the set of main-group atoms involved in the pi-system of the chromophores(?)
# indices are 1-based.
mol2=mol.copy()
mol2.basis = 'MINAO'
mol2.build()
# make a minimal AO basis for our atoms. We load the basis from a library
# just to have access to its shell-composition. Need that to find the indices
# of all atoms, and the AO indices of the px/py/pz functions.
Elements,Coords = Atoms_w_Coords(mol2)
Shells = MakeShells(mol2,Elements)
def AssignTag(iAt, Element, Tag):
assert (Elements[iAt-1] == Element)
Elements[iAt-1] = Element + Tag
# fix type of atom for the donor-acceptor which are exchanging the hydrogens.
# During the process, the hydrogens are moving and the Huckel-Theory
# formal number of contributed pi-electrons may change.
# These settings override the auto-detection of number of pi electrons
# based on atomic connectivity in GetNumPiElec() below.
# AssignTag(50, "N", "1e")
# OrbBasis = mol.basis
C = mf.mo_coeff
S1 = mol.intor_symmetric("int1e_ovlp")
S2 = mol2.intor_symmetric("int1e_ovlp")
S12 = gto.intor_cross('int1e_ovlp', mol, mol2)
SMo = mdot(C.T, S1, C)
print(" MO deviation from orthogonality {:8.2e} \n".format(rmsd(SMo - np.eye(SMo.shape[0]))))
# make arrays of occupation numbers and orbital eigenvalues.
Occ = mf.mo_occ
Eps = mf.mo_energy
nOrb = C.shape[1]
if (mol.spin==0):
nOcc = np.sum(Occ == 2)
else:
nOcc = np.sum(Occ == 2)+np.sum(Occ == 1)
n1=np.sum(Occ == 1)
print(" Number of singly occupied orbitals {} ".format(n1))
nVir = np.sum(Occ == 0)
assert (nOcc + nVir == nOrb)
print(" Number of occupied orbitals {} ".format(nOcc))
print(" Number of unoccupied orbitals {} ".format(nVir))
# Compute Fock matrix from orbital eigenvalues (SCF has to be fully converged)
Fock = mdot(S1, C, np.diag(Eps), C.T, S1.T)
# Rdm = mdot(C, np.diag(Occ), C.T)
COcc = C[:,:nOcc]
CVir = C[:,nOcc:]
# Smh1 = MakeSmh(S1)
# Sh1 = np.dot(S1, Smh1)
# Compute IAO basis (non-orthogonal). Will be used to identify the
# pi-MO-space of the target atom groups.
CIb = MakeIaosRaw(COcc, S1, S2, S12)
# CIbOcc = np.linalg.lstsq(CIb, COcc)[0]
# Err = np.dot(Sh1, np.dot(CIb, CIbOcc) - COcc)
# check orthogonality of IAO basis occupied orbitals
SIb = mdot(CIb.T, S1, CIb)
# SIbOcc = mdot(CIbOcc.T, SIb, CIbOcc)
nIb = SIb.shape[0]
if 0:
# make the a representation of the virtual valence space.
# SmhIb = MakeSmh(SIb)
nIbVir = nIb - nOcc # number of virtual valence orbitals
CTargetIb = CIb
STargetIb = mdot(CTargetIb.T, S1, CTargetIb)
SmhTargetIb = MakeSmh(STargetIb)
STargetIbVir = mdot(SmhTargetIb, CTargetIb.T, S1, CVir)
U,sig,Vt = np.linalg.svd(STargetIbVir, full_matrices=False)
print(" Number of target MINAO basis fn {}\n".format(CTargetIb.shape[1]))
print("SIbVir Singular Values (n={})".format(len(sig)))
print(" [{}]".format(', '.join('{:.4f}'.format(k) for k in sig)))
assert (np.abs(sig[nIbVir-1] - 1.0) < 1e-4)
assert (np.abs(sig[nIbVir] - 0.0) < 1e-4)
# for pi systems: do it like here -^ for the virtuals, but
# for COcc and CVir both. What we should do is:
# - instead of CIb, we use CTargetIb = CIb * (AO-linear-comb-matrix)
# - instead of SmhIb, we use its corresponding target-AO overlap matrix:
# SmhTargetIb = MakeSmh(mdot(CTargetIb.T, S1, CTargetIb))
# - instead of using nOcc/nVir/nIb to determine the target number
# of orbitals, we obtain the target number of pi electrons by counting
# their subset from the selected main group atoms (see CHEM 408 u11).
# From this determine the number of occupied pi-orbitals this system
# is supposed to have, and virtual orbitals as nTargetIb - nPiOcc
# The AoMix matrix is made from the pi-system's inertial tensor to get the
# local z-direction, and then linearly-combining this onto the highest-N p-AOs.
CActOcc = []
CActVir = []
# add pi-HOMOs and pi-LUMOs
CFragOcc, CFragVir,nOccOrbExpected,nVirtOrbExpected = MakePiSystemOrbitals(
"Pi-System", PiAtomsList, None, Elements,Coords, CIb, Shells, S1, S12, S2, Fock, COcc, CVir)
if (nPiOcc is None):
for i in range(1,nOccOrbExpected+1):
CActOcc.append(CFragOcc[:,-i])
else:
for i in range(1,nPiOcc+1):
CActOcc.append(CFragOcc[:,-i])
if (nPiVirt is None):
for j in range(nVirtOrbExpected):
CActVir.append(CFragVir[:,j])
else:
for j in range(nPiVirt):
CActVir.append(CFragVir[:,j])
print("\n -- Joining active spaces")
if (mol.spin==0):
nElec = 2*len(CActOcc)
else:
nElec = 2*len(CActOcc)-n1
CAct = np.array(CActOcc + CActVir).T
if 0:
# orthogonalize and semi-canonicalize
SAct = mdot(CAct.T, S1, CAct)
ew, ev = np.linalg.eigh(SAct)
print(" CAct initial overlap (if all ~approx 1, then the initial "
"active orbitals are near-orthogonal. That is good.)")
print(" [{}]".format(', '.join('{:.4f}'.format(k) for k in ew)))
CAct = np.dot(CAct, MakeSmh(SAct))
CAct = SemiCanonicalize(CAct, Fock, S1, "active")
print(" Number of Active Electrons {} ".format(nElec))
print(" Number of Active Orbitals {} ".format( CAct.shape[1]))
# make new non-active occupied and non-active virtual orbitals,
# in order to rebuild a full MO matrix.
def MakeInactiveSpace(Name, CActList, COrb1):
CAct = np.array(CActList).T
SAct = mdot(CAct.T, S1, CAct)
CAct = np.dot(CAct, MakeSmh(SAct))
SActMo = mdot(COrb1.T, S1.T, CAct, CAct.T, S1, COrb1)
ew, ev = np.linalg.eigh(SActMo) # small ews first (orbs not overlapping with active orbs).
nRest = COrb1.shape[1] - CAct.shape[1]
if 0:
for i in range(len(ew)):
print("{:4} {:15.6f} {}".format(i,ew[i], i == nRest))
assert (np.abs(ew[nRest-1] - 0.0) < 1e-8)
assert (np.abs(ew[nRest] - 1.0) < 1e-8)
CNewOrb1 = np.dot(COrb1, ev[:,:nRest])
return SemiCanonicalize(CNewOrb1, Fock, S1, Name, Print=False)
CNewClo = MakeInactiveSpace("NewClo", CActOcc, COcc)
CNewExt = MakeInactiveSpace("NewVir", CActVir, CVir)
# re-orthogonalize (should be orthogonal already, but just to be sure).
COrbNew = np.hstack([CAct])
if 1:
COrbNew = np.hstack([CNewClo, CAct, CNewExt])
SMo = mdot(COrbNew.T, S1, COrbNew)
COrbNew = np.dot(COrbNew, MakeSmh(SMo))
print(" Number of Core Orbitals {} ".format(CNewClo.shape[1]))
print(" Number of Virtual Orbitals {} ".format(CNewExt.shape[1]))
print(" Total Number of Orbitals {} ".format(COrbNew.shape[1]))
return CNewClo.shape[1],CAct.shape[1],CNewExt.shape[1],nElec, COrbNew
[docs]
def rmsd(a, b = None):
if b is None:
return np.mean(a.flatten()**2)**.5
else:
return rmsd(a - b)
[docs]
def MakeSmh(S):
ew,ev = np.linalg.eigh(S)
assert (np.all(ew > 1e-10))
v = ev * (ew**-0.25)[np.newaxis,:]
return np.dot(v, v.T)
[docs]
def MakeIaosRaw(COcc, S1, S2, S12):
# calculate the molecule-intrinsic atomic orbital (IAO) basis
# ref: [1] Knizia, J. Chem. Theory Comput., http://dx.doi.org/10.1021/ct400687b
# This is the "Simple/2014" version from ibo-ref at sites.psu.edu/knizia/software
assert (S1.shape[0] == S1.shape[1] and S1.shape[0] == S12.shape[0])
assert (S2.shape[0] == S2.shape[1] and S2.shape[0] == S12.shape[1])
assert (COcc.shape[0] == S1.shape[0] and COcc.shape[1] <= S2.shape[0])
P12 = la.solve(S1, S12) # P12 = S1^{-1} S12
COcc2 = mdot(S12.T, COcc) # O(N m^2)
CTil = la.solve(S2, COcc2) # O(m^3)
STil = mdot(COcc2.T, CTil) # O(m^3)
CTil2Bar = la.solve(STil, CTil.T).T # O(m^3)
T4 = COcc - mdot(P12, CTil2Bar) # O(N m^2)
CIb = P12 + mdot(T4, COcc2.T) # O(N m^2)
return CIb
[docs]
def GetPzOrientation(iTargetAtoms, Coords_, Elements_):
# make a xyz vector pointing out of the plane containing the target atoms.
assert (Coords_.shape[1] == 3)
Coords = Coords_[iTargetAtoms,:].copy()
Elements = Elements_[iTargetAtoms]
# get center of mass coordinates
vCom = np.mean(Coords, axis=0)
Coords -= vCom
I = np.zeros((3,3))
for vCoord in Coords:
I += np.outer(vCoord, vCoord)
# diagonalize inertial-like tensor; large eigenvalues first.
ew,ev = np.linalg.eigh(-I)
ew *= -1
print("Target Atom List:")
print("[{}]".format(', '.join('{}'.format(k) for k in iTargetAtoms)))
print("Elements :")
print("[{}]".format(', '.join('{}'.format(k) for k in Elements)))
# direction of smallest spatial extend (last one) should be
# the one indicating the pz-orbital direction.
vPz = ev[:,2]
return vPz
[docs]
def FindValenceAoIndices(iAt, Shells, TargetL):
ipxyz = None
iFn0 = 0
for Atom in range(len(Shells)):
for AtomL in range(len(Shells[Atom])):
if Atom == iAt:
if Shells[Atom][AtomL][0] == TargetL:
# this is a generally contracted p-shell on atom iAt.
assert (ipxyz is None) # should be only one per atom, I think.
nCGTO = Shells[Atom][AtomL][1]
nSphComp = 2*TargetL + 1
iFnHighestP = iFn0 + nSphComp*(nCGTO-1)
ipxyz = np.array([(iFnHighestP + o) for o in range(nSphComp)])
iFn0 += Shells[Atom][AtomL][2]
assert (ipxyz is not None)
return ipxyz
[docs]
def MakePzMinaoVectors(iTargetAtoms, vPz, Shells):
# nIb = len(Shells)*len(Shells[Atom])*len( Shells[AtomShell][AtomL])
nIb=0
for Atom in range(len(Shells)):
for AtomL in range(len(Shells[Atom])):
nIb += Shells[Atom][AtomL][2]
nVec = len(iTargetAtoms)
AoVecs = np.zeros((nIb, nVec))
for (iVec, iAt) in enumerate(iTargetAtoms):
ipxyz = FindValenceAoIndices(iAt, Shells, TargetL = 1)
AoVec = np.zeros(nIb)
AoVec[ipxyz] = vPz
# each atom brings one pz orbital -> iVec = iiAt
AoVecs[:,iVec] = AoVec
return AoVecs
# table of numbers of pi electrons each main group element
# typically brings into a pi-system
_NumPiElec = {
"C": 1,
"B": 0,
"N1e": 1,
"N2e": 2,
"O1e": 1,
"O2e": 2,
"F": 2,
"Si":1,
"P1e": 1,
"P2e": 2,
"S1e": 1,
"S2e": 2,
"C2e": 2,
}
_CovalentRadii = {1: 38.0, 2: 32.0, 3: 134.0, 4: 90.0, 5: 82.0,
6: 77.0, 7: 75.0, 8: 73.0, 9: 71.0, 10: 69.0,
11: 154.0, 12: 130.0, 13: 118.0, 14: 111.0, 15: 106.0,
16: 102.0, 17: 99.0, 18: 97.0, 19: 196.0, 20: 174.0,
21: 144.0, 22: 136.0, 23: 125.0, 24: 127.0, 25: 139.0,
26: 125.0, 27: 126.0, 28: 121.0, 29: 138.0, 30: 131.0,
31: 126.0, 32: 122.0, 33: 119.0, 34: 116.0, 35: 114.0,
36: 110.0, 37: 211.0, 38: 192.0, 39: 162.0, 40: 148.0,
41: 137.0, 42: 145.0, 43: 156.0, 44: 126.0, 45: 135.0,
46: 131.0, 47: 153.0, 48: 148.0, 49: 144.0, 50: 141.0,
51: 138.0, 52: 135.0, 53: 133.0, 54: 130.0, 55: 225.0,
56: 198.0, 57: 169.0, 58: None, 59: None, 60: None,
61: None, 62: None, 63: None, 64: None, 65: None,
66: None, 67: None, 68: None, 69: None, 70: None,
71: 160.0, 72: 150.0, 73: 138.0, 74: 146.0, 75: 159.0,
76: 128.0, 77: 137.0, 78: 128.0, 79: 144.0, 80: 149.0,
81: 148.0, 82: 147.0, 83: 146.0, 84: None, 85: None,
86: 145.0, 87: None, 88: None, 89: None, 90: None,
91: None, 92: None, 93: None, 94: None, 95: None,
96: None, 97: None, 98: None, 99: None, 100: None,
101: None, 102: None, 103: None, 104: None, 105: None,
106: None, 107: None, 108: None, 109: None, 110: None,
111: None, 112: None, 113: None, 114: None, 115: None, 116: None}
# ^- embedded now to remove dependency on finding .txt path.
[docs]
def GetCovalentRadius(At):
# get covalent radius in pm
rcov = _CovalentRadii[At.iElement]
assert (rcov is not None)
# convert to bohr radii
ToAng = 0.52917721092
return (0.01 * rcov)/ToAng
[docs]
def GetNumPiElec(iAt, Elements,Coords):
# deal with C, F, B, Si, and the atoms which have explicit tages assigned.
At = Elements[iAt]
nPiTab = _NumPiElec.get(At, None)
if nPiTab is not None:
return nPiTab
print(" ...determining formal number of pi-electrons for {} (not in table).".format(At.Label))
# find number of bonded partners
iBonded = []
for (jAt, AtJ) in enumerate(Elements):
if iAt == jAt:
continue
rij = np.sum((Coords[At] - Coords[AtJ])**2)**.5
if rij < 1.3 * (GetCovalentRadius(At) + GetCovalentRadius(AtJ)):
iBonded.append(jAt)
nBonds = len(iBonded)
if At in ["N", "P"]: # 5 valence electrons
if nBonds == 2:
# two sigma bonds (1e each) + 1 sigma lone pair (2e) -> 1e left for pi system
return 1
elif nBonds == 3:
# three sigma bonds (1e each) -> 2e left for pi system
return 2
if At in ["O","S"]: # 6 valence electrons
if nBonds == 1:
# one sigma bond (1e) + 2 lone pairs in x/y planes (2e each) -> 1e left for pi system
return 1
elif nBonds == 3:
# two sigma bond (1e each) + 1 lone pairs in x/y planes (2e) -> 2e left for pi system
return 2
raise Exception("GetNumPiElec({}): {} with {} bonds not tabulated :(".format(iAt, At, nBonds))
[docs]
def SemiCanonicalize(COut, Fock, S1, Name, Print=True):
nOrb = COut.shape[1]
SOrb = mdot(COut.T, S1, COut)
fErr = rmsd(SOrb - np.eye(nOrb))
if fErr > 1e-8:
print(" SemiCan: Orbital deviation from orthogonality {:8.2e} \n".format(fErr))
raise Exception("Orbitals not orthogonal.")
FockMo = mdot(COut.T, Fock, COut)
ew,ev = np.linalg.eigh(FockMo)
if Print:
print(" Semi-canonicalized orbital energies ({} subspace)".format(Name))
print(" [{}]".format(', '.join('{:.4f}'.format(k) for k in ew)))
return np.dot(COut, ev)
[docs]
def MakeOverlappingOrbSubspace(Space, Name, COrb, nOrbExpected, CTargetIb, S1, Fock):
print("\n -- Constructing MO subspace space {}/{}.".format(Space, Name))
STargetIb = mdot(CTargetIb.T, S1, CTargetIb)
SmhTargetIb = MakeSmh(STargetIb)
STargetIbVir = mdot(SmhTargetIb, CTargetIb.T, S1, COrb)
U,sig,Vt = np.linalg.svd(STargetIbVir, full_matrices=False)
print(" S[{},{}] singular Values (n={}, nThr={})".format(Space, Name, len(sig), nOrbExpected))
print(" [{}]".format(', '.join('{:.4f}'.format(k) for k in sig)))
print(" Sigma[{}] {:.4f} (should be 1)".format(nOrbExpected-1, sig[nOrbExpected-1]))
if nOrbExpected < len(sig):
print(" Sigma[{}] {:.4f} (should be 0)".format(nOrbExpected, sig[nOrbExpected]))
if sig[nOrbExpected-1] < 0.8 or sig[nOrbExpected] > 0.5:
raise Exception("{} orbital construction okay?".format(Space))
V = Vt.T
COut = np.dot(COrb, V[:,:nOrbExpected])
return SemiCanonicalize(COut, Fock, S1, Name)
# Use these to add 'stuff' to the proto-active space.
[docs]
def MakePiSystemOrbitals(TargetName, iTargetAtomsForPlane_,
iTargetAtomsForBasis_,Elements,Coords, CIb, Shells,
S1, S12, S2, Fock, COcc, CVir):
print("\n *** TARGET: {}\n".format(TargetName))
# convert from 1-based atom indices to 0-based indices.
iTargetAtomsForPlane = np.array(iTargetAtomsForPlane_) - 1
if iTargetAtomsForBasis_ is None:
iTargetAtomsForBasis = iTargetAtomsForPlane
else:
iTargetAtomsForBasis = np.array(iTargetAtomsForBasis_) - 1
vPz = GetPzOrientation(iTargetAtomsForPlane, Coords, Elements)
CAoMix = MakePzMinaoVectors(iTargetAtomsForBasis, vPz, Shells)
nPiElec = 0
for iAt in iTargetAtomsForBasis:
nPiElec += GetNumPiElec(iAt, Elements,Coords)
print(" Formal number of elec in pi-system {} ".format(nPiElec))
if (nPiElec % 2) != 0:
raise Exception("Got {} pi electrons. Shouldn't this be an even number?".format(nPiElec))
# idea: make full pi system for target AOs (via IAOs)
# then we do know how many occupied and virtual pi
# orbitals there are supposed to be (by counting the contributed pi
# electrons from the constituent main group atoms)
# we then isolate those, canonicalize them, and take
# the frontier MO subset.
CTargetIb = np.dot(CIb, CAoMix)
nTargetIb = CAoMix.shape[1]
print(" Number of target MINAO basis fn {} ".format(nTargetIb))
print(" size of CIb1 {} ".format(CIb.shape[0]))
print(" size of CAomix1 {} ".format(CAoMix.shape[0]))
print(" size of CIb2 {} ".format(CIb.shape[1]))
print(" size of CAomix2 {} ".format(CAoMix.shape[1]))
nOccOrbExpected=nPiElec//2
nVirtOrbExpected=nTargetIb - nPiElec//2
CPiOcc = MakeOverlappingOrbSubspace("Pi", "Occ", COcc, nOccOrbExpected, CTargetIb, S1, Fock)
CPiVir = MakeOverlappingOrbSubspace("Pi", "Vir", CVir, nVirtOrbExpected, CTargetIb, S1, Fock)
return CPiOcc, CPiVir,nOccOrbExpected,nVirtOrbExpected