diff --git a/src/bsemol/GNUmakefile b/src/bsemol/GNUmakefile index 3b3051af56..ca007bd758 100644 --- a/src/bsemol/GNUmakefile +++ b/src/bsemol/GNUmakefile @@ -2,7 +2,11 @@ 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_analytic_cis.o bse_analytic_tda.o \ - bse_dipole.o bse_oscstr.o bse_puteris.o bse_get_eia.o + bse_dipole.o bse_oscstr.o bse_puteris.o bse_get_eia.o \ + bse_davidson_diagonal.o bse_davidson_mvec.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 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 c2d8730a74..2b04c0d8f8 100644 --- a/src/bsemol/bse.fh +++ b/src/bsemol/bse.fh @@ -15,7 +15,7 @@ c integer :: ovlo(2),ovhi(2) integer :: vvlo(2),vvhi(2) integer :: llmo(2),nqp(2) - integer :: nroots + integer :: nroots,nmax double precision :: tol2e, eta c integer :: l_mf_evals, k_mf_evals @@ -27,7 +27,7 @@ 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, g_dipa(3), g_dipb(3) + integer :: g_wov(2), g_apb, g_dipa(3), g_dipb(3), g_amb integer :: g_quada(6), g_quadb(6) integer :: g_dipmaga(3), g_dipmagb(3) logical :: noio diff --git a/src/bsemol/bse_analytic.F b/src/bsemol/bse_analytic.F index 34704b7bbb..463c47e0f8 100644 --- a/src/bsemol/bse_analytic.F +++ b/src/bsemol/bse_analytic.F @@ -40,9 +40,6 @@ & dbl_mb(k_dia+lSize),pars%nocc(isp), & pars%nvir(isp),pars%ovlo(isp),pars%ovhi(isp)) enddo - if (pars%me.eq.0) then - call ma_print(dbl_mb(k_dia),pars%mynpoles(1),1,'Diff1') - endif ! obtain screened Coulomb matrix call bse_buildw(pars,dbl_mb(k_dia),dbl_mb(kW),pars%nmo,nri, @@ -61,9 +58,6 @@ & pars%nvir(isp),pars%ovlo(isp), & pars%ovhi(isp)) enddo - if (pars%me.eq.0) then - call ma_print(dbl_mb(k_dia),pars%mynpoles(1),1,'Diff2') - endif ! Build and diagonalize Casida-like matrix if (pars%tda) then diff --git a/src/bsemol/bse_analytic_cis.F b/src/bsemol/bse_analytic_cis.F index c8841fbde5..98d73f3432 100644 --- a/src/bsemol/bse_analytic_cis.F +++ b/src/bsemol/bse_analytic_cis.F @@ -230,7 +230,10 @@ if (.not.ga_destroy(pars%g_erivv(isp))) & call errquit(pname//'could not destroy erivv GA',177,GA_ERR) - call ga_copy(pars%g_apb,g_amb) + call ga_copy_patch('n',pars%g_apb,1+off1, + & pars%npoles(isp)+off1,1+off1,pars%npoles(isp)+off1, + & g_amb,1+off1,pars%npoles(isp)+off1,1+off1, + & pars%npoles(isp)+off1) ! Add Screened Coulomb contribution from B @@ -276,6 +279,7 @@ call ga_add_diagonal(pars%g_apb,g_diag) call ga_add_diagonal(g_amb,g_diag) + ! Destroy ERIs if (.not.ga_destroy(pars%g_eriov(1))) & call errquit(pname//'could not destroy eriov GA',275,GA_ERR) diff --git a/src/bsemol/bse_davidson.F b/src/bsemol/bse_davidson.F index ba4aa6b89a..6c43eb2afc 100644 --- a/src/bsemol/bse_davidson.F +++ b/src/bsemol/bse_davidson.F @@ -2,8 +2,308 @@ implicit none +#include "global.fh" +#include "errquit.fh" +#include "mafdecls.fh" #include "bse.fh" - +#include "cdft.fh" type(bse_params_t) :: pars + character(*), parameter :: pname = 'bse_analytic: ' + + integer lSize,isp + integer l_omega, k_omega + integer l_dia,k_dia + integer lW,kW + integer maxpoles,maxnpoles,totnpoles,nri + + integer indeces(pars%nroots),ntrials,g_trials,npoles,ndim + integer g_x,g_y,g_ym,g_s,g_ks,g_mks,g_w,g_kw,g_mkw,off1 + integer g_wia,g_tmp,info,lwork,kwork,idum,mdim + integer k_x,k_y,k_ym,k_s,k_ks,k_mks,k_w,k_kw,k_mkw,k_r,l_r + integer myld,vkv,l_vkv,myhi,mylo,k_trials,l_ksw,ksw,jdum + double precision minimum,norm(pars%nroots),rdum,factor + double precision,allocatable :: diagonal(:),diagonal2(:) + + integer ipole,jpole,itrial,imo,amo,idavidson + logical found,dohartree,stat,debug + double precision,external :: dasum + + debug = pars%me.eq.0 + debug = .false. + + ! synchronize parallel processes + call ga_sync() + + !initialize useful variables + nri = pars%nri + maxpoles = maxval(pars%mynpoles(1:pars%ipol)) + maxnpoles = maxval(pars%npoles(1:pars%ipol)) + totnpoles = maxnpoles*pars%ipol + npoles = sum(pars%npoles(1:pars%ipol)) + ntrials = pars%nroots + + + allocate(diagonal(npoles),diagonal2(npoles)) + + ! allocation + if(.not.ma_push_get(mt_dbl,totnpoles,'dia',l_dia,k_dia)) + & call errquit(pname//'failed to allocate Delta_ia',0,MA_ERR) + if(.not.ma_push_get(mt_dbl,npoles,'omega',l_omega,k_omega)) + & call errquit(pname//'failed to allocate omega',0,MA_ERR) + if(.not.ma_push_get(mt_dbl,npoles,'residual',l_r,k_r)) + & call errquit(pname//'failed to allocate residual',0,MA_ERR) + if(.not.ma_push_get(mt_dbl,ntrials*ntrials,'vkv',l_vkv,vkv)) + & call errquit(pname//'failed to allocate vkv',0,MA_ERR) + if(.not.ma_push_get(mt_dbl,nri**2,'W',lW,kW)) + & call errquit(pname//'failed to allocate W',0,MA_ERR) + + stat = ga_create(mt_dbl,npoles,ntrials,'x',0,ntrials,g_x).and. + $ ga_create(mt_dbl,npoles,ntrials,'y',0,ntrials,g_y).and. + $ ga_create(mt_dbl,npoles,ntrials,'ym',0,ntrials,g_ym).and. + $ ga_create(mt_dbl,npoles,ntrials,'s',0,ntrials,g_s).and. + $ ga_create(mt_dbl,npoles,ntrials,'ks',0,ntrials,g_ks).and. + $ ga_create(mt_dbl,npoles,ntrials,'mks',0,ntrials,g_mks).and. + $ ga_create(mt_dbl,npoles,ntrials,'w',0,ntrials,g_w).and. + $ ga_create(mt_dbl,npoles,ntrials,'kw',0,ntrials,g_kw).and. + $ ga_create(mt_dbl,npoles,ntrials,'mkw',0,ntrials,g_mkw).and. + $ ga_create(mt_dbl,npoles,ntrials,'t',0,ntrials,g_trials).and. + $ nga_create(mt_dbl,1,npoles,'wia',0,g_wia) + + if(.not.stat) call errquit('ga create failed',66,GA_ERR) + + ! obtain "old" eigenvalue difference to build W + do isp=1,pars%ipol + lSize = (isp-1)*maxpoles + call bse_get_eia(dbl_mb(pars%k_mf_evals+(isp-1)*nbf_ao), + & dbl_mb(k_dia+lSize),pars%nocc(isp), + & pars%nvir(isp),pars%ovlo(isp),pars%ovhi(isp)) + enddo + + ! obtain screened Coulomb matrix + call bse_buildw(pars,dbl_mb(k_dia),dbl_mb(kW),pars%nmo,nri, + & maxpoles,pars%ipol) + + ! transform ERIs using screened Coulomb matrix + call bse_wmn(pars,dbl_mb(kW),pars%nmo,nri,pars%ipol) + if(.not.ma_chop_stack(lW)) + & call errquit(pname//'failed to chop stack',0,MA_ERR) + + ! obtain "GW" eigenvalue difference to build Casida matrix + do isp=1,pars%ipol + lSize = (isp-1)*maxpoles + call bse_get_eia(dbl_mb(pars%k_gw_evals+(isp-1)*nbf_ao), + & dbl_mb(k_dia+lSize),pars%nocc(isp), + & pars%nvir(isp),pars%ovlo(isp), + & pars%ovhi(isp)) + enddo + + ! Get diagonal of Casida matrix + call bse_davidson_diagonal(pars,dbl_mb(k_dia),maxpoles) + + ! Get Eigenvalue differences + do isp=1,pars%ipol + off1 = (isp-1)*pars%npoles(1) + lSize = (isp-1)*maxpoles + call nga_put(g_wia,pars%ovlo(isp)+off1,pars%ovhi(isp)+off1, + $ dbl_mb(k_dia+lSize),pars%ovhi(isp)-pars%ovlo(isp)+1) + enddo + + ! Get guess vectos + if(pars%me.eq.0) + $ call bse_davidson_guess(pars,dbl_mb(k_dia),g_trials,npoles, + $ ntrials) + call ga_brdcst(mt_dbl,diagonal,ma_sizeof(mt_dbl,npoles,mt_byte),0) + + call ga_distribution(g_s,pars%me,mylo,myhi,idum,idum) + myld = myhi - mylo + 1 + + ! K-orthogonalization + call bse_davidson_kvec(pars,g_trials,g_wia,g_x,npoles, + $ ntrials,pars%singlet) + call ga_access(g_x,mylo,myhi,1,ntrials,k_x,myld) + call ga_access(g_trials,mylo,myhi,1,ntrials,k_trials,myld) + call bse_davidson_kortho(pars,dbl_mb(k_trials),dbl_mb(k_x), + $ dbl_mb(k_x),dbl_mb(vkv),dbl_mb(k_x),myld,ntrials,mdim,.true.) + call ga_release_update(g_trials,mylo,myhi,1,ntrials) + call ga_release_update(g_x,mylo,myhi,1,ntrials) + + ! Starting vectors + call bse_davidson_kvec(pars,g_x,g_wia,g_y,npoles, + $ ntrials,pars%singlet) + call bse_davidson_mvec(pars,g_y,g_wia,g_ym,npoles, + $ ntrials) + call ga_copy(g_x,g_s) + call ga_copy(g_y,g_ks) + call ga_copy(g_ym,g_mks) + + if (pars%me.eq.0) then + write(*,*) + write(*,*) "--------------------------------" + write(*,*) " Iter Ntrials Nconv Res " + write(*,*) "--------------------------------" + endif + + do idavidson=1,20 + + ! Get access to local chunks + call ga_access(g_s,mylo,myhi,1,ntrials,k_s,myld) + call ga_access(g_ks,mylo,myhi,1,ntrials,k_ks,myld) + call ga_access(g_mks,mylo,myhi,1,ntrials,k_mks,myld) + call ga_access(g_x,mylo,myhi,1,ntrials,k_x,myld) + call ga_access(g_y,mylo,myhi,1,ntrials,k_y,myld) + call ga_access(g_ym,mylo,myhi,1,ntrials,k_ym,myld) + + ! Build projected matrix + call dgemm('t','n',ntrials,ntrials,myld,1d0,dbl_mb(k_ks), + $ myld,dbl_mb(k_mks),myld,0d0,dbl_mb(vkv),ntrials) + call ga_dgop((/194/),dbl_mb(vkv),ntrials**2,'+') + if (debug) call ma_print(dbl_mb(vkv),ntrials,ntrials,'Proj') + + ! Diagonalize projected matrix + call dsyev('v','l',ntrials,dbl_mb(vkv),ntrials,dbl_mb(k_omega), + $ rdum,-1,info) + if(.not.ma_push_get(mt_dbl,int(rdum),'work',lwork,kwork)) + & call errquit(pname//'failed to allocate work',199,MA_ERR) + call dsyev('v','l',ntrials,dbl_mb(vkv),ntrials,dbl_mb(k_omega), + $ dbl_mb(kwork),int(rdum),info) + if(.not.ma_chop_stack(lwork)) + & call errquit(pname//'failed to allocate work',192,MA_ERR) + if (info.ne.0) + $ call errquit(pname//'Diagonalization failed',206,0) + if (debug) call ma_print(dbl_mb(vkv),ntrials,ntrials,'Eigen') + + + ! Project back vectors + call dgemm('n','n',myld,ntrials,ntrials,1d0,dbl_mb(k_s),myld, + $ dbl_mb(vkv),ntrials,0d0,dbl_mb(k_x),myld) + call dgemm('n','n',myld,ntrials,ntrials,1d0,dbl_mb(k_ks),myld, + $ dbl_mb(vkv),ntrials,0d0,dbl_mb(k_y),myld) + call dgemm('n','n',myld,ntrials,ntrials,1d0,dbl_mb(k_mks),myld, + $ dbl_mb(vkv),ntrials,0d0,dbl_mb(k_ym),myld) + + ! Get residual + call ga_access(g_w,mylo,myhi,1,pars%nroots,k_w,myld) + +#if 1 + call bse_davidson_residual(pars%nroots,ntrials,dbl_mb(k_omega), + $ dbl_mb(k_x),dbl_mb(k_ym),dbl_mb(k_r), + $ dbl_mb(k_dia+mylo-1),dbl_mb(k_w),myld,ndim,norm) + + ! Print info + if (pars%me.eq.0) then + write(*,9010) idavidson,ntrials,pars%nroots-ndim,maxval(norm) + endif + 9010 format(I5,3X,I5,3X,I5,G16.5) +#else + minimum = 0d0 + ndim = 0 + do itrial=1,pars%nroots + call dcopy(myld,dbl_mb(k_ym+(itrial-1)*myld),1,dbl_mb(k_r),1) + call daxpy(myld,-dbl_mb(k_omega+itrial-1), + $ dbl_mb(k_x+(itrial-1)*myld),1,dbl_mb(k_r),1) + norm = maxval(abs(dbl_mb(k_r:k_r+myhi-mylo))) + call ga_dgop((/229/),norm,1,'maxabs') + write(*,*) " ",itrial,norm(1), + $ dsqrt(dbl_mb(k_omega+itrial-1))*27.2114d0 + if (norm(1).gt.1d-5) then + ndim = ndim + 1 + do ipole=mylo,myhi + factor = (dbl_mb(k_dia+ipole-1)-dbl_mb(k_omega+itrial-1)) + dbl_mb(k_r+ipole-mylo) = dbl_mb(k_r+ipole-mylo)/factor + enddo + call ga_put(g_w,mylo,myhi,ndim,ndim,dbl_mb(k_r),myld) + endif + enddo +#endif + + call ga_release_update(g_w,mylo,myhi,1,pars%nroots) + call ga_release_update(g_x,mylo,myhi,1,ntrials) + call ga_release_update(g_y,mylo,myhi,1,ntrials) + call ga_release_update(g_ym,mylo,myhi,1,ntrials) + call ga_release_update(g_s,mylo,myhi,1,ntrials) + call ga_release_update(g_ks,mylo,myhi,1,ntrials) + call ga_release_update(g_mks,mylo,myhi,1,ntrials) + + ! All residuals are below the threshold + if (ndim.eq.0) exit + + ! Restart calculation if dimensions exceeded + if (ndim+ntrials.gt.pars%nmax) then + call bse_davidson_restart(ntrials,pars%nroots,npoles, + $ g_x,g_y,g_ym,g_s,g_ks,g_mks,l_vkv,vkv) + cycle + endif + + ! Project-out current vectors + call ga_access(g_ks,mylo,myhi,1,ntrials,k_ks,myld) + call ga_access(g_s,mylo,myhi,1,ntrials,k_s,myld) + call ga_access(g_w,mylo,myhi,1,pars%nroots,k_w,myld) + if(.not.ma_push_get(mt_dbl,ndim*ntrials,'ksw',l_ksw,ksw)) + $ call errquit(pname//'ksw allocation failed',264,MA_ERR) + call dgemm('t','n',ntrials,ndim,myld,1d0,dbl_mb(k_ks),myld, + $ dbl_mb(k_w),myld,0d0,dbl_mb(ksw),ntrials) + call ga_dgop((/268/),dbl_mb(ksw),ntrials*ndim,'+') + call dgemm('n','n',myld,ndim,ntrials,-1d0,dbl_mb(k_s),myld, + $ dbl_mb(ksw),ntrials,1d0,dbl_mb(k_w),myld) + if(.not.ma_chop_stack(l_ksw)) + $ call errquit(pname//'ksw deallocation failed',272,MA_ERR) + call ga_release(g_ks,mylo,myhi,1,ntrials) + call ga_release(g_s,mylo,myhi,1,ntrials) + call ga_release_update(g_w,mylo,myhi,1,pars%nroots) + + + ! K-orthogonalization + call bse_davidson_kvec(pars,g_w,g_wia,g_kw,npoles, + $ ndim,pars%singlet) + call bse_davidson_mvec(pars,g_kw,g_wia,g_mkw,npoles, + $ ndim) + call ga_access(g_w,mylo,myhi,1,pars%nroots,k_w,myld) + call ga_access(g_kw,mylo,myhi,1,pars%nroots,k_kw,myld) + call ga_access(g_mkw,mylo,myhi,1,pars%nroots,k_mkw,myld) + call ga_access(g_x,mylo,myhi,1,mdim,k_x,myld) + call bse_davidson_kortho(pars,dbl_mb(k_w),dbl_mb(k_kw), + $ dbl_mb(k_mkw),dbl_mb(vkv),dbl_mb(k_x),myld,ndim,mdim, + $ .false.) + call ga_release_update(g_x,mylo,myhi,1,mdim) + call ga_release_update(g_w,mylo,myhi,1,pars%nroots) + call ga_release_update(g_kw,mylo,myhi,1,pars%nroots) + call ga_release_update(g_mkw,mylo,myhi,1,pars%nroots) + + ! All trial vectors were linearly-dependent + ! restart calculation with current best guess + if (mdim.eq.0) then + call bse_davidson_restart(ntrials,pars%nroots,npoles, + $ g_x,g_y,g_ym,g_s,g_ks,g_mks,l_vkv,vkv) + cycle + endif + + ndim = mdim + + ! Extend the size of the current vectors + call bse_davidson_newvecs(ntrials,ndim,npoles,g_x,g_y,g_ym, + $ g_s,g_ks,g_mks,g_w,g_kw,g_mkw, + $ l_vkv,vkv) + + enddo + + stat = ga_destroy(g_x).and.ga_destroy(g_y).and.ga_destroy(g_ym) + stat = stat.and.ga_destroy(g_s).and.ga_destroy(g_ks) + stat = stat.and.ga_destroy(g_mks).and.ga_destroy(g_w) + stat = stat.and.ga_destroy(g_kw).and.ga_destroy(g_mkw) + stat = stat.and.ga_destroy(g_trials).and.ga_destroy(g_wia) + stat = stat.and.ga_destroy(pars%g_apb).and.ga_destroy(pars%g_amb) + do isp=1,pars%ipol + stat = stat.and.ga_destroy(pars%g_erioo(isp)) + stat = stat.and.ga_destroy(pars%g_eriov(isp)) + stat = stat.and.ga_destroy(pars%g_erivv(isp)) + stat = stat.and.ga_destroy(pars%g_wov(isp)) + enddo + if(.not.stat) call errquit('deallocation failed',0,GA_ERR) + stat = ma_chop_stack(l_dia) + if(.not.stat) call errquit('deallocation failed',0,MA_ERR) end subroutine bse_davidson + + + + diff --git a/src/bsemol/bse_davidson_diagonal.F b/src/bsemol/bse_davidson_diagonal.F new file mode 100644 index 0000000000..b8024d043d --- /dev/null +++ b/src/bsemol/bse_davidson_diagonal.F @@ -0,0 +1,98 @@ + subroutine bse_davidson_diagonal(pars,wia,maxpoles) + implicit none +#include "mafdecls.fh" +#include "errquit.fh" +#include "global.fh" +#include "bse.fh" +#include "util.fh" + + type(bse_params_t) :: pars + + integer maxpoles, npoles + double precision wia(maxpoles,pars%ipol) + + character(*),parameter :: pname = 'bse_casida: ' + + integer isp,jsp,nri,off1,off2 + integer ilo,ihi,jlo,jhi + integer nocc2(2),nvir2(2) + 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 ij,jb,llia,ulia,ltmp2,ltmp3,ktmp2,ktmp3 + double precision temp(maxpoles),factor + + double precision, external :: ddot + +#ifdef USE_OPENMP + call util_blas_set_num_threads(pars%iMaxthreads) +#endif + + nri = pars%nri + npoles = pars%npoles(1) + if (pars%ipol.gt.1) npoles = npoles + pars%npoles(2) + if (.not.nga_create(mt_dbl,1,npoles,'apb',0,pars%g_apb)) + $ call errquit(pname//'could not create wia',0,GA_ERR) + if (.not.nga_create(mt_dbl,1,npoles,'amb',0,pars%g_amb)) + $ call errquit(pname//'could not create wia',0,GA_ERR) + if(.not.ma_push_get(mt_dbl,nri,'dia',ltmp2,ktmp2)) + & call errquit(pname//'failed to allocate Delta_ia',0,MA_ERR) + if(.not.ma_push_get(mt_dbl,nri,'dia',ltmp3,ktmp3)) + & call errquit(pname//'failed to allocate Delta_ia',0,MA_ERR) + + factor = 2d0 + if ((pars%singlet) .and. (pars%ipol.eq.1)) factor = 4d0 + if ((pars%triplet) .and. (pars%ipol.eq.1)) factor = 0d0 + if (pars%tda) factor = factor/2d0 + + do isp=1,pars%ipol + off1 = (isp-1)*pars%npoles(1) + + ! Eigenvalue differences + call ga_distribution(pars%g_eriov(isp),pars%me,ilo,ihi,jlo,jhi) + call dcopy(jhi-jlo+1,wia(1,isp),1,temp,1) + + ! Screened Coulomb from A + do ipole=jlo,jhi + imo = (ipole-1)/pars%nvir(isp) + 1 + amo = ipole - (imo-1)*pars%nvir(isp) + imo = imo + (imo-1)*pars%nocc(isp) + amo = amo + (amo-1)*pars%nvir(isp) + call ga_get(pars%g_erioo(isp),1,nri,imo,imo,dbl_mb(ktmp2),nri) + call ga_get(pars%g_erivv(isp),1,nri,amo,amo,dbl_mb(ktmp3),nri) + temp(ipole-jlo+1) = temp(ipole-jlo+1) - + $ ddot(nri,dbl_mb(ktmp2),1,dbl_mb(ktmp3),1) + enddo + call nga_put(pars%g_apb,jlo+off1,jhi+off1,temp,jhi-jlo+1) + call nga_put(pars%g_amb,jlo+off1,jhi+off1,temp,jhi-jlo+1) + + ! Screened Coulomb from B + if (.not.pars%tda) then + call ga_access(pars%g_wov(isp),ilo,ihi,jlo,jhi,ktmp,ld) + do ipole=jlo,jhi + temp(ipole-jlo+1) = ddot(nri,dbl_mb(ktmp+(ipole-jlo)*nri),1, + $ dbl_mb(ktmp+(ipole-jlo)*nri),1) + enddo + call nga_acc(pars%g_apb,jlo+off1,jhi+off1,temp,jhi-jlo+1,-1d0) + call nga_acc(pars%g_amb,jlo+off1,jhi+off1,temp,jhi-jlo+1,1d0) + call ga_release(pars%g_wov(isp),ilo,ihi,jlo,jhi) + endif + + ! Hartree term + call ga_access(pars%g_eriov(isp),ilo,ihi,jlo,jhi,ktmp,ld) + do ipole=jlo,jhi + temp(ipole-jlo+1) = + $ factor*ddot(nri,dbl_mb(ktmp+(ipole-jlo)*nri),1, + $ dbl_mb(ktmp+(ipole-jlo)*nri),1) + enddo + call nga_acc(pars%g_apb,jlo+off1,jhi+off1,temp,jhi-jlo+1,1d0) + call ga_release(pars%g_eriov(isp),ilo,ihi,jlo,jhi) + + enddo + + if(.not.ma_chop_stack(ltmp2)) + $ call errquit('bse_davidson_diagonal: chop stack',0,MA_ERR) + + call ga_sync() + + end subroutine diff --git a/src/bsemol/bse_davidson_guess.F b/src/bsemol/bse_davidson_guess.F new file mode 100644 index 0000000000..c4a3d33ead --- /dev/null +++ b/src/bsemol/bse_davidson_guess.F @@ -0,0 +1,96 @@ + subroutine bse_davidson_guess(pars,diagonal,g_trials,npoles, + $ ntrials) + + implicit none + +#include "bse.fh" +#include "stdio.fh" +#include "mafdecls.fh" +#include "cdft.fh" + + type(bse_params_t) :: pars + integer npoles,ntrials,g_trials + double precision diagonal(npoles) + double precision amb(npoles) + + character*4 oname,vname + logical found + integer amo,imo,itrial,ipole,jpole,indeces(ntrials) + integer isym,jsym,i,j,pol + double precision minimum + + call nga_get(pars%g_apb,1,npoles,diagonal,npoles) + call nga_get(pars%g_amb,1,npoles,amb,npoles) + + diagonal(:) = diagonal(:)*amb(:) + indeces(:) = -1 + + ! Get the indeces of the smallest diagonal elements + do itrial=1,ntrials + minimum = 1d6 + do ipole=1,npoles + if (abs(diagonal(ipole)).lt.minimum) then + found = .true. + do jpole=1,itrial + if (ipole.eq.indeces(jpole)) then + found = .false. + exit + endif + enddo + if (found) then + indeces(itrial) = ipole + minimum = abs(diagonal(ipole)) + endif + endif + enddo + enddo + + ! Get symmetry of ground state + isym = 1 + if (pars%ipol.eq.2) then + do i=1,pars%ipol + do j=1,pars%nocc(i) + isym = ieor(isym-1,int_mb(pars%k_irs(i)+j-1)-1)+1 + enddo + enddo + endif + call sym_irrepname(geom, isym, oname) + + + ! Start the guess vector + write(LuOut,*) + write(LuOut,9000) ntrials + write(LuOut,9020) + write(LuOut,9010) + write(LuOut,9020) + do itrial=1,ntrials + ipole = indeces(itrial) + call ga_put(g_trials,ipole,ipole,itrial,itrial,1d0,1) + + ! Print information + if (ipole.le.pars%npoles(1)) then + imo = (ipole-1)/pars%nvir(1) + 1 + amo = ipole - (imo-1)*pars%nvir(1) + pars%nocc(1) + jsym = ieor(isym-1,int_mb(pars%k_irs(1)+imo-1)-1)+1 + jsym = ieor(jsym-1,int_mb(pars%k_irs(1)+amo-1)-1)+1 + pol = 1 + else + imo = (ipole - pars%npoles(1) - 1)/pars%nvir(2) + 1 + amo = ipole - pars%npoles(1) - (imo-1)*pars%nvir(2) + + $ pars%nocc(2) + jsym = ieor(isym-1,int_mb(pars%k_irs(2)+imo-1)-1)+1 + jsym = ieor(jsym-1,int_mb(pars%k_irs(2)+amo-1)-1)+1 + pol = 2 + endif + call sym_irrepname(geom, jsym, vname) + write(luout,9030) itrial,pol,imo,amo,vname, + $ dsqrt(diagonal(ipole))*27.2114d0 + enddo + + 9000 format(i5,' smallest diagonals (eV) ') + 9010 format(' No. Spin Occ Vir Irrep Diag') + 9020 format('------------------------------------') + 9030 format(4i5,3x,a4,1x,f8.3) + + end subroutine + diff --git a/src/bsemol/bse_davidson_kortho.F b/src/bsemol/bse_davidson_kortho.F new file mode 100644 index 0000000000..821763a600 --- /dev/null +++ b/src/bsemol/bse_davidson_kortho.F @@ -0,0 +1,71 @@ + subroutine bse_davidson_kortho(pars,v,kv,mkv,vkv,work,npoles,n,m, + $ first) + + implicit none +#include "bse.fh" +#include "mafdecls.fh" +#include "errquit.fh" + + type(bse_params_t) :: pars + logical first + integer npoles,n,m + + integer l_work,k_work,info,i,j,k + + double precision rdum + double precision v(npoles,n),kv(npoles,n),mkv(npoles,n) + double precision vkv(n,n),eig(n),work(npoles,n) + + ! V**T K V + call dgemm('t','n',n,n,npoles,1d0,v,npoles,kv,npoles,0d0,vkv,n) + call ga_dgop((/14/),vkv,n*n,'+') + + ! Allocate work field + call dsyev('v','l',n,vkv,n,eig,rdum,-1,info) + if(.not.ma_push_get(mt_dbl,int(rdum),'work',l_work,k_work)) + $ call errquit('bse_kortho: failed to allocate work',199,MA_ERR) + + ! Diagonalize matrix + call dsyev('v','l',n,vkv,n,eig,dbl_mb(k_work),int(rdum),info) + if (info.ne.0) + $ call errquit("bse_davidson_kortho: diago failed",0,0) + + ! Deallocate field + if (.not.ma_chop_stack(l_work)) + $ call errquit('bse_kortho: failed to deallocate work',199,MA_ERR) + + ! Remove small eigenpairs + m = n + k = 1 + 100 continue + do i=k,m + if (eig(i).lt.1d-10) then + do j=i+1,n + eig(j-1) = eig(j) + vkv(:,j-1) = vkv(:,j) + enddo + m = m - 1 + k = i + goto 100 + endif + enddo + + ! Obtain VKV^(1/2) + do i=1,m + vkv(:,i) = vkv(:,i)/dsqrt(eig(i)) + enddo + + if (m.eq.0) return + + ! Obtain orthogonalized vectors + call dgemm('n','n',npoles,m,n,1d0,v,npoles,vkv,n,0d0,work,npoles) + call dcopy(npoles*m,work,1,v,1) + + if (first) return + call dgemm('n','n',npoles,m,n,1d0,kv,npoles,vkv,n,0d0,work,npoles) + call dcopy(npoles*m,work,1,kv,1) + call dgemm('n','n',npoles,m,n,1d0,mkv,npoles,vkv,n,0d0,work, + $ npoles) + call dcopy(npoles*m,work,1,mkv,1) + + end subroutine diff --git a/src/bsemol/bse_davidson_kvec.F b/src/bsemol/bse_davidson_kvec.F new file mode 100644 index 0000000000..90b256e463 --- /dev/null +++ b/src/bsemol/bse_davidson_kvec.F @@ -0,0 +1,158 @@ + subroutine bse_davidson_kvec(pars,vec,wia,kv,npoles,ntrials, + $ dohartree) + implicit none +#include "bse.fh" +#include "global.fh" +#include "mafdecls.fh" +#include "errquit.fh" + + type(bse_params_t) :: pars + + integer npoles,ntrials + integer vec,wia,kv,mv + + logical dohartree + + integer ipol,npoles1,npoles2,uli,lli,ipole + integer nocc(2),nvir(2),amo,imo,itrial,ula,lla,nri + + integer g_tmp1(2),g_tmp2(2),g_tmp3(2) + double precision temp(npoles),factor + + ipol = pars%ipol + nri = pars%nri + nocc = pars%nocc + nvir = pars%nvir + npoles1 = nocc(1)*nvir(1) + npoles2 = nocc(2)*nvir(2) + + call ga_copy(vec,kv) + call ga_scale_rows(kv,wia) + + if (.not.ga_create(mt_dbl,nvir(1),nocc(1),'g_tmp1',0,nocc(1), + $ g_tmp1)) + $ call errquit('bse_davidson_kvec: ga create failed',34,GA_ERR) + if (.not.ga_create(mt_dbl,nri,nocc(1),'g_tmp2',0,nocc(1), + $ g_tmp2)) + $ call errquit('bse_davidson_kvec: ga create failed',37,GA_ERR) + if (pars%ipol.gt.1) then + if (.not.ga_create(mt_dbl,nvir(2),nocc(2),'g_tmp1',0,nocc(2), + $ g_tmp1(2))) + $ call errquit('bse_davidson_kvec: ga create failed',48,GA_ERR) + if (.not.ga_create(mt_dbl,nri,nocc(2),'g_tmp2',0,nocc(2), + $ g_tmp2(2))) + $ call errquit('bse_davidson_kvec: ga create failed',51,GA_ERR) + endif + + do itrial=1,ntrials + ! Screened Coulomb A + call ga_copy_patch('n',vec,1,npoles1,itrial,itrial, + $ g_tmp1,1,nvir(1),1,nocc(1)) + do amo=1,nvir(1) + lla = (amo-1)*nvir(1) + 1 + ula = amo*nvir(1) + call ga_matmul_patch('n','n',1d0,0d0, + $ pars%g_erivv(1),1,nri,lla,ula, + $ g_tmp1,1,nvir(1),1,nocc(1), + $ g_tmp2,1,nri,1,nocc(1)) + do imo=1,nocc(1) + lli = (imo-1)*nocc(1) + 1 + uli = imo*nocc(1) + ipole = (imo-1)*nvir(1) + amo + temp(ipole) = ga_ddot_patch(g_tmp2,'n',1,nri,1,nocc(1), + $ pars%g_erioo(1),'n',1,nri,lli,uli) + enddo + enddo + if (ipol.gt.1) then + call ga_copy_patch('n',vec,1+npoles1,npoles,itrial,itrial, + $ g_tmp1(2),1,nvir(2),1,nocc(2)) + do amo=1,nvir(2) + lla = (amo-1)*nvir(2) + 1 + ula = amo*nvir(2) + call ga_matmul_patch('n','n',1d0,0d0, + $ pars%g_erivv(2),1,nri,lla,ula, + $ g_tmp1(2),1,nvir(2),1,nocc(2), + $ g_tmp2(2),1,nri,1,nocc(2)) + do imo=1,nocc(2) + lli = (imo-1)*nocc(2) + 1 + uli = imo*nocc(2) + ipole = (imo-1)*nvir(2) + amo + npoles1 + temp(ipole) = ga_ddot_patch(g_tmp2(2),'n',1,nri,1,nocc(2), + $ pars%g_erioo(2),'n',1,nri,lli,uli) + enddo + enddo + endif + if (pars%me.eq.0) + $ call ga_acc(kv,1,npoles,itrial,itrial,temp,npoles,-1d0) + enddo + !call ga_copy(kv,mv) + + if(.not. (ga_destroy(g_tmp1(1)) .and. ga_destroy(g_tmp2(1)))) + $ call errquit('bse_davidson_kvec: ga_destroy_failed',93,GA_ERR) + if(pars%ipol.gt.1) then + if(.not. (ga_destroy(g_tmp1(2)) .and. ga_destroy(g_tmp2(2)))) + $ call errquit('bse_davidson_kvec: ga_destroy_failed',96,GA_ERR) + endif + + if (.not.ga_create(mt_dbl,nri,ntrials,'g_tmp1',0,ntrials, + $ g_tmp1)) + $ call errquit('bse_davidson_kvec: ga create failed',102,GA_ERR) + if (.not.ga_create(mt_dbl,npoles,ntrials,'g_tmp2',0,ntrials, + $ g_tmp2)) + $ call errquit('bse_davidson_kvec: ga create failed',105,GA_ERR) + + ! Hartree + if (dohartree) then + factor = 4d0/pars%ipol + if (pars%tda) factor = factor/2d0 + call ga_matmul_patch('n','n',1d0,0d0, + $ pars%g_eriov(1),1,nri,1,npoles1, + $ vec,1,npoles1,1,ntrials, + $ g_tmp1,1,nri,1,ntrials) + if (pars%ipol.gt.1) then + call ga_matmul_patch('n','n',1d0,1d0, + $ pars%g_eriov(2),1,nri,1,npoles2, + $ vec,1+npoles1,npoles,1,ntrials, + $ g_tmp1,1,nri,1,ntrials) + endif + call ga_matmul_patch('t','n',factor,1d0, + $ pars%g_eriov(1),1,npoles1,1,nri, + $ g_tmp1,1,nri,1,ntrials, + $ kv,1,npoles1,1,ntrials) + if (pars%ipol.gt.1) then + call ga_matmul_patch('t','n',factor,1d0, + $ pars%g_eriov(2),1,npoles2,1,nri, + $ g_tmp1,1,nri,1,ntrials, + $ kv,1+npoles1,npoles,1,ntrials) + endif + endif + + ! Screened Coulomb B + call ga_matmul_patch('n','n',1d0,0d0, + $ pars%g_wov(1),1,nri,1,npoles1, + $ vec,1,npoles1,1,ntrials, + $ g_tmp1,1,nri,1,ntrials) + call ga_matmul_patch('t','n',-1d0,0d0, + $ pars%g_wov(1),1,npoles1,1,nri, + $ g_tmp1,1,nri,1,ntrials, + $ g_tmp2,1,npoles1,1,ntrials) + if (pars%ipol.gt.1) then + call ga_matmul_patch('n','n',1d0,0d0, + $ pars%g_wov(2),1,nri,1,npoles2, + $ vec,1+npoles1,npoles,1,ntrials, + $ g_tmp1,1,nri,1,ntrials) + call ga_matmul_patch('t','n',-1d0,0d0, + $ pars%g_wov(2),1,npoles2,1,nri, + $ g_tmp1,1,nri,1,ntrials, + $ g_tmp2,1+npoles1,npoles,1,ntrials) + endif + + call ga_add_patch(1d0,g_tmp2,1,npoles,1,ntrials, + $ 1d0,kv,1,npoles,1,ntrials, + $ kv,1,npoles,1,ntrials) + !call ga_add(-1d0,g_tmp4,1d0,mv,mv) + + if(.not.(ga_destroy(g_tmp1) .and. ga_destroy(g_tmp2))) + $ call errquit('bse_davidson_kvec: ga_destroy failed',156,GA_ERR) + + end diff --git a/src/bsemol/bse_davidson_mvec.F b/src/bsemol/bse_davidson_mvec.F new file mode 100644 index 0000000000..5d50fbfb6e --- /dev/null +++ b/src/bsemol/bse_davidson_mvec.F @@ -0,0 +1,131 @@ + + subroutine bse_davidson_mvec(pars,vec,wia,mv,npoles,ntrials) + implicit none +#include "bse.fh" +#include "global.fh" +#include "mafdecls.fh" +#include "errquit.fh" + + type(bse_params_t) :: pars + + integer npoles,ntrials + integer vec,wia,kv,mv + + integer ipol,npoles1,npoles2,uli,lli,ipole + integer nocc(2),nvir(2),amo,imo,itrial,ula,lla,nri + + integer g_tmp1(2),g_tmp2(2) + double precision temp(npoles),factor + + ipol = pars%ipol + nri = pars%nri + nocc = pars%nocc + nvir = pars%nvir + npoles1 = nocc(1)*nvir(1) + npoles2 = nocc(2)*nvir(2) + + ! Create temporary arrays + if (.not.ga_create(mt_dbl,nvir(1),nocc(1),'g_tmp1',0,nocc(1), + $ g_tmp1)) + $ call errquit('bse_davidson_kvec: ga create failed',34,GA_ERR) + if (.not.ga_create(mt_dbl,nri,nocc(1),'g_tmp2',0,nocc(1), + $ g_tmp2)) + $ call errquit('bse_davidson_kvec: ga create failed',37,GA_ERR) + + if (pars%ipol.gt.1) then + if (.not.ga_create(mt_dbl,nvir(2),nocc(2),'g_tmp1',0,nocc(2), + $ g_tmp1(2))) + $ call errquit('bse_davidson_kvec: ga create failed',48,GA_ERR) + if (.not.ga_create(mt_dbl,nri,nocc(2),'g_tmp2',0,nocc(2), + $ g_tmp2(2))) + $ call errquit('bse_davidson_kvec: ga create failed',51,GA_ERR) + endif + + + call ga_copy(vec,mv) + call ga_scale_rows(mv,wia) + + do itrial=1,ntrials + call ga_copy_patch('n',vec,1,npoles1,itrial,itrial, + $ g_tmp1,1,nvir(1),1,nocc(1)) + ! Screened Coulomb A + do amo=1,nvir(1) + lla = (amo-1)*nvir(1) + 1 + ula = amo*nvir(1) + call ga_matmul_patch('n','n',1d0,0d0, + $ pars%g_erivv(1),1,nri,lla,ula, + $ g_tmp1,1,nvir(1),1,nocc(1), + $ g_tmp2,1,nri,1,nocc(1)) + do imo=1,nocc(1) + lli = (imo-1)*nocc(1) + 1 + uli = imo*nocc(1) + ipole = (imo-1)*nvir(1) + amo + temp(ipole) = ga_ddot_patch(g_tmp2,'n',1,nri,1,nocc(1), + $ pars%g_erioo(1),'n',1,nri,lli,uli) + enddo + enddo + if (ipol.gt.1) then + call ga_copy_patch('n',vec,1+npoles1,npoles,itrial,itrial, + $ g_tmp1,1,nvir(2),1,nocc(2)) + do amo=1,nvir(2) + lla = (amo-1)*nvir(2) + 1 + ula = amo*nvir(2) + call ga_matmul_patch('n','n',1d0,0d0, + $ pars%g_erivv(2),1,nri,lla,ula, + $ g_tmp1,1,nvir(2),1,nocc(2), + $ g_tmp2,1,nri,1,nocc(2)) + do imo=1,nocc(2) + lli = (imo-1)*nocc(2) + 1 + uli = imo*nocc(2) + ipole = (imo-1)*nvir(2) + amo + npoles1 + temp(ipole) = ga_ddot_patch(g_tmp2,'n',1,nri,1,nocc(2), + $ pars%g_erioo(2),'n',1,nri,lli,uli) + enddo + enddo + endif + if (pars%me.eq.0) + $ call ga_acc(mv,1,npoles,itrial,itrial,temp,npoles,-1d0) + enddo + + if(.not. (ga_destroy(g_tmp1(1)) .and. ga_destroy(g_tmp2(1)))) + $ call errquit('bse_davidson_kvec: ga_destroy_failed',93,GA_ERR) + if(pars%ipol.gt.1) then + if(.not. (ga_destroy(g_tmp1(2)) .and. ga_destroy(g_tmp2(2)))) + $ call errquit('bse_davidson_kvec: ga_destroy_failed',96,GA_ERR) + endif + + if (.not.ga_create(mt_dbl,nri,ntrials,'g_tmp1',0,ntrials, + $ g_tmp1)) + $ call errquit('bse_davidson_kvec: ga create failed',102,GA_ERR) + if (.not.ga_create(mt_dbl,npoles,ntrials,'g_tmp2',0,ntrials, + $ g_tmp2)) + $ call errquit('bse_davidson_kvec: ga create failed',105,GA_ERR) + + ! Screened Coulomb B + call ga_matmul_patch('n','n',1d0,0d0, + $ pars%g_wov(1),1,nri,1,npoles1, + $ vec,1,npoles1,1,ntrials, + $ g_tmp1,1,nri,1,ntrials) + call ga_matmul_patch('t','n',-1d0,0d0, + $ pars%g_wov(1),1,npoles1,1,nri, + $ g_tmp1,1,nri,1,ntrials, + $ g_tmp2,1,npoles1,1,ntrials) + if (pars%ipol.gt.1) then + call ga_matmul_patch('n','n',1d0,0d0, + $ pars%g_wov(2),1,nri,1,npoles2, + $ vec,1+npoles1,npoles,1,ntrials, + $ g_tmp1,1,nri,1,ntrials) + call ga_matmul_patch('t','n',-1d0,0d0, + $ pars%g_wov(2),1,npoles2,1,nri, + $ g_tmp1,1,nri,1,ntrials, + $ g_tmp2,1+npoles1,npoles,1,ntrials) + endif + call ga_add_patch(-1d0,g_tmp2,1,npoles,1,ntrials, + $ 1d0,mv,1,npoles,1,ntrials, + $ mv,1,npoles,1,ntrials) + + if(.not.(ga_destroy(g_tmp1) .and. ga_destroy(g_tmp2))) + $ call errquit('bse_davidson_kvec: ga_destroy failed',156,GA_ERR) + + end + diff --git a/src/bsemol/bse_davidson_newvecs.F b/src/bsemol/bse_davidson_newvecs.F new file mode 100644 index 0000000000..f82591bfe2 --- /dev/null +++ b/src/bsemol/bse_davidson_newvecs.F @@ -0,0 +1,73 @@ + subroutine bse_davidson_newvecs(ntrials,ndim,npoles,g_x,g_y, + $ g_ym,g_s,g_ks,g_mks,g_w,g_kw,g_mkw,lvkv,vkv) + + implicit none +#include "global.fh" +#include "mafdecls.fh" +#include "errquit.fh" + + integer ntrials,ndim,npoles + integer g_x,g_y,g_ym,g_s,g_ks,g_mks + integer g_w,g_kw,g_mkw,lvkv,vkv + + logical stat + + ! Save working vectors + call ga_copy(g_s,g_x) + call ga_copy(g_ks,g_y) + call ga_copy(g_mks,g_ym) + + ! Destroy working vectors + stat = ga_destroy(g_s).and.ga_destroy(g_ks).and.ga_destroy(g_mks) + if(.not.stat) + $ call errquit('bse_davidson_newvecs: destroy failed',1,GA_ERR) + + ! Create new working vectors + stat = ga_create(mt_dbl,npoles,ntrials+ndim,'s',0, + $ ntrials+ndim,g_s) .and. + $ ga_create(mt_dbl,npoles,ntrials+ndim,'ks',0, + $ ntrials+ndim,g_ks) .and. + $ ga_create(mt_dbl,npoles,ntrials+ndim,'mks',0, + $ ntrials+ndim,g_mks) + if(.not.stat) + $ call errquit('bse_davidson_newvecs: create failed',1,GA_ERR) + + + ! Copy back old working vectors plus new vectors + call ga_copy_patch('n',g_x,1,npoles,1,ntrials, + $ g_s,1,npoles,1,ntrials) + call ga_copy_patch('n',g_w,1,npoles,1,ndim, + $ g_s,1,npoles,ntrials+1,ntrials+ndim) + + call ga_copy_patch('n',g_y,1,npoles,1,ntrials, + $ g_ks,1,npoles,1,ntrials) + call ga_copy_patch('n',g_kw,1,npoles,1,ndim, + $ g_ks,1,npoles,ntrials+1,ntrials+ndim) + + call ga_copy_patch('n',g_ym,1,npoles,1,ntrials, + $ g_mks,1,npoles,1,ntrials) + call ga_copy_patch('n',g_mkw,1,npoles,1,ndim, + $ g_mks,1,npoles,ntrials+1,ntrials+ndim) + + ! Update ntrials + ntrials = ntrials + ndim + + ! Destroy old solutions + stat = ga_destroy(g_x).and.ga_destroy(g_y).and.ga_destroy(g_ym) + if(.not.stat) + $ call errquit('bse_davidson_newvecs: destroy failed',2,GA_ERR) + + ! Create new solutions + stat = ga_create(mt_dbl,npoles,ntrials,'x',0,ntrials,g_x) .and. + $ ga_create(mt_dbl,npoles,ntrials,'y',0,ntrials,g_y) .and. + $ ga_create(mt_dbl,npoles,ntrials,'ym',0,ntrials,g_ym) + if(.not.stat) + $ call errquit('bse_davidson_newvecs: create failed',2,GA_ERR) + + ! Reallocate VKV + if(.not.ma_chop_stack(lvkv)) + $ call errquit('bse_davidson_restart: alloc failed',3,MA_ERR) + if(.not.ma_push_get(mt_dbl,ntrials*ntrials,'vkv',lvkv,vkv)) + $ call errquit('bse_davidson_restart: dealloc failed',3,MA_ERR) + + end subroutine diff --git a/src/bsemol/bse_davidson_residual.F b/src/bsemol/bse_davidson_residual.F new file mode 100644 index 0000000000..a566fee17d --- /dev/null +++ b/src/bsemol/bse_davidson_residual.F @@ -0,0 +1,29 @@ + subroutine bse_davidson_residual(nroots,ntrials,omega,x,ym,r,dia, + $ w,npoles,ndim,norm) + implicit none + + logical main + integer nroots,npoles,ntrials,ndim + double precision omega(nroots) + double precision ym(npoles,nroots),r(npoles) + double precision x(npoles,nroots),dia(npoles),w(npoles,nroots) + + integer iroot + double precision norm(nroots) + + ndim = 0 + do iroot=1,nroots + r(:) = ym(:,iroot) - omega(iroot)*x(:,iroot) + norm(iroot) = maxval(abs(r)) + call ga_dgop((/1/),norm(iroot),1,"maxabs") + if (norm(iroot).gt.1d-5) then + ndim = ndim + 1 + w(:,ndim) = r(:)/(dia(:) - omega(iroot)) + endif + enddo + + end subroutine + + + + diff --git a/src/bsemol/bse_davidson_restart.F b/src/bsemol/bse_davidson_restart.F new file mode 100644 index 0000000000..f2ba4f9e9e --- /dev/null +++ b/src/bsemol/bse_davidson_restart.F @@ -0,0 +1,68 @@ + subroutine bse_davidson_restart(ntrials,nroots,npoles,g_x,g_y, + $ g_ym,g_s,g_ks,g_mks,lvkv,vkv) + + implicit none +#include "global.fh" +#include "mafdecls.fh" +#include "errquit.fh" + + integer ntrials,nroots,npoles + integer g_x,g_y,g_ym,g_s,g_ks,g_mks,lvkv,vkv + + logical stat + + ! Reset the number of trials + ntrials = nroots + + ! Copy current solutions + call ga_copy_patch('n',g_x,1,npoles,1,ntrials, + $ g_s,1,npoles,1,ntrials) + call ga_copy_patch('n',g_y,1,npoles,1,ntrials, + $ g_ks,1,npoles,1,ntrials) + call ga_copy_patch('n',g_ym,1,npoles,1,ntrials, + $ g_mks,1,npoles,1,ntrials) + + ! Destroy current solutions + stat = ga_destroy(g_x).and.ga_destroy(g_y).and.ga_destroy(g_ym) + if (.not.stat) + $ call errquit('bse_davidson_restar: failed destroy',1,GA_ERR) + + ! Re-create new solutions + stat = ga_create(mt_dbl,npoles,ntrials,'x',0,ntrials,g_x) .and. + $ ga_create(mt_dbl,npoles,ntrials,'y',0,ntrials,g_y) .and. + $ ga_create(mt_dbl,npoles,ntrials,'ym',0,ntrials,g_ym) + if (.not.stat) + $ call errquit('bse_davidson_restar: failed create',1,GA_ERR) + + ! Copy back solutions + call ga_copy_patch('n',g_s,1,npoles,1,ntrials, + $ g_x,1,npoles,1,ntrials) + call ga_copy_patch('n',g_ks,1,npoles,1,ntrials, + $ g_y,1,npoles,1,ntrials) + call ga_copy_patch('n',g_mks,1,npoles,1,ntrials, + $ g_ym,1,npoles,1,ntrials) + + ! Destroy working vectors + stat = ga_destroy(g_s).and.ga_destroy(g_ks).and.ga_destroy(g_mks) + if (.not.stat) + $ call errquit('bse_davidson_restar: failed destroy',2,GA_ERR) + + ! Re-create new solutions + stat = ga_create(mt_dbl,npoles,ntrials,'s',0,ntrials,g_s) .and. + $ ga_create(mt_dbl,npoles,ntrials,'ks',0,ntrials,g_ks) .and. + $ ga_create(mt_dbl,npoles,ntrials,'mks',0,ntrials,g_mks) + if (.not.stat) + $ call errquit('bse_davidson_restar: failed create',2,GA_ERR) + + ! Copy current solutions + call ga_copy(g_x,g_s) + call ga_copy(g_y,g_ks) + call ga_copy(g_ym,g_mks) + + ! Reallocate VKV + if(.not.ma_chop_stack(lvkv)) + $ call errquit('bse_davidson_restart: alloc failed',3,MA_ERR) + if(.not.ma_push_get(mt_dbl,ntrials*ntrials,'vkv',lvkv,vkv)) + $ call errquit('bse_davidson_restart: dealloc failed',3,MA_ERR) + + end subroutine diff --git a/src/bsemol/bse_defaults.F b/src/bsemol/bse_defaults.F index 3e253d9966..167d442f61 100644 --- a/src/bsemol/bse_defaults.F +++ b/src/bsemol/bse_defaults.F @@ -21,6 +21,7 @@ & 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:nmax',mt_int,1,100) .and. & rtdb_put(rtdb,'bse:init',mt_log,1,.true.) if (.not.ok) diff --git a/src/bsemol/bse_init.F b/src/bsemol/bse_init.F index 8c028deadd..7ffc92b036 100644 --- a/src/bsemol/bse_init.F +++ b/src/bsemol/bse_init.F @@ -111,6 +111,14 @@ c if(.not.rtdb_get(pars%rtdb,'bse:nroots',mt_int,1,pars%nroots)) & call errquit(pname//'could not read from rtdb',5,RTDB_ERR) + ! Maximum Number of roots + if(.not.rtdb_get(pars%rtdb,'bse:nmax',mt_int,1,pars%nmax)) + & call errquit(pname//'could not read from rtdb',5,RTDB_ERR) + if (pars%nmax.le.2*pars%nroots) then + call errquit(pname//'maxroots should be at least 2*nroots', + $ 119,INPUT_ERR) + endif + ! Get some info from DFT module call dft_rdinput(pars%rtdb) if (cdfit) then @@ -212,6 +220,10 @@ c pars%nroots = sum(pars%npoles) endif + if (pars%nmax.gt.sum(pars%npoles)) then + pars%nmax = sum(pars%npoles) + endif + ! Obtain NCAP DD shift 1001 format(2X,A12,F8.3,A3) diff --git a/src/bsemol/bse_input.F b/src/bsemol/bse_input.F index 2a04d8de7d..dfad98ab3f 100644 --- a/src/bsemol/bse_input.F +++ b/src/bsemol/bse_input.F @@ -18,10 +18,10 @@ c character*255 test character(*), parameter :: pname = "bse_input: " c - integer, parameter :: num_dirs = 9 + integer, parameter :: num_dirs = 10 character*11 dirs(num_dirs) data dirs / 'method','tda','ncap','ncapr','singlet','triplet', - & 'outwin', 'nroots', 'end' / + & 'outwin', 'nroots', 'maxroots', 'end' / c integer ind, itemp double precision threshold, eta, temp @@ -45,7 +45,7 @@ c call errquit(pname//'unknown directive', 0, INPUT_ERR) end if c - goto (100, 200, 300, 301, 400, 500, 600, 700, 9999) ind + goto (100, 200, 300, 301, 400, 500, 600, 700, 800, 9999) ind call errquit(pname//'unknown directive',ind, INPUT_ERR) c @@ -59,12 +59,12 @@ c if (inp_compare(.false.,'analytic',test)) then if (.not.rtdb_put(rtdb, 'bse:analytic', mt_log, 1, .true.)) & call errquit(pname//'rtdb_put failed', 101, RTDB_ERR) - elseif (inp_compare(.false.,'bse:davidson',test)) then + elseif (inp_compare(.false.,'davidson',test)) then if (.not.rtdb_put(rtdb, 'bse:davidson', mt_log, 1, .true.)) & call errquit(pname//'rtdb_put failed', 102, RTDB_ERR) if (.not.rtdb_put(rtdb, 'bse:analytic', mt_log, 1, .false.)) & call errquit(pname//'rtdb_put failed', 101, RTDB_ERR) - elseif (inp_compare(.false.,'bse:lanczos',test)) then + elseif (inp_compare(.false.,'lanczos',test)) then if (.not.rtdb_put(rtdb, 'bse:lanczos', mt_log, 1, .true.)) & call errquit(pname//'rtdb_put failed', 103, RTDB_ERR) if (.not.rtdb_put(rtdb, 'bse:analytic', mt_log, 1, .false.)) @@ -129,7 +129,14 @@ c if(.not.rtdb_put(rtdb,'bse:nroots',mt_int,1,itemp)) & call errquit(pname//'rtdb_put failed', 700, RTDB_ERR) goto 10 - +c +c maximum number of trial vectors +c + 800 if(.not.inp_i(itemp)) + & call errquit(pname//'could not read maxroots',800,INPUT_ERR) + if(.not.rtdb_put(rtdb,'bse:nmax',mt_int,1,itemp)) + & call errquit(pname//'rtdb_put failed', 800, RTDB_ERR) + goto 10 c c normal termination c