(git:374b731)
Loading...
Searching...
No Matches
qs_loc_utils.F
Go to the documentation of this file.
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 Some utilities for the construction of
10!> the localization environment
11!> \author MI (05-2005)
12! **************************************************************************************************
14
15 USE ai_moments, ONLY: contract_cossin,&
16 cossin
20 USE cell_types, ONLY: cell_type,&
21 pbc
25 USE cp_files, ONLY: close_file,&
32 USE cp_fm_types, ONLY: &
39 USE cp_output_handling, ONLY: cp_p_file,&
44 USE dbcsr_api, ONLY: dbcsr_copy,&
45 dbcsr_get_block_p,&
46 dbcsr_p_type,&
47 dbcsr_set,&
48 dbcsr_type
50 USE input_constants, ONLY: &
57 USE kinds, ONLY: default_path_length,&
59 dp
60 USE mathconstants, ONLY: twopi
63 USE orbital_pointers, ONLY: ncoset
68 USE qs_kind_types, ONLY: get_qs_kind,&
71 USE qs_loc_types, ONLY: get_qs_loc_env,&
78 USE qs_mo_methods, ONLY: make_mo_eig
79 USE qs_mo_types, ONLY: get_mo_set,&
87 USE qs_scf_types, ONLY: ot_method_nr
89#include "./base/base_uses.f90"
90
91 IMPLICIT NONE
92
93 PRIVATE
94
95! *** Global parameters ***
96
97 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_loc_utils'
98
99! *** Public ***
103
104CONTAINS
105
106! **************************************************************************************************
107!> \brief copy old mos to new ones, allocating as necessary
108!> \param mo_loc_history ...
109!> \param mo_loc ...
110! **************************************************************************************************
111 SUBROUTINE retain_history(mo_loc_history, mo_loc)
112
113 TYPE(cp_fm_type), DIMENSION(:), POINTER :: mo_loc_history
114 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_loc
115
116 CHARACTER(len=*), PARAMETER :: routinen = 'retain_history'
117
118 INTEGER :: handle, i, ncol_hist, ncol_loc
119
120 CALL timeset(routinen, handle)
121
122 IF (.NOT. ASSOCIATED(mo_loc_history)) THEN
123 ALLOCATE (mo_loc_history(SIZE(mo_loc)))
124 DO i = 1, SIZE(mo_loc_history)
125 CALL cp_fm_create(mo_loc_history(i), mo_loc(i)%matrix_struct)
126 END DO
127 END IF
128
129 DO i = 1, SIZE(mo_loc_history)
130 CALL cp_fm_get_info(mo_loc_history(i), ncol_global=ncol_hist)
131 CALL cp_fm_get_info(mo_loc(i), ncol_global=ncol_loc)
132 cpassert(ncol_hist == ncol_loc)
133 CALL cp_fm_to_fm(mo_loc(i), mo_loc_history(i))
134 END DO
135
136 CALL timestop(handle)
137
138 END SUBROUTINE retain_history
139
140! **************************************************************************************************
141!> \brief rotate the mo_new, so that the orbitals are as similar
142!> as possible to ones in mo_ref.
143!> \param mo_new ...
144!> \param mo_ref ...
145!> \param matrix_S ...
146! **************************************************************************************************
147 SUBROUTINE rotate_state_to_ref(mo_new, mo_ref, matrix_S)
148
149 TYPE(cp_fm_type), INTENT(IN) :: mo_new, mo_ref
150 TYPE(dbcsr_type), POINTER :: matrix_s
151
152 CHARACTER(len=*), PARAMETER :: routinen = 'rotate_state_to_ref'
153
154 INTEGER :: handle, ncol, ncol_ref, nrow
155 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues
156 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
157 TYPE(cp_fm_type) :: o1, o2, o3, o4, smo
158
159 CALL timeset(routinen, handle)
160
161 CALL cp_fm_get_info(mo_new, nrow_global=nrow, ncol_global=ncol)
162 CALL cp_fm_get_info(mo_ref, ncol_global=ncol_ref)
163 cpassert(ncol == ncol_ref)
164
165 NULLIFY (fm_struct_tmp)
166 CALL cp_fm_create(smo, mo_ref%matrix_struct)
167
168 CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ncol, &
169 ncol_global=ncol, para_env=mo_new%matrix_struct%para_env, &
170 context=mo_new%matrix_struct%context)
171 CALL cp_fm_create(o1, fm_struct_tmp)
172 CALL cp_fm_create(o2, fm_struct_tmp)
173 CALL cp_fm_create(o3, fm_struct_tmp)
174 CALL cp_fm_create(o4, fm_struct_tmp)
175 CALL cp_fm_struct_release(fm_struct_tmp)
176
177 ! o1 = (mo_new)^T matrix_S mo_ref
178 CALL cp_dbcsr_sm_fm_multiply(matrix_s, mo_ref, smo, ncol)
179 CALL parallel_gemm('T', 'N', ncol, ncol, nrow, 1.0_dp, mo_new, smo, 0.0_dp, o1)
180
181 ! o2 = (o1^T o1)
182 CALL parallel_gemm('T', 'N', ncol, ncol, ncol, 1.0_dp, o1, o1, 0.0_dp, o2)
183
184 ! o2 = (o1^T o1)^-1/2
185 ALLOCATE (eigenvalues(ncol))
186 CALL choose_eigv_solver(o2, o3, eigenvalues)
187 CALL cp_fm_to_fm(o3, o4)
188 eigenvalues(:) = 1.0_dp/sqrt(eigenvalues(:))
189 CALL cp_fm_column_scale(o4, eigenvalues)
190 CALL parallel_gemm('N', 'T', ncol, ncol, ncol, 1.0_dp, o3, o4, 0.0_dp, o2)
191
192 ! o3 = o1 (o1^T o1)^-1/2
193 CALL parallel_gemm('N', 'N', ncol, ncol, ncol, 1.0_dp, o1, o2, 0.0_dp, o3)
194
195 ! mo_new o1 (o1^T o1)^-1/2
196 CALL parallel_gemm('N', 'N', nrow, ncol, ncol, 1.0_dp, mo_new, o3, 0.0_dp, smo)
197 CALL cp_fm_to_fm(smo, mo_new)
198
199 ! XXXXXXX testing
200 ! CALL parallel_gemm('N','T',ncol,ncol,ncol,1.0_dp,o3,o3,0.0_dp,o1)
201 ! WRITE(*,*) o1%local_data
202 ! CALL parallel_gemm('T','N',ncol,ncol,ncol,1.0_dp,o3,o3,0.0_dp,o1)
203 ! WRITE(*,*) o1%local_data
204
205 CALL cp_fm_release(o1)
206 CALL cp_fm_release(o2)
207 CALL cp_fm_release(o3)
208 CALL cp_fm_release(o4)
209 CALL cp_fm_release(smo)
210
211 CALL timestop(handle)
212
213 END SUBROUTINE rotate_state_to_ref
214
215! **************************************************************************************************
216!> \brief allocates the data, and initializes the operators
217!> \param qs_loc_env new environment for the localization calculations
218!> \param localized_wfn_control variables and directives for the localization
219!> \param qs_env the qs_env in which the qs_env lives
220!> \param myspin ...
221!> \param do_localize ...
222!> \param loc_coeff ...
223!> \param mo_loc_history ...
224!> \par History
225!> 04.2005 created [MI]
226!> \author MI
227!> \note
228!> similar to the old one, but not quite
229! **************************************************************************************************
230 SUBROUTINE qs_loc_env_init(qs_loc_env, localized_wfn_control, qs_env, myspin, do_localize, &
231 loc_coeff, mo_loc_history)
232
233 TYPE(qs_loc_env_type), POINTER :: qs_loc_env
234 TYPE(localized_wfn_control_type), POINTER :: localized_wfn_control
235 TYPE(qs_environment_type), POINTER :: qs_env
236 INTEGER, INTENT(IN), OPTIONAL :: myspin
237 LOGICAL, INTENT(IN), OPTIONAL :: do_localize
238 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN), &
239 OPTIONAL :: loc_coeff
240 TYPE(cp_fm_type), DIMENSION(:), OPTIONAL, POINTER :: mo_loc_history
241
242 CHARACTER(len=*), PARAMETER :: routinen = 'qs_loc_env_init'
243
244 INTEGER :: dim_op, handle, i, iatom, imo, imoloc, &
245 ispin, j, l_spin, lb, nao, naosub, &
246 natoms, nmo, nmosub, nspins, s_spin, ub
247 REAL(kind=dp) :: my_occ, occ_imo
248 REAL(kind=dp), DIMENSION(:), POINTER :: occupations
249 REAL(kind=dp), DIMENSION(:, :), POINTER :: vecbuffer
250 TYPE(cell_type), POINTER :: cell
251 TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct
252 TYPE(cp_fm_type), DIMENSION(:), POINTER :: moloc_coeff
253 TYPE(cp_fm_type), POINTER :: mat_ptr, mo_coeff
254 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
255 TYPE(distribution_1d_type), POINTER :: local_molecules
256 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
257 TYPE(mp_para_env_type), POINTER :: para_env
258 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
259
260 CALL timeset(routinen, handle)
261
262 NULLIFY (mos, matrix_s, moloc_coeff, particle_set, para_env, cell, local_molecules, occupations, mat_ptr)
263 IF (PRESENT(do_localize)) qs_loc_env%do_localize = do_localize
264 IF (qs_loc_env%do_localize) THEN
265 CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s, cell=cell, &
266 local_molecules=local_molecules, particle_set=particle_set, &
267 para_env=para_env, mos=mos)
268 nspins = SIZE(mos, 1)
269 s_spin = 1
270 l_spin = nspins
271 IF (PRESENT(myspin)) THEN
272 s_spin = myspin
273 l_spin = myspin
274 END IF
275 ALLOCATE (moloc_coeff(s_spin:l_spin))
276 DO ispin = s_spin, l_spin
277 NULLIFY (tmp_fm_struct, mo_coeff)
278 CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, nao=nao, nmo=nmo)
279 nmosub = localized_wfn_control%nloc_states(ispin)
280 CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, &
281 ncol_global=nmosub, para_env=para_env, context=mo_coeff%matrix_struct%context)
282 CALL cp_fm_create(moloc_coeff(ispin), tmp_fm_struct)
283
284 CALL cp_fm_get_info(moloc_coeff(ispin), nrow_global=naosub, &
285 ncol_global=nmosub)
286 cpassert(nao == naosub)
287 IF ((localized_wfn_control%do_homo) .OR. &
288 (localized_wfn_control%set_of_states == state_loc_mixed)) THEN
289 cpassert(nmo >= nmosub)
290 ELSE
291 cpassert(nao - nmo >= nmosub)
292 END IF
293 CALL cp_fm_set_all(moloc_coeff(ispin), 0.0_dp)
294 CALL cp_fm_struct_release(tmp_fm_struct)
295 END DO ! ispin
296 ! Copy the submatrix
297
298 IF (PRESENT(loc_coeff)) ALLOCATE (mat_ptr)
299
300 DO ispin = s_spin, l_spin
301 CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, &
302 occupation_numbers=occupations, nao=nao, nmo=nmo)
303 lb = localized_wfn_control%lu_bound_states(1, ispin)
304 ub = localized_wfn_control%lu_bound_states(2, ispin)
305
306 IF (PRESENT(loc_coeff)) THEN
307 mat_ptr = loc_coeff(ispin)
308 ELSE
309 mat_ptr => mo_coeff
310 END IF
311 IF ((localized_wfn_control%set_of_states == state_loc_list) .OR. &
312 (localized_wfn_control%set_of_states == state_loc_mixed)) THEN
313 ALLOCATE (vecbuffer(1, nao))
314 IF (localized_wfn_control%do_homo) THEN
315 my_occ = occupations(localized_wfn_control%loc_states(1, ispin))
316 END IF
317 nmosub = SIZE(localized_wfn_control%loc_states, 1)
318 cpassert(nmosub > 0)
319 imoloc = 0
320 DO i = lb, ub
321 ! Get the index in the subset
322 imoloc = imoloc + 1
323 ! Get the index in the full set
324 imo = localized_wfn_control%loc_states(i, ispin)
325 IF (localized_wfn_control%do_homo) THEN
326 occ_imo = occupations(imo)
327 IF (abs(occ_imo - my_occ) > localized_wfn_control%eps_occ) THEN
328 IF (localized_wfn_control%localization_method /= do_loc_none) &
329 CALL cp_abort(__location__, &
330 "States with different occupations "// &
331 "cannot be rotated together")
332 END IF
333 END IF
334 ! Take the imo vector from the full set and copy in the imoloc vector of the subset
335 CALL cp_fm_get_submatrix(mat_ptr, vecbuffer, 1, imo, &
336 nao, 1, transpose=.true.)
337 CALL cp_fm_set_submatrix(moloc_coeff(ispin), vecbuffer, 1, imoloc, &
338 nao, 1, transpose=.true.)
339 END DO
340 DEALLOCATE (vecbuffer)
341 ELSE
342 my_occ = occupations(lb)
343 occ_imo = occupations(ub)
344 IF (abs(occ_imo - my_occ) > localized_wfn_control%eps_occ) THEN
345 IF (localized_wfn_control%localization_method /= do_loc_none) &
346 CALL cp_abort(__location__, &
347 "States with different occupations "// &
348 "cannot be rotated together")
349 END IF
350 nmosub = localized_wfn_control%nloc_states(ispin)
351 CALL cp_fm_to_fm(mat_ptr, moloc_coeff(ispin), nmosub, lb, 1)
352 END IF
353
354 ! we have the mo's to be localized now, see if we can rotate them according to the history
355 ! only do that if we have a history of course. The history is filled
356 IF (PRESENT(mo_loc_history)) THEN
357 IF (localized_wfn_control%use_history .AND. ASSOCIATED(mo_loc_history)) THEN
358 CALL rotate_state_to_ref(moloc_coeff(ispin), &
359 mo_loc_history(ispin), matrix_s(1)%matrix)
360 END IF
361 END IF
362
363 END DO
364
365 IF (PRESENT(loc_coeff)) DEALLOCATE (mat_ptr)
366
367 CALL set_qs_loc_env(qs_loc_env=qs_loc_env, cell=cell, local_molecules=local_molecules, &
368 moloc_coeff=moloc_coeff, particle_set=particle_set, para_env=para_env, &
369 localized_wfn_control=localized_wfn_control)
370
371 ! Prepare the operators
372 NULLIFY (tmp_fm_struct, mo_coeff)
373 nmosub = maxval(localized_wfn_control%nloc_states)
374 CALL get_mo_set(mos(1), mo_coeff=mo_coeff)
375 CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nmosub, &
376 ncol_global=nmosub, para_env=para_env, context=mo_coeff%matrix_struct%context)
377
378 IF (localized_wfn_control%operator_type == op_loc_berry) THEN
379 IF (qs_loc_env%cell%orthorhombic) THEN
380 dim_op = 3
381 ELSE
382 dim_op = 6
383 END IF
384 CALL set_qs_loc_env(qs_loc_env=qs_loc_env, dim_op=dim_op)
385 ALLOCATE (qs_loc_env%op_sm_set(2, dim_op))
386 DO i = 1, dim_op
387 DO j = 1, SIZE(qs_loc_env%op_sm_set, 1)
388 NULLIFY (qs_loc_env%op_sm_set(j, i)%matrix)
389 ALLOCATE (qs_loc_env%op_sm_set(j, i)%matrix)
390 CALL dbcsr_copy(qs_loc_env%op_sm_set(j, i)%matrix, matrix_s(1)%matrix, &
391 name="qs_loc_env%op_sm_"//trim(adjustl(cp_to_string(j)))//"-"//trim(adjustl(cp_to_string(i))))
392 CALL dbcsr_set(qs_loc_env%op_sm_set(j, i)%matrix, 0.0_dp)
393 END DO
394 END DO
395
396 ELSEIF (localized_wfn_control%operator_type == op_loc_pipek) THEN
397 natoms = SIZE(qs_loc_env%particle_set, 1)
398 ALLOCATE (qs_loc_env%op_fm_set(natoms, 1))
399 CALL set_qs_loc_env(qs_loc_env=qs_loc_env, dim_op=natoms)
400 DO ispin = 1, SIZE(qs_loc_env%op_fm_set, 2)
401 CALL get_mo_set(mos(ispin), nmo=nmo)
402 DO iatom = 1, natoms
403 CALL cp_fm_create(qs_loc_env%op_fm_set(iatom, ispin), tmp_fm_struct)
404
405 CALL cp_fm_get_info(qs_loc_env%op_fm_set(iatom, ispin), nrow_global=nmosub)
406 cpassert(nmo >= nmosub)
407 CALL cp_fm_set_all(qs_loc_env%op_fm_set(iatom, ispin), 0.0_dp)
408 END DO ! iatom
409 END DO ! ispin
410 ELSE
411 cpabort("Type of operator not implemented")
412 END IF
413 CALL cp_fm_struct_release(tmp_fm_struct)
414
415 IF (localized_wfn_control%operator_type == op_loc_berry) THEN
416
417 CALL initialize_weights(qs_loc_env%cell, qs_loc_env%weights)
418
419 CALL get_berry_operator(qs_loc_env, qs_env)
420
421 ELSEIF (localized_wfn_control%operator_type == op_loc_pipek) THEN
422
423 !! here we don't have to do anything
424 !! CALL get_pipek_mezey_operator ( qs_loc_env, qs_env )
425
426 END IF
427
428 qs_loc_env%molecular_states = .false.
429 qs_loc_env%wannier_states = .false.
430 END IF
431 CALL timestop(handle)
432
433 END SUBROUTINE qs_loc_env_init
434
435! **************************************************************************************************
436!> \brief A wrapper to compute the Berry operator for periodic systems
437!> \param qs_loc_env new environment for the localization calculations
438!> \param qs_env the qs_env in which the qs_env lives
439!> \par History
440!> 04.2005 created [MI]
441!> 04.2018 modified [RZK, ZL]
442!> \author MI
443! **************************************************************************************************
444 SUBROUTINE get_berry_operator(qs_loc_env, qs_env)
445 TYPE(qs_loc_env_type), POINTER :: qs_loc_env
446 TYPE(qs_environment_type), POINTER :: qs_env
447
448 CHARACTER(len=*), PARAMETER :: routinen = 'get_berry_operator'
449
450 INTEGER :: dim_op, handle
451 TYPE(cell_type), POINTER :: cell
452 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: op_sm_set
453
454 CALL timeset(routinen, handle)
455
456 NULLIFY (cell, op_sm_set)
457 CALL get_qs_loc_env(qs_loc_env=qs_loc_env, op_sm_set=op_sm_set, &
458 cell=cell, dim_op=dim_op)
459 CALL compute_berry_operator(qs_env, cell, op_sm_set, dim_op)
460
461 CALL timestop(handle)
462 END SUBROUTINE get_berry_operator
463
464! **************************************************************************************************
465!> \brief Computes the Berry operator for periodic systems
466!> used to define the spread of the MOS
467!> Here the matrix elements of the type <mu|cos(kr)|nu> and <mu|sin(kr)|nu>
468!> are computed, where mu and nu are the contracted basis functions.
469!> Namely the Berry operator is exp(ikr)
470!> k is defined somewhere
471!> the pair lists are exploited and sparse matrixes are constructed
472!> \param qs_env the qs_env in which the qs_env lives
473!> \param cell ...
474!> \param op_sm_set ...
475!> \param dim_op ...
476!> \par History
477!> 04.2005 created [MI]
478!> 04.2018 wrapped old code [RZK, ZL]
479!> \author MI
480!> \note
481!> The intgrals are computed analytically using the primitives GTO
482!> The contraction is performed block-wise
483! **************************************************************************************************
484 SUBROUTINE compute_berry_operator(qs_env, cell, op_sm_set, dim_op)
485 TYPE(qs_environment_type), POINTER :: qs_env
486 TYPE(cell_type), POINTER :: cell
487 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: op_sm_set
488 INTEGER :: dim_op
489
490 CHARACTER(len=*), PARAMETER :: routinen = 'compute_berry_operator'
491
492 INTEGER :: handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, last_jatom, &
493 ldab, ldsa, ldsb, ldwork, maxl, ncoa, ncob, nkind, nrow, nseta, nsetb, sgfa, sgfb
494 INTEGER, DIMENSION(3) :: perd0
495 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
496 npgfb, nsgfa, nsgfb
497 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
498 LOGICAL :: found, new_atom_b
499 REAL(kind=dp) :: dab, kvec(3), rab2, vector_k(3, 6)
500 REAL(kind=dp), DIMENSION(3) :: ra, rab, rb
501 REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
502 REAL(kind=dp), DIMENSION(:, :), POINTER :: cosab, rpgfa, rpgfb, sinab, sphi_a, &
503 sphi_b, work, zeta, zetb
504 TYPE(block_p_type), DIMENSION(:), POINTER :: op_cos, op_sin
505 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
506 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
508 DIMENSION(:), POINTER :: nl_iterator
509 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
510 POINTER :: sab_orb
511 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
512 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
513 TYPE(qs_kind_type), POINTER :: qs_kind
514
515 CALL timeset(routinen, handle)
516 NULLIFY (qs_kind, qs_kind_set)
517 NULLIFY (particle_set)
518 NULLIFY (sab_orb)
519 NULLIFY (cosab, sinab, work)
520 NULLIFY (la_max, la_min, lb_max, lb_min, npgfa, npgfb, nsgfa, nsgfb)
521 NULLIFY (set_radius_a, set_radius_b, rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb)
522
523 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, &
524 particle_set=particle_set, sab_orb=sab_orb)
525
526 nkind = SIZE(qs_kind_set)
527
528 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
529 maxco=ldwork, maxlgto=maxl)
530 ldab = ldwork
531 ALLOCATE (cosab(ldab, ldab))
532 cosab = 0.0_dp
533 ALLOCATE (sinab(ldab, ldab))
534 sinab = 0.0_dp
535 ALLOCATE (work(ldwork, ldwork))
536 work = 0.0_dp
537
538 ALLOCATE (op_cos(dim_op))
539 ALLOCATE (op_sin(dim_op))
540 DO i = 1, dim_op
541 NULLIFY (op_cos(i)%block)
542 NULLIFY (op_sin(i)%block)
543 END DO
544
545 kvec = 0.0_dp
546 vector_k = 0.0_dp
547 vector_k(:, 1) = twopi*cell%h_inv(1, :)
548 vector_k(:, 2) = twopi*cell%h_inv(2, :)
549 vector_k(:, 3) = twopi*cell%h_inv(3, :)
550 vector_k(:, 4) = twopi*(cell%h_inv(1, :) + cell%h_inv(2, :))
551 vector_k(:, 5) = twopi*(cell%h_inv(1, :) + cell%h_inv(3, :))
552 vector_k(:, 6) = twopi*(cell%h_inv(2, :) + cell%h_inv(3, :))
553
554 ! This operator can be used only for periodic systems
555 ! If an isolated system is used the periodicity is overimposed
556 perd0(1:3) = cell%perd(1:3)
557 cell%perd(1:3) = 1
558
559 ALLOCATE (basis_set_list(nkind))
560 DO ikind = 1, nkind
561 qs_kind => qs_kind_set(ikind)
562 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a)
563 IF (ASSOCIATED(basis_set_a)) THEN
564 basis_set_list(ikind)%gto_basis_set => basis_set_a
565 ELSE
566 NULLIFY (basis_set_list(ikind)%gto_basis_set)
567 END IF
568 END DO
569 CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
570 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
571 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
572 iatom=iatom, jatom=jatom, r=rab)
573 basis_set_a => basis_set_list(ikind)%gto_basis_set
574 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
575 basis_set_b => basis_set_list(jkind)%gto_basis_set
576 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
577 ra = pbc(particle_set(iatom)%r, cell)
578 ! basis ikind
579 first_sgfa => basis_set_a%first_sgf
580 la_max => basis_set_a%lmax
581 la_min => basis_set_a%lmin
582 npgfa => basis_set_a%npgf
583 nseta = basis_set_a%nset
584 nsgfa => basis_set_a%nsgf_set
585 rpgfa => basis_set_a%pgf_radius
586 set_radius_a => basis_set_a%set_radius
587 sphi_a => basis_set_a%sphi
588 zeta => basis_set_a%zet
589 ! basis jkind
590 first_sgfb => basis_set_b%first_sgf
591 lb_max => basis_set_b%lmax
592 lb_min => basis_set_b%lmin
593 npgfb => basis_set_b%npgf
594 nsetb = basis_set_b%nset
595 nsgfb => basis_set_b%nsgf_set
596 rpgfb => basis_set_b%pgf_radius
597 set_radius_b => basis_set_b%set_radius
598 sphi_b => basis_set_b%sphi
599 zetb => basis_set_b%zet
600
601 ldsa = SIZE(sphi_a, 1)
602 ldsb = SIZE(sphi_b, 1)
603 IF (inode == 1) last_jatom = 0
604
605 rb = rab + ra
606
607 IF (jatom /= last_jatom) THEN
608 new_atom_b = .true.
609 last_jatom = jatom
610 ELSE
611 new_atom_b = .false.
612 END IF
613
614 IF (new_atom_b) THEN
615 IF (iatom <= jatom) THEN
616 irow = iatom
617 icol = jatom
618 ELSE
619 irow = jatom
620 icol = iatom
621 END IF
622
623 DO i = 1, dim_op
624 NULLIFY (op_cos(i)%block)
625 CALL dbcsr_get_block_p(matrix=op_sm_set(1, i)%matrix, &
626 row=irow, col=icol, block=op_cos(i)%block, found=found)
627 NULLIFY (op_sin(i)%block)
628 CALL dbcsr_get_block_p(matrix=op_sm_set(2, i)%matrix, &
629 row=irow, col=icol, block=op_sin(i)%block, found=found)
630 END DO
631 END IF ! new_atom_b
632
633 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
634 dab = sqrt(rab2)
635
636 nrow = 0
637 DO iset = 1, nseta
638
639 ncoa = npgfa(iset)*ncoset(la_max(iset))
640 sgfa = first_sgfa(1, iset)
641
642 DO jset = 1, nsetb
643
644 ncob = npgfb(jset)*ncoset(lb_max(jset))
645 sgfb = first_sgfb(1, jset)
646
647 IF (set_radius_a(iset) + set_radius_b(jset) >= dab) THEN
648
649! *** Calculate the primitive overlap integrals ***
650 DO i = 1, dim_op
651 kvec(1:3) = vector_k(1:3, i)
652 cosab = 0.0_dp
653 sinab = 0.0_dp
654 CALL cossin(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), &
655 la_min(iset), lb_max(jset), npgfb(jset), zetb(:, jset), &
656 rpgfb(:, jset), lb_min(jset), &
657 ra, rb, kvec, cosab, sinab)
658 CALL contract_cossin(op_cos(i)%block, op_sin(i)%block, &
659 iatom, ncoa, nsgfa(iset), sgfa, sphi_a, ldsa, &
660 jatom, ncob, nsgfb(jset), sgfb, sphi_b, ldsb, &
661 cosab, sinab, ldab, work, ldwork)
662 END DO
663
664 END IF ! >= dab
665
666 END DO ! jset
667
668 nrow = nrow + ncoa
669
670 END DO ! iset
671
672 END DO
673 CALL neighbor_list_iterator_release(nl_iterator)
674
675 ! Set back the correct periodicity
676 cell%perd(1:3) = perd0(1:3)
677
678 DO i = 1, dim_op
679 NULLIFY (op_cos(i)%block)
680 NULLIFY (op_sin(i)%block)
681 END DO
682 DEALLOCATE (op_cos, op_sin)
683
684 DEALLOCATE (cosab, sinab, work, basis_set_list)
685
686 CALL timestop(handle)
687 END SUBROUTINE compute_berry_operator
688
689! **************************************************************************************************
690!> \brief ...
691!> \param qs_loc_env ...
692!> \param section ...
693!> \param mo_array ...
694!> \param coeff_localized ...
695!> \param do_homo ...
696!> \param evals ...
697!> \param do_mixed ...
698! **************************************************************************************************
699 SUBROUTINE loc_write_restart(qs_loc_env, section, mo_array, coeff_localized, &
700 do_homo, evals, do_mixed)
701 TYPE(qs_loc_env_type), POINTER :: qs_loc_env
702 TYPE(section_vals_type), POINTER :: section
703 TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array
704 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: coeff_localized
705 LOGICAL, INTENT(IN) :: do_homo
706 TYPE(cp_1d_r_p_type), DIMENSION(:), OPTIONAL, &
707 POINTER :: evals
708 LOGICAL, INTENT(IN), OPTIONAL :: do_mixed
709
710 CHARACTER(LEN=*), PARAMETER :: routinen = 'loc_write_restart'
711
712 CHARACTER(LEN=default_path_length) :: filename
713 CHARACTER(LEN=default_string_length) :: my_middle
714 INTEGER :: handle, ispin, max_block, nao, nloc, &
715 nmo, output_unit, rst_unit
716 LOGICAL :: my_do_mixed
717 TYPE(cp_logger_type), POINTER :: logger
718 TYPE(section_vals_type), POINTER :: print_key
719
720 CALL timeset(routinen, handle)
721 NULLIFY (logger)
722 logger => cp_get_default_logger()
723 output_unit = cp_logger_get_default_io_unit(logger)
724
725 IF (qs_loc_env%do_localize) THEN
726
727 print_key => section_vals_get_subs_vals(section, "LOC_RESTART")
728 IF (btest(cp_print_key_should_output(logger%iter_info, &
729 section, "LOC_RESTART"), &
730 cp_p_file)) THEN
731
732 ! Open file
733 rst_unit = -1
734
735 my_do_mixed = .false.
736 IF (PRESENT(do_mixed)) my_do_mixed = do_mixed
737 IF (do_homo) THEN
738 my_middle = "LOC_HOMO"
739 ELSEIF (my_do_mixed) THEN
740 my_middle = "LOC_MIXED"
741 ELSE
742 my_middle = "LOC_LUMO"
743 END IF
744
745 rst_unit = cp_print_key_unit_nr(logger, section, "LOC_RESTART", &
746 extension=".wfn", file_status="REPLACE", file_action="WRITE", &
747 file_form="UNFORMATTED", middle_name=trim(my_middle))
748
749 IF (rst_unit > 0) filename = cp_print_key_generate_filename(logger, print_key, &
750 middle_name=trim(my_middle), extension=".wfn", &
751 my_local=.false.)
752
753 IF (output_unit > 0) THEN
754 WRITE (unit=output_unit, fmt="(/,T2,A, A/)") &
755 "LOCALIZATION| Write restart file for the localized MOS : ", &
756 trim(filename)
757 END IF
758
759 IF (rst_unit > 0) THEN
760 WRITE (rst_unit) qs_loc_env%localized_wfn_control%set_of_states
761 WRITE (rst_unit) qs_loc_env%localized_wfn_control%lu_bound_states
762 WRITE (rst_unit) qs_loc_env%localized_wfn_control%nloc_states
763 END IF
764
765 DO ispin = 1, SIZE(coeff_localized)
766 associate(mo_coeff => coeff_localized(ispin))
767 CALL cp_fm_get_info(mo_coeff, nrow_global=nao, ncol_global=nmo, ncol_block=max_block)
768 nloc = qs_loc_env%localized_wfn_control%nloc_states(ispin)
769 IF (rst_unit > 0) THEN
770 WRITE (rst_unit) qs_loc_env%localized_wfn_control%loc_states(1:nloc, ispin)
771 IF (do_homo .OR. my_do_mixed) THEN
772 WRITE (rst_unit) nmo, &
773 mo_array(ispin)%homo, &
774 mo_array(ispin)%lfomo, &
775 mo_array(ispin)%nelectron
776 WRITE (rst_unit) mo_array(ispin)%eigenvalues(1:nmo), &
777 mo_array(ispin)%occupation_numbers(1:nmo)
778 ELSE
779 WRITE (rst_unit) nmo
780 WRITE (rst_unit) evals(ispin)%array(1:nmo)
781 END IF
782 END IF
783
784 CALL cp_fm_write_unformatted(mo_coeff, rst_unit)
785 END associate
786
787 END DO
788
789 CALL cp_print_key_finished_output(rst_unit, logger, section, &
790 "LOC_RESTART")
791 END IF
792
793 END IF
794
795 CALL timestop(handle)
796
797 END SUBROUTINE loc_write_restart
798
799! **************************************************************************************************
800!> \brief ...
801!> \param qs_loc_env ...
802!> \param mos ...
803!> \param mos_localized ...
804!> \param section ...
805!> \param section2 ...
806!> \param para_env ...
807!> \param do_homo ...
808!> \param restart_found ...
809!> \param evals ...
810!> \param do_mixed ...
811! **************************************************************************************************
812 SUBROUTINE loc_read_restart(qs_loc_env, mos, mos_localized, section, section2, para_env, &
813 do_homo, restart_found, evals, do_mixed)
814
815 TYPE(qs_loc_env_type), POINTER :: qs_loc_env
816 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
817 TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT) :: mos_localized
818 TYPE(section_vals_type), POINTER :: section, section2
819 TYPE(mp_para_env_type), POINTER :: para_env
820 LOGICAL, INTENT(IN) :: do_homo
821 LOGICAL, INTENT(INOUT) :: restart_found
822 TYPE(cp_1d_r_p_type), DIMENSION(:), OPTIONAL, &
823 POINTER :: evals
824 LOGICAL, INTENT(IN), OPTIONAL :: do_mixed
825
826 CHARACTER(len=*), PARAMETER :: routinen = 'loc_read_restart'
827
828 CHARACTER(LEN=25) :: fname_key
829 CHARACTER(LEN=default_path_length) :: filename
830 CHARACTER(LEN=default_string_length) :: my_middle
831 INTEGER :: handle, homo_read, i, ispin, lfomo_read, max_nloc, n_rep_val, nao, &
832 nelectron_read, nloc, nmo, nmo_read, nspin, output_unit, rst_unit
833 LOGICAL :: file_exists, my_do_mixed
834 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eig_read, occ_read
835 REAL(kind=dp), DIMENSION(:, :), POINTER :: vecbuffer
836 TYPE(cp_logger_type), POINTER :: logger
837 TYPE(section_vals_type), POINTER :: print_key
838
839 CALL timeset(routinen, handle)
840
841 logger => cp_get_default_logger()
842
843 nspin = SIZE(mos_localized)
844 nao = mos(1)%nao
845 rst_unit = -1
846
847 output_unit = cp_print_key_unit_nr(logger, section2, &
848 "PROGRAM_RUN_INFO", extension=".Log")
849
850 my_do_mixed = .false.
851 IF (PRESENT(do_mixed)) my_do_mixed = do_mixed
852 IF (do_homo) THEN
853 fname_key = "LOCHOMO_RESTART_FILE_NAME"
854 ELSEIF (my_do_mixed) THEN
855 fname_key = "LOCMIXD_RESTART_FILE_NAME"
856 ELSE
857 fname_key = "LOCLUMO_RESTART_FILE_NAME"
858 IF (.NOT. PRESENT(evals)) &
859 cpabort("Missing argument to localize unoccupied states.")
860 END IF
861
862 file_exists = .false.
863 CALL section_vals_val_get(section, fname_key, n_rep_val=n_rep_val)
864 IF (n_rep_val > 0) THEN
865 CALL section_vals_val_get(section, fname_key, c_val=filename)
866 ELSE
867
868 print_key => section_vals_get_subs_vals(section2, "LOC_RESTART")
869 IF (do_homo) THEN
870 my_middle = "LOC_HOMO"
871 ELSEIF (my_do_mixed) THEN
872 my_middle = "LOC_MIXED"
873 ELSE
874 my_middle = "LOC_LUMO"
875 END IF
876 filename = cp_print_key_generate_filename(logger, print_key, &
877 middle_name=trim(my_middle), extension=".wfn", &
878 my_local=.false.)
879 END IF
880
881 IF (para_env%is_source()) INQUIRE (file=filename, exist=file_exists)
882
883 IF (file_exists) THEN
884 IF (para_env%is_source()) THEN
885 CALL open_file(file_name=filename, &
886 file_action="READ", &
887 file_form="UNFORMATTED", &
888 file_status="OLD", &
889 unit_number=rst_unit)
890
891 READ (rst_unit) qs_loc_env%localized_wfn_control%set_of_states
892 READ (rst_unit) qs_loc_env%localized_wfn_control%lu_bound_states
893 READ (rst_unit) qs_loc_env%localized_wfn_control%nloc_states
894 END IF
895 ELSE
896 IF (output_unit > 0) WRITE (output_unit, "(/,T10,A)") &
897 "Restart file not available filename=<"//trim(filename)//'>'
898 END IF
899 CALL para_env%bcast(file_exists)
900
901 IF (file_exists) THEN
902 restart_found = .true.
903
904 CALL para_env%bcast(qs_loc_env%localized_wfn_control%set_of_states)
905 CALL para_env%bcast(qs_loc_env%localized_wfn_control%lu_bound_states)
906 CALL para_env%bcast(qs_loc_env%localized_wfn_control%nloc_states)
907
908 max_nloc = maxval(qs_loc_env%localized_wfn_control%nloc_states(:))
909
910 ALLOCATE (vecbuffer(1, nao))
911 IF (ASSOCIATED(qs_loc_env%localized_wfn_control%loc_states)) THEN
912 DEALLOCATE (qs_loc_env%localized_wfn_control%loc_states)
913 END IF
914 ALLOCATE (qs_loc_env%localized_wfn_control%loc_states(max_nloc, 2))
915 qs_loc_env%localized_wfn_control%loc_states = 0
916
917 DO ispin = 1, nspin
918 IF (do_homo .OR. do_mixed) THEN
919 nmo = mos(ispin)%nmo
920 ELSE
921 nmo = SIZE(evals(ispin)%array, 1)
922 END IF
923 IF (para_env%is_source() .AND. (nmo > 0)) THEN
924 nloc = qs_loc_env%localized_wfn_control%nloc_states(ispin)
925 READ (rst_unit) qs_loc_env%localized_wfn_control%loc_states(1:nloc, ispin)
926 IF (do_homo .OR. do_mixed) THEN
927 READ (rst_unit) nmo_read, homo_read, lfomo_read, nelectron_read
928 ALLOCATE (eig_read(nmo_read), occ_read(nmo_read))
929 eig_read = 0.0_dp
930 occ_read = 0.0_dp
931 READ (rst_unit) eig_read(1:nmo_read), occ_read(1:nmo_read)
932 ELSE
933 READ (rst_unit) nmo_read
934 ALLOCATE (eig_read(nmo_read))
935 eig_read = 0.0_dp
936 READ (rst_unit) eig_read(1:nmo_read)
937 END IF
938 IF (nmo_read < nmo) &
939 CALL cp_warn(__location__, &
940 "The number of MOs on the restart unit is smaller than the number of "// &
941 "the allocated MOs. ")
942 IF (nmo_read > nmo) &
943 CALL cp_warn(__location__, &
944 "The number of MOs on the restart unit is greater than the number of "// &
945 "the allocated MOs. The read MO set will be truncated!")
946
947 nmo = min(nmo, nmo_read)
948 IF (do_homo .OR. do_mixed) THEN
949 mos(ispin)%eigenvalues(1:nmo) = eig_read(1:nmo)
950 mos(ispin)%occupation_numbers(1:nmo) = occ_read(1:nmo)
951 DEALLOCATE (eig_read, occ_read)
952 ELSE
953 evals(ispin)%array(1:nmo) = eig_read(1:nmo)
954 DEALLOCATE (eig_read)
955 END IF
956
957 END IF
958 IF (do_homo .OR. do_mixed) THEN
959 CALL para_env%bcast(mos(ispin)%eigenvalues)
960 CALL para_env%bcast(mos(ispin)%occupation_numbers)
961 ELSE
962 CALL para_env%bcast(evals(ispin)%array)
963 END IF
964
965 DO i = 1, nmo
966 IF (para_env%is_source()) THEN
967 READ (rst_unit) vecbuffer
968 ELSE
969 vecbuffer(1, :) = 0.0_dp
970 END IF
971 CALL para_env%bcast(vecbuffer)
972 CALL cp_fm_set_submatrix(mos_localized(ispin), &
973 vecbuffer, 1, i, nao, 1, transpose=.true.)
974 END DO
975 END DO
976
977 CALL para_env%bcast(qs_loc_env%localized_wfn_control%loc_states)
978
979 DEALLOCATE (vecbuffer)
980
981 END IF
982
983 ! Close restart file
984 IF (para_env%is_source()) THEN
985 IF (file_exists) CALL close_file(unit_number=rst_unit)
986 END IF
987
988 CALL timestop(handle)
989
990 END SUBROUTINE loc_read_restart
991
992! **************************************************************************************************
993!> \brief initializes everything needed for localization of the HOMOs
994!> \param qs_loc_env ...
995!> \param loc_section ...
996!> \param do_homo ...
997!> \param do_mixed ...
998!> \param do_xas ...
999!> \param nloc_xas ...
1000!> \param spin_xas ...
1001!> \par History
1002!> 2009 created
1003! **************************************************************************************************
1004 SUBROUTINE qs_loc_control_init(qs_loc_env, loc_section, do_homo, do_mixed, &
1005 do_xas, nloc_xas, spin_xas)
1006
1007 TYPE(qs_loc_env_type), POINTER :: qs_loc_env
1008 TYPE(section_vals_type), POINTER :: loc_section
1009 LOGICAL, INTENT(IN) :: do_homo
1010 LOGICAL, INTENT(IN), OPTIONAL :: do_mixed, do_xas
1011 INTEGER, INTENT(IN), OPTIONAL :: nloc_xas, spin_xas
1012
1013 LOGICAL :: my_do_mixed
1014 TYPE(localized_wfn_control_type), POINTER :: localized_wfn_control
1015
1016 NULLIFY (localized_wfn_control)
1017
1018 IF (PRESENT(do_mixed)) THEN
1019 my_do_mixed = do_mixed
1020 ELSE
1021 my_do_mixed = .false.
1022 END IF
1023 CALL localized_wfn_control_create(localized_wfn_control)
1024 CALL set_qs_loc_env(qs_loc_env, localized_wfn_control=localized_wfn_control)
1025 CALL localized_wfn_control_release(localized_wfn_control)
1026 CALL get_qs_loc_env(qs_loc_env, localized_wfn_control=localized_wfn_control)
1027 localized_wfn_control%do_homo = do_homo
1028 localized_wfn_control%do_mixed = my_do_mixed
1029 CALL read_loc_section(localized_wfn_control, loc_section, qs_loc_env%do_localize, &
1030 my_do_mixed, do_xas, nloc_xas, spin_xas)
1031
1032 END SUBROUTINE qs_loc_control_init
1033
1034! **************************************************************************************************
1035!> \brief initializes everything needed for localization of the molecular orbitals
1036!> \param qs_env ...
1037!> \param qs_loc_env ...
1038!> \param localize_section ...
1039!> \param mos_localized ...
1040!> \param do_homo ...
1041!> \param do_mo_cubes ...
1042!> \param mo_loc_history ...
1043!> \param evals ...
1044!> \param tot_zeff_corr ...
1045!> \param do_mixed ...
1046! **************************************************************************************************
1047 SUBROUTINE qs_loc_init(qs_env, qs_loc_env, localize_section, mos_localized, &
1048 do_homo, do_mo_cubes, mo_loc_history, evals, &
1049 tot_zeff_corr, do_mixed)
1050
1051 TYPE(qs_environment_type), POINTER :: qs_env
1052 TYPE(qs_loc_env_type), POINTER :: qs_loc_env
1053 TYPE(section_vals_type), POINTER :: localize_section
1054 TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT) :: mos_localized
1055 LOGICAL, OPTIONAL :: do_homo, do_mo_cubes
1056 TYPE(cp_fm_type), DIMENSION(:), OPTIONAL, POINTER :: mo_loc_history
1057 TYPE(cp_1d_r_p_type), DIMENSION(:), OPTIONAL, &
1058 POINTER :: evals
1059 REAL(kind=dp), INTENT(IN), OPTIONAL :: tot_zeff_corr
1060 LOGICAL, OPTIONAL :: do_mixed
1061
1062 CHARACTER(len=*), PARAMETER :: routinen = 'qs_loc_init'
1063
1064 INTEGER :: handle, homo, i, ilast_intocc, ilow, ispin, iup, n_mo(2), n_mos(2), nao, &
1065 nelectron, nextra, nmoloc(2), nocc, npocc, nspin, output_unit
1066 LOGICAL :: my_do_homo, my_do_mixed, my_do_mo_cubes, &
1067 restart_found
1068 REAL(kind=dp) :: maxocc, my_tot_zeff_corr
1069 REAL(kind=dp), DIMENSION(:), POINTER :: mo_eigenvalues, occupation
1070 TYPE(cp_fm_type), POINTER :: mo_coeff
1071 TYPE(cp_logger_type), POINTER :: logger
1072 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_rmpv, mo_derivs
1073 TYPE(dft_control_type), POINTER :: dft_control
1074 TYPE(localized_wfn_control_type), POINTER :: localized_wfn_control
1075 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1076 TYPE(mp_para_env_type), POINTER :: para_env
1077 TYPE(scf_control_type), POINTER :: scf_control
1078 TYPE(section_vals_type), POINTER :: loc_print_section
1079
1080 CALL timeset(routinen, handle)
1081
1082 NULLIFY (mos, mo_coeff, mo_eigenvalues, occupation, ks_rmpv, mo_derivs, scf_control, para_env)
1083 CALL get_qs_env(qs_env, &
1084 mos=mos, &
1085 matrix_ks=ks_rmpv, &
1086 mo_derivs=mo_derivs, &
1087 scf_control=scf_control, &
1088 dft_control=dft_control, &
1089 para_env=para_env)
1090
1091 loc_print_section => section_vals_get_subs_vals(localize_section, "PRINT")
1092
1093 logger => cp_get_default_logger()
1094 output_unit = cp_logger_get_default_io_unit(logger)
1095
1096 nspin = SIZE(mos)
1097 IF (PRESENT(do_homo)) THEN
1098 my_do_homo = do_homo
1099 ELSE
1100 my_do_homo = .true.
1101 END IF
1102 IF (PRESENT(do_mo_cubes)) THEN
1103 my_do_mo_cubes = do_mo_cubes
1104 ELSE
1105 my_do_mo_cubes = .false.
1106 END IF
1107 IF (PRESENT(do_mixed)) THEN
1108 my_do_mixed = do_mixed
1109 ELSE
1110 my_do_mixed = .false.
1111 END IF
1112 IF (PRESENT(tot_zeff_corr)) THEN
1113 my_tot_zeff_corr = tot_zeff_corr
1114 ELSE
1115 my_tot_zeff_corr = 0.0_dp
1116 END IF
1117 restart_found = .false.
1118
1119 IF (qs_loc_env%do_localize) THEN
1120 ! Some setup for MOs to be localized
1121 CALL get_qs_loc_env(qs_loc_env, localized_wfn_control=localized_wfn_control)
1122 IF (localized_wfn_control%loc_restart) THEN
1123 IF (localized_wfn_control%nextra > 0) THEN
1124 ! currently only the occupied guess is read
1125 my_do_homo = .false.
1126 END IF
1127 CALL loc_read_restart(qs_loc_env, mos, mos_localized, localize_section, &
1128 loc_print_section, para_env, my_do_homo, restart_found, evals=evals, &
1129 do_mixed=my_do_mixed)
1130 IF (output_unit > 0) WRITE (output_unit, "(/,T2,A,A)") "LOCALIZATION| ", &
1131 " The orbitals to be localized are read from localization restart file."
1132 nmoloc = localized_wfn_control%nloc_states
1133 localized_wfn_control%nguess = nmoloc
1134 IF (localized_wfn_control%nextra > 0) THEN
1135 ! reset different variables in localized_wfn_control:
1136 ! lu_bound_states, nloc_states, loc_states
1137 localized_wfn_control%loc_restart = restart_found
1138 localized_wfn_control%set_of_states = state_loc_mixed
1139 DO ispin = 1, nspin
1140 CALL get_mo_set(mos(ispin), homo=homo, occupation_numbers=occupation, &
1141 maxocc=maxocc)
1142 nextra = localized_wfn_control%nextra
1143 nocc = homo
1144 DO i = nocc, 1, -1
1145 IF (maxocc - occupation(i) < localized_wfn_control%eps_occ) THEN
1146 ilast_intocc = i
1147 EXIT
1148 END IF
1149 END DO
1150 nocc = ilast_intocc
1151 npocc = homo - nocc
1152 nmoloc(ispin) = nocc + nextra
1153 localized_wfn_control%lu_bound_states(1, ispin) = 1
1154 localized_wfn_control%lu_bound_states(2, ispin) = nmoloc(ispin)
1155 localized_wfn_control%nloc_states(ispin) = nmoloc(ispin)
1156 END DO
1157 my_do_homo = .false.
1158 END IF
1159 END IF
1160 IF (.NOT. restart_found) THEN
1161 nmoloc = 0
1162 DO ispin = 1, nspin
1163 CALL get_mo_set(mos(ispin), nmo=n_mo(ispin), nelectron=nelectron, homo=homo, nao=nao, &
1164 mo_coeff=mo_coeff, eigenvalues=mo_eigenvalues, occupation_numbers=occupation, &
1165 maxocc=maxocc)
1166 ! Get eigenstates (only needed if not already calculated before)
1167 IF ((.NOT. my_do_mo_cubes) &
1168 ! .OR. section_get_ival(dft_section,"PRINT%MO_CUBES%NHOMO")==0)&
1169 .AND. my_do_homo .AND. qs_env%scf_env%method == ot_method_nr .AND. (.NOT. dft_control%restricted)) THEN
1170 CALL make_mo_eig(mos, nspin, ks_rmpv, scf_control, mo_derivs)
1171 END IF
1172 IF (localized_wfn_control%set_of_states == state_loc_all .AND. my_do_homo) THEN
1173 nmoloc(ispin) = nint(nelectron/occupation(1))
1174 IF (n_mo(ispin) > homo) THEN
1175 DO i = nmoloc(ispin), 1, -1
1176 IF (occupation(1) - occupation(i) < localized_wfn_control%eps_occ) THEN
1177 ilast_intocc = i
1178 EXIT
1179 END IF
1180 END DO
1181 ELSE
1182 ilast_intocc = nmoloc(ispin)
1183 END IF
1184 nmoloc(ispin) = ilast_intocc
1185 localized_wfn_control%lu_bound_states(1, ispin) = 1
1186 localized_wfn_control%lu_bound_states(2, ispin) = ilast_intocc
1187 IF (nmoloc(ispin) /= n_mo(ispin)) THEN
1188 IF (output_unit > 0) &
1189 WRITE (output_unit, "(/,T2,A,I4,A,I6,A,/,T15,A,F12.6,A,F12.6,A)") &
1190 "LOCALIZATION| Spin ", ispin, " The first ", &
1191 ilast_intocc, " occupied orbitals are localized,", " with energies from ", &
1192 mo_eigenvalues(1), " to ", mo_eigenvalues(ilast_intocc), " [a.u.]."
1193 END IF
1194 ELSE IF (localized_wfn_control%set_of_states == energy_loc_range .AND. my_do_homo) THEN
1195 ilow = 0
1196 iup = 0
1197 DO i = 1, n_mo(ispin)
1198 IF (mo_eigenvalues(i) >= localized_wfn_control%lu_ene_bound(1)) THEN
1199 ilow = i
1200 EXIT
1201 END IF
1202 END DO
1203 DO i = n_mo(ispin), 1, -1
1204 IF (mo_eigenvalues(i) <= localized_wfn_control%lu_ene_bound(2)) THEN
1205 iup = i
1206 EXIT
1207 END IF
1208 END DO
1209 localized_wfn_control%lu_bound_states(1, ispin) = ilow
1210 localized_wfn_control%lu_bound_states(2, ispin) = iup
1211 localized_wfn_control%nloc_states(ispin) = iup - ilow + 1
1212 nmoloc(ispin) = localized_wfn_control%nloc_states(ispin)
1213 IF (occupation(ilow) - occupation(iup) > localized_wfn_control%eps_occ) THEN
1214 CALL cp_abort(__location__, &
1215 "The selected energy range includes orbitals with different occupation number. "// &
1216 "The localization procedure cannot be applied.")
1217 END IF
1218 IF (output_unit > 0) WRITE (output_unit, "(/,T2,A,I4,A,I6,A)") "LOCALIZATION| Spin ", ispin, " : ", &
1219 nmoloc(ispin), " orbitals in the selected energy range are localized."
1220 ELSE IF (localized_wfn_control%set_of_states == state_loc_all .AND. (.NOT. my_do_homo)) THEN
1221 nmoloc(ispin) = n_mo(ispin) - homo
1222 localized_wfn_control%lu_bound_states(1, ispin) = homo + 1
1223 localized_wfn_control%lu_bound_states(2, ispin) = n_mo(ispin)
1224 IF (output_unit > 0) &
1225 WRITE (output_unit, "(/,T2,A,I4,A,I6,A,/,T15,A,F12.6,A,F12.6,A)") &
1226 "LOCALIZATION| Spin ", ispin, " The first ", &
1227 nmoloc(ispin), " virtual orbitals are localized,", " with energies from ", &
1228 mo_eigenvalues(homo + 1), " to ", mo_eigenvalues(n_mo(ispin)), " [a.u.]."
1229 ELSE IF (localized_wfn_control%set_of_states == state_loc_mixed) THEN
1230 nextra = localized_wfn_control%nextra
1231 nocc = homo
1232 DO i = nocc, 1, -1
1233 IF (maxocc - occupation(i) < localized_wfn_control%eps_occ) THEN
1234 ilast_intocc = i
1235 EXIT
1236 END IF
1237 END DO
1238 nocc = ilast_intocc
1239 npocc = homo - nocc
1240 nmoloc(ispin) = nocc + nextra
1241 localized_wfn_control%lu_bound_states(1, ispin) = 1
1242 localized_wfn_control%lu_bound_states(2, ispin) = nmoloc(ispin)
1243 IF (output_unit > 0) &
1244 WRITE (output_unit, "(/,T2,A,I4,A,I6,A,/,T15,A,I6,/,T15,A,I6,/,T15,A,I6,/,T15,A,F12.6,A)") &
1245 "LOCALIZATION| Spin ", ispin, " The first ", &
1246 nmoloc(ispin), " orbitals are localized.", &
1247 "Number of fully occupied MOs: ", nocc, &
1248 "Number of partially occupied MOs: ", npocc, &
1249 "Number of extra degrees of freedom: ", nextra, &
1250 "Excess charge: ", my_tot_zeff_corr, " electrons"
1251 ELSE
1252 nmoloc(ispin) = min(localized_wfn_control%nloc_states(1), n_mo(ispin))
1253 IF (output_unit > 0 .AND. my_do_homo) WRITE (output_unit, "(/,T2,A,I4,A,I6,A)") "LOCALIZATION| Spin ", ispin, &
1254 " : ", nmoloc(ispin), " occupied orbitals are localized, as given in the input list."
1255 IF (output_unit > 0 .AND. (.NOT. my_do_homo)) WRITE (output_unit, "(/,T2,A,I4,A,I6,A)") "LOCALIZATION| Spin ", &
1256 ispin, " : ", nmoloc(ispin), " unoccupied orbitals are localized, as given in the input list."
1257 IF (n_mo(ispin) > homo .AND. my_do_homo) THEN
1258 ilow = localized_wfn_control%loc_states(1, ispin)
1259 DO i = 2, nmoloc(ispin)
1260 iup = localized_wfn_control%loc_states(i, ispin)
1261 IF (abs(occupation(ilow) - occupation(iup)) > localized_wfn_control%eps_occ) THEN
1262 ! write warning
1263 CALL cp_warn(__location__, &
1264 "User requested the calculation of localized wavefunction from a subset of MOs, "// &
1265 "including MOs with different occupations. Check the selected subset, "// &
1266 "the electronic density is not invariant with "// &
1267 "respect to rotations among orbitals with different occupation numbers!")
1268 END IF
1269 END DO
1270 END IF
1271 END IF
1272 END DO ! ispin
1273 n_mos(:) = nao - n_mo(:)
1274 IF (my_do_homo .OR. my_do_mixed) n_mos = n_mo
1275 CALL set_loc_wfn_lists(localized_wfn_control, nmoloc, n_mos, nspin)
1276 END IF
1277 CALL set_loc_centers(localized_wfn_control, nmoloc, nspin)
1278 IF (my_do_homo .OR. my_do_mixed) THEN
1279 CALL qs_loc_env_init(qs_loc_env, localized_wfn_control, qs_env, &
1280 loc_coeff=mos_localized, mo_loc_history=mo_loc_history)
1281 END IF
1282 ELSE
1283 ! Let's inform in case the section is not present in the input
1284 CALL cp_warn(__location__, &
1285 "User requested the calculation of the localized wavefunction but the section "// &
1286 "LOCALIZE was not specified. Localization will not be performed!")
1287 END IF
1288
1289 CALL timestop(handle)
1290
1291 END SUBROUTINE qs_loc_init
1292
1293! **************************************************************************************************
1294!> \brief read the controlparameter from input, using the new input scheme
1295!> \param localized_wfn_control ...
1296!> \param loc_section ...
1297!> \param localize ...
1298!> \param do_mixed ...
1299!> \param do_xas ...
1300!> \param nloc_xas ...
1301!> \param spin_channel_xas ...
1302!> \par History
1303!> 05.2005 created [MI]
1304! **************************************************************************************************
1305 SUBROUTINE read_loc_section(localized_wfn_control, loc_section, &
1306 localize, do_mixed, do_xas, nloc_xas, spin_channel_xas)
1307
1308 TYPE(localized_wfn_control_type), POINTER :: localized_wfn_control
1309 TYPE(section_vals_type), POINTER :: loc_section
1310 LOGICAL, INTENT(OUT) :: localize
1311 LOGICAL, INTENT(IN), OPTIONAL :: do_mixed, do_xas
1312 INTEGER, INTENT(IN), OPTIONAL :: nloc_xas, spin_channel_xas
1313
1314 INTEGER :: i, ind, ir, n_list, n_rep, n_state, &
1315 nextra, nline, other_spin, &
1316 output_unit, spin_xas
1317 INTEGER, DIMENSION(:), POINTER :: list, loc_list
1318 LOGICAL :: my_do_mixed, my_do_xas
1319 REAL(dp), POINTER :: ene(:)
1320 TYPE(cp_logger_type), POINTER :: logger
1321 TYPE(section_vals_type), POINTER :: loc_print_section
1322
1323 my_do_xas = .false.
1324 spin_xas = 1
1325 IF (PRESENT(do_xas)) THEN
1326 my_do_xas = do_xas
1327 cpassert(PRESENT(nloc_xas))
1328 END IF
1329 IF (PRESENT(spin_channel_xas)) spin_xas = spin_channel_xas
1330 my_do_mixed = .false.
1331 IF (PRESENT(do_mixed)) THEN
1332 my_do_mixed = do_mixed
1333 END IF
1334 cpassert(ASSOCIATED(loc_section))
1335 NULLIFY (logger)
1336 logger => cp_get_default_logger()
1337
1338 CALL section_vals_val_get(loc_section, "_SECTION_PARAMETERS_", l_val=localize)
1339 IF (localize) THEN
1340 loc_print_section => section_vals_get_subs_vals(loc_section, "PRINT")
1341 NULLIFY (list)
1342 NULLIFY (loc_list)
1343 localized_wfn_control%lu_bound_states = 0
1344 localized_wfn_control%lu_ene_bound = 0.0_dp
1345 localized_wfn_control%nloc_states = 0
1346 localized_wfn_control%set_of_states = 0
1347 localized_wfn_control%nextra = 0
1348 n_state = 0
1349
1350 CALL section_vals_val_get(loc_section, "MAX_ITER", &
1351 i_val=localized_wfn_control%max_iter)
1352 CALL section_vals_val_get(loc_section, "MAX_CRAZY_ANGLE", &
1353 r_val=localized_wfn_control%max_crazy_angle)
1354 CALL section_vals_val_get(loc_section, "CRAZY_SCALE", &
1355 r_val=localized_wfn_control%crazy_scale)
1356 CALL section_vals_val_get(loc_section, "EPS_OCCUPATION", &
1357 r_val=localized_wfn_control%eps_occ)
1358 CALL section_vals_val_get(loc_section, "CRAZY_USE_DIAG", &
1359 l_val=localized_wfn_control%crazy_use_diag)
1360 CALL section_vals_val_get(loc_section, "OUT_ITER_EACH", &
1361 i_val=localized_wfn_control%out_each)
1362 CALL section_vals_val_get(loc_section, "EPS_LOCALIZATION", &
1363 r_val=localized_wfn_control%eps_localization)
1364 CALL section_vals_val_get(loc_section, "MIN_OR_MAX", &
1365 i_val=localized_wfn_control%min_or_max)
1366 CALL section_vals_val_get(loc_section, "JACOBI_FALLBACK", &
1367 l_val=localized_wfn_control%jacobi_fallback)
1368 CALL section_vals_val_get(loc_section, "JACOBI_REFINEMENT", &
1369 l_val=localized_wfn_control%jacobi_refinement)
1370 CALL section_vals_val_get(loc_section, "METHOD", &
1371 i_val=localized_wfn_control%localization_method)
1372 CALL section_vals_val_get(loc_section, "OPERATOR", &
1373 i_val=localized_wfn_control%operator_type)
1374 CALL section_vals_val_get(loc_section, "RESTART", &
1375 l_val=localized_wfn_control%loc_restart)
1376 CALL section_vals_val_get(loc_section, "USE_HISTORY", &
1377 l_val=localized_wfn_control%use_history)
1378 CALL section_vals_val_get(loc_section, "NEXTRA", &
1379 i_val=localized_wfn_control%nextra)
1380 CALL section_vals_val_get(loc_section, "CPO_GUESS", &
1381 i_val=localized_wfn_control%coeff_po_guess)
1382 CALL section_vals_val_get(loc_section, "CPO_GUESS_SPACE", &
1383 i_val=localized_wfn_control%coeff_po_guess_mo_space)
1384 CALL section_vals_val_get(loc_section, "CG_PO", &
1385 l_val=localized_wfn_control%do_cg_po)
1386
1387 IF (localized_wfn_control%do_homo) THEN
1388 ! List of States HOMO
1389 CALL section_vals_val_get(loc_section, "LIST", n_rep_val=n_rep)
1390 IF (n_rep > 0) THEN
1391 n_list = 0
1392 DO ir = 1, n_rep
1393 NULLIFY (list)
1394 CALL section_vals_val_get(loc_section, "LIST", i_rep_val=ir, i_vals=list)
1395 IF (ASSOCIATED(list)) THEN
1396 CALL reallocate(loc_list, 1, n_list + SIZE(list))
1397 DO i = 1, SIZE(list)
1398 loc_list(n_list + i) = list(i)
1399 END DO ! i
1400 n_list = n_list + SIZE(list)
1401 END IF
1402 END DO ! ir
1403 IF (n_list /= 0) THEN
1404 localized_wfn_control%set_of_states = state_loc_list
1405 ALLOCATE (localized_wfn_control%loc_states(n_list, 2))
1406 localized_wfn_control%loc_states = 0
1407 localized_wfn_control%loc_states(:, 1) = loc_list(:)
1408 localized_wfn_control%loc_states(:, 2) = loc_list(:)
1409 localized_wfn_control%nloc_states(1) = n_list
1410 localized_wfn_control%nloc_states(2) = n_list
1411 IF (my_do_xas) THEN
1412 other_spin = 2
1413 IF (spin_xas == 2) other_spin = 1
1414 localized_wfn_control%nloc_states(other_spin) = 0
1415 localized_wfn_control%loc_states(:, other_spin) = 0
1416 END IF
1417 DEALLOCATE (loc_list)
1418 END IF
1419 END IF
1420
1421 ELSE
1422 ! List of States LUMO
1423 CALL section_vals_val_get(loc_section, "LIST_UNOCCUPIED", n_rep_val=n_rep)
1424 IF (n_rep > 0) THEN
1425 n_list = 0
1426 DO ir = 1, n_rep
1427 NULLIFY (list)
1428 CALL section_vals_val_get(loc_section, "LIST_UNOCCUPIED", i_rep_val=ir, i_vals=list)
1429 IF (ASSOCIATED(list)) THEN
1430 CALL reallocate(loc_list, 1, n_list + SIZE(list))
1431 DO i = 1, SIZE(list)
1432 loc_list(n_list + i) = list(i)
1433 END DO ! i
1434 n_list = n_list + SIZE(list)
1435 END IF
1436 END DO ! ir
1437 IF (n_list /= 0) THEN
1438 localized_wfn_control%set_of_states = state_loc_list
1439 ALLOCATE (localized_wfn_control%loc_states(n_list, 2))
1440 localized_wfn_control%loc_states = 0
1441 localized_wfn_control%loc_states(:, 1) = loc_list(:)
1442 localized_wfn_control%loc_states(:, 2) = loc_list(:)
1443 localized_wfn_control%nloc_states(1) = n_list
1444 DEALLOCATE (loc_list)
1445 END IF
1446 END IF
1447 END IF
1448
1449 IF (localized_wfn_control%set_of_states == 0) THEN
1450 CALL section_vals_val_get(loc_section, "ENERGY_RANGE", r_vals=ene)
1451 IF (ene(1) /= ene(2)) THEN
1452 localized_wfn_control%set_of_states = energy_loc_range
1453 localized_wfn_control%lu_ene_bound(1) = ene(1)
1454 localized_wfn_control%lu_ene_bound(2) = ene(2)
1455 END IF
1456 END IF
1457
1458 ! All States or XAS specific states
1459 IF (localized_wfn_control%set_of_states == 0) THEN
1460 IF (my_do_xas) THEN
1461 localized_wfn_control%set_of_states = state_loc_range
1462 localized_wfn_control%nloc_states(:) = 0
1463 localized_wfn_control%lu_bound_states(1, :) = 0
1464 localized_wfn_control%lu_bound_states(2, :) = 0
1465 localized_wfn_control%nloc_states(spin_xas) = nloc_xas
1466 localized_wfn_control%lu_bound_states(1, spin_xas) = 1
1467 localized_wfn_control%lu_bound_states(2, spin_xas) = nloc_xas
1468 ELSE IF (my_do_mixed) THEN
1469 localized_wfn_control%set_of_states = state_loc_mixed
1470 nextra = localized_wfn_control%nextra
1471 ELSE
1472 localized_wfn_control%set_of_states = state_loc_all
1473 END IF
1474 END IF
1475
1476 localized_wfn_control%print_centers = &
1477 btest(cp_print_key_should_output(logger%iter_info, loc_print_section, &
1478 "WANNIER_CENTERS"), cp_p_file)
1479 localized_wfn_control%print_spreads = &
1480 btest(cp_print_key_should_output(logger%iter_info, loc_print_section, &
1481 "WANNIER_SPREADS"), cp_p_file)
1482 localized_wfn_control%print_cubes = &
1483 btest(cp_print_key_should_output(logger%iter_info, loc_print_section, &
1484 "WANNIER_CUBES"), cp_p_file)
1485
1486 output_unit = cp_print_key_unit_nr(logger, loc_print_section, "PROGRAM_RUN_INFO", &
1487 extension=".Log")
1488
1489 IF (output_unit > 0) THEN
1490 WRITE (unit=output_unit, fmt="(/,T2,A)") &
1491 "LOCALIZE| The spread relative to a set of orbitals is computed"
1492
1493 SELECT CASE (localized_wfn_control%set_of_states)
1494 CASE (state_loc_all)
1495 WRITE (unit=output_unit, fmt="(T2,A)") &
1496 "LOCALIZE| Orbitals to be localized: All orbitals"
1497 WRITE (unit=output_unit, fmt="(T2,A,/,T12,A,F16.8)") &
1498 "LOCALIZE| If fractional occupation, fully occupied MOs are those ", &
1499 "within occupation tolerance of ", localized_wfn_control%eps_occ
1500 CASE (state_loc_range)
1501 WRITE (unit=output_unit, fmt="(T2,A,T65,I8,A,I8)") &
1502 "LOCALIZE| Orbitals to be localized: Those with index between ", &
1503 localized_wfn_control%lu_bound_states(1, spin_xas), " and ", &
1504 localized_wfn_control%lu_bound_states(2, spin_xas)
1505 CASE (state_loc_list)
1506 WRITE (unit=output_unit, fmt="(T2,A)") &
1507 "LOCALIZE| Orbitals to be localized: Those with index in the following list"
1508 nline = localized_wfn_control%nloc_states(1)/10 + 1
1509 ind = 0
1510 DO i = 1, nline
1511 IF (ind + 10 < localized_wfn_control%nloc_states(1)) THEN
1512 WRITE (unit=output_unit, fmt="(T8,10I7)") localized_wfn_control%loc_states(ind + 1:ind + 10, 1)
1513 ind = ind + 10
1514 ELSE
1515 WRITE (unit=output_unit, fmt="(T8,10I7)") &
1516 localized_wfn_control%loc_states(ind + 1:localized_wfn_control%nloc_states(1), 1)
1517 ind = localized_wfn_control%nloc_states(1)
1518 END IF
1519 END DO
1520 CASE (energy_loc_range)
1521 WRITE (unit=output_unit, fmt="(T2,A,T65,/,f16.6,A,f16.6,A)") &
1522 "LOCALIZE| Orbitals to be localized: Those with energy in the range between ", &
1523 localized_wfn_control%lu_ene_bound(1), " and ", localized_wfn_control%lu_ene_bound(2), " a.u."
1524 CASE (state_loc_mixed)
1525 WRITE (unit=output_unit, fmt="(T2,A,I4,A)") &
1526 "LOCALIZE| Orbitals to be localized: Occupied orbitals + ", nextra, " orbitals"
1527 CASE DEFAULT
1528 WRITE (unit=output_unit, fmt="(T2,A)") &
1529 "LOCALIZE| Orbitals to be localized: None "
1530 END SELECT
1531
1532 SELECT CASE (localized_wfn_control%operator_type)
1533 CASE (op_loc_berry)
1534 WRITE (unit=output_unit, fmt="(T2,A)") &
1535 "LOCALIZE| Spread defined by the Berry phase operator "
1536 CASE (op_loc_boys)
1537 WRITE (unit=output_unit, fmt="(T2,A)") &
1538 "LOCALIZE| Spread defined by the Boys phase operator "
1539 CASE DEFAULT
1540 WRITE (unit=output_unit, fmt="(T2,A)") &
1541 "LOCALIZE| Spread defined by the Pipek phase operator "
1542 END SELECT
1543
1544 SELECT CASE (localized_wfn_control%localization_method)
1545 CASE (do_loc_jacobi)
1546 WRITE (unit=output_unit, fmt="(T2,A)") &
1547 "LOCALIZE| Optimal unitary transformation generated by Jacobi algorithm"
1548 CASE (do_loc_crazy)
1549 WRITE (unit=output_unit, fmt="(T2,A)") &
1550 "LOCALIZE| Optimal unitary transformation generated by Crazy angle algorithm"
1551 WRITE (unit=output_unit, fmt="(T2,A,F16.8)") &
1552 "LOCALIZE| maximum angle: ", localized_wfn_control%max_crazy_angle
1553 WRITE (unit=output_unit, fmt="(T2,A,F16.8)") &
1554 "LOCALIZE| scaling: ", localized_wfn_control%crazy_scale
1555 WRITE (unit=output_unit, fmt="(T2,A,L1)") &
1556 "LOCALIZE| use diag:", localized_wfn_control%crazy_use_diag
1557 CASE (do_loc_gapo)
1558 WRITE (unit=output_unit, fmt="(T2,A)") &
1559 "LOCALIZE| Optimal unitary transformation generated by gradient ascent algorithm "
1560 WRITE (unit=output_unit, fmt="(T2,A)") &
1561 "LOCALIZE| for partially occupied wannier functions"
1562 CASE (do_loc_direct)
1563 WRITE (unit=output_unit, fmt="(T2,A)") &
1564 "LOCALIZE| Optimal unitary transformation generated by direct algorithm"
1565 CASE (do_loc_l1_norm_sd)
1566 WRITE (unit=output_unit, fmt="(T2,A)") &
1567 "LOCALIZE| Optimal unitary transformation generated by "
1568 WRITE (unit=output_unit, fmt="(T2,A)") &
1569 "LOCALIZE| steepest descent algorithm applied on an approximate l1 norm"
1570 CASE (do_loc_none)
1571 WRITE (unit=output_unit, fmt="(T2,A)") &
1572 "LOCALIZE| No unitary transformation is applied"
1573 CASE (do_loc_scdm)
1574 WRITE (unit=output_unit, fmt="(T2,A)") &
1575 "LOCALIZE| Pivoted QR decomposition is used to transform coefficients"
1576 END SELECT
1577
1578 END IF ! process has output_unit
1579
1580 CALL cp_print_key_finished_output(output_unit, logger, loc_print_section, "PROGRAM_RUN_INFO")
1581
1582 ELSE
1583 localized_wfn_control%localization_method = do_loc_none
1584 localized_wfn_control%localization_method = state_loc_none
1585 localized_wfn_control%print_centers = .false.
1586 localized_wfn_control%print_spreads = .false.
1587 localized_wfn_control%print_cubes = .false.
1588 END IF
1589
1590 END SUBROUTINE read_loc_section
1591
1592! **************************************************************************************************
1593!> \brief create the center and spread array and the file names for the output
1594!> \param localized_wfn_control ...
1595!> \param nmoloc ...
1596!> \param nspins ...
1597!> \par History
1598!> 04.2005 created [MI]
1599! **************************************************************************************************
1600 SUBROUTINE set_loc_centers(localized_wfn_control, nmoloc, nspins)
1601
1602 TYPE(localized_wfn_control_type) :: localized_wfn_control
1603 INTEGER, DIMENSION(2), INTENT(IN) :: nmoloc
1604 INTEGER, INTENT(IN) :: nspins
1605
1606 INTEGER :: ispin
1607
1608 DO ispin = 1, nspins
1609 ALLOCATE (localized_wfn_control%centers_set(ispin)%array(6, nmoloc(ispin)))
1610 localized_wfn_control%centers_set(ispin)%array = 0.0_dp
1611 END DO
1612
1613 END SUBROUTINE set_loc_centers
1614
1615! **************************************************************************************************
1616!> \brief create the lists of mos that are taken into account
1617!> \param localized_wfn_control ...
1618!> \param nmoloc ...
1619!> \param nmo ...
1620!> \param nspins ...
1621!> \param my_spin ...
1622!> \par History
1623!> 04.2005 created [MI]
1624! **************************************************************************************************
1625 SUBROUTINE set_loc_wfn_lists(localized_wfn_control, nmoloc, nmo, nspins, my_spin)
1626
1627 TYPE(localized_wfn_control_type) :: localized_wfn_control
1628 INTEGER, DIMENSION(2), INTENT(IN) :: nmoloc, nmo
1629 INTEGER, INTENT(IN) :: nspins
1630 INTEGER, INTENT(IN), OPTIONAL :: my_spin
1631
1632 CHARACTER(len=*), PARAMETER :: routinen = 'set_loc_wfn_lists'
1633
1634 INTEGER :: i, ispin, max_iloc, max_nmoloc, state
1635
1636 CALL timeset(routinen, state)
1637
1638 localized_wfn_control%nloc_states(1:2) = nmoloc(1:2)
1639 max_nmoloc = max(nmoloc(1), nmoloc(2))
1640
1641 SELECT CASE (localized_wfn_control%set_of_states)
1642 CASE (state_loc_list)
1643 ! List
1644 cpassert(ASSOCIATED(localized_wfn_control%loc_states))
1645 DO ispin = 1, nspins
1646 localized_wfn_control%lu_bound_states(1, ispin) = 1
1647 localized_wfn_control%lu_bound_states(2, ispin) = nmoloc(ispin)
1648 IF (nmoloc(ispin) < 1) THEN
1649 localized_wfn_control%lu_bound_states(1, ispin) = 0
1650 localized_wfn_control%loc_states(:, ispin) = 0
1651 END IF
1652 END DO
1653 CASE (state_loc_range)
1654 ! Range
1655 ALLOCATE (localized_wfn_control%loc_states(max_nmoloc, 2))
1656 localized_wfn_control%loc_states = 0
1657 DO ispin = 1, nspins
1658 localized_wfn_control%lu_bound_states(1, ispin) = &
1659 localized_wfn_control%lu_bound_states(1, my_spin)
1660 localized_wfn_control%lu_bound_states(2, ispin) = &
1661 localized_wfn_control%lu_bound_states(1, my_spin) + nmoloc(ispin) - 1
1662 max_iloc = localized_wfn_control%lu_bound_states(2, ispin)
1663 DO i = 1, nmoloc(ispin)
1664 localized_wfn_control%loc_states(i, ispin) = localized_wfn_control%lu_bound_states(1, ispin) + i - 1
1665 END DO
1666 cpassert(max_iloc <= nmo(ispin))
1667 mark_used(nmo)
1668 END DO
1669 CASE (energy_loc_range)
1670 ! Energy
1671 ALLOCATE (localized_wfn_control%loc_states(max_nmoloc, 2))
1672 localized_wfn_control%loc_states = 0
1673 DO ispin = 1, nspins
1674 DO i = 1, nmoloc(ispin)
1675 localized_wfn_control%loc_states(i, ispin) = localized_wfn_control%lu_bound_states(1, ispin) + i - 1
1676 END DO
1677 END DO
1678 CASE (state_loc_all)
1679 ! All
1680 ALLOCATE (localized_wfn_control%loc_states(max_nmoloc, 2))
1681 localized_wfn_control%loc_states = 0
1682
1683 IF (localized_wfn_control%lu_bound_states(1, 1) == 1) THEN
1684 DO ispin = 1, nspins
1685 localized_wfn_control%lu_bound_states(1, ispin) = 1
1686 localized_wfn_control%lu_bound_states(2, ispin) = nmoloc(ispin)
1687 IF (nmoloc(ispin) < 1) localized_wfn_control%lu_bound_states(1, ispin) = 0
1688 DO i = 1, nmoloc(ispin)
1689 localized_wfn_control%loc_states(i, ispin) = i
1690 END DO
1691 END DO
1692 ELSE
1693 DO ispin = 1, nspins
1694 IF (nmoloc(ispin) < 1) localized_wfn_control%lu_bound_states(1, ispin) = 0
1695 DO i = 1, nmoloc(ispin)
1696 localized_wfn_control%loc_states(i, ispin) = &
1697 localized_wfn_control%lu_bound_states(1, ispin) + i - 1
1698 END DO
1699 END DO
1700 END IF
1701 CASE (state_loc_mixed)
1702 ! Mixed
1703 ALLOCATE (localized_wfn_control%loc_states(max_nmoloc, 2))
1704 localized_wfn_control%loc_states = 0
1705 DO ispin = 1, nspins
1706 DO i = 1, nmoloc(ispin)
1707 localized_wfn_control%loc_states(i, ispin) = i
1708 END DO
1709 END DO
1710 END SELECT
1711
1712 CALL timestop(state)
1713
1714 END SUBROUTINE set_loc_wfn_lists
1715
1716END MODULE qs_loc_utils
1717
Calculation of the moment integrals over Cartesian Gaussian-type functions.
Definition ai_moments.F:17
subroutine, public contract_cossin(cos_block, sin_block, iatom, ncoa, nsgfa, sgfa, sphi_a, ldsa, jatom, ncob, nsgfb, sgfb, sphi_b, ldsb, cosab, sinab, ldab, work, ldwork)
...
Definition ai_moments.F:261
subroutine, public cossin(la_max_set, npgfa, zeta, rpgfa, la_min_set, lb_max, npgfb, zetb, rpgfb, lb_min, rac, rbc, kvec, cosab, sinab, dcosab, dsinab)
...
Definition ai_moments.F:339
collect pointers to a block of reals
Handles all functions related to the CELL.
Definition cell_types.F:15
various utilities that regard array of different kinds: output, allocation,... maybe it is not a good...
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
DBCSR operations in CP2K.
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
Utility routines to open and close files. Tracking of preconnections.
Definition cp_files.F:16
subroutine, public open_file(file_name, file_status, file_form, file_action, file_position, file_pad, unit_number, debug, skip_get_unit_number, file_access)
Opens the requested file using a free unit number.
Definition cp_files.F:308
subroutine, public close_file(unit_number, file_status, keep_preconnection)
Close an open file given by its logical unit number. Optionally, keep the file and unit preconnected.
Definition cp_files.F:119
basic linear algebra operations for full matrices
subroutine, public cp_fm_column_scale(matrixa, scaling)
scales column i of matrix a with scaling(i)
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
Definition cp_fm_diag.F:17
subroutine, public choose_eigv_solver(matrix, eigenvectors, eigenvalues, info)
Choose the Eigensolver depending on which library is available ELPA seems to be unstable for small sy...
Definition cp_fm_diag.F:208
represent the structure of a full matrix
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
subroutine, public cp_fm_write_unformatted(fm, unit)
...
subroutine, public cp_fm_set_submatrix(fm, new_values, start_row, start_col, n_rows, n_cols, alpha, beta, transpose)
sets a submatrix of a full matrix fm(start_row:start_row+n_rows,start_col:start_col+n_cols) = alpha*o...
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
subroutine, public cp_fm_get_submatrix(fm, target_m, start_row, start_col, n_rows, n_cols, transpose)
gets a submatrix of a full matrix op(target_m)(1:n_rows,1:n_cols) =fm(start_row:start_row+n_rows,...
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
various routines to log and control the output. The idea is that decisions about where to log should ...
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
character(len=default_path_length) function, public cp_print_key_generate_filename(logger, print_key, middle_name, extension, my_local)
Utility function that returns a unit number to write the print key. Might open a file with a unique f...
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
integer, parameter, public cp_p_file
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public op_loc_pipek
integer, parameter, public do_loc_jacobi
integer, parameter, public do_loc_l1_norm_sd
integer, parameter, public state_loc_all
integer, parameter, public do_loc_none
integer, parameter, public do_loc_scdm
integer, parameter, public energy_loc_range
integer, parameter, public op_loc_berry
integer, parameter, public do_loc_crazy
integer, parameter, public state_loc_mixed
integer, parameter, public do_loc_gapo
integer, parameter, public op_loc_boys
integer, parameter, public state_loc_none
integer, parameter, public state_loc_list
integer, parameter, public state_loc_range
integer, parameter, public do_mixed
integer, parameter, public do_loc_direct
objects that represent the structure of input sections and the data contained in an input section
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
integer, parameter, public default_string_length
Definition kinds.F:57
integer, parameter, public default_path_length
Definition kinds.F:58
An array-based list which grows on demand. When the internal array is full, a new array of twice the ...
Definition list.F:24
Definition of mathematical constants and functions.
real(kind=dp), parameter, public twopi
Utility routines for the memory handling.
Interface to the message passing library MPI.
Provides Cartesian and spherical orbital pointers and indices.
integer, dimension(:), allocatable, public ncoset
basic linear algebra operations for full matrixes
Define the data structure for the particle information.
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Define the quickstep kind type and their sub types.
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_r3d_rs_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, basis_type, total_zeff_corr)
Get attributes of an atomic kind set.
New version of the module for the localization of the molecular orbitals This should be able to use d...
subroutine, public localized_wfn_control_create(localized_wfn_control)
create the localized_wfn_control_type
subroutine, public localized_wfn_control_release(localized_wfn_control)
release the localized_wfn_control_type
subroutine, public get_qs_loc_env(qs_loc_env, cell, local_molecules, localized_wfn_control, moloc_coeff, op_sm_set, op_fm_set, para_env, particle_set, weights, dim_op)
...
subroutine, public set_qs_loc_env(qs_loc_env, cell, local_molecules, localized_wfn_control, moloc_coeff, op_sm_set, op_fm_set, para_env, particle_set, weights, dim_op)
...
Some utilities for the construction of the localization environment.
subroutine, public compute_berry_operator(qs_env, cell, op_sm_set, dim_op)
Computes the Berry operator for periodic systems used to define the spread of the MOS Here the matrix...
subroutine, public set_loc_wfn_lists(localized_wfn_control, nmoloc, nmo, nspins, my_spin)
create the lists of mos that are taken into account
subroutine, public loc_write_restart(qs_loc_env, section, mo_array, coeff_localized, do_homo, evals, do_mixed)
...
subroutine, public qs_loc_env_init(qs_loc_env, localized_wfn_control, qs_env, myspin, do_localize, loc_coeff, mo_loc_history)
allocates the data, and initializes the operators
subroutine, public set_loc_centers(localized_wfn_control, nmoloc, nspins)
create the center and spread array and the file names for the output
subroutine, public qs_loc_control_init(qs_loc_env, loc_section, do_homo, do_mixed, do_xas, nloc_xas, spin_xas)
initializes everything needed for localization of the HOMOs
subroutine, public retain_history(mo_loc_history, mo_loc)
copy old mos to new ones, allocating as necessary
subroutine, public qs_loc_init(qs_env, qs_loc_env, localize_section, mos_localized, do_homo, do_mo_cubes, mo_loc_history, evals, tot_zeff_corr, do_mixed)
initializes everything needed for localization of the molecular orbitals
Localization methods such as 2x2 Jacobi rotations Steepest Decents Conjugate Gradient.
subroutine, public initialize_weights(cell, weights)
...
collects routines that perform operations directly related to MOs
subroutine, public make_mo_eig(mos, nspins, ks_rmpv, scf_control, mo_derivs, admm_env)
Calculate KS eigenvalues starting from OF MOS.
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count)
Get the components of a MO set data structure.
Define the neighbor list data types and the corresponding functionality.
subroutine, public neighbor_list_iterator_create(iterator_set, nl, search, nthread)
Neighbor list iterator functions.
subroutine, public neighbor_list_iterator_release(iterator_set)
...
integer function, public neighbor_list_iterate(iterator_set, mepos)
...
subroutine, public get_iterator_info(iterator_set, mepos, ikind, jkind, nkind, ilist, nlist, inode, nnode, iatom, jatom, r, cell)
...
module that contains the definitions of the scf types
integer, parameter, public ot_method_nr
parameters that control an scf iteration
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
represent a pointer to a 1d array
keeps the information about the structure of a full matrix
represent a full matrix
type of a logger, at the moment it contains just a print level starting at which level it should be l...
structure to store local (to a processor) ordered lists of integers.
stores all the informations relevant to an mpi environment
Provides all information about a quickstep kind.
A type that holds controlling information for the calculation of the spread of wfn and the optimizati...
contains all the info needed by quickstep to calculate the spread of a selected set of orbitals and i...