comp-chem/NWChem
1
0
Fork 0
mirror of https://github.com/nwchemgit/nwchem.git synced 2026-07-28 06:05:44 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 1

Add energy window restriction to trials

Browse source
This commit is contained in:
Daniel Mejia-Rodriguez 2023-08-04 14:22:13 -07:00 • committed by edoapra
parent 6f9a0a8248
commit a39f8e4e36
No known key found for this signature in database
GPG key ID: 9E6A0B70826967BA
6 changed files with 246 additions and 47 deletions
Download patch file Download diff file Expand all files Collapse all files

3
src/bsemol/GNUmakefile
View file

@ -6,7 +6,8 @@ OBJ = bsemol.o bse_input.o bse_init.o bse_analytic.o \
bse_davidson_diagonal.o \
bse_davidson_kvec.o bse_davidson_guess.o \
bse_davidson_kortho.o bse_davidson_restart.o \
bse_davidson_newvecs.o bse_davidson_residual.o
bse_davidson_newvecs.o bse_davidson_residual.o \
bse_s2.o
USES_BLAS = bse_buildw.F bse_ri_init.F bse_wmn.F

2
src/bsemol/bse_davidson.F
View file

@ -515,10 +515,8 @@
enddo
endif
call ga_print(pars%g_apb)
call ga_add(0.5d0,pars%g_apb,0.5d0,pars%g_amb,pars%g_apb)
call ga_add(1d0,pars%g_apb,-1d0,pars%g_amb,pars%g_amb)
call bse_oscstr(pars,dbl_mb(k_omega),npoles)
stat = ma_chop_stack(l_dia)

48
src/bsemol/bse_davidson_guess.F
View file

@ -16,8 +16,10 @@
character*4 oname,vname
logical found
integer amo,imo,itrial,ipole,jpole,indeces(ntrials)
integer isym,jsym,i,j,pol
integer isym,jsym,i,j,pol,idum
integer allindeces(npoles)
double precision minimum
logical swapped
call nga_get(pars%g_apb,1,npoles,diagonal,npoles)
call nga_get(pars%g_amb,1,npoles,amb,npoles)
@ -26,13 +28,52 @@
diagonal(:) = diagonal(:)*amb(:)
endif
indeces(:) = -1
allindeces(:) = -1
do ipole=1,npoles
allindeces(ipole) = ipole
enddo
do ipole=npoles-1,1,-1
swapped = .false.
do jpole=1,ipole
if (diagonal(allindeces(jpole)).gt.
$ diagonal(allindeces(jpole+1))) then
idum = allindeces(jpole)
allindeces(jpole) = allindeces(jpole+1)
allindeces(jpole+1) = idum
swapped = .true.
endif
enddo
if(.not.swapped) exit
enddo
if (pars%dowindow) then
do ipole=1,npoles
if (dsqrt(diagonal(allindeces(ipole))).ge.pars%elower) then
if (ipole+ntrials-1.gt.npoles) then
indeces(1:ntrials) = allindeces(npoles-ntrials+1:npoles)
else
indeces(1:ntrials) = allindeces(ipole:ipole+ntrials-1)
endif
exit
endif
enddo
else
indeces(1:ntrials) = allindeces(1:ntrials)
endif
! Get the indeces of the smallest diagonal elements
! within energy window
#if 0
do itrial=1,ntrials
minimum = 1d6
do ipole=1,npoles
if (abs(diagonal(ipole)).lt.minimum) then
found = .true.
if (pars%dowindow) then
else
found = .true.
endif
do jpole=1,itrial
if (ipole.eq.indeces(jpole)) then
found = .false.
@ -46,6 +87,7 @@
endif
enddo
enddo
#endif
! Get symmetry of ground state
isym = 1
@ -86,7 +128,7 @@
endif
call sym_irrepname(geom, jsym, vname)
write(luout,9030) itrial,pol,imo,amo,vname,
$ dsqrt(diagonal(ipole))*27.2114d0
$ dsqrt(diagonal(ipole))*27.211396d0
enddo
9000 format(i5,' smallest diagonals (eV) ')

2
src/bsemol/bse_input.F
View file

@ -100,7 +100,7 @@ c
& call errquit(pname//'rtdb_put failed', 300, RTDB_ERR)
goto 10
c
c compute singlet transitions
c compute triplet transitions
c
500 if (.not.rtdb_put(rtdb, 'bse:singlet', mt_log, 1, .false.))
& call errquit(pname//'rtdb_put failed', 300, RTDB_ERR)

139
src/bsemol/bse_oscstr.F
View file

@ -21,21 +21,26 @@
integer llpole,ulpole,iroot,ixyz,ipole,mag_transdip(3)
integer isym,jsym,i,j,k,l,pol,occ,vir
integer isym,jsym,i,j,k,l,pol,occ,vir,ioff
integer g_td(2),g_xy(2),g_dipmag, g_ovlp, g_work, g_corr
integer k_corr, l_corr
integer k_corr, l_corr, ulloop
character*4 oname,vname
logical ltransden,stat
double precision origin(3)
double precision na,nb,s2
double precision na,nb,s2,s2_full
data origin/0d0,0d0,0d0/
character*(nw_max_path_len) fn_civecs, fn_transden
logical,external :: dmat_to_file
integer,external :: ga_create_atom_blocked
character*5 spin(2)
double precision x(totpoles),y(totpoles)
double precision thresh ! Threshold for printing
parameter (thresh=5.0d-2)
spin(1)='alpha'
spin(2)='beta '
if(.not.rtdb_get(pars%rtdb,'bse:ltransden',mt_log,1,ltransden))
$ ltransden = .false.
@ -83,7 +88,7 @@
! Get the lower and upper limits
llpole = -1
if (pars%dowindow) then
if (pars%dowindow .and. (.not.pars%davidson)) then
iroot = 0
do ipole=1,npoles
if (omega(ipole).lt.pars%elower) cycle
@ -184,28 +189,29 @@
enddo
endif
if (pars%davidson.and.pars%ipol.eq.2) then
call ga_get(pars%g_apb,1,npoles,ipole,ipole,x,npoles)
call ga_get(pars%g_amb,1,npoles,ipole,ipole,y,npoles)
call tddft_s2(pars%tda, (/0,0/), pars%nocc,
& (/pars%nmo,pars%nmo/), (/0,0/),
& pars%npoles,x,x(pars%npoles(1)+1),y,
& y(pars%npoles(1)+1),dbl_mb(k_corr),s2)
call bse_s2(pars%nocc,pars%nmo,pars%npoles,x,
& x(pars%npoles(1)+1), y, y(pars%npoles(1)+1), dbl_mb(k_corr),
& s2_full)
if (abs(s2_full).lt.0.2d0) then
key = 'singlet'
elseif (abs(abs(s2_full)-2.0d0).lt.0.2) then
key = 'triplet'
else
key = ' '
endif
endif
if (pars%me.eq.0) then
call tddft_oscstr(moms, omega(ipole), mag_transdip,
& oscstr, oscstr_tot)
c do ixyz=1,3
c moms(ixyz) = moms(ixyz) + nga_ddot_patch(pars%g_apb,'n',
c & (/1,ipole/),(/pars%npoles(1),ipole/),
c & pars%g_dipa(ixyz),'n',1,pars%npoles(1))
c oscstr = oscstr + moms(ixyz)**2
c enddo
c if (pars%ipol.gt.1) then
c do ixyz=1,3
c moms(ixyz) = moms(ixyz) +
c & nga_ddot_patch(pars%g_apb,'n',(/pars%npoles(1)+1,ipole/),
c & (/totpoles,ipole/),
c & pars%g_dipb(ixyz),'n',1,pars%npoles(2))
c oscstr = oscstr + moms(ixyz)**2
c enddo
c endif
c oscstr = 2d0/3d0*omega(ipole)*oscstr
l = 0
maxr = 0.0d0
do i=1,pars%ipol
@ -227,28 +233,13 @@ c oscstr = 2d0/3d0*omega(ipole)*oscstr
jsym = ieor(jsym-1,int_mb(pars%k_irs(pol)+vir-1)-1)+1
call sym_irrepname(geom, jsym, vname)
if (pars%ipol.eq.2) then
na = dble(pars%nocc(1))
nb = dble(pars%nocc(2))
s2 = na - (na-nb)/2d0 + ((na-nb)/2d0)**2
do j=1,pars%nocc(1)
do k=1,pars%nocc(2)
enddo
enddo
endif
write(luout,9110) iroot,key,vname,omega(ipole),
& omega(ipole)*ha2ev
if (pars%davidson.and.pars%ipol.eq.2) then
call ga_get(pars%g_apb,1,npoles,ipole,ipole,x,npoles)
call ga_get(pars%g_amb,1,npoles,ipole,ipole,y,npoles)
call tddft_s2(pars%tda, (/0,0/), pars%nocc,
& (/pars%nmo,pars%nmo/), (/0,0/),
& pars%npoles,x,x(pars%npoles(1)+1),y,
& y(pars%npoles(1)+1),dbl_mb(k_corr),s2)
write(luout,*)
write(luout,9180) s2
write(luout,9190) s2_full
endif
if (pars%singlet) then
@ -261,6 +252,73 @@ c oscstr = 2d0/3d0*omega(ipole)*oscstr
write(luout,9210)
endif
if (pars%davidson) then
write(luout,*)
if (ipol.eq.1) then
l=0
do j=1,pars%nocc(1)
do k=pars%nocc(1)+1,pars%nmo
l=l+1
if (pars%tda) then
r=x(l)
call sym_irrepname
1 (geom,int_mb(pars%k_irs(1)+j-1),oname)
call sym_irrepname
1 (geom,int_mb(pars%k_irs(1)+k-1),vname)
if ((dabs(r(1)).gt.thresh))
2 write(LuOut,9120) j,oname,k,vname,r
else
r=x(l)
call sym_irrepname
1 (geom,int_mb(pars%k_irs(1)+j-1),oname)
call sym_irrepname
1 (geom,int_mb(pars%k_irs(1)+k-1),vname)
if ((dabs(r(1)).gt.thresh))
2 write(LuOut,9130) j,oname,k,vname,r
r=y(l)
if ((dabs(r(1)).gt.thresh))
2 write(LuOut,9140) j,oname,k,vname,r
endif
enddo
enddo
else
do i=1,ipol
ioff = (i-1)*pars%npoles(1)
l=0
do j=1,pars%nocc(i)
do k=pars%nocc(i)+1,pars%nmo
l=l+1
if (pars%tda) then
r = x(ioff+l)
call sym_irrepname
1 (geom,int_mb(pars%k_irs(i)+j-1),oname)
call sym_irrepname
1 (geom,int_mb(pars%k_irs(i)+k-1),vname)
if ((dabs(r(1)).gt.thresh))
2 write(LuOut,9150) j,spin(i),oname,
3 k,spin(i),vname,r
else
r = x(ioff+l)
call sym_irrepname
1 (geom,int_mb(pars%k_irs(i)+j-1),oname)
call sym_irrepname
1 (geom,int_mb(pars%k_irs(i)+k-1),vname)
if ((dabs(r(1)).gt.thresh))
2 write(LuOut,9160) j,spin(i),oname,
3 k,spin(i),vname,r
r=y(ioff+l)
if ((dabs(r(1)).gt.thresh))
2 write(LuOut,9170) j,spin(i),oname,
3 k,spin(i),vname,r
endif
enddo
enddo
enddo
endif
call util_flush(LuOut)
endif
endif
@ -301,7 +359,8 @@ c oscstr = 2d0/3d0*omega(ipole)*oscstr
1 i5,1x,a5,1x,a4,f8.5,' X')
9170 format(5x,'Occ.',i5,1x,a5,1x,a4,'--- Virt.',
1 i5,1x,a5,1x,a4,f8.5,' Y')
9180 format(2x,'<S2> = ',f8.4)
9180 format(2x,' <S2> = ',f8.4)
9190 format(2x,' D<S2> = ',f8.4)
end

99
src/bsemol/bse_s2.F Normal file
View file

@ -0,0 +1,99 @@
subroutine bse_s2(nocc, nmo, npoles, xup, xdw, yup, ydw,
$ ovlp, s2)
implicit none
integer nocc(2), npoles(2), nmo
double precision xup(npoles(1)), xdw(npoles(2))
double precision yup(npoles(1)), ydw(npoles(2))
double precision s2
double precision ovlp(nmo,nmo)
double precision ov,xx,yy
integer i,j,a,b,ipole,jpole,nvir(2)
integer llipole,ulipole,lljpole,uljpole
double precision, external :: ddot
s2 = 0d0
nvir(1) = nmo - nocc(1)
nvir(2) = nmo - nocc(2)
! Up-Up
do i=1,nocc(1)
do j=1,nocc(1)
llipole = (i-1)*nvir(1) + 1
ulipole = i*nvir(1)
lljpole = (j-1)*nvir(1) + 1
uljpole = j*nvir(1)
ov = dot_product(ovlp(i,1:nocc(2)),ovlp(j,1:nocc(2)))
xx = dot_product(xup(llipole:ulipole),xup(lljpole:uljpole))
yy = dot_product(yup(llipole:ulipole),yup(lljpole:uljpole))
s2 = s2 + ov*(xx+yy)
enddo
enddo
do a=nocc(1)+1,nmo
do b=nocc(1)+1,nmo
ov = dot_product(ovlp(a,1:nocc(2)),ovlp(b,1:nocc(2)))
llipole = a - nocc(1)
lljpole = b - nocc(1)
s2 = s2 - ov*ddot(nocc(1),xup(llipole),nvir(1),
$ xup(lljpole),nvir(1))
s2 = s2 - ov*ddot(nocc(1),yup(llipole),nvir(1),
$ yup(lljpole),nvir(1))
enddo
enddo
! Dw-Dw
do i=1,nocc(2)
do j=1,nocc(2)
llipole = (i-1)*nvir(2) + 1
ulipole = i*nvir(2)
lljpole = (j-1)*nvir(2) + 1
uljpole = j*nvir(2)
ov = dot_product(ovlp(1:nocc(1),i),ovlp(1:nocc(1),j))
xx = dot_product(xdw(llipole:ulipole),xdw(lljpole:uljpole))
yy = dot_product(ydw(llipole:ulipole),ydw(lljpole:uljpole))
s2 = s2 + ov*(xx+yy)
enddo
enddo
do a=nocc(2)+1,nmo
do b=nocc(2)+1,nmo
ov = dot_product(ovlp(1:nocc(1),a),ovlp(1:nocc(1),b))
llipole = a - nocc(2)
lljpole = b - nocc(2)
s2 = s2 - ov*ddot(nocc(2),xdw(llipole),nvir(2),
$ xdw(lljpole),nvir(2))
s2 = s2 - ov*ddot(nocc(2),ydw(llipole),nvir(2),
$ ydw(lljpole),nvir(2))
enddo
enddo
! Up-Dw
do ipole=1,npoles(1)
i = ipole/nvir(1) + 1
a = ipole - (i-1)*nvir(1) + nocc(1)
do jpole=1,npoles(2)
j = jpole/nvir(2) + 1
b = jpole - (j-1)*nvir(2) + nocc(2)
s2 = s2 -
$ 2d0*xup(ipole)*xdw(jpole)*ovlp(i,j)*ovlp(a,b) -
$ 2d0*yup(ipole)*ydw(jpole)*ovlp(i,j)*ovlp(a,b) +
$ 2d0*xup(ipole)*ydw(jpole)*ovlp(a,j)*ovlp(i,b)
enddo
enddo
! Dw-Up
do ipole=1,npoles(2)
i = ipole/nvir(2) + 1
a = ipole - (i-1)*nvir(2) + nocc(2)
do jpole=1,npoles(1)
j = jpole/nvir(1) + 1
b = jpole - (j-1)*nvir(1) + nocc(1)
s2 = s2 +
$ 2d0*xdw(ipole)*yup(jpole)*ovlp(j,a)*ovlp(b,i)
enddo
enddo
end subroutine