move to 3d ga for DM

This commit is contained in:
edoapra 2023-07-17 22:22:26 -07:00
parent a02ba67801
commit 1f068b59ae
No known key found for this signature in database
GPG key ID: 48E12DA1EDE9E1B0
5 changed files with 123 additions and 79 deletions

View file

@ -68,13 +68,14 @@ c
integer alo(3), ahi(3), blo(3), bhi(3)
integer nfock_xc,i
integer ityp, ndim
logical lcgmin
logical lcgmin,l3d_arg
external xc_gotxc,xc_getipol
integer ga_create_atom_blocked,dft_npatch,dft_dm3d
external ga_create_atom_blocked,dft_dm3d
integer ga_create_atom_blocked,dft_npatch,dft_dm3d,dft_dm3d_calc5
external ga_create_atom_blocked,dft_dm3d,dft_dm3d_calc5
integer g_dens_3d,g_dens_save(2)
c
integer nroots !number of roots for the TDDFT gradient calculations
integer nroots !number of roots for the TDDFT gradient calculations
integer g_test
c
if(.not.xc_gotxc()) return
oprint= util_print('fock_xc',print_debug)
@ -275,6 +276,10 @@ c Use nroots to determine where the ground state density matrix is.
nroots = nfock/ndens
g_dens_scf(1) = int_mb(igxcd+nroots)
if (ndens.eq.2) g_dens_scf(2) = int_mb(igxcd+2*nroots+1)
g_dens_3d=dft_dm3d_calc5(g_dens,
D nbf_ao,geom,ao_bas_han,nfock,ndens)
g_dens_save(1)=g_dens(1)
g_dens(1)=g_dens_3d
else
g_dens_scf(1) = int_mb(igxcd)
if (ndens.eq.2) g_dens_scf(2) = int_mb(igxcd+1)
@ -303,6 +308,7 @@ c
endif
call ga_print(g_dens_scf(2))
endif
if(calc_type.ne.2) then
if(ga_nodeid().eq.0) then
write(luout,*) 'g_dens_xc1',int_mb(igxcd)
call util_flush(luout)
@ -323,6 +329,7 @@ c
endif
call ga_print(int_mb(igxcd+1))
endif
endif
endif
call xc_init_index(geom,ndens,nbf_ao,ao_bas_han,
& g_dens_scf, l_cntoce, k_cntoce,
@ -336,37 +343,31 @@ c Daniel (2-18-13): For calc_type = 5, we don't use a 3D global array
c and therefore don't want l3d to be set equal to true.
c if (.not.lcgmin) then
if (.not.lcgmin.and.calc_type.ne.5) then
call grid_quadv0_gen(rtdb,
cold & int_mb(igxcd),
l3d_arg=.true.
else
l3d_arg=l3d
endif
call grid_quadv0_gen(rtdb,
& g_dens,
& g_xc, nexc,
& rho_n_dum, Exc,
& nfock_xc, calc_type, tdum, .true., triplet)
else
call grid_quadv0_gen(rtdb,
& int_mb(igxcd),
& g_xc, nexc,
& rho_n_dum, Exc,
& nfock_xc, calc_type, tdum, l3d, triplet)
end if
& nfock_xc, calc_type, tdum, l3d_arg, triplet)
c
call dft_store_rhon(rho_n_dum)
c
c Clean up
c
if(l3d) then
if(calc_type.eq.2) then
g_dens(1)=g_dens_save(1)
if (.not. ga_destroy(g_dens_3d))
& call errquit('fock_xc: could not destroy g_dens_3d',0,
if(calc_type.eq.2.or.calc_type.eq.5) then
g_dens(1)=g_dens_save(1)
if (.not. ga_destroy(g_dens_3d))
& call errquit('fock_xc: could not destroy g_dens_3d',0,
& GA_ERR)
else
do i=1,dims(1)
if (.not. ga_destroy(int_mb(igxcd+i-1)))
& call errquit('fock_xc: could not destroy g_dens_xc',0,
& GA_ERR)
else
do i=1,dims(1)
if (.not. ga_destroy(int_mb(igxcd+i-1)))
& call errquit('fock_xc: could not destroy g_dens_xc',0,
& GA_ERR)
enddo
endif
enddo
endif
if (.not.ma_free_heap(lgxcd))
& call errquit('fockxc: cannot pop stack',0, MA_ERR)
@ -536,20 +537,56 @@ c call nga_inquire(g_dens, ityp, ndim, dims)
bhi(1) = 2*nfock_xc+1
call nga_copy_patch('N',g_dens,alo,ahi,dft_dm3d,blo,bhi)
endif
#if 1
alo(1) = 1
ahi(1) = nbf
alo(2) = 1
ahi(2) = nbf
blo(1) = nfock_xc+1
bhi(1) = nfock_xc+1
call nga_copy_patch('N',g_dens_scf,alo,ahi,dft_dm3d,blo,bhi)
call nga_copy_patch('N',g_dens_scf(1),alo,ahi,dft_dm3d,blo,bhi)
if(ndens.eq.2) then
blo(1) = 2*nfock_xc+2
bhi(1) = 2*nfock_xc+2
call nga_copy_patch('N',g_dens_scf(2),alo,ahi,dft_dm3d,blo,bhi)
endif
#endif
return
end
integer function dft_dm3d_calc5(g_dens,
, nbf,geom,basis,nfock,ndens)
implicit none
integer ndens,nfock
integer g_dens(*) ! array of 2-D GA
integer nbf
integer geom,basis
c
integer ga_create_atom_blocked
external ga_create_atom_blocked
integer alo(3), ahi(3), blo(3), bhi(3)
integer dims_3d, i
external ga_nodeid,ga_create3d_atom_blocked
integer ga_nodeid,ga_create3d_atom_blocked
c
dims_3d=ndens+nfock
dft_dm3d_calc5=ga_create3d_atom_blocked(geom, basis,
S 'xcdm3d_calc5', dims_3d)
c2d
alo(1) = 1
ahi(1) = nbf
alo(2) = 1
ahi(2) = nbf
alo(3) = 1
ahi(3) = 1
c3d
blo(2) = 1
bhi(2) = nbf
blo(3) = 1
bhi(3) = nbf
do i=1,nfock+ndens
blo(1)=i
bhi(1)=i
call nga_copy_patch('N',g_dens(i),
D alo,ahi,dft_dm3d_calc5,blo,bhi)
enddo
return
end
subroutine dft_store_rhon(dft_rho_n)

View file

@ -447,7 +447,8 @@ c
c Daniel (1-11-13): I think we can stay similar to calc_type = 2 for
c XC-third derivatives in this portion of the code, since this part
c doesn't depend on having 3-dimensional GAs.
if (calc_type.eq.2) then
if (calc_type.eq.2.or.calc_type.eq.5) then
call dfill(ipol*nctrs*nctrs, 0.0d0, dbl_mb(irdens_atom), 1)
do isp=1,ipol
do ii=1,nmat+1
g_dens_2d =
@ -461,32 +462,6 @@ c DM in last position is the scf DM
& call errquit(pname//'failed destroy',g_dens_2d,GA_ERR)
enddo
enddo
elseif (calc_type.eq.5) then
#if 1
do isp=1,ipol
call util_ga_mat_reduce(nbf_ao, nctrs, int_mb(icetobfr),
A g_dens(1+(isp-1)*(nmat+1)), 1,
A dbl_mb(irdens_atom+(isp-1)*nctrs*nctrs), 'absmax',
& dbl_mb(iscr), nbf_ao_mxnbf_ce,.true.)
do ii=2,(nmat+1)
call util_ga_mat_reduce(nbf_ao,nctrs,int_mb(icetobfr),
A g_dens(ii+(isp-1)*(nmat+1)), 1,
A dbl_mb(irdens_atom+(isp-1)*nctrs*nctrs), 'absmax',
& dbl_mb(iscr), nbf_ao_mxnbf_ce,.false.)
enddo
enddo
#else
do isp=1,ipol
call dfill(nctrs*nctrs, 0.0d0,
c dbl_mb(irdens_atom+(isp-1)*nctrs*nctrs),1)
do ii=1,nmat+1
call util_ga_mat_reduce(nbf_ao, nctrs, int_mb(icetobfr),
A g_dens(ii+(isp-1)*(nmat+1)), 1,
A dbl_mb(irdens_atom+(isp-1)*nctrs*nctrs), 'absmax',
& dbl_mb(iscr), nbf_ao_mxnbf_ce,.false.)
enddo
enddo
#endif
elseif (adft) then
call adft_reduce(nbf_xc, nctrs, int_mb(icetobfr), tmat, ipol,
& dbl_mb(irdens_atom), nbf_xc_mxnbf_ce)

View file

@ -253,6 +253,8 @@ c Functions
c
integer ga_create_atom_blocked
external ga_create_atom_blocked
integer g_dens_3d,dft_dm3d_calc5
external dft_dm3d_calc5
logical grid_reopen
external grid_reopen
logical xc_gotxc ! do we have a density functional?
@ -928,10 +930,25 @@ c the number of perturbed density matrices. This is necessary so that we
c can define a value for looping in xc_3rd_deriv and build the
c perturbed density in xc_rhogen. Since we loop over the roots,
c the number of Fock matrices is always equal to 1.
c
c 3d
c
g_dens_3d=dft_dm3d_calc5(g_dtmp,
D nao,ihdl_geom,ihdl_bfao,2*ipol,0)
if (xc_gotxc()) then
do ip = 1, 2*ipol
if (.not.ga_destroy(g_dtmp(ip)))
1 call errquit(pname//"could not destroy g_dtmp_ao", 0,
2 GA_ERR)
enddo
endif
call dftgh_gridv0(rtdb,ihdl_geom,ihdl_bfao, ipol,nao,
& g_dtmp, dbl_mb(k_frc_dft), dum, idum, calc_type,
& g_dens_3d, dbl_mb(k_frc_dft), dum, idum, calc_type,
& nat, log_mb(k_act), nat, dbl_mb(k_rdens_atom),
& int_mb(k_cetobfr), 1, otriplet, 0d0, .false.)
if (.not.ga_destroy(g_dens_3d))
+ call errquit(pname//"could not destroy g_dens_3d",0,GA_ERR)
c call dftgh_gridv0(rtdb,ihdl_geom,ihdl_bfao, ipol,nao,
c & g_dtmp, dbl_mb(k_frc_dft), dum, idum, calc_type,
c & nat, log_mb(k_act), nat, dbl_mb(k_rdens_atom),
@ -1366,13 +1383,6 @@ c
if (.not.ga_destroy(g_rhs_xc(1)))
+ call errquit(pname//"could not destroy g_rhs_xc",0,GA_ERR)
c
if (xc_gotxc()) then
do ip = 1, 2*ipol
if (.not.ga_destroy(g_dtmp(ip)))
1 call errquit(pname//"could not destroy g_dtmp_ao", 0,
2 GA_ERR)
enddo
endif
c
do ip = 1, ipol
if (.not.ga_destroy(g_d(0+ip)))

View file

@ -127,15 +127,27 @@ c TDDFT gradients.
c Use this for the ground state density matrix (for the gradient of
c the ground state density). Note that g_dens(2) and g_dens(4) store
c the ground state density matrices in TDDFT gradients.
if (calctype.eq.1) then
if (calctype.eq.1) then
#if 1
call xc_get3ddm(g_dens(1),
I 2*imat+ipol*(ispin-1),
% ifirst, ilast, 1, nbf_ao, ppp,nbfia)
#else
call ga_get(g_dens(2*imat+2*(ispin-1)),
1 ifirst, ilast, 1, nbf_ao, PPP, nbfia)
1 ifirst, ilast, 1, nbf_ao, PPP, nbfia)
#endif
c Use this for the perturbed density matrix (for the gradient of the
c perturbed density). Note that g_dens(1) and g_dens(3) store the
c penturbed density matrices in TDDFT gradients.
else if (calctype.eq.2) then
else if (calctype.eq.2) then
#if 1
call xc_get3ddm(g_dens(1),
I imat+ipol*(ispin-1),
% ifirst, ilast, 1, nbf_ao, ppp,nbfia)
#else
call ga_get(g_dens(imat+2*(ispin-1)),
1 ifirst, ilast, 1, nbf_ao, PPP, nbfia)
1 ifirst, ilast, 1, nbf_ao, PPP, nbfia)
#endif
endif
do mu = inizia, ifinia
c mu1 is used for indexing the scratch array hh (which starts at 1,

View file

@ -118,6 +118,7 @@ c
integer ga_nodeid
external ga_nodeid
#endif
integer type,ndim,dims(3)
g_keepd(1)=0
g_keepd(2)=0
call starttimer(monitor_xcrho)
@ -306,16 +307,25 @@ c
else
c3d
if(what.eq.1) then
c write(6,*) ' what ',what,ipert,npert
c write(6,*) ' iniz in',iniz
c write(6,'(a,4(i5))') ' what ',what,ipert,npert,ii
c write(6,*) ' iniz in',iniz,' ipt ',
c B ipert+(ii-1)*npert
c need to handle ttddfrgrad madness tddft_grad_compute_r.F
call nga_inquire(g_dens, type, ndim, dims)
if(ndim.lt.3.and.npert.gt.1) then
write(6,*) ' ndim ',ndim
write(6,*) ' dims ',dims
write(6,*) ' npert ',npert
call errquit(' 3dm fail',npert,0)
endif
call xc_get3ddm(g_dens(1),
I ipert+(ii-1)*npert,ipol*npert,
% ifirst, ilast, jfirst, jlast, Pmat,nbfia)
c write(6,*) ' iniz out',iniz
I ipert+(ii-1)*npert,
% ifirst, ilast, jfirst, jlast, Pmat,nbfia)
else
call ga_get(g_dens(ipert+(ii-1)*npert),
% ifirst, ilast, jfirst, jlast, Pmat,nbfia)
endif
% ifirst, ilast, jfirst, jlast, Pmat,nbfia)
endif
endif
call endtimer(monitor_gaget)
if (wmax*FUNC_MAXI*FUNC_MAXJ*
@ -860,13 +870,13 @@ c
return
end
C>
C> @}
subroutine xc_get3ddm(g_dens,ipert,npert,
C > @}
subroutine xc_get3ddm(g_dens,ipert,
% ifirst, ilast, jfirst, jlast, pmat,nbf)
implicit none
integer g_dens
integer ipert,npert,
% ifirst, ilast, jfirst, jlast
integer ipert
integer ifirst, ilast, jfirst, jlast
integer nbf
double precision pmat(*)
c