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

TDA, printing

Browse source
This commit is contained in:
Daniel Mejia-Rodriguez 2022-08-13 21:13:09 -07:00 • committed by edoapra
parent 1cbb583c2d
commit ca338cc391
No known key found for this signature in database
GPG key ID: 9E6A0B70826967BA
12 changed files with 423 additions and 39 deletions
Download patch file Download diff file Expand all files Collapse all files

4
src/bsemol/GNUmakefile
View file

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

9
src/bsemol/bse.fh
View file

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

12
src/bsemol/bse_analytic.F
View file

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

12
src/bsemol/bse_casida.F → src/bsemol/bse_analytic_cis.F
View file

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

119
src/bsemol/bse_analytic_tda.F Normal file
View file

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

6
src/bsemol/bse_defaults.F
View file

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

85
src/bsemol/bse_dipole.F Normal file
View file

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

15
src/bsemol/bse_finalize.F
View file

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

38
src/bsemol/bse_init.F
View file

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

50
src/bsemol/bse_input.F
View file

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

87
src/bsemol/bse_oscstr.F Normal file
View file

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

25
src/bsemol/bse_wmn.F
View file

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