forked from pyscf/pyscf
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathavas.py
198 lines (174 loc) · 7.07 KB
/
avas.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
#!/usr/bin/env python
#
# Author: Qiming Sun <[email protected]>
# Elvira R. Sayfutyarova <[email protected]>
#
'''
Automated construction of molecular active spaces from atomic valence orbitals.
Ref. arXiv:1701.07862 [physics.chem-ph]
'''
import re
from functools import reduce
import numpy
import scipy.linalg
from pyscf import lib
from pyscf import gto
from pyscf import scf
from pyscf.lib import logger
def kernel(mf, aolabels, threshold=.2, minao='minao', with_iao=False,
openshelloption=2, canonicalize=True, verbose=None):
'''AVAS method to construct mcscf active space.
Ref. arXiv:1701.07862 [physics.chem-ph]
Args:
mf : an :class:`SCF` object
aolabels : string or a list of strings
AO labels for AO active space
Kwargs:
threshold : float
Tructing threshold of the AO-projector above which AOs are kept in
the active space.
minao : str
A reference AOs for AVAS.
with_iao : bool
Whether to use IAO localization to construct the reference active AOs.
openshelloption : int
How to handle singly-occupied orbitals in the active space. The
singly-occupied orbitals are projected as part of alpha orbitals
if openshelloption=2, or completely kept in active space if
openshelloption=3. See Section III.E option 2 or 3 of the
reference paper for more details.
canonicalize : bool
Orbitals defined in AVAS method are local orbitals. Symmetrizing
the core, active and virtual space.
Returns:
active-space-size, #-active-electrons, orbital-initial-guess-for-CASCI/CASSCF
Examples:
>>> from pyscf import gto, scf, mcscf
>>> from pyscf.mcscf import avas
>>> mol = gto.M(atom='Cr 0 0 0; Cr 0 0 1.6', basis='ccpvtz')
>>> mf = scf.RHF(mol).run()
>>> ncas, nelecas, mo = avas.avas(mf, ['Cr 3d', 'Cr 4s'])
>>> mc = mcscf.CASSCF(mf, ncas, nelecas).run(mo)
'''
if isinstance(verbose, logger.Logger):
log = verbose
elif verbose is not None:
log = logger.Logger(mf.stdout, verbose)
else:
log = logger.Logger(mf.stdout, mf.verbose)
mol = mf.mol
log.info('\n** AVAS **')
if isinstance(mf, scf.uhf.UHF):
log.note('UHF/UKS object is found. AVAS takes alpha orbitals only')
mo_coeff = mf.mo_coeff[0]
mo_occ = mf.mo_occ[0]
mo_energy = mf.mo_energy[0]
assert(openshelloption != 1)
else:
mo_coeff = mf.mo_coeff
mo_occ = mf.mo_occ
mo_energy = mf.mo_energy
nocc = numpy.count_nonzero(mo_occ != 0)
ovlp = mol.intor_symmetric('cint1e_ovlp_sph')
log.info(' Total number of HF MOs is equal to %d' ,mo_coeff.shape[1])
log.info(' Number of occupied HF MOs is equal to %d', nocc)
mol = mf.mol
pmol = mol.copy()
pmol.atom = mol._atom
pmol.unit = 'B'
pmol.symmetry = False
pmol.basis = minao
pmol.build(False, False)
if isinstance(aolabels, str):
aolabels = re.sub(' +', ' ', aolabels.strip(), count=1)
baslst = [i for i,s in enumerate(pmol.spherical_labels(1))
if aolabels in s]
elif isinstance(aolabels[0], str):
aolabels = [re.sub(' +', ' ', x.strip(), count=1) for x in aolabels]
baslst = [i for i,s in enumerate(pmol.spherical_labels(1))
if any(x in s for x in aolabels)]
else:
raise RuntimeError
baslst = numpy.asarray(baslst)
log.info('reference AO indices for %s %s: %s', minao, aolabels, baslst)
if with_iao:
from pyscf.lo import iao
c = iao.iao(mol, mo_coeff[:,:nocc], minao)[:,baslst]
s2 = reduce(numpy.dot, (c.T, ovlp, c))
s21 = reduce(numpy.dot, (c.T, ovlp, mo_coeff))
else:
s2 = pmol.intor_symmetric('cint1e_ovlp_sph')[baslst][:,baslst]
s21 = gto.intor_cross('cint1e_ovlp_sph', pmol, mol)[baslst]
s21 = numpy.dot(s21, mo_coeff)
sa = s21.T.dot(scipy.linalg.solve(s2, s21, sym_pos=True))
if openshelloption == 2:
wocc, u = numpy.linalg.eigh(sa[:nocc,:nocc])
log.info('Option 2: threshold %s', threshold)
ncas_occ = (wocc > threshold).sum()
nelecas = mol.nelectron - (wocc < threshold).sum() * 2
mocore = mo_coeff[:,:nocc].dot(u[:,wocc<threshold])
mocas = mo_coeff[:,:nocc].dot(u[:,wocc>threshold])
wvir, u = numpy.linalg.eigh(sa[nocc:,nocc:])
ncas_vir = (wvir > threshold).sum()
mocas = numpy.hstack((mocas, mo_coeff[:,nocc:].dot(u[:,wvir>threshold])))
movir = mo_coeff[:,nocc:].dot(u[:,wvir<threshold])
ncas = mocas.shape[1]
elif openshelloption == 3:
docc = nocc - mol.spin
wocc, u = numpy.linalg.eigh(sa[:docc,:docc])
log.info('Option 3: threshold %s, num open shell %d', threshold, mol.spin)
ncas_occ = (wocc > threshold).sum()
nelecas = mol.nelectron - (wocc < threshold).sum() * 2
mocore = mo_coeff[:,:docc].dot(u[:,wocc<threshold])
mocas = mo_coeff[:,:docc].dot(u[:,wocc>threshold])
wvir, u = numpy.linalg.eigh(sa[nocc:,nocc:])
ncas_vir = (wvir > threshold).sum()
mocas = numpy.hstack((mocas, mo_coeff[:,docc:nocc],
mo_coeff[:,nocc:].dot(u[:,wvir>threshold])))
movir = mo_coeff[:,nocc:].dot(u[:,wvir<threshold])
ncas = mocas.shape[1]
log.debug('projected occ eig %s', wocc[::-1])
log.debug('projected vir eig %s', wvir[::-1])
log.info('Active from occupied = %d , eig %s', ncas_occ, wocc[wocc>threshold][::-1])
log.info('Inactive from occupied = %d', mocore.shape[1])
log.info('Active from unoccupied = %d , eig %s', ncas_vir, wvir[wvir>threshold][::-1])
log.info('Inactive from unoccupied = %d', movir.shape[1])
log.info('Dimensions of active %d', ncas)
nalpha = (nelecas + mol.spin) // 2
nbeta = nelecas - nalpha
log.info('# of alpha electrons %d', nalpha)
log.info('# of beta electrons %d', nbeta)
if canonicalize:
from pyscf.mcscf import dmet_cas
def trans(c):
if c.shape[1] == 0:
return c
else:
csc = reduce(numpy.dot, (c.T, ovlp, mo_coeff))
fock = numpy.dot(csc*mo_energy, csc.T)
e, u = scipy.linalg.eigh(fock)
return dmet_cas.symmetrize(mol, e, numpy.dot(c, u), ovlp, log)
mo = numpy.hstack([trans(mocore), trans(mocas), trans(movir)])
else:
mo = numpy.hstack((mocore, mocas, movir))
return ncas, nelecas, mo
avas = kernel
if __name__ == '__main__':
from pyscf import gto
from pyscf import scf
from pyscf import mcscf
mol = gto.M(
verbose = 0,
atom = '''
H 0.000000, 0.500000, 1.5
O 0.000000, 0.000000, 1.
O 0.000000, 0.000000, -1.
H 0.000000, -0.500000, -1.5''',
basis = 'ccpvdz',
)
mf = scf.UHF(mol)
mf.scf()
ncas, nelecas, mo = avas(mf, 'O 2p', verbose=4)
mc = mcscf.CASSCF(mf, ncas, nelecas).set(verbose=4)
emc = mc.kernel(mo)[0]
print(emc, -150.51496582534054)