Named critical sections, adjusted collocate/integrate, fixed typos.

src/hfx_derivatives.F

@@ -773,7 +773,7 @@ CONTAINS
 !$OMP BARRIER
 !$OMP MASTER
       IF (do_dynamic_load_balancing) THEN
-         ! ** Lets contstruct the task list
+         ! ** Lets construct the task list
          shm_total_bins = 0
          DO i = 1, n_threads
             shm_total_bins = shm_total_bins+SIZE(qs_env%x_data(irep, i)%distribution_forces)
@@ -866,7 +866,7 @@ CONTAINS
                distribution_forces => actual_x_data%distribution_forces(bin)
             END IF
          ELSE
-!$OMP CRITICAL
+!$OMP CRITICAL(hfxderivatives_critical)
             shm_task_counter = shm_task_counter+1
             IF (shm_task_counter <= shm_total_bins) THEN
                my_thread_id = shm_task_list(shm_task_counter)%thread_id
@@ -888,7 +888,7 @@ CONTAINS
                do_it = .FALSE.
                bins_left = .FALSE.
             END IF
-!$OMP END CRITICAL
+!$OMP END CRITICAL(hfxderivatives_critical)
          END IF
 
          IF (.NOT. do_it) CYCLE

src/hfx_energy_potential.F

@@ -975,7 +975,7 @@ CONTAINS
       END IF
 
       !! Since the I/O routines are no thread-safe, i.e. the procedure to get the unit number, put a lock here
-!$OMP CRITICAL
+!$OMP CRITICAL(hfxenergy_io_critical)
       !! If we do disk storage, we have to initialize the containers/caches
       IF (do_disk_storage) THEN
          !! IF new md step -> reinitialize containers
@@ -995,13 +995,13 @@ CONTAINS
             END DO
          END IF
       END IF
-!$OMP END CRITICAL
+!$OMP END CRITICAL(hfxenergy_io_critical)
 
 !$OMP BARRIER
 !$OMP MASTER
 
       IF (do_dynamic_load_balancing) THEN
-         ! ** Lets contstruct the task list
+         ! ** Lets construct the task list
          shm_total_bins = 0
          DO i = 1, n_threads
             shm_total_bins = shm_total_bins+SIZE(x_data(irep, i)%distribution_energy)
@@ -1113,7 +1113,7 @@ CONTAINS
                distribution_energy => actual_x_data%distribution_energy(bin)
             END IF
          ELSE
-!$OMP CRITICAL
+!$OMP CRITICAL(hfxenergy_critical)
             shm_task_counter = shm_task_counter+1
             IF (shm_task_counter <= shm_total_bins) THEN
                my_thread_id = shm_task_list(shm_task_counter)%thread_id
@@ -1135,7 +1135,7 @@ CONTAINS
                do_it = .FALSE.
                bins_left = .FALSE.
             END IF
-!$OMP END CRITICAL
+!$OMP END CRITICAL(hfxenergy_critical)
          END IF
 
          IF (.NOT. do_it) CYCLE
@@ -2121,7 +2121,7 @@ CONTAINS
       END IF
 
       !! Since the I/O routines are no thread-safe, i.e. the procedure to get the unit number, put a lock here
-!$OMP CRITICAL
+!$OMP CRITICAL(hfxenergy_out_critical)
       IF (do_disk_storage) THEN
          !! flush caches if the geometry changed
          IF (my_geo_change) THEN
@@ -2140,7 +2140,7 @@ CONTAINS
                                                memory_parameter%actual_memory_usage_disk, do_disk_storage)
          END DO
       END IF
-!$OMP END CRITICAL
+!$OMP END CRITICAL(hfxenergy_out_critical)
 !$OMP BARRIER
       !! Clean up
       DEALLOCATE (last_sgf_global)

src/motion_utils.F

@@ -219,7 +219,7 @@ CONTAINS
       dof = dof+COUNT(lrot == 1)
       iw = cp_print_key_unit_nr(logger, print_section, "ROTATIONAL_INFO", extension=".vibLog")
       IF (iw > 0) THEN
-         WRITE (iw, '(T2,A,I6)') "ROT| Numer of Rotovibrational vectors:", dof
+         WRITE (iw, '(T2,A,I6)') "ROT| Number of Rotovibrational vectors:", dof
          IF (dof == 5) WRITE (iw, '(T2,A)') "ROT| Linear Molecule detected.."
          IF (dof == 3 .AND. (.NOT. keep_rotations)) WRITE (iw, '(T2,A)') "ROT| Single Atom detected.."
       END IF

src/qs_collocate_density.F

@@ -1489,7 +1489,7 @@ CONTAINS
       IF (PRESENT(soft_valid)) my_soft = soft_valid
 
       ! by default, do not compute the kinetic energy density (tau)
-      ! if compute_tau, all grids referening to rho are actually tau
+      ! if compute_tau, all grids referencing to rho are actually tau
       IF (PRESENT(compute_tau)) THEN
          my_compute_tau = compute_tau
       ELSE
@@ -3642,7 +3642,7 @@ CONTAINS
 
       INTEGER                                            :: i, index_max(3), index_min(3), &
                                                             ipoint(3), j, k
-      REAL(KIND=dp)                                      :: point(3)
+      REAL(KIND=dp)                                      :: inv_ng(3), point(3), primpt
 
 ! still hard-wired (see MODULO)
 
@@ -3650,12 +3650,16 @@ CONTAINS
 
          CALL return_cube_nonortho(cube_info, radius, index_min, index_max, rp)
 
+         inv_ng = 1.0_dp/ng
+
          ! go over the grid, but cycle if the point is not within the radius
          DO k = index_min(3), index_max(3)
          DO j = index_min(2), index_max(2)
          DO i = index_min(1), index_max(1)
             ! point in real space
-            point = MATMUL(cell%hmat, REAL((/i, j, k/), KIND=dp)/ng)
+            point = MATMUL(cell%hmat, REAL((/i, j, k/), KIND=dp)*inv_ng)
+            ! primitive_value of point
+            primpt = primitive_value(point)
             ! skip if outside of the sphere
             IF (SUM((point-rp)**2) > radius**2) CYCLE
             ! point on the grid (including pbc)
@@ -3663,9 +3667,9 @@ CONTAINS
             ! add to grid
             IF (PRESENT(lgrid)) THEN
                ig = ipoint(3)*ng(2)*ng(1)+ipoint(2)*ng(1)+ipoint(1)+1
-               lgrid%r(ig, ithread_l) = lgrid%r(ig, ithread_l)+primitive_value(point)
+               lgrid%r(ig, ithread_l) = lgrid%r(ig, ithread_l)+primpt
             ELSE
-               grid(ipoint(1), ipoint(2), ipoint(3)) = grid(ipoint(1), ipoint(2), ipoint(3))+primitive_value(point)
+               grid(ipoint(1), ipoint(2), ipoint(3)) = grid(ipoint(1), ipoint(2), ipoint(3))+primpt
             ENDIF
          ENDDO
          ENDDO
@@ -3681,7 +3685,8 @@ CONTAINS
       FUNCTION primitive_value(point) RESULT(res)
       REAL(KIND=dp)                                      :: point(3), res
 
-      REAL(KIND=dp)                                      :: dra(3), drap(3), drb(3), drbp(3), myexp
+      REAL(KIND=dp)                                      :: dra(3), drap(3), drb(3), drbp(3), myexp, &
+                                                            pdrap
 
          res = 0.0_dp
          myexp = EXP(-zetp*SUM((point-rp)**2))*prefactor
@@ -3704,10 +3709,11 @@ CONTAINS
                         DO lzb = 1, MAX(lb_min_local-lxb-lyb, 0)
                            drbp(3) = drbp(3)*drb(3)
                         ENDDO
+                        ico = coset(lxa, lya, lza)
+                        pdrap = PRODUCT(drap)
                         DO lzb = MAX(lb_min_local-lxb-lyb, 0), lb_max_local-lxb-lyb
-                           ico = coset(lxa, lya, lza)
                            jco = coset(lxb, lyb, lzb)
-                           res = res+pab_local(ico+o1_local, jco+o2_local)*myexp*PRODUCT(drap)*PRODUCT(drbp)
+                           res = res+pab_local(ico+o1_local, jco+o2_local)*myexp*pdrap*PRODUCT(drbp)
                            drbp(3) = drbp(3)*drb(3)
                         ENDDO
                         drap(3) = drap(3)*dra(3)
@@ -3726,7 +3732,7 @@ CONTAINS
    END SUBROUTINE collocate_pgf_product_rspace
 
 ! **************************************************************************************************
-!> \brief low level collcation of primitive gaussian functions in g-space
+!> \brief low level collocation of primitive gaussian functions in g-space
 !> \param la_max ...
 !> \param zeta ...
 !> \param la_min ...

src/qs_integrate_potential_product.F

@@ -617,7 +617,7 @@ CONTAINS
             rab(1) = dist_ab(1, itask)
             rab(2) = dist_ab(2, itask)
             rab(3) = dist_ab(3, itask)
-            rab2 = rab(1)*rab(1)+rab(2)*rab(2)+rab(3)*rab(3)
+            rab2 = DOT_PRODUCT(rab, rab)
             rb(1) = ra(1)+rab(1)
             rb(2) = ra(2)+rab(2)
             rb(3) = ra(3)+rab(3)