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

BSE K-davidson

Browse source
This commit is contained in:
Daniel Mejia-Rodriguez 2023-07-24 10:54:00 -07:00 • committed by edoapra
parent 00583977fa
commit 396fd02972
No known key found for this signature in database
GPG key ID: 9E6A0B70826967BA
16 changed files with 1063 additions and 17 deletions
Download patch file Download diff file Expand all files Collapse all files

6
src/bsemol/GNUmakefile
View file

@ -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

4
src/bsemol/bse.fh
View file

@ -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

6
src/bsemol/bse_analytic.F
View file

@ -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

6
src/bsemol/bse_analytic_cis.F
View file

@ -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)

302
src/bsemol/bse_davidson.F
View file

@ -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

98
src/bsemol/bse_davidson_diagonal.F Normal file
View file

@ -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

96
src/bsemol/bse_davidson_guess.F Normal file
View file

@ -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

71
src/bsemol/bse_davidson_kortho.F Normal file
View file

@ -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

158
src/bsemol/bse_davidson_kvec.F Normal file
View file

@ -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

131
src/bsemol/bse_davidson_mvec.F Normal file
View file

@ -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

73
src/bsemol/bse_davidson_newvecs.F Normal file
View file

@ -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

29
src/bsemol/bse_davidson_residual.F Normal file
View file

@ -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

68
src/bsemol/bse_davidson_restart.F Normal file
View file

@ -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

1
src/bsemol/bse_defaults.F
View file

@ -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)

12
src/bsemol/bse_init.F
View file

@ -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)

19
src/bsemol/bse_input.F
View file

@ -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