diff --git a/src/grid/grid_integrate.F b/src/grid/grid_integrate.F index da97b92dc7..0d4b4be377 100644 --- a/src/grid/grid_integrate.F +++ b/src/grid/grid_integrate.F @@ -3,7 +3,6 @@ ! Copyright (C) 2000 - 2020 CP2K developers group ! !--------------------------------------------------------------------------------------------------! MODULE grid_integrate - USE ao_util, ONLY: exp_radius_very_extended USE cell_types, ONLY: cell_type USE cube_utils, ONLY: compute_cube_center,& cube_info_type,& @@ -71,16 +70,13 @@ END SUBROUTINE call_to_xyz_to_vab !> \param pab ... !> \param o1 ... !> \param o2 ... -!> \param eps_gvg_rspace ... +!> \param radius ... !> \param calculate_forces ... !> \param hdab ... !> \param hadb ... !> \param force_a ... !> \param force_b ... !> \param compute_tau ... -!> \param map_consistent ... -!> \param collocate_rho0 ... -!> \param rpgf0_s ... !> \param use_virial ... !> \param my_virial_a ... !> \param my_virial_b ... @@ -92,10 +88,10 @@ SUBROUTINE integrate_pgf_product(la_max, zeta, la_min, & lb_max, zetb, lb_min, & ra, rab, rsgrid, cell, & cube_info, hab, pab, o1, o2, & - eps_gvg_rspace, & + radius, & calculate_forces, hdab, hadb, force_a, force_b, & - compute_tau, map_consistent, & - collocate_rho0, rpgf0_s, use_virial, my_virial_a, & + compute_tau, & + use_virial, my_virial_a, & my_virial_b, a_hdab, use_subpatch, subpatch_pattern) INTEGER, INTENT(IN) :: la_max @@ -111,15 +107,13 @@ SUBROUTINE integrate_pgf_product(la_max, zeta, la_min, & REAL(KIND=dp), DIMENSION(:, :), & OPTIONAL, POINTER :: pab INTEGER, INTENT(IN) :: o1, o2 - REAL(KIND=dp), INTENT(IN) :: eps_gvg_rspace + REAL(KIND=dp), INTENT(IN) :: radius LOGICAL, INTENT(IN) :: calculate_forces REAL(KIND=dp), DIMENSION(:, :, :), & OPTIONAL, POINTER :: hdab, hadb REAL(KIND=dp), DIMENSION(3), & INTENT(INOUT), OPTIONAL :: force_a, force_b - LOGICAL, INTENT(IN), OPTIONAL :: compute_tau, map_consistent, & - collocate_rho0 - REAL(dp), INTENT(IN), OPTIONAL :: rpgf0_s + LOGICAL, INTENT(IN), OPTIONAL :: compute_tau LOGICAL, INTENT(IN), OPTIONAL :: use_virial REAL(KIND=dp), DIMENSION(3, 3), OPTIONAL :: my_virial_a, my_virial_b REAL(KIND=dp), DIMENSION(:, :, :, :), OPTIONAL, POINTER :: a_hdab @@ -136,15 +130,11 @@ SUBROUTINE integrate_pgf_product(la_max, zeta, la_min, & INTEGER, DIMENSION(3) :: cubecenter, lb_cube, ng, & ub_cube INTEGER, DIMENSION(:), POINTER :: sphere_bounds - LOGICAL :: my_collocate_rho0, & - my_compute_tau, & - my_map_consistent, & + LOGICAL :: my_compute_tau, & my_use_virial, & subpatch_integrate - REAL(KIND=dp) :: axpm0, cutoff, & - der_a(3), der_b(3), f, ftza, ftzb, pabval, pg, & - prefactor, radius, rpg, zetp - REAL(KIND=dp), DIMENSION(3) :: dr, rap, rb, rbp, roffset, rp + REAL(KIND=dp) :: axpm0, der_a(3), der_b(3), f, ftza, ftzb, pabval, pg, prefactor, rpg, zetp + REAL(KIND=dp), DIMENSION(3) :: dr, rap, rb, roffset, rp REAL(KIND=dp), DIMENSION(:, :, :), & POINTER :: grid @@ -183,19 +173,6 @@ SUBROUTINE integrate_pgf_product(la_max, zeta, la_min, & my_compute_tau = .FALSE. ENDIF - ! use identical radii for integrate and collocate ? - IF (PRESENT(map_consistent)) THEN - my_map_consistent = map_consistent - ELSE - my_map_consistent = .FALSE. - ENDIF - - IF (PRESENT(collocate_rho0) .AND. PRESENT(rpgf0_s)) THEN - my_collocate_rho0 = collocate_rho0 - ELSE - my_collocate_rho0 = .FALSE. - END IF - IF (calculate_forces) THEN la_max_local = la_max + 1 ! needed for the derivative of the gaussian, unimportant which one la_min_local = MAX(la_min - 1, 0) ! just in case the la_min,lb_min is not zero @@ -225,29 +202,9 @@ SUBROUTINE integrate_pgf_product(la_max, zeta, la_min, & prefactor = EXP(-zeta*f*DOT_PRODUCT(rab, rab)) ! *** position of the gaussian product rap(:) = f*rab(:) - rbp(:) = rap(:) - rab(:) rp(:) = ra(:) + rap(:) ! this is the gaussian center in real coordinates rb(:) = ra(:) + rab(:) - IF (my_map_consistent) THEN ! still assumes that eps_gvg_rspace=eps_rho_rspace - cutoff = 1.0_dp - radius = exp_radius_very_extended(la_min, la_max, lb_min, lb_max, ra=ra, rb=rb, rp=rp, & - zetp=zetp, eps=eps_gvg_rspace, prefactor=prefactor, cutoff=cutoff) - ELSE IF (my_collocate_rho0) THEN - cutoff = 0.0_dp - prefactor = 1.0_dp - radius = rpgf0_s - ELSE - cutoff = 1.0_dp - IF (PRESENT(pab)) THEN - radius = exp_radius_very_extended(la_min, la_max, lb_min, lb_max, pab, o1, o2, ra, rb, rp, & - zetp, eps_gvg_rspace, prefactor, cutoff) - ELSE - radius = exp_radius_very_extended(la_min, la_max, lb_min, lb_max, ra=ra, rb=rb, rp=rp, & - zetp=zetp, eps=eps_gvg_rspace, prefactor=prefactor, cutoff=cutoff) - ENDIF - ENDIF - IF (radius == 0.0_dp) THEN RETURN ENDIF diff --git a/src/hirshfeld_methods.F b/src/hirshfeld_methods.F index 37eeee8da7..d3aa13d2f0 100644 --- a/src/hirshfeld_methods.F +++ b/src/hirshfeld_methods.F @@ -543,7 +543,7 @@ SUBROUTINE hirshfeld_integration(qs_env, hirshfeld_env, rfun, fval, fderiv) INTEGER, ALLOCATABLE, DIMENSION(:) :: cores INTEGER, DIMENSION(:), POINTER :: atom_list LOGICAL :: do_force - REAL(KIND=dp) :: alpha, coef, dvol, eps_rho_rspace + REAL(KIND=dp) :: alpha, coef, dvol, eps_rho_rspace, radius REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, ra REAL(KIND=dp), DIMENSION(:, :), POINTER :: hab, pab TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set @@ -610,11 +610,17 @@ SUBROUTINE hirshfeld_integration(qs_env, hirshfeld_env, rfun, fval, fderiv) hab(1, 1) = 0.0_dp force_a(:) = 0.0_dp force_b(:) = 0.0_dp - ! + + radius = exp_radius_very_extended(la_min=0, la_max=0, lb_min=0, lb_max=0, & + ra=ra, rb=ra, rp=ra, & + zetp=alpha, eps=eps_rho_rspace, & + pab=pab, o1=0, o2=0, & ! without map_consistent + prefactor=1.0_dp, cutoff=1.0_dp) + CALL integrate_pgf_product(0, alpha, 0, & 0, 0.0_dp, 0, ra, (/0.0_dp, 0.0_dp, 0.0_dp/), & rs_v, cell, pw_env%cube_info(1), hab, pab=pab, o1=0, o2=0, & - eps_gvg_rspace=eps_rho_rspace, calculate_forces=do_force, & + radius=radius, calculate_forces=do_force, & force_a=force_a, force_b=force_b, use_virial=.FALSE., & use_subpatch=.TRUE., subpatch_pattern=0_int_8) fval(atom_a) = fval(atom_a) + hab(1, 1)*dvol*coef diff --git a/src/mp2_eri_gpw.F b/src/mp2_eri_gpw.F index 8b0739f4b5..54ccb31231 100644 --- a/src/mp2_eri_gpw.F +++ b/src/mp2_eri_gpw.F @@ -9,6 +9,7 @@ !> 07.2019 separated from mp2_integrals.F [Frederick Stein] ! ************************************************************************************************** MODULE mp2_eri_gpw + USE ao_util, ONLY: exp_radius_very_extended USE atomic_kind_types, ONLY: atomic_kind_type USE basis_set_types, ONLY: gto_basis_set_type USE cell_types, ONLY: cell_type,& @@ -182,10 +183,10 @@ SUBROUTINE mp2_eri_2c_integrate_gpw(qs_env, para_env_sub, dimen_RI, mo_coeff, my INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgfa, nsgfa INTEGER, DIMENSION(:, :), POINTER :: first_sgfa LOGICAL :: map_it_here - REAL(KIND=dp) :: cutoff_old, relative_cutoff_old + REAL(KIND=dp) :: cutoff_old, radius, relative_cutoff_old REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: e_cutoff_old, wf_vector REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: I_ab - REAL(KIND=dp), DIMENSION(3) :: ra, rab + REAL(KIND=dp), DIMENSION(3) :: ra REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a REAL(KIND=dp), DIMENSION(:, :), POINTER :: I_tmp2, rpgfa, sphi_a, zeta TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set @@ -260,7 +261,6 @@ SUBROUTINE mp2_eri_2c_integrate_gpw(qs_env, para_env_sub, dimen_RI, mo_coeff, my zeta => basis_set_a%zet ra(:) = pbc(particle_set(iatom)%r, cell) - rab = 0.0_dp DO iset = 1, nseta ncoa = npgfa(iset)*ncoset(la_max(iset)) @@ -303,23 +303,27 @@ SUBROUTINE mp2_eri_2c_integrate_gpw(qs_env, para_env_sub, dimen_RI, mo_coeff, my IF (map_it_here) THEN DO ipgf = 1, npgfa(iset) sgfa = first_sgfa(1, iset) - na1 = (ipgf - 1)*ncoset(la_max(iset)) + 1 na2 = ipgf*ncoset(la_max(iset)) igrid_level = gaussian_gridlevel(pw_env_sub%gridlevel_info, zeta(ipgf, iset)) + radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), & + lb_min=0, lb_max=0, ra=ra, rb=ra, rp=ra, & + zetp=zeta(ipgf, iset), & + eps=dft_control%qs_control%eps_gvg_rspace, & + prefactor=1.0_dp, cutoff=1.0_dp) + CALL integrate_pgf_product( & la_max=la_max(iset), zeta=zeta(ipgf, iset)/2.0_dp, la_min=la_min(iset), & lb_max=0, zetb=zeta(ipgf, iset)/2.0_dp, lb_min=0, & - ra=ra, rab=rab, & + ra=ra, rab=(/0.0_dp, 0.0_dp, 0.0_dp/), & rsgrid=rs_v(igrid_level)%rs_grid, & cell=cell, & cube_info=pw_env_sub%cube_info(igrid_level), & hab=I_tmp2, & o1=na1 - 1, & o2=0, & - map_consistent=.TRUE., & - eps_gvg_rspace=dft_control%qs_control%eps_gvg_rspace, & + radius=radius, & calculate_forces=.FALSE.) END DO diff --git a/src/mp2_ri_grad.F b/src/mp2_ri_grad.F index 88b91e97ec..7ccfb0c8b4 100644 --- a/src/mp2_ri_grad.F +++ b/src/mp2_ri_grad.F @@ -9,6 +9,7 @@ !> 10.2013 created [Mauro Del Ben] ! ************************************************************************************************** MODULE mp2_ri_grad + USE ao_util, ONLY: exp_radius_very_extended USE atomic_kind_types, ONLY: atomic_kind_type,& get_atomic_kind_set USE basis_set_types, ONLY: gto_basis_set_type @@ -171,11 +172,11 @@ SUBROUTINE calc_ri_mp2_nonsep(qs_env, mp2_env, para_env, para_env_sub, cell, par LOGICAL :: alpha_beta, map_it_here, skip_shell, & use_virial REAL(KIND=dp) :: cutoff_old, e_hartree, eps_filter, & - factor, omega, pair_energy, & + factor, omega, pair_energy, radius, & relative_cutoff_old, total_rho REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: e_cutoff_old, wf_vector REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: G_PQ_local, I_ab - REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, ra, rab + REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, ra REAL(KIND=dp), DIMENSION(3, 3) :: h_stress, my_virial_a, my_virial_b REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a REAL(KIND=dp), DIMENSION(:, :), POINTER :: I_tmp2, pab, rpgfa, sphi_a, zeta @@ -507,7 +508,6 @@ SUBROUTINE calc_ri_mp2_nonsep(qs_env, mp2_env, para_env, para_env_sub, cell, par zeta => basis_set_a%zet ra(:) = pbc(particle_set(iatom)%r, cell) - rab = 0.0_dp force_a(:) = 0.0_dp force_b(:) = 0.0_dp @@ -571,9 +571,15 @@ SUBROUTINE calc_ri_mp2_nonsep(qs_env, mp2_env, para_env, para_env_sub, cell, par na2 = ipgf*ncoset(la_max(iset)) igrid_level = gaussian_gridlevel(pw_env_sub%gridlevel_info, zeta(ipgf, iset)) + radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), & + lb_min=0, lb_max=0, ra=ra, rb=ra, rp=ra, & + zetp=zeta(ipgf, iset), & + eps=dft_control%qs_control%eps_gvg_rspace, & + prefactor=1.0_dp, cutoff=1.0_dp) + CALL integrate_pgf_product(la_max=la_max(iset), zeta=zeta(ipgf, iset)/2.0_dp, la_min=la_min(iset), & lb_max=0, zetb=zeta(ipgf, iset)/2.0_dp, lb_min=0, & - ra=ra, rab=rab, & + ra=ra, rab=(/0.0_dp, 0.0_dp, 0.0_dp/), & rsgrid=rs_v(igrid_level)%rs_grid, & cell=cell, & cube_info=pw_env_sub%cube_info(igrid_level), & @@ -581,8 +587,7 @@ SUBROUTINE calc_ri_mp2_nonsep(qs_env, mp2_env, para_env, para_env_sub, cell, par pab=pab, & o1=na1 - 1, & o2=0, & - map_consistent=.TRUE., & - eps_gvg_rspace=dft_control%qs_control%eps_gvg_rspace, & + radius=radius, & calculate_forces=.TRUE., & force_a=force_a, force_b=force_b, & use_virial=use_virial, my_virial_a=my_virial_a, my_virial_b=my_virial_b) @@ -744,7 +749,6 @@ SUBROUTINE calc_ri_mp2_nonsep(qs_env, mp2_env, para_env, para_env_sub, cell, par zeta => basis_set_a%zet ra(:) = pbc(particle_set(iatom)%r, cell) - rab = 0.0_dp force_a(:) = 0.0_dp force_b(:) = 0.0_dp @@ -821,9 +825,15 @@ SUBROUTINE calc_ri_mp2_nonsep(qs_env, mp2_env, para_env, para_env_sub, cell, par na2 = ipgf*ncoset(la_max(iset)) igrid_level = gaussian_gridlevel(pw_env_sub%gridlevel_info, zeta(ipgf, iset)) + radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), & + lb_min=0, lb_max=0, ra=ra, rb=ra, rp=ra, & + zetp=zeta(ipgf, iset), & + eps=dft_control%qs_control%eps_gvg_rspace, & + prefactor=1.0_dp, cutoff=1.0_dp) + CALL integrate_pgf_product(la_max=la_max(iset), zeta=zeta(ipgf, iset)/2.0_dp, la_min=la_min(iset), & lb_max=0, zetb=zeta(ipgf, iset)/2.0_dp, lb_min=0, & - ra=ra, rab=rab, & + ra=ra, rab=(/0.0_dp, 0.0_dp, 0.0_dp/), & rsgrid=rs_v(igrid_level)%rs_grid, & cell=cell, & cube_info=pw_env_sub%cube_info(igrid_level), & @@ -831,8 +841,7 @@ SUBROUTINE calc_ri_mp2_nonsep(qs_env, mp2_env, para_env, para_env_sub, cell, par pab=pab, & o1=na1 - 1, & o2=0, & - map_consistent=.TRUE., & - eps_gvg_rspace=dft_control%qs_control%eps_gvg_rspace, & + radius=radius, & calculate_forces=.TRUE., & force_a=force_a, force_b=force_b, & use_virial=use_virial, my_virial_a=my_virial_a, my_virial_b=my_virial_b) diff --git a/src/qmmm_image_charge.F b/src/qmmm_image_charge.F index ae34671dcf..0a2a3a9d36 100644 --- a/src/qmmm_image_charge.F +++ b/src/qmmm_image_charge.F @@ -10,7 +10,7 @@ !> \author Dorothea Golze ! ************************************************************************************************** MODULE qmmm_image_charge - + USE ao_util, ONLY: exp_radius_very_extended USE cell_types, ONLY: cell_type,& pbc USE cp_control_types, ONLY: dft_control_type @@ -355,7 +355,7 @@ SUBROUTINE integrate_potential_ga_rspace(potential, qmmm_env, qs_env, int_res, & INTEGER :: atom_a, atom_b, atom_ref, handle, iatom, & j, k, natom, npme INTEGER, DIMENSION(:), POINTER :: cores - REAL(KIND=dp) :: eps_rho_rspace + REAL(KIND=dp) :: eps_rho_rspace, radius REAL(KIND=dp), DIMENSION(3) :: ra REAL(KIND=dp), DIMENSION(:, :), POINTER :: hab TYPE(cell_type), POINTER :: cell @@ -425,10 +425,15 @@ SUBROUTINE integrate_potential_ga_rspace(potential, qmmm_env, qs_env, int_res, & hab(1, 1) = 0.0_dp + radius = exp_radius_very_extended(la_min=0, la_max=0, lb_min=0, lb_max=0, & + ra=ra, rb=ra, rp=ra, & + zetp=qmmm_env%image_charge_pot%eta, eps=eps_rho_rspace, & + prefactor=1.0_dp, cutoff=1.0_dp) + CALL integrate_pgf_product(0, qmmm_env%image_charge_pot%eta, 0, & 0, 0.0_dp, 0, ra, (/0.0_dp, 0.0_dp, 0.0_dp/), & rs_v, cell, pw_env%cube_info(1), hab, o1=0, o2=0, & - eps_gvg_rspace=eps_rho_rspace, calculate_forces=.FALSE., & + radius=radius, calculate_forces=.FALSE., & use_subpatch=.TRUE., subpatch_pattern=0_int_8) int_res(iatom) = hab(1, 1) @@ -470,7 +475,7 @@ SUBROUTINE integrate_potential_devga_rspace(potential, coeff, forces, qmmm_env, INTEGER :: atom_a, handle, iatom, j, natom, npme INTEGER, DIMENSION(:), POINTER :: cores LOGICAL :: use_virial - REAL(KIND=dp) :: eps_rho_rspace + REAL(KIND=dp) :: eps_rho_rspace, radius REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, ra REAL(KIND=dp), DIMENSION(:, :), POINTER :: hab, pab TYPE(cell_type), POINTER :: cell @@ -545,10 +550,16 @@ SUBROUTINE integrate_potential_devga_rspace(potential, coeff, forces, qmmm_env, force_a(:) = 0.0_dp force_b(:) = 0.0_dp + radius = exp_radius_very_extended(la_min=0, la_max=0, lb_min=0, lb_max=0, & + ra=ra, rb=ra, rp=ra, & + zetp=qmmm_env%image_charge_pot%eta, eps=eps_rho_rspace, & + pab=pab, o1=0, o2=0, & ! without map_consistent + prefactor=1.0_dp, cutoff=1.0_dp) + CALL integrate_pgf_product(0, qmmm_env%image_charge_pot%eta, 0, & 0, 0.0_dp, 0, ra, (/0.0_dp, 0.0_dp, 0.0_dp/), & rs_v, cell, pw_env%cube_info(1), hab, pab, o1=0, o2=0, & - eps_gvg_rspace=eps_rho_rspace, calculate_forces=.TRUE., & + radius=radius, calculate_forces=.TRUE., & force_a=force_a, force_b=force_b, use_subpatch=.TRUE., & subpatch_pattern=0_int_8) diff --git a/src/qs_integrate_potential_product.F b/src/qs_integrate_potential_product.F index d8e7dee832..c3091114a1 100644 --- a/src/qs_integrate_potential_product.F +++ b/src/qs_integrate_potential_product.F @@ -28,6 +28,7 @@ ! ************************************************************************************************** MODULE qs_integrate_potential_product USE admm_types, ONLY: admm_type + USE ao_util, ONLY: exp_radius_very_extended USE atomic_kind_types, ONLY: atomic_kind_type,& get_atomic_kind_set USE basis_set_types, ONLY: get_gto_basis_set,& @@ -209,9 +210,10 @@ SUBROUTINE integrate_v_rspace_low(v_rspace, hmat, hmat_kp, pmat, pmat_kp, qs_env LOGICAL :: atom_pair_changed, atom_pair_done, distributed_grids, do_kp, found, h_duplicated, & has_threads, my_compute_tau, my_force_adm, my_gapw, new_set_pair_coming, p_duplicated, & pab_required, scatter, use_subpatch, use_virial - REAL(KIND=dp) :: admm_scal_fac, dab, eps_rho_rspace, & - scalef, zetp - REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, ra, rab, rab_inv, rb + REAL(KIND=dp) :: admm_scal_fac, eps_rho_rspace, f, & + prefactor, radius, scalef, zetp + REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, ra, rab, rab_inv, rb, & + rp REAL(KIND=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b, pv_thread REAL(KIND=dp), DIMENSION(3, natom) :: force_thread REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b @@ -465,7 +467,7 @@ SUBROUTINE integrate_v_rspace_low(v_rspace, hmat, hmat_kp, pmat, pmat_kp, qs_env !$OMP PRIVATE(rpgfa,set_radius_a,sphi_a,zeta,first_sgfb,lb_max,lb_min,npgfb), & !$OMP PRIVATE(nsetb,nsgfb,rpgfb,set_radius_b,sphi_b,zetb,found), & !$OMP PRIVATE(force_a,force_b,my_virial_a,my_virial_b,atom_pair_changed,h_block), & -!$OMP PRIVATE(p_block,ncoa,sgfa,ncob,sgfb,rab,ra,rb,zetp,dab,igrid_level), & +!$OMP PRIVATE(p_block,ncoa,sgfa,ncob,sgfb,rab,ra,rb,rp,zetp,f,prefactor,radius,igrid_level), & !$OMP PRIVATE(na1,na2,nb1,nb2,use_subpatch,rab_inv,new_set_pair_coming,atom_pair_done), & !$OMP PRIVATE(iset_new,jset_new,ipgf_new,jpgf_new,scalef), & !$OMP PRIVATE(itask,dist,has_threads) & @@ -640,11 +642,17 @@ SUBROUTINE integrate_v_rspace_low(v_rspace, hmat, hmat_kp, pmat, pmat_kp, qs_env ENDIF rab = tasks(itask)%rab - rb(1) = ra(1) + rab(1) - rb(2) = ra(2) + rab(2) - rb(3) = ra(3) + rab(3) + rb(:) = ra(:) + rab(:) zetp = zeta(ipgf, iset) + zetb(jpgf, jset) - dab = NORM2(rab) + + f = zetb(jpgf, jset)/zetp + rp(:) = ra(:) + f*rab(:) + prefactor = EXP(-zeta(ipgf, iset)*f*DOT_PRODUCT(rab, rab)) + radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), & + lb_min=lb_min(jset), lb_max=lb_max(jset), & + ra=ra, rb=rb, rp=rp, & + zetp=zetp, eps=eps_rho_rspace, & + prefactor=prefactor, cutoff=1.0_dp) na1 = (ipgf - 1)*ncoset(la_max(iset)) + 1 na2 = ipgf*ncoset(la_max(iset)) @@ -671,10 +679,10 @@ SUBROUTINE integrate_v_rspace_low(v_rspace, hmat, hmat_kp, pmat, pmat_kp, qs_env ra, rab, rs_v(igrid_level)%rs_grid, cell, & cube_info(igrid_level), & hab, pab=pab, o1=na1 - 1, o2=nb1 - 1, & - eps_gvg_rspace=eps_rho_rspace, & + radius=radius, & calculate_forces=calculate_forces, & force_a=force_a, force_b=force_b, & - compute_tau=my_compute_tau, map_consistent=.TRUE., & + compute_tau=my_compute_tau, & use_virial=use_virial, my_virial_a=my_virial_a, & my_virial_b=my_virial_b, use_subpatch=use_subpatch, subpatch_pattern=tasks(itask)%subpatch_pattern) ELSE @@ -685,10 +693,10 @@ SUBROUTINE integrate_v_rspace_low(v_rspace, hmat, hmat_kp, pmat, pmat_kp, qs_env rb, rab_inv, rs_v(igrid_level)%rs_grid, cell, & cube_info(igrid_level), & hab, pab=pab, o1=nb1 - 1, o2=na1 - 1, & - eps_gvg_rspace=eps_rho_rspace, & + radius=radius, & calculate_forces=calculate_forces, & force_a=force_b, force_b=force_a, & - compute_tau=my_compute_tau, map_consistent=.TRUE., & + compute_tau=my_compute_tau, & use_virial=use_virial, my_virial_a=my_virial_b, & my_virial_b=my_virial_a, use_subpatch=use_subpatch, subpatch_pattern=tasks(itask)%subpatch_pattern) END IF @@ -700,11 +708,11 @@ SUBROUTINE integrate_v_rspace_low(v_rspace, hmat, hmat_kp, pmat, pmat_kp, qs_env ra, rab, rs_v(igrid_level)%rs_grid, cell, & cube_info(igrid_level), & hab, o1=na1 - 1, o2=nb1 - 1, & - eps_gvg_rspace=eps_rho_rspace, & + radius=radius, & calculate_forces=calculate_forces, & force_a=force_a, force_b=force_b, & compute_tau=my_compute_tau, & - map_consistent=.TRUE., use_subpatch=use_subpatch, subpatch_pattern=tasks(itask)%subpatch_pattern) + use_subpatch=use_subpatch, subpatch_pattern=tasks(itask)%subpatch_pattern) ELSE rab_inv = -rab CALL integrate_pgf_product( & @@ -713,11 +721,11 @@ SUBROUTINE integrate_v_rspace_low(v_rspace, hmat, hmat_kp, pmat, pmat_kp, qs_env rb, rab_inv, rs_v(igrid_level)%rs_grid, cell, & cube_info(igrid_level), & hab, o1=nb1 - 1, o2=na1 - 1, & - eps_gvg_rspace=eps_rho_rspace, & + radius=radius, & calculate_forces=calculate_forces, & force_a=force_b, force_b=force_a, & compute_tau=my_compute_tau, & - map_consistent=.TRUE., use_subpatch=use_subpatch, subpatch_pattern=tasks(itask)%subpatch_pattern) + use_subpatch=use_subpatch, subpatch_pattern=tasks(itask)%subpatch_pattern) END IF END IF diff --git a/src/qs_integrate_potential_single.F b/src/qs_integrate_potential_single.F index 3f019e14da..0b37120f69 100644 --- a/src/qs_integrate_potential_single.F +++ b/src/qs_integrate_potential_single.F @@ -26,6 +26,7 @@ !> \author Matthias Krack (03.04.2001) ! ************************************************************************************************** MODULE qs_integrate_potential_single + USE ao_util, ONLY: exp_radius_very_extended USE atomic_kind_types, ONLY: atomic_kind_type,& get_atomic_kind USE atprop_types, ONLY: atprop_array_init,& @@ -106,7 +107,7 @@ SUBROUTINE integrate_ppl_rspace(rho_rspace, qs_env) n, natom_of_kind, ni, npme INTEGER, DIMENSION(:), POINTER :: atom_list, cores LOGICAL :: use_virial - REAL(KIND=dp) :: alpha, eps_rho_rspace + REAL(KIND=dp) :: alpha, eps_rho_rspace, radius REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, ra REAL(KIND=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b REAL(KIND=dp), DIMENSION(:), POINTER :: cexp_ppl @@ -243,10 +244,16 @@ SUBROUTINE integrate_ppl_rspace(rho_rspace, qs_env) END IF ni = 2*lppl - 2 + radius = exp_radius_very_extended(la_min=0, la_max=ni, lb_min=0, lb_max=0, & + ra=ra, rb=ra, rp=ra, & + zetp=alpha, eps=eps_rho_rspace, & + pab=pab, o1=0, o2=0, & ! without map_consistent + prefactor=1.0_dp, cutoff=1.0_dp) + CALL integrate_pgf_product(ni, alpha, 0, & 0, 0.0_dp, 0, ra, (/0.0_dp, 0.0_dp, 0.0_dp/), & rs_v, cell, pw_env%cube_info(1), hab, pab=pab, o1=0, o2=0, & - eps_gvg_rspace=eps_rho_rspace, & + radius=radius, & calculate_forces=.TRUE., force_a=force_a, & force_b=force_b, use_virial=use_virial, my_virial_a=my_virial_a, & my_virial_b=my_virial_b, use_subpatch=.TRUE., subpatch_pattern=0_int_8) @@ -290,7 +297,7 @@ SUBROUTINE integrate_rho_nlcc(rho_rspace, qs_env) ni, npme, nthread INTEGER, DIMENSION(:), POINTER :: atom_list, cores, nct_nlcc LOGICAL :: nlcc, use_virial - REAL(KIND=dp) :: alpha, eps_rho_rspace + REAL(KIND=dp) :: alpha, eps_rho_rspace, radius REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, ra REAL(KIND=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b REAL(KIND=dp), DIMENSION(:), POINTER :: alpha_nlcc @@ -438,10 +445,16 @@ SUBROUTINE integrate_rho_nlcc(rho_rspace, qs_env) END IF ni = 2*nc - 2 + radius = exp_radius_very_extended(la_min=0, la_max=ni, lb_min=0, lb_max=0, & + ra=ra, rb=ra, rp=ra, & + zetp=1/(2*alpha**2), eps=eps_rho_rspace, & + pab=pab, o1=0, o2=0, & ! without map_consistent + prefactor=1.0_dp, cutoff=1.0_dp) + CALL integrate_pgf_product(ni, 1/(2*alpha**2), 0, & 0, 0.0_dp, 0, ra, (/0.0_dp, 0.0_dp, 0.0_dp/), & rs_v, cell, pw_env%cube_info(1), hab, pab=pab, o1=0, o2=0, & - eps_gvg_rspace=eps_rho_rspace, & + radius=radius, & calculate_forces=.TRUE., force_a=force_a, & force_b=force_b, use_virial=use_virial, my_virial_a=my_virial_a, & my_virial_b=my_virial_b, use_subpatch=.TRUE., subpatch_pattern=0_int_8) @@ -487,7 +500,7 @@ SUBROUTINE integrate_v_core_rspace(v_rspace, qs_env) INTEGER, DIMENSION(:), POINTER :: atom_list, cores LOGICAL :: paw_atom, skip_fcore, use_virial REAL(KIND=dp) :: alpha_core_charge, ccore_charge, & - eps_rho_rspace + eps_rho_rspace, radius REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, ra REAL(KIND=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b REAL(KIND=dp), DIMENSION(:, :), POINTER :: hab, pab @@ -590,10 +603,16 @@ SUBROUTINE integrate_v_core_rspace(v_rspace, qs_env) my_virial_b = 0.0_dp END IF + radius = exp_radius_very_extended(la_min=0, la_max=0, lb_min=0, lb_max=0, & + ra=ra, rb=ra, rp=ra, & + zetp=alpha_core_charge, eps=eps_rho_rspace, & + pab=pab, o1=0, o2=0, & ! without map_consistent + prefactor=1.0_dp, cutoff=1.0_dp) + CALL integrate_pgf_product(0, alpha_core_charge, 0, & 0, 0.0_dp, 0, ra, (/0.0_dp, 0.0_dp, 0.0_dp/), & rs_v, cell, pw_env%cube_info(1), hab, pab=pab, o1=0, o2=0, & - eps_gvg_rspace=eps_rho_rspace, & + radius=radius, & calculate_forces=.TRUE., force_a=force_a, & force_b=force_b, use_virial=use_virial, my_virial_a=my_virial_a, & my_virial_b=my_virial_b, use_subpatch=.TRUE., subpatch_pattern=0_int_8) @@ -654,7 +673,7 @@ SUBROUTINE integrate_v_rspace_one_center(v_rspace, qs_env, int_res, & INTEGER, DIMENSION(:, :), POINTER :: first_sgfa LOGICAL :: use_virial LOGICAL, ALLOCATABLE, DIMENSION(:) :: map_it - REAL(KIND=dp) :: eps_rho_rspace + REAL(KIND=dp) :: eps_rho_rspace, radius REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, ra REAL(KIND=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a @@ -803,6 +822,11 @@ SUBROUTINE integrate_v_rspace_one_center(v_rspace, qs_env, int_res, & igrid_level = gaussian_gridlevel(gridlevel_info, zeta(ipgf, iset)) rs_grid => rs_v(igrid_level)%rs_grid + radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), & + lb_min=0, lb_max=0, ra=ra, rb=ra, rp=ra, & + zetp=zeta(ipgf, iset), eps=eps_rho_rspace, & + prefactor=1.0_dp, cutoff=1.0_dp) + IF (map_it(ipgf)) THEN IF (.NOT. calculate_forces) THEN CALL integrate_pgf_product(la_max=la_max(iset), & @@ -811,9 +835,8 @@ SUBROUTINE integrate_v_rspace_one_center(v_rspace, qs_env, int_res, & ra=ra, rab=(/0.0_dp, 0.0_dp, 0.0_dp/), & rsgrid=rs_grid, cell=cell, & cube_info=pw_env%cube_info(igrid_level), & - hab=hab, o1=na1, o2=0, eps_gvg_rspace=eps_rho_rspace, & - calculate_forces=calculate_forces, & - map_consistent=.TRUE.) + hab=hab, o1=na1, o2=0, radius=radius, & + calculate_forces=calculate_forces) ELSE CALL integrate_pgf_product(la_max=la_max(iset), & zeta=zeta(ipgf, iset), la_min=la_min(iset), & @@ -821,12 +844,11 @@ SUBROUTINE integrate_v_rspace_one_center(v_rspace, qs_env, int_res, & ra=ra, rab=(/0.0_dp, 0.0_dp, 0.0_dp/), & rsgrid=rs_grid, cell=cell, & cube_info=pw_env%cube_info(igrid_level), & - hab=hab, pab=pab, o1=na1, o2=0, eps_gvg_rspace=eps_rho_rspace, & + hab=hab, pab=pab, o1=na1, o2=0, radius=radius, & calculate_forces=calculate_forces, & force_a=force_a, force_b=force_b, & use_virial=use_virial, & - my_virial_a=my_virial_a, my_virial_b=my_virial_b, & - map_consistent=.TRUE.) + my_virial_a=my_virial_a, my_virial_b=my_virial_b) ENDIF ENDIF ENDDO @@ -892,8 +914,8 @@ SUBROUTINE integrate_v_rspace_diagonal(v_rspace, ksmat, pmat, qs_env, calculate_ INTEGER, DIMENSION(:, :), POINTER :: first_sgfa LOGICAL :: found, use_virial LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: map_it2 - REAL(KIND=dp) :: eps_rho_rspace, zetp - REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, ra, rab + REAL(KIND=dp) :: eps_rho_rspace, radius, zetp + REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, ra REAL(KIND=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a REAL(KIND=dp), DIMENSION(:, :), POINTER :: h_block, hab, hmat, p_block, pab, pblk, & @@ -957,7 +979,6 @@ SUBROUTINE integrate_v_rspace_diagonal(v_rspace, ksmat, pmat, qs_env, calculate_ DO iatom = 1, natom_of_kind atom_a = atom_list(iatom) ra(:) = pbc(particle_set(atom_a)%r, cell) - rab(:) = 0.0_dp hmat = 0.0_dp IF (calculate_forces) THEN CALL dbcsr_get_block_p(matrix=pmat, row=atom_a, col=atom_a, BLOCK=p_block, found=found) @@ -1036,32 +1057,39 @@ SUBROUTINE integrate_v_rspace_diagonal(v_rspace, ksmat, pmat, qs_env, calculate_ DO jpgf = 1, npgfa(jset) nb1 = (jpgf - 1)*ncoset(la_max(jset)) nb2 = jpgf*ncoset(la_max(jset)) - zetp = zeta(ipgf, iset) + zeta(jpgf, jset) igrid_level = gaussian_gridlevel(gridlevel_info, zetp) rs_grid => rs_v(igrid_level)%rs_grid + + radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), & + lb_min=la_min(jset), lb_max=la_max(jset), & + ra=ra, rb=ra, rp=ra, & + zetp=zetp, eps=eps_rho_rspace, & + prefactor=1.0_dp, cutoff=1.0_dp) + IF (map_it2(ipgf, jpgf)) THEN IF (calculate_forces) THEN CALL integrate_pgf_product( & la_max(iset), zeta(ipgf, iset), la_min(iset), & la_max(jset), zeta(jpgf, jset), la_min(jset), & - ra, rab, rs_v(igrid_level)%rs_grid, cell, & - pw_env%cube_info(igrid_level), & - hab, pab=pab, o1=na1, o2=nb1, & - eps_gvg_rspace=eps_rho_rspace, & - calculate_forces=calculate_forces, & + ra=ra, rab=(/0.0_dp, 0.0_dp, 0.0_dp/), & + rsgrid=rs_v(igrid_level)%rs_grid, cell=cell, & + cube_info=pw_env%cube_info(igrid_level), & + hab=hab, pab=pab, o1=na1, o2=nb1, & + radius=radius, & + calculate_forces=.TRUE., & force_a=force_a, force_b=force_b, & - use_virial=use_virial, my_virial_a=my_virial_a, my_virial_b=my_virial_b, & - map_consistent=.TRUE.) + use_virial=use_virial, my_virial_a=my_virial_a, my_virial_b=my_virial_b) ELSE CALL integrate_pgf_product( & la_max(iset), zeta(ipgf, iset), la_min(iset), & la_max(jset), zeta(jpgf, jset), la_min(jset), & - ra, rab, rs_v(igrid_level)%rs_grid, cell, & - pw_env%cube_info(igrid_level), & - hab, o1=na1, o2=nb1, & - eps_gvg_rspace=eps_rho_rspace, & - calculate_forces=.FALSE., map_consistent=.TRUE.) + ra=ra, rab=(/0.0_dp, 0.0_dp, 0.0_dp/), & + rsgrid=rs_v(igrid_level)%rs_grid, cell=cell, & + cube_info=pw_env%cube_info(igrid_level), & + hab=hab, o1=na1, o2=nb1, & + radius=radius, & + calculate_forces=.FALSE.) ENDIF ENDIF ENDDO diff --git a/src/qs_rho0_ggrid.F b/src/qs_rho0_ggrid.F index d699e0aff0..75ba275e0e 100644 --- a/src/qs_rho0_ggrid.F +++ b/src/qs_rho0_ggrid.F @@ -581,10 +581,9 @@ SUBROUTINE integrate_vhg0_rspace(qs_env, v_rspace, para_env, calculate_forces, l l0_ikind, zet0, 0, 0, 0.0_dp, 0, & ra, (/0.0_dp, 0.0_dp, 0.0_dp/), rs_v, cell, & cube_info(igrid), hab, pab, o1=0, o2=0, & - eps_gvg_rspace=eps_rho_rspace, & + radius=rpgf0, & calculate_forces=calculate_forces, & hdab=hdab, hadb=hadb, & - collocate_rho0=.TRUE., rpgf0_s=rpgf0, & use_virial=use_virial, my_virial_a=my_virial_a, my_virial_b=my_virial_b, & a_hdab=a_hdab, use_subpatch=.TRUE., subpatch_pattern=0_int_8)