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 methods for deltaSCF calculations
10 : ! **************************************************************************************************
11 : MODULE qs_mom_methods
12 : USE bibliography, ONLY: Barca2018,&
13 : Gilbert2008,&
14 : cite_reference
15 : USE cp_blacs_env, ONLY: cp_blacs_env_type
16 : USE cp_dbcsr_api, ONLY: dbcsr_p_type
17 : USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply
18 : USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale
19 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
20 : cp_fm_struct_release,&
21 : cp_fm_struct_type
22 : USE cp_fm_types, ONLY: cp_fm_create,&
23 : cp_fm_get_info,&
24 : cp_fm_maxabsval,&
25 : cp_fm_to_fm,&
26 : cp_fm_type,&
27 : cp_fm_vectorsnorm,&
28 : cp_fm_vectorssum
29 : USE input_constants, ONLY: momproj_norm,&
30 : momproj_sum,&
31 : momtype_imom,&
32 : momtype_mom
33 : USE input_section_types, ONLY: section_vals_type
34 : USE kinds, ONLY: dp
35 : USE parallel_gemm_api, ONLY: parallel_gemm
36 : USE qs_density_matrices, ONLY: calculate_density_matrix
37 : USE qs_mo_types, ONLY: get_mo_set,&
38 : mo_set_type,&
39 : set_mo_set
40 : USE qs_scf_diagonalization, ONLY: general_eigenproblem
41 : USE qs_scf_types, ONLY: qs_scf_env_type
42 : USE scf_control_types, ONLY: scf_control_type
43 : USE string_utilities, ONLY: integer_to_string
44 : USE util, ONLY: sort,&
45 : sort_unique
46 : #include "./base/base_uses.f90"
47 :
48 : IMPLICIT NONE
49 :
50 : PRIVATE
51 :
52 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_mom_methods'
53 :
54 : PUBLIC :: do_mom_guess, do_mom_diag
55 : PRIVATE :: mom_is_unique_orbital_indices, mom_reoccupy_orbitals
56 :
57 : CONTAINS
58 :
59 : ! **************************************************************************************************
60 : !> \brief check that every molecular orbital index appears only once in each
61 : !> (de-)occupation list supplied by user. Check that all the indices
62 : !> are positive integers and abort if it is not the case.
63 : !> \param iarr list of molecular orbital indices to be checked
64 : !> \return .true. if all the elements are unique or the list contains
65 : !> exactly one 0 element (meaning no excitation)
66 : !> \par History
67 : !> 01.2016 created [Sergey Chulkov]
68 : ! **************************************************************************************************
69 80 : FUNCTION mom_is_unique_orbital_indices(iarr) RESULT(is_unique)
70 : INTEGER, DIMENSION(:), POINTER :: iarr
71 : LOGICAL :: is_unique
72 :
73 : CHARACTER(len=*), PARAMETER :: routineN = 'mom_is_unique_orbital_indices'
74 :
75 : INTEGER :: handle, norbs
76 80 : INTEGER, DIMENSION(:), POINTER :: tmp_iarr
77 :
78 80 : CALL timeset(routineN, handle)
79 :
80 80 : CPASSERT(ASSOCIATED(iarr))
81 80 : norbs = SIZE(iarr)
82 :
83 80 : IF (norbs > 0) THEN
84 240 : ALLOCATE (tmp_iarr(norbs))
85 :
86 320 : tmp_iarr(:) = iarr(:)
87 80 : CALL sort_unique(tmp_iarr, is_unique)
88 :
89 : ! Ensure that all orbital indices are positive integers.
90 : ! A special value '0' means 'disabled keyword',
91 : ! it must appear once to be interpreted in such a way
92 80 : IF (tmp_iarr(1) < 0 .OR. (tmp_iarr(1) == 0 .AND. norbs > 1)) &
93 0 : CPABORT("MOM: all molecular orbital indices must be positive integer numbers")
94 :
95 160 : DEALLOCATE (tmp_iarr)
96 : END IF
97 :
98 : is_unique = .TRUE.
99 :
100 80 : CALL timestop(handle)
101 :
102 80 : END FUNCTION mom_is_unique_orbital_indices
103 :
104 : ! **************************************************************************************************
105 : !> \brief swap occupation numbers between molecular orbitals
106 : !> from occupation and de-occupation lists
107 : !> \param mo_set set of molecular orbitals
108 : !> \param deocc_orb_set list of de-occupied orbital indices
109 : !> \param occ_orb_set list of newly occupied orbital indices
110 : !> \param spin spin component of the molecular orbitals;
111 : !> to be used for diagnostic messages
112 : !> \par History
113 : !> 01.2016 created [Sergey Chulkov]
114 : ! **************************************************************************************************
115 40 : SUBROUTINE mom_reoccupy_orbitals(mo_set, deocc_orb_set, occ_orb_set, spin)
116 : TYPE(mo_set_type), INTENT(INOUT) :: mo_set
117 : INTEGER, DIMENSION(:), POINTER :: deocc_orb_set, occ_orb_set
118 : CHARACTER(len=*), INTENT(in) :: spin
119 :
120 : CHARACTER(len=*), PARAMETER :: routineN = 'mom_reoccupy_orbitals'
121 :
122 : CHARACTER(len=10) :: str_iorb, str_norbs
123 : CHARACTER(len=3) :: str_prefix
124 : INTEGER :: handle, homo, iorb, lfomo, nao, nmo, &
125 : norbs
126 : REAL(kind=dp) :: maxocc
127 40 : REAL(kind=dp), DIMENSION(:), POINTER :: occ_nums
128 :
129 40 : CALL timeset(routineN, handle)
130 :
131 : ! MOM electron excitation should preserve both the number of electrons and
132 : ! multiplicity of the electronic system thus ensuring the following constraint :
133 : ! norbs = SIZE(deocc_orb_set) == SIZE(occ_orb_set)
134 40 : norbs = SIZE(deocc_orb_set)
135 :
136 : ! the following assertion should never raise an exception
137 40 : CPASSERT(SIZE(deocc_orb_set) == SIZE(occ_orb_set))
138 :
139 : ! MOM does not follow aufbau principle producing non-uniformly occupied orbitals
140 40 : CALL set_mo_set(mo_set=mo_set, uniform_occupation=.FALSE.)
141 :
142 40 : IF (deocc_orb_set(1) /= 0 .AND. occ_orb_set(1) /= 0) THEN
143 : CALL get_mo_set(mo_set=mo_set, maxocc=maxocc, &
144 20 : nao=nao, nmo=nmo, occupation_numbers=occ_nums)
145 :
146 20 : IF (deocc_orb_set(norbs) > nao .OR. occ_orb_set(norbs) > nao) THEN
147 : ! STOP: one of the molecular orbital index exceeds the number of atomic basis functions available
148 0 : CALL integer_to_string(nao, str_norbs)
149 :
150 0 : IF (deocc_orb_set(norbs) >= occ_orb_set(norbs)) THEN
151 0 : iorb = deocc_orb_set(norbs)
152 0 : str_prefix = 'de-'
153 : ELSE
154 0 : iorb = occ_orb_set(norbs)
155 0 : str_prefix = ''
156 : END IF
157 0 : CALL integer_to_string(iorb, str_iorb)
158 :
159 : CALL cp_abort(__LOCATION__, "Unable to "//TRIM(str_prefix)//"occupy "// &
160 : TRIM(spin)//" orbital No. "//TRIM(str_iorb)// &
161 : " since its index exceeds the number of atomic orbital functions available ("// &
162 0 : TRIM(str_norbs)//"). Please consider using a larger basis set.")
163 : END IF
164 :
165 20 : IF (deocc_orb_set(norbs) > nmo .OR. occ_orb_set(norbs) > nmo) THEN
166 : ! STOP: one of the molecular orbital index exceeds the number of constructed molecular orbitals
167 0 : IF (deocc_orb_set(norbs) >= occ_orb_set(norbs)) THEN
168 0 : iorb = deocc_orb_set(norbs)
169 : ELSE
170 0 : iorb = occ_orb_set(norbs)
171 : END IF
172 :
173 0 : IF (iorb - nmo > 1) THEN
174 0 : CALL integer_to_string(iorb - nmo, str_iorb)
175 0 : str_prefix = 's'
176 : ELSE
177 0 : str_iorb = 'an'
178 0 : str_prefix = ''
179 : END IF
180 :
181 0 : CALL integer_to_string(nmo, str_norbs)
182 :
183 : CALL cp_abort(__LOCATION__, "The number of molecular orbitals ("//TRIM(str_norbs)// &
184 : ") is not enough to perform MOM calculation. Please add "// &
185 : TRIM(str_iorb)//" extra orbital"//TRIM(str_prefix)// &
186 0 : " using the ADDED_MOS keyword in the SCF section of your input file.")
187 : END IF
188 :
189 40 : DO iorb = 1, norbs
190 : ! swap occupation numbers between two adjoint molecular orbitals
191 20 : IF (occ_nums(deocc_orb_set(iorb)) <= 0.0_dp) THEN
192 0 : CALL integer_to_string(deocc_orb_set(iorb), str_iorb)
193 :
194 : CALL cp_abort(__LOCATION__, "The "//TRIM(spin)//" orbital No. "// &
195 0 : TRIM(str_iorb)//" is not occupied thus it cannot be deoccupied.")
196 : END IF
197 :
198 20 : IF (occ_nums(occ_orb_set(iorb)) > 0.0_dp) THEN
199 0 : CALL integer_to_string(occ_orb_set(iorb), str_iorb)
200 :
201 : CALL cp_abort(__LOCATION__, "The "//TRIM(spin)//" orbital No. "// &
202 0 : TRIM(str_iorb)//" is already occupied thus it cannot be reoccupied.")
203 : END IF
204 :
205 20 : occ_nums(occ_orb_set(iorb)) = occ_nums(deocc_orb_set(iorb))
206 40 : occ_nums(deocc_orb_set(iorb)) = 0.0_dp
207 : END DO
208 :
209 : ! locate the lowest non-maxocc occupied orbital
210 78 : DO lfomo = 1, nmo
211 78 : IF (occ_nums(lfomo) /= maxocc) EXIT
212 : END DO
213 :
214 : ! locate the highest occupied orbital
215 90 : DO homo = nmo, 1, -1
216 90 : IF (occ_nums(homo) > 0.0_dp) EXIT
217 : END DO
218 :
219 20 : CALL set_mo_set(mo_set=mo_set, homo=homo, lfomo=lfomo)
220 :
221 20 : ELSE IF (deocc_orb_set(1) /= 0 .OR. occ_orb_set(1) /= 0) THEN
222 : CALL cp_abort(__LOCATION__, &
223 0 : "Incorrect multiplicity of the MOM reference electronic state")
224 : END IF
225 :
226 40 : CALL timestop(handle)
227 :
228 40 : END SUBROUTINE mom_reoccupy_orbitals
229 :
230 : ! **************************************************************************************************
231 : !> \brief initial guess for the maximum overlap method
232 : !> \param nspins number of spin components
233 : !> \param mos array of molecular orbitals
234 : !> \param scf_control SCF control variables
235 : !> \param p_rmpv density matrix to be computed
236 : !> \par History
237 : !> * 01.2016 created [Sergey Chulkov]
238 : ! **************************************************************************************************
239 20 : SUBROUTINE do_mom_guess(nspins, mos, scf_control, p_rmpv)
240 : INTEGER, INTENT(in) :: nspins
241 : TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
242 : TYPE(scf_control_type), POINTER :: scf_control
243 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_rmpv
244 :
245 : CHARACTER(len=*), PARAMETER :: routineN = 'do_mom_guess'
246 :
247 : CHARACTER(len=10) :: str_iter
248 : INTEGER :: handle, ispin, scf_iter
249 : LOGICAL :: is_mo
250 : REAL(kind=dp) :: maxa
251 : TYPE(cp_fm_type), POINTER :: mo_coeff
252 :
253 20 : CALL timeset(routineN, handle)
254 :
255 : ! we are about to initialise the maximum overlap method,
256 : ! so cite the relevant reference first
257 20 : IF (scf_control%diagonalization%mom_type == momtype_mom) THEN
258 20 : CALL cite_reference(Gilbert2008)
259 0 : ELSE IF (scf_control%diagonalization%mom_type == momtype_imom) THEN
260 0 : CALL cite_reference(Barca2018)
261 : END IF
262 :
263 : ! ensure we do not have duplicated orbital indices
264 20 : IF (.NOT. &
265 : (mom_is_unique_orbital_indices(scf_control%diagonalization%mom_deoccA) .AND. &
266 : mom_is_unique_orbital_indices(scf_control%diagonalization%mom_deoccB) .AND. &
267 : mom_is_unique_orbital_indices(scf_control%diagonalization%mom_occA) .AND. &
268 : mom_is_unique_orbital_indices(scf_control%diagonalization%mom_occB))) &
269 : CALL cp_abort(__LOCATION__, &
270 0 : "Duplicate orbital indices were found in the MOM section")
271 :
272 : ! ignore beta orbitals for spin-unpolarized calculations
273 20 : IF (nspins == 1 .AND. (scf_control%diagonalization%mom_deoccB(1) /= 0 &
274 : .OR. scf_control%diagonalization%mom_occB(1) /= 0)) THEN
275 :
276 : CALL cp_warn(__LOCATION__, "Maximum overlap method will"// &
277 0 : " ignore beta orbitals since neither UKS nor ROKS calculation is performed")
278 : END IF
279 :
280 : ! compute the change in multiplicity and number of electrons
281 : IF (SIZE(scf_control%diagonalization%mom_deoccA) /= &
282 20 : SIZE(scf_control%diagonalization%mom_occA) .OR. &
283 : (nspins > 1 .AND. &
284 : SIZE(scf_control%diagonalization%mom_deoccB) /= &
285 : SIZE(scf_control%diagonalization%mom_occB))) THEN
286 :
287 : CALL cp_abort(__LOCATION__, "Incorrect multiplicity of the MOM reference"// &
288 0 : " electronic state or inconsistent number of electrons")
289 : END IF
290 :
291 20 : is_mo = .FALSE.
292 : ! by default activate MOM at the second SCF iteration as the
293 : ! 'old' molecular orbitals are unavailable from the very beginning
294 20 : scf_iter = 2
295 : ! check if the molecular orbitals are actually there
296 : ! by finding at least one MO coefficient > 0
297 28 : DO ispin = 1, nspins
298 24 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
299 24 : CALL cp_fm_maxabsval(mo_coeff, maxa)
300 : ! is_mo |= maxa > 0.0_dp
301 28 : IF (maxa > 0.0_dp) THEN
302 16 : is_mo = .TRUE.
303 : ! we already have the molecular orbitals (e.g. from a restart file);
304 : ! activate MOM immediately if the input keyword START_ITER is not given
305 16 : scf_iter = 1
306 16 : EXIT
307 : END IF
308 : END DO
309 :
310 : ! proceed alpha orbitals
311 20 : IF (nspins >= 1) &
312 : CALL mom_reoccupy_orbitals(mos(1), &
313 : scf_control%diagonalization%mom_deoccA, &
314 20 : scf_control%diagonalization%mom_occA, 'alpha')
315 :
316 : ! proceed beta orbitals (if any)
317 20 : IF (nspins >= 2) &
318 : CALL mom_reoccupy_orbitals(mos(2), &
319 : scf_control%diagonalization%mom_deoccB, &
320 20 : scf_control%diagonalization%mom_occB, 'beta')
321 :
322 : ! recompute the density matrix if the molecular orbitals are here;
323 : ! otherwise do nothing to prevent zeroing out the density matrix
324 : ! obtained from atomic guess
325 20 : IF (is_mo) THEN
326 48 : DO ispin = 1, nspins
327 48 : CALL calculate_density_matrix(mos(ispin), p_rmpv(ispin)%matrix)
328 : END DO
329 : END IF
330 :
331 : ! adjust the start SCF iteration number if needed
332 20 : IF (scf_control%diagonalization%mom_start < scf_iter) THEN
333 18 : IF (scf_control%diagonalization%mom_start > 0) THEN
334 : ! inappropriate iteration number has been provided through the input file;
335 : ! fix it and issue a warning message
336 0 : CALL integer_to_string(scf_iter, str_iter)
337 : CALL cp_warn(__LOCATION__, &
338 : "The maximum overlap method will be activated at the SCF iteration No. "// &
339 0 : TRIM(str_iter)//" due to the SCF guess method used.")
340 : END IF
341 18 : scf_control%diagonalization%mom_start = scf_iter
342 2 : ELSE IF (scf_control%diagonalization%mom_start > scf_iter .AND. &
343 : (scf_control%diagonalization%mom_occA(1) > 0 .OR. scf_control%diagonalization%mom_occB(1) > 0)) THEN
344 : ! the keyword START_ITER has been provided for an excited state calculation, ignore it
345 2 : CALL integer_to_string(scf_iter, str_iter)
346 : CALL cp_warn(__LOCATION__, &
347 : "The maximum overlap method will be activated at the SCF iteration No. "// &
348 2 : TRIM(str_iter)//" because an excited state calculation has been requested")
349 2 : scf_control%diagonalization%mom_start = scf_iter
350 : END IF
351 :
352 : ! MOM is now initialised properly
353 20 : scf_control%diagonalization%mom_didguess = .TRUE.
354 :
355 20 : CALL timestop(handle)
356 :
357 20 : END SUBROUTINE do_mom_guess
358 :
359 : ! **************************************************************************************************
360 : !> \brief do an SCF iteration, then compute occupation numbers of the new
361 : !> molecular orbitals according to their overlap with the previous ones
362 : !> \param scf_env SCF environment information
363 : !> \param mos array of molecular orbitals
364 : !> \param matrix_ks sparse Kohn-Sham matrix
365 : !> \param matrix_s sparse overlap matrix
366 : !> \param scf_control SCF control variables
367 : !> \param scf_section SCF input section
368 : !> \param diis_step have we done a DIIS step
369 : !> \par History
370 : !> * 07.2014 created [Matt Watkins]
371 : !> * 01.2016 release version [Sergey Chulkov]
372 : !> * 03.2018 initial maximum overlap method [Sergey Chulkov]
373 : ! **************************************************************************************************
374 324 : SUBROUTINE do_mom_diag(scf_env, mos, matrix_ks, matrix_s, scf_control, scf_section, diis_step)
375 : TYPE(qs_scf_env_type), POINTER :: scf_env
376 : TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
377 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
378 : TYPE(scf_control_type), POINTER :: scf_control
379 : TYPE(section_vals_type), POINTER :: scf_section
380 : LOGICAL, INTENT(INOUT) :: diis_step
381 :
382 : CHARACTER(len=*), PARAMETER :: routineN = 'do_mom_diag'
383 :
384 : INTEGER :: handle, homo, iproj, ispin, lfomo, nao, &
385 : nmo, nspins
386 324 : INTEGER, ALLOCATABLE, DIMENSION(:) :: inds
387 : REAL(kind=dp) :: maxocc
388 324 : REAL(kind=dp), DIMENSION(:), POINTER :: occ_nums, proj, tmp_occ_nums
389 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
390 : TYPE(cp_fm_struct_type), POINTER :: ao_mo_fmstruct, mo_mo_fmstruct
391 : TYPE(cp_fm_type), POINTER :: mo_coeff, mo_coeff_ref, overlap, svec
392 :
393 324 : CALL timeset(routineN, handle)
394 :
395 324 : IF (.NOT. scf_control%diagonalization%mom_didguess) &
396 : CALL cp_abort(__LOCATION__, &
397 0 : "The current implementation of the maximum overlap method is incompatible with the initial SCF guess")
398 :
399 : ! number of spins == dft_control%nspins
400 324 : nspins = SIZE(matrix_ks)
401 :
402 : ! copy old molecular orbitals
403 324 : IF (scf_env%iter_count >= scf_control%diagonalization%mom_start) THEN
404 318 : IF (.NOT. ASSOCIATED(scf_env%mom_ref_mo_coeff)) THEN
405 110 : ALLOCATE (scf_env%mom_ref_mo_coeff(nspins))
406 66 : DO ispin = 1, nspins
407 44 : NULLIFY (ao_mo_fmstruct)
408 44 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, occupation_numbers=occ_nums)
409 44 : CALL cp_fm_get_info(mo_coeff, matrix_struct=ao_mo_fmstruct)
410 44 : CALL cp_fm_create(scf_env%mom_ref_mo_coeff(ispin), ao_mo_fmstruct)
411 :
412 : ! Initial Maximum Overlap Method: keep initial molecular orbitals
413 66 : IF (scf_control%diagonalization%mom_type == momtype_imom) THEN
414 0 : CALL cp_fm_to_fm(mo_coeff, scf_env%mom_ref_mo_coeff(ispin))
415 0 : CALL cp_fm_column_scale(scf_env%mom_ref_mo_coeff(ispin), occ_nums)
416 : END IF
417 : END DO
418 : END IF
419 :
420 : ! allocate the molecular orbitals overlap matrix
421 318 : IF (.NOT. ASSOCIATED(scf_env%mom_overlap)) THEN
422 100 : ALLOCATE (scf_env%mom_overlap(nspins))
423 60 : DO ispin = 1, nspins
424 40 : NULLIFY (blacs_env, mo_mo_fmstruct)
425 40 : CALL get_mo_set(mo_set=mos(ispin), nmo=nmo, mo_coeff=mo_coeff)
426 40 : CALL cp_fm_get_info(mo_coeff, context=blacs_env)
427 40 : CALL cp_fm_struct_create(mo_mo_fmstruct, nrow_global=nmo, ncol_global=nmo, context=blacs_env)
428 40 : CALL cp_fm_create(scf_env%mom_overlap(ispin), mo_mo_fmstruct)
429 100 : CALL cp_fm_struct_release(mo_mo_fmstruct)
430 : END DO
431 : END IF
432 :
433 : ! allocate a matrix to store the product S * mo_coeff
434 318 : IF (.NOT. ASSOCIATED(scf_env%mom_s_mo_coeff)) THEN
435 100 : ALLOCATE (scf_env%mom_s_mo_coeff(nspins))
436 60 : DO ispin = 1, nspins
437 40 : NULLIFY (ao_mo_fmstruct)
438 40 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
439 40 : CALL cp_fm_get_info(mo_coeff, matrix_struct=ao_mo_fmstruct)
440 60 : CALL cp_fm_create(scf_env%mom_s_mo_coeff(ispin), ao_mo_fmstruct)
441 : END DO
442 : END IF
443 :
444 : ! Original Maximum Overlap Method: keep orbitals from the previous SCF iteration
445 318 : IF (scf_control%diagonalization%mom_type == momtype_mom) THEN
446 954 : DO ispin = 1, nspins
447 636 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, occupation_numbers=occ_nums)
448 636 : CALL cp_fm_to_fm(mo_coeff, scf_env%mom_ref_mo_coeff(ispin))
449 954 : CALL cp_fm_column_scale(scf_env%mom_ref_mo_coeff(ispin), occ_nums)
450 : END DO
451 : END IF
452 : END IF
453 :
454 : ! solve the eigenproblem
455 324 : CALL general_eigenproblem(scf_env, mos, matrix_ks, matrix_s, scf_control, scf_section, diis_step)
456 :
457 324 : IF (scf_env%iter_count >= scf_control%diagonalization%mom_start) THEN
458 954 : DO ispin = 1, nspins
459 :
460 : ! TO DO: sparse-matrix variant; check if use_mo_coeff_b is set, and if yes use mo_coeff_b instead
461 : CALL get_mo_set(mo_set=mos(ispin), maxocc=maxocc, mo_coeff=mo_coeff, &
462 636 : nao=nao, nmo=nmo, occupation_numbers=occ_nums)
463 :
464 636 : mo_coeff_ref => scf_env%mom_ref_mo_coeff(ispin)
465 636 : overlap => scf_env%mom_overlap(ispin)
466 636 : svec => scf_env%mom_s_mo_coeff(ispin)
467 :
468 : ! svec = S * C(new)
469 636 : CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, mo_coeff, svec, nmo)
470 :
471 : ! overlap = C(reference occupied)^T * S * C(new)
472 636 : CALL parallel_gemm('T', 'N', nmo, nmo, nao, 1.0_dp, mo_coeff_ref, svec, 0.0_dp, overlap)
473 :
474 1908 : ALLOCATE (proj(nmo))
475 1908 : ALLOCATE (inds(nmo))
476 1272 : ALLOCATE (tmp_occ_nums(nmo))
477 :
478 : ! project the new molecular orbitals into the space of the reference occupied orbitals
479 636 : SELECT CASE (scf_control%diagonalization%mom_proj_formula)
480 : CASE (momproj_sum)
481 : ! proj_j = abs( \sum_i overlap(i, j) )
482 0 : CALL cp_fm_vectorssum(overlap, proj)
483 :
484 0 : DO iproj = 1, nmo
485 0 : proj(iproj) = ABS(proj(iproj))
486 : END DO
487 :
488 : CASE (momproj_norm)
489 : ! proj_j = (\sum_i overlap(i, j)**2) ** 0.5
490 636 : CALL cp_fm_vectorsnorm(overlap, proj)
491 :
492 : CASE DEFAULT
493 636 : CPABORT("Unimplemented projection formula")
494 : END SELECT
495 :
496 11688 : tmp_occ_nums(:) = occ_nums(:)
497 : ! sort occupation numbers in ascending order
498 636 : CALL sort(tmp_occ_nums, nmo, inds)
499 : ! sort overlap projection in ascending order
500 636 : CALL sort(proj, nmo, inds)
501 :
502 : ! reorder occupation numbers according to overlap projections
503 5844 : DO iproj = 1, nmo
504 5844 : occ_nums(inds(iproj)) = tmp_occ_nums(iproj)
505 : END DO
506 :
507 636 : DEALLOCATE (tmp_occ_nums)
508 636 : DEALLOCATE (inds)
509 636 : DEALLOCATE (proj)
510 :
511 : ! locate the lowest non-fully occupied orbital
512 2882 : DO lfomo = 1, nmo
513 2882 : IF (occ_nums(lfomo) /= maxocc) EXIT
514 : END DO
515 :
516 : ! locate the highest occupied orbital
517 2918 : DO homo = nmo, 1, -1
518 2918 : IF (occ_nums(homo) > 0.0_dp) EXIT
519 : END DO
520 :
521 1590 : CALL set_mo_set(mo_set=mos(ispin), homo=homo, lfomo=lfomo)
522 : END DO
523 : END IF
524 :
525 : ! recompute density matrix
526 972 : DO ispin = 1, nspins
527 972 : CALL calculate_density_matrix(mos(ispin), scf_env%p_mix_new(ispin, 1)%matrix)
528 : END DO
529 :
530 324 : CALL timestop(handle)
531 :
532 648 : END SUBROUTINE do_mom_diag
533 :
534 : END MODULE qs_mom_methods
|
