Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Utility subroutine for qs energy calculation
10 : !> \par History
11 : !> none
12 : !> \author MK (29.10.2002)
13 : ! **************************************************************************************************
14 : MODULE qs_matrix_w
15 : USE cp_cfm_types, ONLY: cp_cfm_create,&
16 : cp_cfm_get_info,&
17 : cp_cfm_release,&
18 : cp_cfm_to_fm,&
19 : cp_cfm_type
20 : USE cp_control_types, ONLY: dft_control_type
21 : USE cp_dbcsr_api, ONLY: dbcsr_p_type,&
22 : dbcsr_set
23 : USE cp_fm_pool_types, ONLY: cp_fm_pool_p_type,&
24 : fm_pool_get_el_struct
25 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
26 : cp_fm_struct_release,&
27 : cp_fm_struct_type
28 : USE cp_fm_types, ONLY: cp_fm_create,&
29 : cp_fm_release,&
30 : cp_fm_type
31 : USE dft_plus_u, ONLY: lowdin_kp_apply_mo,&
32 : plus_u
33 : USE input_constants, ONLY: plus_u_lowdin
34 : USE kinds, ONLY: dp
35 : USE kpoint_methods, ONLY: kpoint_density_matrices,&
36 : kpoint_density_transform,&
37 : kpoint_ot_energy_weighted_density
38 : USE kpoint_types, ONLY: kpoint_env_type,&
39 : kpoint_type
40 : USE parallel_gemm_api, ONLY: parallel_gemm
41 : USE qs_density_matrices, ONLY: calculate_w_matrix,&
42 : calculate_w_matrix_ot
43 : USE qs_environment_types, ONLY: get_qs_env,&
44 : qs_environment_type
45 : USE qs_kpoint_operators, ONLY: kpoint_operator_context_create,&
46 : kpoint_operator_context_release,&
47 : kpoint_operator_context_type,&
48 : kpoint_operator_get,&
49 : kpoint_operator_prepare,&
50 : kpoint_slot_ks
51 : USE qs_matrix_pools, ONLY: mpools_get
52 : USE qs_mo_types, ONLY: get_mo_set,&
53 : mo_set_type
54 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
55 : USE qs_rho_types, ONLY: qs_rho_get,&
56 : qs_rho_type
57 : USE scf_control_types, ONLY: scf_control_type
58 : #include "./base/base_uses.f90"
59 :
60 : IMPLICIT NONE
61 :
62 : PRIVATE
63 :
64 : ! *** Global parameters ***
65 :
66 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_matrix_w'
67 :
68 : PUBLIC :: compute_matrix_w
69 :
70 : CONTAINS
71 :
72 : ! **************************************************************************************************
73 : !> \brief Refactoring of qs_energies_scf. Moves computation of matrix_w
74 : !> into separate subroutine
75 : !> \param qs_env ...
76 : !> \param calc_forces ...
77 : !> \param include_lowdin_forces include the Lowdin overlap derivative before the W(k) transform
78 : !> \par History
79 : !> 05.2013 created [Florian Schiffmann]
80 : ! **************************************************************************************************
81 :
82 27447 : SUBROUTINE compute_matrix_w(qs_env, calc_forces, include_lowdin_forces)
83 : TYPE(qs_environment_type), POINTER :: qs_env
84 : LOGICAL, INTENT(IN) :: calc_forces
85 : LOGICAL, INTENT(IN), OPTIONAL :: include_lowdin_forces
86 :
87 : CHARACTER(len=*), PARAMETER :: routineN = 'compute_matrix_w'
88 :
89 : INTEGER :: handle, is, ispin, nao, nspin
90 : LOGICAL :: do_kpoints, has_unit_metric
91 27447 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s, matrix_w, &
92 27447 : mo_derivs, rho_ao
93 : TYPE(dft_control_type), POINTER :: dft_control
94 27447 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
95 : TYPE(mo_set_type), POINTER :: mo_set
96 : TYPE(qs_rho_type), POINTER :: rho
97 : TYPE(scf_control_type), POINTER :: scf_control
98 :
99 27447 : CALL timeset(routineN, handle)
100 :
101 : ! if calculate forces, time to compute the w matrix
102 27447 : CALL get_qs_env(qs_env, has_unit_metric=has_unit_metric)
103 :
104 27447 : IF (calc_forces .AND. .NOT. has_unit_metric) THEN
105 6709 : CALL get_qs_env(qs_env, do_kpoints=do_kpoints)
106 :
107 6709 : IF (do_kpoints) THEN
108 474 : BLOCK
109 1422 : TYPE(cp_fm_type), DIMENSION(2) :: fmwork
110 : TYPE(cp_fm_struct_type), POINTER :: ao_ao_fmstruct
111 : TYPE(cp_fm_type), POINTER :: mo_coeff
112 474 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp, matrix_s_kp, &
113 474 : matrix_w_kp
114 : TYPE(kpoint_type), POINTER :: kpoints
115 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
116 474 : POINTER :: sab_nl
117 :
118 : CALL get_qs_env(qs_env, &
119 : matrix_w_kp=matrix_w_kp, &
120 : matrix_ks_kp=matrix_ks_kp, &
121 : matrix_s_kp=matrix_s_kp, &
122 : sab_orb=sab_nl, &
123 : mos=mos, &
124 : dft_control=dft_control, &
125 : kpoints=kpoints, &
126 474 : scf_control=scf_control)
127 :
128 474 : CALL get_mo_set(mos(1), mo_coeff=mo_coeff, nao=nao)
129 : CALL cp_fm_struct_create(fmstruct=ao_ao_fmstruct, nrow_global=nao, ncol_global=nao, &
130 474 : template_fmstruct=mo_coeff%matrix_struct)
131 :
132 1422 : DO is = 1, SIZE(fmwork)
133 1422 : CALL cp_fm_create(fmwork(is), matrix_struct=ao_ao_fmstruct)
134 : END DO
135 474 : CALL cp_fm_struct_release(ao_ao_fmstruct)
136 :
137 : ! energy weighted density matrices in k-space
138 474 : IF (scf_control%use_ot) THEN
139 36 : IF (ALLOCATED(kpoints%lowdin_v)) CALL plus_u(qs_env, matrix_h=matrix_ks_kp)
140 36 : CALL kpoint_ot_energy_weighted_matrices(kpoints, matrix_ks_kp, fmwork(1))
141 : ELSE
142 438 : CALL kpoint_density_matrices(kpoints, energy_weighted=.TRUE.)
143 : END IF
144 : ! Property callers need the ordinary energy-weighted density. Force callers
145 : ! also need -dE_U/dS; transform both together, once, on the overlap support.
146 474 : IF (PRESENT(include_lowdin_forces)) THEN
147 470 : IF (include_lowdin_forces .AND. dft_control%dft_plus_u .AND. &
148 : dft_control%plus_u_method_id == plus_u_lowdin) THEN
149 8 : CALL plus_u(qs_env, matrix_w=matrix_w_kp)
150 : END IF
151 : END IF
152 : ! energy weighted density matrices in real space
153 : CALL kpoint_density_transform(kpoints, matrix_w_kp, .TRUE., &
154 474 : matrix_s_kp(1, 1)%matrix, sab_nl, fmwork)
155 :
156 1896 : DO is = 1, SIZE(fmwork)
157 1422 : CALL cp_fm_release(fmwork(is))
158 : END DO
159 :
160 : END BLOCK
161 : ELSE
162 :
163 6235 : NULLIFY (dft_control, rho_ao)
164 : CALL get_qs_env(qs_env, &
165 : matrix_w=matrix_w, &
166 : matrix_ks=matrix_ks, &
167 : matrix_s=matrix_s, &
168 : mo_derivs=mo_derivs, &
169 : scf_control=scf_control, &
170 : mos=mos, &
171 : rho=rho, &
172 6235 : dft_control=dft_control)
173 :
174 6235 : CALL qs_rho_get(rho, rho_ao=rho_ao)
175 :
176 6235 : nspin = SIZE(mos)
177 13162 : DO ispin = 1, nspin
178 6927 : mo_set => mos(ispin)
179 13162 : IF (dft_control%roks) THEN
180 92 : IF (scf_control%use_ot) THEN
181 40 : IF (ispin > 1) THEN
182 : ! not very elegant, indeed ...
183 20 : CALL dbcsr_set(matrix_w(ispin)%matrix, 0.0_dp)
184 : ELSE
185 : CALL calculate_w_matrix_ot(mo_set, mo_derivs(ispin)%matrix, &
186 20 : matrix_w(ispin)%matrix, matrix_s(1)%matrix)
187 : END IF
188 : ELSE
189 : CALL calculate_w_matrix(mo_set=mo_set, &
190 : matrix_ks=matrix_ks(ispin)%matrix, &
191 : matrix_p=rho_ao(ispin)%matrix, &
192 52 : matrix_w=matrix_w(ispin)%matrix)
193 : END IF
194 : ELSE
195 6835 : IF (scf_control%use_ot) THEN
196 : CALL calculate_w_matrix_ot(mo_set, mo_derivs(ispin)%matrix, &
197 2739 : matrix_w(ispin)%matrix, matrix_s(1)%matrix)
198 : ELSE
199 4096 : CALL calculate_w_matrix(mo_set, matrix_w(ispin)%matrix)
200 : END IF
201 : END IF
202 : END DO
203 :
204 : END IF
205 :
206 : END IF
207 :
208 27447 : CALL timestop(handle)
209 :
210 27447 : END SUBROUTINE compute_matrix_w
211 :
212 : ! **************************************************************************************************
213 : !> \brief Build OT energy-weighted densities using Hamiltonians on the k-point groups.
214 : !> \param kpoint K-point environment
215 : !> \param matrix_ks current real-space Kohn-Sham matrices
216 : !> \param ao_ao_fm square full-matrix template on the global parallel environment
217 : ! **************************************************************************************************
218 36 : SUBROUTINE kpoint_ot_energy_weighted_matrices(kpoint, matrix_ks, ao_ao_fm)
219 :
220 : TYPE(kpoint_type), POINTER :: kpoint
221 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks
222 : TYPE(cp_fm_type), INTENT(IN) :: ao_ao_fm
223 :
224 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_ot_energy_weighted_matrices'
225 : COMPLEX(KIND=dp), PARAMETER :: zone = (1.0_dp, 0.0_dp), &
226 : zzero = (0.0_dp, 0.0_dp)
227 :
228 : INTEGER :: handle, ikpgr, ispin, kplocal, kpmax, &
229 : nao, nmo, nspin
230 36 : REAL(KIND=dp), DIMENSION(:), POINTER :: occupation
231 : TYPE(cp_cfm_type) :: h_k, hc, wmat
232 36 : TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_ao_fm_pools
233 : TYPE(cp_fm_struct_type), POINTER :: h_k_struct
234 : TYPE(kpoint_env_type), POINTER :: kp
235 36 : TYPE(kpoint_operator_context_type) :: op_ctx
236 :
237 36 : CALL timeset(routineN, handle)
238 36 : NULLIFY (ao_ao_fm_pools, h_k_struct, occupation)
239 36 : CPASSERT(.NOT. kpoint%use_real_wfn)
240 36 : CPASSERT(ASSOCIATED(matrix_ks))
241 :
242 36 : nspin = SIZE(matrix_ks, 1)
243 36 : kplocal = SIZE(kpoint%kp_env)
244 102 : kpmax = MAXVAL(kpoint%kp_dist(2, :) - kpoint%kp_dist(1, :) + 1)
245 36 : CALL kpoint_operator_context_create(op_ctx, kpoint, ao_ao_fm, nspin)
246 : ! Capture the current Hamiltonian after SCF, independently of OT's iteration caches.
247 36 : CALL kpoint_operator_prepare(op_ctx, matrix_ks, kpoint_slot_ks)
248 36 : CALL mpools_get(kpoint%mpools, ao_ao_fm_pools=ao_ao_fm_pools)
249 36 : h_k_struct => fm_pool_get_el_struct(ao_ao_fm_pools(1)%pool)
250 36 : CALL cp_cfm_create(h_k, h_k_struct)
251 :
252 94 : DO ikpgr = 1, kpmax
253 174 : DO ispin = 1, nspin
254 : ! Every group calls the collective service; the AO/AO pool layout keeps H(k)
255 : ! valid even on inactive rounds, when fallback assembly may communicate globally.
256 : CALL kpoint_operator_get(op_ctx, ikpgr, ispin, matrix_ks, &
257 : cmat=h_k, &
258 80 : matrix_row=ispin, slot=kpoint_slot_ks)
259 138 : IF (ikpgr <= kplocal) THEN
260 78 : kp => kpoint%kp_env(ikpgr)%kpoint_env
261 : CALL get_mo_set(kp%mos(ispin), nmo=nmo, &
262 78 : occupation_numbers=occupation)
263 78 : CALL cp_cfm_get_info(kp%mos(ispin)%cmo_coeff, nrow_global=nao)
264 78 : CPASSERT(nmo >= 1)
265 78 : CPASSERT(SIZE(occupation) >= nmo)
266 :
267 : ! Smearing occupations are k-point resolved, whereas mo_set%homo is not.
268 78 : CALL cp_cfm_create(hc, kp%mos(ispin)%cmo_coeff%matrix_struct, set_zero=.TRUE.)
269 : CALL parallel_gemm("N", "N", nao, nmo, nao, zone, h_k, &
270 78 : kp%mos(ispin)%cmo_coeff, zzero, hc)
271 78 : IF (ALLOCATED(kpoint%lowdin_v)) THEN
272 : CALL lowdin_kp_apply_mo(kp, kpoint%lowdin_v(:, ispin), &
273 12 : kp%mos(ispin)%cmo_coeff, hc)
274 : END IF
275 :
276 78 : CALL cp_cfm_create(wmat, kp%wmat(1, ispin)%matrix_struct)
277 : CALL kpoint_ot_energy_weighted_density( &
278 78 : kp%mos(ispin)%cmo_coeff, hc, occupation, wmat)
279 : ! wmat remains split-FM storage; keep that conversion at its boundary.
280 78 : CALL cp_cfm_to_fm(wmat, kp%wmat(1, ispin), kp%wmat(2, ispin))
281 :
282 78 : CALL cp_cfm_release(wmat)
283 78 : CALL cp_cfm_release(hc)
284 : END IF
285 : END DO
286 : END DO
287 :
288 36 : CALL cp_cfm_release(h_k)
289 36 : CALL kpoint_operator_context_release(op_ctx)
290 36 : CALL timestop(handle)
291 :
292 72 : END SUBROUTINE kpoint_ot_energy_weighted_matrices
293 :
294 : END MODULE qs_matrix_w
|
