Fix unintialized virial_thread

src/qs_dispersion_nonloc.F

@@ -602,12 +602,8 @@ SUBROUTINE vdW_energy(thetas_g, dispersion_env, vdW_xc_energy, energy_only, viri
       CALL timeset(routineN, handle)
       nqs = dispersion_env%nqs
 
-      IF (PRESENT(virial)) THEN
-         use_virial = .TRUE.
-         virial_thread(:, :) = 0.0_dp
-      ELSE
-         use_virial = .FALSE.
-      END IF
+      use_virial = PRESENT(virial)
+      virial_thread(:, :) = 0.0_dp ! always initialize to avoid floating point exceptions in OMP REDUCTION
 
       vdW_xc_energy = 0._dp
       grid => thetas_g(1)%pw%pw_grid
@@ -769,15 +765,11 @@ SUBROUTINE get_potential(q0, dq0_drho, dq0_dgradrho, total_rho, u_vdW, potential
 
       CALL timeset(routineN, handle)
 
-      IF (PRESENT(virial)) THEN
-         use_virial = .TRUE.
-         virial_thread(:, :) = 0.0_dp
-         CPASSERT(PRESENT(drho))
-         CPASSERT(PRESENT(dvol))
-      ELSE
-         use_virial = .FALSE.
-      END IF
+      use_virial = PRESENT(virial)
+      CPASSERT(.NOT. use_virial .OR. PRESENT(drho))
+      CPASSERT(.NOT. use_virial .OR. PRESENT(dvol))
 
+      virial_thread(:, :) = 0.0_dp ! always initialize to avoid floating point exceptions in OMP REDUCTION
       b_value = dispersion_env%b_value
       const = 1.0_dp/(3.0_dp*b_value**(3.0_dp/2.0_dp)*pi**(5.0_dp/4.0_dp))
       potential = 0.0_dp

src/qs_ks_atom.F

@@ -213,11 +213,12 @@ SUBROUTINE update_ks_atom(qs_env, ksmat, pmat, forces, tddft, rho_atom_external,
          CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind)
          ldCPC = max_gau
          use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
-         IF (use_virial) pv_virial_thread(:, :) = 0.0_dp
       ELSE
          use_virial = .FALSE.
       END IF
 
+      pv_virial_thread(:, :) = 0.0_dp ! always initialize to avoid floating point exceptions in OMP REDUCTION
+
       nkind = SIZE(my_kind_set, 1)
       ! Collect the local potential contributions from all the processors
       mepos = para_env%mepos