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 Methods for preparing and committing k-point orbital states to the QS environment.
10 : ! **************************************************************************************************
11 : MODULE qs_kpoint_state
12 : USE cp_cfm_basic_linalg, ONLY: cp_cfm_gemm,&
13 : cp_cfm_scale
14 : USE cp_cfm_diag, ONLY: cp_cfm_geeig,&
15 : cp_cfm_heevd
16 : USE cp_cfm_types, ONLY: cp_cfm_create,&
17 : cp_cfm_release,&
18 : cp_cfm_to_cfm,&
19 : cp_cfm_to_fm,&
20 : cp_cfm_type,&
21 : cp_fm_to_cfm
22 : USE cp_dbcsr_api, ONLY: dbcsr_p_type
23 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
24 : cp_fm_struct_release,&
25 : cp_fm_struct_type
26 : USE cp_fm_types, ONLY: cp_fm_create,&
27 : cp_fm_get_info,&
28 : cp_fm_release,&
29 : cp_fm_set_all,&
30 : cp_fm_to_fm,&
31 : cp_fm_type,&
32 : cp_fm_vectorsnorm
33 : USE kinds, ONLY: dp
34 : USE kpoint_methods, ONLY: kpoint_density_matrices,&
35 : kpoint_density_transform,&
36 : kpoint_set_mo_occupation
37 : USE kpoint_types, ONLY: get_kpoint_info,&
38 : kpoint_env_type,&
39 : kpoint_type
40 : USE mathconstants, ONLY: z_one,&
41 : z_zero
42 : USE message_passing, ONLY: mp_para_env_type
43 : USE qs_energy_types, ONLY: qs_energy_type
44 : USE qs_environment_types, ONLY: get_qs_env,&
45 : qs_environment_type
46 : USE qs_ks_types, ONLY: qs_ks_did_change
47 : USE qs_mo_types, ONLY: get_mo_set,&
48 : mo_set_type
49 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
50 : USE qs_rho_methods, ONLY: qs_rho_update_rho
51 : USE qs_rho_types, ONLY: qs_rho_get,&
52 : qs_rho_type
53 : USE qs_scf_types, ONLY: qs_scf_env_type
54 : USE scf_control_types, ONLY: scf_control_type
55 : #include "./base/base_uses.f90"
56 :
57 : IMPLICIT NONE
58 : PRIVATE
59 :
60 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_kpoint_state'
61 :
62 : PUBLIC :: qs_kpoint_copy_spin_mos, qs_kpoint_mos_initialized, &
63 : qs_kpoint_state_canonicalize_fixed, qs_kpoint_state_commit, &
64 : qs_kpoint_state_prepare_fixed_density
65 :
66 : CONTAINS
67 :
68 : ! **************************************************************************************************
69 : !> \brief Reconstructs fixed-rank complex k-point orbitals from a density kernel.
70 : !> The density and overlap must first be transferred to kp%ot_hmat and kp%ot_smat. The
71 : !> largest natural-orbital subspace gives a metric-defined fixed-rank projector without
72 : !> replacing the guess by the ground-state projector of its initial Hamiltonian.
73 : !> \param kpoints K-point environment containing the transferred density and overlap matrices.
74 : ! **************************************************************************************************
75 44 : SUBROUTINE qs_kpoint_state_prepare_fixed_density(kpoints)
76 : TYPE(kpoint_type), POINTER :: kpoints
77 :
78 : INTEGER :: ikpoint_local, ispin, nao, nmo, nocc
79 : LOGICAL :: use_real_wfn
80 44 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: density_eigenvalues
81 : TYPE(cp_cfm_type) :: density_metric, density_work, matrix_p, &
82 : matrix_s, matrix_s_work, projected
83 : TYPE(cp_fm_struct_type), POINTER :: active_struct, square_struct
84 : TYPE(cp_fm_type) :: active_im, active_re
85 : TYPE(cp_fm_type), POINTER :: mo_coeff_im, mo_coeff_re
86 : TYPE(kpoint_env_type), POINTER :: kp
87 :
88 44 : NULLIFY (active_struct, kp, mo_coeff_im, mo_coeff_re, square_struct)
89 :
90 0 : CPASSERT(ASSOCIATED(kpoints))
91 44 : CPASSERT(ASSOCIATED(kpoints%kp_env))
92 44 : CALL get_kpoint_info(kpoints, use_real_wfn=use_real_wfn)
93 44 : IF (use_real_wfn) THEN
94 0 : CPABORT("Fixed-density k-point OT initialization requires complex wavefunctions.")
95 : END IF
96 :
97 152 : DO ikpoint_local = 1, SIZE(kpoints%kp_env)
98 108 : kp => kpoints%kp_env(ikpoint_local)%kpoint_env
99 108 : CPASSERT(ASSOCIATED(kp))
100 108 : CPASSERT(ASSOCIATED(kp%mos))
101 108 : CPASSERT(ASSOCIATED(kp%ot_hmat))
102 108 : CPASSERT(ASSOCIATED(kp%ot_smat))
103 108 : CPASSERT(SIZE(kp%ot_hmat, 1) >= 2)
104 108 : CPASSERT(SIZE(kp%ot_hmat, 2) >= SIZE(kp%mos, 2))
105 108 : CPASSERT(SIZE(kp%ot_smat) >= 2)
106 :
107 108 : CALL cp_fm_get_info(kp%ot_hmat(1, 1), matrix_struct=square_struct, nrow_global=nao)
108 108 : CALL cp_cfm_create(matrix_p, square_struct, name="natural-orbital density P(k)")
109 108 : CALL cp_cfm_create(matrix_s, square_struct, name="natural-orbital overlap S(k)")
110 108 : CALL cp_fm_to_cfm(kp%ot_smat(1), kp%ot_smat(2), matrix_s)
111 :
112 236 : DO ispin = 1, SIZE(kp%mos, 2)
113 128 : CALL get_mo_set(kp%mos(1, ispin), homo=nocc, mo_coeff=mo_coeff_re, nmo=nmo)
114 128 : CALL get_mo_set(kp%mos(2, ispin), mo_coeff=mo_coeff_im)
115 128 : CPASSERT(ASSOCIATED(mo_coeff_re))
116 128 : CPASSERT(ASSOCIATED(mo_coeff_im))
117 128 : CPASSERT(nocc >= 1 .AND. nocc <= nmo)
118 :
119 : CALL cp_fm_struct_create(active_struct, template_fmstruct=mo_coeff_re%matrix_struct, &
120 128 : ncol_global=nocc)
121 128 : CALL cp_fm_to_cfm(kp%ot_hmat(1, ispin), kp%ot_hmat(2, ispin), matrix_p)
122 128 : CALL cp_cfm_create(projected, active_struct, name="fixed-density natural orbitals")
123 128 : CALL cp_cfm_create(density_metric, square_struct, name="natural-orbital metric")
124 128 : CALL cp_cfm_create(density_work, square_struct, name="natural-orbital work")
125 128 : CALL cp_cfm_create(matrix_s_work, square_struct, name="natural-orbital overlap work")
126 :
127 : ! The largest eigenvectors of S*P*S*C = S*C*n are the natural orbitals of
128 : ! S**(1/2)*P*S**(1/2). Negating the left-hand side lets the generalized
129 : ! eigensolver request only its lowest nocc eigenpairs when supported.
130 : CALL cp_cfm_gemm('N', 'N', nao, nao, nao, z_one, matrix_p, matrix_s, &
131 128 : z_zero, density_work)
132 : CALL cp_cfm_gemm('N', 'N', nao, nao, nao, z_one, matrix_s, density_work, &
133 128 : z_zero, density_metric)
134 128 : CALL cp_cfm_scale(-z_one, density_metric)
135 128 : CALL cp_cfm_to_cfm(matrix_s, matrix_s_work)
136 384 : ALLOCATE (density_eigenvalues(nocc))
137 : CALL cp_cfm_geeig(density_metric, matrix_s_work, projected, &
138 128 : density_eigenvalues, density_work, lowest_subset=.TRUE.)
139 128 : DEALLOCATE (density_eigenvalues)
140 :
141 128 : CALL cp_fm_create(active_re, active_struct, name="fixed-density orbitals real")
142 128 : CALL cp_fm_create(active_im, active_struct, name="fixed-density orbitals imaginary")
143 128 : CALL cp_cfm_to_fm(projected, active_re, active_im)
144 128 : CALL cp_fm_set_all(mo_coeff_re, 0.0_dp)
145 128 : CALL cp_fm_set_all(mo_coeff_im, 0.0_dp)
146 128 : CALL cp_fm_to_fm(active_re, mo_coeff_re, nocc)
147 128 : CALL cp_fm_to_fm(active_im, mo_coeff_im, nocc)
148 :
149 128 : CALL cp_fm_release(active_im)
150 128 : CALL cp_fm_release(active_re)
151 128 : CALL cp_cfm_release(matrix_s_work)
152 128 : CALL cp_cfm_release(density_work)
153 128 : CALL cp_cfm_release(density_metric)
154 128 : CALL cp_cfm_release(projected)
155 620 : CALL cp_fm_struct_release(active_struct)
156 : END DO
157 :
158 108 : CALL cp_cfm_release(matrix_s)
159 260 : CALL cp_cfm_release(matrix_p)
160 : END DO
161 :
162 88 : END SUBROUTINE qs_kpoint_state_prepare_fixed_density
163 :
164 : ! **************************************************************************************************
165 : !> \brief Rotates a uniformly occupied k-point subspace into Ritz states of the current
166 : !> Hamiltonian and stores their energy labels. The Hamiltonian must first be transferred
167 : !> to kp%ot_hmat. Since the rotation is internal to a uniformly occupied subspace, it does
168 : !> not change the density represented by the orbitals.
169 : !> \param kpoints K-point environment containing the orbitals and transferred Hamiltonians.
170 : ! **************************************************************************************************
171 44 : SUBROUTINE qs_kpoint_state_canonicalize_fixed(kpoints)
172 : TYPE(kpoint_type), POINTER :: kpoints
173 :
174 : INTEGER :: ikpoint_local, ispin, nao, nmo, nocc
175 : LOGICAL :: use_real_wfn
176 44 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ritz_values
177 44 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues_im, eigenvalues_re
178 : TYPE(cp_cfm_type) :: coeff, eigenvectors, hblock, hc, &
179 : matrix_h, rotated
180 : TYPE(cp_fm_struct_type), POINTER :: active_struct, square_struct
181 : TYPE(cp_fm_type) :: active_im, active_re
182 : TYPE(cp_fm_type), POINTER :: mo_coeff_im, mo_coeff_re
183 : TYPE(kpoint_env_type), POINTER :: kp
184 :
185 44 : NULLIFY (active_struct, eigenvalues_im, eigenvalues_re, kp, mo_coeff_im, mo_coeff_re, &
186 44 : square_struct)
187 :
188 0 : CPASSERT(ASSOCIATED(kpoints))
189 44 : CPASSERT(ASSOCIATED(kpoints%kp_env))
190 44 : CALL get_kpoint_info(kpoints, use_real_wfn=use_real_wfn)
191 44 : IF (use_real_wfn) THEN
192 0 : CPABORT("Fixed-density k-point OT canonicalization requires complex wavefunctions.")
193 : END IF
194 :
195 152 : DO ikpoint_local = 1, SIZE(kpoints%kp_env)
196 108 : kp => kpoints%kp_env(ikpoint_local)%kpoint_env
197 108 : CPASSERT(ASSOCIATED(kp))
198 108 : CPASSERT(ASSOCIATED(kp%mos))
199 108 : CPASSERT(ASSOCIATED(kp%ot_hmat))
200 108 : CPASSERT(SIZE(kp%ot_hmat, 1) >= 2)
201 108 : CPASSERT(SIZE(kp%ot_hmat, 2) >= SIZE(kp%mos, 2))
202 :
203 108 : CALL cp_fm_get_info(kp%ot_hmat(1, 1), matrix_struct=square_struct, nrow_global=nao)
204 108 : CALL cp_cfm_create(matrix_h, square_struct, name="fixed-density Ritz Hamiltonian")
205 :
206 236 : DO ispin = 1, SIZE(kp%mos, 2)
207 : CALL get_mo_set(kp%mos(1, ispin), homo=nocc, mo_coeff=mo_coeff_re, nmo=nmo, &
208 128 : eigenvalues=eigenvalues_re)
209 : CALL get_mo_set(kp%mos(2, ispin), mo_coeff=mo_coeff_im, &
210 128 : eigenvalues=eigenvalues_im)
211 128 : CPASSERT(ASSOCIATED(mo_coeff_re))
212 128 : CPASSERT(ASSOCIATED(mo_coeff_im))
213 128 : CPASSERT(ASSOCIATED(eigenvalues_re))
214 128 : CPASSERT(ASSOCIATED(eigenvalues_im))
215 128 : CPASSERT(nocc >= 1 .AND. nocc <= nmo)
216 128 : CPASSERT(SIZE(eigenvalues_re) >= nocc)
217 128 : CPASSERT(SIZE(eigenvalues_im) >= nocc)
218 :
219 : CALL cp_fm_struct_create(active_struct, template_fmstruct=mo_coeff_re%matrix_struct, &
220 128 : ncol_global=nocc)
221 128 : CALL cp_fm_create(active_re, active_struct, name="fixed-density Ritz orbitals real")
222 128 : CALL cp_fm_create(active_im, active_struct, name="fixed-density Ritz orbitals imaginary")
223 128 : CALL cp_fm_to_fm(mo_coeff_re, active_re, nocc)
224 128 : CALL cp_fm_to_fm(mo_coeff_im, active_im, nocc)
225 :
226 128 : CALL cp_cfm_create(coeff, active_struct, name="fixed-density Ritz orbitals")
227 128 : CALL cp_cfm_create(hc, active_struct, name="fixed-density Ritz HC")
228 128 : CALL cp_cfm_create(rotated, active_struct, name="fixed-density Ritz rotated orbitals")
229 : CALL cp_cfm_create(hblock, active_struct, nrow=nocc, ncol=nocc, &
230 128 : name="fixed-density projected Hamiltonian")
231 : CALL cp_cfm_create(eigenvectors, active_struct, nrow=nocc, ncol=nocc, &
232 128 : name="fixed-density Ritz eigenvectors")
233 128 : CALL cp_fm_to_cfm(kp%ot_hmat(1, ispin), kp%ot_hmat(2, ispin), matrix_h)
234 128 : CALL cp_fm_to_cfm(active_re, active_im, coeff)
235 :
236 128 : CALL cp_cfm_gemm('N', 'N', nao, nocc, nao, z_one, matrix_h, coeff, z_zero, hc)
237 128 : CALL cp_cfm_gemm('C', 'N', nocc, nocc, nao, z_one, coeff, hc, z_zero, hblock)
238 384 : ALLOCATE (ritz_values(nocc))
239 128 : CALL cp_cfm_heevd(hblock, eigenvectors, ritz_values)
240 : CALL cp_cfm_gemm('N', 'N', nao, nocc, nocc, z_one, coeff, eigenvectors, &
241 128 : z_zero, rotated)
242 128 : CALL cp_cfm_to_fm(rotated, active_re, active_im)
243 128 : CALL cp_fm_to_fm(active_re, mo_coeff_re, nocc)
244 128 : CALL cp_fm_to_fm(active_im, mo_coeff_im, nocc)
245 488 : eigenvalues_re(1:nocc) = ritz_values
246 488 : eigenvalues_im(1:nocc) = ritz_values
247 128 : DEALLOCATE (ritz_values)
248 :
249 128 : CALL cp_cfm_release(eigenvectors)
250 128 : CALL cp_cfm_release(hblock)
251 128 : CALL cp_cfm_release(rotated)
252 128 : CALL cp_cfm_release(hc)
253 128 : CALL cp_cfm_release(coeff)
254 128 : CALL cp_fm_release(active_im)
255 128 : CALL cp_fm_release(active_re)
256 620 : CALL cp_fm_struct_release(active_struct)
257 : END DO
258 :
259 260 : CALL cp_cfm_release(matrix_h)
260 : END DO
261 :
262 88 : END SUBROUTINE qs_kpoint_state_canonicalize_fixed
263 :
264 : ! **************************************************************************************************
265 : !> \brief Rebuilds the physical density from the current k-point orbitals and occupations.
266 : !> \param qs_env QS environment containing the orbital state to commit.
267 : !> \param update_occupations Recompute occupations from the current MO energies and SCF settings.
268 : !> \param separate_spin_occupations Fill restricted energy channels independently.
269 : !> \param fixed_occupations Restore the fixed occupied rank required by OT.
270 : ! **************************************************************************************************
271 592 : SUBROUTINE qs_kpoint_state_commit(qs_env, update_occupations, &
272 : separate_spin_occupations, fixed_occupations)
273 : TYPE(qs_environment_type), POINTER :: qs_env
274 : LOGICAL, INTENT(IN) :: update_occupations
275 : LOGICAL, INTENT(IN), OPTIONAL :: separate_spin_occupations, &
276 : fixed_occupations
277 :
278 : INTEGER :: ispin
279 : LOGICAL :: fixed, separate_spins
280 592 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s_kp, rho_ao_kp
281 : TYPE(kpoint_type), POINTER :: kpoints
282 592 : TYPE(mo_set_type), DIMENSION(:, :), POINTER :: mos
283 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
284 592 : POINTER :: sab_nl
285 : TYPE(qs_energy_type), POINTER :: energy
286 : TYPE(qs_rho_type), POINTER :: rho
287 : TYPE(qs_scf_env_type), POINTER :: scf_env
288 : TYPE(scf_control_type), POINTER :: scf_control
289 :
290 592 : NULLIFY (energy, kpoints, matrix_s_kp, mos, rho_ao_kp, sab_nl, rho, scf_env, scf_control)
291 :
292 592 : fixed = .FALSE.
293 592 : separate_spins = .FALSE.
294 126 : IF (PRESENT(fixed_occupations)) fixed = fixed_occupations
295 592 : IF (PRESENT(separate_spin_occupations)) separate_spins = separate_spin_occupations
296 592 : IF (fixed .AND. .NOT. update_occupations) THEN
297 0 : CPABORT("Fixed occupations can only be requested while updating occupations.")
298 : END IF
299 :
300 : CALL get_qs_env(qs_env, kpoints=kpoints, matrix_s_kp=matrix_s_kp, &
301 592 : energy=energy, rho=rho, scf_env=scf_env, scf_control=scf_control)
302 592 : CALL get_kpoint_info(kpoints, sab_nl=sab_nl)
303 :
304 592 : IF (update_occupations) THEN
305 : CALL kpoint_set_mo_occupation( &
306 566 : kpoints, scf_control%smear, separate_spin_occupations=separate_spins)
307 566 : IF (fixed) CALL qs_kpoint_set_fixed_occupations(kpoints)
308 : END IF
309 592 : CALL kpoint_density_matrices(kpoints)
310 592 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
311 : CALL kpoint_density_transform(kpoints, rho_ao_kp, .FALSE., &
312 : matrix_s_kp(1, 1)%matrix, sab_nl, scf_env%scf_work1, &
313 592 : overlap_rs=matrix_s_kp)
314 592 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
315 592 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
316 :
317 592 : energy%kTS = 0.0_dp
318 592 : energy%efermi = 0.0_dp
319 592 : mos => kpoints%kp_env(1)%kpoint_env%mos
320 592 : CPASSERT(ASSOCIATED(mos))
321 1230 : DO ispin = 1, SIZE(mos, 2)
322 638 : energy%kTS = energy%kTS + mos(1, ispin)%kTS
323 1230 : energy%efermi = energy%efermi + mos(1, ispin)%mu
324 : END DO
325 592 : energy%efermi = energy%efermi/REAL(SIZE(mos, 2), KIND=dp)
326 :
327 592 : END SUBROUTINE qs_kpoint_state_commit
328 :
329 : ! **************************************************************************************************
330 : !> \brief Checks whether every local k-point channel contains nonzero occupied orbitals.
331 : !> \param kpoints K-point environment.
332 : !> \param require_full_space Require every available orbital instead of the occupied block.
333 : !> \return True if every channel contains a nonzero real or imaginary coefficient.
334 : ! **************************************************************************************************
335 140 : FUNCTION qs_kpoint_mos_initialized(kpoints, require_full_space) RESULT(initialized)
336 : TYPE(kpoint_type), POINTER :: kpoints
337 : LOGICAL, INTENT(IN), OPTIONAL :: require_full_space
338 : LOGICAL :: initialized
339 :
340 : INTEGER :: icomponent, ikpoint_local, ispin, &
341 : nactive, nmo, nocc
342 : LOGICAL :: full_space
343 : REAL(KIND=dp) :: coefficient_norm
344 140 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: column_norm, component_norm
345 : TYPE(cp_fm_type), POINTER :: mo_coeff
346 : TYPE(kpoint_env_type), POINTER :: kp
347 : TYPE(mp_para_env_type), POINTER :: para_env
348 :
349 140 : NULLIFY (kp, mo_coeff, para_env)
350 140 : CPASSERT(ASSOCIATED(kpoints))
351 140 : CPASSERT(ASSOCIATED(kpoints%kp_env))
352 :
353 140 : full_space = .FALSE.
354 140 : IF (PRESENT(require_full_space)) full_space = require_full_space
355 140 : coefficient_norm = HUGE(1.0_dp)
356 374 : DO ikpoint_local = 1, SIZE(kpoints%kp_env)
357 234 : kp => kpoints%kp_env(ikpoint_local)%kpoint_env
358 234 : CPASSERT(ASSOCIATED(kp))
359 234 : CPASSERT(ASSOCIATED(kp%mos))
360 648 : DO ispin = 1, SIZE(kp%mos, 2)
361 274 : CALL get_mo_set(kp%mos(1, ispin), homo=nocc, mo_coeff=mo_coeff, nmo=nmo)
362 274 : CPASSERT(ASSOCIATED(mo_coeff))
363 274 : nactive = MERGE(nmo, nocc, full_space)
364 274 : IF (nactive < 1 .OR. nactive > nmo) THEN
365 0 : coefficient_norm = 0.0_dp
366 0 : CYCLE
367 : END IF
368 1096 : ALLOCATE (column_norm(nmo), component_norm(nmo))
369 274 : column_norm(:) = 0.0_dp
370 822 : DO icomponent = 1, SIZE(kp%mos, 1)
371 548 : CALL get_mo_set(kp%mos(icomponent, ispin), mo_coeff=mo_coeff)
372 548 : CPASSERT(ASSOCIATED(mo_coeff))
373 548 : CALL cp_fm_vectorsnorm(mo_coeff, component_norm)
374 3894 : column_norm(:) = column_norm(:) + component_norm(:)**2
375 : END DO
376 1790 : coefficient_norm = MIN(coefficient_norm, SQRT(MINVAL(column_norm(1:nactive))))
377 782 : DEALLOCATE (column_norm, component_norm)
378 : END DO
379 : END DO
380 :
381 140 : CALL get_kpoint_info(kpoints, para_env=para_env)
382 140 : CPASSERT(ASSOCIATED(para_env))
383 140 : CALL para_env%min(coefficient_norm)
384 140 : initialized = coefficient_norm > 100.0_dp*EPSILON(1.0_dp)
385 :
386 140 : END FUNCTION qs_kpoint_mos_initialized
387 :
388 : ! **************************************************************************************************
389 : !> \brief Copies the first k-point spin channel into the remaining channels.
390 : !> \param kpoints K-point environment.
391 : !> \param nspin Number of spin channels to populate.
392 : ! **************************************************************************************************
393 10 : SUBROUTINE qs_kpoint_copy_spin_mos(kpoints, nspin)
394 : TYPE(kpoint_type), POINTER :: kpoints
395 : INTEGER, INTENT(IN) :: nspin
396 :
397 : INTEGER :: icomponent, ikpoint_local, ispin
398 10 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues_source, eigenvalues_target
399 : TYPE(cp_fm_type), POINTER :: mo_coeff_source, mo_coeff_target
400 : TYPE(kpoint_env_type), POINTER :: kp
401 :
402 10 : NULLIFY (eigenvalues_source, eigenvalues_target, kp, mo_coeff_source, mo_coeff_target)
403 0 : CPASSERT(ASSOCIATED(kpoints))
404 10 : CPASSERT(ASSOCIATED(kpoints%kp_env))
405 10 : CPASSERT(nspin >= 1)
406 :
407 24 : DO ikpoint_local = 1, SIZE(kpoints%kp_env)
408 14 : kp => kpoints%kp_env(ikpoint_local)%kpoint_env
409 14 : CPASSERT(ASSOCIATED(kp))
410 14 : CPASSERT(ASSOCIATED(kp%mos))
411 14 : CPASSERT(SIZE(kp%mos, 2) >= nspin)
412 52 : DO icomponent = 1, SIZE(kp%mos, 1)
413 : CALL get_mo_set(kp%mos(icomponent, 1), mo_coeff=mo_coeff_source, &
414 28 : eigenvalues=eigenvalues_source)
415 28 : CPASSERT(ASSOCIATED(mo_coeff_source))
416 28 : CPASSERT(ASSOCIATED(eigenvalues_source))
417 70 : DO ispin = 2, nspin
418 : CALL get_mo_set(kp%mos(icomponent, ispin), mo_coeff=mo_coeff_target, &
419 28 : eigenvalues=eigenvalues_target)
420 28 : CPASSERT(ASSOCIATED(mo_coeff_target))
421 28 : CPASSERT(ASSOCIATED(eigenvalues_target))
422 28 : CALL cp_fm_to_fm(mo_coeff_source, mo_coeff_target)
423 388 : eigenvalues_target(:) = eigenvalues_source(:)
424 : END DO
425 : END DO
426 : END DO
427 :
428 10 : END SUBROUTINE qs_kpoint_copy_spin_mos
429 :
430 : ! **************************************************************************************************
431 : !> \brief Restores the fixed occupied rank required by k-point OT.
432 : !> \param kpoints K-point environment.
433 : ! **************************************************************************************************
434 44 : SUBROUTINE qs_kpoint_set_fixed_occupations(kpoints)
435 : TYPE(kpoint_type), POINTER :: kpoints
436 :
437 : INTEGER :: icomponent, ikpoint_local, ispin, &
438 : nelectron, nmo, nocc
439 : REAL(KIND=dp) :: maxocc
440 44 : REAL(KIND=dp), DIMENSION(:), POINTER :: occupation_numbers
441 : TYPE(kpoint_env_type), POINTER :: kp
442 :
443 44 : NULLIFY (kp, occupation_numbers)
444 0 : CPASSERT(ASSOCIATED(kpoints))
445 44 : CPASSERT(ASSOCIATED(kpoints%kp_env))
446 :
447 152 : DO ikpoint_local = 1, SIZE(kpoints%kp_env)
448 108 : kp => kpoints%kp_env(ikpoint_local)%kpoint_env
449 108 : CPASSERT(ASSOCIATED(kp))
450 108 : CPASSERT(ASSOCIATED(kp%mos))
451 280 : DO ispin = 1, SIZE(kp%mos, 2)
452 492 : DO icomponent = 1, SIZE(kp%mos, 1)
453 : CALL get_mo_set(kp%mos(icomponent, ispin), homo=nocc, maxocc=maxocc, &
454 256 : nelectron=nelectron, nmo=nmo, occupation_numbers=occupation_numbers)
455 256 : CPASSERT(ASSOCIATED(occupation_numbers))
456 256 : IF (nocc < 1 .OR. nocc > nmo) THEN
457 0 : CPABORT("Fixed-occupation k-point OT requires a nonempty occupied subspace.")
458 : END IF
459 256 : IF (ABS(maxocc*REAL(nocc, KIND=dp) - REAL(nelectron, KIND=dp)) > &
460 : 100.0_dp*EPSILON(1.0_dp)) THEN
461 0 : CPABORT("Fixed-occupation k-point OT requires uniform integer occupations.")
462 : END IF
463 1016 : occupation_numbers(:) = 0.0_dp
464 976 : occupation_numbers(1:nocc) = maxocc
465 640 : kp%mos(icomponent, ispin)%kTS = 0.0_dp
466 : END DO
467 : END DO
468 : END DO
469 :
470 44 : END SUBROUTINE qs_kpoint_set_fixed_occupations
471 :
472 : END MODULE qs_kpoint_state
|
