More parameters, Mulliken charges for xTB

src/mulliken.F

@@ -29,14 +29,18 @@ MODULE mulliken
 
 ! *** Public subroutines ***
 
-   PUBLIC :: mulliken_charges, mulliken_restraint, compute_bond_order
+   PUBLIC :: mulliken_charges, ao_charges, mulliken_restraint, compute_bond_order
 
    INTERFACE mulliken_charges
       MODULE PROCEDURE mulliken_charges_a, mulliken_charges_b, mulliken_charges_c, &
          mulliken_charges_s, &
          mulliken_charges_akp, mulliken_charges_bkp, mulliken_charges_ckp
    END INTERFACE
 
+   INTERFACE ao_charges
+      MODULE PROCEDURE ao_charges_1, ao_charges_kp
+   END INTERFACE
+
 CONTAINS
 
 ! **************************************************************************************************
@@ -681,6 +685,125 @@ SUBROUTINE compute_bond_order(psmat, spmat, bond_order)
 
    END SUBROUTINE compute_bond_order
 
+! **************************************************************************************************
+!> \brief compute the mulliken charges for a single atom (AO resolved)
+!> \param p_matrix , s_matrix, para_env ...
+!> \param s_matrix ...
+!> \param charges ...
+!> \param iatom ...
+!> \param para_env ...
+!> \par History
+!>      06.2004 created [Joost VandeVondele]
+!>      10.2018 adapted [JGH]
+!> \note
+! **************************************************************************************************
+   SUBROUTINE ao_charges_1(p_matrix, s_matrix, charges, iatom, para_env)
+      TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: p_matrix
+      TYPE(dbcsr_type), POINTER                          :: s_matrix
+      REAL(KIND=dp), DIMENSION(:), INTENT(INOUT)         :: charges
+      INTEGER, INTENT(IN)                                :: iatom
+      TYPE(cp_para_env_type), POINTER                    :: para_env
+
+      CHARACTER(len=*), PARAMETER :: routineN = 'ao_charges_1', routineP = moduleN//':'//routineN
+
+      INTEGER                                            :: blk, handle, i, iblock_col, iblock_row, &
+                                                            ispin, j, nspin
+      LOGICAL                                            :: found
+      REAL(kind=dp)                                      :: mult
+      REAL(KIND=dp), DIMENSION(:, :), POINTER            :: p_block, s_block
+      TYPE(dbcsr_iterator_type)                          :: iter
+
+      CALL timeset(routineN, handle)
+
+      nspin = SIZE(p_matrix)
+      charges = 0.0_dp
+      DO ispin = 1, nspin
+         CALL dbcsr_iterator_start(iter, s_matrix)
+         DO WHILE (dbcsr_iterator_blocks_left(iter))
+            NULLIFY (s_block, p_block)
+            CALL dbcsr_iterator_next_block(iter, iblock_row, iblock_col, s_block, blk)
+            CALL dbcsr_get_block_p(matrix=p_matrix(ispin)%matrix, &
+                                   row=iblock_row, col=iblock_col, BLOCK=p_block, found=found)
+
+            ! we can cycle if a block is not present
+            IF (.NOT. found) CYCLE
+            IF (.NOT. (ASSOCIATED(s_block) .AND. ASSOCIATED(p_block))) CYCLE
+
+            IF (iblock_row == iatom) THEN
+               IF (iblock_row == iblock_col) THEN
+                  mult = 0.5_dp ! avoid double counting of diagonal blocks
+               ELSE
+                  mult = 1.0_dp
+               ENDIF
+               DO j = 1, SIZE(p_block, 2)
+                  DO i = 1, SIZE(p_block, 1)
+                     charges(i) = charges(i)+mult*p_block(i, j)*s_block(i, j)
+                  END DO
+               END DO
+            ELSEIF (iblock_col == iatom) THEN
+               DO i = 1, SIZE(p_block, 1)
+                  DO j = 1, SIZE(p_block, 2)
+                     charges(j) = charges(j)+mult*p_block(i, j)*s_block(i, j)
+                  END DO
+               END DO
+            END IF
+         ENDDO
+         CALL dbcsr_iterator_stop(iter)
+      ENDDO
+      CALL mp_sum(charges, para_env%group)
+
+      CALL timestop(handle)
+
+   END SUBROUTINE ao_charges_1
+
+! **************************************************************************************************
+!> \brief ...
+!> \param p_matrix_kp ...
+!> \param s_matrix_kp ...
+!> \param charges ...
+!> \param iatom ...
+!> \param para_env ...
+! **************************************************************************************************
+   SUBROUTINE ao_charges_kp(p_matrix_kp, s_matrix_kp, charges, iatom, para_env)
+
+      TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: p_matrix_kp, s_matrix_kp
+      REAL(KIND=dp), DIMENSION(:), INTENT(INOUT)         :: charges
+      INTEGER, INTENT(IN)                                :: iatom
+      TYPE(cp_para_env_type), POINTER                    :: para_env
+
+      CHARACTER(len=*), PARAMETER :: routineN = 'ao_charges_kp', routineP = moduleN//':'//routineN
+
+      INTEGER                                            :: handle, ic, na
+      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: charge_im
+      TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: p_matrix
+      TYPE(dbcsr_type), POINTER                          :: s_matrix
+
+      CALL timeset(routineN, handle)
+
+      IF (ASSOCIATED(p_matrix_kp) .AND. ASSOCIATED(s_matrix_kp)) THEN
+
+         charges = 0.0_dp
+         na = SIZE(charges, 1)
+         ALLOCATE (charge_im(na))
+
+         DO ic = 1, SIZE(s_matrix_kp, 2)
+            NULLIFY (p_matrix, s_matrix)
+            p_matrix => p_matrix_kp(:, ic)
+            s_matrix => s_matrix_kp(1, ic)%matrix
+            IF (ASSOCIATED(p_matrix) .AND. ASSOCIATED(s_matrix)) THEN
+               CALL ao_charges_1(p_matrix, s_matrix, charge_im, iatom, para_env)
+               charges(:) = charges(:)+charge_im(:)
+            END IF
+         END DO
+
+         DEALLOCATE (charge_im)
+
+      END IF
+
+      CALL timestop(handle)
+
+   END SUBROUTINE ao_charges_kp
+
 ! **************************************************************************************************
 
 END MODULE mulliken