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_control_types, ONLY: dft_control_type
16 : USE cp_dbcsr_api, ONLY: dbcsr_init_p,&
17 : dbcsr_p_type,&
18 : dbcsr_release_p,&
19 : dbcsr_set,&
20 : dbcsr_type
21 : USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply
22 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
23 : cp_fm_struct_release,&
24 : cp_fm_struct_type
25 : USE cp_fm_types, ONLY: cp_fm_create,&
26 : cp_fm_release,&
27 : cp_fm_type
28 : USE kinds, ONLY: dp
29 : USE kpoint_methods, ONLY: kpoint_density_matrices,&
30 : kpoint_density_transform,&
31 : kpoint_ot_energy_weighted_density
32 : USE kpoint_types, ONLY: kpoint_env_type,&
33 : kpoint_type
34 : USE qs_density_matrices, ONLY: calculate_w_matrix,&
35 : calculate_w_matrix_ot
36 : USE qs_environment_types, ONLY: get_qs_env,&
37 : qs_environment_type
38 : USE qs_kpoint_operators, ONLY: kpoint_operator_context_create,&
39 : kpoint_operator_context_release,&
40 : kpoint_operator_context_type,&
41 : kpoint_operator_get,&
42 : kpoint_operator_prepare,&
43 : kpoint_slot_ks
44 : USE qs_mo_types, ONLY: get_mo_set,&
45 : mo_set_type
46 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
47 : USE qs_rho_types, ONLY: qs_rho_get,&
48 : qs_rho_type
49 : USE scf_control_types, ONLY: scf_control_type
50 : #include "./base/base_uses.f90"
51 :
52 : IMPLICIT NONE
53 :
54 : PRIVATE
55 :
56 : ! *** Global parameters ***
57 :
58 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_matrix_w'
59 :
60 : PUBLIC :: compute_matrix_w
61 :
62 : CONTAINS
63 :
64 : ! **************************************************************************************************
65 : !> \brief Refactoring of qs_energies_scf. Moves computation of matrix_w
66 : !> into separate subroutine
67 : !> \param qs_env ...
68 : !> \param calc_forces ...
69 : !> \par History
70 : !> 05.2013 created [Florian Schiffmann]
71 : ! **************************************************************************************************
72 :
73 26979 : SUBROUTINE compute_matrix_w(qs_env, calc_forces)
74 : TYPE(qs_environment_type), POINTER :: qs_env
75 : LOGICAL, INTENT(IN) :: calc_forces
76 :
77 : CHARACTER(len=*), PARAMETER :: routineN = 'compute_matrix_w'
78 :
79 : INTEGER :: handle, is, ispin, nao, nspin
80 : LOGICAL :: do_kpoints, has_unit_metric
81 26979 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s, matrix_w, &
82 26979 : mo_derivs, rho_ao
83 : TYPE(dft_control_type), POINTER :: dft_control
84 26979 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
85 : TYPE(mo_set_type), POINTER :: mo_set
86 : TYPE(qs_rho_type), POINTER :: rho
87 : TYPE(scf_control_type), POINTER :: scf_control
88 :
89 26979 : CALL timeset(routineN, handle)
90 :
91 : ! if calculate forces, time to compute the w matrix
92 26979 : CALL get_qs_env(qs_env, has_unit_metric=has_unit_metric)
93 :
94 26979 : IF (calc_forces .AND. .NOT. has_unit_metric) THEN
95 6747 : CALL get_qs_env(qs_env, do_kpoints=do_kpoints)
96 :
97 6747 : IF (do_kpoints) THEN
98 458 : BLOCK
99 1374 : TYPE(cp_fm_type), DIMENSION(2) :: fmwork
100 : TYPE(cp_fm_struct_type), POINTER :: ao_ao_fmstruct
101 : TYPE(cp_fm_type), POINTER :: mo_coeff
102 458 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp, matrix_s_kp, &
103 458 : matrix_w_kp
104 : TYPE(kpoint_type), POINTER :: kpoints
105 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
106 458 : POINTER :: sab_nl
107 :
108 : CALL get_qs_env(qs_env, &
109 : matrix_w_kp=matrix_w_kp, &
110 : matrix_ks_kp=matrix_ks_kp, &
111 : matrix_s_kp=matrix_s_kp, &
112 : sab_orb=sab_nl, &
113 : mos=mos, &
114 : kpoints=kpoints, &
115 458 : scf_control=scf_control)
116 :
117 458 : CALL get_mo_set(mos(1), mo_coeff=mo_coeff, nao=nao)
118 : CALL cp_fm_struct_create(fmstruct=ao_ao_fmstruct, nrow_global=nao, ncol_global=nao, &
119 458 : template_fmstruct=mo_coeff%matrix_struct)
120 :
121 1374 : DO is = 1, SIZE(fmwork)
122 1374 : CALL cp_fm_create(fmwork(is), matrix_struct=ao_ao_fmstruct)
123 : END DO
124 458 : CALL cp_fm_struct_release(ao_ao_fmstruct)
125 :
126 : ! energy weighted density matrices in k-space
127 458 : IF (scf_control%use_ot) THEN
128 24 : CALL kpoint_ot_energy_weighted_matrices(kpoints, matrix_ks_kp, fmwork(1))
129 : ELSE
130 434 : CALL kpoint_density_matrices(kpoints, energy_weighted=.TRUE.)
131 : END IF
132 : ! energy weighted density matrices in real space
133 : CALL kpoint_density_transform(kpoints, matrix_w_kp, .TRUE., &
134 458 : matrix_s_kp(1, 1)%matrix, sab_nl, fmwork)
135 :
136 1832 : DO is = 1, SIZE(fmwork)
137 1374 : CALL cp_fm_release(fmwork(is))
138 : END DO
139 :
140 : END BLOCK
141 : ELSE
142 :
143 6289 : NULLIFY (dft_control, rho_ao)
144 : CALL get_qs_env(qs_env, &
145 : matrix_w=matrix_w, &
146 : matrix_ks=matrix_ks, &
147 : matrix_s=matrix_s, &
148 : mo_derivs=mo_derivs, &
149 : scf_control=scf_control, &
150 : mos=mos, &
151 : rho=rho, &
152 6289 : dft_control=dft_control)
153 :
154 6289 : CALL qs_rho_get(rho, rho_ao=rho_ao)
155 :
156 6289 : nspin = SIZE(mos)
157 13326 : DO ispin = 1, nspin
158 7037 : mo_set => mos(ispin)
159 13326 : IF (dft_control%roks) THEN
160 168 : IF (scf_control%use_ot) THEN
161 116 : IF (ispin > 1) THEN
162 : ! not very elegant, indeed ...
163 58 : CALL dbcsr_set(matrix_w(ispin)%matrix, 0.0_dp)
164 : ELSE
165 : CALL calculate_w_matrix_ot(mo_set, mo_derivs(ispin)%matrix, &
166 58 : matrix_w(ispin)%matrix, matrix_s(1)%matrix)
167 : END IF
168 : ELSE
169 : CALL calculate_w_matrix(mo_set=mo_set, &
170 : matrix_ks=matrix_ks(ispin)%matrix, &
171 : matrix_p=rho_ao(ispin)%matrix, &
172 52 : matrix_w=matrix_w(ispin)%matrix)
173 : END IF
174 : ELSE
175 6869 : IF (scf_control%use_ot) THEN
176 : CALL calculate_w_matrix_ot(mo_set, mo_derivs(ispin)%matrix, &
177 2773 : matrix_w(ispin)%matrix, matrix_s(1)%matrix)
178 : ELSE
179 4096 : CALL calculate_w_matrix(mo_set, matrix_w(ispin)%matrix)
180 : END IF
181 : END IF
182 : END DO
183 :
184 : END IF
185 :
186 : END IF
187 :
188 26979 : CALL timestop(handle)
189 :
190 26979 : END SUBROUTINE compute_matrix_w
191 :
192 : ! **************************************************************************************************
193 : !> \brief Build OT energy-weighted densities using Hamiltonians on the k-point groups.
194 : !> \param kpoint K-point environment
195 : !> \param matrix_ks current real-space Kohn-Sham matrices
196 : !> \param ao_ao_fm square full-matrix template on the global parallel environment
197 : ! **************************************************************************************************
198 24 : SUBROUTINE kpoint_ot_energy_weighted_matrices(kpoint, matrix_ks, ao_ao_fm)
199 :
200 : TYPE(kpoint_type), POINTER :: kpoint
201 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks
202 : TYPE(cp_fm_type), INTENT(IN) :: ao_ao_fm
203 :
204 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_ot_energy_weighted_matrices'
205 :
206 : INTEGER :: handle, ikpgr, ispin, kplocal, kpmax, &
207 : nmo, nspin
208 24 : REAL(KIND=dp), DIMENSION(:), POINTER :: occupation
209 : TYPE(cp_fm_type) :: hc_im, hc_re
210 : TYPE(cp_fm_type), POINTER :: coeff_im, coeff_re
211 : TYPE(dbcsr_type), POINTER :: matrix_k_im, matrix_k_re
212 : TYPE(kpoint_env_type), POINTER :: kp
213 24 : TYPE(kpoint_operator_context_type) :: op_ctx
214 :
215 24 : CALL timeset(routineN, handle)
216 24 : NULLIFY (coeff_im, coeff_re, matrix_k_im, matrix_k_re, occupation)
217 24 : CPASSERT(.NOT. kpoint%use_real_wfn)
218 24 : CPASSERT(ASSOCIATED(matrix_ks))
219 :
220 24 : nspin = SIZE(matrix_ks, 1)
221 24 : kplocal = SIZE(kpoint%kp_env)
222 58 : kpmax = MAXVAL(kpoint%kp_dist(2, :) - kpoint%kp_dist(1, :) + 1)
223 24 : CALL kpoint_operator_context_create(op_ctx, kpoint, ao_ao_fm, nspin)
224 : ! Capture the current Hamiltonian after SCF, independently of OT's iteration caches.
225 24 : CALL kpoint_operator_prepare(op_ctx, matrix_ks, kpoint_slot_ks)
226 :
227 50 : DO ikpgr = 1, kpmax
228 86 : DO ispin = 1, nspin
229 36 : CALL dbcsr_init_p(matrix_k_re)
230 36 : CALL dbcsr_init_p(matrix_k_im)
231 : ! All groups participate in the operator calls, including inactive groups in the
232 : ! last round. The fallback may communicate over the global parallel environment.
233 : CALL kpoint_operator_get(op_ctx, ikpgr, ispin, matrix_ks, &
234 : mat_re=matrix_k_re, mat_im=matrix_k_im, &
235 36 : matrix_row=ispin, slot=kpoint_slot_ks)
236 36 : IF (ikpgr <= kplocal) THEN
237 34 : kp => kpoint%kp_env(ikpgr)%kpoint_env
238 : CALL get_mo_set(kp%mos(1, ispin), mo_coeff=coeff_re, nmo=nmo, &
239 34 : occupation_numbers=occupation)
240 34 : CALL get_mo_set(kp%mos(2, ispin), mo_coeff=coeff_im)
241 34 : CPASSERT(nmo >= 1)
242 34 : CPASSERT(SIZE(occupation) >= nmo)
243 :
244 : ! Smearing occupations are k-point resolved, whereas mo_set%homo is not.
245 34 : CALL cp_fm_create(hc_re, coeff_re%matrix_struct, set_zero=.TRUE.)
246 34 : CALL cp_fm_create(hc_im, coeff_im%matrix_struct, set_zero=.TRUE.)
247 34 : CALL cp_dbcsr_sm_fm_multiply(matrix_k_re, coeff_re, hc_re, nmo)
248 : CALL cp_dbcsr_sm_fm_multiply(matrix_k_im, coeff_im, hc_re, nmo, &
249 34 : alpha=-1.0_dp, beta=1.0_dp)
250 34 : CALL cp_dbcsr_sm_fm_multiply(matrix_k_re, coeff_im, hc_im, nmo)
251 : CALL cp_dbcsr_sm_fm_multiply(matrix_k_im, coeff_re, hc_im, nmo, &
252 34 : alpha=1.0_dp, beta=1.0_dp)
253 :
254 : CALL kpoint_ot_energy_weighted_density( &
255 : coeff_re, coeff_im, hc_re, hc_im, occupation, &
256 34 : kp%wmat(1, ispin), kp%wmat(2, ispin))
257 :
258 34 : CALL cp_fm_release(hc_im)
259 34 : CALL cp_fm_release(hc_re)
260 : END IF
261 36 : CALL dbcsr_release_p(matrix_k_im)
262 62 : CALL dbcsr_release_p(matrix_k_re)
263 : END DO
264 : END DO
265 :
266 24 : CALL kpoint_operator_context_release(op_ctx)
267 24 : CALL timestop(handle)
268 :
269 48 : END SUBROUTINE kpoint_ot_energy_weighted_matrices
270 :
271 : END MODULE qs_matrix_w
|
