-
-
Notifications
You must be signed in to change notification settings - Fork 550
/
hf.py
2022 lines (1698 loc) · 70.6 KB
/
hf.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
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
#!/usr/bin/env python
# Copyright 2014-2020 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: Qiming Sun <osirpt.sun@gmail.com>
#
'''
Hartree-Fock
'''
import sys
import tempfile
import time
from functools import reduce
import numpy
import scipy.linalg
import h5py
from pyscf import gto
from pyscf import lib
from pyscf.lib import logger
from pyscf.scf import diis
from pyscf.scf import _vhf
from pyscf.scf import chkfile
from pyscf.data import nist
from pyscf import __config__
WITH_META_LOWDIN = getattr(__config__, 'scf_analyze_with_meta_lowdin', True)
PRE_ORTH_METHOD = getattr(__config__, 'scf_analyze_pre_orth_method', 'ANO')
MO_BASE = getattr(__config__, 'MO_BASE', 1)
TIGHT_GRAD_CONV_TOL = getattr(__config__, 'scf_hf_kernel_tight_grad_conv_tol', True)
MUTE_CHKFILE = getattr(__config__, 'scf_hf_SCF_mute_chkfile', False)
# For code compatibility in python-2 and python-3
if sys.version_info >= (3,):
unicode = str
def kernel(mf, conv_tol=1e-10, conv_tol_grad=None,
dump_chk=True, dm0=None, callback=None, conv_check=True, **kwargs):
'''kernel: the SCF driver.
Args:
mf : an instance of SCF class
mf object holds all parameters to control SCF. One can modify its
member functions to change the behavior of SCF. The member
functions which are called in kernel are
| mf.get_init_guess
| mf.get_hcore
| mf.get_ovlp
| mf.get_veff
| mf.get_fock
| mf.get_grad
| mf.eig
| mf.get_occ
| mf.make_rdm1
| mf.energy_tot
| mf.dump_chk
Kwargs:
conv_tol : float
converge threshold.
conv_tol_grad : float
gradients converge threshold.
dump_chk : bool
Whether to save SCF intermediate results in the checkpoint file
dm0 : ndarray
Initial guess density matrix. If not given (the default), the kernel
takes the density matrix generated by ``mf.get_init_guess``.
callback : function(envs_dict) => None
callback function takes one dict as the argument which is
generated by the builtin function :func:`locals`, so that the
callback function can access all local variables in the current
envrionment.
Returns:
A list : scf_conv, e_tot, mo_energy, mo_coeff, mo_occ
scf_conv : bool
True means SCF converged
e_tot : float
Hartree-Fock energy of last iteration
mo_energy : 1D float array
Orbital energies. Depending the eig function provided by mf
object, the orbital energies may NOT be sorted.
mo_coeff : 2D array
Orbital coefficients.
mo_occ : 1D array
Orbital occupancies. The occupancies may NOT be sorted from large
to small.
Examples:
>>> from pyscf import gto, scf
>>> mol = gto.M(atom='H 0 0 0; H 0 0 1.1', basis='cc-pvdz')
>>> conv, e, mo_e, mo, mo_occ = scf.hf.kernel(scf.hf.SCF(mol), dm0=numpy.eye(mol.nao_nr()))
>>> print('conv = %s, E(HF) = %.12f' % (conv, e))
conv = True, E(HF) = -1.081170784378
'''
if 'init_dm' in kwargs:
raise RuntimeError('''
You see this error message because of the API updates in pyscf v0.11.
Keyword argument "init_dm" is replaced by "dm0"''')
cput0 = (time.clock(), time.time())
if conv_tol_grad is None:
conv_tol_grad = numpy.sqrt(conv_tol)
logger.info(mf, 'Set gradient conv threshold to %g', conv_tol_grad)
mol = mf.mol
if dm0 is None:
dm = mf.get_init_guess(mol, mf.init_guess)
else:
dm = dm0
h1e = mf.get_hcore(mol)
vhf = mf.get_veff(mol, dm)
e_tot = mf.energy_tot(dm, h1e, vhf)
logger.info(mf, 'init E= %.15g', e_tot)
scf_conv = False
mo_energy = mo_coeff = mo_occ = None
s1e = mf.get_ovlp(mol)
cond = lib.cond(s1e)
logger.debug(mf, 'cond(S) = %s', cond)
if numpy.max(cond)*1e-17 > conv_tol:
logger.warn(mf, 'Singularity detected in overlap matrix (condition number = %4.3g). '
'SCF may be inaccurate and hard to converge.', numpy.max(cond))
# Skip SCF iterations. Compute only the total energy of the initial density
if mf.max_cycle <= 0:
fock = mf.get_fock(h1e, s1e, vhf, dm) # = h1e + vhf, no DIIS
mo_energy, mo_coeff = mf.eig(fock, s1e)
mo_occ = mf.get_occ(mo_energy, mo_coeff)
return scf_conv, e_tot, mo_energy, mo_coeff, mo_occ
if isinstance(mf.diis, lib.diis.DIIS):
mf_diis = mf.diis
elif mf.diis:
assert issubclass(mf.DIIS, lib.diis.DIIS)
mf_diis = mf.DIIS(mf, mf.diis_file)
mf_diis.space = mf.diis_space
mf_diis.rollback = mf.diis_space_rollback
else:
mf_diis = None
if dump_chk and mf.chkfile:
# Explicit overwrite the mol object in chkfile
# Note in pbc.scf, mf.mol == mf.cell, cell is saved under key "mol"
chkfile.save_mol(mol, mf.chkfile)
# A preprocessing hook before the SCF iteration
mf.pre_kernel(locals())
cput1 = logger.timer(mf, 'initialize scf', *cput0)
for cycle in range(mf.max_cycle):
dm_last = dm
last_hf_e = e_tot
fock = mf.get_fock(h1e, s1e, vhf, dm, cycle, mf_diis)
mo_energy, mo_coeff = mf.eig(fock, s1e)
mo_occ = mf.get_occ(mo_energy, mo_coeff)
dm = mf.make_rdm1(mo_coeff, mo_occ)
# attach mo_coeff and mo_occ to dm to improve DFT get_veff efficiency
dm = lib.tag_array(dm, mo_coeff=mo_coeff, mo_occ=mo_occ)
vhf = mf.get_veff(mol, dm, dm_last, vhf)
e_tot = mf.energy_tot(dm, h1e, vhf)
# Here Fock matrix is h1e + vhf, without DIIS. Calling get_fock
# instead of the statement "fock = h1e + vhf" because Fock matrix may
# be modified in some methods.
fock = mf.get_fock(h1e, s1e, vhf, dm) # = h1e + vhf, no DIIS
norm_gorb = numpy.linalg.norm(mf.get_grad(mo_coeff, mo_occ, fock))
if not TIGHT_GRAD_CONV_TOL:
norm_gorb = norm_gorb / numpy.sqrt(norm_gorb.size)
norm_ddm = numpy.linalg.norm(dm-dm_last)
logger.info(mf, 'cycle= %d E= %.15g delta_E= %4.3g |g|= %4.3g |ddm|= %4.3g',
cycle+1, e_tot, e_tot-last_hf_e, norm_gorb, norm_ddm)
if callable(mf.check_convergence):
scf_conv = mf.check_convergence(locals())
elif abs(e_tot-last_hf_e) < conv_tol and norm_gorb < conv_tol_grad:
scf_conv = True
if dump_chk:
mf.dump_chk(locals())
if callable(callback):
callback(locals())
cput1 = logger.timer(mf, 'cycle= %d'%(cycle+1), *cput1)
if scf_conv:
break
if scf_conv and conv_check:
# An extra diagonalization, to remove level shift
#fock = mf.get_fock(h1e, s1e, vhf, dm) # = h1e + vhf
mo_energy, mo_coeff = mf.eig(fock, s1e)
mo_occ = mf.get_occ(mo_energy, mo_coeff)
dm, dm_last = mf.make_rdm1(mo_coeff, mo_occ), dm
dm = lib.tag_array(dm, mo_coeff=mo_coeff, mo_occ=mo_occ)
vhf = mf.get_veff(mol, dm, dm_last, vhf)
e_tot, last_hf_e = mf.energy_tot(dm, h1e, vhf), e_tot
fock = mf.get_fock(h1e, s1e, vhf, dm)
norm_gorb = numpy.linalg.norm(mf.get_grad(mo_coeff, mo_occ, fock))
if not TIGHT_GRAD_CONV_TOL:
norm_gorb = norm_gorb / numpy.sqrt(norm_gorb.size)
norm_ddm = numpy.linalg.norm(dm-dm_last)
conv_tol = conv_tol * 10
conv_tol_grad = conv_tol_grad * 3
if callable(mf.check_convergence):
scf_conv = mf.check_convergence(locals())
elif abs(e_tot-last_hf_e) < conv_tol or norm_gorb < conv_tol_grad:
scf_conv = True
logger.info(mf, 'Extra cycle E= %.15g delta_E= %4.3g |g|= %4.3g |ddm|= %4.3g',
e_tot, e_tot-last_hf_e, norm_gorb, norm_ddm)
if dump_chk:
mf.dump_chk(locals())
logger.timer(mf, 'scf_cycle', *cput0)
# A post-processing hook before return
mf.post_kernel(locals())
return scf_conv, e_tot, mo_energy, mo_coeff, mo_occ
def energy_elec(mf, dm=None, h1e=None, vhf=None):
r'''Electronic part of Hartree-Fock energy, for given core hamiltonian and
HF potential
... math::
E = \sum_{ij}h_{ij} \gamma_{ji}
+ \frac{1}{2}\sum_{ijkl} \gamma_{ji}\gamma_{lk} \langle ik||jl\rangle
Note this function has side effects which cause mf.scf_summary updated.
Args:
mf : an instance of SCF class
Kwargs:
dm : 2D ndarray
one-partical density matrix
h1e : 2D ndarray
Core hamiltonian
vhf : 2D ndarray
HF potential
Returns:
Hartree-Fock electronic energy and the Coulomb energy
Examples:
>>> from pyscf import gto, scf
>>> mol = gto.M(atom='H 0 0 0; H 0 0 1.1')
>>> mf = scf.RHF(mol)
>>> mf.scf()
>>> dm = mf.make_rdm1()
>>> scf.hf.energy_elec(mf, dm)
(-1.5176090667746334, 0.60917167853723675)
>>> mf.energy_elec(dm)
(-1.5176090667746334, 0.60917167853723675)
'''
if dm is None: dm = mf.make_rdm1()
if h1e is None: h1e = mf.get_hcore()
if vhf is None: vhf = mf.get_veff(mf.mol, dm)
e1 = numpy.einsum('ij,ji->', h1e, dm)
e_coul = numpy.einsum('ij,ji->', vhf, dm) * .5
mf.scf_summary['e1'] = e1.real
mf.scf_summary['e2'] = e_coul.real
logger.debug(mf, 'E1 = %s E_coul = %s', e1, e_coul)
return (e1+e_coul).real, e_coul
def energy_tot(mf, dm=None, h1e=None, vhf=None):
r'''Total Hartree-Fock energy, electronic part plus nuclear repulstion
See :func:`scf.hf.energy_elec` for the electron part
Note this function has side effects which cause mf.scf_summary updated.
'''
nuc = mf.energy_nuc()
e_tot = mf.energy_elec(dm, h1e, vhf)[0] + nuc
mf.scf_summary['nuc'] = nuc.real
return e_tot
def get_hcore(mol):
'''Core Hamiltonian
Examples:
>>> from pyscf import gto, scf
>>> mol = gto.M(atom='H 0 0 0; H 0 0 1.1')
>>> scf.hf.get_hcore(mol)
array([[-0.93767904, -0.59316327],
[-0.59316327, -0.93767904]])
'''
h = mol.intor_symmetric('int1e_kin')
if mol._pseudo:
# Although mol._pseudo for GTH PP is only available in Cell, GTH PP
# may exist if mol is converted from cell object.
from pyscf.gto import pp_int
h += pp_int.get_gth_pp(mol)
else:
h+= mol.intor_symmetric('int1e_nuc')
if len(mol._ecpbas) > 0:
h += mol.intor_symmetric('ECPscalar')
return h
def get_ovlp(mol):
'''Overlap matrix
'''
return mol.intor_symmetric('int1e_ovlp')
def init_guess_by_minao(mol):
'''Generate initial guess density matrix based on ANO basis, then project
the density matrix to the basis set defined by ``mol``
Returns:
Density matrix, 2D ndarray
Examples:
>>> from pyscf import gto, scf
>>> mol = gto.M(atom='H 0 0 0; H 0 0 1.1')
>>> scf.hf.init_guess_by_minao(mol)
array([[ 0.94758917, 0.09227308],
[ 0.09227308, 0.94758917]])
'''
from pyscf.scf import atom_hf
from pyscf.scf import addons
def minao_basis(symb, nelec_ecp):
occ = []
basis_ano = []
if gto.is_ghost_atom(symb):
return occ, basis_ano
stdsymb = gto.mole._std_symbol(symb)
basis_add = gto.basis.load('ano', stdsymb)
# coreshl defines the core shells to be removed in the initial guess
coreshl = gto.ecp.core_configuration(nelec_ecp)
#coreshl = (0,0,0,0) # it keeps all core electrons in the initial guess
for l in range(4):
ndocc, frac = atom_hf.frac_occ(stdsymb, l)
assert ndocc >= coreshl[l]
degen = l * 2 + 1
occ_l = [2,]*(ndocc-coreshl[l]) + [frac,]
occ.append(numpy.repeat(occ_l, degen))
basis_ano.append([l] + [b[:1] + b[1+coreshl[l]:ndocc+2]
for b in basis_add[l][1:]])
occ = numpy.hstack(occ)
if nelec_ecp > 0:
basis4ecp = [[] for i in range(4)]
for bas in mol._basis[symb]:
l = bas[0]
if l < 4:
basis4ecp[l].append(bas)
occ4ecp = []
for l in range(4):
nbas_l = sum((len(bas[1]) - 1) for bas in basis4ecp[l])
ndocc, frac = atom_hf.frac_occ(stdsymb, l)
ndocc -= coreshl[l]
assert ndocc <= nbas_l
occ_l = numpy.zeros(nbas_l)
occ_l[:ndocc] = 2
if frac > 0:
occ_l[ndocc] = frac
occ4ecp.append(numpy.repeat(occ_l, l * 2 + 1))
occ4ecp = numpy.hstack(occ4ecp)
basis4ecp = lib.flatten(basis4ecp)
# Compared to ANO valence basis, to check whether the ECP basis set has
# reasonable AO-character contraction. The ANO valence AO should have
# significant overlap to ECP basis if the ECP basis has AO-character.
atm1 = gto.Mole()
atm2 = gto.Mole()
atom = [[symb, (0.,0.,0.)]]
atm1._atm, atm1._bas, atm1._env = atm1.make_env(atom, {symb:basis4ecp}, [])
atm2._atm, atm2._bas, atm2._env = atm2.make_env(atom, {symb:basis_ano}, [])
atm1._built = True
atm2._built = True
s12 = gto.intor_cross('int1e_ovlp', atm1, atm2)
if abs(numpy.linalg.det(s12[occ4ecp>0][:,occ>0])) > .1:
occ, basis_ano = occ4ecp, basis4ecp
else:
logger.debug(mol, 'Density of valence part of ANO basis '
'will be used as initial guess for %s', symb)
return occ, basis_ano
# Issue 548
if any(gto.charge(mol.atom_symbol(ia)) > 96 for ia in range(mol.natm)):
logger.info(mol, 'MINAO initial guess is not available for super-heavy '
'elements. "atom" initial guess is used.')
return init_guess_by_atom(mol)
nelec_ecp_dic = dict([(mol.atom_symbol(ia), mol.atom_nelec_core(ia))
for ia in range(mol.natm)])
basis = {}
occdic = {}
for symb, nelec_ecp in nelec_ecp_dic.items():
occ_add, basis_add = minao_basis(symb, nelec_ecp)
occdic[symb] = occ_add
basis[symb] = basis_add
occ = []
new_atom = []
for ia in range(mol.natm):
symb = mol.atom_symbol(ia)
if not gto.is_ghost_atom(symb):
occ.append(occdic[symb])
new_atom.append(mol._atom[ia])
occ = numpy.hstack(occ)
pmol = gto.Mole()
pmol._atm, pmol._bas, pmol._env = pmol.make_env(new_atom, basis, [])
pmol._built = True
dm = addons.project_dm_nr2nr(pmol, numpy.diag(occ), mol)
# normalize eletron number
# s = mol.intor_symmetric('int1e_ovlp')
# dm *= mol.nelectron / (dm*s).sum()
return dm
def init_guess_by_1e(mol):
'''Generate initial guess density matrix from core hamiltonian
Returns:
Density matrix, 2D ndarray
'''
mf = RHF(mol)
return mf.init_guess_by_1e(mol)
def init_guess_by_atom(mol):
'''Generate initial guess density matrix from superposition of atomic HF
density matrix. The atomic HF is occupancy averaged RHF
Returns:
Density matrix, 2D ndarray
'''
from pyscf.scf import atom_hf
atm_scf = atom_hf.get_atm_nrhf(mol)
aoslice = mol.aoslice_by_atom()
atm_dms = []
for ia in range(mol.natm):
symb = mol.atom_symbol(ia)
if symb not in atm_scf:
symb = mol.atom_pure_symbol(ia)
if symb in atm_scf:
e_hf, e, c, occ = atm_scf[symb]
dm = numpy.dot(c*occ, c.conj().T)
else: # symb's basis is not specified in the input
nao_atm = aoslice[ia,3] - aoslice[ia,2]
dm = numpy.zeros((nao_atm, nao_atm))
atm_dms.append(dm)
dm = scipy.linalg.block_diag(*atm_dms)
if mol.cart:
cart2sph = mol.cart2sph_coeff(normalized='sp')
dm = reduce(numpy.dot, (cart2sph, dm, cart2sph.T))
for k, v in atm_scf.items():
logger.debug1(mol, 'Atom %s, E = %.12g', k, v[0])
return dm
def init_guess_by_huckel(mol):
'''Generate initial guess density matrix from a Huckel calculation based
on occupancy averaged atomic RHF calculations, doi:10.1021/acs.jctc.8b01089
Returns:
Density matrix, 2D ndarray
'''
mo_energy, mo_coeff = _init_guess_huckel_orbitals(mol)
mo_occ = get_occ(SCF(mol), mo_energy, mo_coeff)
return make_rdm1(mo_coeff, mo_occ)
def _init_guess_huckel_orbitals(mol):
'''Generate initial guess density matrix from a Huckel calculation based
on occupancy averaged atomic RHF calculations, doi:10.1021/acs.jctc.8b01089
Returns:
An 1D array for Huckel orbital energies and an 2D array for orbital coefficients
'''
from pyscf.scf import atom_hf
atm_scf = atom_hf.get_atm_nrhf(mol)
# GWH parameter value
Kgwh = 1.75
# Run atomic SCF calculations to get orbital energies, coefficients and occupations
at_e = []
at_c = []
at_occ = []
for ia in range(mol.natm):
symb = mol.atom_symbol(ia)
if symb not in atm_scf:
symb = mol.atom_pure_symbol(ia)
e_hf, e, c, occ = atm_scf[symb]
at_c.append(c)
at_e.append(e)
at_occ.append(occ)
# Count number of occupied orbitals
nocc = 0
for ia in range(mol.natm):
for iorb in range(len(at_occ[ia])):
if(at_occ[ia][iorb]>0.0):
nocc=nocc+1
# Number of basis functions
nbf = mol.nao_nr()
# Collect AO coefficients and energies
orb_E = numpy.zeros(nocc)
orb_C = numpy.zeros((nbf,nocc))
# Atomic basis info
aoslice = mol.aoslice_by_atom()
iocc = 0
for ia in range(mol.natm):
# First and last bf index
abeg = aoslice[ia, 2]
aend = aoslice[ia, 3]
for iorb in range(len(at_occ[ia])):
if(at_occ[ia][iorb]>0.0):
orb_C[abeg:aend,iocc] = at_c[ia][:,iorb]
orb_E[iocc] = at_e[ia][iorb]
iocc=iocc+1
# Overlap matrix
S = get_ovlp(mol)
# Atomic orbital overlap
orb_S = orb_C.transpose().dot(S).dot(orb_C)
# Build Huckel matrix
orb_H = numpy.zeros((nocc,nocc))
for io in range(nocc):
# Diagonal is just the orbital energies
orb_H[io,io] = orb_E[io]
for jo in range(io):
# Off-diagonal is given by GWH approximation
orb_H[io,jo] = 0.5*Kgwh*orb_S[io,jo]*(orb_E[io]+orb_E[jo])
orb_H[jo,io] = orb_H[io,jo]
# Energies and coefficients in the minimal orbital basis
mo_E, atmo_C = eig(orb_H, orb_S)
# and in the AO basis
mo_C = orb_C.dot(atmo_C)
return mo_E, mo_C
def init_guess_by_chkfile(mol, chkfile_name, project=None):
'''Read the HF results from checkpoint file, then project it to the
basis defined by ``mol``
Returns:
Density matrix, 2D ndarray
'''
from pyscf.scf import uhf
dm = uhf.init_guess_by_chkfile(mol, chkfile_name, project)
return dm[0] + dm[1]
def get_init_guess(mol, key='minao'):
'''Generate density matrix for initial guess
Kwargs:
key : str
One of 'minao', 'atom', 'huckel', 'hcore', '1e', 'chkfile'.
'''
return RHF(mol).get_init_guess(mol, key)
# eigenvalue of d is 1
def level_shift(s, d, f, factor):
r'''Apply level shift :math:`\Delta` to virtual orbitals
.. math::
:nowrap:
\begin{align}
FC &= SCE \\
F &= F + SC \Lambda C^\dagger S \\
\Lambda_{ij} &=
\begin{cases}
\delta_{ij}\Delta & i \in \text{virtual} \\
0 & \text{otherwise}
\end{cases}
\end{align}
Returns:
New Fock matrix, 2D ndarray
'''
dm_vir = s - reduce(numpy.dot, (s, d, s))
return f + dm_vir * factor
def damping(s, d, f, factor):
#dm_vir = s - reduce(numpy.dot, (s,d,s))
#sinv = numpy.linalg.inv(s)
#f0 = reduce(numpy.dot, (dm_vir, sinv, f, d, s))
dm_vir = numpy.eye(s.shape[0]) - numpy.dot(s, d)
f0 = reduce(numpy.dot, (dm_vir, f, d, s))
f0 = (f0+f0.conj().T) * (factor/(factor+1.))
return f - f0
# full density matrix for RHF
def make_rdm1(mo_coeff, mo_occ, **kwargs):
'''One-particle density matrix in AO representation
Args:
mo_coeff : 2D ndarray
Orbital coefficients. Each column is one orbital.
mo_occ : 1D ndarray
Occupancy
'''
mocc = mo_coeff[:,mo_occ>0]
# DO NOT make tag_array for dm1 here because this DM array may be modified and
# passed to functions like get_jk, get_vxc. These functions may take the tags
# (mo_coeff, mo_occ) to compute the potential if tags were found in the DM
# array and modifications to DM array may be ignored.
return numpy.dot(mocc*mo_occ[mo_occ>0], mocc.conj().T)
################################################
# for general DM
# hermi = 0 : arbitary
# hermi = 1 : hermitian
# hermi = 2 : anti-hermitian
################################################
def dot_eri_dm(eri, dm, hermi=0, with_j=True, with_k=True):
'''Compute J, K matrices in terms of the given 2-electron integrals and
density matrix:
J ~ numpy.einsum('pqrs,qp->rs', eri, dm)
K ~ numpy.einsum('pqrs,qr->ps', eri, dm)
Args:
eri : ndarray
8-fold or 4-fold ERIs or complex integral array with N^4 elements
(N is the number of orbitals)
dm : ndarray or list of ndarrays
A density matrix or a list of density matrices
Kwargs:
hermi : int
Whether J, K matrix is hermitian
| 0 : no hermitian or symmetric
| 1 : hermitian
| 2 : anti-hermitian
Returns:
Depending on the given dm, the function returns one J and one K matrix,
or a list of J matrices and a list of K matrices, corresponding to the
input density matrices.
Examples:
>>> from pyscf import gto, scf
>>> from pyscf.scf import _vhf
>>> mol = gto.M(atom='H 0 0 0; H 0 0 1.1')
>>> eri = _vhf.int2e_sph(mol._atm, mol._bas, mol._env)
>>> dms = numpy.random.random((3,mol.nao_nr(),mol.nao_nr()))
>>> j, k = scf.hf.dot_eri_dm(eri, dms, hermi=0)
>>> print(j.shape)
(3, 2, 2)
'''
dm = numpy.asarray(dm)
nao = dm.shape[-1]
if eri.dtype == numpy.complex128 or eri.size == nao**4:
eri = eri.reshape((nao,)*4)
dms = dm.reshape(-1,nao,nao)
vj = vk = None
if with_j:
vj = numpy.einsum('ijkl,xji->xkl', eri, dms)
vj = vj.reshape(dm.shape)
if with_k:
vk = numpy.einsum('ijkl,xjk->xil', eri, dms)
vk = vk.reshape(dm.shape)
else:
vj, vk = _vhf.incore(eri, dm.real, hermi, with_j, with_k)
if dm.dtype == numpy.complex128:
vs = _vhf.incore(eri, dm.imag, 0, with_j, with_k)
if with_j:
vj = vj + vs[0] * 1j
if with_k:
vk = vk + vs[1] * 1j
return vj, vk
def get_jk(mol, dm, hermi=1, vhfopt=None, with_j=True, with_k=True, omega=None):
'''Compute J, K matrices for all input density matrices
Args:
mol : an instance of :class:`Mole`
dm : ndarray or list of ndarrays
A density matrix or a list of density matrices
Kwargs:
hermi : int
Whether J, K matrix is hermitian
| 0 : not hermitian and not symmetric
| 1 : hermitian or symmetric
| 2 : anti-hermitian
vhfopt :
A class which holds precomputed quantities to optimize the
computation of J, K matrices
with_j : boolean
Whether to compute J matrices
with_k : boolean
Whether to compute K matrices
omega : float
Parameter of range-seperated Coulomb operator: erf( omega * r12 ) / r12.
If specified, integration are evaluated based on the long-range
part of the range-seperated Coulomb operator.
Returns:
Depending on the given dm, the function returns one J and one K matrix,
or a list of J matrices and a list of K matrices, corresponding to the
input density matrices.
Examples:
>>> from pyscf import gto, scf
>>> from pyscf.scf import _vhf
>>> mol = gto.M(atom='H 0 0 0; H 0 0 1.1')
>>> dms = numpy.random.random((3,mol.nao_nr(),mol.nao_nr()))
>>> j, k = scf.hf.get_jk(mol, dms, hermi=0)
>>> print(j.shape)
(3, 2, 2)
'''
dm = numpy.asarray(dm, order='C')
dm_shape = dm.shape
dm_dtype = dm.dtype
nao = dm_shape[-1]
if dm_dtype == numpy.complex128:
dm = numpy.vstack((dm.real, dm.imag)).reshape(-1,nao,nao)
hermi = 0
if omega is None:
vj, vk = _vhf.direct(dm, mol._atm, mol._bas, mol._env,
vhfopt, hermi, mol.cart, with_j, with_k)
else:
# The vhfopt of standard Coulomb operator can be used here as an approximate
# integral prescreening conditioner since long-range part Coulomb is always
# smaller than standard Coulomb. It's safe to filter LR integrals with the
# integral estimation from standard Coulomb.
with mol.with_range_coulomb(omega):
vj, vk = _vhf.direct(dm, mol._atm, mol._bas, mol._env,
vhfopt, hermi, mol.cart, with_j, with_k)
if dm_dtype == numpy.complex128:
if with_j:
vj = vj.reshape(2,-1)
vj = vj[0] + vj[1] * 1j
vj = vj.reshape(dm_shape)
if with_k:
vk = vk.reshape(2,-1)
vk = vk[0] + vk[1] * 1j
vk = vk.reshape(dm_shape)
return vj, vk
def get_veff(mol, dm, dm_last=None, vhf_last=None, hermi=1, vhfopt=None):
'''Hartree-Fock potential matrix for the given density matrix
Args:
mol : an instance of :class:`Mole`
dm : ndarray or list of ndarrays
A density matrix or a list of density matrices
Kwargs:
dm_last : ndarray or a list of ndarrays or 0
The density matrix baseline. If not 0, this function computes the
increment of HF potential w.r.t. the reference HF potential matrix.
vhf_last : ndarray or a list of ndarrays or 0
The reference HF potential matrix.
hermi : int
Whether J, K matrix is hermitian
| 0 : no hermitian or symmetric
| 1 : hermitian
| 2 : anti-hermitian
vhfopt :
A class which holds precomputed quantities to optimize the
computation of J, K matrices
Returns:
matrix Vhf = 2*J - K. Vhf can be a list matrices, corresponding to the
input density matrices.
Examples:
>>> import numpy
>>> from pyscf import gto, scf
>>> from pyscf.scf import _vhf
>>> mol = gto.M(atom='H 0 0 0; H 0 0 1.1')
>>> dm0 = numpy.random.random((mol.nao_nr(),mol.nao_nr()))
>>> vhf0 = scf.hf.get_veff(mol, dm0, hermi=0)
>>> dm1 = numpy.random.random((mol.nao_nr(),mol.nao_nr()))
>>> vhf1 = scf.hf.get_veff(mol, dm1, hermi=0)
>>> vhf2 = scf.hf.get_veff(mol, dm1, dm_last=dm0, vhf_last=vhf0, hermi=0)
>>> numpy.allclose(vhf1, vhf2)
True
'''
if dm_last is None:
vj, vk = get_jk(mol, numpy.asarray(dm), hermi, vhfopt)
return vj - vk * .5
else:
ddm = numpy.asarray(dm) - numpy.asarray(dm_last)
vj, vk = get_jk(mol, ddm, hermi, vhfopt)
return vj - vk * .5 + numpy.asarray(vhf_last)
def get_fock(mf, h1e=None, s1e=None, vhf=None, dm=None, cycle=-1, diis=None,
diis_start_cycle=None, level_shift_factor=None, damp_factor=None):
'''F = h^{core} + V^{HF}
Special treatment (damping, DIIS, or level shift) will be applied to the
Fock matrix if diis and cycle is specified (The two parameters are passed
to get_fock function during the SCF iteration)
Kwargs:
h1e : 2D ndarray
Core hamiltonian
s1e : 2D ndarray
Overlap matrix, for DIIS
vhf : 2D ndarray
HF potential matrix
dm : 2D ndarray
Density matrix, for DIIS
cycle : int
Then present SCF iteration step, for DIIS
diis : an object of :attr:`SCF.DIIS` class
DIIS object to hold intermediate Fock and error vectors
diis_start_cycle : int
The step to start DIIS. Default is 0.
level_shift_factor : float or int
Level shift (in AU) for virtual space. Default is 0.
'''
if h1e is None: h1e = mf.get_hcore()
if vhf is None: vhf = mf.get_veff(mf.mol, dm)
f = h1e + vhf
if cycle < 0 and diis is None: # Not inside the SCF iteration
return f
if diis_start_cycle is None:
diis_start_cycle = mf.diis_start_cycle
if level_shift_factor is None:
level_shift_factor = mf.level_shift
if damp_factor is None:
damp_factor = mf.damp
if s1e is None: s1e = mf.get_ovlp()
if dm is None: dm = mf.make_rdm1()
if 0 <= cycle < diis_start_cycle-1 and abs(damp_factor) > 1e-4:
f = damping(s1e, dm*.5, f, damp_factor)
if diis is not None and cycle >= diis_start_cycle:
f = diis.update(s1e, dm, f, mf, h1e, vhf)
if abs(level_shift_factor) > 1e-4:
f = level_shift(s1e, dm*.5, f, level_shift_factor)
return f
def get_occ(mf, mo_energy=None, mo_coeff=None):
'''Label the occupancies for each orbital
Kwargs:
mo_energy : 1D ndarray
Obital energies
mo_coeff : 2D ndarray
Obital coefficients
Examples:
>>> from pyscf import gto, scf
>>> mol = gto.M(atom='H 0 0 0; F 0 0 1.1')
>>> mf = scf.hf.SCF(mol)
>>> energy = numpy.array([-10., -1., 1, -2., 0, -3])
>>> mf.get_occ(energy)
array([2, 2, 0, 2, 2, 2])
'''
if mo_energy is None: mo_energy = mf.mo_energy
e_idx = numpy.argsort(mo_energy)
e_sort = mo_energy[e_idx]
nmo = mo_energy.size
mo_occ = numpy.zeros(nmo)
nocc = mf.mol.nelectron // 2
mo_occ[e_idx[:nocc]] = 2
if mf.verbose >= logger.INFO and nocc < nmo:
if e_sort[nocc-1]+1e-3 > e_sort[nocc]:
logger.warn(mf, 'HOMO %.15g == LUMO %.15g',
e_sort[nocc-1], e_sort[nocc])
else:
logger.info(mf, ' HOMO = %.15g LUMO = %.15g',
e_sort[nocc-1], e_sort[nocc])
if mf.verbose >= logger.DEBUG:
numpy.set_printoptions(threshold=nmo)
logger.debug(mf, ' mo_energy =\n%s', mo_energy)
numpy.set_printoptions(threshold=1000)
return mo_occ
def get_grad(mo_coeff, mo_occ, fock_ao):
'''RHF orbital gradients
Args:
mo_coeff : 2D ndarray
Obital coefficients
mo_occ : 1D ndarray
Orbital occupancy
fock_ao : 2D ndarray
Fock matrix in AO representation
Returns:
Gradients in MO representation. It's a num_occ*num_vir vector.
'''
occidx = mo_occ > 0
viridx = ~occidx
g = reduce(numpy.dot, (mo_coeff[:,viridx].conj().T, fock_ao,
mo_coeff[:,occidx])) * 2
return g.ravel()
def analyze(mf, verbose=logger.DEBUG, with_meta_lowdin=WITH_META_LOWDIN,
**kwargs):
'''Analyze the given SCF object: print orbital energies, occupancies;
print orbital coefficients; Mulliken population analysis; Diople moment.
'''
from pyscf.lo import orth
from pyscf.tools import dump_mat
mo_energy = mf.mo_energy
mo_occ = mf.mo_occ
mo_coeff = mf.mo_coeff
log = logger.new_logger(mf, verbose)
if log.verbose >= logger.NOTE:
mf.dump_scf_summary(log)
log.note('**** MO energy ****')
for i,c in enumerate(mo_occ):
log.note('MO #%-3d energy= %-18.15g occ= %g', i+MO_BASE,
mo_energy[i], c)
ovlp_ao = mf.get_ovlp()
if verbose >= logger.DEBUG:
label = mf.mol.ao_labels()
if with_meta_lowdin:
log.debug(' ** MO coefficients (expansion on meta-Lowdin AOs) **')
orth_coeff = orth.orth_ao(mf.mol, 'meta_lowdin', s=ovlp_ao)
c = reduce(numpy.dot, (orth_coeff.conj().T, ovlp_ao, mo_coeff))
else:
log.debug(' ** MO coefficients (expansion on AOs) **')
c = mo_coeff
dump_mat.dump_rec(mf.stdout, c, label, start=MO_BASE, **kwargs)
dm = mf.make_rdm1(mo_coeff, mo_occ)
if with_meta_lowdin:
return (mf.mulliken_meta(mf.mol, dm, s=ovlp_ao, verbose=log),
mf.dip_moment(mf.mol, dm, verbose=log))
else:
return (mf.mulliken_pop(mf.mol, dm, s=ovlp_ao, verbose=log),
mf.dip_moment(mf.mol, dm, verbose=log))