Parallel Jacobi Rotations for IAO bonds (#2714)

src/iao_analysis.F

@@ -1052,13 +1052,15 @@ SUBROUTINE print_center_spread(moments, nocc, print_section)
       WRITE (filename, '(A18,I1.1)') "IBO_CENTERS_SPREAD"
       iw = cp_print_key_unit_nr(logger, print_section, "", extension=".csp", &
                                 middle_name=TRIM(filename), file_position="REWIND", log_filename=.FALSE.)
-      DO ispin = 1, nspin
-         WRITE (iw, "(A,i1)") "Intrinsic Bond Orbitals: Centers and Spread for Spin ", ispin
-         WRITE (iw, "(A7,T30,A6,T68,A7)") "State", "Center", "Spreads"
-         DO is = 1, nocc(ispin)
-            WRITE (iw, "(i7,3F15.8,8X,2F10.5)") is, moments(1:5, is, ispin)
+      IF (iw > 0) THEN
+         DO ispin = 1, nspin
+            WRITE (iw, "(A,i1)") "Intrinsic Bond Orbitals: Centers and Spread for Spin ", ispin
+            WRITE (iw, "(A7,T30,A6,T68,A7)") "State", "Center", "Spreads"
+            DO is = 1, nocc(ispin)
+               WRITE (iw, "(i7,3F15.8,8X,2F10.5)") is, moments(1:5, is, ispin)
+            END DO
          END DO
-      END DO
+      END IF
       CALL cp_print_key_finished_output(iw, logger, print_section, "")
 
    END SUBROUTINE print_center_spread

src/qs_loc_methods.F

@@ -473,7 +473,6 @@ SUBROUTINE jacobi_rotation_pipek(zij_fm_set, vectors, sweeps)
                   CALL cp_fm_get_element(zij_fm_set(iatom, 1), jstate, jstate, mjj)
                   aij = aij + mij*(mii - mjj)
                   bij = bij + mij*mij - 0.25_dp*(mii - mjj)*(mii - mjj)
-
                END DO
                IF (ABS(bij) > 1.E-10_dp) THEN
                   ratio = -aij/bij
@@ -525,14 +524,19 @@ SUBROUTINE rotate_rmat_real(istate, jstate, st, ct, rmat)
 
       INTEGER                                  :: nstate
 #if defined(__SCALAPACK)
+      INTEGER                                  :: lwork, info
       INTEGER, DIMENSION(9)                    :: desc
+      REAL(dp), DIMENSION(:), ALLOCATABLE      :: work
 #endif
 
       CALL cp_fm_get_info(rmat, nrow_global=nstate)
 #if defined(__SCALAPACK)
+      lwork = 2*nstate + 1
+      ALLOCATE (work(lwork))
       desc(:) = rmat%matrix_struct%descriptor(:)
-      CALL pzrot(nstate, rmat%local_data(1, 1), 1, istate, desc, 1, &
-                 rmat%local_data(1, 1), 1, jstate, desc, 1, ct, st)
+      CALL pdrot(nstate, rmat%local_data(1, 1), 1, istate, desc, 1, &
+                 rmat%local_data(1, 1), 1, jstate, desc, 1, ct, st, work, lwork, info)
+      DEALLOCATE (work)
 #else
       CALL drot(nstate, rmat%local_data(1, istate), 1, &
                 rmat%local_data(1, jstate), 1, ct, st)
@@ -559,23 +563,27 @@ SUBROUTINE rotate_zij_real(istate, jstate, st, ct, zij)
 
       INTEGER                                  :: iatom, natom, nstate
 #if defined(__SCALAPACK)
+      INTEGER                                  :: lwork, info
       INTEGER, DIMENSION(9)                    :: desc
+      REAL(dp), DIMENSION(:), ALLOCATABLE      :: work
 #else
       INTEGER                                  :: stride
 #endif
 
-      CALL cp_fm_get_info(zij(1, 1), nrow_global=nstate)
-
       natom = SIZE(zij, 1)
       DO iatom = 1, natom
+         CALL cp_fm_get_info(zij(iatom, 1), nrow_global=nstate)
 #if defined(__SCALAPACK)
+         lwork = 2*nstate + 1
+         ALLOCATE (work(lwork))
          desc(:) = zij(iatom, 1)%matrix_struct%descriptor(:)
-         CALL pzrot(nstate, zij(iatom, 1)%local_data(1, 1), &
+         CALL pdrot(nstate, zij(iatom, 1)%local_data(1, 1), &
                     1, istate, desc, 1, zij(iatom, 1)%local_data(1, 1), &
-                    1, jstate, desc, nstate, ct, st)
-         CALL pzrot(nstate, zij(iatom, 1)%local_data(1, 1), &
-                    istate, 1, desc, 1, zij(iatom, 1)%local_data(1, 1), &
-                    jstate, 1, desc, nstate, ct, st)
+                    1, jstate, desc, 1, ct, st, work, lwork, info)
+         CALL pdrot(nstate, zij(iatom, 1)%local_data(1, 1), &
+                    istate, 1, desc, nstate, zij(iatom, 1)%local_data(1, 1), &
+                    jstate, 1, desc, nstate, ct, st, work, lwork, info)
+         DEALLOCATE (work)
 #else
          CALL drot(nstate, zij(iatom, 1)%local_data(1, istate), &
                    1, zij(iatom, 1)%local_data(1, jstate), 1, ct, st)
@@ -607,16 +615,16 @@ SUBROUTINE check_tolerance_real(zij_fm_set, tolerance)
       REAL(dp), DIMENSION(:, :), POINTER                 :: diag
       TYPE(cp_fm_type)                                   :: force
 
+      CALL cp_fm_create(force, zij_fm_set(1, 1)%matrix_struct)
+      CALL cp_fm_set_all(force, 0._dp)
+
       NULLIFY (diag, col_indices, row_indices)
       natom = SIZE(zij_fm_set, 1)
       CALL cp_fm_get_info(zij_fm_set(1, 1), nrow_local=nrow_local, &
                           ncol_local=ncol_local, nrow_global=nrow_global, &
                           row_indices=row_indices, col_indices=col_indices)
       ALLOCATE (diag(nrow_global, natom))
 
-      CALL cp_fm_create(force, zij_fm_set(1, 1)%matrix_struct)
-      CALL cp_fm_set_all(force, 0._dp)
-
       DO iatom = 1, natom
          DO istate = 1, nrow_global
             CALL cp_fm_get_element(zij_fm_set(iatom, 1), istate, istate, diag(istate, iatom))
@@ -627,8 +635,8 @@ SUBROUTINE check_tolerance_real(zij_fm_set, tolerance)
          DO jstate = 1, ncol_local
             grad_ij = 0.0_dp
             DO iatom = 1, natom
-               zii = diag(istate, iatom)
-               zjj = diag(jstate, iatom)
+               zii = diag(row_indices(istate), iatom)
+               zjj = diag(col_indices(jstate), iatom)
                zij = zij_fm_set(iatom, 1)%local_data(istate, jstate)
                grad_ij = grad_ij + 4.0_dp*zij*(zjj - zii)
             END DO
@@ -639,7 +647,6 @@ SUBROUTINE check_tolerance_real(zij_fm_set, tolerance)
       DEALLOCATE (diag)
 
       CALL cp_fm_maxabsval(force, tolerance)
-
       CALL cp_fm_release(force)
 
    END SUBROUTINE check_tolerance_real