Skip to content

Commit

Permalink
Bug fix: GAPW TDDFPT Triplets with ADMM/XC-correction kernel (#3018)
Browse files Browse the repository at this point in the history
  • Loading branch information
@juerghutter
juerghutter committed Oct 4, 2023
1 parent cc008e2 commit 21eaee0
Show file tree
Hide file tree
Showing 4 changed files with 39 additions and 25 deletions.
10 changes: 6 additions & 4 deletions src/cp2k_debug.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -90,8 +90,8 @@ SUBROUTINE cp2k_debug_energy_and_forces(force_env)
LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: do_dof_atom_coor
LOGICAL, DIMENSION(3) :: do_dof_dipole
REAL(KIND=dp) :: amplitude, dd, de, derr, difference, dx, &
eps_no_error_check, maxerr, std_value, &
sum_of_differences
eps_no_error_check, errmax, maxerr, &
std_value, sum_of_differences
REAL(KIND=dp), DIMENSION(2) :: numer_energy
REAL(KIND=dp), DIMENSION(3) :: dipole_moment, dipole_numer, err, &
my_maxerr, poldir
Expand Down Expand Up @@ -384,13 +384,15 @@ SUBROUTINE cp2k_debug_energy_and_forces(force_env)
"DEBUG|======================== BEGIN OF SUMMARY ===============================", &
"DEBUG| Atom Coordinate f(numerical) f(analytical) Difference Error [%]"
sum_of_differences = 0.0_dp
errmax = 0.0_dp
DO ip = 1, np
err(1:3) = 0.0_dp
DO k = 1, 3
IF (do_dof_atom_coor(k, ip)) THEN
difference = analyt_forces(ip, k) - numer_forces(ip, k)
IF (ABS(analyt_forces(ip, k)) >= eps_no_error_check) THEN
err(k) = 100_dp*(numer_forces(ip, k) - analyt_forces(ip, k))/analyt_forces(ip, k)
errmax = MAX(errmax, ABS(err(k)))
WRITE (UNIT=iw, FMT="(T2,A,I5,T19,A1,T26,F14.8,T42,F14.8,T57,F12.8,T71,F10.2)") &
"DEBUG|", ip, ACHAR(119 + k), numer_forces(ip, k), analyt_forces(ip, k), difference, err(k)
ELSE
Expand All @@ -401,8 +403,8 @@ SUBROUTINE cp2k_debug_energy_and_forces(force_env)
END IF
END DO
END DO
WRITE (UNIT=iw, FMT="(T2,A,T57,F12.8)") &
"DEBUG| Sum of differences:", sum_of_differences
WRITE (UNIT=iw, FMT="(T2,A,T57,F12.8,T71,F10.2)") &
"DEBUG| Sum of differences:", sum_of_differences, errmax
WRITE (UNIT=iw, FMT="(T2,A)") &
"DEBUG|======================== END OF SUMMARY ================================="
END IF
Expand Down
10 changes: 5 additions & 5 deletions src/qs_rho_types.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -64,13 +64,13 @@ MODULE qs_rho_types
PRIVATE
TYPE(kpoint_transitional_type) :: rho_ao
TYPE(kpoint_transitional_type) :: rho_ao_im
TYPE(pw_type), DIMENSION(:), POINTER :: rho_r => Null(), &
rho_g => Null(), &
tau_r => Null(), &
tau_g => Null()
TYPE(pw_type), DIMENSION(:), POINTER :: rho_r => Null(), &
rho_g => Null(), &
tau_r => Null(), &
tau_g => Null()
TYPE(pw_type), DIMENSION(:, :), POINTER :: drho_r => NULL(), drho_g => NULL()
! Final rho_iter of last SCCS cycle (r-space)
TYPE(pw_type), POINTER :: rho_r_sccs => Null()
TYPE(pw_type), POINTER :: rho_r_sccs => Null()
LOGICAL :: rho_g_valid = .FALSE., &
rho_r_valid = .FALSE., &
drho_r_valid = .FALSE., &
Expand Down
3 changes: 2 additions & 1 deletion src/qs_tddfpt2_fhxc.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -328,7 +328,8 @@ SUBROUTINE fhxc_kernel(Aop_evects, evects, is_rks_triplets, &
rho1_atom_set => work_matrices%local_rho_set_admm%rho_atom_set
CALL calculate_xc_2nd_deriv_atom(rho_atom_set, rho1_atom_set, qs_env, &
kernel_env_admm_aux%xc_section, &
sub_env%para_env, do_tddfpt2=.TRUE., do_triplet=.FALSE., &
sub_env%para_env, do_tddfpt2=.TRUE., &
do_triplet=is_rks_triplets, &
kind_set_external=admm_env%admm_gapw_env%admm_kind_set)
CALL update_ks_atom(qs_env, A_xc_munu_sub, rho_ia_ao_aux_fit, forces=.FALSE., tddft=.TRUE., &
rho_atom_external=rho1_atom_set, &
Expand Down
41 changes: 26 additions & 15 deletions src/response_solver.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1740,35 +1740,46 @@ SUBROUTINE response_force(qs_env, vh_rspace, vxc_rspace, vtau_rspace, vadmm_rspa
CALL get_qs_env(qs_env=qs_env, rho=rho)
CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
IF (debug_stress .AND. use_virial) stdeb = virial%pv_virial
IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
DO ispin = 1, nspins
CALL pw_scale(v_xc(ispin), v_xc(ispin)%pw_grid%dvol)
IF ((.NOT. (gapw)) .AND. (.NOT. gapw_xc)) THEN
END DO
IF ((.NOT. (gapw)) .AND. (.NOT. gapw_xc)) THEN
IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
DO ispin = 1, nspins
CALL pw_axpy(zv_hartree_rspace, v_xc(ispin)) ! Hartree potential of response density
CALL integrate_v_rspace(qs_env=qs_env, &
v_rspace=v_xc(ispin), &
hmat=matrix_hz(ispin), &
pmat=matrix_p(ispin, 1), &
gapw=.FALSE., &
calculate_forces=.TRUE.)
ELSE
IF (myfun /= xc_none) THEN
IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
END DO
IF (debug_forces) THEN
fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
CALL para_env%sum(fodeb)
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pin*dKh*rhoz ", fodeb
END IF
ELSE
IF (myfun /= xc_none) THEN
IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
DO ispin = 1, nspins
CALL integrate_v_rspace(qs_env=qs_env, &
v_rspace=v_xc(ispin), &
hmat=matrix_hz(ispin), &
pmat=matrix_p(ispin, 1), &
gapw=.TRUE., &
calculate_forces=.TRUE.)
IF (debug_forces) THEN
fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
CALL para_env%sum(fodeb)
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pin*dKxc*rhoz ", fodeb
END IF
END IF ! my_fun
! Coulomb T+Dz
IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
END DO
IF (debug_forces) THEN
fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
CALL para_env%sum(fodeb)
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pin*dKxc*rhoz ", fodeb
END IF
END IF ! my_fun
! Coulomb T+Dz
IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
DO ispin = 1, nspins
CALL pw_zero(v_xc(ispin))
IF (gapw) THEN ! Hartree potential of response density
CALL pw_axpy(v_hartree_rspace_t, v_xc(ispin))
Expand All @@ -1781,13 +1792,13 @@ SUBROUTINE response_force(qs_env, vh_rspace, vxc_rspace, vtau_rspace, vadmm_rspa
pmat=matrix_p(ispin, 1), &
gapw=gapw, &
calculate_forces=.TRUE.)
END IF
END DO
IF (debug_forces) THEN
fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
CALL para_env%sum(fodeb)
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pin*dKh*rhoz ", fodeb
END IF
END DO
END IF
IF (gapw .OR. gapw_xc) THEN
! compute hard and soft atomic contributions
Expand Down

0 comments on commit 21eaee0

Please sign in to comment.