Collect further stress tensor components

src/hfx_admm_utils.F

@@ -557,6 +557,7 @@ SUBROUTINE hfx_ks_matrix(qs_env, matrix_ks, rho, energy, calculate_forces, &
          END DO
       END DO
       IF (use_virial .AND. calculate_forces) THEN
+         virial%pv_exx = virial%pv_exx - virial%pv_fock_4c
          virial%pv_virial = virial%pv_virial - virial%pv_fock_4c
          virial%pv_calculate = .FALSE.
       ENDIF

src/lri_forces.F

@@ -389,8 +389,10 @@ SUBROUTINE calculate_v_dadr_sr(lri_env, lri_density, pmatrix, cell_to_index, ato
                   !periodic self-pairs aa' contribute only with factor 0.5
 !$OMP CRITICAL(addvirial)
                   IF (iatom == jatom) THEN
+                     CALL virial_pair_force(virial%pv_lrigpw, 0.5_dp, force_a, rab)
                      CALL virial_pair_force(virial%pv_virial, 0.5_dp, force_a, rab)
                   ELSE
+                     CALL virial_pair_force(virial%pv_lrigpw, 1.0_dp, force_a, rab)
                      CALL virial_pair_force(virial%pv_virial, 1.0_dp, force_a, rab)
                   ENDIF
 !$OMP END CRITICAL(addvirial)
@@ -642,8 +644,10 @@ SUBROUTINE calculate_v_dadr_ff(lri_env, lri_density, pmatrix, cell_to_index, ato
                   !periodic self-pairs aa' contribute only with factor 0.5
 !$OMP CRITICAL(addvirial)
                   IF (iatom == jatom) THEN
+                     CALL virial_pair_force(virial%pv_lrigpw, 0.5_dp, force_a, rab)
                      CALL virial_pair_force(virial%pv_virial, 0.5_dp, force_a, rab)
                   ELSE
+                     CALL virial_pair_force(virial%pv_lrigpw, 1.0_dp, force_a, rab)
                      CALL virial_pair_force(virial%pv_virial, 1.0_dp, force_a, rab)
                   ENDIF
 !$OMP END CRITICAL(addvirial)
@@ -835,6 +839,8 @@ SUBROUTINE calculate_v_dadr_ri(lri_env, lri_density, pmatrix, atomic_kind_set, &
             ! to be debugged
             fscal = 1.0_dp
             IF (iatom == jatom) fscal = 1.0_dp
+            CALL virial_pair_force(virial%pv_lrigpw, 1.0_dp, force_a, rik)
+            CALL virial_pair_force(virial%pv_lrigpw, 1.0_dp, force_b, rjk)
             CALL virial_pair_force(virial%pv_virial, 1.0_dp, force_a, rik)
             CALL virial_pair_force(virial%pv_virial, 1.0_dp, force_b, rjk)
          ENDIF

src/mp2_cphf.F

@@ -849,6 +849,7 @@ SUBROUTINE solve_z_vector_eq(qs_env, mp2_env, para_env, dft_control, &
                h_stress(beta, alpha) = h_stress(alpha, beta)
             END DO
          END DO
+         virial%pv_mp2 = virial%pv_mp2 + h_stress/REAL(para_env%num_pe, dp)
          virial%pv_virial = virial%pv_virial + h_stress/REAL(para_env%num_pe, dp)
 
          ! free stuff

src/mp2_ri_grad.F

@@ -611,6 +611,7 @@ SUBROUTINE calc_ri_mp2_nonsep(qs_env, mp2_env, para_env, para_env_sub, cell, par
                force(ikind)%rho_elec(:, atom_a) = &
                   force(ikind)%rho_elec(:, atom_a) + force_a(:) + force_b
                IF (use_virial) THEN
+                  virial%pv_mp2 = virial%pv_mp2 + my_virial_a + my_virial_b
                   virial%pv_virial = virial%pv_virial + my_virial_a + my_virial_b
                END IF
             END DO
@@ -678,6 +679,7 @@ SUBROUTINE calc_ri_mp2_nonsep(qs_env, mp2_env, para_env, para_env_sub, cell, par
                      h_stress(beta, alpha) = h_stress(alpha, beta)
                   END DO
                END DO
+               virial%pv_mp2 = virial%pv_mp2 + h_stress/REAL(para_env_sub%num_pe, dp)
                virial%pv_virial = virial%pv_virial + h_stress/REAL(para_env_sub%num_pe, dp)
                CALL timestop(handle3)
             END IF
@@ -718,6 +720,7 @@ SUBROUTINE calc_ri_mp2_nonsep(qs_env, mp2_env, para_env, para_env_sub, cell, par
                      h_stress(beta, alpha) = h_stress(alpha, beta)
                   END DO
                END DO
+               virial%pv_mp2 = virial%pv_mp2 + h_stress/REAL(para_env_sub%num_pe, dp)
                virial%pv_virial = virial%pv_virial + h_stress/REAL(para_env_sub%num_pe, dp)
                CALL timestop(handle3)
             END IF
@@ -862,6 +865,7 @@ SUBROUTINE calc_ri_mp2_nonsep(qs_env, mp2_env, para_env, para_env_sub, cell, par
                force(ikind)%rho_elec(:, atom_a) = &
                   force(ikind)%rho_elec(:, atom_a) + force_a(:) + force_b(:)
                IF (use_virial) THEN
+                  virial%pv_mp2 = virial%pv_mp2 + my_virial_a + my_virial_b
                   virial%pv_virial = virial%pv_virial + my_virial_a + my_virial_b
                END IF
             END DO

src/qs_force.F

@@ -412,8 +412,12 @@ SUBROUTINE qs_forces(qs_env)
          CALL mp_sum(virial%pv_ecore_overlap, para_env%group)
          CALL mp_sum(virial%pv_ehartree, para_env%group)
          CALL mp_sum(virial%pv_exc, para_env%group)
-         CALL mp_sum(virial%pv_fock_4c, para_env%group)
+         CALL mp_sum(virial%pv_exx, para_env%group)
          CALL mp_sum(virial%pv_vdw, para_env%group)
+         CALL mp_sum(virial%pv_mp2, para_env%group)
+         CALL mp_sum(virial%pv_nlcc, para_env%group)
+         CALL mp_sum(virial%pv_gapw, para_env%group)
+         CALL mp_sum(virial%pv_lrigpw, para_env%group)
          CALL mp_sum(virial%pv_virial, para_env%group)
          CALL symmetrize_virial(virial)
          ! Add the volume terms of the virial
@@ -422,9 +426,8 @@ SUBROUTINE qs_forces(qs_env)
                      dft_control%qs_control%xtb .OR. &
                      dft_control%qs_control%semi_empirical))) THEN
             DO i = 1, 3
-               virial%pv_virial(i, i) = virial%pv_virial(i, i) - energy%exc - &
-                                        2.0_dp*(energy%hartree + energy%sccs_pol)
                virial%pv_ehartree(i, i) = virial%pv_ehartree(i, i) - 2.0_dp*(energy%hartree + energy%sccs_pol)
+               virial%pv_virial(i, i) = virial%pv_virial(i, i) - energy%exc - 2.0_dp*(energy%hartree + energy%sccs_pol)
                virial%pv_exc(i, i) = virial%pv_exc(i, i) - energy%exc
                IF (dft_control%do_admm) THEN
                   virial%pv_exc(i, i) = virial%pv_exc(i, i) - energy%exc_aux_fit

src/qs_integrate_potential_single.F

@@ -262,8 +262,8 @@ SUBROUTINE integrate_ppl_rspace(rho_rspace, qs_env)
                force(ikind)%gth_ppl(:, iatom) + force_a(:)*rho_rspace%pw%pw_grid%dvol
 
             IF (use_virial) THEN
-               virial%pv_virial = virial%pv_virial + my_virial_a*rho_rspace%pw%pw_grid%dvol
                virial%pv_ppl = virial%pv_ppl + my_virial_a*rho_rspace%pw%pw_grid%dvol
+               virial%pv_virial = virial%pv_virial + my_virial_a*rho_rspace%pw%pw_grid%dvol
                CPABORT("Virial not debuged for CORE_PPL")
             END IF
          END DO
@@ -462,8 +462,8 @@ SUBROUTINE integrate_rho_nlcc(rho_rspace, qs_env)
                   force(ikind)%gth_nlcc(:, iatom) + force_a(:)*rho_rspace%pw%pw_grid%dvol
 
                IF (use_virial) THEN
+                  virial%pv_nlcc = virial%pv_nlcc + my_virial_a*rho_rspace%pw%pw_grid%dvol
                   virial%pv_virial = virial%pv_virial + my_virial_a*rho_rspace%pw%pw_grid%dvol
-                  virial%pv_exc = virial%pv_exc + my_virial_a*rho_rspace%pw%pw_grid%dvol
                END IF
             END DO
 
@@ -862,6 +862,7 @@ SUBROUTINE integrate_v_rspace_one_center(v_rspace, qs_env, int_res, &
             IF (calculate_forces) THEN
                int_res(ikind)%v_dfdr(iatom, :) = int_res(ikind)%v_dfdr(iatom, :) + force_a(:)
                IF (use_virial) THEN
+                  virial%pv_lrigpw = virial%pv_lrigpw + my_virial_a
                   virial%pv_virial = virial%pv_virial + my_virial_a
                ENDIF
             ENDIF
@@ -1113,6 +1114,7 @@ SUBROUTINE integrate_v_rspace_diagonal(v_rspace, ksmat, pmat, qs_env, calculate_
                force(ikind)%rho_elec(:, iatom) = force(ikind)%rho_elec(:, iatom) + 2.0_dp*force_a(:)
                IF (use_virial) THEN
                   IF (use_virial .AND. calculate_forces) THEN
+                     virial%pv_lrigpw = virial%pv_lrigpw + 2.0_dp*my_virial_a
                      virial%pv_virial = virial%pv_virial + 2.0_dp*my_virial_a
                   END IF
                END IF

src/qs_ks_atom.F

@@ -543,6 +543,7 @@ SUBROUTINE update_ks_atom(qs_env, ksmat, pmat, forces, tddft, rho_atom_external,
 !$OMP END PARALLEL
 
       IF (use_virial) THEN
+         virial%pv_gapw(1:3, 1:3) = virial%pv_gapw(1:3, 1:3) + pv_virial_thread(1:3, 1:3)
          virial%pv_virial(1:3, 1:3) = virial%pv_virial(1:3, 1:3) + pv_virial_thread(1:3, 1:3)
       END IF
 

src/qs_rho0_ggrid.F

@@ -678,6 +678,7 @@ SUBROUTINE integrate_vhg0_rspace(qs_env, v_rspace, para_env, calculate_forces, l
                   DO iso = 1, nsoset(l0_ikind)
                      DO ii = 1, 3
                      DO i = 1, 3
+                        virial%pv_gapw(i, ii) = virial%pv_gapw(i, ii) + Qlm(iso)*a_hdab_sph(i, ii, iso)
                         virial%pv_virial(i, ii) = virial%pv_virial(i, ii) + Qlm(iso)*a_hdab_sph(i, ii, iso)
                      END DO
                      END DO

src/subsys/virial_types.F

@@ -31,6 +31,7 @@ MODULE virial_types
                                         pv_kinetic, &
                                         pv_virial, &
                                         pv_xc, &
+                                        pv_fock_4c, &
                                         pv_constraint
       REAL(KIND=dp), DIMENSION(3, 3) :: pv_overlap, &
                                         pv_ekinetic, &
@@ -39,8 +40,12 @@ MODULE virial_types
                                         pv_ecore_overlap, &
                                         pv_ehartree, &
                                         pv_exc, &
-                                        pv_fock_4c, &
-                                        pv_vdw
+                                        pv_exx, &
+                                        pv_vdw, &
+                                        pv_mp2, &
+                                        pv_nlcc, &
+                                        pv_gapw, &
+                                        pv_lrigpw
       LOGICAL                        :: pv_availability, &
                                         pv_calculate, &
                                         pv_numer, &
@@ -69,6 +74,7 @@ SUBROUTINE cp_virial(virial_in, virial_out)
       virial_out%pv_kinetic = virial_in%pv_kinetic
       virial_out%pv_virial = virial_in%pv_virial
       virial_out%pv_xc = virial_in%pv_xc
+      virial_out%pv_fock_4c = virial_in%pv_fock_4c
       virial_out%pv_constraint = virial_in%pv_constraint
 
       virial_out%pv_overlap = virial_in%pv_overlap
@@ -78,8 +84,12 @@ SUBROUTINE cp_virial(virial_in, virial_out)
       virial_out%pv_ecore_overlap = virial_in%pv_ecore_overlap
       virial_out%pv_ehartree = virial_in%pv_ehartree
       virial_out%pv_exc = virial_in%pv_exc
-      virial_out%pv_fock_4c = virial_in%pv_fock_4c
+      virial_out%pv_exx = virial_in%pv_exx
       virial_out%pv_vdw = virial_in%pv_vdw
+      virial_out%pv_mp2 = virial_in%pv_mp2
+      virial_out%pv_nlcc = virial_in%pv_nlcc
+      virial_out%pv_gapw = virial_in%pv_gapw
+      virial_out%pv_lrigpw = virial_in%pv_lrigpw
 
       virial_out%pv_availability = virial_in%pv_availability
       virial_out%pv_calculate = virial_in%pv_calculate
@@ -108,6 +118,8 @@ SUBROUTINE symmetrize_virial(virial)
             virial%pv_virial(i, j) = virial%pv_virial(j, i)
             virial%pv_xc(j, i) = 0.5_dp*(virial%pv_xc(i, j) + virial%pv_xc(j, i))
             virial%pv_xc(i, j) = virial%pv_xc(j, i)
+            virial%pv_fock_4c(j, i) = 0.5_dp*(virial%pv_fock_4c(i, j) + virial%pv_fock_4c(j, i))
+            virial%pv_fock_4c(i, j) = virial%pv_fock_4c(j, i)
             virial%pv_constraint(j, i) = 0.5_dp*(virial%pv_constraint(i, j) + virial%pv_constraint(j, i))
             virial%pv_constraint(i, j) = virial%pv_constraint(j, i)
             ! Virial components
@@ -125,10 +137,18 @@ SUBROUTINE symmetrize_virial(virial)
             virial%pv_ehartree(i, j) = virial%pv_ehartree(j, i)
             virial%pv_exc(j, i) = 0.5_dp*(virial%pv_exc(i, j) + virial%pv_exc(j, i))
             virial%pv_exc(i, j) = virial%pv_exc(j, i)
-            virial%pv_fock_4c(j, i) = 0.5_dp*(virial%pv_fock_4c(i, j) + virial%pv_fock_4c(j, i))
-            virial%pv_fock_4c(i, j) = virial%pv_fock_4c(j, i)
+            virial%pv_exx(j, i) = 0.5_dp*(virial%pv_exx(i, j) + virial%pv_exx(j, i))
+            virial%pv_exx(i, j) = virial%pv_exx(j, i)
             virial%pv_vdw(j, i) = 0.5_dp*(virial%pv_vdw(i, j) + virial%pv_vdw(j, i))
             virial%pv_vdw(i, j) = virial%pv_vdw(j, i)
+            virial%pv_mp2(j, i) = 0.5_dp*(virial%pv_mp2(i, j) + virial%pv_mp2(j, i))
+            virial%pv_mp2(i, j) = virial%pv_mp2(j, i)
+            virial%pv_nlcc(j, i) = 0.5_dp*(virial%pv_nlcc(i, j) + virial%pv_nlcc(j, i))
+            virial%pv_nlcc(i, j) = virial%pv_nlcc(j, i)
+            virial%pv_gapw(j, i) = 0.5_dp*(virial%pv_gapw(i, j) + virial%pv_gapw(j, i))
+            virial%pv_gapw(i, j) = virial%pv_gapw(j, i)
+            virial%pv_lrigpw(j, i) = 0.5_dp*(virial%pv_lrigpw(i, j) + virial%pv_lrigpw(j, i))
+            virial%pv_lrigpw(i, j) = virial%pv_lrigpw(j, i)
          END DO
       END DO
 
@@ -152,6 +172,7 @@ SUBROUTINE zero_virial(virial, reset)
       virial%pv_kinetic = 0.0_dp
       virial%pv_virial = 0.0_dp
       virial%pv_xc = 0.0_dp
+      virial%pv_fock_4c = 0.0_dp
       virial%pv_constraint = 0.0_dp
 
       virial%pv_overlap = 0.0_dp
@@ -161,8 +182,12 @@ SUBROUTINE zero_virial(virial, reset)
       virial%pv_ecore_overlap = 0.0_dp
       virial%pv_ehartree = 0.0_dp
       virial%pv_exc = 0.0_dp
-      virial%pv_fock_4c = 0.0_dp
+      virial%pv_exx = 0.0_dp
       virial%pv_vdw = 0.0_dp
+      virial%pv_mp2 = 0.0_dp
+      virial%pv_nlcc = 0.0_dp
+      virial%pv_gapw = 0.0_dp
+      virial%pv_lrigpw = 0.0_dp
 
       IF (my_reset) THEN
          virial%pv_availability = .FALSE.
@@ -180,6 +205,7 @@ END SUBROUTINE zero_virial
 !> \param pv_kinetic ...
 !> \param pv_virial ...
 !> \param pv_xc ...
+!> \param pv_fock_4c ...
 !> \param pv_constraint ...
 !> \param pv_overlap ...
 !> \param pv_ekinetic ...
@@ -188,29 +214,34 @@ END SUBROUTINE zero_virial
 !> \param pv_ecore_overlap ...
 !> \param pv_ehartree ...
 !> \param pv_exc ...
-!> \param pv_fock_4c ...
+!> \param pv_exx ...
 !> \param pv_vdw ...
+!> \param pv_mp2 ...
+!> \param pv_nlcc ...
+!> \param pv_gapw ...
+!> \param pv_lrigpw ...
 !> \param pv_availability ...
 !> \param pv_calculate ...
 !> \param pv_numer ...
 !> \param pv_diagonal ...
 ! **************************************************************************************************
-   SUBROUTINE virial_set(virial, pv_total, pv_kinetic, pv_virial, pv_xc, pv_constraint, &
+   SUBROUTINE virial_set(virial, pv_total, pv_kinetic, pv_virial, pv_xc, pv_fock_4c, pv_constraint, &
                          pv_overlap, pv_ekinetic, pv_ppl, pv_ppnl, pv_ecore_overlap, pv_ehartree, &
-                         pv_exc, pv_fock_4c, pv_vdw, pv_availability, pv_calculate, pv_numer, &
-                         pv_diagonal)
+                         pv_exc, pv_exx, pv_vdw, pv_mp2, pv_nlcc, pv_gapw, pv_lrigpw, &
+                         pv_availability, pv_calculate, pv_numer, pv_diagonal)
 
       TYPE(virial_type), POINTER                         :: virial
       REAL(KIND=dp), DIMENSION(3, 3), OPTIONAL :: pv_total, pv_kinetic, pv_virial, pv_xc, &
-         pv_constraint, pv_overlap, pv_ekinetic, pv_ppl, pv_ppnl, pv_ecore_overlap, pv_ehartree, &
-         pv_exc, pv_fock_4c, pv_vdw
+         pv_fock_4c, pv_constraint, pv_overlap, pv_ekinetic, pv_ppl, pv_ppnl, pv_ecore_overlap, &
+         pv_ehartree, pv_exc, pv_exx, pv_vdw, pv_mp2, pv_nlcc, pv_gapw, pv_lrigpw
       LOGICAL, OPTIONAL                                  :: pv_availability, pv_calculate, pv_numer, &
                                                             pv_diagonal
 
       IF (PRESENT(pv_total)) virial%pv_total = pv_total
       IF (PRESENT(pv_kinetic)) virial%pv_kinetic = pv_kinetic
       IF (PRESENT(pv_virial)) virial%pv_virial = pv_virial
       IF (PRESENT(pv_xc)) virial%pv_xc = pv_xc
+      IF (PRESENT(pv_fock_4c)) virial%pv_fock_4c = pv_fock_4c
       IF (PRESENT(pv_constraint)) virial%pv_constraint = pv_constraint
 
       IF (PRESENT(pv_overlap)) virial%pv_overlap = pv_overlap
@@ -220,8 +251,12 @@ SUBROUTINE virial_set(virial, pv_total, pv_kinetic, pv_virial, pv_xc, pv_constra
       IF (PRESENT(pv_ecore_overlap)) virial%pv_ecore_overlap = pv_ecore_overlap
       IF (PRESENT(pv_ehartree)) virial%pv_ehartree = pv_ehartree
       IF (PRESENT(pv_exc)) virial%pv_exc = pv_exc
-      IF (PRESENT(pv_fock_4c)) virial%pv_fock_4c = pv_fock_4c
+      IF (PRESENT(pv_exx)) virial%pv_exx = pv_exx
       IF (PRESENT(pv_vdw)) virial%pv_vdw = pv_vdw
+      IF (PRESENT(pv_mp2)) virial%pv_mp2 = pv_mp2
+      IF (PRESENT(pv_nlcc)) virial%pv_nlcc = pv_nlcc
+      IF (PRESENT(pv_gapw)) virial%pv_gapw = pv_gapw
+      IF (PRESENT(pv_lrigpw)) virial%pv_lrigpw = pv_lrigpw
 
       IF (PRESENT(pv_availability)) virial%pv_availability = pv_availability
       IF (PRESENT(pv_calculate)) virial%pv_calculate = pv_calculate

src/virial_methods.F

@@ -197,15 +197,28 @@ SUBROUTINE write_stress_tensor_components(virial, iw, cell)
          pv = fconv*virial%pv_exc
          WRITE (UNIT=iw, FMT=fmt) &
             'STRESS| Exchange-correlation', one_third_sum_diag(pv), det_3x3(pv)
-         pv = -fconv*virial%pv_fock_4c
+         pv = fconv*virial%pv_exx
          WRITE (UNIT=iw, FMT=fmt) &
             'STRESS| Exact exchange (EXX)', one_third_sum_diag(pv), det_3x3(pv)
          pv = fconv*virial%pv_vdw
          WRITE (UNIT=iw, FMT=fmt) &
             'STRESS| vdW correction', one_third_sum_diag(pv), det_3x3(pv)
+         pv = fconv*virial%pv_mp2
+         WRITE (UNIT=iw, FMT=fmt) &
+            'STRESS| Moller-Plesset (MP2)', one_third_sum_diag(pv), det_3x3(pv)
+         pv = fconv*virial%pv_nlcc
+         WRITE (UNIT=iw, FMT=fmt) &
+            'STRESS| Nonlinear core correction (NLCC)', one_third_sum_diag(pv), det_3x3(pv)
+         pv = fconv*virial%pv_gapw
+         WRITE (UNIT=iw, FMT=fmt) &
+            'STRESS| Local atomic parts (GAPW)', one_third_sum_diag(pv), det_3x3(pv)
+         pv = fconv*virial%pv_lrigpw
+         WRITE (UNIT=iw, FMT=fmt) &
+            'STRESS| Resolution-of-the-identity (LRI)', one_third_sum_diag(pv), det_3x3(pv)
          pv = fconv*(virial%pv_overlap + virial%pv_ekinetic + virial%pv_ppl + virial%pv_ppnl + &
-                     virial%pv_ecore_overlap + virial%pv_ehartree + virial%pv_exc - &
-                     virial%pv_fock_4c + virial%pv_vdw)
+                     virial%pv_ecore_overlap + virial%pv_ehartree + virial%pv_exc + &
+                     virial%pv_exx + virial%pv_vdw + virial%pv_mp2 + virial%pv_nlcc + &
+                     virial%pv_gapw + virial%pv_lrigpw)
          WRITE (UNIT=iw, FMT=fmt) &
             'STRESS| Sum of components', one_third_sum_diag(pv), det_3x3(pv)
          pv = fconv*virial%pv_virial