diff --git a/src/bsemol/GNUmakefile b/src/bsemol/GNUmakefile index 3e690a8eaa..c1f2a3f16b 100644 --- a/src/bsemol/GNUmakefile +++ b/src/bsemol/GNUmakefile @@ -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 diff --git a/src/bsemol/bse_davidson.F b/src/bsemol/bse_davidson.F index a944388cea..33e3228452 100644 --- a/src/bsemol/bse_davidson.F +++ b/src/bsemol/bse_davidson.F @@ -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) diff --git a/src/bsemol/bse_davidson_guess.F b/src/bsemol/bse_davidson_guess.F index 8cbd7cef66..735640ef55 100644 --- a/src/bsemol/bse_davidson_guess.F +++ b/src/bsemol/bse_davidson_guess.F @@ -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) ') diff --git a/src/bsemol/bse_input.F b/src/bsemol/bse_input.F index dfad98ab3f..188be807ce 100644 --- a/src/bsemol/bse_input.F +++ b/src/bsemol/bse_input.F @@ -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) diff --git a/src/bsemol/bse_oscstr.F b/src/bsemol/bse_oscstr.F index a227d7f203..161f536db1 100644 --- a/src/bsemol/bse_oscstr.F +++ b/src/bsemol/bse_oscstr.F @@ -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,' = ',f8.4) + 9180 format(2x,' = ',f8.4) + 9190 format(2x,' D = ',f8.4) end diff --git a/src/bsemol/bse_s2.F b/src/bsemol/bse_s2.F new file mode 100644 index 0000000000..03e1a1c7a6 --- /dev/null +++ b/src/bsemol/bse_s2.F @@ -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 +