Enable forces calculation for the DFT+U method Lowdin

cp2k · Apr 26, 2021 · aa88a19 · aa88a19
1 parent 2da74a0
commit aa88a19
Show file tree

Hide file tree

Showing 2 changed files with 11 additions and 6 deletions.
diff --git a/src/dft_plus_u.F b/src/dft_plus_u.F
@@ -73,6 +73,7 @@ MODULE dft_plus_u
    USE qs_rho_types,                    ONLY: qs_rho_get,&
                                               qs_rho_type
    USE qs_scf_types,                    ONLY: qs_scf_env_type
+   USE virial_types,                    ONLY: virial_type
 #include "./base/base_uses.f90"
 
    IMPLICIT NONE
@@ -244,6 +245,7 @@ SUBROUTINE lowdin(qs_env, matrix_h, matrix_w, should_output, output_unit, &
       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
       TYPE(qs_rho_type), POINTER                         :: rho
       TYPE(qs_scf_env_type), POINTER                     :: scf_env
+      TYPE(virial_type), POINTER                         :: virial
 
       CALL timeset(routineN, handle)
 
@@ -293,6 +295,7 @@ SUBROUTINE lowdin(qs_env, matrix_h, matrix_w, should_output, output_unit, &
                       particle_set=particle_set, &
                       rho=rho, &
                       scf_env=scf_env, &
+                      virial=virial, &
                       para_env=para_env, &
                       blacs_env=blacs_env)
 
@@ -567,7 +570,7 @@ SUBROUTINE lowdin(qs_env, matrix_h, matrix_w, should_output, output_unit, &
                                  IF (l(jshell, jset) /= lu) CYCLE
                                  DO jsgf = first_sgf(jshell, jset), last_sgf(jshell, jset)
                                     j = j + 1
-                                    q_matrix(i, j) = q_block(isgf, jsgf)
+                                    IF (isgf == jsgf) q_matrix(i, j) = q_block(isgf, jsgf)
                                  END DO ! next contracted spherical Gaussian function "jsgf"
                               END DO ! next shell "jshell"
                            END DO ! next shell set "jset"
@@ -675,6 +678,7 @@ SUBROUTINE lowdin(qs_env, matrix_h, matrix_w, should_output, output_unit, &
                                        IF (isgf == jsgf) THEN
                                           v_block(isgf, isgf) = u_minus_j*(0.5_dp - fspin*q_matrix(i, i))
                                        ELSE
+                                          CPASSERT(ABS(q_matrix(j, i)) < 1.0E-14_dp)
                                           v_block(isgf, jsgf) = -u_minus_j*fspin*q_matrix(j, i)
                                        END IF
                                     END DO ! next contracted spherical Gaussian function "jsgf"
@@ -784,7 +788,8 @@ SUBROUTINE lowdin(qs_env, matrix_h, matrix_w, should_output, output_unit, &
 
          IF (ASSOCIATED(sm_w)) THEN
 
-            CPABORT("Forces are not yet fully implemented for DFT+U method LOWDIN")
+            CPWARN("This is an experimental version of the forces calculation for the DFT+U method LOWDIN")
+            IF (virial%pv_calculate) CPABORT("The stress tensor is not implemented for the DFT+U method LOWDIN")
 
          END IF ! W matrix update requested
 

diff --git a/tests/QS/regtest-plus_u/TEST_FILES b/tests/QS/regtest-plus_u/TEST_FILES
@@ -6,10 +6,10 @@ H2O-rks-diag-mulliken.inp                             24      2e-11          0.0
 H2O-rks-otcg-mulliken.inp                             24      9e-08          0.0401616150180899990
 H2O-uks-diag-mulliken.inp                             24      5e-13          0.0434047728284000000
 H2O-uks-otcg-mulliken.inp                             24      2e-08          0.0434080727848500030
-H2O-rks-diag-lowdin.inp                               24      1E-12          0.0306388860861400000
-H2O-rks-otcg-lowdin.inp                               24      1E-12          0.0306388956677500000
-H2O-uks-diag-lowdin.inp                               24      1E-12          0.0327909187785700000
-H2O-uks-otcg-lowdin.inp                               24      1e-08          0.0327909985079399970
+H2O-rks-diag-lowdin.inp                               24      1E-12          0.07674874415138
+H2O-rks-otcg-lowdin.inp                               24      1E-12          0.07674870925125
+H2O-uks-diag-lowdin.inp                               24      1E-12          0.08364867485442
+H2O-uks-otcg-lowdin.inp                               24      1e-08          0.08365011953223
 # U ramping
 H2O-rks-u_ramping.inp                                 24      1e-06          0.0401615789381200020
 H2O-uks-u_ramping.inp                                 24      3e-07          0.0401667259122000030