Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 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_add_on_diag, dbcsr_dot, dbcsr_get_info, dbcsr_iterator_blocks_left, &
20 : dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
21 : dbcsr_p_type, dbcsr_scale, dbcsr_set, dbcsr_type
22 : USE kinds, ONLY: dp
23 : USE message_passing, ONLY: mp_para_env_type
24 : USE qs_kind_types, ONLY: qs_kind_type
25 : #include "./base/base_uses.f90"
26 :
27 : IMPLICIT NONE
28 :
29 : PRIVATE
30 :
31 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_atomic_block'
32 :
33 : PUBLIC :: calculate_atomic_block_dm
34 :
35 : TYPE atom_matrix_type
36 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: mat => NULL()
37 : END TYPE atom_matrix_type
38 :
39 : CONTAINS
40 :
41 : ! **************************************************************************************************
42 : !> \brief returns a block diagonal density matrix. Blocks correspond to the atomic densities.
43 : !> \param pmatrix ...
44 : !> \param matrix_s ...
45 : !> \param atomic_kind_set ...
46 : !> \param qs_kind_set ...
47 : !> \param nspin ...
48 : !> \param nelectron_spin ...
49 : !> \param ounit ...
50 : !> \param para_env ...
51 : ! **************************************************************************************************
52 4575 : SUBROUTINE calculate_atomic_block_dm(pmatrix, matrix_s, atomic_kind_set, &
53 4575 : qs_kind_set, nspin, nelectron_spin, ounit, para_env)
54 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT) :: pmatrix
55 : TYPE(dbcsr_type), INTENT(INOUT) :: matrix_s
56 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
57 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
58 : INTEGER, INTENT(IN) :: nspin
59 : INTEGER, DIMENSION(:), INTENT(IN) :: nelectron_spin
60 : INTEGER, INTENT(IN) :: ounit
61 : TYPE(mp_para_env_type) :: para_env
62 :
63 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_atomic_block_dm'
64 :
65 : INTEGER :: blk, handle, icol, ikind, irow, ispin, &
66 : nc, nkind, nocc(2)
67 4575 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
68 4575 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: nok
69 4575 : REAL(dp), DIMENSION(:, :), POINTER :: pdata
70 : REAL(KIND=dp) :: rds, rscale, trps1
71 4575 : TYPE(atom_matrix_type), ALLOCATABLE, DIMENSION(:) :: pmat
72 : TYPE(atomic_kind_type), POINTER :: atomic_kind
73 : TYPE(dbcsr_iterator_type) :: iter
74 : TYPE(dbcsr_type), POINTER :: matrix_p
75 : TYPE(qs_kind_type), POINTER :: qs_kind
76 :
77 4575 : CALL timeset(routineN, handle)
78 :
79 4575 : nkind = SIZE(atomic_kind_set)
80 4575 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
81 21945 : ALLOCATE (pmat(nkind))
82 13725 : ALLOCATE (nok(2, nkind))
83 :
84 : ! precompute the atomic blocks corresponding to spherical atoms
85 12795 : DO ikind = 1, nkind
86 8220 : atomic_kind => atomic_kind_set(ikind)
87 8220 : qs_kind => qs_kind_set(ikind)
88 8220 : NULLIFY (pmat(ikind)%mat)
89 8220 : IF (ounit > 0) THEN
90 : WRITE (UNIT=ounit, FMT="(/,T2,A)") &
91 1968 : "Guess for atomic kind: "//TRIM(atomic_kind%name)
92 : END IF
93 : CALL calculate_atomic_orbitals(atomic_kind, qs_kind, iunit=ounit, &
94 8220 : pmat=pmat(ikind)%mat, nocc=nocc)
95 29235 : nok(1:2, ikind) = nocc(1:2)
96 : END DO
97 :
98 4575 : rscale = 1.0_dp
99 4575 : IF (nspin == 2) rscale = 0.5_dp
100 :
101 10109 : DO ispin = 1, nspin
102 5534 : IF ((ounit > 0) .AND. (nspin > 1)) THEN
103 502 : WRITE (UNIT=ounit, FMT="(/,T2,A,I0)") "Spin ", ispin
104 : END IF
105 :
106 5534 : matrix_p => pmatrix(ispin)%matrix
107 5534 : CALL dbcsr_set(matrix_p, 0.0_dp)
108 :
109 5534 : nocc(ispin) = 0
110 5534 : CALL dbcsr_iterator_start(iter, matrix_p)
111 82087 : DO WHILE (dbcsr_iterator_blocks_left(iter))
112 76553 : CALL dbcsr_iterator_next_block(iter, irow, icol, pdata, blk)
113 76553 : ikind = kind_of(irow)
114 82087 : IF (icol .EQ. irow) THEN
115 11648 : IF (ispin == 1) THEN
116 : pdata(:, :) = pmat(ikind)%mat(:, :, 1)*rscale + &
117 1008167 : pmat(ikind)%mat(:, :, 2)*rscale
118 : ELSE
119 : pdata(:, :) = pmat(ikind)%mat(:, :, 1)*rscale - &
120 175531 : pmat(ikind)%mat(:, :, 2)*rscale
121 : END IF
122 11648 : nocc(ispin) = nocc(ispin) + nok(ispin, ikind)
123 : END IF
124 : END DO
125 5534 : CALL dbcsr_iterator_stop(iter)
126 :
127 5534 : CALL dbcsr_dot(matrix_p, matrix_s, trps1)
128 5534 : rds = 0.0_dp
129 : ! could be a ghost-atoms-only simulation
130 5534 : IF (nelectron_spin(ispin) > 0) THEN
131 5406 : rds = REAL(nelectron_spin(ispin), dp)/trps1
132 : END IF
133 5534 : CALL dbcsr_scale(matrix_p, rds)
134 :
135 5534 : IF (ounit > 0) THEN
136 1363 : IF (nspin > 1) THEN
137 : WRITE (UNIT=ounit, FMT="(T2,A,I1)") &
138 502 : "Re-scaling the density matrix to get the right number of electrons for spin ", ispin
139 : ELSE
140 : WRITE (UNIT=ounit, FMT="(T2,A)") &
141 861 : "Re-scaling the density matrix to get the right number of electrons"
142 : END IF
143 1363 : WRITE (ounit, '(T19,A,T44,A,T67,A)') "# Electrons", "Trace(P)", "Scaling factor"
144 1363 : WRITE (ounit, '(T20,I10,T40,F12.3,T67,F14.3)') nelectron_spin(ispin), trps1, rds
145 : END IF
146 :
147 15643 : IF (nspin > 1) THEN
148 1918 : CALL para_env%sum(nocc)
149 1918 : IF (nelectron_spin(ispin) > nocc(ispin)) THEN
150 2 : rds = 0.99_dp
151 2 : CALL dbcsr_scale(matrix_p, rds)
152 2 : rds = (1.0_dp - rds)*nelectron_spin(ispin)
153 2 : CALL dbcsr_get_info(matrix_p, nfullcols_total=nc)
154 2 : rds = rds/REAL(nc, KIND=dp)
155 2 : CALL dbcsr_add_on_diag(matrix_p, rds)
156 2 : IF (ounit > 0) THEN
157 : WRITE (UNIT=ounit, FMT="(T4,A,/,T4,A,T59,F20.12)") &
158 1 : "More MOs than initial guess orbitals detected", &
159 2 : "Add constant to diagonal elements ", rds
160 : END IF
161 : END IF
162 : END IF
163 :
164 : END DO
165 :
166 12795 : DO ikind = 1, nkind
167 12795 : IF (ASSOCIATED(pmat(ikind)%mat)) THEN
168 8220 : DEALLOCATE (pmat(ikind)%mat)
169 : END IF
170 : END DO
171 4575 : DEALLOCATE (pmat)
172 :
173 4575 : DEALLOCATE (kind_of, nok)
174 :
175 4575 : CALL timestop(handle)
176 :
177 4575 : END SUBROUTINE calculate_atomic_block_dm
178 :
179 0 : END MODULE qs_atomic_block
|