Fix computation of stress tensor components

cp2k · Oct 19, 2020 · 9ca2177 · 9ca2177
1 parent 965dda1
commit 9ca2177
Show file tree

Hide file tree

Showing 14 changed files with 160 additions and 154 deletions.
diff --git a/src/core_ae.F b/src/core_ae.F
@@ -112,7 +112,7 @@ SUBROUTINE build_core_ae(matrix_h, matrix_p, force, virial, calculate_forces, us
       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: habd, work
       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: hab, pab, verf, vnuc
       REAL(KIND=dp), DIMENSION(3)                        :: force_a, force_b, rab, rac, rbc
-      REAL(KIND=dp), DIMENSION(3, 3)                     :: pv_loc, pv_thread
+      REAL(KIND=dp), DIMENSION(3, 3)                     :: pv_thread
       TYPE(neighbor_list_iterator_p_type), &
          DIMENSION(:), POINTER                           :: ap_iterator
       TYPE(gto_basis_set_type), POINTER                  :: basis_set_a, basis_set_b
@@ -153,11 +153,8 @@ SUBROUTINE build_core_ae(matrix_h, matrix_p, force, virial, calculate_forces, us
             END DO
          END IF
       END IF
-      force_thread = 0.0_dp
 
-      IF (use_virial) THEN
-         pv_loc = virial%pv_virial
-      END IF
+      force_thread = 0.0_dp
       pv_thread = 0.0_dp
 
       ALLOCATE (basis_set_list(nkind))
@@ -479,8 +476,8 @@ SUBROUTINE build_core_ae(matrix_h, matrix_p, force, virial, calculate_forces, us
       END IF
 
       IF (calculate_forces .AND. use_virial) THEN
+         virial%pv_ppl = virial%pv_ppl + pv_thread
          virial%pv_virial = virial%pv_virial + pv_thread
-         virial%pv_ppl = virial%pv_virial - pv_loc
       END IF
 
       CALL timestop(handle)

diff --git a/src/core_ppl.F b/src/core_ppl.F
@@ -123,7 +123,7 @@ SUBROUTINE build_core_ppl(matrix_h, matrix_p, force, virial, calculate_forces, u
       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :)  :: hab, pab
       REAL(KIND=dp), DIMENSION(1:10)                     :: aloc, bloc
       REAL(KIND=dp), DIMENSION(3)                        :: force_a, force_b, rab, rac, rbc
-      REAL(KIND=dp), DIMENSION(3, 3)                     :: pv_loc, pv_thread
+      REAL(KIND=dp), DIMENSION(3, 3)                     :: pv_thread
       TYPE(neighbor_list_iterator_p_type), &
          DIMENSION(:), POINTER                           :: ap_iterator
       TYPE(gto_basis_set_type), POINTER                  :: basis_set_a, basis_set_b
@@ -191,9 +191,6 @@ SUBROUTINE build_core_ppl(matrix_h, matrix_p, force, virial, calculate_forces, u
          END IF
       END DO
 
-      IF (use_virial) THEN
-         pv_loc = virial%pv_virial
-      END IF
       pv_thread = 0.0_dp
 
       nthread = 1
@@ -537,8 +534,8 @@ SUBROUTINE build_core_ppl(matrix_h, matrix_p, force, virial, calculate_forces, u
       END IF
 
       IF (calculate_forces .AND. use_virial) THEN
+         virial%pv_ppl = virial%pv_ppl + pv_thread
          virial%pv_virial = virial%pv_virial + pv_thread
-         virial%pv_ppl = virial%pv_virial - pv_loc
       END IF
 
       CALL timestop(handle)
@@ -800,8 +797,8 @@ SUBROUTINE build_core_ppl_ri(lri_ppl_coef, force, virial, calculate_forces, use_
       DEALLOCATE (atom_of_kind, kind_of)
 
       IF (calculate_forces .AND. use_virial) THEN
+         virial%pv_ppl = virial%pv_ppl + pv_thread
          virial%pv_virial = virial%pv_virial + pv_thread
-         virial%pv_ppl = pv_thread
       END IF
 
       DEALLOCATE (basis_set_list)

diff --git a/src/core_ppnl.F b/src/core_ppnl.F
@@ -118,7 +118,7 @@ SUBROUTINE build_core_ppnl(matrix_h, matrix_p, force, virial, calculate_forces,
       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: ai_work
       REAL(KIND=dp), DIMENSION(1)                        :: rprjc, zetc
       REAL(KIND=dp), DIMENSION(3)                        :: fa, fb, rab, rac, rbc
-      REAL(KIND=dp), DIMENSION(3, 3)                     :: pv_loc, pv_thread
+      REAL(KIND=dp), DIMENSION(3, 3)                     :: pv_thread
       TYPE(gto_basis_set_type), POINTER                  :: orb_basis_set
       TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER  :: basis_set
       TYPE(gth_potential_type), POINTER                  :: gth_potential
@@ -168,10 +168,6 @@ SUBROUTINE build_core_ppnl(matrix_h, matrix_p, force, virial, calculate_forces,
             END IF
          END IF
 
-         IF (use_virial) THEN
-            pv_loc = virial%pv_virial
-         END IF
-
          maxder = ncoset(nder)
 
          CALL get_qs_kind_set(qs_kind_set, &
@@ -428,6 +424,7 @@ SUBROUTINE build_core_ppnl(matrix_h, matrix_p, force, virial, calculate_forces,
          CALL sap_sort(sap_int)
          ! *** All integrals needed have been calculated and stored in sap_int
          ! *** We now calculate the Hamiltonian matrix elements
+
          force_thread = 0.0_dp
          pv_thread = 0.0_dp
 
@@ -598,8 +595,8 @@ SUBROUTINE build_core_ppnl(matrix_h, matrix_p, force, virial, calculate_forces,
          END IF
 
          IF (calculate_forces .AND. use_virial) THEN
+            virial%pv_ppnl = virial%pv_ppnl + pv_thread
             virial%pv_virial = virial%pv_virial + pv_thread
-            virial%pv_ppnl = virial%pv_virial - pv_loc
          END IF
 
          IF (calculate_forces) THEN

diff --git a/src/cp2k_debug.F b/src/cp2k_debug.F
@@ -276,7 +276,7 @@ SUBROUTINE cp2k_debug_energy_and_forces(force_env)
                particles(ip)%r(k) = std_value
                numer_forces(ip, k) = -0.5_dp*(numer_energy(1) - numer_energy(2))/dx
                IF (iw > 0) THEN
-                  WRITE (UNIT=iw, FMT="(T2,A,T17,A,F7.4,A,T34,A,F7.4,A,T52,A,T68,A)") &
+                  WRITE (UNIT=iw, FMT="(/,T2,A,T17,A,F7.4,A,T34,A,F7.4,A,T52,A,T68,A)") &
                      "DEBUG| Atom", "E("//ACHAR(119 + k)//" +", dx, ")", &
                      "E("//ACHAR(119 + k)//" -", dx, ")", &
                      "f(numerical)", "f(analytical)"

diff --git a/src/force_env_methods.F b/src/force_env_methods.F
@@ -165,7 +165,7 @@ MODULE force_env_methods
    USE virial_methods,                  ONLY: write_stress_tensor,&
                                               write_stress_tensor_components
    USE virial_types,                    ONLY: cp_virial,&
-                                              sym_virial,&
+                                              symmetrize_virial,&
                                               virial_create,&
                                               virial_p_type,&
                                               virial_release,&
@@ -301,7 +301,7 @@ RECURSIVE SUBROUTINE force_env_calc_energy_force(force_env, calc_force, &
          ELSE
             IF (calculate_forces) THEN
                ! Symmetrize analytical stress tensor
-               CALL sym_virial(virial)
+               CALL symmetrize_virial(virial)
             ELSE
                IF (calculate_stress_tensor) THEN
                   CALL cp_warn(__LOCATION__, "The calculation of the stress tensor "// &

diff --git a/src/qs_core_energies.F b/src/qs_core_energies.F
@@ -320,8 +320,8 @@ SUBROUTINE calculate_ecore_overlap(qs_env, para_env, calculate_forces, molecular
          DEALLOCATE (atom_of_kind)
       END IF
       IF (calculate_forces .AND. use_virial) THEN
+         virial%pv_ecore_overlap = virial%pv_ecore_overlap + pv_loc
          virial%pv_virial = virial%pv_virial + pv_loc
-         virial%pv_hartree = virial%pv_hartree + pv_loc
       END IF
 
       CALL mp_sum(ecore_overlap, group)

diff --git a/src/qs_dispersion_pairpot.F b/src/qs_dispersion_pairpot.F
@@ -2104,7 +2104,7 @@ SUBROUTINE calculate_dispersion_pairpot(qs_env, dispersion_env, energy, calculat
       END IF
 
       IF (calculate_forces .AND. use_virial) THEN
-         virial%pv_vdw = virial%pv_virial - pv_loc
+         virial%pv_vdw = virial%pv_vdw + (virial%pv_virial - pv_loc)
       END IF
 
       IF (ASSOCIATED(dispersion_env%dftd_section)) THEN

diff --git a/src/qs_force.F b/src/qs_force.F
@@ -78,7 +78,8 @@ MODULE qs_force
                                               rt_admm_force
    USE se_core_core,                    ONLY: se_core_core_interaction
    USE se_core_matrix,                  ONLY: build_se_core_matrix
-   USE virial_types,                    ONLY: virial_type
+   USE virial_types,                    ONLY: symmetrize_virial,&
+                                              virial_type
    USE xtb_matrices,                    ONLY: build_xtb_matrices
 #include "./base/base_uses.f90"
 
@@ -131,7 +132,7 @@ SUBROUTINE qs_forces(qs_env)
 
       CHARACTER(len=*), PARAMETER                        :: routineN = 'qs_forces'
 
-      INTEGER                                            :: after, dir, handle, i, iatom, ic, ikind, &
+      INTEGER                                            :: after, handle, i, iatom, ic, ikind, &
                                                             ispin, iw, natom, nkind, nspin, &
                                                             output_unit
       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atom_of_kind, kind_of, natom_of_kind
@@ -403,23 +404,35 @@ SUBROUTINE qs_forces(qs_env)
 
       NULLIFY (virial, energy)
       CALL get_qs_env(qs_env=qs_env, virial=virial, energy=energy)
-      !   *** distribute virial ***
       IF (virial%pv_availability) THEN
+         CALL mp_sum(virial%pv_overlap, para_env%group)
+         CALL mp_sum(virial%pv_ekinetic, para_env%group)
+         CALL mp_sum(virial%pv_ppl, para_env%group)
+         CALL mp_sum(virial%pv_ppnl, para_env%group)
+         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_vdw, para_env%group)
          CALL mp_sum(virial%pv_virial, para_env%group)
-         !  *** add the volume terms of the virial ***
+         CALL symmetrize_virial(virial)
+         ! Add the volume terms of the virial
          IF ((.NOT. virial%pv_numer) .AND. &
              (.NOT. (dft_control%qs_control%dftb .OR. &
                      dft_control%qs_control%xtb .OR. &
                      dft_control%qs_control%semi_empirical))) THEN
-            DO dir = 1, 3
-               virial%pv_virial(dir, dir) = virial%pv_virial(dir, dir) - energy%exc &
-                                            - 2.0_dp*(energy%hartree + energy%sccs_pol)
+            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_exc(i, i) = virial%pv_exc(i, i) - energy%exc
                IF (dft_control%do_admm) THEN
-                  virial%pv_virial(dir, dir) = virial%pv_virial(dir, dir) - energy%exc_aux_fit
+                  virial%pv_exc(i, i) = virial%pv_exc(i, i) - energy%exc_aux_fit
+                  virial%pv_virial(i, i) = virial%pv_virial(i, i) - energy%exc_aux_fit
                END IF
                ! The factor 2 is a hack. It compensates the plus sign in h_stress/pw_poisson_solve.
                ! The sign in pw_poisson_solve is correct for FIST, but not for QS.
-               ! There should be a more elegant solution to that...
+               ! There should be a more elegant solution to that ...
             END DO
          END IF
       END IF

diff --git a/src/qs_integrate_potential_single.F b/src/qs_integrate_potential_single.F
@@ -619,8 +619,8 @@ SUBROUTINE integrate_v_core_rspace(v_rspace, qs_env)
                   force(ikind)%rho_core(:, iatom) = force(ikind)%rho_core(:, iatom) + force_a(:)
                END IF
                IF (use_virial) THEN
+                  virial%pv_ehartree = virial%pv_ehartree + my_virial_a
                   virial%pv_virial = virial%pv_virial + my_virial_a
-                  virial%pv_hartree = virial%pv_hartree + my_virial_a
                END IF
                IF (ASSOCIATED(atprop)) THEN
                   atprop%ateb(atom_a) = atprop%ateb(atom_a) + 0.5_dp*hab(1, 1)*pab(1, 1)

diff --git a/src/qs_kinetic.F b/src/qs_kinetic.F
@@ -166,7 +166,7 @@ SUBROUTINE build_kinetic_matrix_low(ks_env, matrix_t, matrixkp_t, matrix_name, b
       REAL(KIND=dp)                                      :: f0, ff, rab2, tab
       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: kab, pab, qab
       REAL(KIND=dp), DIMENSION(3)                        :: force_a, rab
-      REAL(KIND=dp), DIMENSION(3, 3)                     :: pv_loc, pv_thread
+      REAL(KIND=dp), DIMENSION(3, 3)                     :: pv_thread
       REAL(KIND=dp), DIMENSION(3, natom)                 :: force_thread
       REAL(KIND=dp), DIMENSION(:), POINTER               :: set_radius_a, set_radius_b
       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: k_block, p_block, rpgfa, rpgfb, scon_a, &
@@ -238,14 +238,14 @@ SUBROUTINE build_kinetic_matrix_low(ks_env, matrix_t, matrixkp_t, matrix_name, b
          ! if forces -> maybe virial too
          CALL get_ks_env(ks_env=ks_env, force=force, virial=virial)
          use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
-         pv_loc = virial%pv_virial
          ! we need density matrix for forces
          IF (dokp) THEN
             CPASSERT(PRESENT(matrixkp_p))
          ELSE
             CPASSERT(PRESENT(matrix_p))
          END IF
       END IF
+
       force_thread = 0.0_dp
       pv_thread = 0.0_dp
 
@@ -448,8 +448,8 @@ SUBROUTINE build_kinetic_matrix_low(ks_env, matrix_t, matrixkp_t, matrix_name, b
          DEALLOCATE (atom_of_kind, kind_of)
       END IF
       IF (do_forces .AND. use_virial) THEN
+         virial%pv_ekinetic = virial%pv_ekinetic + pv_thread
          virial%pv_virial = virial%pv_virial + pv_thread
-         virial%pv_ekin = virial%pv_virial - pv_loc
       END IF
 
       IF (dokp) THEN

diff --git a/src/qs_ks_methods.F b/src/qs_ks_methods.F
@@ -215,8 +215,8 @@ SUBROUTINE qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces, just_energy, &
       CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_ks_build_kohn_sham_matrix'
 
       CHARACTER(len=default_string_length)               :: name
-      INTEGER                                            :: handle, iatom, idir, img, ispin, &
-                                                            nimages, ns, nspins
+      INTEGER                                            :: handle, iatom, img, ispin, nimages, ns, &
+                                                            nspins
       LOGICAL :: do_adiabatic_rescaling, do_ddapc, do_hfx, do_ppl, gapw, gapw_xc, &
          hfx_treat_lsd_in_core, just_energy_xc, lrigpw, my_print, rigpw, use_virial
       REAL(KIND=dp)                                      :: ecore_ppl, edisp, ee_ener, ekin_mol, &
@@ -419,8 +419,8 @@ SUBROUTINE qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces, just_energy, &
          IF (use_virial .AND. calculate_forces) THEN
             CALL sccs(qs_env, rho_tot_gspace%pw, v_hartree_gspace%pw, v_sccs_rspace%pw, &
                       h_stress=h_stress)
+            virial%pv_ehartree = virial%pv_ehartree + h_stress/REAL(para_env%num_pe, dp)
             virial%pv_virial = virial%pv_virial + h_stress/REAL(para_env%num_pe, dp)
-            virial%pv_hartree = virial%pv_hartree + h_stress/REAL(para_env%num_pe, dp)
          ELSE
             CALL sccs(qs_env, rho_tot_gspace%pw, v_hartree_gspace%pw, v_sccs_rspace%pw)
          END IF
@@ -432,8 +432,8 @@ SUBROUTINE qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces, just_energy, &
             CALL pw_poisson_solve(poisson_env, rho_tot_gspace%pw, energy%hartree, &
                                   v_hartree_gspace%pw, h_stress=h_stress, &
                                   rho_core=rho_core)
+            virial%pv_ehartree = virial%pv_ehartree + h_stress/REAL(para_env%num_pe, dp)
             virial%pv_virial = virial%pv_virial + h_stress/REAL(para_env%num_pe, dp)
-            virial%pv_hartree = virial%pv_hartree + h_stress/REAL(para_env%num_pe, dp)
          ELSE
             CALL pw_poisson_solve(poisson_env, rho_tot_gspace%pw, energy%hartree, &
                                   v_hartree_gspace%pw, rho_core=rho_core)
@@ -577,8 +577,9 @@ SUBROUTINE qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces, just_energy, &
          NULLIFY (rho_struct)
 
          IF (use_virial .AND. calculate_forces) THEN
+            virial%pv_exc = virial%pv_exc - virial%pv_xc
             virial%pv_virial = virial%pv_virial - virial%pv_xc
-            ! ** virial%pv_xc will be zeroed in the xc routines
+            ! virial%pv_xc will be zeroed in the xc routines
          END IF
          xc_section => admm_env%xc_section_primary
       ELSE
@@ -628,18 +629,9 @@ SUBROUTINE qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces, just_energy, &
 
       NULLIFY (rho_struct)
       IF (use_virial .AND. calculate_forces) THEN
+         virial%pv_exc = virial%pv_exc - virial%pv_xc
          virial%pv_virial = virial%pv_virial - virial%pv_xc
-         virial%pv_exc = -virial%pv_xc
-         DO idir = 1, 3
-            virial%pv_exc(idir, idir) = virial%pv_exc(idir, idir) - energy%exc
-            virial%pv_hartree(idir, idir) = virial%pv_hartree(idir, idir) - 2.0_dp*energy%hartree
-         END DO
-         IF (dft_control%do_admm) THEN
-            DO idir = 1, 3
-               virial%pv_exc(idir, idir) = virial%pv_exc(idir, idir) - energy%exc_aux_fit
-            END DO
-         END IF
-      ENDIF
+      END IF
 
       ! *** Add Hartree-Fock contribution if required ***
       IF (do_hfx) THEN
@@ -714,7 +706,7 @@ SUBROUTINE qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces, just_energy, &
                                    cdft_control, calculate_forces)
 
          IF (use_virial .AND. calculate_forces) THEN
-            virial%pv_hartree = virial%pv_hartree + (virial%pv_virial - pv_loc)
+            virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
          END IF
          IF (dft_control%qs_control%do_kg) THEN
             CPASSERT(.NOT. (gapw .OR. gapw_xc))
@@ -725,6 +717,7 @@ SUBROUTINE qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces, just_energy, &
             energy%exc = energy%exc - ekin_mol
             ! virial corrections
             IF (use_virial .AND. calculate_forces) THEN
+               virial%pv_exc = virial%pv_exc + virial%pv_xc
                virial%pv_virial = virial%pv_virial + virial%pv_xc
                virial%pv_xc = 0.0_dp
             END IF

diff --git a/src/qs_overlap.F b/src/qs_overlap.F
@@ -197,7 +197,7 @@ SUBROUTINE build_overlap_matrix_low(ks_env, matrix_s, matrixkp_s, matrix_name, n
       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: owork, pmat
       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: oint
       REAL(KIND=dp), DIMENSION(3)                        :: force_a, rab
-      REAL(KIND=dp), DIMENSION(3, 3)                     :: pv_loc, pv_thread
+      REAL(KIND=dp), DIMENSION(3, 3)                     :: pv_thread
       REAL(KIND=dp), DIMENSION(3, natom)                 :: force_thread
       REAL(KIND=dp), DIMENSION(:), POINTER               :: set_radius_a, set_radius_b
       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: p_block, rpgfa, rpgfb, scon_a, scon_b, &
@@ -282,8 +282,8 @@ SUBROUTINE build_overlap_matrix_low(ks_env, matrix_s, matrixkp_s, matrix_name, n
       IF (do_forces) THEN
          CALL get_ks_env(ks_env=ks_env, force=force, virial=virial)
          use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
-         pv_loc = virial%pv_virial
       END IF
+
       force_thread = 0.0_dp
       pv_thread = 0.0_dp
 
@@ -516,8 +516,8 @@ SUBROUTINE build_overlap_matrix_low(ks_env, matrix_s, matrixkp_s, matrix_name, n
          DEALLOCATE (atom_of_kind, kind_of)
       END IF
       IF (do_forces .AND. use_virial) THEN
+         virial%pv_overlap = virial%pv_overlap + pv_thread
          virial%pv_virial = virial%pv_virial + pv_thread
-         virial%pv_overlap = virial%pv_virial - pv_loc
       ENDIF
 
       IF (dokp) THEN