MP2: Add ADMM and add forces and stress tensor for other operators

Refactor calculation of non-separable part
Add more operators for GPW plus stress tensor terms
Add terms if potential != metric
Set default value for EPS_SVD to 0.0
Rebuild MP2 gradient code and add ADMM stress tensors and numerical XC derivatives

src/admm_methods.F

@@ -37,7 +37,8 @@ MODULE admm_methods
                                               cp_fm_cholesky_reduce,&
                                               cp_fm_cholesky_restore
    USE cp_fm_diag,                      ONLY: cp_fm_syevd
-   USE cp_fm_types,                     ONLY: cp_fm_p_type,&
+   USE cp_fm_types,                     ONLY: cp_fm_get_info,&
+                                              cp_fm_p_type,&
                                               cp_fm_set_all,&
                                               cp_fm_set_element,&
                                               cp_fm_to_fm,&
@@ -125,7 +126,8 @@ MODULE admm_methods
              calc_aux_mo_derivs_none, &
              scale_dm, &
              admm_fit_mo_coeffs, &
-             admm_update_ks_atom
+             admm_update_ks_atom, &
+             admm_aux_reponse_density
 
    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'admm_methods'
 
@@ -2448,5 +2450,49 @@ FUNCTION Heaviside(x)
          Heaviside = 1.0_dp
       END IF
    END FUNCTION Heaviside
+ 
+! **************************************************************************************************
+!> \brief Calculate ADMM auxiliary response density
+!> \param qs_env ...
+!> \param dm ...
+!> \param dm_admm ...
+! **************************************************************************************************
+   SUBROUTINE admm_aux_reponse_density(qs_env, dm, dm_admm)
+      TYPE(qs_environment_type), INTENT(IN), POINTER     :: qs_env
+      TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN)       :: dm
+      TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT)    :: dm_admm
+
+      CHARACTER(LEN=*), PARAMETER :: routineN = 'admm_aux_reponse_density'
+
+      INTEGER                                            :: handle, ispin, nao, nao_aux, ncol, nspins
+      TYPE(admm_type), POINTER                           :: admm_env
+      TYPE(dft_control_type), POINTER                    :: dft_control
+
+      CALL timeset(routineN, handle)
+
+      CALL get_qs_env(qs_env, admm_env=admm_env, dft_control=dft_control)
+
+      nspins = dft_control%nspins
+
+      CPASSERT(ASSOCIATED(admm_env%A))
+      CPASSERT(ASSOCIATED(admm_env%work_orb_orb))
+      CPASSERT(ASSOCIATED(admm_env%work_aux_orb))
+      CPASSERT(ASSOCIATED(admm_env%work_aux_aux))
+      CALL cp_fm_get_info(admm_env%A, nrow_global=nao_aux, ncol_global=nao)
+
+      ! P1 -> AUX BASIS
+      CALL cp_fm_get_info(admm_env%work_orb_orb, nrow_global=nao, ncol_global=ncol)
+      DO ispin = 1, nspins
+         CALL copy_dbcsr_to_fm(dm(ispin)%matrix, admm_env%work_orb_orb)
+         CALL cp_gemm('N', 'N', nao_aux, ncol, nao, 1.0_dp, admm_env%A, &
+                      admm_env%work_orb_orb, 0.0_dp, admm_env%work_aux_orb)
+         CALL cp_gemm('N', 'T', nao_aux, nao_aux, nao, 1.0_dp, admm_env%A, &
+                      admm_env%work_aux_orb, 0.0_dp, admm_env%work_aux_aux)
+         CALL copy_fm_to_dbcsr(admm_env%work_aux_aux, dm_admm(ispin)%matrix, keep_sparsity=.TRUE.)
+      ENDDO
+
+      CALL timestop(handle)
+
+   END SUBROUTINE admm_aux_reponse_density
 
 END MODULE admm_methods

src/common/mathlib.F

@@ -1482,9 +1482,10 @@ END FUNCTION transpose_3d
 !>      05.2007 Created
 !> \author Manuel Guidon (adapted from Numerical recipies)
 ! **************************************************************************************************
-   FUNCTION expint(n, x)
-      INTEGER                                            :: n
-      REAL(dp)                                           :: x, expint
+   ELEMENTAL IMPURE FUNCTION expint(n, x)
+      INTEGER, INTENT(IN)                                :: n
+      REAL(dp), INTENT(IN)                               :: x
+      REAL(dp)                                           :: expint
 
       INTEGER, PARAMETER                                 :: maxit = 100
       REAL(dp), PARAMETER :: eps = 6.e-14_dp, euler = 0.5772156649015328606065120_dp, &

src/ec_env_types.F

@@ -89,7 +89,7 @@ MODULE ec_env_types
       TYPE(dbcsr_p_type), DIMENSION(:), POINTER        :: matrix_hz
       ! potentials from input density
       TYPE(pw_p_type), POINTER                         :: vh_rspace
-      TYPE(pw_p_type), DIMENSION(:), POINTER           :: vxc_rspace, vtau_rspace
+      TYPE(pw_p_type), DIMENSION(:), POINTER           :: vxc_rspace, vtau_rspace, vadmm_rspace
       ! efield
       TYPE(efield_berry_type), POINTER                 :: efield => NULL()
    END TYPE energy_correction_type
@@ -164,6 +164,12 @@ SUBROUTINE ec_env_release(ec_env)
             END DO
             DEALLOCATE (ec_env%vtau_rspace)
          END IF
+         IF (ASSOCIATED(ec_env%vadmm_rspace)) THEN
+            DO iab = 1, SIZE(ec_env%vadmm_rspace)
+               CALL pw_release(ec_env%vadmm_rspace(iab)%pw)
+            END DO
+            DEALLOCATE (ec_env%vadmm_rspace)
+         END IF
          CALL efield_berry_release(ec_env%efield)
 
          DEALLOCATE (ec_env)

src/ec_environment.F

@@ -116,6 +116,7 @@ SUBROUTINE init_ec_env(qs_env, ec_env, dft_section)
       NULLIFY (ec_env%vh_rspace)
       NULLIFY (ec_env%vxc_rspace)
       NULLIFY (ec_env%vtau_rspace)
+      NULLIFY (ec_env%vadmm_rspace)
       ec_env%should_update = .TRUE.
       ec_env%mao = .FALSE.
 

src/energy_corrections.F

@@ -156,12 +156,8 @@ MODULE energy_corrections
                                               xc_prep_2nd_deriv
    USE xc_derivative_set_types,         ONLY: xc_derivative_set_type,&
                                               xc_dset_release
-   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_release,&
-                                              xc_rho_set_type,&
-                                              xc_rho_set_update
+   USE xc_rho_set_types,                ONLY: xc_rho_set_release,&
+                                              xc_rho_set_type
 #include "./base/base_uses.f90"
 
    IMPLICIT NONE
@@ -310,7 +306,7 @@ SUBROUTINE harris_energy(qs_env, ec_env, calculate_forces, unit_nr)
          CALL ec_build_neighborlist(qs_env, ec_env)
          !
          CALL ec_build_core_hamiltonian(qs_env, ec_env)
-         CALL ks_ref_potential(qs_env, ec_env%vh_rspace, ec_env%vxc_rspace, ec_env%vtau_rspace, &
+         CALL ks_ref_potential(qs_env, ec_env%vh_rspace, ec_env%vxc_rspace, ec_env%vtau_rspace, ec_env%vadmm_rspace, &
                                ec_env%ehartree, exc)
          CALL ec_build_ks_matrix(qs_env, ec_env)
          !
@@ -336,7 +332,7 @@ SUBROUTINE harris_energy(qs_env, ec_env, calculate_forces, unit_nr)
          CALL response_equation(qs_env, ec_env%p_env, ec_env%cpmos)
          !
          CALL response_force(qs_env, ec_env%p_env, &
-                             ec_env%vh_rspace, ec_env%vxc_rspace, ec_env%vtau_rspace, &
+                             ec_env%vh_rspace, ec_env%vxc_rspace, ec_env%vtau_rspace, ec_env%vadmm_rspace, &
                              ec_env%matrix_hz)
          !
          CALL ec_properties(qs_env, ec_env)
@@ -896,13 +892,6 @@ SUBROUTINE ec_build_ks_matrix_force(qs_env, ec_env)
       ! Pulay force from Tr P_in (V_H(drho)+ Fxc(rho_in)*drho)
       ! RHS of CPKS equations: (V_H(drho)+ Fxc(rho_in)*drho)*C0
       ! Fxc*drho term
-      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
       CALL get_qs_env(qs_env, input=input)
       xc_section => ec_env%xc_section
       CALL calc_ks_response_potentials(v_xc, xc_section, qs_env, rhoout_r, rhoout_g)
@@ -1065,58 +1054,47 @@ END SUBROUTINE ec_build_ks_matrix_force
 !> \param rho1_r ...
 !> \param rho1_g ...
 ! **************************************************************************************************
-   SUBROUTINE calc_ks_response_potentials(v_xc, xc_section, qs_env, rho1_r, rho1_g, compute_virial, virial)
+   SUBROUTINE calc_ks_response_potentials(v_xc, xc_section, qs_env, rho1_r, rho1_g, compute_virial, virial, do_admm)
       TYPE(pw_p_type), DIMENSION(:), INTENT(IN), POINTER :: v_xc
       TYPE(section_vals_type), INTENT(IN), POINTER       :: xc_section
       TYPE(qs_environment_type), POINTER                 :: qs_env
       TYPE(pw_p_type), DIMENSION(:), INTENT(IN), POINTER :: rho1_r, rho1_g
-      LOGICAL, INTENT(IN), OPTIONAL :: compute_virial
+      LOGICAL, INTENT(IN), OPTIONAL :: compute_virial, do_admm
       REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT), OPTIONAL :: virial
 
       CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_ks_response_potentials'
 
-      INTEGER                                            :: handle, nspins
-      INTEGER, DIMENSION(2, 3)                           :: bo
-      LOGICAL                                            :: lsd
+      INTEGER                                            :: handle
+      LOGICAL                                            :: my_do_admm
       TYPE(pw_env_type), POINTER                         :: pw_env
       TYPE(pw_p_type), DIMENSION(:), POINTER             :: rho_r, tau_r
       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
       TYPE(qs_ks_env_type), POINTER                      :: ks_env
       TYPE(qs_rho_type), POINTER                         :: rho
-      TYPE(section_vals_type), POINTER                   :: xc_fun_section
       TYPE(xc_derivative_set_type), POINTER              :: deriv_set
-      TYPE(xc_rho_cflags_type)                           :: needs
-      TYPE(xc_rho_set_type), POINTER                     :: rho1_set, rho_set
+      TYPE(xc_rho_set_type), POINTER                     :: rho_set
 
       CALL timeset(routineN, handle)
 
-      CALL get_qs_env(qs_env, ks_env=ks_env, rho=rho, pw_env=pw_env)
+      my_do_admm = .FALSE.
+      IF (PRESENT(do_admm)) my_do_admm = do_admm
+
+      IF (.NOT. my_do_admm) THEN
+         CALL get_qs_env(qs_env, ks_env=ks_env, rho=rho, pw_env=pw_env)
+      ELSE
+         CALL get_qs_env(qs_env, ks_env=ks_env, rho_aux_fit=rho, pw_env=pw_env)
+      END IF
       CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
       CALL qs_rho_get(rho, rho_r=rho_r, tau_r=tau_r)
 
-      nspins = SIZE(v_xc, 1)
-      lsd = (nspins == 2)
-      NULLIFY (deriv_set, rho_set, rho1_set)
+      NULLIFY (deriv_set, rho_set)
       CALL xc_prep_2nd_deriv(deriv_set, rho_set, rho_r, auxbas_pw_pool, xc_section=xc_section)
-      bo = rho1_r(1)%pw%pw_grid%bounds_local
-      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, rho1_g, tau_r, needs, &
-                             section_get_ival(xc_section, "XC_GRID%XC_DERIV"), &
-                             section_get_ival(xc_section, "XC_GRID%XC_SMOOTH_RHO"), &
-                             auxbas_pw_pool)
+
       CALL xc_calc_2nd_deriv(v_xc, deriv_set, rho_set, &
-                             rho1_set, auxbas_pw_pool, xc_section=xc_section, compute_virial=compute_virial, virial_xc=virial)
+                             rho1_r, rho1_g, auxbas_pw_pool, xc_section=xc_section, gapw=.FALSE., &
+                             compute_virial=compute_virial, virial_xc=virial)
       CALL xc_dset_release(deriv_set)
       CALL xc_rho_set_release(rho_set)
-      CALL xc_rho_set_release(rho1_set)
 
       CALL timestop(handle)
 
@@ -1130,7 +1108,6 @@ END SUBROUTINE calc_ks_response_potentials
 !>      03.2014 created [JGH]
 !> \author JGH
 ! **************************************************************************************************
-
    SUBROUTINE ec_ks_solver(qs_env, ec_env)
 
       TYPE(qs_environment_type), POINTER                 :: qs_env

src/fm/cp_fm_basic_linalg.F

@@ -38,6 +38,7 @@ MODULE cp_fm_basic_linalg
 
    PUBLIC :: cp_fm_scale, & ! scale a matrix
              cp_fm_scale_and_add, & ! scale and add two matrices
+             cp_fm_geadd, & ! general addition
              cp_fm_column_scale, & ! scale columns of a matrix
              cp_fm_row_scale, & ! scale rows of a matrix
              cp_fm_trace, & ! trace of the transpose(A)*B

src/hfx_admm_utils.F

@@ -453,7 +453,11 @@ SUBROUTINE hfx_ks_matrix(qs_env, matrix_ks, rho, energy, calculate_forces, &
                CALL scale_dm(qs_env, rho_ao_orb, scale_back=.FALSE.)
             END IF
             IF (ASSOCIATED(qs_env%mp2_env)) THEN
-               CALL get_qs_env(qs_env, matrix_p_mp2=rho_ao_resp)
+               IF (dft_control%do_admm) THEN
+                  CALL get_qs_env(qs_env, matrix_p_mp2_admm=rho_ao_resp)
+               ELSE
+                  CALL get_qs_env(qs_env, matrix_p_mp2=rho_ao_resp)
+               END IF
             ELSE
                NULLIFY (rho_ao_resp)
             END IF
@@ -653,9 +657,9 @@ SUBROUTINE create_admm_xc_section(qs_env, xc_section, admm_env)
 
       !in case of no admm exchange corr., no auxiliary exchange functional needed
       IF (admm_env%aux_exch_func == do_admm_aux_exch_func_none) THEN
-         hfx_fraction = 0.0_dp
          CALL section_vals_val_set(xc_fun_section, "_SECTION_PARAMETERS_", &
-                                   i_val=xc_none)
+                                i_val=xc_none)
+         hfx_fraction = 0.0_dp
       ELSE IF (admm_env%aux_exch_func == do_admm_aux_exch_func_default) THEN
          ! default PBE Functional
          !! ** Add functionals evaluated with auxiliary basis

src/hfx_energy_potential.F

@@ -469,7 +469,6 @@ SUBROUTINE integrate_four_center(qs_env, x_data, ks_matrix, ehfx, rho_ao, hfx_se
       do_p_screening = screening_parameter%do_initial_p_screening
       ! Special treatment for MP2 with initial density screening
       IF (do_p_screening) THEN
-         nspins = dft_control%nspins
          IF (ASSOCIATED(qs_env%mp2_env)) THEN
             IF ((qs_env%mp2_env%ri_mp2%free_hfx_buffer)) THEN
                do_p_screening = ((qs_env%mp2_env%p_screen) .AND. (qs_env%mp2_env%not_last_hfx))

src/input_cp2k_mp2.F

@@ -1302,25 +1302,30 @@ SUBROUTINE create_mp2_potential(section)
          description="Which interaction potential should be used "// &
          "(Coulomb, TShPSC operator).", &
          usage="POTENTIAL_TYPE TSHPSC", &
-         enum_c_vals=s2a("COULOMB", "TShPSC", "LONGRANGE"), &
+         enum_c_vals=s2a("COULOMB", "TShPSC", "LONGRANGE", "SHORTRANGE", "TRUNCATED"), &
          enum_i_vals=(/do_potential_coulomb, &
                        do_potential_TShPSC, &
-                       do_potential_long/), &
+                       do_potential_long, &
+                       do_potential_short, &
+                       do_potential_truncated/), &
          enum_desc=s2a("Coulomb potential: 1/r", &
                        "TShPSC:<ul><li>1/x - s/Rc for x &le; Rc</li>"// &
                        "<li>(1 - s)/Rc - (x - Rc)/Rc^2 + (x - Rc)^2/Rc^3 - "// &
                        "(2*n^2 - 7*n + 9 - 4*s)*(x - Rc)^3/(Rc^4*(n^2 - 2*n + 1)*(n - 1)) + "// &
                        "(6-3*s - 4*n + n^2)*(x - Rc)^4/(Rc^5*(n^4 - 4*n^3 + 6*n^2 - 4*n + 1)) "// &
                        "for Rc &lt; x &le; n*Rc (4th order polynomial)</li>"// &
                        "<li>0 for x &gt; n*Rc</li></ul>", &
-                       "Longrange Coulomg potential: erf(w*r)/r)"), &
+                       "Longrange Coulomb potential: erf(w*r)/r)", &
+                       "Shortrange Coulomb potential: erfc(w*r)/r", &
+                       "Truncated Coulomb potential"), &
          default_i_val=do_potential_coulomb)
       CALL section_add_keyword(section, keyword)
       CALL keyword_release(keyword)
 
       CALL keyword_create(keyword, __LOCATION__, name="TRUNCATION_RADIUS", &
-                          description="Determines truncation radius for the truncated TShPSC potential. "// &
-                          "Only valid when doing truncated calculation", &
+                          variants=(/"CUTOFF_RADIUS"/), &
+                          description="Determines truncation radius for the truncated potentials. "// &
+                          "Only valid when doing truncated calculations", &
                           usage="TRUNCATION_RADIUS 10.0", type_of_var=real_t, &
                           default_r_val=10.0_dp, &
                           unit_str="angstrom")
@@ -1329,9 +1334,10 @@ SUBROUTINE create_mp2_potential(section)
 
       CALL keyword_create( &
          keyword, __LOCATION__, &
-         name="TShPSC_DATA", &
-         description="Location of the file TShPSC.dat that contains the data for the "// &
-         "evaluation of the TShPSC G0 ", &
+         name="POTENTIAL_DATA", &
+         variants=s2a("TShPSC_DATA", "T_C_G_DATA"), &
+         description="Location of the file TShPSC.dat or t_c_g.dat that contains the data for the "// &
+         "evaluation of the evaluation of the truncated potentials", &
          usage="TShPSC_DATA t_sh_p_s_c.dat", &
          default_c_val="t_sh_p_s_c.dat")
       CALL section_add_keyword(section, keyword)
@@ -1340,7 +1346,8 @@ SUBROUTINE create_mp2_potential(section)
       CALL keyword_create( &
          keyword, __LOCATION__, &
          name="OMEGA", &
-         description="Range separation parameter for the longrange potential. Only valid when longrange potential is requested.", &
+         description="Range separation parameter for the longrange or shortrange potential. "// &
+                     "Only valid when longrange or shortrange potential is requested.", &
          usage="OMEGA 0.5", type_of_var=real_t, &
          default_r_val=0.5_dp)
       CALL section_add_keyword(section, keyword)
@@ -1419,14 +1426,6 @@ SUBROUTINE create_ri_section(section)
       CALL section_add_keyword(section, keyword)
       CALL keyword_release(keyword)
 
-      CALL keyword_create(keyword, __LOCATION__, name="EPS_SVD", &
-                          description="Include all singular vectors with a singular value smaller than EPS_SVD. "// &
-                          "Is used to provide info on smallest eigenvalues.", &
-                          usage="EPS_SVD  1.0e-5", &
-                          default_r_val=1.0E-3_dp)
-      CALL section_add_keyword(section, keyword)
-      CALL keyword_release(keyword)
-
       CALL keyword_create(keyword, __LOCATION__, name="ERI_BLKSIZE", &
                           description="block sizes for tensors (only used if ERI_METHOD=MME). First value "// &
                           "is the block size for ORB basis, second value is the block size for RI_AUX basis.", &

src/input_cp2k_xc.F

@@ -1251,6 +1251,13 @@ SUBROUTINE create_xc_section(section)
       CALL section_add_keyword(section, keyword)
       CALL keyword_release(keyword)
 
+      CALL keyword_create(keyword, __LOCATION__, name="STEP_SIZE", &
+                          description="Step size in terms of the first order potential for the numerical "// &
+                          "evaluation of XC derivatives", &
+                          usage="STEP_SIZE 1.0E-3", default_r_val=1e-3_dp)
+      CALL section_add_keyword(section, keyword)
+      CALL keyword_release(keyword)
+
       CALL section_add_subsection(section, subsection)
       CALL section_release(subsection)