diff --git a/src/bsemol/GNUmakefile b/src/bsemol/GNUmakefile index 63bbd26bf6..bf8113b5ea 100644 --- a/src/bsemol/GNUmakefile +++ b/src/bsemol/GNUmakefile @@ -1,6 +1,8 @@ OBJ = bsemol.o bse_input.o bse_init.o bse_analytic.o \ bse_davidson.o bse_lanczos.o bse_defaults.o \ - bse_finalize.o bse_ri_init.o bse_wmn.o bse_casida.o + bse_finalize.o bse_ri_init.o bse_wmn.o \ + bse_analytic_cis.o bse_analytic_tda.o \ + bse_dipole.o bse_oscstr.o USES_BLAS = bse_buildw.F bse_ri_init.F bse_wmn.F diff --git a/src/bsemol/bse.fh b/src/bsemol/bse.fh index 9f177f729f..4c62604599 100644 --- a/src/bsemol/bse.fh +++ b/src/bsemol/bse.fh @@ -10,13 +10,12 @@ c integer :: natoms integer :: nmo, nri - integer :: nocc(2),nvir(2),npoles(2) - integer :: noqp(2),nvqp(2),nstates(2),mynpoles(2) + integer :: nocc(2),nvir(2),npoles(2),mynpoles(2) integer :: oolo(2),oohi(2) integer :: ovlo(2),ovhi(2) integer :: vvlo(2),vvhi(2) integer :: llmo(2),nqp(2) - integer :: ngl + integer :: nroots double precision :: tol2e, eta c integer :: l_mf_evals, k_mf_evals @@ -27,11 +26,13 @@ c c integer :: g_movecs(2), g_moeris, g_sols integer :: g_erioo(2), g_eriov(2), g_erivv(2), g_erim - integer :: g_wov(2), g_apb + integer :: g_wov(2), g_apb, g_dipa(3), g_dipb(3) logical :: noio c logical :: ncap, davidson, lanczos, analytic, tda + logical :: singlet, triplet, dowindow double precision :: thresh, fdstep, vxddp(2), vxddm(2) + double precision :: elower, eupper end type bse_params_t diff --git a/src/bsemol/bse_analytic.F b/src/bsemol/bse_analytic.F index 85294a6e5c..273a19e22c 100644 --- a/src/bsemol/bse_analytic.F +++ b/src/bsemol/bse_analytic.F @@ -59,12 +59,16 @@ enddo ! Build and diagonalize Casida-like matrix - call bse_casida(pars,dbl_mb(k_dia),dbl_mb(k_omega),maxpoles, - & totnpoles) + if (pars%tda) then + call bse_analytic_tda(pars,dbl_mb(k_dia),dbl_mb(k_omega), + & maxpoles,totnpoles,pars%singlet) + else + call bse_analytic_cis(pars,dbl_mb(k_dia),dbl_mb(k_omega), + & maxpoles,totnpoles,pars%singlet) + endif ! Compute oscillator strengths - - ! Print output + call bse_oscstr(pars,dbl_mb(k_omega),totnpoles) ! Destroy GAs if (.not.ga_destroy(pars%g_apb)) diff --git a/src/bsemol/bse_casida.F b/src/bsemol/bse_analytic_cis.F similarity index 95% rename from src/bsemol/bse_casida.F rename to src/bsemol/bse_analytic_cis.F index 53849aa2ed..7d73fa5452 100644 --- a/src/bsemol/bse_casida.F +++ b/src/bsemol/bse_analytic_cis.F @@ -1,4 +1,5 @@ - subroutine bse_casida(pars,wia,omega,maxpoles,totpoles) + subroutine bse_analytic_cis(pars,wia,omega,maxpoles,totpoles, + $ singlet) implicit none #include "mafdecls.fh" #include "errquit.fh" @@ -12,6 +13,7 @@ type(bse_params_t) :: pars + logical singlet integer maxpoles, totpoles double precision wia(maxpoles,pars%ipol) double precision omega(totpoles) @@ -93,9 +95,11 @@ if (.not.ga_destroy(g_tmp)) & call errquit(pname//'could not destroy temp GA',1,GA_ERR) + if (singlet) then ! Hartree contribution - call ga_dgemm('t','n',totpoles,totpoles,pars%nri,4d0, - & pars%g_eriov(1),pars%g_eriov(1),1d0,pars%g_apb) + call ga_dgemm('t','n',totpoles,totpoles,pars%nri,4d0, + & pars%g_eriov(1),pars%g_eriov(1),1d0,pars%g_apb) + endif ! Diagonal term call ga_add_diagonal(pars%g_apb,g_diag) @@ -157,8 +161,6 @@ enddo !$omp end parallel do - write(*,*) omega(1:20) - ! Transform eigenvectors with (A-B)^(1/2) if(.not.ga_create(mt_dbl,totpoles,totpoles,'tmp',0,totpoles, & g_tmp)) diff --git a/src/bsemol/bse_analytic_tda.F b/src/bsemol/bse_analytic_tda.F new file mode 100644 index 0000000000..029eb6dda9 --- /dev/null +++ b/src/bsemol/bse_analytic_tda.F @@ -0,0 +1,119 @@ + subroutine bse_analytic_tda(pars,wia,omega,maxpoles,totpoles, + $ singlet) + implicit none +#include "mafdecls.fh" +#include "errquit.fh" +#include "global.fh" +#include "bse.fh" +#ifdef SCALAPACK +#include "scaleig.fh" + integer ga_cholesky,ga_llt_i + external ga_cholesky,ga_llt_i +#endif + + type(bse_params_t) :: pars + + integer maxpoles, totpoles + double precision wia(maxpoles,pars%ipol) + double precision omega(totpoles) + logical singlet + + character(*),parameter :: pname = 'bse_analytic_tda: ' + + integer ilo,ihi,jlo,jhi + integer nocc2(3),nvir2(3) + integer ipole,imo,jmo,amo,bmo + integer k,l,kproc,lproc,klocal,llocal,kglobal,lglobal + integer g_diag, g_tmp, g_amb, ktmp, ld, info, i, j + integer modtn,tovern + double precision alpha + +#ifdef USE_OPENMP + call util_blas_set_num_threads(pars%iMaxthreads) +#endif + + if(.not.ga_create(mt_dbl,totpoles,totpoles,'apb',0,totpoles, + & pars%g_apb)) + & call errquit(pname//'could not create A+B GA',0,GA_ERR) + if (.not.nga_create(mt_dbl,1,totpoles,'wia',0,g_diag)) + $ call errquit(pname//'could not create wia',0,GA_ERR) + + if (pars%ipol.eq.1) then + nocc2(1) = pars%nocc(1)**2 + nvir2(1) = pars%nvir(1)**2 + + call ga_distribution(pars%g_apb,pars%me,ilo,ihi,jlo,jhi) + if (jhi.ge.jlo) call nga_put(g_diag,jlo,jhi,wia,jhi-jlo+1) + + ! Screened Coulomb contribution + if(.not.ga_create(mt_dbl,nocc2(1),nvir2(1),'tmp',nocc2(1),0, + & g_tmp)) + & call errquit(pname//'could not create temp GA',1,GA_ERR) + call ga_dgemm('t','n',nocc2(1),nvir2(1),pars%nri,-1d0, + % pars%g_erioo(1),pars%g_erivv(1),0d0,g_tmp) + call ga_distribution(g_tmp,pars%me,ilo,ihi,jlo,jhi) + call ga_access(g_tmp,ilo,ihi,jlo,jhi,ktmp,ld) + + ! Need to take care about some crazy indexing + tovern = totpoles/pars%nprocs + modtn = mod(totpoles,pars%nprocs) + do i=1,nocc2(1) + jmo = (i-1)/pars%nocc(1) + 1 + imo = i - (jmo-1)*pars%nocc(1) + do j=jlo,jhi + bmo = (j-1)/pars%nvir(1) + 1 + amo = j - (bmo-1)*pars%nvir(1) + k = amo + (imo-1)*pars%nvir(1) + l = bmo + (jmo-1)*pars%nvir(1) + kproc = mod(k-1,pars%nprocs) + lproc = mod(l-1,pars%nprocs) + klocal = (k-1)/pars%nprocs + 1 + llocal = (l-1)/pars%nprocs + 1 + kglobal = klocal + kproc*tovern + min(kproc, modtn) + lglobal = llocal + lproc*tovern + min(lproc, modtn) + call ga_put(pars%g_apb,kglobal,kglobal,lglobal,lglobal, + & dbl_mb(ktmp+(j-jlo)*nocc2(1)+i-1),1) + enddo + enddo + call ga_release(g_tmp,ilo,ihi,jlo,jhi) + if (.not.ga_destroy(g_tmp)) + & call errquit(pname//'could not destroy temp GA',1,GA_ERR) + + ! Hartree contribution + if (singlet) then + call ga_dgemm('t','n',totpoles,totpoles,pars%nri,2d0, + & pars%g_eriov(1),pars%g_eriov(1),1d0,pars%g_apb) + endif + + ! Diagonal term + call ga_add_diagonal(pars%g_apb,g_diag) + + ! Destroy ERIs + if (.not.ga_destroy(pars%g_erioo(1))) + & call errquit(pname//'could not destroy erioo GA',1,GA_ERR) + if (.not.ga_destroy(pars%g_erivv(1))) + & call errquit(pname//'could not destroy erivv GA',1,GA_ERR) + if (.not.ga_destroy(pars%g_eriov(1))) + & call errquit(pname//'could not destroy eriov GA',1,GA_ERR) + + + ! Diagonalize (A+B) matrix +#ifdef PARALLEL_DIAG +#ifdef SCALAPACK + call dft_scaleig(pars%g_apb,pars%g_apb,omega,0) +#else + call ga_diag_std(pars%g_apb,pars%g_apb,omega) +#endif +#else + call ga_diag_std_seq(pars%g_apb,pars%g_apb,omega) +#endif + + ! Destroy unused GAs + if (.not.ga_destroy(g_diag)) + & call errquit(pname//'could not destroy temp GA',1,GA_ERR) + + else + + endif + + end subroutine diff --git a/src/bsemol/bse_defaults.F b/src/bsemol/bse_defaults.F index dc62ab0eb0..226718a77b 100644 --- a/src/bsemol/bse_defaults.F +++ b/src/bsemol/bse_defaults.F @@ -14,6 +14,12 @@ & rtdb_put(rtdb,'bse:lanczos',mt_log,1,.false.) .and. & rtdb_put(rtdb,'bse:tda',mt_log,1,.false.) .and. & rtdb_put(rtdb,'bse:ncap',mt_log,1,.false.) .and. + & rtdb_put(rtdb,'bse:singlet',mt_log,1,.true.) .and. + & rtdb_put(rtdb,'bse:triplet',mt_log,1,.false.) .and. + & rtdb_put(rtdb,'bse:dowindow',mt_log,1,.false.) .and. + & rtdb_put(rtdb,'bse:elower',mt_dbl,1,0d0) .and. + & rtdb_put(rtdb,'bse:upper',mt_dbl,1,0d0) .and. + & rtdb_put(rtdb,'bse:nroots',mt_int,1,10) .and. & rtdb_put(rtdb,'bse:init',mt_log,1,.true.) if (.not.ok) diff --git a/src/bsemol/bse_dipole.F b/src/bsemol/bse_dipole.F new file mode 100644 index 0000000000..421414fb97 --- /dev/null +++ b/src/bsemol/bse_dipole.F @@ -0,0 +1,85 @@ + subroutine bse_dipole(pars) + implicit none +#include "global.fh" +#include "mafdecls.fh" +#include "errquit.fh" +#include "bse.fh" +#include "cdft.fh" + type(bse_params_t) :: pars + + character(*),parameter :: pname = 'bse_dipole: ' + + integer ixyz,isp,ilo,ihi,jlo,jhi,amo,imo + integer g_dip(3),handle(3),g_tmp1,g_tmp2 + integer k,klocal,kproc,kglobal,tovern,modtn,ktmp,ld + integer ga_create_atom_blocked + external ga_create_atom_blocked + + ! Create dipole matrices + do ixyz=1,3 + g_dip(ixyz) = ga_create_atom_blocked(geom,ao_bas_han,'Dipoles') + if(.not.nga_create(mt_dbl,1,pars%npoles(1),'dip a',0, + & pars%g_dipa(ixyz))) + & call errquit(pname//'could not create g_dipa',0,GA_ERR) + if (pars%ipol.gt.1) then + if(.not.nga_create(mt_dbl,1,pars%npoles(2),'dip b',0, + & pars%g_dipb(ixyz))) + & call errquit(pname//'could not create g_dipb',0,GA_ERR) + endif + enddo + + ! Compute dipole AO integrals for the 3 components + call int_dip_ga(ao_bas_han, ao_bas_han, g_dip(1), g_dip(2), + & g_dip(3)) + + ! Transform dipole matrices + do isp=1,pars%ipol + if (.not.ga_create(mt_dbl,nbf_ao,pars%nocc(isp),'tmp1', + & 0,pars%nocc(isp),g_tmp1)) + & call errquit(pname//'could not create g_tmp1',0,GA_ERR) + if (.not.ga_create(mt_dbl,pars%nvir(isp),pars%nocc(isp),'tmp2', + & pars%nvir(isp),0,g_tmp2)) + & call errquit(pname//'could not create g_tmp2',0,GA_ERR) + + if (isp.eq.1) then + handle(:) = pars%g_dipa(:) + else + handle(:) = pars%g_dipb(:) + endif + + tovern = pars%npoles(isp)/pars%nprocs + modtn = mod(pars%npoles(isp),pars%nprocs) + call ga_distribution(g_tmp2,pars%me,ilo,ihi,jlo,jhi) + do ixyz=1,3 + call ga_matmul_patch('n','n',1d0,0d0, + & g_dip(ixyz),1,nbf_ao,1,nbf_ao, + & pars%g_movecs(isp),1,nbf_ao,1,pars%nocc(isp), + & g_tmp1,1,nbf_ao,1,pars%nocc(isp)) + call ga_matmul_patch('t','n',1d0,0d0, + & pars%g_movecs(isp),1,pars%nvir(isp),1,nbf_ao, + & g_tmp1,1,nbf_ao,1,pars%nocc(isp), + & g_tmp2,1,pars%nvir(isp),1,pars%nocc(isp)) + + call ga_access(g_tmp2,ilo,ihi,jlo,jhi,ktmp,ld) + do imo=jlo,jhi + do amo=ilo,ihi + k = amo + (imo-1)*pars%nvir(isp) + kproc = mod(k-1,pars%nprocs) + klocal = (k-1)/pars%nprocs + 1 + kglobal = klocal + tovern*kproc + min(kproc,modtn) + call nga_put(handle(ixyz),kglobal,kglobal, + & dbl_mb(ktmp+(imo-1)*(ihi-ilo+1)+amo-ilo),1) + enddo + enddo + call ga_release(g_tmp2,ilo,ihi,jlo,jhi) + enddo + if ((.not.ga_destroy(g_tmp1)).and.(.not.ga_destroy(g_tmp2))) + & call errquit(pname//'could not destroy GAs',0,GA_ERR) + enddo + + do ixyz=1,3 + if (.not.ga_destroy(g_dip(ixyz))) + & call errquit(pname//'could not destroy GAs',0,GA_ERR) + enddo + + end diff --git a/src/bsemol/bse_finalize.F b/src/bsemol/bse_finalize.F index dd48873aef..a0b0299cb5 100644 --- a/src/bsemol/bse_finalize.F +++ b/src/bsemol/bse_finalize.F @@ -17,12 +17,19 @@ do i=1,ipol if (.not.ga_destroy(pars%g_movecs(i))) $ call errquit(pname//'failed to destroy g_movecs',0,GA_ERR) -c if (pars%npoles(i).gt.0) then -c if (.not.ga_destroy(pars%g_eriov(i))) -c $ call errquit(pname//'failed to destroy g_moeris',0,GA_ERR) -c endif enddo + do i=1,3 + if(.not.ga_destroy(pars%g_dipa(i))) + $ call errquit(pname//'failed to destroy g_dipa',0,GA_ERR) + enddo + if (ipol.gt.1) then + do i=1,3 + if(.not.ga_destroy(pars%g_dipb(i))) + $ call errquit(pname//'failed to destroy g_dipa',0,GA_ERR) + enddo + endif + ok = ma_chop_stack(pars%l_mf_evals) if (.not.ok) & call errquit(pname//'could not chop stack',0,MA_ERR) diff --git a/src/bsemol/bse_init.F b/src/bsemol/bse_init.F index fd893ffe4c..f962bf7f06 100644 --- a/src/bsemol/bse_init.F +++ b/src/bsemol/bse_init.F @@ -78,8 +78,32 @@ c & (pars%lanczos.and.(pars%analytic.or.pars%davidson))) & call errquit(pname//'more than one method',4,INPUT_ERR) + ! Check wether we are computing singlet or triplet excitations + if(.not.rtdb_get(pars%rtdb,'bse:singlet',mt_log,1,pars%singlet)) + & call errquit(pname//'could not read from rtdb',1,RTDB_ERR) + if(.not.rtdb_get(pars%rtdb,'bse:triplet',mt_log,1,pars%triplet)) + & call errquit(pname//'could not read from rtdb',1,RTDB_ERR) + if((pars%singlet .and. pars%triplet).and.(pars%analytic)) then + call errquit(pname//'select either triplet or singlet',0, + & INPUT_ERR) + endif + ! Check if we are doing Tamm-Dancoff approximation - if (.not.rtdb_get(pars%rtdb,'bse:tda',mt_log,1,pars%tda)) + if(.not.rtdb_get(pars%rtdb,'bse:tda',mt_log,1,pars%tda)) + & call errquit(pname//'could not read from rtdb',5,RTDB_ERR) + + ! Check if Output is restricted to an energy window + if(.not.rtdb_get(pars%rtdb,'bse:dowindow',mt_log,1,pars%dowindow)) + & call errquit(pname//'could not read from rtdb',5,RTDB_ERR) + if(pars%dowindow) then + if(.not.rtdb_get(pars%rtdb,'bse:elower',mt_dbl,1,pars%elower)) + & call errquit(pname//'could not read from rtdb',5,RTDB_ERR) + if(.not.rtdb_get(pars%rtdb,'bse:eupper',mt_dbl,1,pars%eupper)) + & call errquit(pname//'could not read from rtdb',5,RTDB_ERR) + endif + + ! Number of roots + if(.not.rtdb_get(pars%rtdb,'bse:nroots',mt_int,1,pars%nroots)) & call errquit(pname//'could not read from rtdb',5,RTDB_ERR) ! Get some info from DFT module @@ -188,13 +212,14 @@ c aq2 = 2d0/9d0 * (1d0-zeta)**2 do i=1,ipol ehomo = dbl_mb(pars%k_mf_evals+(i-1)*nbf_ao+pars%nocc(i)-1) - vxddp = -0.5d0*aq2*(1d0 + sqrt(1d0-4d0*ehomo/aq2)) - vxddm = -0.5d0*aq2*(1d0 - sqrt(1d0-4d0*ehomo/aq2)) + vxddp = -0.5d0*aq2*(1d0 + sqrt(1d0-4d0*ehomo/aq2))*0.8d0 + vxddm = -0.5d0*aq2*(1d0 - sqrt(1d0-4d0*ehomo/aq2))*0.8d0 pars%vxddp(i) = vxddp pars%vxddm(i) = vxddm if (pars%me.eq.0) then - write(luout,1001) ' v_x^{DD+}: ',vxddp*ha2ev, ' eV' - write(luout,1001) ' v_x^(DD-): ',vxddm*ha2ev, ' eV' + write(luout,*) ' Shifting eigvals by 80% NCAP DD' + write(luout,1001) ' v_x^{DD-}: ',vxddp*ha2ev, ' eV' + write(luout,1001) ' v_x^(DD+): ',vxddm*ha2ev, ' eV' endif call dcopy(nbf_ao,dbl_mb(pars%k_mf_evals+(i-1)*nbf_ao),1, & dbl_mb(pars%k_gw_evals+(i-1)*nbf_ao),1) @@ -229,6 +254,9 @@ c call schwarz_init(geom, ao_bas_han) call scf_get_fock_param(pars%rtdb, pars%tol2e) + ! Get dipole in CIS basis + call bse_dipole(pars) + !Clean previous grid call grid_cleanup(.false.) diff --git a/src/bsemol/bse_input.F b/src/bsemol/bse_input.F index 716c48739a..457cfea91b 100644 --- a/src/bsemol/bse_input.F +++ b/src/bsemol/bse_input.F @@ -18,12 +18,13 @@ c character*255 test character(*), parameter :: pname = "bse_input: " c - integer, parameter :: num_dirs = 4 + integer, parameter :: num_dirs = 8 character*11 dirs(num_dirs) - data dirs / 'method', 'tda', 'ncap', 'end' / + data dirs / 'method', 'tda', 'ncap', 'singlet', 'triplet', + & 'outwin', 'nroots', 'end' / c - integer ind, qpiter,eviter,noqp(2),nvqp(2), ngl, first - double precision threshold, eta + integer ind, itemp + double precision threshold, eta, temp double precision, parameter :: ha2ev = 27.211386245988d0 c c Set default values and put them in RTDB @@ -44,7 +45,7 @@ c call errquit(pname//'unknown directive', 0, INPUT_ERR) end if c - goto (100, 200, 300, 9999) ind + goto (100, 200, 300, 400, 500, 600, 700, 9999) ind call errquit(pname//'unknown directive',ind, INPUT_ERR) c @@ -85,6 +86,45 @@ c & call errquit(pname//'rtdb_put failed', 300, RTDB_ERR) goto 10 c +c compute singlet transitions +c + 400 if (.not.rtdb_put(rtdb, 'bse:singlet', mt_log, 1, .true.)) + & call errquit(pname//'rtdb_put failed', 300, RTDB_ERR) + if (.not.rtdb_put(rtdb, 'bse:triplet', mt_log, 1, .false.)) + & call errquit(pname//'rtdb_put failed', 300, RTDB_ERR) + goto 10 +c +c compute singlet transitions +c + 500 if (.not.rtdb_put(rtdb, 'bse:singlet', mt_log, 1, .false.)) + & call errquit(pname//'rtdb_put failed', 300, RTDB_ERR) + if (.not.rtdb_put(rtdb, 'bse:triplet', mt_log, 1, .true.)) + & call errquit(pname//'rtdb_put failed', 300, RTDB_ERR) + goto 10 +c +c output energy window +c + 600 if (.not.rtdb_put(rtdb, 'bse:dowindow', mt_log, 1, .true.)) + & call errquit(pname//'rtdb_put failed', 600, RTDB_ERR) + if(.not.inp_f(temp)) + & call errquit(pname//'could not read lower energy',600,INPUT_ERR) + if(.not.rtdb_put(rtdb,'bse:elower',mt_dbl,1,temp)) + & call errquit(pname//'rtdb_put failed', 600, RTDB_ERR) + if(.not.inp_f(temp)) + & call errquit(pname//'could not read upper energy',600,INPUT_ERR) + if(.not.rtdb_put(rtdb,'bse:eupper',mt_dbl,1,temp)) + & call errquit(pname//'rtdb_put failed', 600, RTDB_ERR) + goto 10 +c +c number of roots +c + 700 if(.not.inp_i(itemp)) + & call errquit(pname//'could not read no. of roots',700,INPUT_ERR) + if(.not.rtdb_put(rtdb,'bse:eupper',mt_int,1,itemp)) + & call errquit(pname//'rtdb_put failed', 700, RTDB_ERR) + goto 10 + +c c normal termination c 9999 return diff --git a/src/bsemol/bse_oscstr.F b/src/bsemol/bse_oscstr.F new file mode 100644 index 0000000000..1740c8302e --- /dev/null +++ b/src/bsemol/bse_oscstr.F @@ -0,0 +1,87 @@ + subroutine bse_oscstr(pars,omega,totpoles) + implicit none +#include "errquit.fh" +#include "stdio.fh" +#include "global.fh" +#include "mafdecls.fh" +#include "bse.fh" + type(bse_params_t) :: pars + + character(*), parameter :: pname = 'bse_oscstr: ' + character(len=7) :: key + + integer totpoles + double precision omega(*) + double precision oscstr,tmp,moms(3) + + integer llpole,ulpole,iroot,ixyz,ipole + + if (pars%singlet) then + key = 'singlet' + else + key = 'triplet' + endif + + ! Get the lower and upper limits + llpole = -1 + if (pars%dowindow) then + iroot = 0 + do ipole=1,totpoles + if (omega(ipole).lt.pars%elower) cycle + if (llpole.eq.-1) llpole = ipole + iroot = iroot + 1 + if (omega(ipole).gt.pars%eupper) then + ulpole = ipole - 1 + exit + endif + if (iroot.eq.pars%nroots) then + ulpole = ipole + exit + endif + enddo + else + llpole = 1 + ulpole = pars%nroots + endif + + + ! Compute the oscillator strengths around an energy window + iroot = 0 + do ipole=llpole,ulpole + iroot = iroot + 1 + oscstr = 0d0 + moms(:) = 0d0 + do ixyz=1,3 + moms(ixyz) = moms(ixyz) + nga_ddot_patch(pars%g_apb,'n', + & (/1,ipole/),(/pars%npoles(1),ipole/), + & pars%g_dipa(ixyz),'n',1,pars%npoles(1)) + oscstr = oscstr + moms(ixyz)**2 + enddo + if (pars%ipol.gt.1) then + do ixyz=1,3 + moms(ixyz) = moms(ixyz) + + & nga_ddot_patch(pars%g_apb,'n',(/pars%npoles(1)+1,ipole/), + & (/totpoles,ipole/), + & pars%g_dipb(ixyz),'n',1,pars%npoles(2)) + oscstr = oscstr + moms(ixyz)**2 + enddo + endif + oscstr = 2d0/3d0*omega(ipole)*oscstr + + if (pars%me.eq.0) then + write(luout,9110) iroot,key,'a1 ',omega(ipole), + & omega(ipole)*ha2ev + write(luout,9200) moms(1),moms(2),moms(3),oscstr + endif + enddo + + 9110 format(2x, + 1 '--------------------------------------------------------------- + 3-------------', + 2 /,2x,'Root',i4,1x,a7,1x,a4,f22.9,' a.u.',f22.4,' eV ',/,2x, + 1 '--------------------------------------------------------------- + 3-------------') + 9200 format(5x,'Transition Moments X',f9.5,' Y',f9.5,' Z',f9.5 + 5 ,/,5x,'Dipole Oscillator Strength',18x,f14.10) + + end diff --git a/src/bsemol/bse_wmn.F b/src/bsemol/bse_wmn.F index 0ed4e09241..61b7e488fb 100644 --- a/src/bsemol/bse_wmn.F +++ b/src/bsemol/bse_wmn.F @@ -15,17 +15,20 @@ integer isp,ilo,ihi,jlo,jhi,keri,ld do isp=1,ipol - ! copy and transform occ-vir block - if(.not.ga_duplicate(pars%g_eriov(isp),pars%g_wov(isp),'Wov')) - & call errquit(pname//'could not duplicate ERI OV',0,GA_ERR) - call ga_copy(pars%g_eriov(isp),pars%g_wov(isp)) - if(pars%mynpoles(isp).eq.0) cycle - call ga_access(pars%g_wov(isp),1,nri,pars%ovlo(isp), - & pars%ovhi(isp),keri,ld) - call dtfsm('n','l','l','n','n',nri,pars%mynpoles(isp),1d0, - & w,dbl_mb(keri),nri) - call ga_release_update(pars%g_wov(isp),1,nri,pars%ovlo(isp), - & pars%ovhi(isp)) + + if (.not.pars%tda) then + ! copy and transform occ-vir block + if(.not.ga_duplicate(pars%g_eriov(isp),pars%g_wov(isp),'Wov')) + & call errquit(pname//'could not duplicate ERI OV',0,GA_ERR) + call ga_copy(pars%g_eriov(isp),pars%g_wov(isp)) + if(pars%mynpoles(isp).eq.0) cycle + call ga_access(pars%g_wov(isp),1,nri,pars%ovlo(isp), + & pars%ovhi(isp),keri,ld) + call dtfsm('n','l','l','n','n',nri,pars%mynpoles(isp),1d0, + & w,dbl_mb(keri),nri) + call ga_release_update(pars%g_wov(isp),1,nri,pars%ovlo(isp), + & pars%ovhi(isp)) + endif ! transform occ-occ and vir-vir blocks in place call ga_distribution(pars%g_erioo(isp),pars%me,ilo,ihi,jlo,jhi)