GAPW TDDFPT (cp2k#2932)

src/qs_environment.F

@@ -1146,7 +1146,7 @@ SUBROUTINE qs_init_subsys(qs_env, para_env, subsys, cell, cell_ref, use_ref_cell
          IF (dft_control%qs_control%method_id /= do_method_gapw_xc) THEN
             CALL get_qs_env(qs_env=qs_env, local_rho_set=local_rho_set, natom=natom)
             !       *** Allocate and initialize the compensation density rho0 ***
-            CALL init_rho0(local_rho_set, qs_env, gapw_control, .FALSE.)
+            CALL init_rho0(local_rho_set, qs_env, gapw_control)
             !       *** Allocate and Initialize the local coulomb term ***
             CALL init_coulomb_local(qs_env%hartree_local, natom)
          END IF

src/qs_fxc.F

@@ -203,7 +203,9 @@ SUBROUTINE qs_fxc_fdiff(ks_env, rho0_struct, rho1_struct, xc_section, accuracy,
             NULLIFY (vxc00, v_tau_rspace)
             IF (is_triplet) THEN
                CPASSERT(nspins == 1)
+               ! rhoin = (0.5 rho0, 0.5 rho0)
                CALL qs_rho_copy(rho0_struct, rhoin, auxbas_pw_pool, 2)
+               ! rhoin = (0.5 rho0 + 0.5 rho1, 0.5 rho0)
                CALL qs_rho_scale_and_add(rhoin, rho1_struct, alpha, 0.5_dp*beta)
                CALL qs_vxc_create(ks_env=ks_env, rho_struct=rhoin, xc_section=xc_section, &
                                   vxc_rho=vxc00, vxc_tau=v_tau_rspace, exc=exc, just_energy=.FALSE.)

src/qs_ks_atom.F

@@ -94,6 +94,8 @@ MODULE qs_ks_atom
 !> \param oce_external ...
 !> \param sab_external ...
 !> \param kscale ...
+!> \param kintegral ...
+!> \param kforce ...
 !> \param fscale ...
 !> \par History
 !>      created [MI]
@@ -103,7 +105,8 @@ MODULE qs_ks_atom
 !>      Allow for external kind_set, rho_atom_set, oce, sab 12.2019 (A. Bussy)
 ! **************************************************************************************************
    SUBROUTINE update_ks_atom(qs_env, ksmat, pmat, forces, tddft, rho_atom_external, &
-                             kind_set_external, oce_external, sab_external, kscale, fscale)
+                             kind_set_external, oce_external, sab_external, kscale, &
+                             kintegral, kforce, fscale)
 
       TYPE(qs_environment_type), POINTER                 :: qs_env
       TYPE(dbcsr_p_type), DIMENSION(*), INTENT(INOUT)    :: ksmat, pmat
@@ -116,7 +119,8 @@ SUBROUTINE update_ks_atom(qs_env, ksmat, pmat, forces, tddft, rho_atom_external,
       TYPE(oce_matrix_type), OPTIONAL, POINTER           :: oce_external
       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
          OPTIONAL, POINTER                               :: sab_external
-      REAL(KIND=dp), INTENT(IN), OPTIONAL                :: kscale, fscale(2)
+      REAL(KIND=dp), INTENT(IN), OPTIONAL                :: kscale, kintegral, kforce
+      REAL(KIND=dp), DIMENSION(2), INTENT(IN), OPTIONAL  :: fscale
 
       CHARACTER(len=*), PARAMETER                        :: routineN = 'update_ks_atom'
 
@@ -134,10 +138,11 @@ SUBROUTINE update_ks_atom(qs_env, ksmat, pmat, forces, tddft, rho_atom_external,
       REAL(dp), ALLOCATABLE, DIMENSION(:, :)             :: a_matrix, p_matrix
       REAL(dp), DIMENSION(3)                             :: rac, rbc
       REAL(dp), DIMENSION(3, 3)                          :: force_tmp
-      REAL(kind=dp)                                      :: eps_cpc, factor1, factor2, force_fac(2)
+      REAL(kind=dp)                                      :: eps_cpc, factor1, factor2
       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: C_int_h, C_int_s, coc
       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: dCPC_h, dCPC_s
       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: PC_h, PC_s
+      REAL(KIND=dp), DIMENSION(2)                        :: force_fac
       REAL(KIND=dp), DIMENSION(3, 3)                     :: pv_virial_thread
       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: C_coeff_hh_a, C_coeff_hh_b, &
                                                             C_coeff_ss_a, C_coeff_ss_b
@@ -185,11 +190,12 @@ SUBROUTINE update_ks_atom(qs_env, ksmat, pmat, forces, tddft, rho_atom_external,
       nspins = dft_control%nspins
       nimages = dft_control%nimages
 
+      factor1 = 1.0_dp
+      factor2 = 1.0_dp
+
       !deal with externals
       my_tddft = .FALSE.
       IF (PRESENT(tddft)) my_tddft = tddft
-      factor1 = 1.0_dp
-      factor2 = 1.0_dp
       IF (my_tddft) THEN
          IF (nspins == 1) factor1 = 2.0_dp
          CPASSERT(nimages == 1)
@@ -198,7 +204,8 @@ SUBROUTINE update_ks_atom(qs_env, ksmat, pmat, forces, tddft, rho_atom_external,
          factor1 = factor1*kscale
          factor2 = factor2*kscale
       END IF
-
+      IF (PRESENT(kintegral)) factor1 = kintegral
+      IF (PRESENT(kforce)) factor2 = kforce
       force_fac = 1.0_dp
       IF (PRESENT(fscale)) force_fac(:) = fscale(:)
 
@@ -376,6 +383,7 @@ SUBROUTINE update_ks_atom(qs_env, ksmat, pmat, forces, tddft, rho_atom_external,
                DO ispin = 1, nspins
                   NULLIFY (mat_h(ispin)%array, mat_p(ispin)%array)
 
+                  found = .FALSE.
                   IF (iatom <= jatom) THEN
                      CALL dbcsr_get_block_p(matrix=ksmat(nspins*(img - 1) + ispin)%matrix, &
                                             row=iatom, col=jatom, &
@@ -385,8 +393,10 @@ SUBROUTINE update_ks_atom(qs_env, ksmat, pmat, forces, tddft, rho_atom_external,
                                             row=jatom, col=iatom, &
                                             BLOCK=mat_h(ispin)%array, found=found)
                   END IF
+                  CPASSERT(found)
 
                   IF (forces) THEN
+                     found = .FALSE.
                      IF (iatom <= jatom) THEN
                         CALL dbcsr_get_block_p(matrix=pmat(nspins*(img - 1) + ispin)%matrix, &
                                                row=iatom, col=jatom, &
@@ -396,6 +406,7 @@ SUBROUTINE update_ks_atom(qs_env, ksmat, pmat, forces, tddft, rho_atom_external,
                                                row=jatom, col=iatom, &
                                                BLOCK=mat_p(ispin)%array, found=found)
                      END IF
+                     CPASSERT(found)
                   END IF
                END DO
 

src/qs_ks_reference.F

@@ -19,8 +19,7 @@ MODULE qs_ks_reference
    USE dbcsr_api,                       ONLY: dbcsr_p_type
    USE hartree_local_methods,           ONLY: Vh_1c_gg_integrals,&
                                               init_coulomb_local
-   USE hartree_local_types,             ONLY: ecoul_1center_type,&
-                                              hartree_local_create,&
+   USE hartree_local_types,             ONLY: hartree_local_create,&
                                               hartree_local_release,&
                                               hartree_local_type
    USE input_constants,                 ONLY: do_admm_aux_exch_func_none
@@ -58,7 +57,8 @@ MODULE qs_ks_reference
                                               local_rho_type
    USE qs_neighbor_list_types,          ONLY: neighbor_list_set_p_type
    USE qs_oce_types,                    ONLY: oce_matrix_type
-   USE qs_rho0_ggrid,                   ONLY: rho0_s_grid_create
+   USE qs_rho0_ggrid,                   ONLY: integrate_vhg0_rspace,&
+                                              rho0_s_grid_create
    USE qs_rho0_methods,                 ONLY: init_rho0
    USE qs_rho_atom_methods,             ONLY: allocate_rho_atom_internals,&
                                               calculate_rho_atom_coeff
@@ -301,13 +301,15 @@ END SUBROUTINE ks_ref_potential
 !> \param qs_env ...
 !> \param local_rho_set ...
 !> \param local_rho_set_admm ...
+!> \param v_hartree_rspace ...
 !> \par History
 !>      07.2022 created [JGH]
 !> \author JGH
 ! **************************************************************************************************
-   SUBROUTINE ks_ref_potential_atom(qs_env, local_rho_set, local_rho_set_admm)
+   SUBROUTINE ks_ref_potential_atom(qs_env, local_rho_set, local_rho_set_admm, v_hartree_rspace)
       TYPE(qs_environment_type), POINTER                 :: qs_env
       TYPE(local_rho_type), POINTER                      :: local_rho_set, local_rho_set_admm
+      TYPE(pw_type), INTENT(IN)                          :: v_hartree_rspace
 
       CHARACTER(LEN=*), PARAMETER :: routineN = 'ks_ref_potential_atom'
 
@@ -318,7 +320,6 @@ SUBROUTINE ks_ref_potential_atom(qs_env, local_rho_set, local_rho_set_admm)
       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: rho_ao_aux, rho_ao_kp
       TYPE(dft_control_type), POINTER                    :: dft_control
-      TYPE(ecoul_1center_type), DIMENSION(:), POINTER    :: ecoul_1c
       TYPE(hartree_local_type), POINTER                  :: hartree_local
       TYPE(mp_para_env_type), POINTER                    :: para_env
       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
@@ -349,7 +350,7 @@ SUBROUTINE ks_ref_potential_atom(qs_env, local_rho_set, local_rho_set_admm)
                                           qs_kind_set, dft_control, para_env)
          IF (gapw) THEN
             CALL get_qs_env(qs_env, natom=natom)
-            CALL init_rho0(local_rho_set, qs_env, dft_control%qs_control%gapw_control, .FALSE.)
+            CALL init_rho0(local_rho_set, qs_env, dft_control%qs_control%gapw_control)
             CALL rho0_s_grid_create(pw_env, local_rho_set%rho0_mpole)
             CALL hartree_local_create(hartree_local)
             CALL init_coulomb_local(hartree_local, natom)
@@ -362,9 +363,9 @@ SUBROUTINE ks_ref_potential_atom(qs_env, local_rho_set, local_rho_set_admm)
          CALL prepare_gapw_den(qs_env, local_rho_set, do_rho0=gapw)
 
          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., &
-                                    core_2nd=.FALSE.)
+            CALL Vh_1c_gg_integrals(qs_env, eh1c, hartree_local%ecoul_1c, local_rho_set, para_env, .FALSE.)
+            CALL integrate_vhg0_rspace(qs_env, v_hartree_rspace, para_env, calculate_forces=.FALSE., &
+                                       local_rho_set=local_rho_set)
          END IF
          IF (dft_control%do_admm) THEN
             CALL get_qs_env(qs_env, admm_env=admm_env)
@@ -386,7 +387,7 @@ SUBROUTINE ks_ref_potential_atom(qs_env, local_rho_set, local_rho_set_admm)
                CALL qs_rho_get(rho, rho_ao_kp=rho_ao_aux)
                CALL calculate_rho_atom_coeff(qs_env, rho_ao_aux, local_rho_set_admm%rho_atom_set, &
                                              admm_env%admm_gapw_env%admm_kind_set, oce, sab, para_env)
-               CALL prepare_gapw_den(qs_env, local_rho_set=admm_env%admm_gapw_env%local_rho_set, &
+               CALL prepare_gapw_den(qs_env, local_rho_set=local_rho_set_admm, &
                                      do_rho0=.FALSE., kind_set_external=admm_env%admm_gapw_env%admm_kind_set)
                !compute the potential due to atomic densities
                xc_section => admm_env%xc_section_aux

src/qs_p_env_methods.F

@@ -262,7 +262,8 @@ SUBROUTINE p_env_create(p_env, qs_env, p1_option, p1_admm_option, &
          CALL allocate_rho_atom_internals(p_env%local_rho_set%rho_atom_set, atomic_kind_set, &
                                           qs_kind_set, dft_control, para_env)
 
-         CALL init_rho0(p_env%local_rho_set, qs_env, dft_control%qs_control%gapw_control, .TRUE.)
+         CALL init_rho0(p_env%local_rho_set, qs_env, dft_control%qs_control%gapw_control, &
+                        zcore=0.0_dp)
          CALL rho0_s_grid_create(pw_env, p_env%local_rho_set%rho0_mpole)
          CALL hartree_local_create(p_env%hartree_local)
          CALL init_coulomb_local(p_env%hartree_local, natom)

src/qs_rho0_ggrid.F

@@ -319,16 +319,18 @@ END SUBROUTINE rho0_s_grid_create
 !> \param local_rho_set ...
 !> \param local_rho_set_2nd ...
 !> \param atener ...
+!> \param kforce ...
 ! **************************************************************************************************
    SUBROUTINE integrate_vhg0_rspace(qs_env, v_rspace, para_env, calculate_forces, local_rho_set, &
-                                    local_rho_set_2nd, atener)
+                                    local_rho_set_2nd, atener, kforce)
 
       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, local_rho_set_2nd
       REAL(KIND=dp), DIMENSION(:), OPTIONAL              :: atener
+      REAL(KIND=dp), INTENT(IN), OPTIONAL                :: kforce
 
       CHARACTER(LEN=*), PARAMETER :: routineN = 'integrate_vhg0_rspace'
 
@@ -340,8 +342,8 @@ SUBROUTINE integrate_vhg0_rspace(qs_env, v_rspace, para_env, calculate_forces, l
       INTEGER, ALLOCATABLE, DIMENSION(:, :, :)           :: cg_list
       INTEGER, DIMENSION(:), POINTER                     :: atom_list, lmax, lmin, npgf
       LOGICAL                                            :: grid_distributed, paw_atom, use_virial
-      REAL(KIND=dp)                                      :: eps_rho_rspace, force_tmp(3), ra(3), &
-                                                            rpgf0, scale, zet0
+      REAL(KIND=dp)                                      :: eps_rho_rspace, force_tmp(3), fscale, &
+                                                            ra(3), rpgf0, zet0
       REAL(KIND=dp), DIMENSION(3, 3)                     :: my_virial_a, my_virial_b
       REAL(KIND=dp), DIMENSION(:), POINTER               :: hab_sph, norm_l, Qlm
       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: hab, hdab_sph, intloc, pab
@@ -455,8 +457,8 @@ SUBROUTINE integrate_vhg0_rspace(qs_env, v_rspace, para_env, calculate_forces, l
          CALL pw_pool_give_back_pw(pw_aux, coeff_gaux)
          CALL pw_pool_create_pw(pw_aux, coeff_raux, use_data=REALDATA3D, &
                                 in_space=REALSPACE)
-         scale = coeff_rspace%pw_grid%dvol/coeff_raux%pw_grid%dvol
-         coeff_rspace%cr3d = scale*coeff_rspace%cr3d
+         fscale = coeff_rspace%pw_grid%dvol/coeff_raux%pw_grid%dvol
+         coeff_rspace%cr3d = fscale*coeff_rspace%cr3d
          CALL pw_pool_give_back_pw(pw_aux, coeff_raux)
       ELSE
 
@@ -496,6 +498,11 @@ SUBROUTINE integrate_vhg0_rspace(qs_env, v_rspace, para_env, calculate_forces, l
 
       grid_distributed = rs_v%desc%distributed
 
+      fscale = 1.0_dp
+      IF (PRESENT(kforce)) THEN
+         fscale = kforce
+      END IF
+
       DO ikind = 1, SIZE(atomic_kind_set, 1)
          NULLIFY (basis_1c_set, atom_list, harmonics)
          CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=nat)
@@ -655,16 +662,16 @@ SUBROUTINE integrate_vhg0_rspace(qs_env, v_rspace, para_env, calculate_forces, l
                      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
-                  force(ikind)%g0s_Vh_elec(1:3, iat) = force(ikind)%g0s_Vh_elec(1:3, iat) + force_tmp(1:3)
+                  force(ikind)%g0s_Vh_elec(1:3, iat) = force(ikind)%g0s_Vh_elec(1:3, iat) + fscale*force_tmp(1:3)
                END IF
                IF (use_virial) THEN
                   my_virial_a = 0.0_dp
                   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)
+                        virial%pv_gapw(i, ii) = virial%pv_gapw(i, ii) + fscale*Qlm(iso)*a_hdab_sph(i, ii, iso)
+                        virial%pv_virial(i, ii) = virial%pv_virial(i, ii) + fscale*Qlm(iso)*a_hdab_sph(i, ii, iso)
                      END DO
                      END DO
                   END DO

src/qs_rho0_methods.F

@@ -342,14 +342,14 @@ END SUBROUTINE calculate_rho0_atom
 !> \param local_rho_set ...
 !> \param qs_env ...
 !> \param gapw_control ...
-!> \param tddft ...
+!> \param zcore ...
 ! **************************************************************************************************
-   SUBROUTINE init_rho0(local_rho_set, qs_env, gapw_control, tddft)
+   SUBROUTINE init_rho0(local_rho_set, qs_env, gapw_control, zcore)
 
       TYPE(local_rho_type), POINTER                      :: local_rho_set
       TYPE(qs_environment_type), POINTER                 :: qs_env
       TYPE(gapw_control_type), POINTER                   :: gapw_control
-      LOGICAL, INTENT(in)                                :: tddft
+      REAL(KIND=dp), INTENT(IN), OPTIONAL                :: zcore
 
       CHARACTER(len=*), PARAMETER                        :: routineN = 'init_rho0'
 
@@ -430,7 +430,7 @@ SUBROUTINE init_rho0(local_rho_set, qs_env, gapw_control, tddft)
                           basis_set=basis_1c, basis_type="GAPW_1C")
 
          ! Set charge distribution of ionic cores to zero when computing the response-density
-         IF (tddft) zeff = 0.0_dp
+         IF (PRESENT(zcore)) zeff = zcore
 
          CALL get_gto_basis_set(gto_basis_set=basis_1c, &
                                 maxl=maxl, &