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