mirror of
https://github.com/nwchemgit/nwchem.git
synced 2026-07-28 14:15:30 -04:00
CS singlet
This commit is contained in:
Daniel Mejia-Rodriguez
edoapra
parent
3c22431c69
commit
1cbb583c2d
10 changed files with 768 additions and 11 deletions
|
|
@ -1,9 +1,10 @@
|
|||
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_finalize.o bse_ri_init.o bse_wmn.o bse_casida.o
|
||||
|
||||
USES_BLAS =
|
||||
OBJ_OPTIMIZE =
|
||||
USES_BLAS = bse_buildw.F bse_ri_init.F bse_wmn.F
|
||||
|
||||
OBJ_OPTIMIZE = bse_buildw.o
|
||||
|
||||
LIBRARY = libbsemol.a
|
||||
|
||||
|
|
|
|||
|
|
@ -9,7 +9,7 @@
|
|||
logical :: main, odft
|
||||
c
|
||||
integer :: natoms
|
||||
integer :: nmo
|
||||
integer :: nmo, nri
|
||||
integer :: nocc(2),nvir(2),npoles(2)
|
||||
integer :: noqp(2),nvqp(2),nstates(2),mynpoles(2)
|
||||
integer :: oolo(2),oohi(2)
|
||||
|
|
@ -27,6 +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
|
||||
logical :: noio
|
||||
c
|
||||
logical :: ncap, davidson, lanczos, analytic, tda
|
||||
|
|
|
|||
|
|
@ -2,8 +2,76 @@
|
|||
|
||||
implicit none
|
||||
|
||||
#include "global.fh"
|
||||
#include "errquit.fh"
|
||||
#include "mafdecls.fh"
|
||||
#include "bse.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
|
||||
|
||||
! 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
|
||||
|
||||
! allocation
|
||||
if(.not.ma_push_get(mt_dbl,totnpoles,'omega',l_omega,k_omega))
|
||||
& call errquit(pname//'failed to allocate omega',0,MA_ERR)
|
||||
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,nri**2,'W',lW,kW))
|
||||
& call errquit(pname//'failed to allocate W',0,MA_ERR)
|
||||
|
||||
! obtain "old" eigenvalue difference to build W
|
||||
do isp=1,pars%ipol
|
||||
lSize = (isp-1)*pars%npoles(1)
|
||||
call gw_get_eia(dbl_mb(pars%k_mf_evals+(isp-1)*pars%nmo),
|
||||
& dbl_mb(k_dia+lSize),pars%nocc(isp),
|
||||
& pars%nvir(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)*pars%npoles(1)
|
||||
call gw_get_eia(dbl_mb(pars%k_gw_evals+(isp-1)*pars%nmo),
|
||||
& dbl_mb(k_dia+lSize),pars%nocc(isp),
|
||||
& pars%nvir(isp))
|
||||
enddo
|
||||
|
||||
! Build and diagonalize Casida-like matrix
|
||||
call bse_casida(pars,dbl_mb(k_dia),dbl_mb(k_omega),maxpoles,
|
||||
& totnpoles)
|
||||
|
||||
! Compute oscillator strengths
|
||||
|
||||
! Print output
|
||||
|
||||
! Destroy GAs
|
||||
if (.not.ga_destroy(pars%g_apb))
|
||||
& call errquit(pname//'could not destroy temp GA',1,GA_ERR)
|
||||
|
||||
! chop stack
|
||||
if(.not.ma_chop_stack(l_omega))
|
||||
$ call errquit(pname//'failed to chop stack',0,MA_ERR)
|
||||
|
||||
end subroutine bse_analytic
|
||||
|
|
|
|||
75
src/bsemol/bse_buildw.F
Normal file
75
src/bsemol/bse_buildw.F
Normal file
|
|
@ -0,0 +1,75 @@
|
|||
subroutine bse_buildW(pars,wia,w,nmo,nri,maxpoles,ipol)
|
||||
implicit none
|
||||
#include "errquit.fh"
|
||||
#include "mafdecls.fh"
|
||||
#include "global.fh"
|
||||
#include "stdio.fh"
|
||||
#include "bse.fh"
|
||||
|
||||
type(bse_params_t) :: pars
|
||||
integer nmo, nri, maxpoles, ipol
|
||||
|
||||
character(*), parameter :: pname = 'bse_buildw: '
|
||||
|
||||
double precision wia(maxpoles,ipol)
|
||||
double precision w(nri,nri)
|
||||
|
||||
integer :: ovlo,ovhi,kov(2),isp
|
||||
integer :: ldpi, sizepi
|
||||
integer :: alloc, info, ld, ipole
|
||||
double precision,allocatable :: factor(:,:)
|
||||
|
||||
allocate(factor(maxpoles,ipol),stat=alloc)
|
||||
if (alloc.ne.0)
|
||||
& call errquit(pname//'allocation failed',0,MA_ERR)
|
||||
C
|
||||
C Decide the leading dimension of the \Pi array according to
|
||||
C the rectangular full-packed (RFP) format
|
||||
C
|
||||
if (mod(nri,2).eq.0) then
|
||||
ldpi = nri + 1
|
||||
else
|
||||
ldpi = nri
|
||||
endif
|
||||
sizepi = (nri*(nri+1))/2
|
||||
|
||||
! Get pointer to local integrals
|
||||
do isp=1,pars%ipol
|
||||
kov(isp) = 1
|
||||
if (pars%mynpoles(isp).lt.1) cycle
|
||||
|
||||
ovlo = pars%ovlo(isp)
|
||||
ovhi = pars%ovhi(isp)
|
||||
call ga_access(pars%g_eriov(1),1,nri,ovlo,ovhi,kov(isp),ld)
|
||||
|
||||
!$omp parallel do
|
||||
do ipole=1,pars%mynpoles(isp)
|
||||
factor(ipole,isp) = dsqrt(1.0d0/wia(ipole,isp))
|
||||
enddo
|
||||
!$omp end parallel do
|
||||
enddo
|
||||
|
||||
! Build the Polarizability matrix and transform it to the
|
||||
! dielectric matrix \epsilon = 1 - \pi
|
||||
call gw_cdgw_buildpi('w',dbl_mb(kov(1)),dbl_mb(kov(2)),w,factor,
|
||||
$ ldpi,nri,pars%mynpoles,pars%me.eq.0,
|
||||
$ ipol,maxpoles)
|
||||
call ga_dgop(1038,w,sizepi,'+')
|
||||
|
||||
!Factorize dieletric matrix for further use
|
||||
call dpftrf('n','l',nri,w,info)
|
||||
|
||||
!Deallocate
|
||||
deallocate(factor,stat=alloc)
|
||||
if (alloc.ne.0)
|
||||
& call errquit(pname//'deallocation failed',0,MA_ERR)
|
||||
|
||||
!Release GAs
|
||||
do isp=1,ipol
|
||||
if (pars%mynpoles(isp).lt.1) cycle
|
||||
ovlo = pars%ovlo(isp)
|
||||
ovhi = pars%ovhi(isp)
|
||||
call ga_release_update(pars%g_eriov(isp),1,nri,ovlo,ovhi)
|
||||
enddo
|
||||
|
||||
end
|
||||
191
src/bsemol/bse_casida.F
Normal file
191
src/bsemol/bse_casida.F
Normal file
|
|
@ -0,0 +1,191 @@
|
|||
subroutine bse_casida(pars,wia,omega,maxpoles,totpoles)
|
||||
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)
|
||||
|
||||
character(*),parameter :: pname = 'bse_casida: '
|
||||
|
||||
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
|
||||
|
||||
#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.ga_create(mt_dbl,totpoles,totpoles,'amb',0,totpoles,
|
||||
& g_amb))
|
||||
& 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 from A
|
||||
if(.not.ga_create(mt_dbl,totpoles,totpoles,'tmp',0,totpoles,
|
||||
& g_tmp))
|
||||
& call errquit(pname//'could not create temp GA',1,GA_ERR)
|
||||
call ga_dgemm('t','n',totpoles,totpoles,pars%nri,-1d0,
|
||||
& pars%g_wov(1),pars%g_wov(1),0d0,g_tmp)
|
||||
call ga_copy(g_tmp,pars%g_apb)
|
||||
call ga_zero(g_amb)
|
||||
call ga_add(-1d0,g_tmp,1d0,g_amb,g_amb)
|
||||
if (.not.ga_destroy(g_tmp))
|
||||
& call errquit(pname//'could not destroy temp GA',1,GA_ERR)
|
||||
|
||||
! Screened Coulomb contribution from B
|
||||
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_acc(pars%g_apb,kglobal,kglobal,lglobal,lglobal,
|
||||
& dbl_mb(ktmp+(j-jlo)*nocc2(1)+i-1),1,1d0)
|
||||
call ga_acc(g_amb,kglobal,kglobal,lglobal,lglobal,
|
||||
& dbl_mb(ktmp+(j-jlo)*nocc2(1)+i-1),1,1d0)
|
||||
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
|
||||
call ga_dgemm('t','n',totpoles,totpoles,pars%nri,4d0,
|
||||
& pars%g_eriov(1),pars%g_eriov(1),1d0,pars%g_apb)
|
||||
|
||||
! Diagonal term
|
||||
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_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)
|
||||
if (.not.ga_destroy(pars%g_wov(1)))
|
||||
& call errquit(pname//'could not destroy wov GA',1,GA_ERR)
|
||||
|
||||
|
||||
! First step is to get (A-B)^(1/2)
|
||||
call ga_sync()
|
||||
#if defined(PARALLEL_DIAG)
|
||||
#ifdef SCALAPACK
|
||||
info= ga_cholesky('L',g_amb)
|
||||
#else
|
||||
call ga_chol(g_amb, g_amb, info)
|
||||
#endif
|
||||
#else
|
||||
call ga_chol_seq(g_amb, g_amb, info)
|
||||
#endif
|
||||
|
||||
! Transform (A+B) matrix with Cholesky factors
|
||||
if(.not.ga_create(mt_dbl,totpoles,totpoles,'tmp',0,totpoles,
|
||||
& g_tmp))
|
||||
& call errquit(pname//'could not create temp GA',1,GA_ERR)
|
||||
call ga_dgemm('n','n',totpoles,totpoles,totpoles,1d0,
|
||||
& pars%g_apb,g_amb,0d0,g_tmp)
|
||||
call ga_dgemm('t','n',totpoles,totpoles,totpoles,1d0,
|
||||
& g_amb,g_tmp,0d0,pars%g_apb)
|
||||
if (.not.ga_destroy(g_tmp))
|
||||
& call errquit(pname//'could not destroy temp 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
|
||||
|
||||
!
|
||||
! Since we transformed the eigenvalue equation to an Hermitian
|
||||
! form, we obtain the square of the neutral excitations
|
||||
! (Equation (32)).
|
||||
!
|
||||
!$omp parallel do
|
||||
do ipole=1,totpoles
|
||||
omega(ipole) = dsqrt(omega(ipole))
|
||||
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))
|
||||
& call errquit(pname//'could not create temp GA',1,GA_ERR)
|
||||
call ga_dgemm('n','n',totpoles,totpoles,totpoles,1d0,
|
||||
& g_amb,pars%g_apb,0d0,g_tmp)
|
||||
|
||||
! Gather neutral excitations vectors
|
||||
call nga_distribution(g_diag,pars%me,ilo,ihi)
|
||||
if(ihi.ge.ilo) then
|
||||
call nga_put(g_diag,ilo,ihi,1d0/dsqrt(omega(ilo:ihi)),
|
||||
& ihi-ilo+1)
|
||||
endif
|
||||
call ga_sync()
|
||||
call ga_scale_cols(g_tmp,g_diag)
|
||||
|
||||
! Destroy unused GAs
|
||||
call ga_copy(g_tmp,pars%g_apb)
|
||||
if (.not.ga_destroy(g_tmp))
|
||||
& call errquit(pname//'could not destroy temp GA',1,GA_ERR)
|
||||
if (.not.ga_destroy(g_amb))
|
||||
& call errquit(pname//'could not destroy temp GA',1,GA_ERR)
|
||||
if (.not.ga_destroy(g_diag))
|
||||
& call errquit(pname//'could not destroy temp GA',1,GA_ERR)
|
||||
|
||||
else
|
||||
|
||||
endif
|
||||
|
||||
end subroutine
|
||||
|
|
@ -21,7 +21,7 @@
|
|||
type(bse_params_t) :: pars
|
||||
character(*), parameter :: pname = 'bse_init: '
|
||||
|
||||
character(len=255) :: basisname, scftype
|
||||
character(len=255) :: basisname, scftype, gwevals
|
||||
integer i, ilo, ihi, j
|
||||
integer nbf_temp, ipol_temp, nocc, nvir, nmo_temp(2)
|
||||
logical ldum
|
||||
|
|
@ -39,6 +39,8 @@ c
|
|||
character(len=80) names
|
||||
double precision t0, ehomo, vxddp, vxddm, aq2, zeta
|
||||
|
||||
integer,parameter :: unitno = 66
|
||||
|
||||
#ifdef USE_OPENMP
|
||||
integer,external :: omp_get_max_threads, omp_get_thread_num
|
||||
integer,external :: omp_get_num_threads
|
||||
|
|
@ -119,6 +121,7 @@ c
|
|||
if (bas_is_spherical(ao_bas_han).and.
|
||||
$ (.not.bas_is_spherical(cd_bas_han)))
|
||||
$ call int_app_set_no_texas(pars%rtdb)
|
||||
pars%nri = nbf_cd
|
||||
|
||||
!Spin multiplicity, occupations
|
||||
if (.not.rtdb_get(pars%rtdb,'dft:ipol',mt_int,1,ipol))
|
||||
|
|
@ -193,9 +196,33 @@ c
|
|||
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)
|
||||
dbl_mb(pars%k_gw_evals+(i-1)*nbf_ao:
|
||||
& pars%k_gw_evals+(i-1)*nbf_ao+pars%nocc(i)-1) =
|
||||
& dbl_mb(pars%k_gw_evals+(i-1)*nbf_ao:
|
||||
& pars%k_gw_evals+(i-1)*nbf_ao+pars%nocc(i)-1) +
|
||||
& vxddp
|
||||
dbl_mb(pars%k_gw_evals+(i-1)*nbf_ao+pars%nocc(i):
|
||||
& pars%k_gw_evals+(i-1)*nbf_ao+nbf_ao-1) =
|
||||
& dbl_mb(pars%k_gw_evals+(i-1)*nbf_ao+pars%nocc(i):
|
||||
& pars%k_gw_evals+(i-1)*nbf_ao+nbf_ao-1) +
|
||||
& vxddm
|
||||
enddo
|
||||
endif
|
||||
|
||||
! Read GW QP energies
|
||||
if (.not.pars%ncap) then
|
||||
call util_file_name('gwevals',.false.,.false.,gwevals)
|
||||
call util_file_name_resolve(gwevals,.false.)
|
||||
open(unit=unitno,status='old',form='unformatted',file=gwevals)
|
||||
do i=1,pars%ipol
|
||||
call sread(unitno, dbl_mb(pars%k_gw_evals+(i-1)*nbf_ao),
|
||||
& pars%nmo)
|
||||
enddo
|
||||
close(unit=unitno)
|
||||
endif
|
||||
|
||||
!3-center ERIs
|
||||
call int_init(pars%rtdb, 2, (/ao_bas_han, cd_bas_han/))
|
||||
call print_integrals((/ao_bas_han,cd_bas_han/),.false.)
|
||||
|
|
@ -205,5 +232,9 @@ c
|
|||
!Clean previous grid
|
||||
call grid_cleanup(.false.)
|
||||
|
||||
!Print expected memory requirements
|
||||
|
||||
!Compute three-center integrals
|
||||
call bse_ri_init(pars)
|
||||
|
||||
end subroutine bse_init
|
||||
|
|
|
|||
326
src/bsemol/bse_ri_init.F
Normal file
326
src/bsemol/bse_ri_init.F
Normal file
|
|
@ -0,0 +1,326 @@
|
|||
subroutine bse_ri_init(pars)
|
||||
c
|
||||
implicit none
|
||||
#include "errquit.fh"
|
||||
#include "mafdecls.fh"
|
||||
#include "cdft.fh"
|
||||
#include "bas.fh"
|
||||
#include "bse.fh"
|
||||
#include "stdio.fh"
|
||||
#include "util.fh"
|
||||
#include "global.fh"
|
||||
#ifdef SCALAPACK
|
||||
#include "scaleig.fh"
|
||||
#endif
|
||||
c
|
||||
type(bse_params_t) :: pars
|
||||
c
|
||||
character(len=*), parameter :: pname = "bse_ri_init: "
|
||||
c
|
||||
logical iolgc
|
||||
integer g_2ceris, g_2cinv, info, lmiss
|
||||
c
|
||||
integer ilo,ihi,jlo,jhi
|
||||
integer lbuf,lscr,ltmp,leri,lvecs
|
||||
integer l_buf,k_buf,l_scr,k_scr,l_tmp,k_tmp,k_eri,k_mo,l_mo
|
||||
integer l_vecs,k_vecs,l_idx,k_idx
|
||||
integer l_shlbuf, k_shlbuf, l_mobuf, k_mobuf
|
||||
integer me, isp, addr
|
||||
integer g_3ceris,shlo,shhi,offset,ld
|
||||
integer nmo,nvir,rihi,rilo,nri,nri_me,nocc,offk
|
||||
integer ish,jsh,ksh,ni,nj,nk,nij
|
||||
integer ifirst,jfirst,kfirst
|
||||
integer ilast,jlast,klast,imaxthreads
|
||||
integer aopairs,mopairs(2),idum,totmopairs
|
||||
logical IeqJ,ok,oactive
|
||||
double precision eri_est,timer,timing
|
||||
c
|
||||
double precision schwarz_shell
|
||||
external schwarz_shell
|
||||
c
|
||||
#ifdef USE_OPENMP
|
||||
integer,external :: omp_get_max_threads
|
||||
#endif
|
||||
|
||||
timer = util_wallsec()
|
||||
c
|
||||
c Create Global Arrays
|
||||
c
|
||||
if (.not.ga_create(mt_dbl,nbf_cd,nbf_cd,'2c ERI',
|
||||
$ 0, nbf_cd, g_2ceris))
|
||||
$ call errquit(pname//'Error creating 2c ERI',0,GA_ERR)
|
||||
if (.not.ga_create(mt_dbl,nbf_cd,nbf_cd,'2c INV',
|
||||
$ 0, nbf_cd, g_2cinv))
|
||||
$ call errquit(pname//'Error creating 2c ERI',0,GA_ERR)
|
||||
call ga_zero(g_2ceris)
|
||||
c
|
||||
c Compute 2-center ERIs
|
||||
c
|
||||
if (pars%me.eq.0) write(luout,9000)
|
||||
9000 format(/,10x,'Computing 2-center integrals')
|
||||
c
|
||||
call dft_get2eri(cd_bas_han, g_2ceris, .false.)
|
||||
c
|
||||
c Obtain inverse Cholesky factor
|
||||
c
|
||||
if (pars%me.eq.0) write(luout,9010)
|
||||
9010 format(10x,'Computing Inverse Cholesky factor')
|
||||
c
|
||||
call ga_sync()
|
||||
call ga_chol_inv_seq(g_2ceris, g_2cinv, info)
|
||||
if (info.ne.0) then
|
||||
if (pars%me.eq.0) then
|
||||
write(luout,*)
|
||||
write(luout,*)' Problem in performing Cholesky'
|
||||
write(luout,*)' Obtaining square root via truncated '
|
||||
write(luout,*)' eigenvalue decomposition '
|
||||
write(luout,*)
|
||||
endif
|
||||
call dft_invdiag(g_2ceris, g_2cinv, nbf_cd)
|
||||
endif
|
||||
c
|
||||
if (.not. ga_destroy(g_2ceris))
|
||||
$ call errquit(pname//'failed to destrou g_2ceri',0,GA_ERR)
|
||||
c
|
||||
c Allocate storage
|
||||
c
|
||||
nri = nbf_cd
|
||||
nmo = pars%nmo
|
||||
aopairs = (nbf_ao*(nbf_ao+1))/2
|
||||
|
||||
do isp=1,pars%ipol
|
||||
mopairs(isp) = pars%nmo*pars%nmo
|
||||
enddo
|
||||
totmopairs = sum(mopairs(1:pars%ipol))
|
||||
|
||||
if (.not.ga_create(mt_dbl,1,nri,'3c ERI',aopairs,0,
|
||||
$ g_3ceris))
|
||||
$ call errquit(pname//'can''t create 3c ERIs',0,GA_ERR)
|
||||
call ga_distribution(g_3ceris,pars%me,ilo,ihi,rilo,rihi)
|
||||
nri_me = rihi - rilo + 1
|
||||
ok = ga_destroy(g_3ceris)
|
||||
|
||||
c call ga_access(g_3ceris,ilo,ihi,rilo,rihi,k_eri,ld)
|
||||
c call dfill(aopairs*nri_me,0d0,dbl_mb(k_eri),1)
|
||||
|
||||
shlo = 0
|
||||
do ksh=1,nshells_cd
|
||||
ok = bas_cn2bfr(cd_bas_han,ksh,kfirst,klast)
|
||||
if (klast.lt.rilo) cycle
|
||||
if (shlo.eq.0) then
|
||||
offset = (rilo-kfirst)
|
||||
shlo = ksh
|
||||
endif
|
||||
shhi = ksh
|
||||
if (klast.ge.rihi) exit
|
||||
enddo
|
||||
|
||||
oactive = nri_me.gt.0
|
||||
|
||||
call int_mem_2e3c(lbuf, lscr)
|
||||
if (.not.ma_alloc_get(mt_dbl,lbuf, 'RI buffer',l_buf, k_buf))
|
||||
$ call errquit(pname//'can''t get buffer space',0,MA_ERR)
|
||||
lscr = max(lscr,nbf_cd*nbf_ao_mxnbf_cn**2)
|
||||
if (.not.ma_alloc_get(mt_dbl,lscr, 'RI scratch', l_scr, k_scr))
|
||||
$ call errquit(pname//'can''t get scratch space',0,MA_ERR)
|
||||
|
||||
if (.not.ma_alloc_get(mt_dbl,nri_me*nbf_ao_mxnbf_cn**2,
|
||||
$ 'Shell buffer',l_shlbuf, k_shlbuf))
|
||||
$ call errquit(pname//'can''t allocate shell buffer',0,MA_ERR)
|
||||
if (.not.ma_alloc_get(mt_dbl,
|
||||
$ nri_me*nbf_ao_mxnbf_cn*nmo*pars%ipol,
|
||||
$ 'MO buffer',l_mobuf, k_mobuf))
|
||||
$ call errquit(pname//'can''t allocate shell buffer',0,MA_ERR)
|
||||
c
|
||||
c Compute three-center ERIs in MO representation
|
||||
c
|
||||
call ga_sync()
|
||||
c
|
||||
if (pars%me.eq.0) write(luout,9020)
|
||||
9020 format(10x,'Computing 3-center integrals')
|
||||
|
||||
lvecs = nbf_ao*nmo*pars%ipol
|
||||
if (.not.ma_alloc_get(mt_dbl,lvecs,'MOVECS',l_vecs,k_vecs))
|
||||
$ call errquit(pname//'can''t get movecs space',0,MA_ERR)
|
||||
if (.not.ma_alloc_get(mt_dbl,totmopairs*nri_me,'MO Eris',
|
||||
$ l_mo,k_mo))
|
||||
$ call errquit(pname//'can''t get MO eris space',0,MA_ERR)
|
||||
call dfill(totmopairs*nri_me,0.0d0,dbl_mb(k_mo),1)
|
||||
|
||||
do isp=1,pars%ipol
|
||||
call ga_get(pars%g_movecs(isp),1,nbf_ao,1,nmo,
|
||||
$ dbl_mb(k_vecs+(isp-1)*nbf_ao*nmo),nbf_ao)
|
||||
enddo
|
||||
c
|
||||
#ifdef USE_OPENMP
|
||||
iMaxThreads = omp_get_max_threads()
|
||||
call util_blas_set_num_threads(iMaxThreads)
|
||||
#endif
|
||||
|
||||
do ish=1,nshells_ao
|
||||
ok = bas_cn2bfr(ao_bas_han,ish,ifirst,ilast)
|
||||
ni = ilast - ifirst + 1
|
||||
call dfill(ni*nmo*pars%ipol*nri_me,0d0,dbl_mb(k_mobuf),1)
|
||||
do jsh=1,nshells_ao
|
||||
eri_est = schwarz_shell(ish,jsh)
|
||||
if (eri_est.lt.pars%tol2e) cycle
|
||||
IeqJ = ish.eq.jsh
|
||||
ok = bas_cn2bfr(ao_bas_han,jsh,jfirst,jlast)
|
||||
nj = jlast - jfirst + 1
|
||||
nij = ni*nj
|
||||
do ksh=shlo,shhi
|
||||
ok = bas_cn2bfr(cd_bas_han, ksh, kfirst, klast)
|
||||
call int_2e3c(cd_bas_han,ksh,ao_bas_han,ish,jsh,lscr,
|
||||
$ dbl_mb(k_scr),lbuf,dbl_mb(k_buf))
|
||||
|
||||
if (ksh.eq.shlo) then
|
||||
kfirst = rilo
|
||||
offk = k_buf + offset*nij
|
||||
else
|
||||
offk = k_buf
|
||||
endif
|
||||
|
||||
if (ksh.eq.shhi) klast = rihi
|
||||
|
||||
nk = klast - kfirst + 1
|
||||
kfirst = kfirst - rilo + 1
|
||||
call dcopy(ni*nj*nk,dbl_mb(offk),1,
|
||||
$ dbl_mb(k_shlbuf+(kfirst-1)*ni*nj),1)
|
||||
enddo
|
||||
do isp=1,pars%ipol
|
||||
call gw_mo2(dbl_mb(k_shlbuf),
|
||||
$ dbl_mb(k_mobuf+(isp-1)*nmo*ni*nri_me),
|
||||
$ dbl_mb(k_vecs+(isp-1)*nbf_ao*nmo),nri_me,
|
||||
$ nbf_ao,ni,nj,nmo,nmo,ifirst,
|
||||
$ jfirst,ieqj)
|
||||
enddo
|
||||
enddo
|
||||
do isp=1,pars%ipol
|
||||
call gw_mo(dbl_mb(k_mobuf+(isp-1)*nmo*ni*nri_me),
|
||||
$ dbl_mb(k_mo+(isp-1)*mopairs(1)*nri_me),
|
||||
$ dbl_mb(k_vecs+(isp-1)*nbf_ao*nmo),nri_me,
|
||||
$ nbf_ao,ni,nmo,nmo,ifirst)
|
||||
enddo
|
||||
enddo
|
||||
c
|
||||
ok = ma_free_heap(l_buf) .and. ma_free_heap(l_scr) .and.
|
||||
$ ma_free_heap(l_shlbuf) .and.
|
||||
$ ma_free_heap(l_mobuf)
|
||||
if (.not.ok)
|
||||
$ call errquit(pname//'can''t free heap',165,MA_ERR)
|
||||
c
|
||||
c
|
||||
#ifdef GWDEBUG
|
||||
if (pars%me.eq.0) write(*,*) ' Transform ERIs'
|
||||
#endif
|
||||
|
||||
call ga_sync()
|
||||
|
||||
ok = ma_free_heap(l_vecs)
|
||||
if (.not.ok)
|
||||
$ call errquit(pname//'can''t free heap',176,MA_ERR)
|
||||
c
|
||||
call ga_sync()
|
||||
c
|
||||
c Distribute and orthonormalize ERIs
|
||||
c
|
||||
timer = util_wallsec()
|
||||
#ifdef GWDEBUG
|
||||
if (pars%me.eq.0) write(*,*) ' Distribute ERIs'
|
||||
#endif
|
||||
|
||||
addr = k_mo
|
||||
do isp=1,pars%ipol
|
||||
if (nmo.eq.0) cycle
|
||||
if (isp.eq.2) addr = k_mo + mopairs(1)*nri_me
|
||||
|
||||
! OO block
|
||||
call gw_puteris(dbl_mb(addr),nmo,nmo,nri,nri_me,
|
||||
$ rilo,1,pars%nocc(isp),1,pars%nocc(isp),
|
||||
$ pars%g_erioo(isp),'oo',pars%oolo(isp),
|
||||
$ pars%oohi(isp),idum)
|
||||
|
||||
! OV block
|
||||
call gw_puteris(dbl_mb(addr),nmo,nmo,nri,nri_me,
|
||||
$ rilo,1,pars%nocc(isp),pars%nocc(isp)+1,
|
||||
$ nmo,pars%g_eriov(isp),'ov',pars%ovlo(isp),
|
||||
$ pars%ovhi(isp),pars%mynpoles(isp))
|
||||
|
||||
! VV block
|
||||
call gw_puteris(dbl_mb(addr),nmo,nmo,nri,nri_me,
|
||||
$ rilo,pars%nocc(isp)+1,nmo,
|
||||
$ pars%nocc(isp)+1,nmo,pars%g_erivv(isp),
|
||||
$ 'vv',pars%vvlo(isp),pars%vvhi(isp),idum)
|
||||
|
||||
enddo
|
||||
|
||||
ok = ma_free_heap(l_mo)
|
||||
#ifdef GWDEBUG
|
||||
if (pars%me.eq.0) write(*,*) ' Read inverse'
|
||||
#endif
|
||||
c
|
||||
if (.not.ma_alloc_get(mt_dbl,nri**2,'RI temporary',l_tmp,k_tmp))
|
||||
$ call errquit(pname//'can''t get temp space',0,MA_ERR)
|
||||
if (pars%me.eq.0) then
|
||||
call ga_get(g_2cinv,1,nbf_cd,1,nbf_cd,dbl_mb(k_tmp),nbf_cd)
|
||||
endif
|
||||
call ga_sync()
|
||||
if (.not. ga_destroy(g_2cinv))
|
||||
$ call errquit(pname//'failed to destroy g_2cinv',0,GA_ERR)
|
||||
call ga_brdcst(1038,dbl_mb(k_tmp),ma_sizeof(mt_dbl,nri**2,
|
||||
$ mt_byte),0)
|
||||
|
||||
#ifdef GWDEBUG
|
||||
if (pars%me.eq.0) write(*,*) ' Orthogonalize ERIs'
|
||||
#endif
|
||||
do isp=1,pars%ipol
|
||||
if (nmo.eq.0) cycle
|
||||
c
|
||||
c Transform OO block
|
||||
c
|
||||
call ga_distribution(pars%g_erioo(isp),pars%me,ilo,ihi,jlo,jhi)
|
||||
if (jlo.gt.jhi) goto 101
|
||||
call ga_access(pars%g_erioo(isp),ilo,ihi,jlo,jhi,k_eri,ld)
|
||||
nij = jhi-jlo+1
|
||||
call dtrmm('l','l','n','n',nbf_cd,nij,1.0d0,dbl_mb(k_tmp),
|
||||
$ nbf_cd,dbl_mb(k_eri),nbf_cd)
|
||||
call ga_release_update(pars%g_erioo(isp),ilo,ihi,jlo,jhi)
|
||||
c
|
||||
c Transform OV block
|
||||
c
|
||||
101 continue
|
||||
call ga_distribution(pars%g_eriov(isp),pars%me,ilo,ihi,jlo,jhi)
|
||||
if (jlo.gt.jhi) goto 102
|
||||
call ga_access(pars%g_eriov(isp),ilo,ihi,jlo,jhi,k_eri,ld)
|
||||
nij = jhi-jlo+1
|
||||
call dtrmm('l','l','n','n',nbf_cd,nij,1.0d0,dbl_mb(k_tmp),
|
||||
$ nbf_cd,dbl_mb(k_eri),nbf_cd)
|
||||
call ga_release_update(pars%g_eriov(isp),ilo,ihi,jlo,jhi)
|
||||
c
|
||||
c Transform VV block
|
||||
c
|
||||
102 continue
|
||||
if (nmo.gt.pars%nocc(isp)) then
|
||||
call ga_distribution(pars%g_erivv(isp),pars%me,ilo,ihi,
|
||||
$ jlo,jhi)
|
||||
if (jlo.gt.jhi) goto 103
|
||||
call ga_access(pars%g_erivv(isp),ilo,ihi,jlo,jhi,k_eri,ld)
|
||||
nij = jhi-jlo+1
|
||||
call dtrmm('l','l','n','n',nbf_cd,nij,1.0d0,dbl_mb(k_tmp),
|
||||
$ nbf_cd,dbl_mb(k_eri),nbf_cd)
|
||||
call ga_release_update(pars%g_erivv(isp),ilo,ihi,jlo,jhi)
|
||||
endif
|
||||
103 continue
|
||||
enddo
|
||||
c
|
||||
c Free heap
|
||||
c
|
||||
ok = ma_free_heap(l_tmp)
|
||||
if (.not.ok)
|
||||
$ call errquit(pname//'failed to free heap',224,MA_ERR)
|
||||
c
|
||||
#ifdef USE_OPENMP
|
||||
call util_blas_set_num_threads(1)
|
||||
#endif
|
||||
return
|
||||
end subroutine
|
||||
47
src/bsemol/bse_wmn.F
Normal file
47
src/bsemol/bse_wmn.F
Normal file
|
|
@ -0,0 +1,47 @@
|
|||
subroutine bse_wmn(pars,w,nmo,nri,ipol)
|
||||
implicit none
|
||||
#include "errquit.fh"
|
||||
#include "mafdecls.fh"
|
||||
#include "global.fh"
|
||||
#include "bse.fh"
|
||||
|
||||
type(bse_params_t) :: pars
|
||||
integer nmo, nri, ipol
|
||||
|
||||
character(*), parameter :: pname = 'bse_wmn: '
|
||||
|
||||
double precision w(*)
|
||||
|
||||
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))
|
||||
|
||||
! transform occ-occ and vir-vir blocks in place
|
||||
call ga_distribution(pars%g_erioo(isp),pars%me,ilo,ihi,jlo,jhi)
|
||||
if (jlo.gt.jhi) goto 101
|
||||
call ga_access(pars%g_erioo(isp),ilo,ihi,jlo,jhi,keri,ld)
|
||||
call dtfsm('n','l','l','n','n',nri,jhi-jlo+1,1.0d0,w,
|
||||
& dbl_mb(keri),nri)
|
||||
call ga_release_update(pars%g_erioo(isp),1,nri,jlo,jhi)
|
||||
|
||||
101 continue
|
||||
call ga_distribution(pars%g_erivv(isp),pars%me,ilo,ihi,jlo,jhi)
|
||||
if (jlo.gt.jhi) cycle
|
||||
call ga_access(pars%g_erivv(isp),ilo,ihi,jlo,jhi,keri,ld)
|
||||
call dtfsm('n','l','l','n','n',nri,jhi-jlo+1,1.0d0,w,
|
||||
& dbl_mb(keri),nri)
|
||||
call ga_release_update(pars%g_erivv(isp),1,nri,jlo,jhi)
|
||||
enddo
|
||||
|
||||
end
|
||||
|
|
@ -12,6 +12,7 @@ c
|
|||
type(gw_params_t) :: pars
|
||||
c
|
||||
character(len=13), parameter :: pname = 'gw_analytic: '
|
||||
character*255 gwevals
|
||||
c
|
||||
integer,dimension(2) :: nocc,nvir,npoles,nstates,mynpoles
|
||||
|
||||
|
|
@ -45,6 +46,8 @@ c
|
|||
integer maxpoles,maxnpoles,isp,lSize
|
||||
|
||||
logical ok
|
||||
|
||||
integer,parameter :: unitno = 66
|
||||
c
|
||||
c -----------------------------------------------------------------
|
||||
c
|
||||
|
|
@ -185,15 +188,16 @@ c
|
|||
endif
|
||||
|
||||
c
|
||||
if (pars%evgw0.or.pars%evgw) then
|
||||
call gw_scissor(dbl_mb(k_Enew),dbl_mb(k_Eold),pars%nmo,
|
||||
call gw_scissor(dbl_mb(k_Enew),dbl_mb(k_Eold),pars%nmo,
|
||||
$ pars%nocc,pars%nvir,pars%noqp,pars%nvqp,
|
||||
$ pars%ipol,pars%me)
|
||||
|
||||
call ga_brdcst(1038,dbl_mb(k_enew),
|
||||
$ ma_sizeof(mt_dbl,nmo,mt_byte),0)
|
||||
call ga_brdcst(1038,dbl_mb(k_sigma),
|
||||
$ ma_sizeof(mt_dbl,nmo,mt_byte),0)
|
||||
call ga_brdcst(1038,dbl_mb(k_enew),
|
||||
$ ma_sizeof(mt_dbl,nmo,mt_byte),0)
|
||||
call ga_brdcst(1038,dbl_mb(k_sigma),
|
||||
$ ma_sizeof(mt_dbl,nmo,mt_byte),0)
|
||||
|
||||
if (pars%evgw0.or.pars%evgw) then
|
||||
|
||||
!
|
||||
! Update Fermi energy
|
||||
|
|
@ -214,6 +218,18 @@ c
|
|||
endif
|
||||
endif
|
||||
|
||||
!
|
||||
! Save eigenvalues
|
||||
!
|
||||
call gw_shift(dbl_mb(k_enew),dbl_mb(k_eold),-efermi,nmo,nmo,ipol)
|
||||
call util_file_name('gwevals',.false.,.false.,gwevals)
|
||||
call util_file_name_resolve(gwevals,.false.)
|
||||
open(unit=unitno,status='unknown',form='unformatted',file=gwevals)
|
||||
do i=1,ipol
|
||||
call swrite(unitno, dbl_mb(k_eold+(i-1)*nmo),nmo)
|
||||
enddo
|
||||
close(unit=unitno)
|
||||
|
||||
9050 format(10x,I3,6x,F8.3,6x,F6.3,6x,F8.3)
|
||||
c
|
||||
if (.not.ma_chop_stack(l_sigma))
|
||||
|
|
|
|||
|
|
@ -484,6 +484,7 @@ c if (.not.ga_create(mt_dbl,nri,npairs,erilabel,nri,0,eri))
|
|||
enddo
|
||||
|
||||
ok = ma_chop_stack(l_tmp)
|
||||
ok = ma_chop_stack(lmap)
|
||||
|
||||
call ga_distribution(eri,me,ld,ld,polelo,polehi)
|
||||
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue