1#!/usr/bin/env python 2# Copyright 2014-2021 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''' 20FCI solver for Singlet state 21 22Different FCI solvers are implemented to support different type of symmetry. 23 Symmetry 24File Point group Spin singlet Real hermitian* Alpha/beta degeneracy 25direct_spin0_symm Yes Yes Yes Yes 26direct_spin1_symm Yes No Yes Yes 27direct_spin0 No Yes Yes Yes 28direct_spin1 No No Yes Yes 29direct_uhf No No Yes No 30direct_nosym No No No** Yes 31 32* Real hermitian Hamiltonian implies (ij|kl) = (ji|kl) = (ij|lk) = (ji|lk) 33** Hamiltonian is real but not hermitian, (ij|kl) != (ji|kl) ... 34 35direct_spin0 solver is specified for singlet state. However, calling this 36solver sometimes ends up with the error "State not singlet x.xxxxxxe-06" due 37to numerical issues. Calling direct_spin1 for singlet state is slightly 38slower but more robust than direct_spin0 especially when combining to energy 39penalty method (:func:`fix_spin_`) 40''' 41 42import ctypes 43import numpy 44from pyscf import lib 45from pyscf import ao2mo 46from pyscf.lib import logger 47from pyscf.fci import cistring 48from pyscf.fci import rdm 49from pyscf.fci import direct_spin1 50from pyscf.fci.spin_op import contract_ss 51 52libfci = lib.load_library('libfci') 53 54@lib.with_doc(direct_spin1.contract_1e.__doc__) 55def contract_1e(f1e, fcivec, norb, nelec, link_index=None): 56 fcivec = numpy.asarray(fcivec, order='C') 57 link_index = _unpack(norb, nelec, link_index) 58 na, nlink = link_index.shape[:2] 59 assert(fcivec.size == na**2) 60 ci1 = numpy.empty_like(fcivec) 61 f1e_tril = lib.pack_tril(f1e) 62 libfci.FCIcontract_1e_spin0(f1e_tril.ctypes.data_as(ctypes.c_void_p), 63 fcivec.ctypes.data_as(ctypes.c_void_p), 64 ci1.ctypes.data_as(ctypes.c_void_p), 65 ctypes.c_int(norb), ctypes.c_int(na), 66 ctypes.c_int(nlink), 67 link_index.ctypes.data_as(ctypes.c_void_p)) 68# no *.5 because FCIcontract_2e_spin0 only compute half of the contraction 69 return lib.transpose_sum(ci1, inplace=True).reshape(fcivec.shape) 70 71# Note eri is NOT the 2e hamiltonian matrix, the 2e hamiltonian is 72# h2e = eri_{pq,rs} p^+ q r^+ s 73# = (pq|rs) p^+ r^+ s q - (pq|rs) \delta_{qr} p^+ s 74# so eri is defined as 75# eri_{pq,rs} = (pq|rs) - (1/Nelec) \sum_q (pq|qs) 76# to restore the symmetry between pq and rs, 77# eri_{pq,rs} = (pq|rs) - (.5/Nelec) [\sum_q (pq|qs) + \sum_p (pq|rp)] 78# Please refer to the treatment in direct_spin1.absorb_h1e 79# the input fcivec should be symmetrized 80@lib.with_doc(direct_spin1.contract_2e.__doc__) 81def contract_2e(eri, fcivec, norb, nelec, link_index=None): 82 fcivec = numpy.asarray(fcivec, order='C') 83 eri = ao2mo.restore(4, eri, norb) 84 lib.transpose_sum(eri, inplace=True) 85 eri *= .5 86 link_index = _unpack(norb, nelec, link_index) 87 na, nlink = link_index.shape[:2] 88 assert(fcivec.size == na**2) 89 ci1 = numpy.empty((na,na)) 90 91 libfci.FCIcontract_2e_spin0(eri.ctypes.data_as(ctypes.c_void_p), 92 fcivec.ctypes.data_as(ctypes.c_void_p), 93 ci1.ctypes.data_as(ctypes.c_void_p), 94 ctypes.c_int(norb), ctypes.c_int(na), 95 ctypes.c_int(nlink), 96 link_index.ctypes.data_as(ctypes.c_void_p)) 97# no *.5 because FCIcontract_2e_spin0 only compute half of the contraction 98 return lib.transpose_sum(ci1, inplace=True).reshape(fcivec.shape) 99 100absorb_h1e = direct_spin1.absorb_h1e 101 102@lib.with_doc(direct_spin1.make_hdiag.__doc__) 103def make_hdiag(h1e, eri, norb, nelec): 104 hdiag = direct_spin1.make_hdiag(h1e, eri, norb, nelec) 105 na = int(numpy.sqrt(hdiag.size)) 106# symmetrize hdiag to reduce numerical error 107 hdiag = lib.transpose_sum(hdiag.reshape(na,na), inplace=True) * .5 108 return hdiag.ravel() 109 110pspace = direct_spin1.pspace 111 112# be careful with single determinant initial guess. It may lead to the 113# eigvalue of first davidson iter being equal to hdiag 114def kernel(h1e, eri, norb, nelec, ci0=None, level_shift=1e-3, tol=1e-10, 115 lindep=1e-14, max_cycle=50, max_space=12, nroots=1, 116 davidson_only=False, pspace_size=400, orbsym=None, wfnsym=None, 117 ecore=0, **kwargs): 118 e, c = direct_spin1._kfactory(FCISolver, h1e, eri, norb, nelec, ci0, level_shift, 119 tol, lindep, max_cycle, max_space, nroots, 120 davidson_only, pspace_size, ecore=ecore, **kwargs) 121 return e, c 122 123# dm[p,q] = <|q^+ p|> 124@lib.with_doc(direct_spin1.make_rdm1.__doc__) 125def make_rdm1(fcivec, norb, nelec, link_index=None): 126 rdm1 = rdm.make_rdm1('FCImake_rdm1a', fcivec, fcivec, 127 norb, nelec, link_index) 128 return rdm1 * 2 129 130# alpha and beta 1pdm 131@lib.with_doc(direct_spin1.make_rdm1s.__doc__) 132def make_rdm1s(fcivec, norb, nelec, link_index=None): 133 rdm1 = rdm.make_rdm1('FCImake_rdm1a', fcivec, fcivec, 134 norb, nelec, link_index) 135 return rdm1, rdm1 136 137# Chemist notation 138@lib.with_doc(direct_spin1.make_rdm12.__doc__) 139def make_rdm12(fcivec, norb, nelec, link_index=None, reorder=True): 140 #dm1, dm2 = rdm.make_rdm12('FCIrdm12kern_spin0', fcivec, fcivec, 141 # norb, nelec, link_index, 1) 142 143 # NOT use FCIrdm12kern_spin0 because for small system, the kernel may call 144 # direct diagonalization, which may not fulfil fcivec = fcivet.T 145 dm1, dm2 = rdm.make_rdm12('FCIrdm12kern_sf', fcivec, fcivec, 146 norb, nelec, link_index, 1) 147 if reorder: 148 dm1, dm2 = rdm.reorder_rdm(dm1, dm2, True) 149 return dm1, dm2 150 151# dm[p,q] = <I|q^+ p|J> 152@lib.with_doc(direct_spin1.trans_rdm1s.__doc__) 153def trans_rdm1s(cibra, ciket, norb, nelec, link_index=None): 154 if link_index is None: 155 if isinstance(nelec, (int, numpy.number)): 156 neleca = nelec//2 157 else: 158 neleca, nelecb = nelec 159 assert(neleca == nelecb) 160 link_index = cistring.gen_linkstr_index(range(norb), neleca) 161 rdm1a = rdm.make_rdm1('FCItrans_rdm1a', cibra, ciket, 162 norb, nelec, link_index) 163 rdm1b = rdm.make_rdm1('FCItrans_rdm1b', cibra, ciket, 164 norb, nelec, link_index) 165 return rdm1a, rdm1b 166 167@lib.with_doc(direct_spin1.trans_rdm1.__doc__) 168def trans_rdm1(cibra, ciket, norb, nelec, link_index=None): 169 rdm1a, rdm1b = trans_rdm1s(cibra, ciket, norb, nelec, link_index) 170 return rdm1a + rdm1b 171 172# dm[p,q,r,s] = <I|p^+ q r^+ s|J> 173@lib.with_doc(direct_spin1.trans_rdm12.__doc__) 174def trans_rdm12(cibra, ciket, norb, nelec, link_index=None, reorder=True): 175 dm1, dm2 = rdm.make_rdm12('FCItdm12kern_sf', cibra, ciket, 176 norb, nelec, link_index, 2) 177 if reorder: 178 dm1, dm2 = rdm.reorder_rdm(dm1, dm2, True) 179 return dm1, dm2 180 181def energy(h1e, eri, fcivec, norb, nelec, link_index=None): 182 h2e = direct_spin1.absorb_h1e(h1e, eri, norb, nelec, .5) 183 ci1 = contract_2e(h2e, fcivec, norb, nelec, link_index) 184 return numpy.dot(fcivec.ravel(), ci1.ravel()) 185 186def get_init_guess(norb, nelec, nroots, hdiag): 187 if isinstance(nelec, (int, numpy.number)): 188 nelecb = nelec//2 189 neleca = nelec - nelecb 190 else: 191 neleca, nelecb = nelec 192 na = cistring.num_strings(norb, neleca) 193 nb = cistring.num_strings(norb, nelecb) 194 195 init_strs = [] 196 iroot = 0 197 for addr in numpy.argsort(hdiag): 198 addra = addr // nb 199 addrb = addr % nb 200 if (addrb,addra) not in init_strs: # avoid initial guess linear dependency 201 init_strs.append((addra,addrb)) 202 iroot += 1 203 if iroot >= nroots: 204 break 205 ci0 = [] 206 for addra,addrb in init_strs: 207 x = numpy.zeros((na,nb)) 208 if addra == addrb: 209 x[addra,addrb] = 1 210 else: 211 x[addra,addrb] = x[addrb,addra] = numpy.sqrt(.5) 212 ci0.append(x.ravel()) 213 214 # Add noise 215 ci0[0][0 ] += 1e-5 216 ci0[0][-1] -= 1e-5 217 return ci0 218 219 220############################################################### 221# direct-CI driver 222############################################################### 223 224def kernel_ms0(fci, h1e, eri, norb, nelec, ci0=None, link_index=None, 225 tol=None, lindep=None, max_cycle=None, max_space=None, 226 nroots=None, davidson_only=None, pspace_size=None, 227 max_memory=None, verbose=None, ecore=0, **kwargs): 228 if nroots is None: nroots = fci.nroots 229 if davidson_only is None: davidson_only = fci.davidson_only 230 if pspace_size is None: pspace_size = fci.pspace_size 231 if max_memory is None: 232 max_memory = fci.max_memory - lib.current_memory()[0] 233 log = logger.new_logger(fci, verbose) 234 235 assert(fci.spin is None or fci.spin == 0) 236 assert(0 <= numpy.sum(nelec) <= norb*2) 237 238 link_index = _unpack(norb, nelec, link_index) 239 h1e = numpy.ascontiguousarray(h1e) 240 eri = numpy.ascontiguousarray(eri) 241 na = link_index.shape[0] 242 243 if max_memory < na**2*6*8e-6: 244 log.warn('Not enough memory for FCI solver. ' 245 'The minimal requirement is %.0f MB', na**2*60e-6) 246 247 hdiag = fci.make_hdiag(h1e, eri, norb, nelec) 248 nroots = min(hdiag.size, nroots) 249 250 try: 251 addr, h0 = fci.pspace(h1e, eri, norb, nelec, hdiag, max(pspace_size,nroots)) 252 if pspace_size > 0: 253 pw, pv = fci.eig(h0) 254 else: 255 pw = pv = None 256 257 if pspace_size >= na*na and ci0 is None and not davidson_only: 258 # The degenerated wfn can break symmetry. The davidson iteration with proper 259 # initial guess doesn't have this issue 260 if na*na == 1: 261 return pw[0]+ecore, pv[:,0].reshape(1,1) 262 elif nroots > 1: 263 civec = numpy.empty((nroots,na*na)) 264 civec[:,addr] = pv[:,:nroots].T 265 civec = civec.reshape(nroots,na,na) 266 try: 267 return pw[:nroots]+ecore, [_check_(ci) for ci in civec] 268 except ValueError: 269 pass 270 elif abs(pw[0]-pw[1]) > 1e-12: 271 civec = numpy.empty((na*na)) 272 civec[addr] = pv[:,0] 273 civec = civec.reshape(na,na) 274 civec = lib.transpose_sum(civec) * .5 275 # direct diagonalization may lead to triplet ground state 276 277 #TODO: optimize initial guess. Using pspace vector as initial guess may have 278 # spin problems. The 'ground state' of psapce vector may have different spin 279 # state to the true ground state. 280 try: 281 return pw[0]+ecore, _check_(civec.reshape(na,na)) 282 except ValueError: 283 pass 284 except NotImplementedError: 285 addr = [0] 286 pw = pv = None 287 288 precond = fci.make_precond(hdiag, pw, pv, addr) 289 290 h2e = fci.absorb_h1e(h1e, eri, norb, nelec, .5) 291 def hop(c): 292 hc = fci.contract_2e(h2e, c.reshape(na,na), norb, nelec, link_index) 293 return hc.ravel() 294 295#TODO: check spin of initial guess 296 if ci0 is None: 297 if callable(getattr(fci, 'get_init_guess', None)): 298 ci0 = lambda: fci.get_init_guess(norb, nelec, nroots, hdiag) 299 else: 300 def ci0(): 301 x0 = [] 302 for i in range(nroots): 303 x = numpy.zeros((na,na)) 304 addra = addr[i] // na 305 addrb = addr[i] % na 306 if addra == addrb: 307 x[addra,addrb] = 1 308 else: 309 x[addra,addrb] = x[addrb,addra] = numpy.sqrt(.5) 310 x0.append(x.ravel()) 311 return x0 312 elif not callable(ci0): 313 if isinstance(ci0, numpy.ndarray) and ci0.size == na*na: 314 ci0 = [ci0.ravel()] 315 else: 316 ci0 = [x.ravel() for x in ci0] 317 318 if tol is None: tol = fci.conv_tol 319 if lindep is None: lindep = fci.lindep 320 if max_cycle is None: max_cycle = fci.max_cycle 321 if max_space is None: max_space = fci.max_space 322 tol_residual = getattr(fci, 'conv_tol_residual', None) 323 324 with lib.with_omp_threads(fci.threads): 325 #e, c = lib.davidson(hop, ci0, precond, tol=fci.conv_tol, lindep=fci.lindep) 326 e, c = fci.eig(hop, ci0, precond, tol=tol, lindep=lindep, 327 max_cycle=max_cycle, max_space=max_space, nroots=nroots, 328 max_memory=max_memory, verbose=log, follow_state=True, 329 tol_residual=tol_residual, **kwargs) 330 if nroots > 1: 331 return e+ecore, [_check_(ci.reshape(na,na)) for ci in c] 332 else: 333 return e+ecore, _check_(c.reshape(na,na)) 334 335def _check_(c): 336 c = lib.transpose_sum(c, inplace=True) 337 c *= .5 338 norm = numpy.linalg.norm(c) 339 if abs(norm-1) > 1e-6: 340 raise ValueError('State not singlet %g' % (norm - 1)) 341 return c/norm 342 343 344class FCISolver(direct_spin1.FCISolver): 345 346 def make_hdiag(self, h1e, eri, norb, nelec): 347 return make_hdiag(h1e, eri, norb, nelec) 348 349 def contract_1e(self, f1e, fcivec, norb, nelec, link_index=None, **kwargs): 350 return contract_1e(f1e, fcivec, norb, nelec, link_index, **kwargs) 351 352 def contract_2e(self, eri, fcivec, norb, nelec, link_index=None, **kwargs): 353 return contract_2e(eri, fcivec, norb, nelec, link_index, **kwargs) 354 355 def get_init_guess(self, norb, nelec, nroots, hdiag): 356 return get_init_guess(norb, nelec, nroots, hdiag) 357 358 def kernel(self, h1e, eri, norb, nelec, ci0=None, 359 tol=None, lindep=None, max_cycle=None, max_space=None, 360 nroots=None, davidson_only=None, pspace_size=None, 361 orbsym=None, wfnsym=None, ecore=0, **kwargs): 362 if self.verbose >= logger.WARN: 363 self.check_sanity() 364 self.norb = norb 365 self.nelec = nelec 366 self.eci, self.ci = \ 367 kernel_ms0(self, h1e, eri, norb, nelec, ci0, None, 368 tol, lindep, max_cycle, max_space, nroots, 369 davidson_only, pspace_size, ecore=ecore, **kwargs) 370 return self.eci, self.ci 371 372 def energy(self, h1e, eri, fcivec, norb, nelec, link_index=None): 373 h2e = self.absorb_h1e(h1e, eri, norb, nelec, .5) 374 ci1 = self.contract_2e(h2e, fcivec, norb, nelec, link_index) 375 return numpy.dot(fcivec.reshape(-1), ci1.reshape(-1)) 376 377 def make_rdm1s(self, fcivec, norb, nelec, link_index=None): 378 return make_rdm1s(fcivec, norb, nelec, link_index) 379 380 def make_rdm1(self, fcivec, norb, nelec, link_index=None): 381 return make_rdm1(fcivec, norb, nelec, link_index) 382 383 @lib.with_doc(make_rdm12.__doc__) 384 def make_rdm12(self, fcivec, norb, nelec, link_index=None, reorder=True): 385 return make_rdm12(fcivec, norb, nelec, link_index, reorder) 386 387 def trans_rdm1s(self, cibra, ciket, norb, nelec, link_index=None): 388 return trans_rdm1s(cibra, ciket, norb, nelec, link_index) 389 390 def trans_rdm1(self, cibra, ciket, norb, nelec, link_index=None): 391 return trans_rdm1(cibra, ciket, norb, nelec, link_index) 392 393 @lib.with_doc(trans_rdm12.__doc__) 394 def trans_rdm12(self, cibra, ciket, norb, nelec, link_index=None, 395 reorder=True): 396 return trans_rdm12(cibra, ciket, norb, nelec, link_index, reorder) 397 398 def gen_linkstr(self, norb, nelec, tril=True, spin=None): 399 if isinstance(nelec, (int, numpy.number)): 400 neleca = nelec//2 401 else: 402 neleca, nelecb = nelec 403 assert(neleca == nelecb) 404 if tril: 405 link_index = cistring.gen_linkstr_index_trilidx(range(norb), neleca) 406 else: 407 link_index = cistring.gen_linkstr_index(range(norb), neleca) 408 return link_index 409 410FCI = FCISolver 411 412 413def _unpack(norb, nelec, link_index): 414 if link_index is None: 415 if isinstance(nelec, (int, numpy.number)): 416 neleca = nelec//2 417 else: 418 neleca, nelecb = nelec 419 assert(neleca == nelecb) 420 return cistring.gen_linkstr_index_trilidx(range(norb), neleca) 421 else: 422 return link_index 423 424 425if __name__ == '__main__': 426 from functools import reduce 427 from pyscf import gto 428 from pyscf import scf 429 430 mol = gto.Mole() 431 mol.verbose = 0 432 mol.output = None#"out_h2o" 433 mol.atom = [ 434 ['H', ( 1.,-1. , 0. )], 435 ['H', ( 0.,-1. ,-1. )], 436 ['H', ( 1.,-0.5 ,-1. )], 437 ['H', ( 0.,-0.5 ,-1. )], 438 ['H', ( 0.,-0.5 ,-0. )], 439 ['H', ( 0.,-0. ,-1. )], 440 ['H', ( 1.,-0.5 , 0. )], 441 ['H', ( 0., 1. , 1. )], 442 ] 443 444 mol.basis = {'H': 'sto-3g'} 445 mol.build() 446 447 m = scf.RHF(mol) 448 ehf = m.scf() 449 450 cis = FCISolver(mol) 451 norb = m.mo_coeff.shape[1] 452 nelec = mol.nelectron 453 h1e = reduce(numpy.dot, (m.mo_coeff.T, m.get_hcore(), m.mo_coeff)) 454 eri = ao2mo.incore.general(m._eri, (m.mo_coeff,)*4, compact=False) 455 e, c = cis.kernel(h1e, eri, norb, nelec) 456 print(e - -15.9977886375) 457 print('t',logger.process_clock()) 458 459