kpoint (cp2k#3609)

hfp · Aug 6, 2024 · 729dda6 · 729dda6
1 parent 615b0eb
commit 729dda6
Show file tree

Hide file tree

Showing 4 changed files with 43 additions and 10 deletions.
diff --git a/src/aobasis/orbital_transformation_matrices.F b/src/aobasis/orbital_transformation_matrices.F
@@ -327,8 +327,6 @@ SUBROUTINE calculate_rotmat(orbrotmat, rotmat, lval)
       REAL(KIND=dp), DIMENSION(-2:2, -2:2)               :: r2
       REAL(KIND=dp), DIMENSION(3, 3)                     :: t
 
-      MARK_USED(rotmat)
-
       CALL release_rotmat(orbrotmat)
 
       ALLOCATE (orbrotmat(0:lval))

diff --git a/src/kpoint_methods.F b/src/kpoint_methods.F
@@ -231,6 +231,8 @@ SUBROUTINE kpoint_initialize(kpoint, particle_set, cell)
                   NULLIFY (kpoint%kind_rotmat(i, j)%rmat)
                END DO
             END DO
+            ALLOCATE (kpoint%ibrot(ns))
+            kpoint%ibrot(1:ns) = crys_sym%ibrot(1:ns)
          END IF
 
          CALL release_csym_type(crys_sym)
@@ -1152,8 +1154,8 @@ SUBROUTINE kpoint_density_transform(kpoint, denmat, wtype, tempmat, sab_nl, fmwo
 
       CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_density_transform'
 
-      INTEGER                                            :: handle, ic, ik, ikk, indx, ir, is, &
-                                                            ispin, nc, nimg, nkp, nspin
+      INTEGER                                            :: handle, ic, ik, ikk, indx, ir, ira, is, &
+                                                            ispin, jr, nc, nimg, nkp, nspin
       INTEGER, DIMENSION(:, :, :), POINTER               :: cell_to_index
       LOGICAL                                            :: aux_fit, do_ext, do_symmetric, my_kpgrp, &
                                                             real_only
@@ -1292,7 +1294,12 @@ SUBROUTINE kpoint_density_transform(kpoint, denmat, wtype, tempmat, sab_nl, fmwo
                wkpx = wkp(ik)/REAL(kpsym%nwght, KIND=dp)
                DO is = 1, kpsym%nwght
                   ir = ABS(kpsym%rotp(is))
-                  kind_rot => kpoint%kind_rotmat(ir, :)
+                  ira = 0
+                  DO jr = 1, SIZE(kpoint%ibrot)
+                     IF (ir == kpoint%ibrot(jr)) ira = jr
+                  END DO
+                  CPASSERT(ira > 0)
+                  kind_rot => kpoint%kind_rotmat(ira, :)
                   IF (real_only) THEN
                      CALL symtrans(srpmat, rpmat, kind_rot, kpsym%rot(1:3, 1:3, is), &
                                    kpsym%f0(:, is), kpoint%atype, symmetric=.TRUE.)
@@ -1582,5 +1589,29 @@ SUBROUTINE symtrans(smat, pmat, kmat, rot, f0, atype, symmetric, antisymmetric)
    END SUBROUTINE symtrans
 
 ! **************************************************************************************************
+!> \brief ...
+!> \param mat ...
+! **************************************************************************************************
+   SUBROUTINE matprint(mat)
+      TYPE(dbcsr_type), POINTER                          :: mat
+
+      INTEGER                                            :: i, icol, irow
+      REAL(KIND=dp), DIMENSION(:, :), POINTER            :: mblock
+      TYPE(dbcsr_iterator_type)                          :: iter
+
+      CALL dbcsr_iterator_start(iter, mat)
+      DO WHILE (dbcsr_iterator_blocks_left(iter))
+         CALL dbcsr_iterator_next_block(iter, irow, icol, mblock)
+         !
+         WRITE (6, '(A,2I4)') 'BLOCK  ', irow, icol
+         DO i = 1, SIZE(mblock, 1)
+            WRITE (6, '(8F12.6)') mblock(i, :)
+         END DO
+         !
+      END DO
+      CALL dbcsr_iterator_stop(iter)
+
+   END SUBROUTINE matprint
+! **************************************************************************************************
 
 END MODULE kpoint_methods
diff --git a/src/kpoint_types.F b/src/kpoint_types.F
@@ -184,6 +184,7 @@ MODULE kpoint_types
       TYPE(kpoint_sym_p_type), DIMENSION(:), &
          POINTER                              :: kp_sym => Null()
       INTEGER, DIMENSION(:), POINTER          :: atype => Null()
+      INTEGER, DIMENSION(:), POINTER          :: ibrot => Null()
       TYPE(kind_rotmat_type), DIMENSION(:, :), &
          POINTER                              :: kind_rotmat => Null()
       ! pools
@@ -299,6 +300,7 @@ SUBROUTINE kpoint_release(kpoint)
          END IF
 
          IF (ASSOCIATED(kpoint%atype)) DEALLOCATE (kpoint%atype)
+         IF (ASSOCIATED(kpoint%ibrot)) DEALLOCATE (kpoint%ibrot)
 
          IF (ASSOCIATED(kpoint%kind_rotmat)) THEN
             DO i = 1, SIZE(kpoint%kind_rotmat, 1)

diff --git a/src/qs_basis_rotation_methods.F b/src/qs_basis_rotation_methods.F
@@ -50,8 +50,8 @@ SUBROUTINE qs_basis_rotation(qs_env, kpoints)
       TYPE(qs_environment_type), POINTER                 :: qs_env
       TYPE(kpoint_type), POINTER                         :: kpoints
 
-      INTEGER                                            :: ik, ikind, ir, ira, irot, lval, nkind, &
-                                                            nrot
+      INTEGER                                            :: ik, ikind, ir, ira, irot, jr, lval, &
+                                                            nkind, nrot
       REAL(KIND=dp), DIMENSION(3, 3)                     :: rotmat
       TYPE(cell_type), POINTER                           :: cell
       TYPE(dft_control_type), POINTER                    :: dft_control
@@ -84,12 +84,14 @@ SUBROUTINE qs_basis_rotation(qs_env, kpoints)
             IF (kpsym%apply_symmetry) THEN
                DO irot = 1, SIZE(kpsym%rotp)
                   ir = kpsym%rotp(irot)
-                  ira = ABS(ir)
-                  IF (ira > 0 .AND. ira <= nrot) THEN
+                  ira = 0
+                  DO jr = 1, SIZE(kpoints%ibrot)
+                     IF (ir == kpoints%ibrot(jr)) ira = jr
+                  END DO
+                  IF (ira > 0) THEN
                      IF (.NOT. ASSOCIATED(kpoints%kind_rotmat(ira, 1)%rmat)) THEN
                         rotmat(1:3, 1:3) = MATMUL(cell%h_inv, &
                                                   MATMUL(kpsym%rot(:, :, irot), cell%hmat))
-                        NULLIFY (orbrot)
                         CALL calculate_rotmat(orbrot, rotmat, lval)
                         IF (dft_control%qs_control%method_id == do_method_dftb) THEN
                            CPABORT("ROTMAT")