Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Routine to return block diagonal density matrix. Blocks correspond to the atomic densities
10 : !> \par History
11 : !> 2006.03 Moved here from qs_scf.F [Joost VandeVondele]
12 : !> 2022.05 split from qs_initial_guess.F to break circular dependency [Harald Forbert]
13 : ! **************************************************************************************************
14 : MODULE qs_atomic_block
15 : USE atom_kind_orbitals, ONLY: calculate_atomic_orbitals
16 : USE atomic_kind_types, ONLY: atomic_kind_type,&
17 : get_atomic_kind_set
18 : USE cp_dbcsr_api, ONLY: &
19 : dbcsr_get_info, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
20 : dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_p_type, dbcsr_scale, &
21 : dbcsr_set, dbcsr_type
22 : USE cp_dbcsr_contrib, ONLY: dbcsr_add_on_diag,&
23 : dbcsr_dot
24 : USE kinds, ONLY: dp
25 : USE message_passing, ONLY: mp_para_env_type
26 : USE qs_kind_types, ONLY: qs_kind_type
27 : #include "./base/base_uses.f90"
28 :
29 : IMPLICIT NONE
30 :
31 : PRIVATE
32 :
33 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_atomic_block'
34 :
35 : PUBLIC :: calculate_atomic_block_dm
36 :
37 : TYPE atom_matrix_type
38 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: mat => NULL()
39 : END TYPE atom_matrix_type
40 :
41 : CONTAINS
42 :
43 : ! **************************************************************************************************
44 : !> \brief returns a block diagonal density matrix. Blocks correspond to the atomic densities.
45 : !> \param pmatrix ...
46 : !> \param matrix_s ...
47 : !> \param atomic_kind_set ...
48 : !> \param qs_kind_set ...
49 : !> \param nspin ...
50 : !> \param nelectron_spin ...
51 : !> \param ounit ...
52 : !> \param para_env ...
53 : ! **************************************************************************************************
54 4647 : SUBROUTINE calculate_atomic_block_dm(pmatrix, matrix_s, atomic_kind_set, &
55 4647 : qs_kind_set, nspin, nelectron_spin, ounit, para_env)
56 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT) :: pmatrix
57 : TYPE(dbcsr_type), INTENT(INOUT) :: matrix_s
58 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
59 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
60 : INTEGER, INTENT(IN) :: nspin
61 : INTEGER, DIMENSION(:), INTENT(IN) :: nelectron_spin
62 : INTEGER, INTENT(IN) :: ounit
63 : TYPE(mp_para_env_type) :: para_env
64 :
65 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_atomic_block_dm'
66 :
67 : INTEGER :: handle, icol, ikind, irow, ispin, nc, &
68 : nkind, nocc(2)
69 4647 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
70 4647 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: nok
71 4647 : REAL(dp), DIMENSION(:, :), POINTER :: pdata
72 : REAL(KIND=dp) :: rds, rscale, trps1
73 4647 : TYPE(atom_matrix_type), ALLOCATABLE, DIMENSION(:) :: pmat
74 : TYPE(atomic_kind_type), POINTER :: atomic_kind
75 : TYPE(dbcsr_iterator_type) :: iter
76 : TYPE(dbcsr_type), POINTER :: matrix_p
77 : TYPE(qs_kind_type), POINTER :: qs_kind
78 :
79 4647 : CALL timeset(routineN, handle)
80 :
81 4647 : nkind = SIZE(atomic_kind_set)
82 4647 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
83 22305 : ALLOCATE (pmat(nkind))
84 13941 : ALLOCATE (nok(2, nkind))
85 :
86 : ! precompute the atomic blocks corresponding to spherical atoms
87 13011 : DO ikind = 1, nkind
88 8364 : atomic_kind => atomic_kind_set(ikind)
89 8364 : qs_kind => qs_kind_set(ikind)
90 8364 : NULLIFY (pmat(ikind)%mat)
91 8364 : IF (ounit > 0) THEN
92 : WRITE (UNIT=ounit, FMT="(/,T2,A)") &
93 2031 : "Guess for atomic kind: "//TRIM(atomic_kind%name)
94 : END IF
95 : CALL calculate_atomic_orbitals(atomic_kind, qs_kind, iunit=ounit, &
96 8364 : pmat=pmat(ikind)%mat, nocc=nocc)
97 29739 : nok(1:2, ikind) = nocc(1:2)
98 : END DO
99 :
100 4647 : rscale = 1.0_dp
101 4647 : IF (nspin == 2) rscale = 0.5_dp
102 :
103 10259 : DO ispin = 1, nspin
104 5612 : IF ((ounit > 0) .AND. (nspin > 1)) THEN
105 512 : WRITE (UNIT=ounit, FMT="(/,T2,A,I0)") "Spin ", ispin
106 : END IF
107 :
108 5612 : matrix_p => pmatrix(ispin)%matrix
109 5612 : CALL dbcsr_set(matrix_p, 0.0_dp)
110 :
111 5612 : nocc(ispin) = 0
112 5612 : CALL dbcsr_iterator_start(iter, matrix_p)
113 83756 : DO WHILE (dbcsr_iterator_blocks_left(iter))
114 78144 : CALL dbcsr_iterator_next_block(iter, irow, icol, pdata)
115 78144 : ikind = kind_of(irow)
116 83756 : IF (icol .EQ. irow) THEN
117 11858 : IF (ispin == 1) THEN
118 : pdata(:, :) = pmat(ikind)%mat(:, :, 1)*rscale + &
119 1024765 : pmat(ikind)%mat(:, :, 2)*rscale
120 : ELSE
121 : pdata(:, :) = pmat(ikind)%mat(:, :, 1)*rscale - &
122 175319 : pmat(ikind)%mat(:, :, 2)*rscale
123 : END IF
124 11858 : nocc(ispin) = nocc(ispin) + nok(ispin, ikind)
125 : END IF
126 : END DO
127 5612 : CALL dbcsr_iterator_stop(iter)
128 :
129 5612 : CALL dbcsr_dot(matrix_p, matrix_s, trps1)
130 5612 : rds = 0.0_dp
131 : ! could be a ghost-atoms-only simulation
132 5612 : IF (nelectron_spin(ispin) > 0) THEN
133 5480 : rds = REAL(nelectron_spin(ispin), dp)/trps1
134 : END IF
135 5612 : CALL dbcsr_scale(matrix_p, rds)
136 :
137 5612 : IF (ounit > 0) THEN
138 1404 : IF (nspin > 1) THEN
139 : WRITE (UNIT=ounit, FMT="(T2,A,I1)") &
140 512 : "Re-scaling the density matrix to get the right number of electrons for spin ", ispin
141 : ELSE
142 : WRITE (UNIT=ounit, FMT="(T2,A)") &
143 892 : "Re-scaling the density matrix to get the right number of electrons"
144 : END IF
145 1404 : WRITE (ounit, '(T19,A,T44,A,T67,A)') "# Electrons", "Trace(P)", "Scaling factor"
146 1404 : WRITE (ounit, '(T20,I10,T40,F12.3,T67,F14.3)') nelectron_spin(ispin), trps1, rds
147 : END IF
148 :
149 21483 : IF (nspin > 1) THEN
150 1930 : CALL para_env%sum(nocc)
151 1930 : IF (nelectron_spin(ispin) > nocc(ispin)) THEN
152 2 : rds = 0.99_dp
153 2 : CALL dbcsr_scale(matrix_p, rds)
154 2 : rds = (1.0_dp - rds)*nelectron_spin(ispin)
155 2 : CALL dbcsr_get_info(matrix_p, nfullcols_total=nc)
156 2 : rds = rds/REAL(nc, KIND=dp)
157 2 : CALL dbcsr_add_on_diag(matrix_p, rds)
158 2 : IF (ounit > 0) THEN
159 : WRITE (UNIT=ounit, FMT="(T4,A,/,T4,A,T59,F20.12)") &
160 1 : "More MOs than initial guess orbitals detected", &
161 2 : "Add constant to diagonal elements ", rds
162 : END IF
163 : END IF
164 : END IF
165 :
166 : END DO
167 :
168 13011 : DO ikind = 1, nkind
169 13011 : IF (ASSOCIATED(pmat(ikind)%mat)) THEN
170 8364 : DEALLOCATE (pmat(ikind)%mat)
171 : END IF
172 : END DO
173 4647 : DEALLOCATE (pmat)
174 :
175 4647 : DEALLOCATE (kind_of, nok)
176 :
177 4647 : CALL timestop(handle)
178 :
179 4647 : END SUBROUTINE calculate_atomic_block_dm
180 :
181 0 : END MODULE qs_atomic_block
|
