From 1cbb583c2daefb7e89675ebc1b69b68398dc70c0 Mon Sep 17 00:00:00 2001 From: Daniel Mejia-Rodriguez Date: Sat, 13 Aug 2022 01:05:12 -0700 Subject: [PATCH] CS singlet --- src/bsemol/GNUmakefile | 7 +- src/bsemol/bse.fh | 3 +- src/bsemol/bse_analytic.F | 68 ++++++++ src/bsemol/bse_buildw.F | 75 +++++++++ src/bsemol/bse_casida.F | 191 ++++++++++++++++++++++ src/bsemol/bse_init.F | 33 +++- src/bsemol/bse_ri_init.F | 326 ++++++++++++++++++++++++++++++++++++++ src/bsemol/bse_wmn.F | 47 ++++++ src/gwmol/gw_analytic.F | 28 +++- src/gwmol/gw_ri_init.F | 1 + 10 files changed, 768 insertions(+), 11 deletions(-) create mode 100644 src/bsemol/bse_buildw.F create mode 100644 src/bsemol/bse_casida.F create mode 100644 src/bsemol/bse_ri_init.F create mode 100644 src/bsemol/bse_wmn.F diff --git a/src/bsemol/GNUmakefile b/src/bsemol/GNUmakefile index 6f695705a1..63bbd26bf6 100644 --- a/src/bsemol/GNUmakefile +++ b/src/bsemol/GNUmakefile @@ -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 diff --git a/src/bsemol/bse.fh b/src/bsemol/bse.fh index 19d3f5f2a3..9f177f729f 100644 --- a/src/bsemol/bse.fh +++ b/src/bsemol/bse.fh @@ -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 diff --git a/src/bsemol/bse_analytic.F b/src/bsemol/bse_analytic.F index 764ceae421..85294a6e5c 100644 --- a/src/bsemol/bse_analytic.F +++ b/src/bsemol/bse_analytic.F @@ -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 diff --git a/src/bsemol/bse_buildw.F b/src/bsemol/bse_buildw.F new file mode 100644 index 0000000000..47d48e3532 --- /dev/null +++ b/src/bsemol/bse_buildw.F @@ -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 diff --git a/src/bsemol/bse_casida.F b/src/bsemol/bse_casida.F new file mode 100644 index 0000000000..53849aa2ed --- /dev/null +++ b/src/bsemol/bse_casida.F @@ -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 diff --git a/src/bsemol/bse_init.F b/src/bsemol/bse_init.F index ee0353b6bc..fd893ffe4c 100644 --- a/src/bsemol/bse_init.F +++ b/src/bsemol/bse_init.F @@ -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 diff --git a/src/bsemol/bse_ri_init.F b/src/bsemol/bse_ri_init.F new file mode 100644 index 0000000000..385e78fa2e --- /dev/null +++ b/src/bsemol/bse_ri_init.F @@ -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 diff --git a/src/bsemol/bse_wmn.F b/src/bsemol/bse_wmn.F new file mode 100644 index 0000000000..0ed4e09241 --- /dev/null +++ b/src/bsemol/bse_wmn.F @@ -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 diff --git a/src/gwmol/gw_analytic.F b/src/gwmol/gw_analytic.F index 23e497990c..e558cc8e0d 100644 --- a/src/gwmol/gw_analytic.F +++ b/src/gwmol/gw_analytic.F @@ -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)) diff --git a/src/gwmol/gw_ri_init.F b/src/gwmol/gw_ri_init.F index 08b6e443f3..5621749465 100644 --- a/src/gwmol/gw_ri_init.F +++ b/src/gwmol/gw_ri_init.F @@ -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)