Skip to content

Commit

Permalink
GAPW excited-state forces
Browse files Browse the repository at this point in the history
  • Loading branch information
Anna Hehn committed Jul 16, 2023
1 parent b8a1131 commit f80ece7
Show file tree
Hide file tree
Showing 16 changed files with 737 additions and 181 deletions.
88 changes: 63 additions & 25 deletions src/hartree_local_methods.F
View file
Expand Up @@ -203,36 +203,41 @@ END SUBROUTINE calculate_Vh_1center
!> \param local_rho_set ...
!> \param para_env ...
!> \param tddft ...
!> \param local_rho_set_2nd ...
!> \param core_2nd ...
!> \par History
!> 05.2012 JGH refactoring
!> \author ??
! **************************************************************************************************
SUBROUTINE Vh_1c_gg_integrals(qs_env, energy_hartree_1c, ecoul_1c, local_rho_set, para_env, tddft)
SUBROUTINE Vh_1c_gg_integrals(qs_env, energy_hartree_1c, ecoul_1c, local_rho_set, para_env, tddft, local_rho_set_2nd, &
core_2nd)

TYPE(qs_environment_type), POINTER :: qs_env
REAL(kind=dp), INTENT(out) :: energy_hartree_1c
TYPE(ecoul_1center_type), DIMENSION(:), POINTER :: ecoul_1c
TYPE(local_rho_type), POINTER :: local_rho_set
TYPE(mp_para_env_type), POINTER :: para_env
LOGICAL, INTENT(IN) :: tddft
TYPE(local_rho_type), OPTIONAL, POINTER :: local_rho_set_2nd
LOGICAL, INTENT(IN), OPTIONAL :: core_2nd

CHARACTER(LEN=*), PARAMETER :: routineN = 'Vh_1c_gg_integrals'

INTEGER :: bo(2), handle, iat, iatom, ikind, ipgf1, is1, iset1, iso, l_ang, llmax, lmax0, &
lmax_0, m1, max_iso, max_iso_not0, max_s_harm, maxl, maxso, mepos, n1, nat, natom, &
nchan_0, nkind, nr, nset, nsotot, nspins, num_pe
lmax0_2nd, lmax_0, m1, max_iso, max_iso_not0, max_s_harm, maxl, maxso, mepos, n1, nat, &
natom, nchan_0, nkind, nr, nset, nsotot, nspins, num_pe
INTEGER, ALLOCATABLE, DIMENSION(:) :: cg_n_list
INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cg_list
INTEGER, DIMENSION(:), POINTER :: atom_list, lmax, lmin, npgf
LOGICAL :: atenergy, core_charge, my_periodic, &
paw_atom
LOGICAL :: atenergy, core_charge, l_2nd_local_rho, &
my_core_2nd, my_periodic, paw_atom
REAL(dp) :: back_ch, factor
REAL(dp), ALLOCATABLE, DIMENSION(:) :: gexp, sqrtwr
REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: aVh1b_00, aVh1b_hh, aVh1b_ss, g0_h_w
REAL(dp), DIMENSION(:), POINTER :: rrad_z, vrrad_z
REAL(dp), DIMENSION(:, :), POINTER :: g0_h, gsph, rrad_0, Vh1_h, Vh1_s, &
vrrad_0, zet
REAL(dp), DIMENSION(:, :, :), POINTER :: my_CG, Qlm_gg
REAL(dp), DIMENSION(:, :), POINTER :: g0_h, g0_h_2nd, gsph, rrad_0, Vh1_h, &
Vh1_s, vrrad_0, zet
REAL(dp), DIMENSION(:, :, :), POINTER :: my_CG, Qlm_gg, Qlm_gg_2nd
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(atprop_type), POINTER :: atprop
TYPE(cell_type), POINTER :: cell
Expand All @@ -244,11 +249,11 @@ SUBROUTINE Vh_1c_gg_integrals(qs_env, energy_hartree_1c, ecoul_1c, local_rho_set
TYPE(pw_poisson_type), POINTER :: poisson_env
TYPE(qs_charges_type), POINTER :: qs_charges
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(rho0_atom_type), DIMENSION(:), POINTER :: rho0_atom_set
TYPE(rho0_mpole_type), POINTER :: rho0_mpole
TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
TYPE(rho0_atom_type), DIMENSION(:), POINTER :: rho0_atom_set, rho0_atom_set_2nd
TYPE(rho0_mpole_type), POINTER :: rho0_mpole, rho0_mpole_2nd
TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set, rho_atom_set_2nd
TYPE(rho_atom_type), POINTER :: rho_atom
TYPE(rhoz_type), DIMENSION(:), POINTER :: rhoz_set
TYPE(rhoz_type), DIMENSION(:), POINTER :: rhoz_set, rhoz_set_2nd

CALL timeset(routineN, handle)

Expand All @@ -271,12 +276,22 @@ SUBROUTINE Vh_1c_gg_integrals(qs_env, energy_hartree_1c, ecoul_1c, local_rho_set
back_ch = qs_charges%background*cell%deth

! rhoz_set is not accessed in TDDFT
CALL get_local_rho(local_rho_set, rho_atom_set, rho0_atom_set, rho0_mpole, rhoz_set)
CALL get_local_rho(local_rho_set, rho_atom_set, rho0_atom_set, rho0_mpole, rhoz_set) ! for integral space

! for forces we need a second local_rho_set
l_2nd_local_rho = .FALSE.
IF (PRESENT(local_rho_set_2nd)) THEN ! for potential
l_2nd_local_rho = .TRUE.
NULLIFY (rho_atom_set_2nd, rho0_atom_set_2nd, rhoz_set_2nd) ! for potential
CALL get_local_rho(local_rho_set_2nd, rho_atom_set_2nd, rho0_atom_set_2nd, rho0_mpole_2nd, rhoz_set=rhoz_set_2nd)
END IF

nkind = SIZE(atomic_kind_set, 1)
nspins = dft_control%nspins

core_charge = .NOT. tddft
core_charge = .NOT. tddft ! for forces mixed version
my_core_2nd = .TRUE.
IF (PRESENT(core_2nd)) my_core_2nd = .NOT. core_2nd ! if my_core_2nd true, include core charge

! The aim of the following code was to return immediately if the subroutine
! was called for triplet excited states in spin-restricted case. This check
Expand Down Expand Up @@ -339,8 +354,14 @@ SUBROUTINE Vh_1c_gg_integrals(qs_env, energy_hartree_1c, ecoul_1c, local_rho_set
NULLIFY (Qlm_gg, g0_h)
CALL get_rho0_mpole(rho0_mpole=rho0_mpole, ikind=ikind, &
l0_ikind=lmax0, &
Qlm_gg=Qlm_gg, g0_h=g0_h)

Qlm_gg=Qlm_gg, g0_h=g0_h) ! Qlm_gg of density

IF (PRESENT(local_rho_set_2nd)) THEN ! for potential
NULLIFY (Qlm_gg_2nd, g0_h_2nd)
CALL get_rho0_mpole(rho0_mpole=rho0_mpole_2nd, ikind=ikind, &
l0_ikind=lmax0_2nd, &
Qlm_gg=Qlm_gg_2nd, g0_h=g0_h_2nd) ! Qlm_gg of density
END IF
nchan_0 = nsoset(lmax0)

IF (nchan_0 > max_iso_not0) CPABORT("channels for rho0 > # max of spherical harmonics")
Expand Down Expand Up @@ -384,29 +405,46 @@ SUBROUTINE Vh_1c_gg_integrals(qs_env, energy_hartree_1c, ecoul_1c, local_rho_set

NULLIFY (rrad_z, vrrad_z, rrad_0, vrrad_0)
IF (core_charge) THEN
rrad_z => rhoz_set(ikind)%r_coef
vrrad_z => rhoz_set(ikind)%vr_coef
rrad_z => rhoz_set(ikind)%r_coef ! for density
END IF
IF (my_core_2nd) THEN
IF (l_2nd_local_rho) THEN
vrrad_z => rhoz_set_2nd(ikind)%vr_coef ! for potential
ELSE
vrrad_z => rhoz_set(ikind)%vr_coef ! for potential
END IF
END IF
rrad_0 => rho0_atom_set(iatom)%rho0_rad_h%r_coef
rrad_0 => rho0_atom_set(iatom)%rho0_rad_h%r_coef ! for density
vrrad_0 => rho0_atom_set(iatom)%vrho0_rad_h%r_coef
IF (l_2nd_local_rho) THEN
rho_atom => rho_atom_set_2nd(iatom)
vrrad_0 => rho0_atom_set_2nd(iatom)%vrho0_rad_h%r_coef ! for potential
END IF
IF (my_periodic .AND. back_ch .GT. 1.E-3_dp) THEN
factor = -2.0_dp*pi/3.0_dp*SQRT(fourpi)*qs_charges%background
ELSE
factor = 0._dp
END IF

CALL Vh_1c_atom_potential(rho_atom, vrrad_0, &
grid_atom, core_charge, vrrad_z, Vh1_h, Vh1_s, &
grid_atom, my_core_2nd, vrrad_z, Vh1_h, Vh1_s, & ! core charge for potential (2nd)
nchan_0, nspins, max_iso_not0, factor)

IF (l_2nd_local_rho) rho_atom => rho_atom_set(iatom) ! rho_atom for density

CALL Vh_1c_atom_energy(energy_hartree_1c, ecoul_1c, rho_atom, rrad_0, &
grid_atom, iatom, core_charge, rrad_z, Vh1_h, Vh1_s, &
grid_atom, iatom, core_charge, rrad_z, Vh1_h, Vh1_s, & ! core charge for density
nchan_0, nspins, max_iso_not0, atenergy, atprop%ate1c)

CALL Vh_1c_atom_integrals(rho_atom, &
IF (l_2nd_local_rho) rho_atom => rho_atom_set_2nd(iatom) ! rho_atom for potential (2nd)

CALL Vh_1c_atom_integrals(rho_atom, & ! results (int_local_h and int_local_s) written on rho_atom_2nd
! int_local_h and int_local_s are used in update_ks_atom
! on int_local_h mixed core / non-core
aVh1b_hh, aVh1b_ss, aVh1b_00, Vh1_h, Vh1_s, max_iso_not0, &
max_s_harm, llmax, cg_list, cg_n_list, &
nset, npgf, lmin, lmax, nsotot, maxso, nspins, nchan_0, gsph, g0_h_w, my_CG, Qlm_gg)
nset, npgf, lmin, lmax, nsotot, maxso, nspins, nchan_0, gsph, &
g0_h_w, my_CG, Qlm_gg) ! Qlm_gg for density from local_rho_set

END DO ! iat

Expand Down Expand Up @@ -607,7 +645,8 @@ END SUBROUTINE Vh_1c_atom_energy
SUBROUTINE Vh_1c_atom_integrals(rho_atom, &
aVh1b_hh, aVh1b_ss, aVh1b_00, Vh1_h, Vh1_s, max_iso_not0, &
max_s_harm, llmax, cg_list, cg_n_list, &
nset, npgf, lmin, lmax, nsotot, maxso, nspins, nchan_0, gsph, g0_h_w, my_CG, Qlm_gg)
nset, npgf, lmin, lmax, nsotot, maxso, nspins, nchan_0, gsph, &
g0_h_w, my_CG, Qlm_gg)

TYPE(rho_atom_type), POINTER :: rho_atom
REAL(dp), DIMENSION(:, :) :: aVh1b_hh, aVh1b_ss, aVh1b_00
Expand Down Expand Up @@ -702,7 +741,6 @@ SUBROUTINE Vh_1c_atom_integrals(rho_atom, &
END DO ! iset2
m1 = m1 + maxso
END DO !iset1

DO ispin = 1, nspins
CALL daxpy(nsotot*nsotot, 1.0_dp, aVh1b_hh, 1, int_local_h(ispin)%r_coef, 1)
CALL daxpy(nsotot*nsotot, 1.0_dp, aVh1b_ss, 1, int_local_s(ispin)%r_coef, 1)
Expand Down
2 changes: 1 addition & 1 deletion src/qs_kpp1_env_methods.F
View file
Expand Up @@ -527,7 +527,7 @@ SUBROUTINE calc_kpp1(rho1_xc, rho1, xc_section, do_tddft, lsd_singlets, lrigpw,
CALL Vh_1c_gg_integrals(qs_env, energy_hartree_1c, &
p_env%hartree_local%ecoul_1c, &
p_env%local_rho_set, &
para_env, tddft=.TRUE.)
para_env, tddft=.TRUE., core_2nd=.TRUE.)
CALL integrate_vhg0_rspace(qs_env, v_hartree_rspace, para_env, &
calculate_forces=my_calc_forces, &
local_rho_set=p_env%local_rho_set)
Expand Down
3 changes: 2 additions & 1 deletion src/qs_ks_methods.F
View file
Expand Up @@ -436,7 +436,8 @@ SUBROUTINE qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces, just_energy, &

IF (gapw) THEN
CALL get_qs_env(qs_env, ecoul_1c=ecoul_1c, local_rho_set=local_rho_set)
CALL Vh_1c_gg_integrals(qs_env, energy%hartree_1c, ecoul_1c, local_rho_set, para_env, tddft=.FALSE.)
CALL Vh_1c_gg_integrals(qs_env, energy%hartree_1c, ecoul_1c, local_rho_set, para_env, tddft=.FALSE., &
core_2nd=.FALSE.)
END IF

! Check if CDFT constraint is needed
Expand Down
3 changes: 2 additions & 1 deletion src/qs_ks_reference.F
View file
Expand Up @@ -363,7 +363,8 @@ SUBROUTINE ks_ref_potential_atom(qs_env, local_rho_set, local_rho_set_admm)

IF (gapw) THEN
CALL get_qs_env(qs_env, ecoul_1c=ecoul_1c)
CALL Vh_1c_gg_integrals(qs_env, eh1c, ecoul_1c, local_rho_set, para_env, tddft=.FALSE.)
CALL Vh_1c_gg_integrals(qs_env, eh1c, ecoul_1c, local_rho_set, para_env, tddft=.FALSE., &
core_2nd=.FALSE.)
END IF
IF (dft_control%do_admm) THEN
CALL get_qs_env(qs_env, admm_env=admm_env)
Expand Down
2 changes: 1 addition & 1 deletion src/qs_linres_kernel.F
View file
Expand Up @@ -495,7 +495,7 @@ SUBROUTINE apply_op_2_dft(qs_env, p_env)
CALL Vh_1c_gg_integrals(qs_env, energy_hartree_1c, &
p_env%hartree_local%ecoul_1c, &
p_env%local_rho_set, &
para_env, tddft=.TRUE.)
para_env, tddft=.TRUE., core_2nd=.TRUE.)

CALL integrate_vhg0_rspace(qs_env, v_hartree_rspace, para_env, &
calculate_forces=.FALSE., &
Expand Down
26 changes: 17 additions & 9 deletions src/qs_rho0_ggrid.F
View file
Expand Up @@ -312,14 +312,16 @@ END SUBROUTINE rho0_s_grid_create
!> \param para_env ...
!> \param calculate_forces ...
!> \param local_rho_set ...
!> \param local_rho_set_2nd ...
! **************************************************************************************************
SUBROUTINE integrate_vhg0_rspace(qs_env, v_rspace, para_env, calculate_forces, local_rho_set)
SUBROUTINE integrate_vhg0_rspace(qs_env, v_rspace, para_env, calculate_forces, local_rho_set, &
local_rho_set_2nd)

TYPE(qs_environment_type), POINTER :: qs_env
TYPE(pw_type), INTENT(IN) :: v_rspace
TYPE(mp_para_env_type), POINTER :: para_env
LOGICAL, INTENT(IN) :: calculate_forces
TYPE(local_rho_type), OPTIONAL, POINTER :: local_rho_set
TYPE(local_rho_type), OPTIONAL, POINTER :: local_rho_set, local_rho_set_2nd

CHARACTER(LEN=*), PARAMETER :: routineN = 'integrate_vhg0_rspace'

Expand Down Expand Up @@ -368,9 +370,9 @@ SUBROUTINE integrate_vhg0_rspace(qs_env, v_rspace, para_env, calculate_forces, l
! attribute 'OPTIONAL' from the argument 'local_rho_set'.

! CPASSERT(.NOT. (calculate_forces .AND. PRESENT(local_rho_set)))
IF (calculate_forces .AND. PRESENT(local_rho_set)) THEN
CPWARN("Forces and External Density")
END IF
! IF (calculate_forces .AND. PRESENT(local_rho_set)) THEN
! CPWARN("Forces and External Density!")
! END IF

NULLIFY (atomic_kind_set, qs_kind_set, dft_control, particle_set)
NULLIFY (cell, force, pw_env, rho0_mpole, rho_atom_set)
Expand Down Expand Up @@ -405,7 +407,12 @@ SUBROUTINE integrate_vhg0_rspace(qs_env, v_rspace, para_env, calculate_forces, l

IF (PRESENT(local_rho_set)) &
CALL get_local_rho(local_rho_set, rho0_mpole=rho0_mpole, rho_atom_set=rho_atom_set)

! Q from rho0_mpole of local_rho_set
! for TDDFT forces we need mixed potential / integral space
! potential stored on local_rho_set_2nd
IF (PRESENT(local_rho_set_2nd)) THEN
CALL get_local_rho(local_rho_set_2nd, rho_atom_set=rho_atom_set)
END IF
CALL get_rho0_mpole(rho0_mpole=rho0_mpole, lmax_0=lmax0, &
zet0_h=zet0, igrid_zet0_s=igrid, &
norm_g0l_h=norm_l)
Expand Down Expand Up @@ -494,7 +501,7 @@ SUBROUTINE integrate_vhg0_rspace(qs_env, v_rspace, para_env, calculate_forces, l

NULLIFY (Qlm_gg, lmax, npgf)
CALL get_rho0_mpole(rho0_mpole=rho0_mpole, ikind=ikind, &
l0_ikind=l0_ikind, Qlm_gg=Qlm_gg, &
l0_ikind=l0_ikind, Qlm_gg=Qlm_gg, & ! Qs different
rpgf0_s=rpgf0)

CALL get_gto_basis_set(gto_basis_set=basis_1c_set, &
Expand Down Expand Up @@ -602,7 +609,7 @@ SUBROUTINE integrate_vhg0_rspace(qs_env, v_rspace, para_env, calculate_forces, l
ig1 = iso1 + n1*(ipgf1 - 1) + m1
ig2 = iso2 + n2*(ipgf2 - 1) + m2

intloc(ig1, ig2) = intloc(ig1, ig2) + Qlm_gg(ig1, ig2, iso)*hab_sph(iso)
intloc(ig1, ig2) = intloc(ig1, ig2) + Qlm_gg(ig1, ig2, iso)*hab_sph(iso) ! potential times Q

END DO ! icg
END DO ! iso
Expand Down Expand Up @@ -631,7 +638,7 @@ SUBROUTINE integrate_vhg0_rspace(qs_env, v_rspace, para_env, calculate_forces, l
IF (calculate_forces) THEN
force_tmp(1:3) = 0.0_dp
DO iso = 1, nsoset(l0_ikind)
force_tmp(1) = force_tmp(1) + Qlm(iso)*hdab_sph(1, iso)
force_tmp(1) = force_tmp(1) + Qlm(iso)*hdab_sph(1, iso) ! Q here is from local_rho_set
force_tmp(2) = force_tmp(2) + Qlm(iso)*hdab_sph(2, iso)
force_tmp(3) = force_tmp(3) + Qlm(iso)*hdab_sph(3, iso)
END DO
Expand All @@ -642,6 +649,7 @@ SUBROUTINE integrate_vhg0_rspace(qs_env, v_rspace, para_env, calculate_forces, l
DO iso = 1, nsoset(l0_ikind)
DO ii = 1, 3
DO i = 1, 3
! Q from local_rho_set
virial%pv_gapw(i, ii) = virial%pv_gapw(i, ii) + Qlm(iso)*a_hdab_sph(i, ii, iso)
virial%pv_virial(i, ii) = virial%pv_virial(i, ii) + Qlm(iso)*a_hdab_sph(i, ii, iso)
END DO
Expand Down
31 changes: 30 additions & 1 deletion src/qs_rho_types.F
View file
Expand Up @@ -24,6 +24,8 @@ MODULE qs_rho_types
kpoint_transitional_type,&
set_1d_pointer,&
set_2d_pointer
USE pw_pool_types, ONLY: pw_pool_give_back_pw,&
pw_pool_type
USE pw_types, ONLY: pw_release,&
pw_type
#include "./base/base_uses.f90"
Expand All @@ -35,7 +37,7 @@ MODULE qs_rho_types
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_rho_types'

PUBLIC :: qs_rho_p_type, qs_rho_type
PUBLIC :: qs_rho_release, &
PUBLIC :: qs_rho_release, qs_rho_clear_pwpool, &
qs_rho_get, qs_rho_set, qs_rho_clear, qs_rho_create, qs_rho_unset_rho_ao

! **************************************************************************************************
Expand Down Expand Up @@ -337,5 +339,32 @@ SUBROUTINE qs_rho_set(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rh
IF (PRESENT(complex_rho_ao)) rho_struct%complex_rho_ao = complex_rho_ao
END SUBROUTINE qs_rho_set
! **************************************************************************************************
!> \brief ...
!> \param rho_struct ...
!> \param auxbas_pw_pool ...
! **************************************************************************************************
SUBROUTINE qs_rho_clear_pwpool(rho_struct, auxbas_pw_pool)
TYPE(qs_rho_type), INTENT(INOUT) :: rho_struct
TYPE(pw_pool_type), INTENT(IN), POINTER :: auxbas_pw_pool
INTEGER :: i
IF (ASSOCIATED(rho_struct%rho_r)) THEN
DO i = 1, SIZE(rho_struct%rho_r)
CALL pw_pool_give_back_pw(auxbas_pw_pool, rho_struct%rho_r(i))
END DO
DEALLOCATE (rho_struct%rho_r)
NULLIFY (rho_struct%rho_r)
END IF
IF (ASSOCIATED(rho_struct%rho_g)) THEN
DO i = 1, SIZE(rho_struct%rho_g)
CALL pw_pool_give_back_pw(auxbas_pw_pool, rho_struct%rho_g(i))
END DO
DEALLOCATE (rho_struct%rho_g)
NULLIFY (rho_struct%rho_g)
END IF
END SUBROUTINE qs_rho_clear_pwpool
END MODULE qs_rho_types

0 comments on commit f80ece7

Please sign in to comment.