Python numpy.einsum() Examples
The following are 30
code examples of numpy.einsum().
You can vote up the ones you like or vote down the ones you don't like,
and go to the original project or source file by following the links above each example.
You may also want to check out all available functions/classes of the module
numpy
, or try the search function
.
Example #1
Source File: 30-force_on_mm_particles.py From pyscf with Apache License 2.0 | 6 votes |
def force(dm): # The interaction between QM atoms and MM particles # \sum_K d/dR (1/|r_K-R|) = \sum_K (r_K-R)/|r_K-R|^3 qm_coords = mol.atom_coords() qm_charges = mol.atom_charges() dr = qm_coords[:,None,:] - coords r = numpy.linalg.norm(dr, axis=2) g = numpy.einsum('r,R,rRx,rR->Rx', qm_charges, charges, dr, r**-3) # The interaction between electron density and MM particles # d/dR <i| (1/|r-R|) |j> = <i| d/dR (1/|r-R|) |j> = <i| -d/dr (1/|r-R|) |j> # = <d/dr i| (1/|r-R|) |j> + <i| (1/|r-R|) |d/dr j> for i, q in enumerate(charges): with mol.with_rinv_origin(coords[i]): v = mol.intor('int1e_iprinv') f =(numpy.einsum('ij,xji->x', dm, v) + numpy.einsum('ij,xij->x', dm, v.conj())) * -q g[i] += f # Force = -d/dR return -g # The force from HF electron density # Be careful with the unit of the MM particle coordinates. The gradients are # computed in the atomic unit.
Example #2
Source File: sfx2c1e_hess.py From pyscf with Apache License 2.0 | 6 votes |
def hcore_hess_generator(x2cobj, mol=None): '''nuclear gradients of 1-component X2c hcore Hamiltonian (spin-free part only) ''' if mol is None: mol = x2cobj.mol xmol, contr_coeff = x2cobj.get_xmol(mol) if x2cobj.basis is not None: s22 = xmol.intor_symmetric('int1e_ovlp') s21 = gto.intor_cross('int1e_ovlp', xmol, mol) contr_coeff = lib.cho_solve(s22, s21) get_h1_xmol = gen_sf_hfw(xmol, x2cobj.approx) def hcore_deriv(ia, ja): h1 = get_h1_xmol(ia, ja) if contr_coeff is not None: h1 = lib.einsum('pi,xypq,qj->xyij', contr_coeff, h1, contr_coeff) return numpy.asarray(h1) return hcore_deriv
Example #3
Source File: test_x2c_grad.py From pyscf with Apache License 2.0 | 6 votes |
def get_x1(mol, ia): h0, s0 = get_h0_s0(mol) h1, s1 = get_h1_s1(mol, ia) e0, c0 = scipy.linalg.eigh(h0, s0) nao = mol.nao_nr() cl0 = c0[:nao,nao:] cs0 = c0[nao:,nao:] x0 = scipy.linalg.solve(cl0.T, cs0.T).T h1 = numpy.einsum('pi,xpq,qj->xij', c0.conj(), h1, c0[:,nao:]) s1 = numpy.einsum('pi,xpq,qj->xij', c0.conj(), s1, c0[:,nao:]) epi = e0[:,None] - e0[nao:] degen_mask = abs(epi) < 1e-7 epi[degen_mask] = 1e200 c1 = (h1 - s1 * e0[nao:]) / -epi c1[:,degen_mask] = -.5 * s1[:,degen_mask] c1 = numpy.einsum('pq,xqi->xpi', c0, c1) cl1 = c1[:,:nao] cs1 = c1[:,nao:] x1 = [scipy.linalg.solve(cl0.T, (cs1[i] - x0.dot(cl1[i])).T).T for i in range(3)] return numpy.asarray(x1)
Example #4
Source File: test_x2c_hess.py From pyscf with Apache License 2.0 | 6 votes |
def _invsqrt2(a0, a1i, a1j, a2ij): '''Solving first order derivative of x^2 = a^{-1}''' w, v = scipy.linalg.eigh(a0) w = 1./numpy.sqrt(w) a1i = reduce(numpy.dot, (v.conj().T, a1i, v)) x1i = numpy.einsum('i,ij,j->ij', w**2, a1i, w**2) / (w[:,None] + w) a1j = reduce(numpy.dot, (v.conj().T, a1j, v)) x1j = numpy.einsum('i,ij,j->ij', w**2, a1j, w**2) / (w[:,None] + w) a2ij = reduce(numpy.dot, (v.conj().T, a2ij, v)) tmp = (a1i*w**2).dot(a1j) a2ij -= tmp + tmp.conj().T a2ij = -numpy.einsum('i,ij,j->ij', w**2, a2ij, w**2) tmp = x1i.dot(x1j) a2ij -= tmp + tmp.conj().T x2 = a2ij / (w[:,None] + w) x2 = reduce(numpy.dot, (v, x2, v.conj().T)) return x2
Example #5
Source File: ModelingCloth.py From Modeling-Cloth with MIT License | 6 votes |
def generate_inflate(cloth): """Blow it up baby!""" tri_nor = cloth.normals #* cloth.ob.modeling_cloth_inflate # non-unit calculated by tri_normals_in_place() per each triangle #tri_nor /= np.einsum("ij, ij->i", tri_nor, tri_nor)[:, nax] # reshape for add.at shape = cloth.inflate.shape cloth.inflate += tri_nor[:, nax] * cloth.ob.modeling_cloth_inflate# * cloth.tri_mix cloth.inflate.shape = (shape[0] * 3, 3) cloth.inflate *= cloth.tri_mix np.add.at(cloth.vel, cloth.tridex.ravel(), cloth.inflate) cloth.inflate.shape = shape cloth.inflate *= 0
Example #6
Source File: hf.py From pyscf with Apache License 2.0 | 6 votes |
def normalize_dm_(mf, dm): ''' Scale density matrix to make it produce the correct number of electrons. ''' cell = mf.cell if isinstance(dm, np.ndarray) and dm.ndim == 2: ne = np.einsum('ij,ji->', dm, mf.get_ovlp(cell)).real else: ne = np.einsum('xij,ji->', dm, mf.get_ovlp(cell)).real if abs(ne - cell.nelectron).sum() > 1e-7: logger.debug(mf, 'Big error detected in the electron number ' 'of initial guess density matrix (Ne/cell = %g)!\n' ' This can cause huge error in Fock matrix and ' 'lead to instability in SCF for low-dimensional ' 'systems.\n DM is normalized wrt the number ' 'of electrons %s', ne, cell.nelectron) dm *= cell.nelectron / ne return dm
Example #7
Source File: kuhf.py From pyscf with Apache License 2.0 | 6 votes |
def spin_square(self, mo_coeff=None, s=None): '''Treating the k-point sampling wfn as a giant Slater determinant, the spin_square value is the <S^2> of the giant determinant. ''' nkpts = len(self.kpts) if mo_coeff is None: mo_a = [self.mo_coeff[0][k][:,self.mo_occ[0][k]>0] for k in range(nkpts)] mo_b = [self.mo_coeff[1][k][:,self.mo_occ[1][k]>0] for k in range(nkpts)] else: mo_a, mo_b = mo_coeff if s is None: s = self.get_ovlp() nelec_a = sum([mo_a[k].shape[1] for k in range(nkpts)]) nelec_b = sum([mo_b[k].shape[1] for k in range(nkpts)]) ssxy = (nelec_a + nelec_b) * .5 for k in range(nkpts): sij = reduce(np.dot, (mo_a[k].T.conj(), s[k], mo_b[k])) ssxy -= np.einsum('ij,ij->', sij.conj(), sij).real ssz = (nelec_b-nelec_a)**2 * .25 ss = ssxy + ssz s = np.sqrt(ss+.25) - .5 return ss, s*2+1
Example #8
Source File: khf.py From pyscf with Apache License 2.0 | 6 votes |
def energy_elec(mf, dm_kpts=None, h1e_kpts=None, vhf_kpts=None): '''Following pyscf.scf.hf.energy_elec() ''' if dm_kpts is None: dm_kpts = mf.make_rdm1() if h1e_kpts is None: h1e_kpts = mf.get_hcore() if vhf_kpts is None: vhf_kpts = mf.get_veff(mf.cell, dm_kpts) nkpts = len(dm_kpts) e1 = 1./nkpts * np.einsum('kij,kji', dm_kpts, h1e_kpts) e_coul = 1./nkpts * np.einsum('kij,kji', dm_kpts, vhf_kpts) * 0.5 mf.scf_summary['e1'] = e1.real mf.scf_summary['e2'] = e_coul.real logger.debug(mf, 'E1 = %s E_coul = %s', e1, e_coul) if CHECK_COULOMB_IMAG and abs(e_coul.imag > mf.cell.precision*10): logger.warn(mf, "Coulomb energy has imaginary part %s. " "Coulomb integrals (e-e, e-N) may not converge !", e_coul.imag) return (e1+e_coul).real, e_coul.real
Example #9
Source File: 22-density.py From pyscf with Apache License 2.0 | 6 votes |
def tda_denisty_matrix(td, state_id): ''' Taking the TDA amplitudes as the CIS coefficients, calculate the density matrix (in AO basis) of the excited states ''' cis_t1 = td.xy[state_id][0] dm_oo =-np.einsum('ia,ka->ik', cis_t1.conj(), cis_t1) dm_vv = np.einsum('ia,ic->ac', cis_t1, cis_t1.conj()) # The ground state density matrix in mo_basis mf = td._scf dm = np.diag(mf.mo_occ) # Add CIS contribution nocc = cis_t1.shape[0] dm[:nocc,:nocc] += dm_oo * 2 dm[nocc:,nocc:] += dm_vv * 2 # Transform density matrix to AO basis mo = mf.mo_coeff dm = np.einsum('pi,ij,qj->pq', mo, dm, mo.conj()) return dm # Density matrix for the 3rd excited state
Example #10
Source File: sfx2c1e_grad.py From pyscf with Apache License 2.0 | 6 votes |
def hcore_grad_generator(x2cobj, mol=None): '''nuclear gradients of 1-component X2c hcore Hamiltonian (spin-free part only) ''' if mol is None: mol = x2cobj.mol xmol, contr_coeff = x2cobj.get_xmol(mol) if x2cobj.basis is not None: s22 = xmol.intor_symmetric('int1e_ovlp') s21 = gto.intor_cross('int1e_ovlp', xmol, mol) contr_coeff = lib.cho_solve(s22, s21) get_h1_xmol = gen_sf_hfw(xmol, x2cobj.approx) def hcore_deriv(atm_id): h1 = get_h1_xmol(atm_id) if contr_coeff is not None: h1 = lib.einsum('pi,xpq,qj->xij', contr_coeff, h1, contr_coeff) return numpy.asarray(h1) return hcore_deriv
Example #11
Source File: from_arrays.py From QCElemental with BSD 3-Clause "New" or "Revised" License | 6 votes |
def validate_and_fill_geometry(geom=None, tooclose=0.1, copy=True): """Check `geom` for overlapping atoms. Return flattened""" npgeom = np.array(geom, copy=copy, dtype=np.float).reshape((-1, 3)) # Upper triangular metric = tooclose ** 2 tooclose_inds = [] for x in range(npgeom.shape[0]): diffs = npgeom[x] - npgeom[x + 1 :] dists = np.einsum("ij,ij->i", diffs, diffs) # Record issues if np.any(dists < metric): indices = np.where(dists < metric)[0] tooclose_inds.extend([(x, y, dist) for y, dist in zip(indices + x + 1, dists[indices] ** 0.5)]) if tooclose_inds: raise ValidationError( """Following atoms are too close: {}""".format([(i, j, dist) for i, j, dist in tooclose_inds]) ) return {"geom": npgeom.reshape((-1))}
Example #12
Source File: kmp2.py From pyscf with Apache License 2.0 | 6 votes |
def _gamma1_intermediates(mp, t2=None): # Memory optimization should be here if t2 is None: t2 = mp.t2 if t2 is None: raise NotImplementedError("Run kmp2.kernel with `with_t2=True`") nmo = mp.nmo nocc = mp.nocc nvir = nmo - nocc nkpts = mp.nkpts dtype = t2.dtype dm1occ = np.zeros((nkpts, nocc, nocc), dtype=dtype) dm1vir = np.zeros((nkpts, nvir, nvir), dtype=dtype) for ki in range(nkpts): for kj in range(nkpts): for ka in range(nkpts): kb = mp.khelper.kconserv[ki, ka, kj] dm1vir[kb] += einsum('ijax,ijay->yx', t2[ki][kj][ka].conj(), t2[ki][kj][ka]) * 2 -\ einsum('ijax,ijya->yx', t2[ki][kj][ka].conj(), t2[ki][kj][kb]) dm1occ[kj] += einsum('ixab,iyab->xy', t2[ki][kj][ka].conj(), t2[ki][kj][ka]) * 2 -\ einsum('ixab,iyba->xy', t2[ki][kj][ka].conj(), t2[ki][kj][kb]) return -dm1occ, dm1vir
Example #13
Source File: 11-ump2.py From pyscf with Apache License 2.0 | 6 votes |
def myump2(mol, mo_energy, mo_coeff, mo_occ): o = numpy.hstack((mo_coeff[0][:,mo_occ[0]>0] ,mo_coeff[1][:,mo_occ[1]>0])) v = numpy.hstack((mo_coeff[0][:,mo_occ[0]==0],mo_coeff[1][:,mo_occ[1]==0])) eo = numpy.hstack((mo_energy[0][mo_occ[0]>0] ,mo_energy[1][mo_occ[1]>0])) ev = numpy.hstack((mo_energy[0][mo_occ[0]==0],mo_energy[1][mo_occ[1]==0])) no = o.shape[1] nv = v.shape[1] noa = sum(mo_occ[0]>0) nva = sum(mo_occ[0]==0) eri = ao2mo.general(mol, (o,v,o,v)).reshape(no,nv,no,nv) eri[:noa,nva:] = eri[noa:,:nva] = eri[:,:,:noa,nva:] = eri[:,:,noa:,:nva] = 0 g = eri - eri.transpose(0,3,2,1) eov = eo.reshape(-1,1) - ev.reshape(-1) de = 1/(eov.reshape(-1,1) + eov.reshape(-1)).reshape(g.shape) emp2 = .25 * numpy.einsum('iajb,iajb,iajb->', g, g, de) return emp2
Example #14
Source File: 20-soc_ao_integrals.py From pyscf with Apache License 2.0 | 6 votes |
def ss(mol): n = mol.nao_nr() mat1 = mol.intor('int2e_ip1ip2_sph').reshape(3,3,n,n,n,n) # <nabla1 nabla2 | 1 2> mat2 =-mat1.transpose(0,1,2,3,5,4) # <nabla1 2 | 1 nabla2> mat3 =-mat2.transpose(1,0,3,2,4,5) # <1 nabla2 | nabla1 2> mat4 = mat1.transpose(0,1,3,2,5,4) # <1 2 | nabla1 nabla2> mat = mat1 - mat2 - mat3 + mat4 # Fermi contact term h_fc = mol.intor('int4c1e').reshape(nao,nao,nao,nao) * (4*numpy.pi/3) mat[0,0] -= h_fc mat[1,1] -= h_fc mat[2,2] -= h_fc s = lib.PauliMatrices * .5 # wxyz are the spin indices, ijkl are the AO indicies alpha = 137.036 fac = alpha ** 2 / 2 mat = numpy.einsum('swx,tyz,stijkl->wxyzijkl', s[:,0,0], s[:,0,0], mat) * fac return mat
Example #15
Source File: eom_kccsd_ghf.py From pyscf with Apache License 2.0 | 5 votes |
def get_kconserv_ee_r2(self, kshift=0): r'''Get the momentum conservation array for a set of k-points. Given k-point indices (k, l, m) the array kconserv_r2[k,l,m] returns the index n that satisfies momentum conservation, (k(k) - k(l) + k(m) - k(n) - kshift) \dot a = 2n\pi This is used for symmetry of 2p-2h excitation operator vector R_{k k_k, m k_m}^{l k_l n k_n} is zero unless n satisfies the above. Note that this method is adapted from `kpts_helper.get_kconserv()`. ''' cell = self._cc._scf.cell kpts = self.kpts nkpts = kpts.shape[0] a = cell.lattice_vectors() / (2 * np.pi) kconserv_r2 = np.zeros((nkpts, nkpts, nkpts), dtype=int) kvKLM = kpts[:, None, None, :] - kpts[:, None, :] + kpts # Apply k shift kvKLM = kvKLM - kpts[kshift] for N, kvN in enumerate(kpts): kvKLMN = np.einsum('wx,klmx->wklm', a, kvKLM - kvN) # check whether (1/(2pi) k_{KLMN} dot a) is an integer kvKLMN_int = np.rint(kvKLMN) mask = np.einsum('wklm->klm', abs(kvKLMN - kvKLMN_int)) < 1e-9 kconserv_r2[mask] = N return kconserv_r2
Example #16
Source File: khf.py From pyscf with Apache License 2.0 | 5 votes |
def get_init_guess(self, cell=None, key='minao'): if cell is None: cell = self.cell dm_kpts = None key = key.lower() if key == '1e' or key == 'hcore': dm_kpts = self.init_guess_by_1e(cell) elif getattr(cell, 'natm', 0) == 0: logger.info(self, 'No atom found in cell. Use 1e initial guess') dm_kpts = self.init_guess_by_1e(cell) elif key == 'atom': dm = self.init_guess_by_atom(cell) elif key[:3] == 'chk': try: dm_kpts = self.from_chk() except (IOError, KeyError): logger.warn(self, 'Fail to read %s. Use MINAO initial guess', self.chkfile) dm = self.init_guess_by_minao(cell) else: dm = self.init_guess_by_minao(cell) if dm_kpts is None: dm_kpts = lib.asarray([dm]*len(self.kpts)) ne = np.einsum('kij,kji->', dm_kpts, self.get_ovlp(cell)).real # FIXME: consider the fractional num_electron or not? This maybe # relate to the charged system. nkpts = len(self.kpts) nelectron = float(self.cell.tot_electrons(nkpts)) if abs(ne - nelectron) > 1e-7*nkpts: logger.debug(self, 'Big error detected in the electron number ' 'of initial guess density matrix (Ne/cell = %g)!\n' ' This can cause huge error in Fock matrix and ' 'lead to instability in SCF for low-dimensional ' 'systems.\n DM is normalized wrt the number ' 'of electrons %s', ne/nkpts, nelectron/nkpts) dm_kpts *= (nelectron / ne).reshape(-1,1,1) return dm_kpts
Example #17
Source File: multigrid.py From pyscf with Apache License 2.0 | 5 votes |
def get_j_kpts(mydf, dm_kpts, hermi=1, kpts=numpy.zeros((1,3)), kpts_band=None): '''Get the Coulomb (J) AO matrix at sampled k-points. Args: dm_kpts : (nkpts, nao, nao) ndarray or a list of (nkpts,nao,nao) ndarray Density matrix at each k-point. If a list of k-point DMs, eg, UHF alpha and beta DM, the alpha and beta DMs are contracted separately. kpts : (nkpts, 3) ndarray Kwargs: kpts_band : (3,) ndarray or (*,3) ndarray A list of arbitrary "band" k-points at which to evalute the matrix. Returns: vj : (nkpts, nao, nao) ndarray or list of vj if the input dm_kpts is a list of DMs ''' cell = mydf.cell dm_kpts = numpy.asarray(dm_kpts) rhoG = _eval_rhoG(mydf, dm_kpts, hermi, kpts, deriv=0) coulG = tools.get_coulG(cell, mesh=cell.mesh) #:vG = numpy.einsum('ng,g->ng', rhoG[:,0], coulG) vG = rhoG[:,0] vG *= coulG kpts_band, input_band = _format_kpts_band(kpts_band, kpts), kpts_band vj_kpts = _get_j_pass2(mydf, vG, kpts_band) return _format_jks(vj_kpts, dm_kpts, input_band, kpts)
Example #18
Source File: eom_kccsd_ghf.py From pyscf with Apache License 2.0 | 5 votes |
def eaccsd_diag(eom, kshift, imds=None): if imds is None: imds = eom.make_imds() t1, t2 = imds.t1, imds.t2 nkpts, nocc, nvir = t1.shape kconserv = imds.kconserv Hr1 = np.diag(imds.Fvv[kshift]) Hr2 = np.zeros((nkpts,nkpts,nocc,nvir,nvir), dtype=t1.dtype) if eom.partition == 'mp': # This case is untested foo = eom.eris.fock[:,:nocc,:nocc] fvv = eom.eris.fock[:,nocc:,nocc:] for kj in range(nkpts): for ka in range(nkpts): kb = kconserv[kshift,ka,kj] Hr2[kj,ka] -= foo[kj].diagonal()[:,None,None] Hr2[kj,ka] -= fvv[ka].diagonal()[None,:,None] Hr2[kj,ka] += fvv[kb].diagonal()[None,None,:] else: for kj in range(nkpts): for ka in range(nkpts): kb = kconserv[kshift,ka,kj] Hr2[kj,ka] -= imds.Foo[kj].diagonal()[:,None,None] Hr2[kj,ka] += imds.Fvv[ka].diagonal()[None,:,None] Hr2[kj,ka] += imds.Fvv[kb].diagonal()[None,None,:] Hr2[kj,ka] += np.einsum('jbbj->jb', imds.Wovvo[kj,kb,kb])[:, None, :] Hr2[kj,ka] += np.einsum('jaaj->ja', imds.Wovvo[kj,ka,ka])[:, :, None] if ka == kconserv[ka,kb,kb]: Hr2[kj,ka] += np.einsum('abab->ab', imds.Wvvvv[ka,kb,ka])[None,:,:] Hr2[kj,ka] -= np.einsum('kjab,kjab->jab',imds.Woovv[kshift,kj,ka],imds.t2[kshift,kj,ka]) vector = amplitudes_to_vector_ea(Hr1, Hr2, kshift, kconserv) return vector
Example #19
Source File: krohf.py From pyscf with Apache License 2.0 | 5 votes |
def eig(self, fock, s): e, c = khf.KSCF.eig(self, fock, s) if getattr(fock, 'focka', None) is not None: for k, mo in enumerate(c): fa, fb = fock.focka[k], fock.fockb[k] mo_ea = np.einsum('pi,pi->i', mo.conj(), fa.dot(mo)).real mo_eb = np.einsum('pi,pi->i', mo.conj(), fb.dot(mo)).real e[k] = lib.tag_array(e[k], mo_ea=mo_ea, mo_eb=mo_eb) return e, c
Example #20
Source File: krohf.py From pyscf with Apache License 2.0 | 5 votes |
def get_init_guess(self, cell=None, key='minao'): if cell is None: cell = self.cell dm_kpts = None key = key.lower() if key == '1e' or key == 'hcore': dm_kpts = self.init_guess_by_1e(cell) elif getattr(cell, 'natm', 0) == 0: logger.info(self, 'No atom found in cell. Use 1e initial guess') dm_kpts = self.init_guess_by_1e(cell) elif key == 'atom': dm = self.init_guess_by_atom(cell) elif key[:3] == 'chk': try: dm_kpts = self.from_chk() except (IOError, KeyError): logger.warn(self, 'Fail to read %s. Use MINAO initial guess', self.chkfile) dm = self.init_guess_by_minao(cell) else: dm = self.init_guess_by_minao(cell) if dm_kpts is None: nkpts = len(self.kpts) # dm[spin,nao,nao] at gamma point -> dm_kpts[spin,nkpts,nao,nao] dm_kpts = np.repeat(dm[:,None,:,:], nkpts, axis=1) ne = np.einsum('xkij,kji->', dm_kpts, self.get_ovlp(cell)).real # FIXME: consider the fractional num_electron or not? This maybe # relates to the charged system. nkpts = len(self.kpts) nelec = float(sum(self.nelec)) if np.any(abs(ne - nelec) > 1e-7*nkpts): logger.debug(self, 'Big error detected in the electron number ' 'of initial guess density matrix (Ne/cell = %g)!\n' ' This can cause huge error in Fock matrix and ' 'lead to instability in SCF for low-dimensional ' 'systems.\n DM is normalized wrt the number ' 'of electrons %g', ne/nkpts, nelec/nkpts) dm_kpts *= nelec / ne return dm_kpts
Example #21
Source File: eom_kccsd_ghf.py From pyscf with Apache License 2.0 | 5 votes |
def eaccsd_matvec(eom, vector, kshift, imds=None, diag=None): '''2hp operators are of the form s_{ j}^{ab}, i.e. 'jb' indices are coupled.''' if imds is None: imds = eom.make_imds() nocc = eom.nocc nmo = eom.nmo nkpts = eom.nkpts kconserv = imds.kconserv r1, r2 = vector_to_amplitudes_ea(vector, kshift, nkpts, nmo, nocc, kconserv) Hr1 = np.einsum('ac,c->a', imds.Fvv[kshift], r1) for kl in range(nkpts): Hr1 += np.einsum('ld,lad->a', imds.Fov[kl], r2[kl, kshift]) for kc in range(nkpts): Hr1 += 0.5*np.einsum('alcd,lcd->a', imds.Wvovv[kshift,kl,kc], r2[kl,kc]) Hr2 = np.zeros_like(r2) for kj, ka in itertools.product(range(nkpts), repeat=2): kb = kconserv[kshift,ka,kj] Hr2[kj,ka] += np.einsum('abcj,c->jab', imds.Wvvvo[ka,kb,kshift], r1) Hr2[kj,ka] += lib.einsum('ac,jcb->jab', imds.Fvv[ka], r2[kj,ka]) Hr2[kj,ka] -= lib.einsum('bc,jca->jab', imds.Fvv[kb], r2[kj,kb]) Hr2[kj,ka] -= lib.einsum('lj,lab->jab', imds.Foo[kj], r2[kj,ka]) for kd in range(nkpts): kl = kconserv[kj, kb, kd] Hr2[kj, ka] += lib.einsum('lbdj,lad->jab', imds.Wovvo[kl, kb, kd], r2[kl, ka]) # P(ab) kl = kconserv[kj, ka, kd] Hr2[kj, ka] -= lib.einsum('ladj,lbd->jab', imds.Wovvo[kl, ka, kd], r2[kl, kb]) kc = kconserv[ka, kd, kb] Hr2[kj, ka] += 0.5 * lib.einsum('abcd,jcd->jab', imds.Wvvvv[ka, kb, kc], r2[kj, kc]) tmp = lib.einsum('xyklcd,xylcd->k', imds.Woovv[kshift, :, :], r2[:, :]) # contract_{kl, kc} Hr2[:, :] -= 0.5*lib.einsum('k,xykjab->xyjab', tmp, imds.t2[kshift, :, :]) # sum_{kj, ka] vector = eom.amplitudes_to_vector(Hr1, Hr2, kshift) return vector
Example #22
Source File: krohf.py From pyscf with Apache License 2.0 | 5 votes |
def canonicalize(mf, mo_coeff_kpts, mo_occ_kpts, fock=None): '''Canonicalization diagonalizes the ROHF Fock matrix within occupied, virtual subspaces separatedly (without change occupancy). ''' if fock is None: dm = mf.make_rdm1(mo_coeff_kpts, mo_occ_kpts) fock = mf.get_fock(dm=dm) mo_coeff = [] mo_energy = [] for k, mo in enumerate(mo_coeff_kpts): mo1 = np.empty_like(mo) mo_e = np.empty_like(mo_occ_kpts[k]) coreidx = mo_occ_kpts[k] == 2 openidx = mo_occ_kpts[k] == 1 viridx = mo_occ_kpts[k] == 0 for idx in (coreidx, openidx, viridx): if np.count_nonzero(idx) > 0: orb = mo[:,idx] f1 = reduce(np.dot, (orb.T.conj(), fock[k], orb)) e, c = scipy.linalg.eigh(f1) mo1[:,idx] = np.dot(orb, c) mo_e[idx] = e if getattr(fock, 'focka', None) is not None: fa, fb = fock.focka[k], fock.fockb[k] mo_ea = np.einsum('pi,pi->i', mo1.conj(), fa.dot(mo1)).real mo_eb = np.einsum('pi,pi->i', mo1.conj(), fb.dot(mo1)).real mo_e = lib.tag_array(mo_e, mo_ea=mo_ea, mo_eb=mo_eb) mo_coeff.append(mo1) mo_energy.append(mo_e) return mo_energy, mo_coeff
Example #23
Source File: eom_kccsd_ghf.py From pyscf with Apache License 2.0 | 5 votes |
def ipccsd_diag(eom, kshift, imds=None): if imds is None: imds = eom.make_imds() t1, t2 = imds.t1, imds.t2 nkpts, nocc, nvir = t1.shape kconserv = imds.kconserv Hr1 = -np.diag(imds.Foo[kshift]) Hr2 = np.zeros((nkpts,nkpts,nocc,nocc,nvir), dtype=t1.dtype) if eom.partition == 'mp': foo = eom.eris.fock[:,:nocc,:nocc] fvv = eom.eris.fock[:,nocc:,nocc:] for ki in range(nkpts): for kj in range(nkpts): ka = kconserv[ki,kshift,kj] Hr2[ki,kj] -= foo[ki].diagonal()[:,None,None] Hr2[ki,kj] -= foo[kj].diagonal()[None,:,None] Hr2[ki,kj] += fvv[ka].diagonal()[None,None,:] else: for ki in range(nkpts): for kj in range(nkpts): ka = kconserv[ki,kshift,kj] Hr2[ki,kj] -= imds.Foo[ki].diagonal()[:,None,None] Hr2[ki,kj] -= imds.Foo[kj].diagonal()[None,:,None] Hr2[ki,kj] += imds.Fvv[ka].diagonal()[None,None,:] if ki == kconserv[ki,kj,kj]: Hr2[ki,kj] += np.einsum('ijij->ij', imds.Woooo[ki, kj, ki])[:,:,None] Hr2[ki, kj] += lib.einsum('iaai->ia', imds.Wovvo[ki, ka, ka])[:,None,:] Hr2[ki, kj] += lib.einsum('jaaj->ja', imds.Wovvo[kj, ka, ka])[None,:,:] Hr2[ki, kj] += lib.einsum('ijea,jiea->ija',imds.Woovv[ki,kj,kshift], imds.t2[kj,ki,kshift]) vector = amplitudes_to_vector_ip(Hr1, Hr2, kshift, kconserv) return vector
Example #24
Source File: test_x2c.py From pyscf with Apache License 2.0 | 5 votes |
def test_hf(self): with lib.light_speed(2) as c: mf = scf.RHF(cell).sfx2c1e() mf.with_df = df.AFTDF(cell) dm = mf.get_init_guess() h1 = mf.get_hcore() self.assertAlmostEqual(numpy.einsum('ij,ji', dm, h1), -0.2458227312351979+0j, 8) kpts = cell.make_kpts([3,1,1]) h1 = mf.get_hcore(kpt=kpts[1]) self.assertAlmostEqual(numpy.einsum('ij,ji', dm, h1), -0.04113247191600125+0j, 8)
Example #25
Source File: sfx2c1e.py From pyscf with Apache License 2.0 | 5 votes |
def picture_change(self, even_operator=(None, None), odd_operator=None): '''Picture change for even_operator + odd_operator even_operator has two terms at diagonal blocks [ v 0 ] [ 0 w ] odd_operator has the term at off-diagonal blocks [ 0 p ] [ p^T 0 ] v, w, and p can be strings (integral name) or matrices. ''' mol = self.mol xmol, c = self.get_xmol(mol) pc_mat = self._picture_change(xmol, even_operator, odd_operator) if self.basis is not None: s22 = xmol.intor_symmetric('int1e_ovlp') s21 = gto.mole.intor_cross('int1e_ovlp', xmol, mol) c = lib.cho_solve(s22, s21) elif self.xuncontract: pass else: return pc_mat if pc_mat.ndim == 2: return lib.einsum('pi,pq,qj->ij', c, pc_mat, c) else: return lib.einsum('pi,xpq,qj->xij', c, pc_mat, c)
Example #26
Source File: test_x2c_grad.py From pyscf with Apache License 2.0 | 5 votes |
def _invsqrt1(a0, a1): '''Solving first order of x^2 = a^{-1}''' w, v = scipy.linalg.eigh(a0) w = 1./numpy.sqrt(w) a1 = -reduce(numpy.dot, (v.conj().T, a1, v)) x1 = numpy.einsum('i,ij,j->ij', w**2, a1, w**2) / (w[:,None] + w) x1 = reduce(numpy.dot, (v, x1, v.conj().T)) return x1
Example #27
Source File: test_x2c_hess.py From pyscf with Apache License 2.0 | 5 votes |
def _invsqrt1(a0, a1): '''Solving first order derivative of x^2 = a^{-1}''' w, v = scipy.linalg.eigh(a0) w = 1./numpy.sqrt(w) a1 = -reduce(numpy.dot, (v.conj().T, a1, v)) x1 = numpy.einsum('i,ij,j->ij', w**2, a1, w**2) / (w[:,None] + w) x1 = reduce(numpy.dot, (v, x1, v.conj().T)) return x1
Example #28
Source File: sfx2c1e_grad.py From pyscf with Apache License 2.0 | 5 votes |
def _get_r1(s0_roots, s_nesc0, s1, s_nesc1, r0_roots): # See JCP 135, 084114 (2011); DOI:10.1063/1.3624397, Eq (34) w_sqrt, v_s = s0_roots w_invsqrt = 1. / w_sqrt wr0_sqrt, vr0 = r0_roots wr0_invsqrt = 1. / wr0_sqrt s1 = reduce(numpy.dot, (v_s.T, s1, v_s)) s_nesc1 = reduce(numpy.dot, (v_s.T, s_nesc1, v_s)) s1_sqrt = s1 / (w_sqrt[:,None] + w_sqrt) s1_invsqrt = (numpy.einsum('i,ij,j->ij', w_invsqrt**2, s1, w_invsqrt**2) / -(w_invsqrt[:,None] + w_invsqrt)) R1_mid = numpy.dot(s1_invsqrt, s_nesc0) * w_invsqrt R1_mid = R1_mid + R1_mid.T R1_mid += numpy.einsum('i,ij,j->ij', w_invsqrt, s_nesc1, w_invsqrt) R1_mid = reduce(numpy.dot, (vr0.T, R1_mid, vr0)) R1_mid /= -(wr0_invsqrt[:,None] + wr0_invsqrt) R1_mid = numpy.einsum('i,ij,j->ij', wr0_invsqrt**2, R1_mid, wr0_invsqrt**2) vr0_wr0_sqrt = vr0 * wr0_invsqrt vr0_s0_sqrt = vr0.T * w_sqrt vr0_s0_invsqrt = vr0.T * w_invsqrt R1 = reduce(numpy.dot, (vr0_s0_invsqrt.T, R1_mid, vr0_s0_sqrt)) R1 += reduce(numpy.dot, (s1_invsqrt, vr0_wr0_sqrt, vr0_s0_sqrt)) R1 += reduce(numpy.dot, (vr0_s0_invsqrt.T, vr0_wr0_sqrt.T, s1_sqrt)) R1 = reduce(numpy.dot, (v_s, R1, v_s.T)) return R1
Example #29
Source File: sfx2c1e_grad.py From pyscf with Apache License 2.0 | 5 votes |
def gen_sf_hfw(mol, approx='1E'): approx = approx.upper() c = lib.param.LIGHT_SPEED h0, s0 = _get_h0_s0(mol) e0, c0 = scipy.linalg.eigh(h0, s0) aoslices = mol.aoslice_by_atom() nao = mol.nao_nr() if 'ATOM' in approx: x0 = numpy.zeros((nao,nao)) for ia in range(mol.natm): ish0, ish1, p0, p1 = aoslices[ia] shls_slice = (ish0, ish1, ish0, ish1) t1 = mol.intor('int1e_kin', shls_slice=shls_slice) s1 = mol.intor('int1e_ovlp', shls_slice=shls_slice) with mol.with_rinv_at_nucleus(ia): z = -mol.atom_charge(ia) v1 = z * mol.intor('int1e_rinv', shls_slice=shls_slice) w1 = z * mol.intor('int1e_prinvp', shls_slice=shls_slice) x0[p0:p1,p0:p1] = x2c._x2c1e_xmatrix(t1, v1, w1, s1, c) else: cl0 = c0[:nao,nao:] cs0 = c0[nao:,nao:] x0 = scipy.linalg.solve(cl0.T, cs0.T).T s_nesc0 = s0[:nao,:nao] + reduce(numpy.dot, (x0.T, s0[nao:,nao:], x0)) R0 = x2c._get_r(s0[:nao,:nao], s_nesc0) c_fw0 = numpy.vstack((R0, numpy.dot(x0, R0))) h0_fw_half = numpy.dot(h0, c_fw0) get_h1_etc = _gen_first_order_quantities(mol, e0, c0, x0, approx) def hcore_deriv(ia): h1_ao, s1_ao, e1, c1, x1, s_nesc1, R1, c_fw1 = get_h1_etc(ia) hfw1 = lib.einsum('xpi,pj->xij', c_fw1, h0_fw_half) hfw1 = hfw1 + hfw1.transpose(0,2,1) hfw1+= lib.einsum('pi,xpq,qj->xij', c_fw0, h1_ao, c_fw0) return hfw1 return hcore_deriv
Example #30
Source File: molecular_example_odd_qubits.py From OpenFermion-Cirq with Apache License 2.0 | 5 votes |
def make_h3_2_5() -> Tuple[RestrictedHartreeFockObjective, of.MolecularData, np. ndarray, np.ndarray, np.ndarray]: # load the molecule from moelcular data h3_2_5_path = os.path.join( hfvqe.__path__[0], 'molecular_data/hydrogen_chains/h_3_p_sto-3g/bond_distance_2.5') molfile = os.path.join(h3_2_5_path, 'H3_plus_sto-3g_singlet_linear_r-2.5.hdf5') molecule = of.MolecularData(filename=molfile) molecule.load() S = np.load(os.path.join(h3_2_5_path, 'overlap.npy')) Hcore = np.load(os.path.join(h3_2_5_path, 'h_core.npy')) TEI = np.load(os.path.join(h3_2_5_path, 'tei.npy')) _, X = sp.linalg.eigh(Hcore, S) obi = of.general_basis_change(Hcore, X, (1, 0)) tbi = np.einsum('psqr', of.general_basis_change(TEI, X, (1, 0, 1, 0))) molecular_hamiltonian = generate_hamiltonian(obi, tbi, molecule.nuclear_repulsion) rhf_objective = RestrictedHartreeFockObjective(molecular_hamiltonian, molecule.n_electrons) scipy_result = rhf_minimization(rhf_objective) return rhf_objective, molecule, scipy_result.x, obi, tbi