from __future__ import print_function, division
import numpy as np
import time, socket, getpass
from pynao.m_siesta_units import siesta_conv_coefficients
HARTREE2EV = siesta_conv_coefficients["ha2ev"]
[docs]def report_gw(self):
"""
Prints the energy levels of mean-field and G0W0
"""
import sys, re
if not hasattr(self, 'mo_energy_gw'):
if (self.__class__.__name__ == 'gw'):
self.kernel_gw()
elif (self.__class__.__name__ == 'gw_iter'):
self.kernel_gw_iter()
emfev = self.mo_energy[0].T * HARTREE2EV
egwev = self.mo_energy_gw[0].T * HARTREE2EV
file_name= ''.join(self.get_symbols())
with open('report_'+file_name+'.out', 'w') as out_file:
out_file.write('Job started on {}, @ {}-{}.\n'.format(time.strftime("%c",
time.localtime(self.start_time)), socket.gethostname(), getpass.getuser()))
if self.nspin == 1:
print('\n','='*50,'\n','='*12,'|G0W0 eigenvalues (eV)|','='*13,'\n','='*50)
out_file.write('='*15+'|G0W0 eigenvalues (eV)|'+'='*15+'\n')
line = "\n n %14s %14s %7s " % ("E_mf", "E_gw", "occ")
print(line)
out_file.write(line + "\n")
EFERMI, EHOMO, ELUMO = get_energy_levels_FHL(self, egwev)
for ie, (emf, egw, f) in enumerate(zip(emfev,egwev, self.mo_occ[0].T)):
line = "%5d %14.7g %14.7g %7.2g" % (ie, emf[0], egw[0], f[0])
print(line)
out_file.write(line + "\n")
line = '\nFermi energy (eV):%16.7g'%(EFERMI)
print(line)
out_file.write(line + "\n")
line = 'G0W0 HOMO energy (eV):%16.7g' % (EHOMO[0])
print(line)
out_file.write(line + "\n")
line = 'G0W0 LUMO energy (eV):%16.7g' % (ELUMO[0])
print(line)
out_file.write(line + "\n")
line = 'G0W0 HOMO-LUMO gap (eV):%16.7g' %(ELUMO[0] - EHOMO[0])
print(line)
out_file.write(line + "\n")
elif self.nspin == 2:
print('\n','='*85,'\n','='*30,'|G0W0 eigenvalues (eV)|','='*30,'\n','='*85)
out_file.write('='*31+'|G0W0 eigenvalues (eV)|'+'='*31+'\n')
line = "\n n %14s %14s %7s | %14s %14s %7s" % ("E_mf_up",
"E_gw_up",
"occ_up",
"E_mf_down",
"E_gw_down",
"occ_down")
print(line)
out_file.write(line + "\n")
EFERMI, EHOMO, ELUMO = get_energy_levels_FHL(self, egwev)
for ie, (emf, egw, f) in enumerate(zip(emfev, egwev, self.mo_occ[0].T)):
line = "%5d %14.7g %14.7g %7.2g | %14.7g %14.7g %7.2g" % (ie,
emf[0],
egw[0],
f[0],
emf[1],
egw[1],
f[1])
print(line)
out_file.write (line + "\n")
line = '\nFermi energy (eV):%16.7g'%(EFERMI)
print(line)
out_file.write(line + "\n")
line = 'G0W0 HOMO energy (eV):%16.7f %16.7g'%(EHOMO[0], EHOMO[1])
print(line)
out_file.write(line + "\n")
line = 'G0W0 LUMO energy (eV):%16.7f %16.7g' % (ELUMO[0], ELUMO[1])
print(line)
out_file.write(line + "\n")
line = 'G0W0 HOMO-LUMO gap (eV):%16.7f %16.7g' % (ELUMO[0] - EHOMO[0],
ELUMO[1] - EHOMO[1])
print(line)
out_file.write(line + "\n")
else:
raise RuntimeError('not implemented...')
line = 'G0W0 Total energy (eV):%16.7g' % (self.etot_gw*HARTREE2EV)
print(line)
out_file.write(line + "\n")
for s in range(self.nspin):
swap = self.argsort[s][self.start_st[s]:self.finish_st[s]]
if (np.allclose(swap,np.sort(swap))==False):
line = "Warning: Spin {}, swapping in QP orbital energies has happened!".format(s+1)
print(line)
out_file.write(line + "\n")
line = 'Energy-sorted MO indices: \t {}'.format(swap)
print(line)
out_file.write(line + "\n")
out_file.write('\n=====| Timings of main algorithms |=====')
timing, steps = report_gw_t(self,ret=True)
for i in zip(timing, steps):
if (round(i[0],3) != 0):
out_file.write('\n{:20.19s}:{:14.2f} secs'.format(i[1],i[0]))
elapsed_time = self.time_gw[-1]-self.time_gw[0]
out_file.write('\n'+'-'*20+'\nTotal spent time :{:14.2f} secs\n'.format (elapsed_time))
finish_t = elapsed_time + self.start_time
line = '\nJOB DONE! \t {}'.format(time.strftime("%c", time.localtime(finish_t)))
print(line)
out_file.write(line + "\n")
[docs]def report_mfx(self, dm1=None):
"""
This collects the h_core, Hartree (K) and Exchange(J) and Fock expectation value with given density matrix.
"""
#import re
#import sys
#file_name= ''.join(self.get_symbols())
#sys.stdout=open('report_mf_'+file_name+'.out','w')
if hasattr(self, 'mf'): report_mocc(self)
if dm1 is None: dm1 = self.make_rdm1()
if dm1.ndim==5 : dm1=dm1[0,...,0]
assert dm1.shape == (self.nspin, self.norbs, self.norbs), "Given density matrix has wrong dimension"
H = self.get_hcore()
ecore = (H*dm1).sum()
co = self.get_j()
fock = self.get_fock()
exp_h = np.zeros((self.nspin, self.norbs))
exp_co = np.zeros((self.nspin, self.norbs))
exp_x = np.zeros((self.nspin, self.norbs))
exp_f = np.zeros((self.nspin, self.norbs))
if not hasattr(self, 'kmat'):
self.kmat = self.get_k()
if self.nspin==1:
x = -0.5* self.kmat
mat_h = np.dot(self.mo_coeff[0,0,:,:,0], H)
exp_h = np.einsum('nb,nb->n', mat_h, self.mo_coeff[0,0,:,:,0])
mat_co = np.dot(self.mo_coeff[0,0,:,:,0], co)
exp_co = np.einsum('nb,nb->n', mat_co, self.mo_coeff[0,0,:,:,0])
mat_x = np.dot(self.mo_coeff[0,0,:,:,0], x)
exp_x = np.einsum('nb,nb->n', mat_x, self.mo_coeff[0,0,:,:,0])
mat_f = np.dot(self.mo_coeff[0,0,:,:,0], fock)
exp_f = np.einsum('nb,nb->n', mat_f, self.mo_coeff[0,0,:,:,0])
print('\n','='*20,'| Restricted HF expectation values (eV) |','='*20)
print('%2s %13s %13s %13s %13s %13s %3s'%('no.','<H_core>','<K> ','<Sigma_x>','Fock ','MF energy ','occ'))
for i, (a,b,c,d,e,f) in enumerate(zip(exp_h.T*HARTREE2EV,exp_co.T*HARTREE2EV, exp_x.T*HARTREE2EV,exp_f.T*HARTREE2EV,self.mo_energy.T*HARTREE2EV, self.mo_occ[0].T)):
if (i==self.nfermi[0]): print('-'*83)
print(' %3d %13.6f %13.6f %13.6f %13.6f %13.6f %3d'%(i, a,b,c,d,e,f))
Vha = 0.5*(co*dm1).sum()
EX = 0.5*(x*dm1).sum()
elif self.nspin==2:
x = -self.kmat
cou = co[0]+co[1]
Vha = 0.5*(cou*dm1).sum()
for s in range(self.nspin):
mat_h = np.dot(self.mo_coeff[0,s,:,:,0], H)
exp_h[s] = np.einsum('nb,nb->n', mat_h, self.mo_coeff[0,s,:,:,0])
mat_co = np.dot(self.mo_coeff[0,s,:,:,0], cou)
exp_co[s] = np.einsum('nb,nb->n', mat_co, self.mo_coeff[0,s,:,:,0])
mat_x = np.dot(self.mo_coeff[0,s,:,:,0], x[s])
exp_x[s] = np.einsum('nb,nb->n', mat_x, self.mo_coeff[0,s,:,:,0])
mat_f = np.dot(self.mo_coeff[0,s,:,:,0], fock[s])
exp_f[s] = np.einsum('nb,nb->n', mat_f, self.mo_coeff[0,s,:,:,0])
print('\n','='*59,'| Unrestricted HF expectation values (eV) |','='*60)
print('%2s %13s %13s %13s %13s %13s %3s |%13s %13s %13s %13s %13s %3s '%('no.','<H_core>','<K> ','<Sigma_x>','Fock ','MF energy ','occ','<H_core>','<K> ','<Sigma_x>','Fock ','MF energy','occ'))
for i , (a,b,c,d,e,f) in enumerate(zip(exp_h.T*HARTREE2EV,exp_co.T*HARTREE2EV, exp_x.T*HARTREE2EV,exp_f.T*HARTREE2EV,self.mo_energy.T*HARTREE2EV, self.mo_occ[0].T)):
if (i==self.nfermi[0] or i==self.nfermi[1]): print('-'*163)
print(' %3d %13.6f %13.6f %13.6f %13.6f %13.6f %3d | %13.6f %13.6f %13.6f %13.6f %13.6f %3d'%(i, a[0],b[0],c[0],d[0],e[0],f[0],a[1],b[1],c[1],d[1],e[1],f[1]))
EX = 0.5*(x*dm1).sum()
if hasattr(self, 'mf'):
print('\nmean-field Nucleus-Nucleus (Ha):%16.6f'%(self.energy_nuc()))
print('mean-field Core energy (Ha):%16.6f'%(ecore))
print('mean-field Exchange energy (Ha):%16.6f'%(EX))
print('mean-field Hartree energy (Ha):%16.6f'%(Vha))
print('mean-field Total energy (Ha):%16.6f'%(self.mf.e_tot))
if (self.nspin==2):
sp = self.spin/2
s_ref = sp*(sp+1)
ss = self.mf.spin_square()
print('<S^2> and 2S+1 :%16.7f %16.7f'%(ss[0],ss[1]))
print('Instead of :%16.7f %16.7f'%(s_ref, 2*sp+1))
elapsed_time = time.time() - self.start_time
print('\nMean-field running time is: {}'.format(time.strftime("%H:%M:%S", time.gmtime(elapsed_time))))
#sys.stdout.close()
[docs]def report_mocc (self):
'''Lists indeces of occupied and virtual orbitals in a limited range'''
m, l = min(self.start_st),max(self.finish_st)
if (self.nspin==1):
print('==|Occupations|==\nno. occ')
for i in range(m,l,1): print('{:3d} {:3.1f}'.format(i,self.mo_occ[0,0,i]))
elif (self.nspin==2):
print('===|Occupations|===\nno. occ occ')
for i in range(m,l,1): print('{:3d} {:3.1f} {:3.1f}'.format(i,self.mo_occ[0,0,i],self.mo_occ[0,1,i]))
[docs]def sigma_xc(self):
"""
This calculates the Exchange expectation value and correlation part of self energy, when:
self.get_k() = Exchange operator/energy
mat1 is product of this operator and molecular coefficients and it will be diagonalized in expval by einsum
Sigma_c = E_GW - E_HF
"""
if not hasattr(self, 'kmat'):
self.kmat = self.get_k()
if self.nspin==1:
mat = -0.5*self.kmat
mat1 = np.dot(self.mo_coeff[0,0,:,:,0], mat)
expval = np.einsum('nb,nb->n', mat1, self.mo_coeff[0,0,:,:,0]).reshape((1,self.norbs))
print('===| Expectationvalues of Exchange energy(eV) |===\n %3s %16s %3s'%('no.','<Sigma_x> ','occ'))
for i, (a,b) in enumerate(zip(expval.T*HARTREE2EV,self.mo_occ[0].T)): #self.h0_vh_x_expval[0,:self.nfermi[0]+5] to limit the virual states
if (i==self.nfermi[0]): print('-'*50)
print (' %3d %16.6f %3d'%(i,a[0], b[0]))
elif self.nspin==2:
mat = -self.kmat
expval = np.zeros((self.nspin, self.norbs))
for s in range(self.nspin):
mat1 = np.dot(self.mo_coeff[0,s,:,:,0], mat[s])
expval[s] = np.einsum('nb,nb->n', mat1, self.mo_coeff[0,s,:,:,0])
print('========| Expectationvalues of Exchange energy(eV) |========\n %3s %16s %3s | %13s %4s'%('no.','<Sigma_x>','occ','<Sigma_x>','occ'))
for i , (a,b) in enumerate(zip(expval.T* HARTREE2EV,self.mo_occ[0].T)):
if (i==self.nfermi[0] or i==self.nfermi[1]): print('-'*60)
print(' %3d %16.6f %3d | %13.6f %3d'%(i, a[0],b[0],a[1], b[1]))
if hasattr(self,'mo_energy_gw'):
sigma_gw_c = self.mo_energy_gw - self.mo_energy
#sigma_gw_c= np.asanyarray([gw.mo_energy_gw[0,s,nn] - gw.mo_energy[0,s,nn] for s,nn in enumerate(gw.nn) ]) #Only corrected by GW not scisorres
if self.nspin==1:
print('\n===| Correlation contribution of GW@HF (eV) |===\n %3s %16s %3s'%('no.','<Sigma_c> ','occ'))
for i, (a,b) in enumerate(zip(sigma_gw_c.T*HARTREE2EV,self.mo_occ[0].T)): #self.h0_vh_x_expval[0,:self.nfermi[0]+5] to limit the virual states
if (i==self.nfermi[0]): print('-'*48)
print (' %3d %16.6f %3d'%(i,a[0], b[0]))
elif self.nspin==2:
print('\n========| Correlation contribution of GW@HF (eV) |=========\n %3s %16s %3s | %13s %4s'%('no.','<Sigma_c>','occ','<Sigma_c>','occ'))
for i , (a,b) in enumerate(zip(sigma_gw_c.T* HARTREE2EV,self.mo_occ[0].T)):
if (i==self.nfermi[0] or i==self.nfermi[1]): print('-'*60)
print(' %3d %16.6f %3d | %13.6f %3d'%(i, a[0],b[0],a[1], b[1]))
[docs]def report_gw_t(self, ret=None):
"""Lists spent time within main GW calculation"""
steps = ['initialize NAO', #01
'get_h0_vh_x_expval',' 1-get_K', ' 2-get_J',#23-45-67
'G0W0_eigvals',' 1-snmw2sf',' 1-1-XVX',' 1-2-chi_0', #89-1011-1213-1415
' 2-corr_int_tot', #1617
' 3-corr_res_tot',' 3-1-chi_0', #1819-2021
' 4-veffmatvec', #2223
'scissor', #2425
'etot_gw'] #2627
timing = list(self.time_gw[1::2]- self.time_gw[0:-1:2])
print('\n=====| Timings of main algorithms |=====')
for i in zip(timing, steps):
if (round(i[0],3) != 0): print('{:20.19s}:{:14.2f} secs'.format(i[1],i[0]))
print('-'*20,'\nTotal spent time :{:14.2f} secs'.format (self.time_gw[-1]-self.time_gw[0]),'\n')
if ret:
return timing, steps
[docs]def get_energy_levels_FHL(self, egwev):
"""
Return the values of the Fermi, HOMO and LUMO energies in eV
"""
EFermi = self.fermi_energy*HARTREE2EV
EHOMO = np.zeros((self.nspin), dtype=self.dtype)
ELUMO = np.zeros((self.nspin), dtype=self.dtype)
for spin in range(self.nspin):
for ilevel, El in enumerate(egwev[:, spin]):
if El > EFermi:
break
EHOMO[spin] = egwev[ilevel-1, spin]
ELUMO[spin] = egwev[ilevel, spin]
return EFermi, EHOMO, ELUMO
#
# Example of reporting expectation values of mean-field calculations.
#
if __name__=='__main__':
import numpy as np
from pyscf import gto, scf
from pynao import gw as gw_c
mol = gto.M(verbose=0, atom='O 0.0, 0.0, 0.622978 ; O 0.0, 0.0, -0.622978',
basis='cc-pvdz', spin=2, charge=0)
gto_mf = scf.UHF(mol)
gto_mf.kernel()
gw = gw_c(mf=gto_mf, gto=mol, verbosity=3, niter_max_ev=1, kmat_algo='sm0_sum')
#gw.report_mf() # prints the energy levels of mean-field components
gw.report() # gives G0W0 spectra