TDDFT Forces (beta version) (#1670)

* Excited state forces: input and flow logic Collect all density matrix type routines into qs_density_matrices TDDFPT with kernel = NONE * TDA dipole moments Kernel forces LSD Triplet case (forces) * Further work on TDDFPT forces (ADMM part) * Hybrid functinal TDDFPT forces with regtests * HFX/ADMM (ground state) TDDFT forces + regtests HFX/ADMM kernel forces (w/o correction term) + regtests * ADMM X correction TDA forces + regtests * Towards xTB/sTDA forces + regtests * TB excited state dipoles * xTB/sTDA excited state dipoles + regtests * xTB/sTDA Ewald forces + regtests * HFX/ADMM functional simplifications, many new regtests * Hybrid functionals polarizability tests activated * Resolve merge conflicts Reduce regtest times * Merge codes Harris forces and TDDFT forces
cp2k · Oct 11, 2021 · a938171 · a938171
1 parent 15378df
commit a938171
Show file tree

Hide file tree

Showing 303 changed files with 24,127 additions and 1,560 deletions.
diff --git a/src/common/bibliography.F b/src/common/bibliography.F
@@ -1982,8 +1982,8 @@ SUBROUTINE add_all_references()
 
       CALL add_reference(key=Kolafa2004, ISI_record=s2a( &
                          "AU Kolafa, J", &
-                         "TI Time-reversible always stable predictor-corrector method for "// &
-                         "molecular dynamics of polarizable molecules", &
+                         "TI Time-reversible always stable predictor-corrector method for ", &
+                         "   molecular dynamics of polarizable molecules", &
                          "SO JOURNAL OF COMPUTATIONAL CHEMISTRY", &
                          "SN 0192-8651", &
                          "J9 J COMPUT CHEM", &

diff --git a/src/cp2k_debug.F b/src/cp2k_debug.F
@@ -407,14 +407,27 @@ SUBROUTINE cp2k_debug_energy_and_forces(force_env)
                      WRITE (UNIT=iw, FMT="(T13,A1,T21,F16.8,T38,F16.8,T56,G12.3,T72,F9.3)") &
                         ACHAR(119 + k), dipole_numer(k), dipole_moment(k), dd, derr
                   ELSE
+                     derr = 0.0_dp
                      WRITE (UNIT=iw, FMT="(T13,A1,T21,F16.8,T38,F16.8,T56,G12.3)") &
                         ACHAR(119 + k), dipole_numer(k), dipole_moment(k), dd
                   END IF
+                  err(k) = derr
                END DO
                WRITE (UNIT=iw, FMT="((T2,A))") &
                   "DEBUG|========================================================================="
-               WRITE (UNIT=iw, FMT="(T2,A,T61,E20.12)") ' DIPOLE : CheckSum  =', SUM(dipole_moment)
+               WRITE (UNIT=iw, FMT="(T2,A,T61,E20.12)") 'DIPOLE : CheckSum  =', SUM(dipole_moment)
+               IF (ANY(ABS(err(1:3)) > 0.1_dp)) THEN
+                  message = "A mismatch between analytical and numerical dipoles "// &
+                            "has been detected. Check the implementation of the "// &
+                            "analytical dipole calculation"
+                  IF (stop_on_mismatch) THEN
+                     CPABORT(message)
+                  ELSE
+                     CPWARN(message)
+                  END IF
+               END IF
             END IF
+
          ELSE
             CALL cp_warn(__LOCATION__, "Debug of dipole moments only for Quickstep code available")
          END IF

diff --git a/src/cp_control_types.F b/src/cp_control_types.F
@@ -117,6 +117,7 @@ MODULE cp_control_types
       LOGICAL                              :: xb_interaction
       LOGICAL                              :: do_nonbonded
       LOGICAL                              :: coulomb_interaction
+      LOGICAL                              :: coulomb_lr
       LOGICAL                              :: tb3_interaction
       LOGICAL                              :: check_atomic_charges
       !
@@ -441,6 +442,9 @@ MODULE cp_control_types
       INTEGER                              :: nprocs
       !> type of kernel function/approximation to use
       INTEGER                              :: kernel
+      !> fro full kernel, do we have HFX/ADMM
+      LOGICAL                              :: do_hfx
+      LOGICAL                              :: do_admm
       !> options used in sTDA calculation (Kernel)
       TYPE(stda_control_type)              :: stda_control
       !> algorithm to correct orbital energies
@@ -458,6 +462,8 @@ MODULE cp_control_types
       LOGICAL                              :: is_restart
       !> compute triplet excited states using spin-unpolarised molecular orbitals
       LOGICAL                              :: rks_triplets
+      !> use symmetric definition of ADMM Kernel correction
+      LOGICAL                              :: admm_symm
       !
       ! DIPOLE_MOMENTS subsection
       !

diff --git a/src/cp_control_utils.F b/src/cp_control_utils.F
@@ -1157,6 +1157,8 @@ SUBROUTINE read_qs_section(qs_control, qs_section)
          ! For debug purposes
          CALL section_vals_val_get(xtb_section, "COULOMB_INTERACTION", &
                                    l_val=qs_control%xtb_control%coulomb_interaction)
+         CALL section_vals_val_get(xtb_section, "COULOMB_LR", &
+                                   l_val=qs_control%xtb_control%coulomb_lr)
          CALL section_vals_val_get(xtb_section, "TB3_INTERACTION", &
                                    l_val=qs_control%xtb_control%tb3_interaction)
          ! Check for bad atomic charges
@@ -1265,6 +1267,7 @@ SUBROUTINE read_tddfpt2_control(t_control, t_section, qs_control)
 
       CALL section_vals_val_get(t_section, "RESTART", l_val=t_control%is_restart)
       CALL section_vals_val_get(t_section, "RKS_TRIPLETS", l_val=t_control%rks_triplets)
+      CALL section_vals_val_get(t_section, "ADMM_KERNEL_CORRECTION_SYMMETRIC", l_val=t_control%admm_symm)
 
       IF (t_control%conv < 0) &
          t_control%conv = ABS(t_control%conv)

diff --git a/src/ec_env_types.F b/src/ec_env_types.F
@@ -13,6 +13,8 @@
 ! **************************************************************************************************
 MODULE ec_env_types
    USE cp_dbcsr_operations,             ONLY: dbcsr_deallocate_matrix_set
+   USE cp_fm_types,                     ONLY: cp_fm_p_type,&
+                                              cp_fm_release
    USE dbcsr_api,                       ONLY: dbcsr_p_type
    USE dm_ls_scf_types,                 ONLY: ls_scf_env_type,&
                                               ls_scf_release
@@ -22,6 +24,8 @@ MODULE ec_env_types
                                               pw_release
    USE qs_dispersion_types,             ONLY: qs_dispersion_release,&
                                               qs_dispersion_type
+   USE qs_force_types,                  ONLY: deallocate_qs_force,&
+                                              qs_force_type
    USE qs_neighbor_list_types,          ONLY: neighbor_list_set_p_type,&
                                               release_neighbor_list_sets
    USE qs_p_env_types,                  ONLY: p_env_release,&
@@ -66,6 +70,8 @@ MODULE ec_env_types
       REAL(KIND=dp)                                    :: etotal
       REAL(KIND=dp)                                    :: eband, exc, ehartree, vhxc
       REAL(KIND=dp)                                    :: edispersion, efield_nuclear
+      ! forces
+      TYPE(qs_force_type), DIMENSION(:), POINTER       :: force => Null()
       ! full neighbor lists and corresponding task list
       TYPE(neighbor_list_set_p_type), &
          DIMENSION(:), POINTER                         :: sab_orb, sac_ppl, sap_ppnl
@@ -85,6 +91,7 @@ MODULE ec_env_types
       TYPE(dbcsr_p_type), DIMENSION(:), POINTER        :: mao_coef
       ! CP equations
       TYPE(qs_p_env_type), POINTER                     :: p_env
+      TYPE(cp_fm_p_type), DIMENSION(:), POINTER        :: cpmos
       TYPE(dbcsr_p_type), DIMENSION(:), POINTER        :: matrix_hz, matrix_z, matrix_wz, z_admm
       ! potentials from input density
       TYPE(pw_p_type), POINTER                         :: vh_rspace
@@ -115,6 +122,8 @@ SUBROUTINE ec_env_release(ec_env)
          CALL release_neighbor_list_sets(ec_env%sab_orb)
          CALL release_neighbor_list_sets(ec_env%sac_ppl)
          CALL release_neighbor_list_sets(ec_env%sap_ppnl)
+         ! forces
+         IF (ASSOCIATED(ec_env%force)) CALL deallocate_qs_force(ec_env%force)
          ! operator matrices
          IF (ASSOCIATED(ec_env%matrix_ks)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_ks)
          IF (ASSOCIATED(ec_env%matrix_h)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_h)
@@ -131,6 +140,14 @@ SUBROUTINE ec_env_release(ec_env)
          IF (ASSOCIATED(ec_env%dispersion_env)) THEN
             CALL qs_dispersion_release(ec_env%dispersion_env)
          END IF
+         ! CP env
+         IF (ASSOCIATED(ec_env%cpmos)) THEN
+            DO iab = 1, SIZE(ec_env%cpmos)
+               CALL cp_fm_release(ec_env%cpmos(iab)%matrix)
+            END DO
+            DEALLOCATE (ec_env%cpmos)
+            NULLIFY (ec_env%cpmos)
+         END IF
 
          IF (ASSOCIATED(ec_env%matrix_z)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_z)
          IF (ASSOCIATED(ec_env%matrix_hz)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_hz)