Finite difference for Linres Kernels (#591)

src/input_cp2k_xc.F

@@ -1230,6 +1230,20 @@ SUBROUTINE create_xc_section(section)
       CALL section_add_keyword(subsection, keyword)
       CALL keyword_release(keyword)
 
+      CALL keyword_create(keyword, __LOCATION__, name="2ND_DERIV_ANALYTICAL", &
+                          description="Use analytical formulas or finite differences for 2nd derivatives of XC", &
+                          usage="2ND_DERIV_ANALYTICAL logical", default_l_val=.TRUE., &
+                          lone_keyword_l_val=.TRUE.)
+      CALL section_add_keyword(section, keyword)
+      CALL keyword_release(keyword)
+
+      CALL keyword_create(keyword, __LOCATION__, name="3RD_DERIV_ANALYTICAL", &
+                          description="Use analytical formulas or finite differences for 3rd derivatives of XC", &
+                          usage="3RD_DERIV_ANALYTICAL logical", default_l_val=.TRUE., &
+                          lone_keyword_l_val=.TRUE.)
+      CALL section_add_keyword(section, keyword)
+      CALL keyword_release(keyword)
+
       CALL section_add_subsection(section, subsection)
       CALL section_release(subsection)
 

src/qs_linres_methods.F

@@ -57,6 +57,7 @@ MODULE qs_linres_methods
                                               ot_precond_solver_default,&
                                               state_loc_all
    USE input_section_types,             ONLY: section_get_ival,&
+                                              section_get_lval,&
                                               section_get_rval,&
                                               section_vals_get_subs_vals,&
                                               section_vals_type,&
@@ -129,10 +130,12 @@ MODULE qs_linres_methods
    USE qs_vxc_atom,                     ONLY: calculate_xc_2nd_deriv_atom
    USE string_utilities,                ONLY: xstring
    USE xc,                              ONLY: xc_calc_2nd_deriv,&
-                                              xc_prep_2nd_deriv
+                                              xc_prep_2nd_deriv,&
+                                              xc_vxc_pw_create
    USE xc_derivatives,                  ONLY: xc_functionals_get_needs
    USE xc_rho_cflags_types,             ONLY: xc_rho_cflags_type
    USE xc_rho_set_types,                ONLY: xc_rho_set_create,&
+                                              xc_rho_set_get,&
                                               xc_rho_set_release,&
                                               xc_rho_set_type,&
                                               xc_rho_set_update
@@ -840,12 +843,17 @@ SUBROUTINE apply_op_2_dft(qs_env, p_env, c0, v, Av, chc)
       TYPE(cp_fm_p_type), DIMENSION(:), POINTER          :: c0, v, Av, chc
 
       CHARACTER(len=*), PARAMETER :: routineN = 'apply_op_2_dft', routineP = moduleN//':'//routineN
+      REAL(KIND=dp), PARAMETER                           :: h = 0.001_dp
 
       INTEGER                                            :: handle, ikind, ispin, ncol, nkind, ns, &
                                                             nspins
       INTEGER, DIMENSION(2, 3)                           :: bo
-      LOGICAL                                            :: gapw, gapw_xc, lr_triplet, lrigpw, lsd
-      REAL(KIND=dp)                                      :: energy_hartree, energy_hartree_1c, fac
+      LOGICAL                                            :: deriv2_analytic, gapw, gapw_xc, &
+                                                            lr_triplet, lrigpw, lsd
+      REAL(KIND=dp)                                      :: energy_hartree, energy_hartree_1c, exc, &
+                                                            fac
+      REAL(KIND=dp), DIMENSION(3, 3)                     :: virial_xc
+      REAL(kind=dp), DIMENSION(:, :, :), POINTER         :: rho3, rhoa, rhob
       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
       TYPE(cp_logger_type), POINTER                      :: logger
       TYPE(cp_para_env_type), POINTER                    :: para_env
@@ -859,8 +867,8 @@ SUBROUTINE apply_op_2_dft(qs_env, p_env, c0, v, Av, chc)
       TYPE(pw_env_type), POINTER                         :: pw_env
       TYPE(pw_p_type)                                    :: rho1_tot_gspace, v_hartree_gspace, &
                                                             v_hartree_rspace
-      TYPE(pw_p_type), DIMENSION(:), POINTER             :: rho1_g, rho1_g_pw, rho1_r, rho1_r_pw, &
-                                                            tau_pw, v_rspace_new, v_xc
+      TYPE(pw_p_type), DIMENSION(:), POINTER :: rho1_g, rho1_g_pw, rho1_r, rho1_r_pw, rho_g, &
+         rho_r, tau, tau_pw, v_rspace_new, v_xc, vxc_rho_1, vxc_rho_2, vxc_rho_3, vxc_rho_4
       TYPE(pw_poisson_type), POINTER                     :: poisson_env
       TYPE(pw_pool_p_type), DIMENSION(:), POINTER        :: pw_pools
       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
@@ -993,52 +1001,159 @@ SUBROUTINE apply_op_2_dft(qs_env, p_env, c0, v, Av, chc)
 
       NULLIFY (tau_pw)
 
-      !------!
-      ! rho1 !
-      !------!
-      bo = rho1_r(1)%pw%pw_grid%bounds_local
-      ! create the place where to store the argument for the functionals
-      CALL xc_rho_set_create(rho1_set, bo, &
-                             rho_cutoff=section_get_rval(xc_section, "DENSITY_CUTOFF"), &
-                             drho_cutoff=section_get_rval(xc_section, "GRADIENT_CUTOFF"), &
-                             tau_cutoff=section_get_rval(xc_section, "TAU_CUTOFF"))
-
-      xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
-      needs = xc_functionals_get_needs(xc_fun_section, lsd, .TRUE.)
-
-      ! calculate the arguments needed by the functionals
-      CALL xc_rho_set_update(rho1_set, rho1_r_pw, rho1_g_pw, tau_pw, needs, &
-                             section_get_ival(xc_section, "XC_GRID%XC_DERIV"), &
-                             section_get_ival(xc_section, "XC_GRID%XC_SMOOTH_RHO"), &
-                             auxbas_pw_pool)
-
-      ALLOCATE (v_xc(nspins))
-      DO ispin = 1, nspins
-         NULLIFY (v_xc(ispin)%pw)
-         CALL pw_pool_create_pw(auxbas_pw_pool, v_xc(ispin)%pw, &
-                                use_data=REALDATA3D, &
-                                in_space=REALSPACE)
-         CALL pw_zero(v_xc(ispin)%pw)
-      END DO
+      deriv2_analytic = section_get_lval(xc_section, "2ND_DERIV_ANALYTICAL")
+
+      IF (deriv2_analytic) THEN
+         !------!
+         ! rho1 !
+         !------!
+         bo = rho1_r(1)%pw%pw_grid%bounds_local
+         ! create the place where to store the argument for the functionals
+         CALL xc_rho_set_create(rho1_set, bo, &
+                                rho_cutoff=section_get_rval(xc_section, "DENSITY_CUTOFF"), &
+                                drho_cutoff=section_get_rval(xc_section, "GRADIENT_CUTOFF"), &
+                                tau_cutoff=section_get_rval(xc_section, "TAU_CUTOFF"))
+
+         xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
+         needs = xc_functionals_get_needs(xc_fun_section, lsd, .TRUE.)
+
+         ! calculate the arguments needed by the functionals
+         CALL xc_rho_set_update(rho1_set, rho1_r_pw, rho1_g_pw, tau_pw, needs, &
+                                section_get_ival(xc_section, "XC_GRID%XC_DERIV"), &
+                                section_get_ival(xc_section, "XC_GRID%XC_SMOOTH_RHO"), &
+                                auxbas_pw_pool)
+
+         ALLOCATE (v_xc(nspins))
+         DO ispin = 1, nspins
+            NULLIFY (v_xc(ispin)%pw)
+            CALL pw_pool_create_pw(auxbas_pw_pool, v_xc(ispin)%pw, &
+                                   use_data=REALDATA3D, &
+                                   in_space=REALSPACE)
+            CALL pw_zero(v_xc(ispin)%pw)
+         END DO
 
-      fac = 0._dp
-      IF (nspins == 1) THEN
-         IF (lr_triplet) fac = -1.0_dp
-      END IF
+         fac = 0._dp
+         IF (nspins == 1) THEN
+            IF (lr_triplet) fac = -1.0_dp
+         END IF
 
-      CALL xc_calc_2nd_deriv(v_xc, p_env%kpp1_env%deriv_set, p_env%kpp1_env%rho_set, &
-                             rho1_set, auxbas_pw_pool, xc_section=xc_section, &
-                             tddfpt_fac=fac)
+         CALL xc_calc_2nd_deriv(v_xc, p_env%kpp1_env%deriv_set, p_env%kpp1_env%rho_set, &
+                                rho1_set, auxbas_pw_pool, xc_section=xc_section, &
+                                tddfpt_fac=fac)
 
-      DO ispin = 1, nspins
-         v_rspace_new(ispin)%pw => v_xc(ispin)%pw
-      END DO
-      DEALLOCATE (v_xc)
+         DO ispin = 1, nspins
+            v_rspace_new(ispin)%pw => v_xc(ispin)%pw
+         END DO
+         DEALLOCATE (v_xc)
 
-      IF (gapw) CALL calculate_xc_2nd_deriv_atom(p_env, qs_env, xc_section, do_tddft=.FALSE., &
-                                                 do_triplet=lr_triplet)
+         IF (gapw) CALL calculate_xc_2nd_deriv_atom(p_env, qs_env, xc_section, do_tddft=.FALSE., &
+                                                    do_triplet=lr_triplet)
 
-      CALL xc_rho_set_release(rho1_set)
+         CALL xc_rho_set_release(rho1_set)
+
+      ELSE
+         NULLIFY (tau)
+         ALLOCATE (v_xc(nspins))
+         DO ispin = 1, nspins
+            NULLIFY (v_xc(ispin)%pw)
+            CALL pw_pool_create_pw(auxbas_pw_pool, v_xc(ispin)%pw, &
+                                   use_data=REALDATA3D, &
+                                   in_space=REALSPACE)
+            CALL pw_zero(v_xc(ispin)%pw)
+         END DO
+         ! finite differences
+         CPASSERT(.NOT. gapw)
+         IF (nspins == 2) THEN
+            NULLIFY (vxc_rho_1, vxc_rho_2, rho_g)
+            ALLOCATE (rho_r(2))
+            DO ispin = 1, nspins
+               NULLIFY (rho_r(ispin)%pw)
+               CALL pw_pool_create_pw(auxbas_pw_pool, rho_r(ispin)%pw, in_space=REALSPACE, use_data=REALDATA3D)
+            END DO
+            CALL xc_rho_set_get(p_env%kpp1_env%rho_set, rhoa=rhoa, rhob=rhob)
+            rho_r(1)%pw%cr3d(:, :, :) = rhoa(:, :, :)+0.5_dp*h*rho1_r(1)%pw%cr3d(:, :, :)
+            rho_r(2)%pw%cr3d(:, :, :) = rhob(:, :, :)+0.5_dp*h*rho1_r(2)%pw%cr3d(:, :, :)
+            CALL xc_vxc_pw_create(vxc_rho_1, tau_pw, exc, rho_r, rho_g, tau, xc_section, &
+                                  auxbas_pw_pool, .FALSE., virial_xc)
+            rho_r(1)%pw%cr3d(:, :, :) = rhoa(:, :, :)-0.5_dp*h*rho1_r(1)%pw%cr3d(:, :, :)
+            rho_r(2)%pw%cr3d(:, :, :) = rhob(:, :, :)-0.5_dp*h*rho1_r(2)%pw%cr3d(:, :, :)
+            CALL xc_vxc_pw_create(vxc_rho_2, tau_pw, exc, rho_r, rho_g, tau, xc_section, &
+                                  auxbas_pw_pool, .FALSE., virial_xc)
+            v_xc(1)%pw%cr3d(:, :, :) = (vxc_rho_1(1)%pw%cr3d(:, :, :)-vxc_rho_2(1)%pw%cr3d(:, :, :))/h
+            v_xc(2)%pw%cr3d(:, :, :) = (vxc_rho_1(2)%pw%cr3d(:, :, :)-vxc_rho_2(2)%pw%cr3d(:, :, :))/h
+            DO ispin = 1, nspins
+               CALL pw_pool_give_back_pw(auxbas_pw_pool, rho_r(ispin)%pw)
+            END DO
+            DEALLOCATE (rho_r)
+            DO ispin = 1, nspins
+               CALL pw_release(vxc_rho_1(ispin)%pw)
+               CALL pw_release(vxc_rho_2(ispin)%pw)
+            END DO
+            DEALLOCATE (vxc_rho_1, vxc_rho_2)
+         ELSE IF (nspins == 1 .AND. (lr_triplet)) THEN
+            NULLIFY (vxc_rho_1, vxc_rho_2, vxc_rho_3, vxc_rho_4, rho_g)
+            ALLOCATE (rho_r(2))
+            DO ispin = 1, 2
+               NULLIFY (rho_r(ispin)%pw)
+               CALL pw_pool_create_pw(auxbas_pw_pool, rho_r(ispin)%pw, in_space=REALSPACE, use_data=REALDATA3D)
+            END DO
+            CALL xc_rho_set_get(p_env%kpp1_env%rho_set, rhoa=rhoa, rhob=rhob)
+            ! K(alpha,alpha)
+            rho_r(1)%pw%cr3d(:, :, :) = rhoa(:, :, :)+0.5_dp*h*rho1_r(1)%pw%cr3d(:, :, :)
+            rho_r(2)%pw%cr3d(:, :, :) = rhob(:, :, :)
+            CALL xc_vxc_pw_create(vxc_rho_1, tau_pw, exc, rho_r, rho_g, tau, xc_section, &
+                                  auxbas_pw_pool, .FALSE., virial_xc)
+            rho_r(1)%pw%cr3d(:, :, :) = rhoa(:, :, :)-0.5_dp*h*rho1_r(1)%pw%cr3d(:, :, :)
+            rho_r(2)%pw%cr3d(:, :, :) = rhob(:, :, :)
+            CALL xc_vxc_pw_create(vxc_rho_2, tau_pw, exc, rho_r, rho_g, tau, xc_section, &
+                                  auxbas_pw_pool, .FALSE., virial_xc)
+            v_xc(1)%pw%cr3d(:, :, :) = (vxc_rho_1(1)%pw%cr3d(:, :, :)-vxc_rho_2(1)%pw%cr3d(:, :, :))/h
+            ! K(alpha,beta)
+            rho_r(1)%pw%cr3d(:, :, :) = rhoa(:, :, :)
+            rho_r(2)%pw%cr3d(:, :, :) = rhob(:, :, :)+0.5_dp*h*rho1_r(1)%pw%cr3d(:, :, :)
+            CALL xc_vxc_pw_create(vxc_rho_3, tau_pw, exc, rho_r, rho_g, tau, xc_section, &
+                                  auxbas_pw_pool, .FALSE., virial_xc)
+            rho_r(1)%pw%cr3d(:, :, :) = rhoa(:, :, :)
+            rho_r(2)%pw%cr3d(:, :, :) = rhob(:, :, :)-0.5_dp*h*rho1_r(1)%pw%cr3d(:, :, :)
+            CALL xc_vxc_pw_create(vxc_rho_4, tau_pw, exc, rho_r, rho_g, tau, xc_section, &
+                                  auxbas_pw_pool, .FALSE., virial_xc)
+            v_xc(1)%pw%cr3d(:, :, :) = v_xc(1)%pw%cr3d(:, :, :)- &
+                                       (vxc_rho_3(1)%pw%cr3d(:, :, :)-vxc_rho_4(1)%pw%cr3d(:, :, :))/h
+            DO ispin = 1, 2
+               CALL pw_pool_give_back_pw(auxbas_pw_pool, rho_r(ispin)%pw)
+            END DO
+            DEALLOCATE (rho_r)
+            DO ispin = 1, 2
+               CALL pw_release(vxc_rho_1(ispin)%pw)
+               CALL pw_release(vxc_rho_2(ispin)%pw)
+               CALL pw_release(vxc_rho_3(ispin)%pw)
+               CALL pw_release(vxc_rho_4(ispin)%pw)
+            END DO
+            DEALLOCATE (vxc_rho_1, vxc_rho_2, vxc_rho_3, vxc_rho_4)
+         ELSE
+            NULLIFY (vxc_rho_1, vxc_rho_2, rho_r, rho_g)
+            ALLOCATE (rho_r(1))
+            NULLIFY (rho_r(1)%pw)
+            CALL pw_pool_create_pw(auxbas_pw_pool, rho_r(1)%pw, in_space=REALSPACE, use_data=REALDATA3D)
+            CALL xc_rho_set_get(p_env%kpp1_env%rho_set, rho=rho3)
+            rho_r(1)%pw%cr3d(:, :, :) = rho3(:, :, :)+0.5_dp*h*rho1_r(1)%pw%cr3d(:, :, :)
+            CALL xc_vxc_pw_create(vxc_rho_1, tau_pw, exc, rho_r, rho_g, tau, xc_section, &
+                                  auxbas_pw_pool, .FALSE., virial_xc)
+            rho_r(1)%pw%cr3d(:, :, :) = rho3(:, :, :)-0.5_dp*h*rho1_r(1)%pw%cr3d(:, :, :)
+            CALL xc_vxc_pw_create(vxc_rho_2, tau_pw, exc, rho_r, rho_g, tau, xc_section, &
+                                  auxbas_pw_pool, .FALSE., virial_xc)
+            v_xc(1)%pw%cr3d(:, :, :) = (vxc_rho_1(1)%pw%cr3d(:, :, :)-vxc_rho_2(1)%pw%cr3d(:, :, :))/h
+            CALL pw_pool_give_back_pw(auxbas_pw_pool, rho_r(1)%pw)
+            DEALLOCATE (rho_r)
+            CALL pw_release(vxc_rho_1(1)%pw)
+            CALL pw_release(vxc_rho_2(1)%pw)
+            DEALLOCATE (vxc_rho_1, vxc_rho_2)
+         END IF
+         DO ispin = 1, nspins
+            v_rspace_new(ispin)%pw => v_xc(ispin)%pw
+         END DO
+         DEALLOCATE (v_xc)
+      END IF
 
       DO ispin = 1, SIZE(rho1_r_pw)
          CALL pw_release(rho1_r_pw(ispin)%pw)

tests/QS/regtest-debug/TEST_FILES

@@ -10,4 +10,5 @@ h2o_periodic.inp                                      87      1e-05
 h2o_gga.inp                                           87      1e-05             0.163373158129E+02
 h2o_gapw.inp                                          87      1e-05             0.163123128994E+02
 h2o_lri.inp                                           87      1e-05             0.216331245542E+02
+h2o_pade_fd.inp                                       87      1e-05             0.164539007391E+02
 #EOF

tests/QS/regtest-debug/h2o_pade_fd.inp

@@ -0,0 +1,85 @@
+&FORCE_EVAL
+  METHOD Quickstep
+  &PROPERTIES
+    &LINRES
+       PRECONDITIONER FULL_ALL
+       EPS 1.e-09
+       &POLAR
+          DO_RAMAN T
+          PERIODIC_DIPOLE_OPERATOR F
+       &END
+    &END
+  &END
+  &DFT
+    LSD
+    &QS
+      METHOD GPW
+      EPS_DEFAULT 1.e-14
+    &END QS
+    &EFIELD
+    &END
+    &SCF
+      SCF_GUESS ATOMIC 
+      &OT
+         PRECONDITIONER FULL_SINGLE_INVERSE
+         MINIMIZER DIIS
+      &END
+      &OUTER_SCF
+         MAX_SCF  10
+         EPS_SCF 1.0E-7
+      &END
+      MAX_SCF  10
+      EPS_SCF 1.0E-7
+    &END SCF
+    &XC
+       &XC_FUNCTIONAL PADE
+       &END XC_FUNCTIONAL
+       2ND_DERIV_ANALYTICAL F
+    &END XC
+    &PRINT
+      &MOMENTS ON
+         PERIODIC .FALSE.
+         REFERENCE COM
+      &END
+    &END
+  &END DFT
+  &SUBSYS
+    &CELL
+      ABC [angstrom] 6.0 6.0 6.0
+      PERIODIC NONE
+    &END
+    &COORD
+    O   0.000000    0.000000   -0.065587
+    H   0.000000   -0.757136    0.520545
+    H   0.000000    0.757136    0.520545
+    &END COORD
+    &TOPOLOGY
+     &CENTER_COORDINATES
+     &END
+    &END
+    &KIND H
+      BASIS_SET DZV-GTH-PADE
+      POTENTIAL GTH-PADE-q1
+    &END KIND
+    &KIND O
+      BASIS_SET DZVP-GTH-PADE
+      POTENTIAL GTH-PADE-q6
+    &END KIND
+  &END SUBSYS
+&END FORCE_EVAL
+&GLOBAL
+  PRINT_LEVEL LOW
+  PROJECT dipole
+  RUN_TYPE DEBUG
+&END GLOBAL
+
+&DEBUG
+   DEBUG_FORCES .FALSE.
+   DEBUG_STRESS_TENSOR .FALSE.
+   DEBUG_DIPOLE .FALSE.
+   DEBUG_POLARIZABILITY .TRUE.
+   DE 0.002
+   EPS_NO_ERROR_CHECK 1.e-4
+&END
+
+