Bug fixes for TDDFPT-Forces + kernel ADMM + X correction (#3141)

cp2k · Dec 1, 2023 · 05a2402 · 05a2402
1 parent ddec052
commit 05a2402
Show file tree

Hide file tree

Showing 6 changed files with 236 additions and 27 deletions.
diff --git a/src/qs_tddfpt2_fhxc_forces.F b/src/qs_tddfpt2_fhxc_forces.F
@@ -791,17 +791,17 @@ SUBROUTINE fhxc_force(qs_env, ex_env, gs_mos, full_kernel, debug_forces)
 
                IF (debug_forces) fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1)
                IF (nspins == 1) THEN
-                  CALL update_ks_atom(qs_env, mfx, matrix_px1_admm, forces=.TRUE., tddft=.TRUE., &
+                  CALL update_ks_atom(qs_env, mfx, matrix_px1_admm, forces=.TRUE., &
                                       rho_atom_external=rho_atom_set_f, &
                                       kind_set_external=admm_env%admm_gapw_env%admm_kind_set, &
                                       oce_external=admm_env%admm_gapw_env%oce, sab_external=sab_aux_fit, &
-                                      kintegral=1.0_dp, kforce=0.25_dp)
+                                      kintegral=2.0_dp, kforce=0.5_dp)
                ELSE
-                  CALL update_ks_atom(qs_env, mfx, matrix_px1_admm, forces=.TRUE., tddft=.TRUE., &
+                  CALL update_ks_atom(qs_env, mfx, matrix_px1_admm, forces=.TRUE., &
                                       rho_atom_external=rho_atom_set_f, &
                                       kind_set_external=admm_env%admm_gapw_env%admm_kind_set, &
                                       oce_external=admm_env%admm_gapw_env%oce, sab_external=sab_aux_fit, &
-                                      kintegral=1.0_dp, kforce=0.5_dp)
+                                      kintegral=2.0_dp, kforce=1.0_dp)
                END IF
                IF (debug_forces) THEN
                   fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1) - fodeb(1:3)
@@ -810,17 +810,17 @@ SUBROUTINE fhxc_force(qs_env, ex_env, gs_mos, full_kernel, debug_forces)
                END IF
                IF (debug_forces) fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1)
                IF (nspins == 1) THEN
-                  CALL update_ks_atom(qs_env, mgx, matrix_p_admm, forces=.TRUE., tddft=.TRUE., &
+                  CALL update_ks_atom(qs_env, mgx, matrix_p, forces=.TRUE., &
                                       rho_atom_external=rho_atom_set_g, &
                                       kind_set_external=admm_env%admm_gapw_env%admm_kind_set, &
                                       oce_external=admm_env%admm_gapw_env%oce, sab_external=sab_aux_fit, &
-                                      kintegral=0.5_dp, kforce=0.25_dp)
+                                      kintegral=1.0_dp, kforce=0.5_dp)
                ELSE
-                  CALL update_ks_atom(qs_env, mgx, matrix_p_admm, forces=.TRUE., &
+                  CALL update_ks_atom(qs_env, mgx, matrix_p, forces=.TRUE., &
                                       rho_atom_external=rho_atom_set_g, &
                                       kind_set_external=admm_env%admm_gapw_env%admm_kind_set, &
                                       oce_external=admm_env%admm_gapw_env%oce, sab_external=sab_aux_fit, &
-                                      kintegral=0.5_dp, kforce=0.5_dp)
+                                      kintegral=1.0_dp, kforce=1.0_dp)
                END IF
                IF (debug_forces) THEN
                   fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1) - fodeb(1:3)

diff --git a/src/response_solver.F b/src/response_solver.F
@@ -2052,7 +2052,7 @@ SUBROUTINE response_force(qs_env, vh_rspace, vxc_rspace, vtau_rspace, vadmm_rspa
                CALL get_admm_env(admm_env, sab_aux_fit=sab_aux_fit)
                IF (debug_forces) fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1)
                CALL update_ks_atom(qs_env, mhx(:, 1), mpz(:, 1), forces=.TRUE., tddft=.FALSE., &
-                                   rho_atom_external=ex_env%local_rho_set_admm%rho_atom_set, &
+                                   rho_atom_external=admm_env%admm_gapw_env%local_rho_set%rho_atom_set, &
                                    kind_set_external=admm_env%admm_gapw_env%admm_kind_set, &
                                    oce_external=admm_env%admm_gapw_env%oce, &
                                    sab_external=sab_aux_fit)
@@ -2077,8 +2077,7 @@ SUBROUTINE response_force(qs_env, vh_rspace, vxc_rspace, vtau_rspace, vadmm_rspa
             DO ispin = 1, nspins
                CALL calculate_rho_elec(ks_env=ks_env, matrix_p=mpz(ispin, 1)%matrix, &
                                        rho=rhoz_r_aux(ispin), rho_gspace=rhoz_g_aux(ispin), &
-                                       basis_type="AUX_FIT", &
-                                       task_list_external=task_list_aux_fit)
+                                       basis_type=basis_type, task_list_external=task_list)
             END DO
             !
             ! Add ADMM volume contribution to stress tensor
@@ -2139,8 +2138,8 @@ SUBROUTINE response_force(qs_env, vh_rspace, vxc_rspace, vtau_rspace, vadmm_rspa
                CALL pw_scale(v_xc(ispin), v_xc(ispin)%pw_grid%dvol)
                CALL integrate_v_rspace(qs_env=qs_env, v_rspace=v_xc(ispin), &
                                        hmat=mhy(ispin, 1), pmat=matrix_p(ispin, 1), &
-                                       calculate_forces=.TRUE., basis_type="AUX_FIT", &
-                                       task_list_external=task_list_aux_fit)
+                                       calculate_forces=.TRUE., &
+                                       basis_type=basis_type, task_list_external=task_list)
             END DO
             IF (debug_forces) THEN
                fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
@@ -2174,7 +2173,8 @@ SUBROUTINE response_force(qs_env, vh_rspace, vxc_rspace, vtau_rspace, vadmm_rspa
                CALL prepare_gapw_den(qs_env, local_rho_set=local_rhoz_set_admm, &
                                      do_rho0=.FALSE., kind_set_external=admm_env%admm_gapw_env%admm_kind_set)
                !compute the potential due to atomic densities
-               CALL calculate_xc_2nd_deriv_atom(ex_env%local_rho_set_admm%rho_atom_set, local_rhoz_set_admm%rho_atom_set, &
+               CALL calculate_xc_2nd_deriv_atom(admm_env%admm_gapw_env%local_rho_set%rho_atom_set, &
+                                                local_rhoz_set_admm%rho_atom_set, &
                                                 qs_env, xc_section, para_env, do_tddft=.FALSE., do_triplet=.FALSE., &
                                                 kind_set_external=admm_env%admm_gapw_env%admm_kind_set)
                !

diff --git a/tests/QS/regtest-tddfpt-force-gapw/TEST_FILES b/tests/QS/regtest-tddfpt-force-gapw/TEST_FILES
@@ -25,4 +25,6 @@ bnt2.inp                                 0
 bnt3.inp                                 0
 #
 ft1.inp                                  0
+ft2.inp                                  0
+ft3.inp                                  0
 #EOF
diff --git a/tests/QS/regtest-tddfpt-force-gapw/ft1.inp b/tests/QS/regtest-tddfpt-force-gapw/ft1.inp
@@ -21,7 +21,7 @@
     &END
     &MGRID
       CUTOFF 180
-      REL_CUTOFF 40
+      REL_CUTOFF 30
     &END MGRID
     &QS
       METHOD GAPW
@@ -33,17 +33,8 @@
     &END EXCITED_STATES
     &SCF
       SCF_GUESS ATOMIC
-      &OT
-         PRECONDITIONER FULL_ALL
-         MINIMIZER DIIS
-         STEPSIZE 0.01
-      &END
-      &OUTER_SCF
-         MAX_SCF  10
-         EPS_SCF 1.0E-7
-      &END
-      MAX_SCF  10
-      EPS_SCF 1.0E-7
+      MAX_SCF  100
+      EPS_SCF 1.0E-6
     &END SCF
     &XC
       &XC_FUNCTIONAL PBE0 
@@ -67,7 +58,7 @@
   &END DFT
   &SUBSYS
     &CELL
-      ABC [angstrom]  5.0  5.0  5.0
+      ABC [angstrom]  4.5  4.5  4.5
       PERIODIC NONE
     &END
     &COORD

diff --git a/tests/QS/regtest-tddfpt-force-gapw/ft2.inp b/tests/QS/regtest-tddfpt-force-gapw/ft2.inp
@@ -0,0 +1,108 @@
+&FORCE_EVAL
+  METHOD Quickstep
+  &PROPERTIES
+    &TDDFPT
+       KERNEL FULL
+       RKS_TRIPLETS F
+       NSTATES      1
+       MAX_ITER    50
+       CONVERGENCE [eV] 1.0e-06
+       ADMM_KERNEL_CORRECTION_SYMMETRIC T
+       ADMM_KERNEL_XC_CORRECTION        F
+    &END TDDFPT
+  &END PROPERTIES
+  &DFT
+    BASIS_SET_FILE_NAME BASIS_SET
+    BASIS_SET_FILE_NAME BASIS_ADMM
+    POTENTIAL_FILE_NAME GTH_POTENTIALS
+    &AUXILIARY_DENSITY_MATRIX_METHOD
+        ADMM_PURIFICATION_METHOD NONE
+        METHOD BASIS_PROJECTION
+        EXCH_SCALING_MODEL NONE
+        EXCH_CORRECTION_FUNC LDA_X_LIBXC
+    &END
+    &MGRID
+      CUTOFF 180
+      REL_CUTOFF 30
+    &END MGRID
+    &QS
+      METHOD GAPW
+      EPS_DEFAULT 1.E-10
+    &END QS
+    &EXCITED_STATES T
+       STATE 1
+       DEBUG_FORCES T
+    &END EXCITED_STATES
+    &SCF
+      SCF_GUESS ATOMIC
+      MAX_SCF  20
+      EPS_SCF 1.0E-6
+    &END SCF
+    &XC
+      &XC_FUNCTIONAL PBE0 
+      &END XC_FUNCTIONAL
+      &HF
+         &SCREENING
+            SCREEN_P_FORCES F
+         &END SCREENING
+      &END HF
+    &END XC
+    &PRINT
+      &MOMENTS ON
+         PERIODIC .FALSE.
+         REFERENCE COM
+      &END
+    &END
+    &POISSON
+       PERIODIC NONE
+       POISSON_SOLVER MT
+    &END
+  &END DFT
+  &SUBSYS
+    &CELL
+      ABC [angstrom]  4.5  4.5  4.5
+      PERIODIC NONE
+    &END
+    &COORD
+C          0.0000000000        0.0000000000       -0.5265526741
+O          0.0000000000        0.0000000000        0.6555124750
+H          0.0000000000       -0.9325664988       -1.1133424527
+H          0.0000000000        0.9325664988       -1.1133424527
+    &END COORD
+    &TOPOLOGY
+     &CENTER_COORDINATES
+     &END
+    &END
+    &KIND H
+      BASIS_SET DZVP-GTH-PADE
+      BASIS_SET AUX_FIT cFIT3
+      POTENTIAL GTH-PADE-q1
+    &END KIND
+    &KIND C
+      BASIS_SET DZVP-GTH-PADE
+      BASIS_SET AUX_FIT cFIT3
+      POTENTIAL GTH-PADE-q4
+    &END KIND
+    &KIND O
+      BASIS_SET DZVP-GTH-PADE
+      BASIS_SET AUX_FIT cFIT3
+      POTENTIAL GTH-PADE-q6
+    &END KIND
+  &END SUBSYS
+&END FORCE_EVAL
+&GLOBAL
+  PRINT_LEVEL LOW
+  PROJECT cc02
+  RUN_TYPE DEBUG
+&END GLOBAL
+
+&DEBUG
+   DEBUG_FORCES .TRUE.
+   DEBUG_STRESS_TENSOR .FALSE.
+   DEBUG_DIPOLE .FALSE.
+   DEBUG_POLARIZABILITY .FALSE.
+   CHECK_ATOM_FORCE   1 z
+   STOP_ON_MISMATCH  T
+   DE 0.001
+   MAX_RELATIVE_ERROR  0.05
+&END
diff --git a/tests/QS/regtest-tddfpt-force-gapw/ft3.inp b/tests/QS/regtest-tddfpt-force-gapw/ft3.inp
@@ -0,0 +1,108 @@
+&FORCE_EVAL
+  METHOD Quickstep
+  &PROPERTIES
+    &TDDFPT
+       KERNEL FULL
+       RKS_TRIPLETS F
+       NSTATES      1
+       MAX_ITER    50
+       CONVERGENCE [eV] 1.0e-06
+       ADMM_KERNEL_CORRECTION_SYMMETRIC T
+       ADMM_KERNEL_XC_CORRECTION        T
+    &END TDDFPT
+  &END PROPERTIES
+  &DFT
+    BASIS_SET_FILE_NAME BASIS_SET
+    BASIS_SET_FILE_NAME BASIS_ADMM
+    POTENTIAL_FILE_NAME GTH_POTENTIALS
+    &AUXILIARY_DENSITY_MATRIX_METHOD
+        ADMM_PURIFICATION_METHOD NONE
+        METHOD BASIS_PROJECTION
+        EXCH_SCALING_MODEL NONE
+        EXCH_CORRECTION_FUNC LDA_X_LIBXC
+    &END
+    &MGRID
+      CUTOFF 180
+      REL_CUTOFF 40
+    &END MGRID
+    &QS
+      METHOD GAPW
+      EPS_DEFAULT 1.E-10
+    &END QS
+    &EXCITED_STATES T
+       STATE 1
+       DEBUG_FORCES T
+    &END EXCITED_STATES
+    &SCF
+      SCF_GUESS ATOMIC
+      MAX_SCF  100
+      EPS_SCF 1.0E-6
+    &END SCF
+    &XC
+      &XC_FUNCTIONAL PBE0 
+      &END XC_FUNCTIONAL
+      &HF
+         &SCREENING
+            SCREEN_P_FORCES F
+         &END SCREENING
+      &END HF
+    &END XC
+    &PRINT
+      &MOMENTS ON
+         PERIODIC .FALSE.
+         REFERENCE COM
+      &END
+    &END
+    &POISSON
+       PERIODIC NONE
+       POISSON_SOLVER MT
+    &END
+  &END DFT
+  &SUBSYS
+    &CELL
+      ABC [angstrom]  5.0  5.0  5.0
+      PERIODIC NONE
+    &END
+    &COORD
+C          0.0000000000        0.0000000000       -0.5265526741
+O          0.0000000000        0.0000000000        0.6555124750
+H          0.0000000000       -0.9325664988       -1.1133424527
+H          0.0000000000        0.9325664988       -1.1133424527
+    &END COORD
+    &TOPOLOGY
+     &CENTER_COORDINATES
+     &END
+    &END
+    &KIND H
+      BASIS_SET DZVP-GTH-PADE
+      BASIS_SET AUX_FIT cFIT3
+      POTENTIAL GTH-PADE-q1
+    &END KIND
+    &KIND C
+      BASIS_SET DZVP-GTH-PADE
+      BASIS_SET AUX_FIT cFIT3
+      POTENTIAL GTH-PADE-q4
+    &END KIND
+    &KIND O
+      BASIS_SET DZVP-GTH-PADE
+      BASIS_SET AUX_FIT cFIT3
+      POTENTIAL GTH-PADE-q6
+    &END KIND
+  &END SUBSYS
+&END FORCE_EVAL
+&GLOBAL
+  PRINT_LEVEL LOW
+  PROJECT cc02
+  RUN_TYPE DEBUG
+&END GLOBAL
+
+&DEBUG
+   DEBUG_FORCES .TRUE.
+   DEBUG_STRESS_TENSOR .FALSE.
+   DEBUG_DIPOLE .FALSE.
+   DEBUG_POLARIZABILITY .FALSE.
+   CHECK_ATOM_FORCE   1 z
+   STOP_ON_MISMATCH  T
+   DE 0.001
+   MAX_RELATIVE_ERROR  0.05
+&END