From 1f068b59ae442d074348f107a49c3193f7bc2ea2 Mon Sep 17 00:00:00 2001 From: edoapra Date: Mon, 17 Jul 2023 22:22:26 -0700 Subject: [PATCH] move to 3d ga for DM --- src/ddscf/fock_xc.F | 95 +++++++++++++------ src/nwdft/grid/grid_quadv0.F | 29 +----- .../lr_tddft_grad/tddft_grad_compute_g.F | 26 +++-- src/nwdft/lr_tddft_grad/tddft_grad_quadv0b.F | 20 +++- src/nwdft/xc/xc_rhogen.F | 32 ++++--- 5 files changed, 123 insertions(+), 79 deletions(-) diff --git a/src/ddscf/fock_xc.F b/src/ddscf/fock_xc.F index 920582cedd..7a5c87b137 100644 --- a/src/ddscf/fock_xc.F +++ b/src/ddscf/fock_xc.F @@ -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) diff --git a/src/nwdft/grid/grid_quadv0.F b/src/nwdft/grid/grid_quadv0.F index 02790b84ec..6632e960c1 100644 --- a/src/nwdft/grid/grid_quadv0.F +++ b/src/nwdft/grid/grid_quadv0.F @@ -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) diff --git a/src/nwdft/lr_tddft_grad/tddft_grad_compute_g.F b/src/nwdft/lr_tddft_grad/tddft_grad_compute_g.F index bf6b74d0fb..8aa1386eba 100644 --- a/src/nwdft/lr_tddft_grad/tddft_grad_compute_g.F +++ b/src/nwdft/lr_tddft_grad/tddft_grad_compute_g.F @@ -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))) diff --git a/src/nwdft/lr_tddft_grad/tddft_grad_quadv0b.F b/src/nwdft/lr_tddft_grad/tddft_grad_quadv0b.F index d106fcdba8..8722d02500 100644 --- a/src/nwdft/lr_tddft_grad/tddft_grad_quadv0b.F +++ b/src/nwdft/lr_tddft_grad/tddft_grad_quadv0b.F @@ -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, diff --git a/src/nwdft/xc/xc_rhogen.F b/src/nwdft/xc/xc_rhogen.F index 4ac0d9966b..c6b10409c8 100644 --- a/src/nwdft/xc/xc_rhogen.F +++ b/src/nwdft/xc/xc_rhogen.F @@ -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