diff --git a/.flake8 b/.flake8 index ed8ac06aa1..3e1f237f91 100644 --- a/.flake8 +++ b/.flake8 @@ -12,8 +12,8 @@ ignore = # Imports: E401, E402, # Other: - E701, E731, E741, E275, - F401, C901, W391, W503, W504 + E701, E713, E721, E731, E741, E275, + F401, F403, C901, W391, W503, W504 exclude = test, .git, __pycache__, build, dist, __init__.py, .eggs, *.egg max-line-length = 120 diff --git a/.github/workflows/lint.yml b/.github/workflows/lint.yml index 755e768b41..a3153bb291 100644 --- a/.github/workflows/lint.yml +++ b/.github/workflows/lint.yml @@ -3,6 +3,17 @@ name: Lint on: [push, pull_request] jobs: + ruff: + runs-on: ubuntu-latest + steps: + - name: Checkout + uses: actions/checkout@v4 + - name: Install ruff + run: pip install ruff + - name: Check style + run: ruff check --config .ruff.toml pyscf + - name: Check NumPy + run: ruff check --select NPY --ignore NPY002 pyscf flake: runs-on: ubuntu-latest steps: diff --git a/.ruff.toml b/.ruff.toml new file mode 100644 index 0000000000..b49fb70faf --- /dev/null +++ b/.ruff.toml @@ -0,0 +1,32 @@ +exclude = [ + "test", +] + +respect-gitignore = true +line-length = 120 + +[lint] +# Enable Pyflakes (`F`) and a subset of the pycodestyle (`E`) codes by default. +#select = ["E4", "E7", "E9", "F"] +ignore = [ +# Whitespaces: +"E201", "E202", "E203", "E211", "E221", "E222", "E225", "E226", "E228", "E231", "E241", "E251", +# Comments: +"E261", "E262", "E265", "E266", +# Blank lines: +"E301", "E302", "E303", "E305", "E306", +# Imports: +"E401", "E402", +# Other: +"E701", "E713", "E721", "E731", "E741", "E275", +"F401", "F403", "C901", "W391" +] + +[format] +# Like Black, indent with spaces, rather than tabs. +indent-style = "space" + +# Like Black, automatically detect the appropriate line ending. +line-ending = "auto" + +quote-style = "single" diff --git a/examples/1-advanced/019-state_average_dmrg_casscf.py b/examples/1-advanced/019-state_average_dmrg_casscf.py deleted file mode 120000 index e742614c54..0000000000 --- a/examples/1-advanced/019-state_average_dmrg_casscf.py +++ /dev/null @@ -1 +0,0 @@ -../dmrg/10-state_average.py \ No newline at end of file diff --git a/examples/1-advanced/020-dmrg_casscf_nevpt2_for_FeS.py b/examples/1-advanced/020-dmrg_casscf_nevpt2_for_FeS.py deleted file mode 120000 index 4e369af616..0000000000 --- a/examples/1-advanced/020-dmrg_casscf_nevpt2_for_FeS.py +++ /dev/null @@ -1 +0,0 @@ -../dmrg/32-dmrg_casscf_nevpt2_for_FeS.py \ No newline at end of file diff --git a/examples/local_orb/nlocal.py b/examples/local_orb/nlocal.py deleted file mode 100644 index 04b9b78580..0000000000 --- a/examples/local_orb/nlocal.py +++ /dev/null @@ -1,213 +0,0 @@ -# Copyright 2014-2018 The PySCF Developers. All Rights Reserved. -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -# -# Localization based on UNO from UHF/UKS check files -# -from functools import reduce -import numpy -import scipy.linalg -import h5py -from pyscf import tools,gto,scf,dft -from pyscf.tools import molden -import pmloc -import ulocal - -def sqrtm(s): - e, v = numpy.linalg.eigh(s) - return numpy.dot(v*numpy.sqrt(e), v.T.conj()) - -def lowdin(s): - e, v = numpy.linalg.eigh(s) - return numpy.dot(v/numpy.sqrt(e), v.T.conj()) - -def dumpLUNO(fname,thresh=0.01): - chkfile = fname+'.chk' - outfile = fname+'_cmo.molden' - tools.molden.from_chkfile(outfile, chkfile) - #============================= - # Natural orbitals - # Lowdin basis X=S{-1/2} - # psi = chi * C - # = chi' * C' - # = chi*X*(X{-1}C') - #============================= - mol,mf = scf.chkfile.load_scf(chkfile) - mo_coeff = mf["mo_coeff"] - ova=mol.intor_symmetric("cint1e_ovlp_sph") - nb = mo_coeff.shape[1] - # Check overlap - diff = reduce(numpy.dot,(mo_coeff[0].T,ova,mo_coeff[0])) - numpy.identity(nb) - print numpy.linalg.norm(diff) - diff = reduce(numpy.dot,(mo_coeff[1].T,ova,mo_coeff[1])) - numpy.identity(nb) - print numpy.linalg.norm(diff) - # UHF-alpha/beta - ma = mo_coeff[0] - mb = mo_coeff[1] - nalpha = (mol.nelectron+mol.spin)/2 - nbeta = (mol.nelectron-mol.spin)/2 - # Spin-averaged DM - pTa = numpy.dot(ma[:,:nalpha],ma[:,:nalpha].T) - pTb = numpy.dot(mb[:,:nbeta],mb[:,:nbeta].T) - pT = 0.5*(pTa+pTb) - # Lowdin basis - s12 = sqrtm(ova) - s12inv = lowdin(ova) - pTOAO = reduce(numpy.dot,(s12,pT,s12)) - eig,coeff = scipy.linalg.eigh(-pTOAO) - eig = -2.0*eig - eig[eig<0.0]=0.0 - eig[abs(eig)<1.e-14]=0.0 - ifplot = False #True - if ifplot: - import matplotlib.pyplot as plt - plt.plot(range(nb),eig,'ro') - plt.show() - # Back to AO basis - coeff = numpy.dot(s12inv,coeff) - diff = reduce(numpy.dot,(coeff.T,ova,coeff)) - numpy.identity(nb) - print 'CtSC-I',numpy.linalg.norm(diff) - # - # Averaged Fock - # - enorb = mf["mo_energy"] - fa = reduce(numpy.dot,(ma,numpy.diag(enorb[0]),ma.T)) - fb = reduce(numpy.dot,(mb,numpy.diag(enorb[1]),mb.T)) - # Non-orthogonal cases: FC=SCE - # Fao = SC*e*C{-1} = S*C*e*Ct*S - fav = 0.5*(fa+fb) - # Expectation value of natural orbitals - fexpt = reduce(numpy.dot,(coeff.T,ova,fav,ova,coeff)) - enorb = numpy.diag(fexpt) - nocc = eig.copy() - # - # Reordering and define active space according to thresh - # - idx = 0 - active=[] - for i in range(nb): - if nocc[i]<=2.0-thresh and nocc[i]>=thresh: - active.append(True) - else: - active.append(False) - print '\nNatural orbitals:' - for i in range(nb): - print 'orb:',i,active[i],nocc[i],enorb[i] - active = numpy.array(active) - actIndices = list(numpy.argwhere(active==True).flatten()) - cOrbs = coeff[:,:actIndices[0]] - aOrbs = coeff[:,actIndices] - vOrbs = coeff[:,actIndices[-1]+1:] - nb = cOrbs.shape[0] - nc = cOrbs.shape[1] - na = aOrbs.shape[1] - nv = vOrbs.shape[1] - print 'core orbs:',cOrbs.shape - print 'act orbs:',aOrbs.shape - print 'vir orbs:',vOrbs.shape - assert nc+na+nv == nb - # dump UNO - with open(fname+'_uno.molden','w') as thefile: - molden.header(mol,thefile) - molden.orbital_coeff(mol,thefile,coeff) - #===================== - # Population analysis - #===================== - from pyscf import lo - aux = lo.orth_ao(mol,method='meta_lowdin') - #clmo = ulocal.scdm(cOrbs,ova,aux) - #almo = ulocal.scdm(aOrbs,ova,aux) - clmo = cOrbs - almo = aOrbs - ierr,uc = pmloc.loc(mol,clmo) - ierr,ua = pmloc.loc(mol,almo) - clmo = clmo.dot(uc) - almo = almo.dot(ua) - vlmo = ulocal.scdm(vOrbs,ova,aux) - # P-SORT - mo_c,n_c,e_c = ulocal.psort(ova,fav,pT,clmo) - mo_o,n_o,e_o = ulocal.psort(ova,fav,pT,almo) - mo_v,n_v,e_v = ulocal.psort(ova,fav,pT,vlmo) - lmo = numpy.hstack((mo_c,mo_o,mo_v)).copy() - enorb = numpy.hstack([e_c,e_o,e_v]) - occ = numpy.hstack([n_c,n_o,n_v]) - # CHECK - diff = reduce(numpy.dot,(lmo.T,ova,lmo)) - numpy.identity(nb) - print 'diff=',numpy.linalg.norm(diff) - ulocal.lowdinPop(mol,lmo,ova,enorb,occ) - ulocal.dumpLMO(mol,fname,lmo) - print 'nalpha,nbeta,mol.spin,nb:',\ - nalpha,nbeta,mol.spin,nb - return mol,ova,fav,pT,nb,nalpha,nbeta,nc,na,nv,lmo,enorb,occ - -def dumpAct(fname,info,actlst,base=1): - actlst2 = [i-base for i in actlst] - mol,ova,fav,pT,nb,nalpha,nbeta,nc,na,nv,lmo,enorb,occ = info - corb = set(range(nc)) - aorb = set(range(nc,nc+na)) - vorb = set(range(nc+na,nc+na+nv)) - print '[dumpAct]' - print ' corb=',corb - print ' aorb=',aorb - print ' vorb=',vorb - sorb = set(actlst2) - rcorb = corb.difference(corb.intersection(sorb)) - #assuming act in actlst - #raorb = aorb.difference(aorb.intersection(sorb)) - rvorb = vorb.difference(vorb.intersection(sorb)) - corb = list(rcorb) - aorb = list(sorb) - vorb = list(rvorb) - print ' corb=',corb - print ' aorb=',aorb - print ' vorb=',vorb - clmo = lmo[:,corb].copy() - almo = lmo[:,aorb].copy() - vlmo = lmo[:,vorb].copy() - ierr,ua = pmloc.loc(mol,almo) - almo = almo.dot(ua) - #>>> DUMP <<<# - # P-SORT - mo_c = clmo - mo_v = vlmo - e_c = enorb[corb].copy() - e_v = enorb[vorb].copy() - n_c = occ[corb].copy() - n_v = occ[vorb].copy() - mo_o,n_o,e_o = ulocal.psort(ova,fav,pT,almo) - lmo2 = numpy.hstack((mo_c,mo_o,mo_v)) - enorb = numpy.hstack([e_c,e_o,e_v]) - occ = numpy.hstack([n_c,n_o,n_v]) - assert len(enorb)==nb - assert len(occ)==nb - # CHECK - diff = reduce(numpy.dot,(lmo2.T,ova,lmo2)) - numpy.identity(nb) - print 'diff=',numpy.linalg.norm(diff) - ulocal.lowdinPop(mol,lmo,ova,enorb,occ) - ulocal.dumpLMO(mol,fname+'_new',lmo2) - print 'nalpha,nbeta,mol.spin,nb:',\ - nalpha,nbeta,mol.spin,nb - print 'diff(LMO2-LMO)=',numpy.linalg.norm(lmo2-lmo) - nc = len(e_c) - na = len(e_o) - nv = len(e_v) - assert na == len(actlst) - assert nc+na+nv == nb - print 'nc,na,nv,nb=',nc,na,nv,nb - return lmo2,nc,na,nv - -if __name__ == '__main__': - fname = 'hs_bp86' - info = dumpLUNO(fname) - actlst = [117,118,119,120,125,126]+range(127,137) - dumpAct(fname,info,actlst,base=1) diff --git a/examples/local_orb/pmloc.py b/examples/local_orb/pmloc.py deleted file mode 100644 index 6a9626e796..0000000000 --- a/examples/local_orb/pmloc.py +++ /dev/null @@ -1,286 +0,0 @@ -# Copyright 2014-2018 The PySCF Developers. All Rights Reserved. -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -# -# Localization function: pmloc.loc(mol,mocoeff) -# -# -ZL@20151227- Add Boys. Note that to view LMOs -# in jmol, 'orbital energies' need -# to be defined and used to reorder -# LMOs. -# -ZL@20151225- Pipek-Mezey localizaiton -# -from functools import reduce -import math -import numpy -import scipy.linalg -from pyscf.tools import molden - -#------------------------------------------------ -# Initialization -#------------------------------------------------ -def sqrtm(s): - e, v = numpy.linalg.eigh(s) - return numpy.dot(v*numpy.sqrt(e), v.T.conj()) - -def lowdin(s): - e, v = numpy.linalg.eigh(s) - return numpy.dot(v/numpy.sqrt(e), v.T.conj()) - -#def scdmU(coeff,ova): -# aux = numpy.identity(ova.shape[0]) -# #aux = lowdin(ova) -# no = coeff.shape[1] -# ova = reduce(numpy.dot,(coeff.T,ova,aux)) -# # ova = no*nb -# q,r,piv = scipy.linalg.qr(ova, pivoting=True) -# # In fact, it is just "Lowdin-orthonormalized PAO". -# bc = ova[:,piv[:no]] -# ova = numpy.dot(bc.T,bc) -# s12inv = lowdin(ova) -# u = numpy.dot(bc,s12inv) -# return u - -#------------------------------------------------ -# Boys/PM-Localization -#------------------------------------------------ -def loc(mol,mocoeff,tol=1.e-6,maxcycle=1000,iop=0): - partition = genAtomPartition(mol) - ova = mol.intor_symmetric("cint1e_ovlp_sph") - ierr,u = loc_kernel(mol,mocoeff,ova,partition,tol,maxcycle,iop) - #ierr = 0 - #u = scdmU(mocoeff,ova) - #if iop <= 2: - # mocoeff2 = mocoeff.dot(u) - # ierr,u2 = loc_kernel(mol,mocoeff2,ova,partition,tol,maxcycle,iop) - # u = u.dot(u2) - return ierr,u - -def loc_kernel(mol,mocoeff,ova,partition,tol,maxcycle,iop): - debug = False - print - print '[pm_loc_kernel]' - print ' mocoeff.shape=',mocoeff.shape - print ' tol=',tol - print ' maxcycle=',maxcycle - print ' partition=',len(partition),'\n',partition - k = mocoeff.shape[0] - n = mocoeff.shape[1] - natom = len(partition) - - def genPaij(mol,mocoeff,ova,partition,iop): - c = mocoeff.copy() - # Mulliken matrix - if iop == 0: - cts = c.T.dot(ova) - natom = len(partition) - pija = numpy.zeros((natom,n,n)) - for iatom in range(natom): - idx = partition[iatom] - tmp = numpy.dot(cts[:,idx],c[idx,:]) - pija[iatom] = 0.5*(tmp+tmp.T) - # Lowdin - elif iop == 1: - s12 = sqrtm(ova) - s12c = s12.T.dot(c) - natom = len(partition) - pija = numpy.zeros((natom,n,n)) - for iatom in range(natom): - idx = partition[iatom] - pija[iatom] = numpy.dot(s12c[idx,:].T,s12c[idx,:]) - # Boys - elif iop == 2: - rmat = mol.intor_symmetric('cint1e_r_sph',3) - pija = numpy.zeros((3,n,n)) - for icart in range(3): - pija[icart] = reduce(numpy.dot,(c.T,rmat[icart],c)) - # P[i,j,a] - pija = pija.transpose(1,2,0).copy() - return pija - - ## Initial from random unitary - #iden = numpy.identity(n) - #iden += 1.e-2*numpy.random.uniform(-1,1,size=n*n).reshape(n,n) - #u,r = scipy.linalg.qr(iden) - #mou = mocoeff.dot(u) - #pija = genPaij(mol,mou,ova,partition,iop) - u = numpy.identity(n) - pija = genPaij(mol,mocoeff,ova,partition,iop) - if debug: pija0 = pija.copy() - - # Start - def funval(pija): - return numpy.einsum('iia,iia',pija,pija) - - fun = funval(pija) - print ' initial funval = ',fun - # - # Iteration - # - for icycle in range(maxcycle): - delta = 0.0 - # i>j - ijdx = [] - for i in range(n-1): - for j in range(i+1,n): - bij = abs(numpy.sum(pija[i,j]*(pija[i,i]-pija[j,j]))) - ijdx.append((i,j,bij)) - # Without pivoting - # 6-31g: icycle= 184 delta= 5.6152945523e-09 fun= 54.4719738182 - # avtz : icycle= 227 delta= 4.43639097128e-09 fun= 3907.60402435 - # With pivoting - significantly accelerated! - # 6-31g: icycle= 71 delta= 7.3566217445e-09 fun= 54.4719739144 - # avdz : icycle= 28 delta= 2.31739493914e-10 fun= 3907.60402153 - # The delta value generally decay monotonically (adjoint diagonalization) - ijdx = sorted(ijdx,key=lambda x:x[2],reverse=True) - for i,j,bij in ijdx: - # - # determine angle - # - vij = pija[i,i]-pija[j,j] - aij = numpy.dot(pija[i,j],pija[i,j]) \ - - 0.25*numpy.dot(vij,vij) - bij = numpy.dot(pija[i,j],vij) - if abs(aij)<1.e-10 and abs(bij)<1.e-10: continue - p1 = math.sqrt(aij**2+bij**2) - cos4a = -aij/p1 - sin4a = bij/p1 - cos2a = math.sqrt((1+cos4a)*0.5) - sin2a = math.sqrt((1-cos4a)*0.5) - cosa = math.sqrt((1+cos2a)*0.5) - sina = math.sqrt((1-cos2a)*0.5) - # Why? Because we require alpha in [0,pi/2] - if sin4a < 0.0: - cos2a = -cos2a - sina,cosa = cosa,sina - # stationary condition - if abs(cosa-1.0)<1.e-10: continue - if abs(sina-1.0)<1.e-10: continue - # incremental value - delta += p1*(1-cos4a) - # - # Transformation - # - if debug: - g = numpy.identity(n) - g[i,i] = cosa - g[j,j] = cosa - g[i,j] = -sina - g[j,i] = sina - ug = u.dot(g) - pijag = numpy.einsum('ik,jl,ija->kla',ug,ug,pija0) - # Urot - ui = u[:,i]*cosa+u[:,j]*sina - uj = -u[:,i]*sina+u[:,j]*cosa - u[:,i] = ui.copy() - u[:,j] = uj.copy() - # Bra-transformation of Integrals - tmp_ip = pija[i,:,:]*cosa+pija[j,:,:]*sina - tmp_jp = -pija[i,:,:]*sina+pija[j,:,:]*cosa - pija[i,:,:] = tmp_ip.copy() - pija[j,:,:] = tmp_jp.copy() - # Ket-transformation of Integrals - tmp_ip = pija[:,i,:]*cosa+pija[:,j,:]*sina - tmp_jp = -pija[:,i,:]*sina+pija[:,j,:]*cosa - pija[:,i,:] = tmp_ip.copy() - pija[:,j,:] = tmp_jp.copy() - if debug: - diff1 = numpy.linalg.norm(u-ug) - diff2 = numpy.linalg.norm(pija-pijag) - cu = numpy.dot(mocoeff,u) - pija2 = genPaij(cu,ova,partition) - fun2 = funval(pija2) - diff = abs(fun+delta-fun2) - print 'diff(u/p/f)=',diff1,diff2,diff - if diff1>1.e-6: - print 'Error: ug',diff1 - exit() - if diff2>1.e-6: - print 'Error: pijag',diff2 - exit() - if diff>1.e-6: - print 'Error: inconsistency in PMloc: fun/fun2=',fun+delta,fun2 - exit() - - fun = fun+delta - print 'icycle=',icycle,'delta=',delta,'fun=',fun - if delta -def psort(ova,fav,pT,coeff): - # Compute expectation value - pTnew = 2.0*reduce(numpy.dot,(coeff.T,ova,pT,ova,coeff)) - nocc = numpy.diag(pTnew) - index = numpy.argsort(-nocc) - ncoeff = coeff[:,index] - nocc = nocc[index] - enorb = numpy.diag(reduce(numpy.dot,(coeff.T,ova,fav,ova,coeff))) - enorb = enorb[index] - return ncoeff,nocc,enorb - -def lowdinPop(mol,coeff,ova,enorb,occ): - print '\nLowdin population for LMOs:' - nb,nc = coeff.shape - s12 = sqrtm(ova) - lcoeff = s12.dot(coeff) - diff = reduce(numpy.dot,(lcoeff.T,lcoeff)) - numpy.identity(nc) - print 'diff=',numpy.linalg.norm(diff) - pthresh = 0.05 - labels = mol.ao_labels(None) - nelec = 0.0 - for iorb in range(nc): - vec = lcoeff[:,iorb]**2 - idx = list(numpy.argwhere(vec>pthresh)) - print ' iorb=',iorb,' occ=',occ[iorb],' =',enorb[iorb] - for iao in idx: - print ' iao=',labels[iao],' pop=',vec[iao] - nelec += occ[iorb] - print 'nelec=',nelec - return 0 - -def scdm(coeff,ova,aux): - no = coeff.shape[1] - ova = reduce(numpy.dot,(coeff.T,ova,aux)) - # ova = no*nb - q,r,piv = scipy.linalg.qr(ova,pivoting=True) - bc = ova[:,piv[:no]] - ova2 = numpy.dot(bc.T,bc) - s12inv = lowdin(ova2) - cnew = reduce(numpy.dot,(coeff,bc,s12inv)) - return cnew - -def dumpLMO(mol,fname,lmo): - print 'Dump into '+fname+'.h5' - f = h5py.File(fname+'.h5','w') - f.create_dataset("lmo",data=lmo) - f.close() - print 'Dump into '+fname+'_lmo.molden' - with open(fname+'_lmo.molden','w') as thefile: - molden.header(mol,thefile) - molden.orbital_coeff(mol,thefile,lmo) - return 0 - -#============================= -# DUMP from chkfile to molden -#============================= -def dumpLocal(fname): - chkfile = fname+'.chk' - outfile = fname+'_cmo.molden' - tools.molden.from_chkfile(outfile, chkfile) - - mol,mf = scf.chkfile.load_scf(chkfile) - mo_coeff = mf["mo_coeff"] - ova=mol.intor_symmetric("cint1e_ovlp_sph") - nb = mo_coeff.shape[1] - nalpha = (mol.nelectron+mol.spin)/2 - nbeta = (mol.nelectron-mol.spin)/2 - print 'nalpha,nbeta,mol.spin,nb:',\ - nalpha,nbeta,mol.spin,nb - # UHF-alpha/beta - ma = mo_coeff[0] - mb = mo_coeff[1] - - #============================= - # Localization - #============================= - ma_c = ma[:,:nalpha].copy() - ma_v = ma[:,nalpha:].copy() - #-------------------- - # Occupied space: PM - #-------------------- - import pmloc - ierr,uc = pmloc.loc(mol,ma_c) - mc = numpy.dot(ma_c,uc) - #-------------------- - # Virtual space: PAO - #-------------------- - from pyscf import lo - aux = lo.orth_ao(mol,method='meta_lowdin') - mv = scdm(ma_v,ova,aux) - - # P[dm] - pa = numpy.dot(ma[:,:nalpha],ma[:,:nalpha].T) - pb = numpy.dot(mb[:,:nbeta],mb[:,:nbeta].T) - pT = 0.5*(pa+pb) - # E-SORT - enorb = mf["mo_energy"] - fa = reduce(numpy.dot,(ma,numpy.diag(enorb[0]),ma.T)) - fb = reduce(numpy.dot,(mb,numpy.diag(enorb[1]),mb.T)) - fav = 0.5*(fa+fb) - mc,occ_c,ec = psort(ova,fav,pT,mc) - mv,occ_v,ev = psort(ova,fav,pT,mv) - #---Check--- - tij = reduce(numpy.dot,(mc.T,ova,ma_c)) - sig = scipy.linalg.svd(tij,compute_uv=False) - print 'nc=',nalpha,numpy.sum(sig**2) - assert abs(nalpha-numpy.sum(sig**2))<1.e-8 - tij = reduce(numpy.dot,(mv.T,ova,ma_v)) - sig = scipy.linalg.svd(tij,compute_uv=False) - print 'nv=',nb-nalpha,numpy.sum(sig**2) - assert abs(nb-nalpha-numpy.sum(sig**2))<1.e-8 - - lmo = numpy.hstack([mc,mv]) - enorb = numpy.hstack([ec,ev]) - occ = numpy.hstack([occ_c,occ_v]) - lowdinPop(mol,lmo,ova,enorb,occ) - dumpLMO(mol,fname,lmo) - print 'nalpha,nbeta,mol.spin,nb:',\ - nalpha,nbeta,mol.spin,nb - return 0 - -if __name__ == '__main__': - fname = 'hs_bp86' - dumpLocal(fname) diff --git a/examples/mcscf/33-make_init_guess/32-dmrg_casscf_nevpt2_for_FeS.py b/examples/mcscf/33-make_init_guess/32-dmrg_casscf_nevpt2_for_FeS.py deleted file mode 120000 index 59b169a4b3..0000000000 --- a/examples/mcscf/33-make_init_guess/32-dmrg_casscf_nevpt2_for_FeS.py +++ /dev/null @@ -1 +0,0 @@ -../../dmrg/32-dmrg_casscf_nevpt2_for_FeS.py \ No newline at end of file diff --git a/examples/mcscf/43-dmet_cas.py b/examples/mcscf/43-dmet_cas.py index 8c9e6011c1..8b48c4ecc5 100644 --- a/examples/mcscf/43-dmet_cas.py +++ b/examples/mcscf/43-dmet_cas.py @@ -6,6 +6,7 @@ from pyscf import scf from pyscf import gto from pyscf import mcscf, fci +from pyscf.mcscf import dmet_cas ''' Triplet and quintet energy gap of Iron-Porphyrin molecule @@ -115,7 +116,7 @@ mf = scf.fast_newton(mf) idx3d4d = mol.search_ao_label(['Fe 3d', 'Fe 4d']) ncas, nelecas, mo = dmet_cas.guess_cas(mf, mf.make_rdm1(), idx3d) -mc = mcscf.approx_hessian(mcscf.CASSCF(mf, ncas, nelecas) +mc = mcscf.approx_hessian(mcscf.CASSCF(mf, ncas, nelecas)) mc.kernel(mo) e_t = mc.e_tot # -2244.81493852189 diff --git a/examples/mrpt/03-dmrg_nevpt2.py b/examples/mrpt/03-dmrg_nevpt2.py deleted file mode 120000 index 06a4f8cf72..0000000000 --- a/examples/mrpt/03-dmrg_nevpt2.py +++ /dev/null @@ -1 +0,0 @@ -../../pyscf/dmrgscf/example/nevpt2.py \ No newline at end of file diff --git a/examples/pbc/44-k_point_eom_ccsd_ghf.py b/examples/pbc/44-k_point_eom_ccsd_ghf.py index 22887e543d..beeefb7bf6 100644 --- a/examples/pbc/44-k_point_eom_ccsd_ghf.py +++ b/examples/pbc/44-k_point_eom_ccsd_ghf.py @@ -96,11 +96,11 @@ # Molecular EOM-GCCSD myeomee = mol_eom_gccsd.EOMEE(mycc) -eee_mol, vee_mol = myeomee.kernel(nroots=nroots_test*np.product(nmp)) +eee_mol, vee_mol = myeomee.kernel(nroots=nroots_test*np.prod(nmp)) print("PBC KRHF Energy:", ekrhf) print("PBC RHF Energy :", erhf) print("PBC KGCCSD Energy :", ekgcc) -print("Mol GCCSD Energy per cell:", egcc / np.product(nmp)) +print("Mol GCCSD Energy per cell:", egcc / np.prod(nmp)) print("PBC EOMEE roots:", eee) print("Mol EOMEE roots:", eee_mol) diff --git a/examples/pbc/48-k_point_adc.py b/examples/pbc/48-k_point_adc.py index b19d43ad1d..6bf397ac0d 100644 --- a/examples/pbc/48-k_point_adc.py +++ b/examples/pbc/48-k_point_adc.py @@ -76,15 +76,15 @@ # Molecular IP-ADC myadc.method_type = 'ip' myadc.method = 'adc(3)' -mol_e_ip,mol_v_ip,mol_p_ip,mol_x_ip = myadc.kernel(nroots=nroots_test*np.product(nmp)) +mol_e_ip,mol_v_ip,mol_p_ip,mol_x_ip = myadc.kernel(nroots=nroots_test*np.prod(nmp)) # Molecular EA-ADC myadc.method_type = 'ea' myadc.method = 'adc(3)' -mol_e_ea,mol_v_ea,mol_p_ea,mol_x_ea = myadc.kernel(nroots=nroots_test*np.product(nmp)) +mol_e_ea,mol_v_ea,mol_p_ea,mol_x_ea = myadc.kernel(nroots=nroots_test*np.prod(nmp)) print("PBC KRHF Energy:", ekrhf) -print("PBC RHF Energy :", erhf / np.product(nmp)) +print("PBC RHF Energy :", erhf / np.prod(nmp)) print("PBC IP-ADC(3) roots:", k_e_ip) print("Mol IP-ADC(3) roots:", mol_e_ip) print("PBC EA-ADC(3) roots:", k_e_ea) diff --git a/pyscf/adc/__init__.py b/pyscf/adc/__init__.py index 230b763d9f..fdfa047839 100644 --- a/pyscf/adc/__init__.py +++ b/pyscf/adc/__init__.py @@ -70,8 +70,6 @@ def UADC(mf, frozen=None, mo_coeff=None, mo_occ=None): UADC.__doc__ = uadc.UADC.__doc__ def RADC(mf, frozen=None, mo_coeff=None, mo_occ=None): - __doc__ = radc.RADC.__doc__ - if not (frozen is None or frozen == 0): raise NotImplementedError diff --git a/pyscf/geomopt/__init__.py b/pyscf/geomopt/__init__.py index 80200c53c2..9b8c3125f4 100644 --- a/pyscf/geomopt/__init__.py +++ b/pyscf/geomopt/__init__.py @@ -22,6 +22,6 @@ def optimize(method, *args, **kwargs): except ImportError as e1: try: from . import berny_solver as geom - except ImportError as e2: + except ImportError: raise e1 return geom.optimize(method, *args, **kwargs) diff --git a/pyscf/gto/__init__.py b/pyscf/gto/__init__.py index 60be5e5d35..0dd37df489 100644 --- a/pyscf/gto/__init__.py +++ b/pyscf/gto/__init__.py @@ -24,5 +24,4 @@ from pyscf.gto.eval_gto import eval_gto from pyscf.gto import ecp -parse = basis.parse #import pyscf.gto.mole.cmd_args diff --git a/pyscf/gto/basis/__init__.py b/pyscf/gto/basis/__init__.py index dccda21443..87f0f68be7 100644 --- a/pyscf/gto/basis/__init__.py +++ b/pyscf/gto/basis/__init__.py @@ -531,7 +531,6 @@ def _convert_contraction(contr_string): def _truncate(basis, contr_scheme, symb, split_name): # keep only first n_keep contractions for each l contr_b = [] - b_index = 0 for l, n_keep in enumerate(contr_scheme): n_saved = 0 if n_keep > 0: diff --git a/pyscf/gto/mole.py b/pyscf/gto/mole.py index 2b72898cd4..429ca4fd7f 100644 --- a/pyscf/gto/mole.py +++ b/pyscf/gto/mole.py @@ -3840,7 +3840,7 @@ def _parse_nuc_mod(str_or_int_or_fn): def _update_from_cmdargs_(mol): try: # Detect whether in Ipython shell - __IPYTHON__ # noqa: + __IPYTHON__ # noqa: F821 return except Exception: pass diff --git a/pyscf/mp/ump2.py b/pyscf/mp/ump2.py index c6af263f50..547fdf95dd 100644 --- a/pyscf/mp/ump2.py +++ b/pyscf/mp/ump2.py @@ -303,7 +303,7 @@ def make_fno(mp, thresh=1e-6, pct_occ=None, nvir_act=None, t2=None, eris=None): masks = _mo_splitter(mp) if nvir_act is not None: - if isinstance(nvir_act, (int,numpy.int)): + if isinstance(nvir_act, (int, numpy.integer)): nvir_act = [nvir_act]*2 no_frozen = [] diff --git a/pyscf/pbc/adc/__init__.py b/pyscf/pbc/adc/__init__.py index 4e0d2c69c0..32db3be369 100644 --- a/pyscf/pbc/adc/__init__.py +++ b/pyscf/pbc/adc/__init__.py @@ -19,7 +19,6 @@ from pyscf.pbc.adc import kadc_rhf_ea def KRADC(mf, frozen=None, mo_coeff=None, mo_occ=None): - from pyscf.pbc.adc import kadc_rhf if not isinstance(mf, scf.khf.KRHF): mf = mf.to_rhf() return kadc_rhf.RADC(mf, frozen, mo_coeff, mo_occ) diff --git a/pyscf/pbc/dft/kgks.py b/pyscf/pbc/dft/kgks.py index d1187a8b33..e968a73af7 100644 --- a/pyscf/pbc/dft/kgks.py +++ b/pyscf/pbc/dft/kgks.py @@ -35,6 +35,7 @@ from pyscf.pbc.dft import rks from pyscf.pbc.dft import multigrid from pyscf.pbc.dft.numint2c import KNumInt2C +from pyscf.dft import gks as mol_ks from pyscf import __config__ def get_veff(ks, cell=None, dm=None, dm_last=0, vhf_last=0, hermi=1, @@ -117,6 +118,8 @@ class KGKS(rks.KohnShamDFT, kghf.KGHF): '''GKS class adapted for PBCs with k-point sampling. ''' + collinear = mol_ks.GKS.collinear + spin_samples = mol_ks.GKS.spin_samples get_veff = get_veff energy_elec = krks.energy_elec get_rho = krks.get_rho diff --git a/pyscf/pbc/dft/numint.py b/pyscf/pbc/dft/numint.py index 056349f88e..7332b45240 100644 --- a/pyscf/pbc/dft/numint.py +++ b/pyscf/pbc/dft/numint.py @@ -1264,7 +1264,7 @@ def _gen_rho_evaluator(self, cell, dms, hermi=0, with_lapl=False, grids=None): if isinstance(dms[0], numpy.ndarray) and dms[0].ndim == 2: mo_coeff = [mo_coeff] mo_occ = [mo_occ] - nao = cell.nao_nr() + nao = dms[0].shape[-1] ndms = len(mo_occ) def make_rho(idm, ao, non0tab, xctype): return self.eval_rho2(cell, ao, mo_coeff[idm], mo_occ[idm], diff --git a/pyscf/pbc/dft/numint2c.py b/pyscf/pbc/dft/numint2c.py index c645f66e57..5c8f8eedfc 100644 --- a/pyscf/pbc/dft/numint2c.py +++ b/pyscf/pbc/dft/numint2c.py @@ -335,6 +335,7 @@ def __init__(self, kpts=np.zeros((1,3))): collinear_samples = numint2c.NumInt2C.collinear_samples make_mask = lib.invalid_method('make_mask') + eval_ao = staticmethod(pnumint.eval_ao_kpts) def eval_rho(self, cell, ao_kpts, dm_kpts, non0tab=None, xctype='LDA', hermi=0, with_lapl=True, verbose=None): diff --git a/pyscf/pbc/dft/test/test_kgks.py b/pyscf/pbc/dft/test/test_kgks.py index 6ebd1165ea..89a3a8f681 100644 --- a/pyscf/pbc/dft/test/test_kgks.py +++ b/pyscf/pbc/dft/test/test_kgks.py @@ -133,11 +133,11 @@ def test_collinear_kgks_gga(self): cell.basis = [[0, [2, 1]], [1, [.5, 1]]] cell.spin = 1 cell.build() - mf = cell.KGKS(kpts=cell.make_kpts([3,1,1])) + mf = cell.KGKS(kpts=cell.make_kpts([2,1,1])) mf.xc = 'pbe' mf.collinear = 'col' mf.run() - self.assertAlmostEqual(mf.e_tot, -1.5775787238220018, 7) + self.assertAlmostEqual(mf.e_tot, -1.6373456924395708, 7) @unittest.skipIf(mcfun is None, "mcfun library not found.") def test_mcol_kgks_gga(self): @@ -150,12 +150,12 @@ def test_mcol_kgks_gga(self): cell.basis = [[0, [2, 1]], [1, [.5, 1]]] cell.spin = 1 cell.build() - mf = cell.KGKS(kpts=cell.make_kpts([3,1,1])).density_fit() + mf = cell.KGKS(kpts=cell.make_kpts([2,1,1])).density_fit() mf.xc = 'pbe' mf.collinear = 'mcol' mf._numint.spin_samples = 6 mf.run() - self.assertAlmostEqual(mf.e_tot, -1.5896974092061769, 7) + self.assertAlmostEqual(mf.e_tot, -1.6312141891350893, 6) def test_ncol_x2c_kgks_lda(self): from pyscf.pbc import gto @@ -167,11 +167,11 @@ def test_ncol_x2c_kgks_lda(self): cell.basis = [[0, [2, 1]], [1, [.5, 1]]] cell.spin = 1 cell.build() - mf = cell.KGKS(kpts=cell.make_kpts([3,1,1])).x2c() + mf = cell.KGKS(kpts=cell.make_kpts([2,1,1])).x2c() mf.xc = 'lda,vwn' mf.collinear = 'ncol' mf.run() - self.assertAlmostEqual(mf.e_tot, -1.4910121442258883, 7) + self.assertAlmostEqual(mf.e_tot, -1.533809531603591, 6) @unittest.skipIf(mcfun is None, "mcfun library not found.") def test_mcol_x2c_kgks_lda(self): @@ -184,11 +184,12 @@ def test_mcol_x2c_kgks_lda(self): cell.basis = [[0, [2, 1]], [1, [.5, 1]]] cell.spin = 1 cell.build() - mf = cell.KGKS(kpts=cell.make_kpts([3,1,1])).x2c() + mf = cell.KGKS(kpts=cell.make_kpts([2,1,1])).x2c() mf.xc = 'lda,vwn' + mf.collinear = 'mcol' mf._numint.spin_samples = 50 mf.run() - self.assertAlmostEqual(mf.e_tot, -1.4910121442258883, 6) + self.assertAlmostEqual(mf.e_tot, -1.533809531603591, 6) def test_to_hf(self): mf = dft.KGKS(cell).density_fit() diff --git a/pyscf/pbc/lib/kpts.py b/pyscf/pbc/lib/kpts.py index cd3de233e8..0e3c3e4d94 100644 --- a/pyscf/pbc/lib/kpts.py +++ b/pyscf/pbc/lib/kpts.py @@ -410,7 +410,7 @@ def symmetrize_wavefunction(kpts, psiR_k, mesh): ''' raise RuntimeError('need verification') psiR_k = np.asarray(psiR_k, order='C') - is_complex = psiR_k.dtype == np.complex + is_complex = np.iscomplexobj(psiR_k.dtype) nao = psiR_k.shape[1] nG = psiR_k.shape[2] psiR = np.zeros([kpts.nkpts,nao,nG], dtype = psiR_k.dtype, order='C') diff --git a/pyscf/pbc/scf/ghf.py b/pyscf/pbc/scf/ghf.py index e5f12dc901..fc5830006e 100644 --- a/pyscf/pbc/scf/ghf.py +++ b/pyscf/pbc/scf/ghf.py @@ -100,6 +100,7 @@ def get_jk(mf, cell=None, dm=None, hermi=0, kpt=None, kpts_band=None, class GHF(pbchf.SCF): '''GHF class for PBCs. ''' + _keys = {'with_soc'} def __init__(self, cell, kpt=np.zeros(3), exxdiv=getattr(__config__, 'pbc_scf_SCF_exxdiv', 'ewald')): diff --git a/pyscf/pbc/scf/kghf.py b/pyscf/pbc/scf/kghf.py index 2ce492067c..89259e4ea8 100644 --- a/pyscf/pbc/scf/kghf.py +++ b/pyscf/pbc/scf/kghf.py @@ -186,6 +186,8 @@ def fn_init_guess(mf, cell=None, kpts=None): class KGHF(khf.KSCF): '''GHF class for PBCs. ''' + _keys = {'with_soc'} + def __init__(self, cell, kpts=np.zeros((1,3)), exxdiv=getattr(__config__, 'pbc_scf_SCF_exxdiv', 'ewald')): khf.KSCF.__init__(self, cell, kpts, exxdiv) diff --git a/pyscf/solvent/__init__.py b/pyscf/solvent/__init__.py index 40cdc46106..eeba899549 100644 --- a/pyscf/solvent/__init__.py +++ b/pyscf/solvent/__init__.py @@ -12,7 +12,7 @@ # See the License for the specific language governing permissions and # limitations under the License. -from pyscf.solvent import ddcosmo, pcm +from pyscf.solvent import ddcosmo def ddCOSMO(method_or_mol, solvent_obj=None, dm=None): '''Initialize ddCOSMO model.