1#!/usr/bin/env python
2# Copyright 2014-2019 The PySCF Developers. All Rights Reserved.
3#
4# Licensed under the Apache License, Version 2.0 (the "License");
5# you may not use this file except in compliance with the License.
6# You may obtain a copy of the License at
7#
8#     http://www.apache.org/licenses/LICENSE-2.0
9#
10# Unless required by applicable law or agreed to in writing, software
11# distributed under the License is distributed on an "AS IS" BASIS,
12# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13# See the License for the specific language governing permissions and
14# limitations under the License.
15#
16# Author: Qiming Sun <osirpt.sun@gmail.com>
17#
18
19'''
20Non-relativistic generalized Hartree-Fock
21'''
22
23from functools import reduce
24import numpy
25import scipy.linalg
26from pyscf import lib
27from pyscf import gto
28from pyscf.lib import logger
29from pyscf.scf import hf
30from pyscf.scf import uhf
31from pyscf.scf import chkfile
32from pyscf import __config__
33
34WITH_META_LOWDIN = getattr(__config__, 'scf_analyze_with_meta_lowdin', True)
35PRE_ORTH_METHOD = getattr(__config__, 'scf_analyze_pre_orth_method', 'ANO')
36MO_BASE = getattr(__config__, 'MO_BASE', 1)
37
38
39def init_guess_by_chkfile(mol, chkfile_name, project=None):
40    '''Read SCF chkfile and make the density matrix for GHF initial guess.
41
42    Kwargs:
43        project : None or bool
44            Whether to project chkfile's orbitals to the new basis.  Note when
45            the geometry of the chkfile and the given molecule are very
46            different, this projection can produce very poor initial guess.
47            In PES scanning, it is recommended to swith off project.
48
49            If project is set to None, the projection is only applied when the
50            basis sets of the chkfile's molecule are different to the basis
51            sets of the given molecule (regardless whether the geometry of
52            the two molecules are different).  Note the basis sets are
53            considered to be different if the two molecules are derived from
54            the same molecule with different ordering of atoms.
55    '''
56    from pyscf.scf import addons
57    chk_mol, scf_rec = chkfile.load_scf(chkfile_name)
58    if project is None:
59        project = not gto.same_basis_set(chk_mol, mol)
60
61    # Check whether the two molecules are similar
62    if abs(mol.inertia_moment() - chk_mol.inertia_moment()).sum() > 0.5:
63        logger.warn(mol, "Large deviations found between the input "
64                    "molecule and the molecule from chkfile\n"
65                    "Initial guess density matrix may have large error.")
66
67    if project:
68        s = hf.get_ovlp(mol)
69
70    def fproj(mo):
71        if project:
72            mo = addons.project_mo_nr2nr(chk_mol, mo, mol)
73            norm = numpy.einsum('pi,pi->i', mo.conj(), s.dot(mo))
74            mo /= numpy.sqrt(norm)
75        return mo
76
77    nao = chk_mol.nao_nr()
78    mo = scf_rec['mo_coeff']
79    mo_occ = scf_rec['mo_occ']
80    if getattr(mo[0], 'ndim', None) == 1:  # RHF/GHF/DHF
81        if nao*2 == mo.shape[0]:  # GHF or DHF
82            if project:
83                raise NotImplementedError('Project initial guess from '
84                                          'different geometry')
85            else:
86                dm = hf.make_rdm1(mo, mo_occ)
87        else:  # RHF
88            mo_coeff = fproj(mo)
89            mo_occa = (mo_occ>1e-8).astype(numpy.double)
90            mo_occb = mo_occ - mo_occa
91            dma, dmb = uhf.make_rdm1([mo_coeff]*2, (mo_occa, mo_occb))
92            dm = scipy.linalg.block_diag(dma, dmb)
93    else: #UHF
94        if getattr(mo[0][0], 'ndim', None) == 2:  # KUHF
95            logger.warn(mol, 'k-point UHF results are found.  Density matrix '
96                        'at Gamma point is used for the molecular SCF initial guess')
97            mo = mo[0]
98        dma, dmb = uhf.make_rdm1([fproj(mo[0]),fproj(mo[1])], mo_occ)
99        dm = scipy.linalg.block_diag(dma, dmb)
100    return dm
101
102
103@lib.with_doc(hf.get_jk.__doc__)
104def get_jk(mol, dm, hermi=0,
105           with_j=True, with_k=True, jkbuild=hf.get_jk, omega=None):
106
107    dm = numpy.asarray(dm)
108    nso = dm.shape[-1]
109    nao = nso // 2
110    dms = dm.reshape(-1,nso,nso)
111    n_dm = dms.shape[0]
112
113    dmaa = dms[:,:nao,:nao]
114    dmab = dms[:,:nao,nao:]
115    dmbb = dms[:,nao:,nao:]
116    dms = numpy.vstack((dmaa, dmbb, dmab))
117    if dm.dtype == numpy.complex128:
118        dms = numpy.vstack((dms.real, dms.imag))
119
120    hermi = 0  # The off-diagonal blocks are not hermitian
121    j1, k1 = jkbuild(mol, dms, hermi, with_j, with_k, omega)
122    if with_j: j1 = j1.reshape(-1,n_dm,nao,nao)
123    if with_k: k1 = k1.reshape(-1,n_dm,nao,nao)
124
125    if dm.dtype == numpy.complex128:
126        if with_j: j1 = j1[:3] + j1[3:] * 1j
127        if with_k: k1 = k1[:3] + k1[3:] * 1j
128
129    vj = vk = None
130    if with_j:
131        vj = numpy.zeros((n_dm,nso,nso), dm.dtype)
132        vj[:,:nao,:nao] = vj[:,nao:,nao:] = j1[0] + j1[1]
133        vj = vj.reshape(dm.shape)
134
135    if with_k:
136        vk = numpy.zeros((n_dm,nso,nso), dm.dtype)
137        vk[:,:nao,:nao] = k1[0]
138        vk[:,nao:,nao:] = k1[1]
139        vk[:,:nao,nao:] = k1[2]
140        vk[:,nao:,:nao] = k1[2].transpose(0,2,1).conj()
141        vk = vk.reshape(dm.shape)
142
143    return vj, vk
144
145def get_occ(mf, mo_energy=None, mo_coeff=None):
146    if mo_energy is None: mo_energy = mf.mo_energy
147    e_idx = numpy.argsort(mo_energy)
148    e_sort = mo_energy[e_idx]
149    nmo = mo_energy.size
150    mo_occ = numpy.zeros(nmo)
151    nocc = mf.mol.nelectron
152    mo_occ[e_idx[:nocc]] = 1
153    if mf.verbose >= logger.INFO and nocc < nmo:
154        if e_sort[nocc-1]+1e-3 > e_sort[nocc]:
155            logger.warn(mf, 'HOMO %.15g == LUMO %.15g',
156                        e_sort[nocc-1], e_sort[nocc])
157        else:
158            logger.info(mf, '  HOMO = %.15g  LUMO = %.15g',
159                        e_sort[nocc-1], e_sort[nocc])
160
161    if mf.verbose >= logger.DEBUG:
162        numpy.set_printoptions(threshold=nmo)
163        logger.debug(mf, '  mo_energy =\n%s', mo_energy)
164        numpy.set_printoptions(threshold=1000)
165
166    if mo_coeff is not None and mf.verbose >= logger.DEBUG:
167        ss, s = mf.spin_square(mo_coeff[:,mo_occ>0], mf.get_ovlp())
168        logger.debug(mf, 'multiplicity <S^2> = %.8g  2S+1 = %.8g', ss, s)
169    return mo_occ
170
171# mo_a and mo_b are occupied orbitals
172def spin_square(mo, s=1):
173    r'''Spin of the GHF wavefunction
174
175    .. math::
176
177        S^2 = \frac{1}{2}(S_+ S_-  +  S_- S_+) + S_z^2
178
179    where :math:`S_+ = \sum_i S_{i+}` is effective for all beta occupied
180    orbitals; :math:`S_- = \sum_i S_{i-}` is effective for all alpha occupied
181    orbitals.
182
183    1. There are two possibilities for :math:`S_+ S_-`
184        1) same electron :math:`S_+ S_- = \sum_i s_{i+} s_{i-}`,
185
186        .. math::
187
188            \sum_i \langle UHF|s_{i+} s_{i-}|UHF\rangle
189             = \sum_{pq}\langle p|s_+s_-|q\rangle \gamma_{qp} = n_\alpha
190
191        2) different electrons :math:`S_+ S_- = \sum s_{i+} s_{j-},  (i\neq j)`.
192        There are in total :math:`n(n-1)` terms.  As a two-particle operator,
193
194        .. math::
195
196            \langle S_+ S_- \rangle
197            =\sum_{ij}(\langle i^\alpha|i^\beta\rangle \langle j^\beta|j^\alpha\rangle
198            - \langle i^\alpha|j^\beta\rangle \langle j^\beta|i^\alpha\rangle)
199
200    2. Similarly, for :math:`S_- S_+`
201        1) same electron
202
203        .. math::
204
205           \sum_i \langle s_{i-} s_{i+}\rangle = n_\beta
206
207        2) different electrons
208
209        .. math::
210
211            \langle S_- S_+ \rangle
212            =\sum_{ij}(\langle i^\beta|i^\alpha\rangle \langle j^\alpha|j^\beta\rangle
213            - \langle i^\beta|j^\alpha\rangle \langle j^\alpha|i^\beta\rangle)
214
215    3. For :math:`S_z^2`
216        1) same electron
217
218        .. math::
219
220            \langle s_z^2\rangle = \frac{1}{4}(n_\alpha + n_\beta)
221
222        2) different electrons
223
224        .. math::
225
226            &\sum_{ij}(\langle ij|s_{z1}s_{z2}|ij\rangle
227                      -\langle ij|s_{z1}s_{z2}|ji\rangle) \\
228            &=\frac{1}{4}\sum_{ij}(\langle i^\alpha|i^\alpha\rangle \langle j^\alpha|j^\alpha\rangle
229             - \langle i^\alpha|i^\alpha\rangle \langle j^\beta|j^\beta\rangle
230             - \langle i^\beta|i^\beta\rangle \langle j^\alpha|j^\alpha\rangle
231             + \langle i^\beta|i^\beta\rangle \langle j^\beta|j^\beta\rangle) \\
232            &-\frac{1}{4}\sum_{ij}(\langle i^\alpha|j^\alpha\rangle \langle j^\alpha|i^\alpha\rangle
233             - \langle i^\alpha|j^\alpha\rangle \langle j^\beta|i^\beta\rangle
234             - \langle i^\beta|j^\beta\rangle \langle j^\alpha|i^\alpha\rangle
235             + \langle i^\beta|j^\beta\rangle\langle j^\beta|i^\beta\rangle) \\
236            &=\frac{1}{4}\sum_{ij}|\langle i^\alpha|i^\alpha\rangle - \langle i^\beta|i^\beta\rangle|^2
237             -\frac{1}{4}\sum_{ij}|\langle i^\alpha|j^\alpha\rangle - \langle i^\beta|j^\beta\rangle|^2 \\
238            &=\frac{1}{4}(n_\alpha - n_\beta)^2
239             -\frac{1}{4}\sum_{ij}|\langle i^\alpha|j^\alpha\rangle - \langle i^\beta|j^\beta\rangle|^2
240
241    Args:
242        mo : a list of 2 ndarrays
243            Occupied alpha and occupied beta orbitals
244
245    Kwargs:
246        s : ndarray
247            AO overlap
248
249    Returns:
250        A list of two floats.  The first is the expectation value of S^2.
251        The second is the corresponding 2S+1
252
253    Examples:
254
255    >>> mol = gto.M(atom='O 0 0 0; H 0 0 1; H 0 1 0', basis='ccpvdz', charge=1, spin=1, verbose=0)
256    >>> mf = scf.UHF(mol)
257    >>> mf.kernel()
258    -75.623975516256706
259    >>> mo = (mf.mo_coeff[0][:,mf.mo_occ[0]>0], mf.mo_coeff[1][:,mf.mo_occ[1]>0])
260    >>> print('S^2 = %.7f, 2S+1 = %.7f' % spin_square(mo, mol.intor('int1e_ovlp_sph')))
261    S^2 = 0.7570150, 2S+1 = 2.0070027
262    '''
263    nao = mo.shape[0] // 2
264    if isinstance(s, numpy.ndarray):
265        assert(s.size == nao**2 or numpy.allclose(s[:nao,:nao], s[nao:,nao:]))
266        s = s[:nao,:nao]
267    mo_a = mo[:nao]
268    mo_b = mo[nao:]
269    saa = reduce(numpy.dot, (mo_a.conj().T, s, mo_a))
270    sbb = reduce(numpy.dot, (mo_b.conj().T, s, mo_b))
271    sab = reduce(numpy.dot, (mo_a.conj().T, s, mo_b))
272    sba = sab.conj().T
273    nocc_a = saa.trace()
274    nocc_b = sbb.trace()
275    ssxy = (nocc_a+nocc_b) * .5
276    ssxy+= sba.trace() * sab.trace() - numpy.einsum('ij,ji->', sba, sab)
277    ssz  = (nocc_a+nocc_b) * .25
278    ssz += (nocc_a-nocc_b)**2 * .25
279    tmp  = saa - sbb
280    ssz -= numpy.einsum('ij,ji', tmp, tmp) * .25
281    ss = (ssxy + ssz).real
282    s = numpy.sqrt(ss+.25) - .5
283    return ss, s*2+1
284
285def analyze(mf, verbose=logger.DEBUG, with_meta_lowdin=WITH_META_LOWDIN,
286            **kwargs):
287    '''Analyze the given SCF object:  print orbital energies, occupancies;
288    print orbital coefficients; Mulliken population analysis; Dipole moment
289    '''
290    log = logger.new_logger(mf, verbose)
291    mo_energy = mf.mo_energy
292    mo_occ = mf.mo_occ
293    mo_coeff = mf.mo_coeff
294
295    log.note('**** MO energy ****')
296    for i,c in enumerate(mo_occ):
297        log.note('MO #%-3d energy= %-18.15g occ= %g', i+MO_BASE, mo_energy[i], c)
298    ovlp_ao = mf.get_ovlp()
299    dm = mf.make_rdm1(mo_coeff, mo_occ)
300    dip = mf.dip_moment(mf.mol, dm, verbose=log)
301    if with_meta_lowdin:
302        pop_and_chg = mf.mulliken_meta(mf.mol, dm, s=ovlp_ao, verbose=log)
303    else:
304        pop_and_chg = mf.mulliken_pop(mf.mol, dm, s=ovlp_ao, verbose=log)
305    return pop_and_chg, dip
306
307def mulliken_pop(mol, dm, s=None, verbose=logger.DEBUG):
308    '''Mulliken population analysis
309    '''
310    nao = mol.nao_nr()
311    dma = dm[:nao,:nao]
312    dmb = dm[nao:,nao:]
313    if s is not None:
314        assert(s.size == nao**2 or numpy.allclose(s[:nao,:nao], s[nao:,nao:]))
315        s = s[:nao,:nao]
316    return uhf.mulliken_pop(mol, (dma,dmb), s, verbose)
317
318def mulliken_meta(mol, dm_ao, verbose=logger.DEBUG,
319                  pre_orth_method=PRE_ORTH_METHOD, s=None):
320    '''Mulliken population analysis, based on meta-Lowdin AOs.
321    '''
322    nao = mol.nao_nr()
323    dma = dm_ao[:nao,:nao]
324    dmb = dm_ao[nao:,nao:]
325    if s is not None:
326        assert(s.size == nao**2 or numpy.allclose(s[:nao,:nao], s[nao:,nao:]))
327        s = s[:nao,:nao]
328    return uhf.mulliken_meta(mol, (dma,dmb), verbose, pre_orth_method, s)
329
330def det_ovlp(mo1, mo2, occ1, occ2, ovlp):
331    r''' Calculate the overlap between two different determinants. It is the product
332    of single values of molecular orbital overlap matrix.
333
334    Return:
335        A list:
336            the product of single values: float
337            x_a: :math:`\mathbf{U} \mathbf{\Lambda}^{-1} \mathbf{V}^\dagger`
338            They are used to calculate asymmetric density matrix
339    '''
340
341    if numpy.sum(occ1) != numpy.sum(occ2):
342        raise RuntimeError('Electron numbers are not equal. Electronic coupling does not exist.')
343    s = reduce(numpy.dot, (mo1[:,occ1>0].T.conj(), ovlp, mo2[:,occ2>0]))
344    u, s, vt = numpy.linalg.svd(s)
345    x = numpy.dot(u/s, vt)
346    return numpy.prod(s), x
347
348def dip_moment(mol, dm, unit_symbol='Debye', verbose=logger.NOTE):
349    nao = mol.nao_nr()
350    dma = dm[:nao,:nao]
351    dmb = dm[nao:,nao:]
352    return hf.dip_moment(mol, dma+dmb, unit_symbol, verbose)
353
354canonicalize = hf.canonicalize
355
356def guess_orbspin(mo_coeff):
357    '''Guess the orbital spin (alpha 0, beta 1, unknown -1) based on the
358    orbital coefficients
359    '''
360    nao, nmo = mo_coeff.shape
361    mo_a = mo_coeff[:nao//2]
362    mo_b = mo_coeff[nao//2:]
363    # When all coefficients on alpha AOs are close to 0, it's a beta orbital
364    bidx = numpy.all(abs(mo_a) < 1e-14, axis=0)
365    aidx = numpy.all(abs(mo_b) < 1e-14, axis=0)
366    orbspin = numpy.empty(nmo, dtype=int)
367    orbspin[:] = -1
368    orbspin[aidx] = 0
369    orbspin[bidx] = 1
370    return orbspin
371
372class GHF(hf.SCF):
373    __doc__ = hf.SCF.__doc__ + '''
374
375    Attributes for GHF method
376        GHF orbital coefficients are 2D array.  Let nao be the number of spatial
377        AOs, mo_coeff[:nao] are the coefficients of AO with alpha spin;
378        mo_coeff[nao:nao*2] are the coefficients of AO with beta spin.
379    '''
380
381    def __init__(self, mol):
382        hf.SCF.__init__(self, mol)
383        self.with_soc = None
384        self._keys = self._keys.union(['with_soc'])
385
386    def get_hcore(self, mol=None):
387        if mol is None: mol = self.mol
388        hcore = hf.get_hcore(mol)
389        hcore = scipy.linalg.block_diag(hcore, hcore)
390
391        if self.with_soc and mol.has_ecp_soc():
392            # The ECP SOC contribution = <|1j * s * U_SOC|>
393            s = .5 * lib.PauliMatrices
394            ecpso = numpy.einsum('sxy,spq->xpyq', -1j * s, mol.intor('ECPso'))
395            hcore = hcore + ecpso.reshape(hcore.shape)
396        return hcore
397
398    def get_ovlp(self, mol=None):
399        if mol is None: mol = self.mol
400        s = hf.get_ovlp(mol)
401        return scipy.linalg.block_diag(s, s)
402
403    get_occ = get_occ
404
405    def get_grad(self, mo_coeff, mo_occ, fock=None):
406        if fock is None:
407            dm1 = self.make_rdm1(mo_coeff, mo_occ)
408            fock = self.get_hcore(self.mol) + self.get_veff(self.mol, dm1)
409        occidx = mo_occ > 0
410        viridx = ~occidx
411        g = reduce(numpy.dot, (mo_coeff[:,occidx].T.conj(), fock,
412                               mo_coeff[:,viridx]))
413        return g.conj().T.ravel()
414
415    @lib.with_doc(hf.SCF.init_guess_by_minao.__doc__)
416    def init_guess_by_minao(self, mol=None):
417        return _from_rhf_init_dm(hf.SCF.init_guess_by_minao(self, mol))
418
419    @lib.with_doc(hf.SCF.init_guess_by_atom.__doc__)
420    def init_guess_by_atom(self, mol=None):
421        return _from_rhf_init_dm(hf.SCF.init_guess_by_atom(self, mol))
422
423    @lib.with_doc(hf.SCF.init_guess_by_huckel.__doc__)
424    def init_guess_by_huckel(self, mol=None):
425        if mol is None: mol = self.mol
426        logger.info(self, 'Initial guess from on-the-fly Huckel, doi:10.1021/acs.jctc.8b01089.')
427        return _from_rhf_init_dm(hf.init_guess_by_huckel(mol))
428
429    @lib.with_doc(hf.SCF.init_guess_by_chkfile.__doc__)
430    def init_guess_by_chkfile(self, chkfile=None, project=None):
431        if chkfile is None: chkfile = self.chkfile
432        return init_guess_by_chkfile(self.mol, chkfile, project)
433
434    @lib.with_doc(hf.get_jk.__doc__)
435    def get_jk(self, mol=None, dm=None, hermi=0, with_j=True, with_k=True,
436               omega=None):
437        if mol is None: mol = self.mol
438        if dm is None: dm = self.make_rdm1()
439        nao = mol.nao
440        dm = numpy.asarray(dm)
441
442        def jkbuild(mol, dm, hermi, with_j, with_k, omega=None):
443            if (not omega and
444                (self._eri is not None or mol.incore_anyway or self._is_mem_enough())):
445                if self._eri is None:
446                    self._eri = mol.intor('int2e', aosym='s8')
447                return hf.dot_eri_dm(self._eri, dm, hermi, with_j, with_k)
448            else:
449                return hf.SCF.get_jk(self, mol, dm, hermi, with_j, with_k, omega)
450
451        if nao == dm.shape[-1]:
452            vj, vk = jkbuild(mol, dm, hermi, with_j, with_k, omega)
453        else:  # GHF density matrix, shape (2N,2N)
454            vj, vk = get_jk(mol, dm, hermi, with_j, with_k, jkbuild, omega)
455        return vj, vk
456
457    def get_veff(self, mol=None, dm=None, dm_last=0, vhf_last=0, hermi=1):
458        if mol is None: mol = self.mol
459        if dm is None: dm = self.make_rdm1()
460        if self._eri is not None or not self.direct_scf:
461            vj, vk = self.get_jk(mol, dm, hermi)
462            vhf = vj - vk
463        else:
464            ddm = numpy.asarray(dm) - numpy.asarray(dm_last)
465            vj, vk = self.get_jk(mol, ddm, hermi)
466            vhf = vj - vk + numpy.asarray(vhf_last)
467        return vhf
468
469    def analyze(self, verbose=None, **kwargs):
470        if verbose is None: verbose = self.verbose
471        return analyze(self, verbose, **kwargs)
472
473    def mulliken_pop(self, mol=None, dm=None, s=None, verbose=logger.DEBUG):
474        if mol is None: mol = self.mol
475        if dm is None: dm = self.make_rdm1()
476        if s is None: s = self.get_ovlp(mol)
477        return mulliken_pop(mol, dm, s=s, verbose=verbose)
478
479    def mulliken_meta(self, mol=None, dm=None, verbose=logger.DEBUG,
480                      pre_orth_method=PRE_ORTH_METHOD, s=None):
481        if mol is None: mol = self.mol
482        if dm is None: dm = self.make_rdm1()
483        if s is None: s = self.get_ovlp(mol)
484        return mulliken_meta(mol, dm, s=s, verbose=verbose,
485                             pre_orth_method=pre_orth_method)
486
487    @lib.with_doc(spin_square.__doc__)
488    def spin_square(self, mo_coeff=None, s=None):
489        if mo_coeff is None: mo_coeff = self.mo_coeff[:,self.mo_occ>0]
490        if s is None: s = self.get_ovlp()
491        return spin_square(mo_coeff, s)
492
493    @lib.with_doc(det_ovlp.__doc__)
494    def det_ovlp(self, mo1, mo2, occ1, occ2, ovlp=None):
495        if ovlp is None: ovlp = self.get_ovlp()
496        return det_ovlp(mo1, mo2, occ1, occ2, ovlp)
497
498    @lib.with_doc(dip_moment.__doc__)
499    def dip_moment(self, mol=None, dm=None, unit_symbol='Debye',
500                   verbose=logger.NOTE):
501        if mol is None: mol = self.mol
502        if dm is None: dm = self.make_rdm1()
503        return dip_moment(mol, dm, unit_symbol, verbose=verbose)
504
505    def _finalize(self):
506        ss, s = self.spin_square()
507
508        if self.converged:
509            logger.note(self, 'converged SCF energy = %.15g  '
510                        '<S^2> = %.8g  2S+1 = %.8g', self.e_tot, ss, s)
511        else:
512            logger.note(self, 'SCF not converged.')
513            logger.note(self, 'SCF energy = %.15g after %d cycles  '
514                        '<S^2> = %.8g  2S+1 = %.8g',
515                        self.e_tot, self.max_cycle, ss, s)
516        return self
517
518    def convert_from_(self, mf):
519        '''Create GHF object based on the RHF/UHF object'''
520        from pyscf.scf import addons
521        return addons.convert_to_ghf(mf, out=self)
522
523    def stability(self, internal=None, external=None, verbose=None):
524        from pyscf.scf.stability import ghf_stability
525        return ghf_stability(self, verbose)
526
527    def nuc_grad_method(self):
528        raise NotImplementedError
529
530def _from_rhf_init_dm(dm, breaksym=True):
531    dma = dm * .5
532    dm = scipy.linalg.block_diag(dma, dma)
533    if breaksym:
534        nao = dma.shape[0]
535        idx, idy = numpy.diag_indices(nao)
536        dm[idx+nao,idy] = dm[idx,idy+nao] = dma.diagonal() * .05
537    return dm
538
539del(PRE_ORTH_METHOD)
540
541
542if __name__ == '__main__':
543    mol = gto.Mole()
544    mol.verbose = 3
545    mol.atom = 'H 0 0 0; H 0 0 1; O .5 .6 .2'
546    mol.basis = 'ccpvdz'
547    mol.build()
548
549    mf = GHF(mol)
550    mf.kernel()
551
552    dm = mf.init_guess_by_1e(mol)
553    dm = dm + 0j
554    nao = mol.nao_nr()
555    numpy.random.seed(12)
556    dm[:nao,nao:] = numpy.random.random((nao,nao)) * .1j
557    dm[nao:,:nao] = dm[:nao,nao:].T.conj()
558    mf.kernel(dm)
559    mf.canonicalize(mf.mo_coeff, mf.mo_occ)
560    mf.analyze()
561    print(mf.spin_square())
562    print(mf.e_tot - -75.9125824421352)
563