1#!/usr/bin/env python 2# Copyright 2014-2020 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 19import sys 20 21from functools import reduce 22import numpy 23from pyscf import lib 24from pyscf.lib import logger 25from pyscf import gto 26from pyscf import scf 27from pyscf import ao2mo 28from pyscf import fci 29from pyscf.mcscf import addons 30from pyscf import __config__ 31 32WITH_META_LOWDIN = getattr(__config__, 'mcscf_analyze_with_meta_lowdin', True) 33LARGE_CI_TOL = getattr(__config__, 'mcscf_analyze_large_ci_tol', 0.1) 34PENALTY = getattr(__config__, 'mcscf_casci_CASCI_fix_spin_shift', 0.2) 35FRAC_OCC_THRESHOLD = 1e-6 36 37if sys.version_info < (3,): 38 RANGE_TYPE = list 39else: 40 RANGE_TYPE = range 41 42 43def h1e_for_cas(casci, mo_coeff=None, ncas=None, ncore=None): 44 '''CAS sapce one-electron hamiltonian 45 46 Args: 47 casci : a CASSCF/CASCI object or RHF object 48 49 Returns: 50 A tuple, the first is the effective one-electron hamiltonian defined in CAS space, 51 the second is the electronic energy from core. 52 ''' 53 if mo_coeff is None: mo_coeff = casci.mo_coeff 54 if ncas is None: ncas = casci.ncas 55 if ncore is None: ncore = casci.ncore 56 mo_core = mo_coeff[:,:ncore] 57 mo_cas = mo_coeff[:,ncore:ncore+ncas] 58 59 hcore = casci.get_hcore() 60 energy_core = casci.energy_nuc() 61 if mo_core.size == 0: 62 corevhf = 0 63 else: 64 core_dm = numpy.dot(mo_core, mo_core.conj().T) * 2 65 corevhf = casci.get_veff(casci.mol, core_dm) 66 energy_core += numpy.einsum('ij,ji', core_dm, hcore).real 67 energy_core += numpy.einsum('ij,ji', core_dm, corevhf).real * .5 68 h1eff = reduce(numpy.dot, (mo_cas.conj().T, hcore+corevhf, mo_cas)) 69 return h1eff, energy_core 70 71def analyze(casscf, mo_coeff=None, ci=None, verbose=None, 72 large_ci_tol=LARGE_CI_TOL, with_meta_lowdin=WITH_META_LOWDIN, 73 **kwargs): 74 from pyscf.lo import orth 75 from pyscf.tools import dump_mat 76 from pyscf.mcscf import addons 77 log = logger.new_logger(casscf, verbose) 78 79 if mo_coeff is None: mo_coeff = casscf.mo_coeff 80 if ci is None: ci = casscf.ci 81 nelecas = casscf.nelecas 82 ncas = casscf.ncas 83 ncore = casscf.ncore 84 nocc = ncore + ncas 85 mocore = mo_coeff[:,:ncore] 86 mocas = mo_coeff[:,ncore:nocc] 87 88 label = casscf.mol.ao_labels() 89 if (isinstance(ci, (list, tuple, RANGE_TYPE)) and 90 not isinstance(casscf.fcisolver, addons.StateAverageFCISolver)): 91 log.warn('Mulitple states found in CASCI/CASSCF solver. Density ' 92 'matrix of the first state is generated in .analyze() function.') 93 civec = ci[0] 94 else: 95 civec = ci 96 if getattr(casscf.fcisolver, 'make_rdm1s', None): 97 casdm1a, casdm1b = casscf.fcisolver.make_rdm1s(civec, ncas, nelecas) 98 casdm1 = casdm1a + casdm1b 99 dm1b = numpy.dot(mocore, mocore.conj().T) 100 dm1a = dm1b + reduce(numpy.dot, (mocas, casdm1a, mocas.conj().T)) 101 dm1b += reduce(numpy.dot, (mocas, casdm1b, mocas.conj().T)) 102 dm1 = dm1a + dm1b 103 if log.verbose >= logger.DEBUG2: 104 log.info('alpha density matrix (on AO)') 105 dump_mat.dump_tri(log.stdout, dm1a, label, **kwargs) 106 log.info('beta density matrix (on AO)') 107 dump_mat.dump_tri(log.stdout, dm1b, label, **kwargs) 108 else: 109 casdm1 = casscf.fcisolver.make_rdm1(civec, ncas, nelecas) 110 dm1a = (numpy.dot(mocore, mocore.conj().T) * 2 + 111 reduce(numpy.dot, (mocas, casdm1, mocas.conj().T))) 112 dm1b = None 113 dm1 = dm1a 114 115 if log.verbose >= logger.INFO: 116 ovlp_ao = casscf._scf.get_ovlp() 117 # note the last two args of ._eig for mc1step_symm 118 occ, ucas = casscf._eig(-casdm1, ncore, nocc) 119 log.info('Natural occ %s', str(-occ)) 120 mocas = numpy.dot(mocas, ucas) 121 if with_meta_lowdin: 122 log.info('Natural orbital (expansion on meta-Lowdin AOs) in CAS space') 123 orth_coeff = orth.orth_ao(casscf.mol, 'meta_lowdin', s=ovlp_ao) 124 mocas = reduce(numpy.dot, (orth_coeff.conj().T, ovlp_ao, mocas)) 125 else: 126 log.info('Natural orbital (expansion on AOs) in CAS space') 127 dump_mat.dump_rec(log.stdout, mocas, label, start=1, **kwargs) 128 if log.verbose >= logger.DEBUG2: 129 if not casscf.natorb: 130 log.debug2('NOTE: mc.mo_coeff in active space is different to ' 131 'the natural orbital coefficients printed in above.') 132 if with_meta_lowdin: 133 c = reduce(numpy.dot, (orth_coeff.conj().T, ovlp_ao, mo_coeff)) 134 log.debug2('MCSCF orbital (expansion on meta-Lowdin AOs)') 135 else: 136 c = mo_coeff 137 log.debug2('MCSCF orbital (expansion on AOs)') 138 dump_mat.dump_rec(log.stdout, c, label, start=1, **kwargs) 139 140 if casscf._scf.mo_coeff is not None: 141 addons.map2hf(casscf, casscf._scf.mo_coeff) 142 143 if (ci is not None and 144 (getattr(casscf.fcisolver, 'large_ci', None) or 145 getattr(casscf.fcisolver, 'states_large_ci', None))): 146 log.info('** Largest CI components **') 147 if isinstance(ci, (list, tuple, RANGE_TYPE)): 148 if hasattr(casscf.fcisolver, 'states_large_ci'): 149 # defined in state_average_mix_ mcscf object 150 res = casscf.fcisolver.states_large_ci(ci, casscf.ncas, casscf.nelecas, 151 large_ci_tol, return_strs=False) 152 else: 153 res = [casscf.fcisolver.large_ci(civec, casscf.ncas, casscf.nelecas, 154 large_ci_tol, return_strs=False) 155 for civec in ci] 156 for i, civec in enumerate(ci): 157 log.info(' [alpha occ-orbitals] [beta occ-orbitals] state %-3d CI coefficient', i) 158 for c,ia,ib in res[i]: 159 log.info(' %-20s %-30s %.12f', ia, ib, c) 160 else: 161 log.info(' [alpha occ-orbitals] [beta occ-orbitals] CI coefficient') 162 res = casscf.fcisolver.large_ci(ci, casscf.ncas, casscf.nelecas, 163 large_ci_tol, return_strs=False) 164 for c,ia,ib in res: 165 log.info(' %-20s %-30s %.12f', ia, ib, c) 166 167 if with_meta_lowdin: 168 casscf._scf.mulliken_meta(casscf.mol, dm1, s=ovlp_ao, verbose=log) 169 else: 170 casscf._scf.mulliken_pop(casscf.mol, dm1, s=ovlp_ao, verbose=log) 171 return dm1a, dm1b 172 173def get_fock(mc, mo_coeff=None, ci=None, eris=None, casdm1=None, verbose=None): 174 r''' 175 Effective one-electron Fock matrix in AO representation 176 f = \sum_{pq} E_{pq} F_{pq} 177 F_{pq} = h_{pq} + \sum_{rs} [(pq|rs)-(ps|rq)] DM_{sr} 178 179 Ref. 180 Theor. Chim. Acta., 91, 31 181 Chem. Phys. 48, 157 182 183 For state-average CASCI/CASSCF object, the effective fock matrix is based 184 on the state-average density matrix. To obtain Fock matrix of a specific 185 state in the state-average calculations, you can pass "casdm1" of the 186 specific state to this function. 187 188 Args: 189 mc: a CASSCF/CASCI object or RHF object 190 191 Kwargs: 192 mo_coeff (ndarray): orbitals that span the core, active and external 193 space. 194 ci (ndarray): CI coefficients (or objects to represent the CI 195 wavefunctions in DMRG/QMC-MCSCF calculations). 196 eris: Integrals for the MCSCF object. Input this object to reduce the 197 overhead of computing integrals. It can be generated by 198 :func:`mc.ao2mo` method. 199 casdm1 (ndarray): 1-particle density matrix in active space. Without 200 input casdm1, the density matrix is computed with the input ci 201 coefficients/object. If neither ci nor casdm1 were given, density 202 matrix is computed by :func:`mc.fcisolver.make_rdm1` method. For 203 state-average CASCI/CASCF calculation, this results in the 204 effective Fock matrix based on the state-average density matrix. 205 To obtain the effective Fock matrix for one particular state, you 206 can assign the density matrix of that state to the kwarg casdm1. 207 208 Returns: 209 Fock matrix 210 ''' 211 212 if ci is None: ci = mc.ci 213 if mo_coeff is None: mo_coeff = mc.mo_coeff 214 nmo = mo_coeff.shape[1] 215 ncore = mc.ncore 216 ncas = mc.ncas 217 nocc = ncore + ncas 218 nelecas = mc.nelecas 219 220 if casdm1 is None: 221 casdm1 = mc.fcisolver.make_rdm1(ci, ncas, nelecas) 222 if getattr(eris, 'ppaa', None) is not None: 223 vj = numpy.empty((nmo,nmo)) 224 vk = numpy.empty((nmo,nmo)) 225 for i in range(nmo): 226 vj[i] = numpy.einsum('ij,qij->q', casdm1, eris.ppaa[i]) 227 vk[i] = numpy.einsum('ij,iqj->q', casdm1, eris.papa[i]) 228 mo_inv = numpy.dot(mo_coeff.conj().T, mc._scf.get_ovlp()) 229 fock = (mc.get_hcore() + 230 reduce(numpy.dot, (mo_inv.conj().T, eris.vhf_c+vj-vk*.5, mo_inv))) 231 else: 232 dm_core = numpy.dot(mo_coeff[:,:ncore]*2, mo_coeff[:,:ncore].conj().T) 233 mocas = mo_coeff[:,ncore:nocc] 234 dm = dm_core + reduce(numpy.dot, (mocas, casdm1, mocas.conj().T)) 235 vj, vk = mc._scf.get_jk(mc.mol, dm) 236 fock = mc.get_hcore() + vj-vk*.5 237 return fock 238 239def cas_natorb(mc, mo_coeff=None, ci=None, eris=None, sort=False, 240 casdm1=None, verbose=None, with_meta_lowdin=WITH_META_LOWDIN): 241 '''Transform active orbitals to natrual orbitals, and update the CI wfn 242 accordingly 243 244 Args: 245 mc : a CASSCF/CASCI object or RHF object 246 247 Kwargs: 248 sort : bool 249 Sort natural orbitals wrt the occupancy. 250 251 Returns: 252 A tuple, the first item is natural orbitals, the second is updated CI 253 coefficients, the third is the natural occupancy associated to the 254 natural orbitals. 255 ''' 256 from pyscf.lo import orth 257 from pyscf.tools import dump_mat 258 from pyscf.tools.mo_mapping import mo_1to1map 259 if mo_coeff is None: mo_coeff = mc.mo_coeff 260 if ci is None: ci = mc.ci 261 log = logger.new_logger(mc, verbose) 262 ncore = mc.ncore 263 ncas = mc.ncas 264 nocc = ncore + ncas 265 nelecas = mc.nelecas 266 nmo = mo_coeff.shape[1] 267 if casdm1 is None: 268 casdm1 = mc.fcisolver.make_rdm1(ci, ncas, nelecas) 269 # orbital symmetry is reserved in this _eig call 270 cas_occ, ucas = mc._eig(-casdm1, ncore, nocc) 271 if sort: 272 casorb_idx = numpy.argsort(cas_occ.round(9), kind='mergesort') 273 cas_occ = cas_occ[casorb_idx] 274 ucas = ucas[:,casorb_idx] 275 276 cas_occ = -cas_occ 277 mo_occ = numpy.zeros(mo_coeff.shape[1]) 278 mo_occ[:ncore] = 2 279 mo_occ[ncore:nocc] = cas_occ 280 281 mo_coeff1 = mo_coeff.copy() 282 mo_coeff1[:,ncore:nocc] = numpy.dot(mo_coeff[:,ncore:nocc], ucas) 283 if getattr(mo_coeff, 'orbsym', None) is not None: 284 orbsym = numpy.copy(mo_coeff.orbsym) 285 if sort: 286 orbsym[ncore:nocc] = orbsym[ncore:nocc][casorb_idx] 287 mo_coeff1 = lib.tag_array(mo_coeff1, orbsym=orbsym) 288 else: 289 orbsym = numpy.zeros(nmo, dtype=int) 290 291 # When occupancies of active orbitals equal to 2 or 0, these orbitals 292 # need to be canonicalized along with inactive(core or virtual) orbitals 293 # using general Fock matrix. Because they are strongly coupled with 294 # inactive orbitals, the 0th order Hamiltonian of MRPT methods can be 295 # strongly affected. Numerical uncertainty may be found in the perturbed 296 # correlation energy. 297 # See issue https://github.com/pyscf/pyscf/issues/1041 298 occ2_idx = numpy.where(2 - cas_occ < FRAC_OCC_THRESHOLD)[0] 299 occ0_idx = numpy.where(cas_occ < FRAC_OCC_THRESHOLD)[0] 300 if occ2_idx.size > 0 or occ0_idx.size > 0: 301 fock_ao = mc.get_fock(mo_coeff, ci, eris, casdm1, verbose) 302 303 def _diag_subfock_(idx): 304 c = mo_coeff1[:,idx] 305 fock = reduce(numpy.dot, (c.conj().T, fock_ao, c)) 306 w, c = mc._eig(fock, None, None, orbsym[idx]) 307 mo_coeff1[:,idx] = mo_coeff1[:,idx].dot(c) 308 309 if occ2_idx.size > 0: 310 log.warn('Active orbitals %s (occs = %s) are canonicalized with core orbitals', 311 occ2_idx, cas_occ[occ2_idx]) 312 full_occ2_idx = numpy.append(numpy.arange(ncore), ncore + occ2_idx) 313 _diag_subfock_(full_occ2_idx) 314 if occ0_idx.size > 0: 315 log.warn('Active orbitals %s (occs = %s) are canonicalized with external orbitals', 316 occ0_idx, cas_occ[occ0_idx]) 317 full_occ0_idx = numpy.append(ncore + occ0_idx, numpy.arange(nocc, nmo)) 318 _diag_subfock_(full_occ0_idx) 319 320 # Rotate CI according to the unitary coefficients ucas if applicable 321 fcivec = None 322 if getattr(mc.fcisolver, 'transform_ci_for_orbital_rotation', None): 323 if isinstance(ci, numpy.ndarray): 324 fcivec = mc.fcisolver.transform_ci_for_orbital_rotation(ci, ncas, nelecas, ucas) 325 elif (isinstance(ci, (list, tuple)) and 326 all(isinstance(x[0], numpy.ndarray) for x in ci)): 327 fcivec = [mc.fcisolver.transform_ci_for_orbital_rotation(x, ncas, nelecas, ucas) 328 for x in ci] 329 elif getattr(mc.fcisolver, 'states_transform_ci_for_orbital_rotation', None): 330 fcivec = mc.fcisolver.states_transform_ci_for_orbital_rotation(ci, ncas, nelecas, ucas) 331 332 # Rerun fcisolver to get wavefunction if it cannot be transformed from 333 # existed one. 334 if fcivec is None: 335 log.info('FCI vector not available, call CASCI to update wavefunction') 336 mocas = mo_coeff1[:,ncore:nocc] 337 hcore = mc.get_hcore() 338 dm_core = numpy.dot(mo_coeff1[:,:ncore]*2, mo_coeff1[:,:ncore].conj().T) 339 ecore = mc.energy_nuc() 340 ecore+= numpy.einsum('ij,ji', hcore, dm_core) 341 h1eff = reduce(numpy.dot, (mocas.conj().T, hcore, mocas)) 342 if getattr(eris, 'ppaa', None) is not None: 343 ecore += eris.vhf_c[:ncore,:ncore].trace() 344 h1eff += reduce(numpy.dot, (ucas.conj().T, eris.vhf_c[ncore:nocc,ncore:nocc], ucas)) 345 aaaa = ao2mo.restore(4, eris.ppaa[ncore:nocc,ncore:nocc,:,:], ncas) 346 aaaa = ao2mo.incore.full(aaaa, ucas) 347 else: 348 if getattr(mc, 'with_df', None): 349 aaaa = mc.with_df.ao2mo(mocas) 350 else: 351 aaaa = ao2mo.kernel(mc.mol, mocas) 352 corevhf = mc.get_veff(mc.mol, dm_core) 353 ecore += numpy.einsum('ij,ji', dm_core, corevhf) * .5 354 h1eff += reduce(numpy.dot, (mocas.conj().T, corevhf, mocas)) 355 356 357 # See label_symmetry_ function in casci_symm.py which initialize the 358 # orbital symmetry information in fcisolver. This orbital symmetry 359 # labels should be reordered to match the sorted active space orbitals. 360 if sort and getattr(mo_coeff1, 'orbsym', None) is not None: 361 mc.fcisolver.orbsym = mo_coeff1.orbsym[ncore:nocc] 362 363 max_memory = max(400, mc.max_memory-lib.current_memory()[0]) 364 e, fcivec = mc.fcisolver.kernel(h1eff, aaaa, ncas, nelecas, ecore=ecore, 365 max_memory=max_memory, verbose=log) 366 log.debug('In Natural orbital, CASCI energy = %s', e) 367 368 if log.verbose >= logger.INFO: 369 ovlp_ao = mc._scf.get_ovlp() 370 # where_natorb gives the new locations of the natural orbitals 371 where_natorb = mo_1to1map(ucas) 372 log.debug('where_natorb %s', str(where_natorb)) 373 log.info('Natural occ %s', str(cas_occ)) 374 if with_meta_lowdin: 375 log.info('Natural orbital (expansion on meta-Lowdin AOs) in CAS space') 376 label = mc.mol.ao_labels() 377 orth_coeff = orth.orth_ao(mc.mol, 'meta_lowdin', s=ovlp_ao) 378 mo_cas = reduce(numpy.dot, (orth_coeff.conj().T, ovlp_ao, mo_coeff1[:,ncore:nocc])) 379 else: 380 log.info('Natural orbital (expansion on AOs) in CAS space') 381 label = mc.mol.ao_labels() 382 mo_cas = mo_coeff1[:,ncore:nocc] 383 dump_mat.dump_rec(log.stdout, mo_cas, label, start=1) 384 385 if mc._scf.mo_coeff is not None: 386 s = reduce(numpy.dot, (mo_coeff1[:,ncore:nocc].conj().T, 387 mc._scf.get_ovlp(), mc._scf.mo_coeff)) 388 idx = numpy.argwhere(abs(s)>.4) 389 for i,j in idx: 390 log.info('<CAS-nat-orb|mo-hf> %d %d %12.8f', 391 ncore+i+1, j+1, s[i,j]) 392 return mo_coeff1, fcivec, mo_occ 393 394def canonicalize(mc, mo_coeff=None, ci=None, eris=None, sort=False, 395 cas_natorb=False, casdm1=None, verbose=logger.NOTE, 396 with_meta_lowdin=WITH_META_LOWDIN): 397 '''Canonicalized CASCI/CASSCF orbitals of effecitive Fock matrix and 398 update CI coefficients accordingly. 399 400 Effective Fock matrix is built with one-particle density matrix (see 401 also :func:`mcscf.casci.get_fock`). For state-average CASCI/CASSCF object, 402 the canonicalized orbitals are based on the state-average density matrix. 403 To obtain canonicalized orbitals for an individual state, you need to pass 404 "casdm1" of the specific state to this function. 405 406 Args: 407 mc: a CASSCF/CASCI object or RHF object 408 409 Kwargs: 410 mo_coeff (ndarray): orbitals that span the core, active and external 411 space. 412 ci (ndarray): CI coefficients (or objects to represent the CI 413 wavefunctions in DMRG/QMC-MCSCF calculations). 414 eris: Integrals for the MCSCF object. Input this object to reduce the 415 overhead of computing integrals. It can be generated by 416 :func:`mc.ao2mo` method. 417 sort (bool): Whether the canonicalized orbitals are sorted based on 418 the orbital energy (diagonal part of the effective Fock matrix) 419 within each subspace (core, active, external). If point group 420 symmetry is not available in the system, orbitals are always 421 sorted. When point group symmetry is available, sort=False will 422 preserve the symmetry label of input orbitals and only sort the 423 orbitals in each symmetry sector. sort=True will reorder all 424 orbitals over all symmetry sectors in each subspace and the 425 symmetry labels may be changed. 426 cas_natorb (bool): Whether to transform active orbitals to natual 427 orbitals. If enabled, the output orbitals in active space are 428 transformed to natural orbitals and CI coefficients are updated 429 accordingly. 430 casdm1 (ndarray): 1-particle density matrix in active space. This 431 density matrix is used to build effective fock matrix. Without 432 input casdm1, the density matrix is computed with the input ci 433 coefficients/object. If neither ci nor casdm1 were given, density 434 matrix is computed by :func:`mc.fcisolver.make_rdm1` method. For 435 state-average CASCI/CASCF calculation, this results in a set of 436 canonicalized orbitals of state-average effective Fock matrix. 437 To canonicalize the orbitals for one particular state, you can 438 assign the density matrix of that state to the kwarg casdm1. 439 440 Returns: 441 A tuple, (natural orbitals, CI coefficients, orbital energies) 442 The orbital energies are the diagonal terms of effective Fock matrix. 443 ''' 444 from pyscf.mcscf import addons 445 log = logger.new_logger(mc, verbose) 446 447 if mo_coeff is None: mo_coeff = mc.mo_coeff 448 if ci is None: ci = mc.ci 449 if casdm1 is None: 450 if (isinstance(ci, (list, tuple, RANGE_TYPE)) and 451 not isinstance(mc.fcisolver, addons.StateAverageFCISolver)): 452 log.warn('Mulitple states found in CASCI solver. First state is ' 453 'used to compute the natural orbitals in active space.') 454 casdm1 = mc.fcisolver.make_rdm1(ci[0], mc.ncas, mc.nelecas) 455 else: 456 casdm1 = mc.fcisolver.make_rdm1(ci, mc.ncas, mc.nelecas) 457 458 ncore = mc.ncore 459 nocc = ncore + mc.ncas 460 nmo = mo_coeff.shape[1] 461 fock_ao = mc.get_fock(mo_coeff, ci, eris, casdm1, verbose) 462 463 if cas_natorb: 464 mo_coeff1, ci, mc.mo_occ = mc.cas_natorb(mo_coeff, ci, eris, sort, casdm1, 465 verbose, with_meta_lowdin) 466 else: 467 # Keep the active space unchanged by default. The rotation in active space 468 # may cause problem for external CI solver eg DMRG. 469 mo_coeff1 = mo_coeff.copy() 470 log.info('Density matrix diagonal elements %s', casdm1.diagonal()) 471 472 mo_energy = numpy.einsum('pi,pi->i', mo_coeff1.conj(), fock_ao.dot(mo_coeff1)) 473 474 if getattr(mo_coeff, 'orbsym', None) is not None: 475 orbsym = mo_coeff.orbsym 476 else: 477 orbsym = numpy.zeros(nmo, dtype=int) 478 479 def _diag_subfock_(idx): 480 if idx.size > 1: 481 c = mo_coeff1[:,idx] 482 fock = reduce(numpy.dot, (c.conj().T, fock_ao, c)) 483 # note the last argument orbysm is needed by mc1step_symm._eig 484 w, c = mc._eig(fock, None, None, orbsym[idx]) 485 486 if sort: 487 sub_order = numpy.argsort(w.round(9), kind='mergesort') 488 w = w[sub_order] 489 c = c[:,sub_order] 490 orbsym[idx] = orbsym[idx][sub_order] 491 492 mo_coeff1[:,idx] = mo_coeff1[:,idx].dot(c) 493 mo_energy[idx] = w 494 495 mask = numpy.ones(nmo, dtype=bool) 496 frozen = getattr(mc, 'frozen', None) 497 if frozen is not None: 498 if isinstance(frozen, (int, numpy.integer)): 499 mask[:frozen] = False 500 else: 501 mask[frozen] = False 502 core_idx = numpy.where(mask[:ncore])[0] 503 vir_idx = numpy.where(mask[nocc:])[0] + nocc 504 _diag_subfock_(core_idx) 505 _diag_subfock_(vir_idx) 506 507 # orbsym is required only for symmetry-adapted methods. Here to use 508 # mo_coeff.orbsym to test if a symmetry-adapted calculation. 509 if getattr(mo_coeff, 'orbsym', None) is not None: 510 mo_coeff1 = lib.tag_array(mo_coeff1, orbsym=orbsym) 511 512 if log.verbose >= logger.DEBUG: 513 for i in range(nmo): 514 log.debug('i = %d <i|F|i> = %12.8f', i+1, mo_energy[i]) 515# still return ci coefficients, in case the canonicalization funciton changed 516# cas orbitals, the ci coefficients should also be updated. 517 return mo_coeff1, ci, mo_energy 518 519 520def kernel(casci, mo_coeff=None, ci0=None, verbose=logger.NOTE): 521 '''CASCI solver 522 ''' 523 if mo_coeff is None: mo_coeff = casci.mo_coeff 524 log = logger.new_logger(casci, verbose) 525 t0 = (logger.process_clock(), logger.perf_counter()) 526 log.debug('Start CASCI') 527 528 ncas = casci.ncas 529 nelecas = casci.nelecas 530 531 # 2e 532 eri_cas = casci.get_h2eff(mo_coeff) 533 t1 = log.timer('integral transformation to CAS space', *t0) 534 535 # 1e 536 h1eff, energy_core = casci.get_h1eff(mo_coeff) 537 log.debug('core energy = %.15g', energy_core) 538 t1 = log.timer('effective h1e in CAS space', *t1) 539 540 if h1eff.shape[0] != ncas: 541 raise RuntimeError('Active space size error. nmo=%d ncore=%d ncas=%d' % 542 (mo_coeff.shape[1], casci.ncore, ncas)) 543 544 # FCI 545 max_memory = max(400, casci.max_memory-lib.current_memory()[0]) 546 e_tot, fcivec = casci.fcisolver.kernel(h1eff, eri_cas, ncas, nelecas, 547 ci0=ci0, verbose=log, 548 max_memory=max_memory, 549 ecore=energy_core) 550 551 t1 = log.timer('FCI solver', *t1) 552 e_cas = e_tot - energy_core 553 return e_tot, e_cas, fcivec 554 555 556def as_scanner(mc): 557 '''Generating a scanner for CASCI PES. 558 559 The returned solver is a function. This function requires one argument 560 "mol" as input and returns total CASCI energy. 561 562 The solver will automatically use the results of last calculation as the 563 initial guess of the new calculation. All parameters of MCSCF object 564 are automatically applied in the solver. 565 566 Note scanner has side effects. It may change many underlying objects 567 (_scf, with_df, with_x2c, ...) during calculation. 568 569 Examples: 570 571 >>> from pyscf import gto, scf, mcscf 572 >>> mf = scf.RHF(gto.Mole().set(verbose=0)) 573 >>> mc_scanner = mcscf.CASCI(mf, 4, 4).as_scanner() 574 >>> mc_scanner(gto.M(atom='N 0 0 0; N 0 0 1.1')) 575 >>> mc_scanner(gto.M(atom='N 0 0 0; N 0 0 1.5')) 576 ''' 577 if isinstance(mc, lib.SinglePointScanner): 578 return mc 579 580 logger.info(mc, 'Create scanner for %s', mc.__class__) 581 582 class CASCI_Scanner(mc.__class__, lib.SinglePointScanner): 583 def __init__(self, mc): 584 self.__dict__.update(mc.__dict__) 585 self._scf = mc._scf.as_scanner() 586 587 def __call__(self, mol_or_geom, mo_coeff=None, ci0=None): 588 if isinstance(mol_or_geom, gto.Mole): 589 mol = mol_or_geom 590 else: 591 mol = self.mol.set_geom_(mol_or_geom, inplace=False) 592 593 # These properties can be updated when calling mf_scanner(mol) if 594 # they are shared with mc._scf. In certain scenario the properties 595 # may be created for mc separately, e.g. when mcscf.approx_hessian is 596 # called. For safety, the code below explicitly resets these 597 # properties. 598 for key in ('with_df', 'with_x2c', 'with_solvent', 'with_dftd3'): 599 sub_mod = getattr(self, key, None) 600 if sub_mod: 601 sub_mod.reset(mol) 602 603 if mo_coeff is None: 604 mf_scanner = self._scf 605 mf_scanner(mol) 606 mo_coeff = mf_scanner.mo_coeff 607 if ci0 is None: 608 ci0 = self.ci 609 self.mol = mol 610 e_tot = self.kernel(mo_coeff, ci0)[0] 611 return e_tot 612 return CASCI_Scanner(mc) 613 614 615class CASCI(lib.StreamObject): 616 '''CASCI 617 618 Args: 619 mf_or_mol : SCF object or Mole object 620 SCF or Mole to define the problem size. 621 ncas : int 622 Number of active orbitals. 623 nelecas : int or a pair of int 624 Number of electrons in active space. 625 626 Kwargs: 627 ncore : int 628 Number of doubly occupied core orbitals. If not presented, this 629 parameter can be automatically determined. 630 631 Attributes: 632 verbose : int 633 Print level. Default value equals to :class:`Mole.verbose`. 634 max_memory : float or int 635 Allowed memory in MB. Default value equals to :class:`Mole.max_memory`. 636 ncas : int 637 Active space size. 638 nelecas : tuple of int 639 Active (nelec_alpha, nelec_beta) 640 ncore : int or tuple of int 641 Core electron number. In UHF-CASSCF, it's a tuple to indicate the different core eletron numbers. 642 natorb : bool 643 Whether to transform natural orbitals in active space. 644 Note: when CASCI/CASSCF are combined with DMRG solver or selected 645 CI solver, enabling this parameter may slightly change the total energy. 646 False by default. 647 canonicalization : bool 648 Whether to canonicalize orbitals in core and external space 649 against the general Fock matrix. 650 The orbitals in active space are NOT transformed by default. To 651 get the natural orbitals in active space, the attribute .natorb 652 needs to be enabled. 653 True by default. 654 sorting_mo_energy : bool 655 Whether to sort the orbitals based on the diagonal elements of the 656 general Fock matrix. Default is False. 657 fcisolver : an instance of :class:`FCISolver` 658 The pyscf.fci module provides several FCISolver for different scenario. Generally, 659 fci.direct_spin1.FCISolver can be used for all RHF-CASSCF. However, a proper FCISolver 660 can provide better performance and better numerical stability. One can either use 661 :func:`fci.solver` function to pick the FCISolver by the program or manually assigen 662 the FCISolver to this attribute, e.g. 663 664 >>> from pyscf import fci 665 >>> mc = mcscf.CASSCF(mf, 4, 4) 666 >>> mc.fcisolver = fci.solver(mol, singlet=True) 667 >>> mc.fcisolver = fci.direct_spin1.FCISolver(mol) 668 669 You can control FCISolver by setting e.g.:: 670 671 >>> mc.fcisolver.max_cycle = 30 672 >>> mc.fcisolver.conv_tol = 1e-7 673 674 For more details of the parameter for FCISolver, See :mod:`fci`. 675 676 Saved results 677 678 e_tot : float 679 Total MCSCF energy (electronic energy plus nuclear repulsion) 680 e_cas : float 681 CAS space FCI energy 682 ci : ndarray 683 CAS space FCI coefficients 684 mo_coeff : ndarray 685 When canonicalization is specified, the orbitals are canonical 686 orbitals which make the general Fock matrix (Fock operator on top 687 of MCSCF 1-particle density matrix) diagonalized within each 688 subspace (core, active, external). If natorb (natural orbitals in 689 active space) is specified, the active segment of the mo_coeff is 690 natural orbitls. 691 mo_energy : ndarray 692 Diagonal elements of general Fock matrix (in mo_coeff 693 representation). 694 mo_occ : ndarray 695 Occupation numbers of natural orbitals if natorb is specified. 696 697 Examples: 698 699 >>> from pyscf import gto, scf, mcscf 700 >>> mol = gto.M(atom='N 0 0 0; N 0 0 1', basis='ccpvdz', verbose=0) 701 >>> mf = scf.RHF(mol) 702 >>> mf.scf() 703 >>> mc = mcscf.CASCI(mf, 6, 6) 704 >>> mc.kernel()[0] 705 -108.980200816243354 706 ''' 707 708 natorb = getattr(__config__, 'mcscf_casci_CASCI_natorb', False) 709 canonicalization = getattr(__config__, 'mcscf_casci_CASCI_canonicalization', True) 710 sorting_mo_energy = getattr(__config__, 'mcscf_casci_CASCI_sorting_mo_energy', False) 711 712 def __init__(self, mf_or_mol, ncas, nelecas, ncore=None): 713 if isinstance(mf_or_mol, gto.Mole): 714 mf = scf.RHF(mf_or_mol) 715 else: 716 mf = mf_or_mol 717 718 mol = mf.mol 719 self.mol = mol 720 self._scf = mf 721 self.verbose = mol.verbose 722 self.stdout = mol.stdout 723 self.max_memory = mf.max_memory 724 self.ncas = ncas 725 if isinstance(nelecas, (int, numpy.integer)): 726 nelecb = (nelecas-mol.spin)//2 727 neleca = nelecas - nelecb 728 self.nelecas = (neleca, nelecb) 729 else: 730 self.nelecas = (nelecas[0],nelecas[1]) 731 self.ncore = ncore 732 singlet = (getattr(__config__, 'mcscf_casci_CASCI_fcisolver_direct_spin0', False) 733 and self.nelecas[0] == self.nelecas[1]) # leads to direct_spin1 734 self.fcisolver = fci.solver(mol, singlet, symm=False) 735# CI solver parameters are set in fcisolver object 736 self.fcisolver.lindep = getattr(__config__, 737 'mcscf_casci_CASCI_fcisolver_lindep', 1e-10) 738 self.fcisolver.max_cycle = getattr(__config__, 739 'mcscf_casci_CASCI_fcisolver_max_cycle', 200) 740 self.fcisolver.conv_tol = getattr(__config__, 741 'mcscf_casci_CASCI_fcisolver_conv_tol', 1e-8) 742 743################################################## 744# don't modify the following attributes, they are not input options 745 self.e_tot = 0 746 self.e_cas = None 747 self.ci = None 748 self.mo_coeff = mf.mo_coeff 749 self.mo_energy = mf.mo_energy 750 self.mo_occ = None 751 self.converged = False 752 753 keys = set(('natorb', 'canonicalization', 'sorting_mo_energy')) 754 self._keys = set(self.__dict__.keys()).union(keys) 755 756 @property 757 def ncore(self): 758 if self._ncore is None: 759 ncorelec = self.mol.nelectron - sum(self.nelecas) 760 assert ncorelec % 2 == 0 761 assert ncorelec >= 0 762 return ncorelec // 2 763 else: 764 return self._ncore 765 @ncore.setter 766 def ncore(self, x): 767 assert x is None or isinstance(x, (int, numpy.integer)) 768 assert x is None or x >= 0 769 self._ncore = x 770 771 def dump_flags(self, verbose=None): 772 log = logger.new_logger(self, verbose) 773 log.info('') 774 log.info('******** CASCI flags ********') 775 ncore = self.ncore 776 ncas = self.ncas 777 nvir = self.mo_coeff.shape[1] - ncore - ncas 778 log.info('CAS (%de+%de, %do), ncore = %d, nvir = %d', 779 self.nelecas[0], self.nelecas[1], ncas, ncore, nvir) 780 log.info('natorb = %s', self.natorb) 781 log.info('canonicalization = %s', self.canonicalization) 782 log.info('sorting_mo_energy = %s', self.sorting_mo_energy) 783 log.info('max_memory %d (MB)', self.max_memory) 784 if getattr(self.fcisolver, 'dump_flags', None): 785 self.fcisolver.dump_flags(log.verbose) 786 if self.mo_coeff is None: 787 log.error('Orbitals for CASCI are not specified. The relevant SCF ' 788 'object may not be initialized.') 789 790 if (getattr(self._scf, 'with_solvent', None) and 791 not getattr(self, 'with_solvent', None)): 792 log.warn('''Solvent model %s was found at SCF level but not applied to the CASCI object. 793The SCF solvent model will not be applied to the current CASCI calculation. 794To enable the solvent model for CASCI, the following code needs to be called 795 from pyscf import solvent 796 mc = mcscf.CASCI(...) 797 mc = solvent.ddCOSMO(mc) 798''', 799 self._scf.with_solvent.__class__) 800 return self 801 802 def check_sanity(self): 803 assert self.ncas > 0 804 ncore = self.ncore 805 nvir = self.mo_coeff.shape[1] - ncore - self.ncas 806 assert ncore >= 0 807 assert nvir >= 0 808 assert ncore * 2 + sum(self.nelecas) == self.mol.nelectron 809 assert 0 <= self.nelecas[0] <= self.ncas 810 assert 0 <= self.nelecas[1] <= self.ncas 811 return self 812 813 def reset(self, mol=None): 814 if mol is not None: 815 self.mol = mol 816 self.fcisolver.mol = mol 817 self._scf.reset(mol) 818 return self 819 820 def energy_nuc(self): 821 return self._scf.energy_nuc() 822 823 def get_hcore(self, mol=None): 824 return self._scf.get_hcore(mol) 825 826 @lib.with_doc(scf.hf.get_jk.__doc__) 827 def get_jk(self, mol, dm, hermi=1, with_j=True, with_k=True, omega=None): 828 return self._scf.get_jk(mol, dm, hermi, 829 with_j=with_j, with_k=with_k, omega=omega) 830 831 @lib.with_doc(scf.hf.get_veff.__doc__) 832 def get_veff(self, mol=None, dm=None, hermi=1): 833 if mol is None: mol = self.mol 834 if dm is None: 835 mocore = self.mo_coeff[:,:self.ncore] 836 dm = numpy.dot(mocore, mocore.conj().T) * 2 837# don't call self._scf.get_veff because _scf might be DFT object 838 vj, vk = self.get_jk(mol, dm, hermi) 839 return vj - vk * .5 840 841 def _eig(self, h, *args): 842 return scf.hf.eig(h, None) 843 844 def get_h2cas(self, mo_coeff=None): 845 '''Compute the active space two-particle Hamiltonian. 846 847 Note It is different to get_h2eff when df.approx_hessian is applied, 848 in which get_h2eff function returns the DF integrals while get_h2cas 849 returns the regular 2-electron integrals. 850 ''' 851 return self.ao2mo(mo_coeff) 852 853 def get_h2eff(self, mo_coeff=None): 854 '''Compute the active space two-particle Hamiltonian. 855 856 Note It is different to get_h2cas when df.approx_hessian is applied. 857 in which get_h2eff function returns the DF integrals while get_h2cas 858 returns the regular 2-electron integrals. 859 ''' 860 return self.ao2mo(mo_coeff) 861 862 def ao2mo(self, mo_coeff=None): 863 '''Compute the active space two-particle Hamiltonian. 864 ''' 865 ncore = self.ncore 866 ncas = self.ncas 867 nocc = ncore + ncas 868 if mo_coeff is None: 869 ncore = self.ncore 870 mo_coeff = self.mo_coeff[:,ncore:nocc] 871 elif mo_coeff.shape[1] != ncas: 872 mo_coeff = mo_coeff[:,ncore:nocc] 873 874 if self._scf._eri is not None: 875 eri = ao2mo.full(self._scf._eri, mo_coeff, 876 max_memory=self.max_memory) 877 else: 878 eri = ao2mo.full(self.mol, mo_coeff, verbose=self.verbose, 879 max_memory=self.max_memory) 880 return eri 881 882 get_h1cas = h1e_for_cas = h1e_for_cas 883 884 def get_h1eff(self, mo_coeff=None, ncas=None, ncore=None): 885 return self.h1e_for_cas(mo_coeff, ncas, ncore) 886 get_h1eff.__doc__ = h1e_for_cas.__doc__ 887 888 def casci(self, mo_coeff=None, ci0=None, verbose=None): 889 return self.kernel(mo_coeff, ci0, verbose) 890 def kernel(self, mo_coeff=None, ci0=None, verbose=None): 891 ''' 892 Returns: 893 Five elements, they are 894 total energy, 895 active space CI energy, 896 the active space FCI wavefunction coefficients or DMRG wavefunction ID, 897 the MCSCF canonical orbital coefficients, 898 the MCSCF canonical orbital coefficients. 899 900 They are attributes of mcscf object, which can be accessed by 901 .e_tot, .e_cas, .ci, .mo_coeff, .mo_energy 902 ''' 903 if mo_coeff is None: 904 mo_coeff = self.mo_coeff 905 else: 906 self.mo_coeff = mo_coeff 907 if ci0 is None: 908 ci0 = self.ci 909 log = logger.new_logger(self, verbose) 910 911 self.check_sanity() 912 self.dump_flags(log) 913 914 self.e_tot, self.e_cas, self.ci = \ 915 kernel(self, mo_coeff, ci0=ci0, verbose=log) 916 917 if self.canonicalization: 918 self.canonicalize_(mo_coeff, self.ci, 919 sort=self.sorting_mo_energy, 920 cas_natorb=self.natorb, verbose=log) 921 elif self.natorb: 922 # FIXME (pyscf-2.0): Whether to transform natural orbitals in 923 # active space when this flag is enabled? 924 log.warn('The attribute .natorb of mcscf object affects only the ' 925 'orbital canonicalization.\n' 926 'If you would like to get natural orbitals in active space ' 927 'without touching core and external orbitals, an explicit ' 928 'call to mc.cas_natorb_() is required') 929 930 if getattr(self.fcisolver, 'converged', None) is not None: 931 self.converged = numpy.all(self.fcisolver.converged) 932 if self.converged: 933 log.info('CASCI converged') 934 else: 935 log.info('CASCI not converged') 936 else: 937 self.converged = True 938 self._finalize() 939 return self.e_tot, self.e_cas, self.ci, self.mo_coeff, self.mo_energy 940 941 def _finalize(self): 942 log = logger.Logger(self.stdout, self.verbose) 943 if log.verbose >= logger.NOTE and getattr(self.fcisolver, 'spin_square', None): 944 if isinstance(self.e_cas, (float, numpy.number)): 945 ss = self.fcisolver.spin_square(self.ci, self.ncas, self.nelecas) 946 log.note('CASCI E = %.15g E(CI) = %.15g S^2 = %.7f', 947 self.e_tot, self.e_cas, ss[0]) 948 else: 949 for i, e in enumerate(self.e_cas): 950 ss = self.fcisolver.spin_square(self.ci[i], self.ncas, self.nelecas) 951 log.note('CASCI state %d E = %.15g E(CI) = %.15g S^2 = %.7f', 952 i, self.e_tot[i], e, ss[0]) 953 else: 954 if isinstance(self.e_cas, (float, numpy.number)): 955 log.note('CASCI E = %.15g E(CI) = %.15g', self.e_tot, self.e_cas) 956 else: 957 for i, e in enumerate(self.e_cas): 958 log.note('CASCI state %d E = %.15g E(CI) = %.15g', 959 i, self.e_tot[i], e) 960 return self 961 962 as_scanner = as_scanner 963 964 @lib.with_doc(cas_natorb.__doc__) 965 def cas_natorb(self, mo_coeff=None, ci=None, eris=None, sort=False, 966 casdm1=None, verbose=None, with_meta_lowdin=WITH_META_LOWDIN): 967 return cas_natorb(self, mo_coeff, ci, eris, sort, casdm1, verbose, 968 with_meta_lowdin) 969 @lib.with_doc(cas_natorb.__doc__) 970 def cas_natorb_(self, mo_coeff=None, ci=None, eris=None, sort=False, 971 casdm1=None, verbose=None, with_meta_lowdin=WITH_META_LOWDIN): 972 self.mo_coeff, self.ci, self.mo_occ = cas_natorb(self, mo_coeff, ci, eris, 973 sort, casdm1, verbose) 974 return self.mo_coeff, self.ci, self.mo_occ 975 976 def get_fock(self, mo_coeff=None, ci=None, eris=None, casdm1=None, 977 verbose=None): 978 return get_fock(self, mo_coeff, ci, eris, casdm1, verbose) 979 980 canonicalize = canonicalize 981 982 @lib.with_doc(canonicalize.__doc__) 983 def canonicalize_(self, mo_coeff=None, ci=None, eris=None, sort=False, 984 cas_natorb=False, casdm1=None, verbose=None, 985 with_meta_lowdin=WITH_META_LOWDIN): 986 self.mo_coeff, ci, self.mo_energy = \ 987 canonicalize(self, mo_coeff, ci, eris, 988 sort, cas_natorb, casdm1, verbose, with_meta_lowdin) 989 if cas_natorb: # When active space is changed, the ci solution needs to be updated 990 self.ci = ci 991 return self.mo_coeff, ci, self.mo_energy 992 993 analyze = analyze 994 995 @lib.with_doc(addons.sort_mo.__doc__) 996 def sort_mo(self, caslst, mo_coeff=None, base=1): 997 if mo_coeff is None: mo_coeff = self.mo_coeff 998 return addons.sort_mo(self, mo_coeff, caslst, base) 999 1000 @lib.with_doc(addons.state_average.__doc__) 1001 def state_average_(self, weights=(0.5,0.5)): 1002 addons.state_average_(self, weights) 1003 return self 1004 @lib.with_doc(addons.state_average.__doc__) 1005 def state_average(self, weights=(0.5,0.5)): 1006 return addons.state_average(self, weights) 1007 1008 @lib.with_doc(addons.state_specific_.__doc__) 1009 def state_specific_(self, state=1): 1010 addons.state_specific(self, state) 1011 return self 1012 1013 def make_rdm1s(self, mo_coeff=None, ci=None, ncas=None, nelecas=None, 1014 ncore=None, **kwargs): 1015 '''One-particle density matrices for alpha and beta spin on AO basis 1016 ''' 1017 if mo_coeff is None: mo_coeff = self.mo_coeff 1018 if ci is None: ci = self.ci 1019 if ncas is None: ncas = self.ncas 1020 if nelecas is None: nelecas = self.nelecas 1021 if ncore is None: ncore = self.ncore 1022 1023 casdm1a, casdm1b = self.fcisolver.make_rdm1s(ci, ncas, nelecas) 1024 mocore = mo_coeff[:,:ncore] 1025 mocas = mo_coeff[:,ncore:ncore+ncas] 1026 dm1b = numpy.dot(mocore, mocore.conj().T) 1027 dm1a = dm1b + reduce(numpy.dot, (mocas, casdm1a, mocas.conj().T)) 1028 dm1b += reduce(numpy.dot, (mocas, casdm1b, mocas.conj().T)) 1029 return dm1a, dm1b 1030 1031 def make_rdm1(self, mo_coeff=None, ci=None, ncas=None, nelecas=None, 1032 ncore=None, **kwargs): 1033 '''One-particle density matrix in AO representation 1034 ''' 1035 if mo_coeff is None: mo_coeff = self.mo_coeff 1036 if ci is None: ci = self.ci 1037 if ncas is None: ncas = self.ncas 1038 if nelecas is None: nelecas = self.nelecas 1039 if ncore is None: ncore = self.ncore 1040 1041 casdm1 = self.fcisolver.make_rdm1(ci, ncas, nelecas) 1042 mocore = mo_coeff[:,:ncore] 1043 mocas = mo_coeff[:,ncore:ncore+ncas] 1044 dm1 = numpy.dot(mocore, mocore.conj().T) * 2 1045 dm1 = dm1 + reduce(numpy.dot, (mocas, casdm1, mocas.conj().T)) 1046 return dm1 1047 1048 def fix_spin_(self, shift=PENALTY, ss=None): 1049 r'''Use level shift to control FCI solver spin. 1050 1051 .. math:: 1052 1053 (H + shift*S^2) |\Psi\rangle = E |\Psi\rangle 1054 1055 Kwargs: 1056 shift : float 1057 Energy penalty for states which have wrong spin 1058 ss : number 1059 S^2 expection value == s*(s+1) 1060 ''' 1061 fci.addons.fix_spin_(self.fcisolver, shift, ss) 1062 return self 1063 fix_spin = fix_spin_ 1064 1065 def density_fit(self, auxbasis=None, with_df=None): 1066 from pyscf.mcscf import df 1067 return df.density_fit(self, auxbasis, with_df) 1068 1069 def sfx2c1e(self): 1070 from pyscf.x2c import sfx2c1e 1071 self._scf = sfx2c1e.sfx2c1e(self._scf).run() 1072 self.mo_coeff = self._scf.mo_coeff 1073 self.mo_energy = self._scf.mo_energy 1074 return self 1075 x2c = x2c1e = sfx2c1e 1076 1077 def nuc_grad_method(self): 1078 from pyscf.grad import casci 1079 return casci.Gradients(self) 1080 1081scf.hf.RHF.CASCI = scf.rohf.ROHF.CASCI = lib.class_as_method(CASCI) 1082scf.uhf.UHF.CASCI = None 1083 1084del(WITH_META_LOWDIN, LARGE_CI_TOL, PENALTY) 1085 1086 1087if __name__ == '__main__': 1088 from pyscf import mcscf 1089 mol = gto.Mole() 1090 mol.verbose = 0 1091 mol.output = None#"out_h2o" 1092 mol.atom = [ 1093 ['O', ( 0., 0. , 0. )], 1094 ['H', ( 0., -0.757, 0.587)], 1095 ['H', ( 0., 0.757 , 0.587)],] 1096 1097 mol.basis = {'H': 'sto-3g', 1098 'O': '6-31g',} 1099 mol.build() 1100 1101 m = scf.RHF(mol) 1102 ehf = m.scf() 1103 mc = mcscf.CASCI(m, 4, 4) 1104 mc.fcisolver = fci.solver(mol) 1105 mc.natorb = 1 1106 emc = mc.kernel()[0] 1107 print(ehf, emc, emc-ehf) 1108 #-75.9577817425 -75.9624554777 -0.00467373522233 1109 print(emc+75.9624554777) 1110 1111# mc = CASCI(m, 4, (3,1)) 1112# #mc.fcisolver = fci.direct_spin1 1113# mc.fcisolver = fci.solver(mol, False) 1114# emc = mc.casci()[0] 1115# print(emc - -75.439016172976) 1116# 1117# mol = gto.Mole() 1118# mol.verbose = 0 1119# mol.output = "out_casci" 1120# mol.atom = [ 1121# ["C", (-0.65830719, 0.61123287, -0.00800148)], 1122# ["C", ( 0.73685281, 0.61123287, -0.00800148)], 1123# ["C", ( 1.43439081, 1.81898387, -0.00800148)], 1124# ["C", ( 0.73673681, 3.02749287, -0.00920048)], 1125# ["C", (-0.65808819, 3.02741487, -0.00967948)], 1126# ["C", (-1.35568919, 1.81920887, -0.00868348)], 1127# ["H", (-1.20806619, -0.34108413, -0.00755148)], 1128# ["H", ( 1.28636081, -0.34128013, -0.00668648)], 1129# ["H", ( 2.53407081, 1.81906387, -0.00736748)], 1130# ["H", ( 1.28693681, 3.97963587, -0.00925948)], 1131# ["H", (-1.20821019, 3.97969587, -0.01063248)], 1132# ["H", (-2.45529319, 1.81939187, -0.00886348)],] 1133# 1134# mol.basis = {'H': 'sto-3g', 1135# 'C': 'sto-3g',} 1136# mol.build() 1137# 1138# m = scf.RHF(mol) 1139# ehf = m.scf() 1140# mc = CASCI(m, 9, 8) 1141# mc.fcisolver = fci.solver(mol) 1142# emc = mc.casci()[0] 1143# print(ehf, emc, emc-ehf) 1144# print(emc - -227.948912536) 1145# 1146# mc = CASCI(m, 9, (5,3)) 1147# #mc.fcisolver = fci.direct_spin1 1148# mc.fcisolver = fci.solver(mol, False) 1149# mc.fcisolver.nroots = 3 1150# emc = mc.casci()[0] 1151# print(emc[0] - -227.7674519720) 1152