Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Some utilities for the construction of
10 : !> the localization environment
11 : !> \author MI (05-2005)
12 : ! **************************************************************************************************
13 : MODULE qs_loc_utils
14 :
15 : USE ai_moments, ONLY: contract_cossin,&
16 : cossin
17 : USE basis_set_types, ONLY: gto_basis_set_p_type,&
18 : gto_basis_set_type
19 : USE block_p_types, ONLY: block_p_type
20 : USE cell_types, ONLY: cell_type,&
21 : pbc
22 : USE cp_array_utils, ONLY: cp_1d_r_p_type
23 : USE cp_control_types, ONLY: dft_control_type
24 : USE cp_dbcsr_api, ONLY: dbcsr_copy,&
25 : dbcsr_get_block_p,&
26 : dbcsr_p_type,&
27 : dbcsr_set,&
28 : dbcsr_type
29 : USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply
30 : USE cp_files, ONLY: close_file,&
31 : open_file
32 : USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale
33 : USE cp_fm_diag, ONLY: choose_eigv_solver
34 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
35 : cp_fm_struct_release,&
36 : cp_fm_struct_type
37 : USE cp_fm_types, ONLY: &
38 : cp_fm_create, cp_fm_get_info, cp_fm_get_submatrix, cp_fm_release, cp_fm_set_all, &
39 : cp_fm_set_submatrix, cp_fm_to_fm, cp_fm_type, cp_fm_write_unformatted
40 : USE cp_log_handling, ONLY: cp_get_default_logger,&
41 : cp_logger_get_default_io_unit,&
42 : cp_logger_type,&
43 : cp_to_string
44 : USE cp_output_handling, ONLY: cp_p_file,&
45 : cp_print_key_finished_output,&
46 : cp_print_key_generate_filename,&
47 : cp_print_key_should_output,&
48 : cp_print_key_unit_nr
49 : USE distribution_1d_types, ONLY: distribution_1d_type
50 : USE input_constants, ONLY: &
51 : do_loc_crazy, do_loc_direct, do_loc_gapo, do_loc_jacobi, do_loc_l1_norm_sd, do_loc_none, &
52 : do_loc_scdm, energy_loc_range, op_loc_berry, op_loc_boys, op_loc_pipek, state_loc_all, &
53 : state_loc_list, state_loc_mixed, state_loc_none, state_loc_range
54 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
55 : section_vals_type,&
56 : section_vals_val_get
57 : USE kinds, ONLY: default_path_length,&
58 : default_string_length,&
59 : dp
60 : USE mathconstants, ONLY: twopi
61 : USE memory_utilities, ONLY: reallocate
62 : USE message_passing, ONLY: mp_para_env_type
63 : USE orbital_pointers, ONLY: ncoset
64 : USE parallel_gemm_api, ONLY: parallel_gemm
65 : USE particle_types, ONLY: particle_type
66 : USE qs_environment_types, ONLY: get_qs_env,&
67 : qs_environment_type
68 : USE qs_kind_types, ONLY: get_qs_kind,&
69 : get_qs_kind_set,&
70 : qs_kind_type
71 : USE qs_loc_types, ONLY: get_qs_loc_env,&
72 : localized_wfn_control_create,&
73 : localized_wfn_control_release,&
74 : localized_wfn_control_type,&
75 : qs_loc_env_type,&
76 : set_qs_loc_env
77 : USE qs_localization_methods, ONLY: initialize_weights
78 : USE qs_mo_methods, ONLY: make_mo_eig
79 : USE qs_mo_types, ONLY: get_mo_set,&
80 : mo_set_type
81 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
82 : neighbor_list_iterate,&
83 : neighbor_list_iterator_create,&
84 : neighbor_list_iterator_p_type,&
85 : neighbor_list_iterator_release,&
86 : neighbor_list_set_p_type
87 : USE qs_scf_types, ONLY: ot_method_nr
88 : USE scf_control_types, ONLY: scf_control_type
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 ***
100 : PUBLIC :: qs_loc_env_init, loc_write_restart, &
101 : retain_history, qs_loc_init, compute_berry_operator, &
102 : set_loc_centers, set_loc_wfn_lists, qs_loc_control_init
103 :
104 : CONTAINS
105 :
106 : ! **************************************************************************************************
107 : !> \brief copy old mos to new ones, allocating as necessary
108 : !> \param mo_loc_history ...
109 : !> \param mo_loc ...
110 : ! **************************************************************************************************
111 10 : 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 10 : CALL timeset(routineN, handle)
121 :
122 10 : IF (.NOT. ASSOCIATED(mo_loc_history)) THEN
123 8 : ALLOCATE (mo_loc_history(SIZE(mo_loc)))
124 4 : DO i = 1, SIZE(mo_loc_history)
125 4 : CALL cp_fm_create(mo_loc_history(i), mo_loc(i)%matrix_struct)
126 : END DO
127 : END IF
128 :
129 20 : DO i = 1, SIZE(mo_loc_history)
130 10 : CALL cp_fm_get_info(mo_loc_history(i), ncol_global=ncol_hist)
131 10 : CALL cp_fm_get_info(mo_loc(i), ncol_global=ncol_loc)
132 10 : CPASSERT(ncol_hist == ncol_loc)
133 30 : CALL cp_fm_to_fm(mo_loc(i), mo_loc_history(i))
134 : END DO
135 :
136 10 : CALL timestop(handle)
137 :
138 10 : 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 8 : 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 8 : CALL timeset(routineN, handle)
160 :
161 8 : CALL cp_fm_get_info(mo_new, nrow_global=nrow, ncol_global=ncol)
162 8 : CALL cp_fm_get_info(mo_ref, ncol_global=ncol_ref)
163 8 : CPASSERT(ncol == ncol_ref)
164 :
165 8 : NULLIFY (fm_struct_tmp)
166 8 : 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 8 : context=mo_new%matrix_struct%context)
171 8 : CALL cp_fm_create(o1, fm_struct_tmp)
172 8 : CALL cp_fm_create(o2, fm_struct_tmp)
173 8 : CALL cp_fm_create(o3, fm_struct_tmp)
174 8 : CALL cp_fm_create(o4, fm_struct_tmp)
175 8 : CALL cp_fm_struct_release(fm_struct_tmp)
176 :
177 : ! o1 = (mo_new)^T matrix_S mo_ref
178 8 : CALL cp_dbcsr_sm_fm_multiply(matrix_S, mo_ref, smo, ncol)
179 8 : 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 8 : 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 24 : ALLOCATE (eigenvalues(ncol))
186 8 : CALL choose_eigv_solver(o2, o3, eigenvalues)
187 8 : CALL cp_fm_to_fm(o3, o4)
188 72 : eigenvalues(:) = 1.0_dp/SQRT(eigenvalues(:))
189 8 : CALL cp_fm_column_scale(o4, eigenvalues)
190 8 : 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 8 : 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 8 : CALL parallel_gemm('N', 'N', nrow, ncol, ncol, 1.0_dp, mo_new, o3, 0.0_dp, smo)
197 8 : 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 8 : CALL cp_fm_release(o1)
206 8 : CALL cp_fm_release(o2)
207 8 : CALL cp_fm_release(o3)
208 8 : CALL cp_fm_release(o4)
209 8 : CALL cp_fm_release(smo)
210 :
211 8 : CALL timestop(handle)
212 :
213 32 : 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 968 : SUBROUTINE qs_loc_env_init(qs_loc_env, localized_wfn_control, qs_env, myspin, do_localize, &
231 484 : 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 : LOGICAL :: loc_coeff_spin_resolved
248 : REAL(KIND=dp) :: my_occ, occ_imo
249 484 : REAL(KIND=dp), DIMENSION(:), POINTER :: occupations
250 484 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: vecbuffer
251 : TYPE(cell_type), POINTER :: cell
252 : TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct
253 484 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: moloc_coeff
254 : TYPE(cp_fm_type), POINTER :: mat_ptr, mo_coeff
255 484 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
256 : TYPE(distribution_1d_type), POINTER :: local_molecules
257 484 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
258 : TYPE(mp_para_env_type), POINTER :: para_env
259 484 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
260 :
261 484 : CALL timeset(routineN, handle)
262 :
263 484 : NULLIFY (mos, matrix_s, moloc_coeff, particle_set, para_env, cell, &
264 484 : local_molecules, occupations, mat_ptr)
265 484 : IF (PRESENT(do_localize)) qs_loc_env%do_localize = do_localize
266 484 : IF (qs_loc_env%do_localize) THEN
267 : CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s, cell=cell, &
268 : local_molecules=local_molecules, particle_set=particle_set, &
269 484 : para_env=para_env, mos=mos)
270 484 : nspins = SIZE(mos, 1)
271 484 : loc_coeff_spin_resolved = .FALSE.
272 484 : IF (PRESENT(loc_coeff)) THEN
273 308 : loc_coeff_spin_resolved = nspins*2 == SIZE(loc_coeff)
274 : END IF
275 484 : s_spin = 1
276 484 : l_spin = nspins
277 484 : IF (PRESENT(myspin)) THEN
278 162 : s_spin = myspin
279 162 : l_spin = myspin
280 : END IF
281 484 : IF (loc_coeff_spin_resolved) THEN
282 166 : ALLOCATE (moloc_coeff(s_spin:s_spin + 2*(l_spin - s_spin) + 1))
283 : ELSE
284 1948 : ALLOCATE (moloc_coeff(s_spin:l_spin))
285 : END IF
286 1102 : DO ispin = s_spin, l_spin
287 618 : NULLIFY (tmp_fm_struct, mo_coeff)
288 618 : CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, nao=nao, nmo=nmo)
289 618 : nmosub = localized_wfn_control%nloc_states(ispin)
290 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, &
291 618 : ncol_global=nmosub, para_env=para_env, context=mo_coeff%matrix_struct%context)
292 618 : IF (loc_coeff_spin_resolved) THEN
293 44 : CALL cp_fm_create(moloc_coeff(2*ispin - 1), tmp_fm_struct)
294 44 : CALL cp_fm_create(moloc_coeff(2*ispin), tmp_fm_struct)
295 : ELSE
296 574 : CALL cp_fm_create(moloc_coeff(ispin), tmp_fm_struct)
297 : END IF
298 :
299 : CALL cp_fm_get_info(moloc_coeff(ispin), nrow_global=naosub, &
300 618 : ncol_global=nmosub)
301 618 : CPASSERT(nao == naosub)
302 618 : IF ((localized_wfn_control%do_homo) .OR. &
303 : (localized_wfn_control%set_of_states == state_loc_mixed)) THEN
304 606 : CPASSERT(nmo >= nmosub)
305 : ELSE
306 12 : CPASSERT(nao - nmo >= nmosub)
307 : END IF
308 618 : CALL cp_fm_set_all(moloc_coeff(ispin), 0.0_dp)
309 2338 : CALL cp_fm_struct_release(tmp_fm_struct)
310 : END DO ! ispin
311 : ! Copy the submatrix
312 :
313 484 : IF (PRESENT(loc_coeff)) ALLOCATE (mat_ptr)
314 :
315 1102 : DO ispin = s_spin, l_spin
316 : CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, &
317 618 : occupation_numbers=occupations, nao=nao, nmo=nmo)
318 618 : lb = localized_wfn_control%lu_bound_states(1, ispin)
319 618 : ub = localized_wfn_control%lu_bound_states(2, ispin)
320 :
321 618 : IF (PRESENT(loc_coeff)) THEN
322 428 : mat_ptr = loc_coeff(ispin)
323 : ELSE
324 190 : mat_ptr => mo_coeff
325 : END IF
326 618 : IF ((localized_wfn_control%set_of_states == state_loc_list) .OR. &
327 : (localized_wfn_control%set_of_states == state_loc_mixed)) THEN
328 444 : ALLOCATE (vecbuffer(1, nao))
329 148 : IF (localized_wfn_control%do_homo) THEN
330 134 : my_occ = occupations(localized_wfn_control%loc_states(1, ispin))
331 : END IF
332 148 : nmosub = SIZE(localized_wfn_control%loc_states, 1)
333 148 : CPASSERT(nmosub > 0)
334 148 : imoloc = 0
335 934 : DO i = lb, ub
336 : ! Get the index in the subset
337 786 : imoloc = imoloc + 1
338 : ! Get the index in the full set
339 786 : imo = localized_wfn_control%loc_states(i, ispin)
340 786 : IF (localized_wfn_control%do_homo) THEN
341 652 : occ_imo = occupations(imo)
342 652 : IF (ABS(occ_imo - my_occ) > localized_wfn_control%eps_occ) THEN
343 0 : IF (localized_wfn_control%localization_method /= do_loc_none) THEN
344 : CALL cp_abort(__LOCATION__, &
345 : "States with different occupations "// &
346 0 : "cannot be rotated together")
347 : END IF
348 : END IF
349 : END IF
350 : ! Take the imo vector from the full set and copy in the imoloc vector of the subset
351 : CALL cp_fm_get_submatrix(mat_ptr, vecbuffer, 1, imo, &
352 786 : nao, 1, transpose=.TRUE.)
353 : CALL cp_fm_set_submatrix(moloc_coeff(ispin), vecbuffer, 1, imoloc, &
354 934 : nao, 1, transpose=.TRUE.)
355 : END DO
356 148 : DEALLOCATE (vecbuffer)
357 : ELSE
358 470 : my_occ = occupations(lb)
359 470 : occ_imo = occupations(ub)
360 470 : IF (ABS(occ_imo - my_occ) > localized_wfn_control%eps_occ) THEN
361 0 : IF (localized_wfn_control%localization_method /= do_loc_none) THEN
362 : CALL cp_abort(__LOCATION__, &
363 : "States with different occupations "// &
364 0 : "cannot be rotated together")
365 : END IF
366 : END IF
367 470 : nmosub = localized_wfn_control%nloc_states(ispin)
368 :
369 470 : IF (loc_coeff_spin_resolved) THEN
370 44 : CALL cp_fm_to_fm(loc_coeff(2*ispin - 1), moloc_coeff(2*ispin - 1))
371 44 : CALL cp_fm_to_fm(loc_coeff(2*ispin), moloc_coeff(2*ispin))
372 : ELSE
373 426 : CALL cp_fm_to_fm(mat_ptr, moloc_coeff(ispin), nmosub, lb, 1)
374 : END IF
375 : END IF
376 :
377 : ! we have the mo's to be localized now, see if we can rotate them according to the history
378 : ! only do that if we have a history of course. The history is filled
379 1720 : IF (PRESENT(mo_loc_history)) THEN
380 104 : IF (localized_wfn_control%use_history .AND. ASSOCIATED(mo_loc_history)) THEN
381 : CALL rotate_state_to_ref(moloc_coeff(ispin), &
382 8 : mo_loc_history(ispin), matrix_s(1)%matrix)
383 : END IF
384 : END IF
385 :
386 : END DO
387 :
388 484 : IF (PRESENT(loc_coeff)) DEALLOCATE (mat_ptr)
389 :
390 : CALL set_qs_loc_env(qs_loc_env=qs_loc_env, cell=cell, local_molecules=local_molecules, &
391 : moloc_coeff=moloc_coeff, particle_set=particle_set, para_env=para_env, &
392 484 : localized_wfn_control=localized_wfn_control)
393 :
394 : ! Prepare the operators
395 484 : NULLIFY (tmp_fm_struct, mo_coeff)
396 1452 : nmosub = MAXVAL(localized_wfn_control%nloc_states)
397 484 : CALL get_mo_set(mos(1), mo_coeff=mo_coeff)
398 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nmosub, &
399 484 : ncol_global=nmosub, para_env=para_env, context=mo_coeff%matrix_struct%context)
400 :
401 484 : IF (localized_wfn_control%operator_type == op_loc_berry) THEN
402 478 : IF (qs_loc_env%cell%orthorhombic) THEN
403 466 : dim_op = 3
404 : ELSE
405 12 : dim_op = 6
406 : END IF
407 478 : CALL set_qs_loc_env(qs_loc_env=qs_loc_env, dim_op=dim_op)
408 5844 : ALLOCATE (qs_loc_env%op_sm_set(2, dim_op))
409 1948 : DO i = 1, dim_op
410 4888 : DO j = 1, SIZE(qs_loc_env%op_sm_set, 1)
411 2940 : NULLIFY (qs_loc_env%op_sm_set(j, i)%matrix)
412 2940 : ALLOCATE (qs_loc_env%op_sm_set(j, i)%matrix)
413 : CALL dbcsr_copy(qs_loc_env%op_sm_set(j, i)%matrix, matrix_s(1)%matrix, &
414 2940 : name="qs_loc_env%op_sm_"//TRIM(ADJUSTL(cp_to_string(j)))//"-"//TRIM(ADJUSTL(cp_to_string(i))))
415 4410 : CALL dbcsr_set(qs_loc_env%op_sm_set(j, i)%matrix, 0.0_dp)
416 : END DO
417 : END DO
418 :
419 6 : ELSE IF (localized_wfn_control%operator_type == op_loc_pipek) THEN
420 6 : natoms = SIZE(qs_loc_env%particle_set, 1)
421 96 : ALLOCATE (qs_loc_env%op_fm_set(natoms, 1))
422 6 : CALL set_qs_loc_env(qs_loc_env=qs_loc_env, dim_op=natoms)
423 12 : DO ispin = 1, SIZE(qs_loc_env%op_fm_set, 2)
424 6 : CALL get_mo_set(mos(ispin), nmo=nmo)
425 84 : DO iatom = 1, natoms
426 72 : CALL cp_fm_create(qs_loc_env%op_fm_set(iatom, ispin), tmp_fm_struct)
427 :
428 72 : CALL cp_fm_get_info(qs_loc_env%op_fm_set(iatom, ispin), nrow_global=nmosub)
429 72 : CPASSERT(nmo >= nmosub)
430 150 : CALL cp_fm_set_all(qs_loc_env%op_fm_set(iatom, ispin), 0.0_dp)
431 : END DO ! iatom
432 : END DO ! ispin
433 : ELSE
434 0 : CPABORT("Type of operator not implemented")
435 : END IF
436 484 : CALL cp_fm_struct_release(tmp_fm_struct)
437 :
438 484 : IF (localized_wfn_control%operator_type == op_loc_berry) THEN
439 :
440 478 : CALL initialize_weights(qs_loc_env%cell, qs_loc_env%weights)
441 :
442 478 : CALL get_berry_operator(qs_loc_env, qs_env)
443 :
444 : ELSE IF (localized_wfn_control%operator_type == op_loc_pipek) THEN
445 :
446 : !! here we don't have to do anything
447 : !! CALL get_pipek_mezey_operator ( qs_loc_env, qs_env )
448 :
449 : END IF
450 :
451 484 : qs_loc_env%molecular_states = .FALSE.
452 484 : qs_loc_env%wannier_states = .FALSE.
453 : END IF
454 484 : CALL timestop(handle)
455 :
456 484 : END SUBROUTINE qs_loc_env_init
457 :
458 : ! **************************************************************************************************
459 : !> \brief A wrapper to compute the Berry operator for periodic systems
460 : !> \param qs_loc_env new environment for the localization calculations
461 : !> \param qs_env the qs_env in which the qs_env lives
462 : !> \par History
463 : !> 04.2005 created [MI]
464 : !> 04.2018 modified [RZK, ZL]
465 : !> \author MI
466 : ! **************************************************************************************************
467 478 : SUBROUTINE get_berry_operator(qs_loc_env, qs_env)
468 : TYPE(qs_loc_env_type), POINTER :: qs_loc_env
469 : TYPE(qs_environment_type), POINTER :: qs_env
470 :
471 : CHARACTER(len=*), PARAMETER :: routineN = 'get_berry_operator'
472 :
473 : INTEGER :: dim_op, handle
474 : TYPE(cell_type), POINTER :: cell
475 478 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: op_sm_set
476 :
477 478 : CALL timeset(routineN, handle)
478 :
479 478 : NULLIFY (cell, op_sm_set)
480 : CALL get_qs_loc_env(qs_loc_env=qs_loc_env, op_sm_set=op_sm_set, &
481 478 : cell=cell, dim_op=dim_op)
482 478 : CALL compute_berry_operator(qs_env, cell, op_sm_set, dim_op)
483 :
484 478 : CALL timestop(handle)
485 478 : END SUBROUTINE get_berry_operator
486 :
487 : ! **************************************************************************************************
488 : !> \brief Computes the Berry operator for periodic systems
489 : !> used to define the spread of the MOS
490 : !> Here the matrix elements of the type <mu|cos(kr)|nu> and <mu|sin(kr)|nu>
491 : !> are computed, where mu and nu are the contracted basis functions.
492 : !> Namely the Berry operator is exp(ikr)
493 : !> k is defined somewhere
494 : !> the pair lists are exploited and sparse matrixes are constructed
495 : !> \param qs_env the qs_env in which the qs_env lives
496 : !> \param cell ...
497 : !> \param op_sm_set ...
498 : !> \param dim_op ...
499 : !> \par History
500 : !> 04.2005 created [MI]
501 : !> 04.2018 wrapped old code [RZK, ZL]
502 : !> \author MI
503 : !> \note
504 : !> The intgrals are computed analytically using the primitives GTO
505 : !> The contraction is performed block-wise
506 : ! **************************************************************************************************
507 504 : SUBROUTINE compute_berry_operator(qs_env, cell, op_sm_set, dim_op)
508 : TYPE(qs_environment_type), POINTER :: qs_env
509 : TYPE(cell_type), POINTER :: cell
510 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: op_sm_set
511 : INTEGER :: dim_op
512 :
513 : CHARACTER(len=*), PARAMETER :: routineN = 'compute_berry_operator'
514 :
515 : INTEGER :: handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, last_jatom, &
516 : ldab, ldsa, ldsb, ldwork, maxl, ncoa, ncob, nkind, nrow, nseta, nsetb, sgfa, sgfb
517 : INTEGER, DIMENSION(3) :: perd0
518 504 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
519 504 : npgfb, nsgfa, nsgfb
520 504 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
521 : LOGICAL :: found, new_atom_b
522 : REAL(KIND=dp) :: dab, kvec(3), rab2, vector_k(3, 6)
523 : REAL(KIND=dp), DIMENSION(3) :: ra, rab, rb
524 504 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
525 504 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: cosab, rpgfa, rpgfb, sinab, sphi_a, &
526 504 : sphi_b, work, zeta, zetb
527 504 : TYPE(block_p_type), DIMENSION(:), POINTER :: op_cos, op_sin
528 504 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
529 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
530 : TYPE(neighbor_list_iterator_p_type), &
531 504 : DIMENSION(:), POINTER :: nl_iterator
532 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
533 504 : POINTER :: sab_orb
534 504 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
535 504 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
536 : TYPE(qs_kind_type), POINTER :: qs_kind
537 :
538 504 : CALL timeset(routineN, handle)
539 504 : NULLIFY (qs_kind, qs_kind_set)
540 504 : NULLIFY (particle_set)
541 504 : NULLIFY (sab_orb)
542 : NULLIFY (cosab, sinab, work)
543 504 : NULLIFY (la_max, la_min, lb_max, lb_min, npgfa, npgfb, nsgfa, nsgfb)
544 504 : NULLIFY (set_radius_a, set_radius_b, rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb)
545 :
546 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, &
547 504 : particle_set=particle_set, sab_orb=sab_orb)
548 :
549 504 : nkind = SIZE(qs_kind_set)
550 :
551 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
552 504 : maxco=ldwork, maxlgto=maxl)
553 504 : ldab = ldwork
554 2016 : ALLOCATE (cosab(ldab, ldab))
555 136612 : cosab = 0.0_dp
556 1512 : ALLOCATE (sinab(ldab, ldab))
557 136612 : sinab = 0.0_dp
558 1512 : ALLOCATE (work(ldwork, ldwork))
559 136612 : work = 0.0_dp
560 :
561 3060 : ALLOCATE (op_cos(dim_op))
562 2556 : ALLOCATE (op_sin(dim_op))
563 2052 : DO i = 1, dim_op
564 1548 : NULLIFY (op_cos(i)%block)
565 2052 : NULLIFY (op_sin(i)%block)
566 : END DO
567 :
568 504 : kvec = 0.0_dp
569 504 : vector_k = 0.0_dp
570 2016 : vector_k(:, 1) = twopi*cell%h_inv(1, :)
571 2016 : vector_k(:, 2) = twopi*cell%h_inv(2, :)
572 2016 : vector_k(:, 3) = twopi*cell%h_inv(3, :)
573 2016 : vector_k(:, 4) = twopi*(cell%h_inv(1, :) + cell%h_inv(2, :))
574 2016 : vector_k(:, 5) = twopi*(cell%h_inv(1, :) + cell%h_inv(3, :))
575 2016 : vector_k(:, 6) = twopi*(cell%h_inv(2, :) + cell%h_inv(3, :))
576 :
577 : ! This operator can be used only for periodic systems
578 : ! If an isolated system is used the periodicity is overimposed
579 2016 : perd0(1:3) = cell%perd(1:3)
580 2016 : cell%perd(1:3) = 1
581 :
582 2424 : ALLOCATE (basis_set_list(nkind))
583 1416 : DO ikind = 1, nkind
584 912 : qs_kind => qs_kind_set(ikind)
585 912 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a)
586 1416 : IF (ASSOCIATED(basis_set_a)) THEN
587 912 : basis_set_list(ikind)%gto_basis_set => basis_set_a
588 : ELSE
589 0 : NULLIFY (basis_set_list(ikind)%gto_basis_set)
590 : END IF
591 : END DO
592 504 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
593 70256 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
594 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
595 69752 : iatom=iatom, jatom=jatom, r=rab)
596 69752 : basis_set_a => basis_set_list(ikind)%gto_basis_set
597 69752 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
598 69752 : basis_set_b => basis_set_list(jkind)%gto_basis_set
599 69752 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
600 69752 : ra = pbc(particle_set(iatom)%r, cell)
601 : ! basis ikind
602 69752 : first_sgfa => basis_set_a%first_sgf
603 69752 : la_max => basis_set_a%lmax
604 69752 : la_min => basis_set_a%lmin
605 69752 : npgfa => basis_set_a%npgf
606 69752 : nseta = basis_set_a%nset
607 69752 : nsgfa => basis_set_a%nsgf_set
608 69752 : rpgfa => basis_set_a%pgf_radius
609 69752 : set_radius_a => basis_set_a%set_radius
610 69752 : sphi_a => basis_set_a%sphi
611 69752 : zeta => basis_set_a%zet
612 : ! basis jkind
613 69752 : first_sgfb => basis_set_b%first_sgf
614 69752 : lb_max => basis_set_b%lmax
615 69752 : lb_min => basis_set_b%lmin
616 69752 : npgfb => basis_set_b%npgf
617 69752 : nsetb = basis_set_b%nset
618 69752 : nsgfb => basis_set_b%nsgf_set
619 69752 : rpgfb => basis_set_b%pgf_radius
620 69752 : set_radius_b => basis_set_b%set_radius
621 69752 : sphi_b => basis_set_b%sphi
622 69752 : zetb => basis_set_b%zet
623 :
624 69752 : ldsa = SIZE(sphi_a, 1)
625 69752 : ldsb = SIZE(sphi_b, 1)
626 69752 : IF (inode == 1) last_jatom = 0
627 :
628 279008 : rb = rab + ra
629 :
630 69752 : IF (jatom /= last_jatom) THEN
631 : new_atom_b = .TRUE.
632 : last_jatom = jatom
633 : ELSE
634 : new_atom_b = .FALSE.
635 : END IF
636 :
637 : IF (new_atom_b) THEN
638 18941 : IF (iatom <= jatom) THEN
639 9900 : irow = iatom
640 9900 : icol = jatom
641 : ELSE
642 9041 : irow = jatom
643 9041 : icol = iatom
644 : END IF
645 :
646 76232 : DO i = 1, dim_op
647 57291 : NULLIFY (op_cos(i)%block)
648 : CALL dbcsr_get_block_p(matrix=op_sm_set(1, i)%matrix, &
649 57291 : row=irow, col=icol, block=op_cos(i)%block, found=found)
650 57291 : NULLIFY (op_sin(i)%block)
651 : CALL dbcsr_get_block_p(matrix=op_sm_set(2, i)%matrix, &
652 76232 : row=irow, col=icol, block=op_sin(i)%block, found=found)
653 : END DO
654 : END IF ! new_atom_b
655 :
656 69752 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
657 69752 : dab = SQRT(rab2)
658 :
659 69752 : nrow = 0
660 212568 : DO iset = 1, nseta
661 :
662 142312 : ncoa = npgfa(iset)*ncoset(la_max(iset))
663 142312 : sgfa = first_sgfa(1, iset)
664 :
665 461162 : DO jset = 1, nsetb
666 :
667 318850 : ncob = npgfb(jset)*ncoset(lb_max(jset))
668 318850 : sgfb = first_sgfb(1, jset)
669 :
670 461162 : IF (set_radius_a(iset) + set_radius_b(jset) >= dab) THEN
671 :
672 : ! *** Calculate the primitive overlap integrals ***
673 608842 : DO i = 1, dim_op
674 1844340 : kvec(1:3) = vector_k(1:3, i)
675 162735237 : cosab = 0.0_dp
676 162735237 : sinab = 0.0_dp
677 : CALL cossin(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), &
678 : la_min(iset), lb_max(jset), npgfb(jset), zetb(:, jset), &
679 : rpgfb(:, jset), lb_min(jset), &
680 461085 : ra, rb, kvec, cosab, sinab)
681 : CALL contract_cossin(op_cos(i)%block, op_sin(i)%block, &
682 : iatom, ncoa, nsgfa(iset), sgfa, sphi_a, ldsa, &
683 : jatom, ncob, nsgfb(jset), sgfb, sphi_b, ldsb, &
684 608842 : cosab, sinab, ldab, work, ldwork)
685 : END DO
686 :
687 : END IF ! >= dab
688 :
689 : END DO ! jset
690 :
691 212064 : nrow = nrow + ncoa
692 :
693 : END DO ! iset
694 :
695 : END DO
696 504 : CALL neighbor_list_iterator_release(nl_iterator)
697 :
698 : ! Set back the correct periodicity
699 2016 : cell%perd(1:3) = perd0(1:3)
700 :
701 2052 : DO i = 1, dim_op
702 1548 : NULLIFY (op_cos(i)%block)
703 2052 : NULLIFY (op_sin(i)%block)
704 : END DO
705 504 : DEALLOCATE (op_cos, op_sin)
706 :
707 504 : DEALLOCATE (cosab, sinab, work, basis_set_list)
708 :
709 504 : CALL timestop(handle)
710 1008 : END SUBROUTINE compute_berry_operator
711 :
712 : ! **************************************************************************************************
713 : !> \brief ...
714 : !> \param qs_loc_env ...
715 : !> \param section ...
716 : !> \param mo_array ...
717 : !> \param coeff_localized ...
718 : !> \param do_homo ...
719 : !> \param evals ...
720 : !> \param do_mixed ...
721 : ! **************************************************************************************************
722 298 : SUBROUTINE loc_write_restart(qs_loc_env, section, mo_array, coeff_localized, &
723 : do_homo, evals, do_mixed)
724 : TYPE(qs_loc_env_type), POINTER :: qs_loc_env
725 : TYPE(section_vals_type), POINTER :: section
726 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array
727 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: coeff_localized
728 : LOGICAL, INTENT(IN) :: do_homo
729 : TYPE(cp_1d_r_p_type), DIMENSION(:), OPTIONAL, &
730 : POINTER :: evals
731 : LOGICAL, INTENT(IN), OPTIONAL :: do_mixed
732 :
733 : CHARACTER(LEN=*), PARAMETER :: routineN = 'loc_write_restart'
734 :
735 : CHARACTER(LEN=default_path_length) :: filename
736 : CHARACTER(LEN=default_string_length) :: my_middle
737 : INTEGER :: handle, ispin, max_block, nao, nloc, &
738 : nmo, output_unit, rst_unit
739 : LOGICAL :: my_do_mixed
740 : TYPE(cp_logger_type), POINTER :: logger
741 : TYPE(section_vals_type), POINTER :: print_key
742 :
743 298 : CALL timeset(routineN, handle)
744 298 : NULLIFY (logger)
745 298 : logger => cp_get_default_logger()
746 298 : output_unit = cp_logger_get_default_io_unit(logger)
747 :
748 298 : IF (qs_loc_env%do_localize) THEN
749 :
750 282 : print_key => section_vals_get_subs_vals(section, "LOC_RESTART")
751 282 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
752 : section, "LOC_RESTART"), &
753 : cp_p_file)) THEN
754 :
755 : ! Open file
756 : rst_unit = -1
757 :
758 30 : my_do_mixed = .FALSE.
759 30 : IF (PRESENT(do_mixed)) my_do_mixed = do_mixed
760 30 : IF (do_homo) THEN
761 30 : my_middle = "LOC_HOMO"
762 0 : ELSE IF (my_do_mixed) THEN
763 0 : my_middle = "LOC_MIXED"
764 : ELSE
765 0 : my_middle = "LOC_LUMO"
766 : END IF
767 :
768 : rst_unit = cp_print_key_unit_nr(logger, section, "LOC_RESTART", &
769 : extension=".wfn", file_status="REPLACE", file_action="WRITE", &
770 30 : file_form="UNFORMATTED", middle_name=TRIM(my_middle))
771 :
772 30 : IF (rst_unit > 0) filename = cp_print_key_generate_filename(logger, print_key, &
773 : middle_name=TRIM(my_middle), extension=".wfn", &
774 15 : my_local=.FALSE.)
775 :
776 30 : IF (output_unit > 0) THEN
777 : WRITE (UNIT=output_unit, FMT="(/,T2,A, A/)") &
778 15 : "LOCALIZATION| Write restart file for the localized MOS : ", &
779 30 : TRIM(filename)
780 : END IF
781 :
782 30 : IF (rst_unit > 0) THEN
783 15 : WRITE (rst_unit) qs_loc_env%localized_wfn_control%set_of_states
784 105 : WRITE (rst_unit) qs_loc_env%localized_wfn_control%lu_bound_states
785 45 : WRITE (rst_unit) qs_loc_env%localized_wfn_control%nloc_states
786 : END IF
787 :
788 70 : DO ispin = 1, SIZE(coeff_localized)
789 30 : ASSOCIATE (mo_coeff => coeff_localized(ispin))
790 40 : CALL cp_fm_get_info(mo_coeff, nrow_global=nao, ncol_global=nmo, ncol_block=max_block)
791 40 : nloc = qs_loc_env%localized_wfn_control%nloc_states(ispin)
792 40 : IF (rst_unit > 0) THEN
793 198 : WRITE (rst_unit) qs_loc_env%localized_wfn_control%loc_states(1:nloc, ispin)
794 20 : IF (do_homo .OR. my_do_mixed) THEN
795 20 : WRITE (rst_unit) nmo, &
796 20 : mo_array(ispin)%homo, &
797 20 : mo_array(ispin)%lfomo, &
798 40 : mo_array(ispin)%nelectron
799 456 : WRITE (rst_unit) mo_array(ispin)%eigenvalues(1:nmo), &
800 476 : mo_array(ispin)%occupation_numbers(1:nmo)
801 : ELSE
802 0 : WRITE (rst_unit) nmo
803 0 : WRITE (rst_unit) evals(ispin)%array(1:nmo)
804 : END IF
805 : END IF
806 :
807 80 : CALL cp_fm_write_unformatted(mo_coeff, rst_unit)
808 : END ASSOCIATE
809 :
810 : END DO
811 :
812 : CALL cp_print_key_finished_output(rst_unit, logger, section, &
813 30 : "LOC_RESTART")
814 : END IF
815 :
816 : END IF
817 :
818 298 : CALL timestop(handle)
819 :
820 298 : END SUBROUTINE loc_write_restart
821 :
822 : ! **************************************************************************************************
823 : !> \brief ...
824 : !> \param qs_loc_env ...
825 : !> \param mos ...
826 : !> \param mos_localized ...
827 : !> \param section ...
828 : !> \param section2 ...
829 : !> \param para_env ...
830 : !> \param do_homo ...
831 : !> \param restart_found ...
832 : !> \param evals ...
833 : !> \param do_mixed ...
834 : ! **************************************************************************************************
835 6 : SUBROUTINE loc_read_restart(qs_loc_env, mos, mos_localized, section, section2, para_env, &
836 : do_homo, restart_found, evals, do_mixed)
837 :
838 : TYPE(qs_loc_env_type), POINTER :: qs_loc_env
839 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
840 : TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT) :: mos_localized
841 : TYPE(section_vals_type), POINTER :: section, section2
842 : TYPE(mp_para_env_type), POINTER :: para_env
843 : LOGICAL, INTENT(IN) :: do_homo
844 : LOGICAL, INTENT(INOUT) :: restart_found
845 : TYPE(cp_1d_r_p_type), DIMENSION(:), OPTIONAL, &
846 : POINTER :: evals
847 : LOGICAL, INTENT(IN), OPTIONAL :: do_mixed
848 :
849 : CHARACTER(len=*), PARAMETER :: routineN = 'loc_read_restart'
850 :
851 : CHARACTER(LEN=25) :: fname_key
852 : CHARACTER(LEN=default_path_length) :: filename
853 : CHARACTER(LEN=default_string_length) :: my_middle
854 : INTEGER :: handle, homo_read, i, ispin, lfomo_read, max_nloc, n_rep_val, nao, &
855 : nelectron_read, nloc, nmo, nmo_read, nspin, output_unit, rst_unit
856 : LOGICAL :: file_exists, my_do_mixed
857 6 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eig_read, occ_read
858 6 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: vecbuffer
859 : TYPE(cp_logger_type), POINTER :: logger
860 : TYPE(section_vals_type), POINTER :: print_key
861 :
862 6 : CALL timeset(routineN, handle)
863 :
864 6 : logger => cp_get_default_logger()
865 :
866 6 : nspin = SIZE(mos_localized)
867 6 : nao = mos(1)%nao
868 6 : rst_unit = -1
869 :
870 : output_unit = cp_print_key_unit_nr(logger, section2, &
871 6 : "PROGRAM_RUN_INFO", extension=".Log")
872 :
873 6 : my_do_mixed = .FALSE.
874 6 : IF (PRESENT(do_mixed)) my_do_mixed = do_mixed
875 6 : IF (do_homo) THEN
876 6 : fname_key = "LOCHOMO_RESTART_FILE_NAME"
877 0 : ELSE IF (my_do_mixed) THEN
878 0 : fname_key = "LOCMIXD_RESTART_FILE_NAME"
879 : ELSE
880 0 : fname_key = "LOCLUMO_RESTART_FILE_NAME"
881 0 : IF (.NOT. PRESENT(evals)) THEN
882 0 : CPABORT("Missing argument to localize unoccupied states.")
883 : END IF
884 : END IF
885 :
886 6 : file_exists = .FALSE.
887 6 : CALL section_vals_val_get(section, fname_key, n_rep_val=n_rep_val)
888 6 : IF (n_rep_val > 0) THEN
889 0 : CALL section_vals_val_get(section, fname_key, c_val=filename)
890 : ELSE
891 :
892 6 : print_key => section_vals_get_subs_vals(section2, "LOC_RESTART")
893 6 : IF (do_homo) THEN
894 6 : my_middle = "LOC_HOMO"
895 0 : ELSE IF (my_do_mixed) THEN
896 0 : my_middle = "LOC_MIXED"
897 : ELSE
898 0 : my_middle = "LOC_LUMO"
899 : END IF
900 : filename = cp_print_key_generate_filename(logger, print_key, &
901 : middle_name=TRIM(my_middle), extension=".wfn", &
902 6 : my_local=.FALSE.)
903 : END IF
904 :
905 6 : IF (para_env%is_source()) INQUIRE (FILE=filename, exist=file_exists)
906 :
907 6 : IF (file_exists) THEN
908 2 : IF (para_env%is_source()) THEN
909 : CALL open_file(file_name=filename, &
910 : file_action="READ", &
911 : file_form="UNFORMATTED", &
912 : file_status="OLD", &
913 2 : unit_number=rst_unit)
914 :
915 2 : READ (rst_unit) qs_loc_env%localized_wfn_control%set_of_states
916 14 : READ (rst_unit) qs_loc_env%localized_wfn_control%lu_bound_states
917 6 : READ (rst_unit) qs_loc_env%localized_wfn_control%nloc_states
918 : END IF
919 : ELSE
920 4 : IF (output_unit > 0) THEN
921 : WRITE (output_unit, "(/,T10,A)") &
922 1 : "Restart file not available filename=<"//TRIM(filename)//'>'
923 : END IF
924 : END IF
925 6 : CALL para_env%bcast(file_exists)
926 :
927 6 : IF (file_exists) THEN
928 4 : restart_found = .TRUE.
929 :
930 4 : CALL para_env%bcast(qs_loc_env%localized_wfn_control%set_of_states)
931 4 : CALL para_env%bcast(qs_loc_env%localized_wfn_control%lu_bound_states)
932 4 : CALL para_env%bcast(qs_loc_env%localized_wfn_control%nloc_states)
933 :
934 12 : max_nloc = MAXVAL(qs_loc_env%localized_wfn_control%nloc_states(:))
935 :
936 12 : ALLOCATE (vecbuffer(1, nao))
937 4 : IF (ASSOCIATED(qs_loc_env%localized_wfn_control%loc_states)) THEN
938 2 : DEALLOCATE (qs_loc_env%localized_wfn_control%loc_states)
939 : END IF
940 12 : ALLOCATE (qs_loc_env%localized_wfn_control%loc_states(max_nloc, 2))
941 56 : qs_loc_env%localized_wfn_control%loc_states = 0
942 :
943 10 : DO ispin = 1, nspin
944 6 : IF (do_homo .OR. do_mixed) THEN
945 6 : nmo = mos(ispin)%nmo
946 : ELSE
947 0 : nmo = SIZE(evals(ispin)%array, 1)
948 : END IF
949 6 : IF (para_env%is_source() .AND. (nmo > 0)) THEN
950 3 : nloc = qs_loc_env%localized_wfn_control%nloc_states(ispin)
951 21 : READ (rst_unit) qs_loc_env%localized_wfn_control%loc_states(1:nloc, ispin)
952 3 : IF (do_homo .OR. do_mixed) THEN
953 3 : READ (rst_unit) nmo_read, homo_read, lfomo_read, nelectron_read
954 12 : ALLOCATE (eig_read(nmo_read), occ_read(nmo_read))
955 3 : eig_read = 0.0_dp
956 3 : occ_read = 0.0_dp
957 3 : READ (rst_unit) eig_read(1:nmo_read), occ_read(1:nmo_read)
958 : ELSE
959 0 : READ (rst_unit) nmo_read
960 0 : ALLOCATE (eig_read(nmo_read))
961 0 : eig_read = 0.0_dp
962 0 : READ (rst_unit) eig_read(1:nmo_read)
963 : END IF
964 3 : IF (nmo_read < nmo) THEN
965 : CALL cp_warn(__LOCATION__, &
966 : "The number of MOs on the restart unit is smaller than the number of "// &
967 0 : "the allocated MOs. ")
968 : END IF
969 3 : IF (nmo_read > nmo) THEN
970 : CALL cp_warn(__LOCATION__, &
971 : "The number of MOs on the restart unit is greater than the number of "// &
972 0 : "the allocated MOs. The read MO set will be truncated!")
973 : END IF
974 :
975 3 : nmo = MIN(nmo, nmo_read)
976 3 : IF (do_homo .OR. do_mixed) THEN
977 79 : mos(ispin)%eigenvalues(1:nmo) = eig_read(1:nmo)
978 79 : mos(ispin)%occupation_numbers(1:nmo) = occ_read(1:nmo)
979 3 : DEALLOCATE (eig_read, occ_read)
980 : ELSE
981 0 : evals(ispin)%array(1:nmo) = eig_read(1:nmo)
982 0 : DEALLOCATE (eig_read)
983 : END IF
984 :
985 : END IF
986 6 : IF (do_homo .OR. do_mixed) THEN
987 310 : CALL para_env%bcast(mos(ispin)%eigenvalues)
988 310 : CALL para_env%bcast(mos(ispin)%occupation_numbers)
989 : ELSE
990 0 : CALL para_env%bcast(evals(ispin)%array)
991 : END IF
992 :
993 162 : DO i = 1, nmo
994 152 : IF (para_env%is_source()) THEN
995 15236 : READ (rst_unit) vecbuffer
996 : ELSE
997 7656 : vecbuffer(1, :) = 0.0_dp
998 : END IF
999 60792 : CALL para_env%bcast(vecbuffer)
1000 : CALL cp_fm_set_submatrix(mos_localized(ispin), &
1001 158 : vecbuffer, 1, i, nao, 1, transpose=.TRUE.)
1002 : END DO
1003 : END DO
1004 :
1005 108 : CALL para_env%bcast(qs_loc_env%localized_wfn_control%loc_states)
1006 :
1007 4 : DEALLOCATE (vecbuffer)
1008 :
1009 : END IF
1010 :
1011 : ! Close restart file
1012 6 : IF (para_env%is_source()) THEN
1013 3 : IF (file_exists) CALL close_file(unit_number=rst_unit)
1014 : END IF
1015 :
1016 6 : CALL timestop(handle)
1017 :
1018 6 : END SUBROUTINE loc_read_restart
1019 :
1020 : ! **************************************************************************************************
1021 : !> \brief initializes everything needed for localization of the HOMOs
1022 : !> \param qs_loc_env ...
1023 : !> \param loc_section ...
1024 : !> \param do_homo ...
1025 : !> \param do_mixed ...
1026 : !> \param do_xas ...
1027 : !> \param nloc_xas ...
1028 : !> \param spin_xas ...
1029 : !> \par History
1030 : !> 2009 created
1031 : ! **************************************************************************************************
1032 432 : SUBROUTINE qs_loc_control_init(qs_loc_env, loc_section, do_homo, do_mixed, &
1033 : do_xas, nloc_xas, spin_xas)
1034 :
1035 : TYPE(qs_loc_env_type), POINTER :: qs_loc_env
1036 : TYPE(section_vals_type), POINTER :: loc_section
1037 : LOGICAL, INTENT(IN) :: do_homo
1038 : LOGICAL, INTENT(IN), OPTIONAL :: do_mixed, do_xas
1039 : INTEGER, INTENT(IN), OPTIONAL :: nloc_xas, spin_xas
1040 :
1041 : LOGICAL :: my_do_mixed
1042 : TYPE(localized_wfn_control_type), POINTER :: localized_wfn_control
1043 :
1044 432 : NULLIFY (localized_wfn_control)
1045 :
1046 432 : IF (PRESENT(do_mixed)) THEN
1047 2 : my_do_mixed = do_mixed
1048 : ELSE
1049 430 : my_do_mixed = .FALSE.
1050 : END IF
1051 432 : CALL localized_wfn_control_create(localized_wfn_control)
1052 432 : CALL set_qs_loc_env(qs_loc_env, localized_wfn_control=localized_wfn_control)
1053 432 : CALL localized_wfn_control_release(localized_wfn_control)
1054 432 : CALL get_qs_loc_env(qs_loc_env, localized_wfn_control=localized_wfn_control)
1055 432 : localized_wfn_control%do_homo = do_homo
1056 432 : localized_wfn_control%do_mixed = my_do_mixed
1057 : CALL read_loc_section(localized_wfn_control, loc_section, qs_loc_env%do_localize, &
1058 432 : my_do_mixed, do_xas, nloc_xas, spin_xas)
1059 :
1060 432 : END SUBROUTINE qs_loc_control_init
1061 :
1062 : ! **************************************************************************************************
1063 : !> \brief initializes everything needed for localization of the molecular orbitals
1064 : !> \param qs_env ...
1065 : !> \param qs_loc_env ...
1066 : !> \param localize_section ...
1067 : !> \param mos_localized ...
1068 : !> \param do_homo ...
1069 : !> \param do_mo_cubes ...
1070 : !> \param mo_loc_history ...
1071 : !> \param evals ...
1072 : !> \param tot_zeff_corr ...
1073 : !> \param do_mixed ...
1074 : ! **************************************************************************************************
1075 324 : SUBROUTINE qs_loc_init(qs_env, qs_loc_env, localize_section, mos_localized, &
1076 : do_homo, do_mo_cubes, mo_loc_history, evals, &
1077 : tot_zeff_corr, do_mixed)
1078 :
1079 : TYPE(qs_environment_type), POINTER :: qs_env
1080 : TYPE(qs_loc_env_type), POINTER :: qs_loc_env
1081 : TYPE(section_vals_type), POINTER :: localize_section
1082 : TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT) :: mos_localized
1083 : LOGICAL, OPTIONAL :: do_homo, do_mo_cubes
1084 : TYPE(cp_fm_type), DIMENSION(:), OPTIONAL, POINTER :: mo_loc_history
1085 : TYPE(cp_1d_r_p_type), DIMENSION(:), OPTIONAL, &
1086 : POINTER :: evals
1087 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: tot_zeff_corr
1088 : LOGICAL, OPTIONAL :: do_mixed
1089 :
1090 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_loc_init'
1091 :
1092 : INTEGER :: handle, homo, i, ilast_intocc, ilow, ispin, iup, n_mo(2), n_mos(2), nao, &
1093 : nelectron, nextra, nmoloc(2), nocc, npocc, nspin, output_unit
1094 : LOGICAL :: my_do_homo, my_do_mixed, my_do_mo_cubes, &
1095 : restart_found
1096 : REAL(KIND=dp) :: maxocc, my_tot_zeff_corr
1097 324 : REAL(KIND=dp), DIMENSION(:), POINTER :: mo_eigenvalues, occupation
1098 : TYPE(cp_fm_type), POINTER :: mo_coeff
1099 : TYPE(cp_logger_type), POINTER :: logger
1100 324 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_rmpv, mo_derivs
1101 : TYPE(dft_control_type), POINTER :: dft_control
1102 : TYPE(localized_wfn_control_type), POINTER :: localized_wfn_control
1103 324 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1104 : TYPE(mp_para_env_type), POINTER :: para_env
1105 : TYPE(scf_control_type), POINTER :: scf_control
1106 : TYPE(section_vals_type), POINTER :: loc_print_section
1107 :
1108 324 : CALL timeset(routineN, handle)
1109 :
1110 324 : NULLIFY (mos, mo_coeff, mo_eigenvalues, occupation, ks_rmpv, mo_derivs, scf_control, para_env)
1111 : CALL get_qs_env(qs_env, &
1112 : mos=mos, &
1113 : matrix_ks=ks_rmpv, &
1114 : mo_derivs=mo_derivs, &
1115 : scf_control=scf_control, &
1116 : dft_control=dft_control, &
1117 324 : para_env=para_env)
1118 :
1119 324 : loc_print_section => section_vals_get_subs_vals(localize_section, "PRINT")
1120 :
1121 324 : logger => cp_get_default_logger()
1122 324 : output_unit = cp_logger_get_default_io_unit(logger)
1123 :
1124 324 : nspin = SIZE(mos)
1125 324 : IF (PRESENT(do_homo)) THEN
1126 324 : my_do_homo = do_homo
1127 : ELSE
1128 0 : my_do_homo = .TRUE.
1129 : END IF
1130 324 : IF (PRESENT(do_mo_cubes)) THEN
1131 134 : my_do_mo_cubes = do_mo_cubes
1132 : ELSE
1133 : my_do_mo_cubes = .FALSE.
1134 : END IF
1135 324 : IF (PRESENT(do_mixed)) THEN
1136 2 : my_do_mixed = do_mixed
1137 : ELSE
1138 322 : my_do_mixed = .FALSE.
1139 : END IF
1140 324 : IF (PRESENT(tot_zeff_corr)) THEN
1141 2 : my_tot_zeff_corr = tot_zeff_corr
1142 : ELSE
1143 322 : my_tot_zeff_corr = 0.0_dp
1144 : END IF
1145 324 : restart_found = .FALSE.
1146 :
1147 324 : IF (qs_loc_env%do_localize) THEN
1148 : ! Some setup for MOs to be localized
1149 308 : CALL get_qs_loc_env(qs_loc_env, localized_wfn_control=localized_wfn_control)
1150 308 : IF (localized_wfn_control%loc_restart) THEN
1151 6 : IF (localized_wfn_control%nextra > 0) THEN
1152 : ! currently only the occupied guess is read
1153 0 : my_do_homo = .FALSE.
1154 : END IF
1155 : CALL loc_read_restart(qs_loc_env, mos, mos_localized, localize_section, &
1156 : loc_print_section, para_env, my_do_homo, restart_found, evals=evals, &
1157 6 : do_mixed=my_do_mixed)
1158 9 : IF (output_unit > 0) WRITE (output_unit, "(/,T2,A,A)") "LOCALIZATION| ", &
1159 6 : " The orbitals to be localized are read from localization restart file."
1160 18 : nmoloc = localized_wfn_control%nloc_states
1161 18 : localized_wfn_control%nguess = nmoloc
1162 6 : IF (localized_wfn_control%nextra > 0) THEN
1163 : ! reset different variables in localized_wfn_control:
1164 : ! lu_bound_states, nloc_states, loc_states
1165 0 : localized_wfn_control%loc_restart = restart_found
1166 0 : localized_wfn_control%set_of_states = state_loc_mixed
1167 0 : DO ispin = 1, nspin
1168 : CALL get_mo_set(mos(ispin), homo=homo, occupation_numbers=occupation, &
1169 0 : maxocc=maxocc)
1170 0 : nextra = localized_wfn_control%nextra
1171 0 : nocc = homo
1172 0 : DO i = nocc, 1, -1
1173 0 : IF (maxocc - occupation(i) < localized_wfn_control%eps_occ) THEN
1174 0 : ilast_intocc = i
1175 0 : EXIT
1176 : END IF
1177 : END DO
1178 0 : nocc = ilast_intocc
1179 0 : npocc = homo - nocc
1180 0 : nmoloc(ispin) = nocc + nextra
1181 0 : localized_wfn_control%lu_bound_states(1, ispin) = 1
1182 0 : localized_wfn_control%lu_bound_states(2, ispin) = nmoloc(ispin)
1183 0 : localized_wfn_control%nloc_states(ispin) = nmoloc(ispin)
1184 : END DO
1185 0 : my_do_homo = .FALSE.
1186 : END IF
1187 : END IF
1188 308 : IF (.NOT. restart_found) THEN
1189 304 : nmoloc = 0
1190 726 : DO ispin = 1, nspin
1191 : CALL get_mo_set(mos(ispin), nmo=n_mo(ispin), nelectron=nelectron, homo=homo, nao=nao, &
1192 : mo_coeff=mo_coeff, eigenvalues=mo_eigenvalues, occupation_numbers=occupation, &
1193 422 : maxocc=maxocc)
1194 : ! Get eigenstates (only needed if not already calculated before)
1195 : IF ((.NOT. my_do_mo_cubes) &
1196 : .AND. my_do_homo .AND. ASSOCIATED(qs_env%scf_env) &
1197 422 : .AND. qs_env%scf_env%method == ot_method_nr .AND. (.NOT. dft_control%restricted)) THEN
1198 30 : CALL make_mo_eig(mos, nspin, ks_rmpv, scf_control, mo_derivs)
1199 : END IF
1200 1148 : IF (localized_wfn_control%set_of_states == state_loc_all .AND. my_do_homo) THEN
1201 382 : nmoloc(ispin) = NINT(nelectron/occupation(1))
1202 382 : IF (n_mo(ispin) > homo) THEN
1203 28 : DO i = nmoloc(ispin), 1, -1
1204 28 : IF (occupation(1) - occupation(i) < localized_wfn_control%eps_occ) THEN
1205 14 : ilast_intocc = i
1206 14 : EXIT
1207 : END IF
1208 : END DO
1209 : ELSE
1210 368 : ilast_intocc = nmoloc(ispin)
1211 : END IF
1212 382 : nmoloc(ispin) = ilast_intocc
1213 382 : localized_wfn_control%lu_bound_states(1, ispin) = 1
1214 382 : localized_wfn_control%lu_bound_states(2, ispin) = ilast_intocc
1215 382 : IF (nmoloc(ispin) /= n_mo(ispin)) THEN
1216 14 : IF (output_unit > 0) THEN
1217 : WRITE (output_unit, "(/,T2,A,I4,A,I6,A,/,T15,A,F12.6,A,F12.6,A)") &
1218 7 : "LOCALIZATION| Spin ", ispin, " The first ", &
1219 7 : ilast_intocc, " occupied orbitals are localized,", " with energies from ", &
1220 14 : mo_eigenvalues(1), " to ", mo_eigenvalues(ilast_intocc), " [a.u.]."
1221 : END IF
1222 : END IF
1223 40 : ELSE IF (localized_wfn_control%set_of_states == energy_loc_range .AND. my_do_homo) THEN
1224 12 : ilow = 0
1225 12 : iup = 0
1226 20 : DO i = 1, n_mo(ispin)
1227 20 : IF (mo_eigenvalues(i) >= localized_wfn_control%lu_ene_bound(1)) THEN
1228 : ilow = i
1229 : EXIT
1230 : END IF
1231 : END DO
1232 306 : DO i = n_mo(ispin), 1, -1
1233 306 : IF (mo_eigenvalues(i) <= localized_wfn_control%lu_ene_bound(2)) THEN
1234 : iup = i
1235 : EXIT
1236 : END IF
1237 : END DO
1238 12 : localized_wfn_control%lu_bound_states(1, ispin) = ilow
1239 12 : localized_wfn_control%lu_bound_states(2, ispin) = iup
1240 12 : localized_wfn_control%nloc_states(ispin) = iup - ilow + 1
1241 12 : nmoloc(ispin) = localized_wfn_control%nloc_states(ispin)
1242 12 : IF (occupation(ilow) - occupation(iup) > localized_wfn_control%eps_occ) THEN
1243 : CALL cp_abort(__LOCATION__, &
1244 : "The selected energy range includes orbitals with different occupation number. "// &
1245 0 : "The localization procedure cannot be applied.")
1246 : END IF
1247 18 : IF (output_unit > 0) WRITE (output_unit, "(/,T2,A,I4,A,I6,A)") "LOCALIZATION| Spin ", ispin, " : ", &
1248 12 : nmoloc(ispin), " orbitals in the selected energy range are localized."
1249 28 : ELSE IF (localized_wfn_control%set_of_states == state_loc_all .AND. (.NOT. my_do_homo)) THEN
1250 0 : nmoloc(ispin) = n_mo(ispin) - homo
1251 0 : localized_wfn_control%lu_bound_states(1, ispin) = homo + 1
1252 0 : localized_wfn_control%lu_bound_states(2, ispin) = n_mo(ispin)
1253 0 : IF (output_unit > 0) THEN
1254 : WRITE (output_unit, "(/,T2,A,I4,A,I6,A,/,T15,A,F12.6,A,F12.6,A)") &
1255 0 : "LOCALIZATION| Spin ", ispin, " The first ", &
1256 0 : nmoloc(ispin), " virtual orbitals are localized,", " with energies from ", &
1257 0 : mo_eigenvalues(homo + 1), " to ", mo_eigenvalues(n_mo(ispin)), " [a.u.]."
1258 : END IF
1259 28 : ELSE IF (localized_wfn_control%set_of_states == state_loc_mixed) THEN
1260 2 : nextra = localized_wfn_control%nextra
1261 2 : nocc = homo
1262 6 : DO i = nocc, 1, -1
1263 6 : IF (maxocc - occupation(i) < localized_wfn_control%eps_occ) THEN
1264 2 : ilast_intocc = i
1265 2 : EXIT
1266 : END IF
1267 : END DO
1268 2 : nocc = ilast_intocc
1269 2 : npocc = homo - nocc
1270 2 : nmoloc(ispin) = nocc + nextra
1271 2 : localized_wfn_control%lu_bound_states(1, ispin) = 1
1272 2 : localized_wfn_control%lu_bound_states(2, ispin) = nmoloc(ispin)
1273 2 : IF (output_unit > 0) THEN
1274 : WRITE (output_unit, "(/,T2,A,I4,A,I6,A,/,T15,A,I6,/,T15,A,I6,/,T15,A,I6,/,T15,A,F12.6,A)") &
1275 1 : "LOCALIZATION| Spin ", ispin, " The first ", &
1276 1 : nmoloc(ispin), " orbitals are localized.", &
1277 1 : "Number of fully occupied MOs: ", nocc, &
1278 1 : "Number of partially occupied MOs: ", npocc, &
1279 1 : "Number of extra degrees of freedom: ", nextra, &
1280 2 : "Excess charge: ", my_tot_zeff_corr, " electrons"
1281 : END IF
1282 : ELSE
1283 26 : nmoloc(ispin) = MIN(localized_wfn_control%nloc_states(1), n_mo(ispin))
1284 33 : IF (output_unit > 0 .AND. my_do_homo) WRITE (output_unit, "(/,T2,A,I4,A,I6,A)") "LOCALIZATION| Spin ", ispin, &
1285 14 : " : ", nmoloc(ispin), " occupied orbitals are localized, as given in the input list."
1286 19 : IF (output_unit > 0 .AND. (.NOT. my_do_homo)) WRITE (output_unit, "(/,T2,A,I4,A,I6,A)") "LOCALIZATION| Spin ", &
1287 12 : ispin, " : ", nmoloc(ispin), " unoccupied orbitals are localized, as given in the input list."
1288 26 : IF (n_mo(ispin) > homo .AND. my_do_homo) THEN
1289 8 : ilow = localized_wfn_control%loc_states(1, ispin)
1290 56 : DO i = 2, nmoloc(ispin)
1291 48 : iup = localized_wfn_control%loc_states(i, ispin)
1292 56 : IF (ABS(occupation(ilow) - occupation(iup)) > localized_wfn_control%eps_occ) THEN
1293 : ! write warning
1294 : CALL cp_warn(__LOCATION__, &
1295 : "User requested the calculation of localized wavefunction from a subset of MOs, "// &
1296 : "including MOs with different occupations. Check the selected subset, "// &
1297 : "the electronic density is not invariant with "// &
1298 0 : "respect to rotations among orbitals with different occupation numbers!")
1299 : END IF
1300 : END DO
1301 : END IF
1302 : END IF
1303 : END DO ! ispin
1304 912 : n_mos(:) = nao - n_mo(:)
1305 304 : IF (my_do_homo .OR. my_do_mixed) n_mos = n_mo
1306 304 : CALL set_loc_wfn_lists(localized_wfn_control, nmoloc, n_mos, nspin)
1307 : END IF
1308 308 : CALL set_loc_centers(localized_wfn_control, nmoloc, nspin)
1309 308 : IF (my_do_homo .OR. my_do_mixed) THEN
1310 : CALL qs_loc_env_init(qs_loc_env, localized_wfn_control, qs_env, &
1311 302 : loc_coeff=mos_localized, mo_loc_history=mo_loc_history)
1312 : END IF
1313 : ELSE
1314 : ! Let's inform in case the section is not present in the input
1315 : CALL cp_warn(__LOCATION__, &
1316 : "User requested the calculation of the localized wavefunction but the section "// &
1317 16 : "LOCALIZE was not specified. Localization will not be performed!")
1318 : END IF
1319 :
1320 324 : CALL timestop(handle)
1321 :
1322 324 : END SUBROUTINE qs_loc_init
1323 :
1324 : ! **************************************************************************************************
1325 : !> \brief read the controlparameter from input, using the new input scheme
1326 : !> \param localized_wfn_control ...
1327 : !> \param loc_section ...
1328 : !> \param localize ...
1329 : !> \param do_mixed ...
1330 : !> \param do_xas ...
1331 : !> \param nloc_xas ...
1332 : !> \param spin_channel_xas ...
1333 : !> \par History
1334 : !> 05.2005 created [MI]
1335 : ! **************************************************************************************************
1336 864 : SUBROUTINE read_loc_section(localized_wfn_control, loc_section, &
1337 : localize, do_mixed, do_xas, nloc_xas, spin_channel_xas)
1338 :
1339 : TYPE(localized_wfn_control_type), POINTER :: localized_wfn_control
1340 : TYPE(section_vals_type), POINTER :: loc_section
1341 : LOGICAL, INTENT(OUT) :: localize
1342 : LOGICAL, INTENT(IN), OPTIONAL :: do_mixed, do_xas
1343 : INTEGER, INTENT(IN), OPTIONAL :: nloc_xas, spin_channel_xas
1344 :
1345 : INTEGER :: i, ind, ir, n_list, n_rep, n_state, &
1346 : nextra, nline, other_spin, &
1347 : output_unit, spin_xas
1348 432 : INTEGER, DIMENSION(:), POINTER :: list, loc_list
1349 : LOGICAL :: my_do_mixed, my_do_xas
1350 432 : REAL(dp), POINTER :: ene(:)
1351 : TYPE(cp_logger_type), POINTER :: logger
1352 : TYPE(section_vals_type), POINTER :: loc_print_section
1353 :
1354 432 : my_do_xas = .FALSE.
1355 432 : spin_xas = 1
1356 432 : IF (PRESENT(do_xas)) THEN
1357 108 : my_do_xas = do_xas
1358 108 : CPASSERT(PRESENT(nloc_xas))
1359 : END IF
1360 432 : IF (PRESENT(spin_channel_xas)) spin_xas = spin_channel_xas
1361 432 : my_do_mixed = .FALSE.
1362 432 : IF (PRESENT(do_mixed)) THEN
1363 432 : my_do_mixed = do_mixed
1364 : END IF
1365 432 : CPASSERT(ASSOCIATED(loc_section))
1366 432 : NULLIFY (logger)
1367 432 : logger => cp_get_default_logger()
1368 :
1369 432 : CALL section_vals_val_get(loc_section, "_SECTION_PARAMETERS_", l_val=localize)
1370 432 : IF (localize) THEN
1371 384 : loc_print_section => section_vals_get_subs_vals(loc_section, "PRINT")
1372 384 : NULLIFY (list)
1373 384 : NULLIFY (loc_list)
1374 2688 : localized_wfn_control%lu_bound_states = 0
1375 1152 : localized_wfn_control%lu_ene_bound = 0.0_dp
1376 1152 : localized_wfn_control%nloc_states = 0
1377 384 : localized_wfn_control%set_of_states = 0
1378 384 : localized_wfn_control%nextra = 0
1379 384 : n_state = 0
1380 :
1381 : CALL section_vals_val_get(loc_section, "MAX_ITER", &
1382 384 : i_val=localized_wfn_control%max_iter)
1383 : CALL section_vals_val_get(loc_section, "MAX_CRAZY_ANGLE", &
1384 384 : r_val=localized_wfn_control%max_crazy_angle)
1385 : CALL section_vals_val_get(loc_section, "CRAZY_SCALE", &
1386 384 : r_val=localized_wfn_control%crazy_scale)
1387 : CALL section_vals_val_get(loc_section, "EPS_OCCUPATION", &
1388 384 : r_val=localized_wfn_control%eps_occ)
1389 : CALL section_vals_val_get(loc_section, "CRAZY_USE_DIAG", &
1390 384 : l_val=localized_wfn_control%crazy_use_diag)
1391 : CALL section_vals_val_get(loc_section, "OUT_ITER_EACH", &
1392 384 : i_val=localized_wfn_control%out_each)
1393 : CALL section_vals_val_get(loc_section, "EPS_LOCALIZATION", &
1394 384 : r_val=localized_wfn_control%eps_localization)
1395 : CALL section_vals_val_get(loc_section, "MIN_OR_MAX", &
1396 384 : i_val=localized_wfn_control%min_or_max)
1397 : CALL section_vals_val_get(loc_section, "JACOBI_FALLBACK", &
1398 384 : l_val=localized_wfn_control%jacobi_fallback)
1399 : CALL section_vals_val_get(loc_section, "JACOBI_REFINEMENT", &
1400 384 : l_val=localized_wfn_control%jacobi_refinement)
1401 : CALL section_vals_val_get(loc_section, "METHOD", &
1402 384 : i_val=localized_wfn_control%localization_method)
1403 : CALL section_vals_val_get(loc_section, "OPERATOR", &
1404 384 : i_val=localized_wfn_control%operator_type)
1405 : CALL section_vals_val_get(loc_section, "RESTART", &
1406 384 : l_val=localized_wfn_control%loc_restart)
1407 : CALL section_vals_val_get(loc_section, "USE_HISTORY", &
1408 384 : l_val=localized_wfn_control%use_history)
1409 : CALL section_vals_val_get(loc_section, "NEXTRA", &
1410 384 : i_val=localized_wfn_control%nextra)
1411 : CALL section_vals_val_get(loc_section, "CPO_GUESS", &
1412 384 : i_val=localized_wfn_control%coeff_po_guess)
1413 : CALL section_vals_val_get(loc_section, "CPO_GUESS_SPACE", &
1414 384 : i_val=localized_wfn_control%coeff_po_guess_mo_space)
1415 : CALL section_vals_val_get(loc_section, "CG_PO", &
1416 384 : l_val=localized_wfn_control%do_cg_po)
1417 :
1418 384 : IF (localized_wfn_control%do_homo) THEN
1419 : ! List of States HOMO
1420 376 : CALL section_vals_val_get(loc_section, "LIST", n_rep_val=n_rep)
1421 376 : IF (n_rep > 0) THEN
1422 14 : n_list = 0
1423 28 : DO ir = 1, n_rep
1424 14 : NULLIFY (list)
1425 14 : CALL section_vals_val_get(loc_section, "LIST", i_rep_val=ir, i_vals=list)
1426 28 : IF (ASSOCIATED(list)) THEN
1427 14 : CALL reallocate(loc_list, 1, n_list + SIZE(list))
1428 90 : DO i = 1, SIZE(list)
1429 90 : loc_list(n_list + i) = list(i)
1430 : END DO ! i
1431 14 : n_list = n_list + SIZE(list)
1432 : END IF
1433 : END DO ! ir
1434 14 : IF (n_list /= 0) THEN
1435 14 : localized_wfn_control%set_of_states = state_loc_list
1436 42 : ALLOCATE (localized_wfn_control%loc_states(n_list, 2))
1437 194 : localized_wfn_control%loc_states = 0
1438 180 : localized_wfn_control%loc_states(:, 1) = loc_list(:)
1439 180 : localized_wfn_control%loc_states(:, 2) = loc_list(:)
1440 14 : localized_wfn_control%nloc_states(1) = n_list
1441 14 : localized_wfn_control%nloc_states(2) = n_list
1442 14 : IF (my_do_xas) THEN
1443 4 : other_spin = 2
1444 4 : IF (spin_xas == 2) other_spin = 1
1445 4 : localized_wfn_control%nloc_states(other_spin) = 0
1446 22 : localized_wfn_control%loc_states(:, other_spin) = 0
1447 : END IF
1448 14 : DEALLOCATE (loc_list)
1449 : END IF
1450 : END IF
1451 :
1452 : ELSE
1453 : ! List of States LUMO
1454 8 : CALL section_vals_val_get(loc_section, "LIST_UNOCCUPIED", n_rep_val=n_rep)
1455 8 : IF (n_rep > 0) THEN
1456 6 : n_list = 0
1457 12 : DO ir = 1, n_rep
1458 6 : NULLIFY (list)
1459 6 : CALL section_vals_val_get(loc_section, "LIST_UNOCCUPIED", i_rep_val=ir, i_vals=list)
1460 12 : IF (ASSOCIATED(list)) THEN
1461 6 : CALL reallocate(loc_list, 1, n_list + SIZE(list))
1462 46 : DO i = 1, SIZE(list)
1463 46 : loc_list(n_list + i) = list(i)
1464 : END DO ! i
1465 6 : n_list = n_list + SIZE(list)
1466 : END IF
1467 : END DO ! ir
1468 6 : IF (n_list /= 0) THEN
1469 6 : localized_wfn_control%set_of_states = state_loc_list
1470 18 : ALLOCATE (localized_wfn_control%loc_states(n_list, 2))
1471 98 : localized_wfn_control%loc_states = 0
1472 92 : localized_wfn_control%loc_states(:, 1) = loc_list(:)
1473 92 : localized_wfn_control%loc_states(:, 2) = loc_list(:)
1474 6 : localized_wfn_control%nloc_states(1) = n_list
1475 6 : DEALLOCATE (loc_list)
1476 : END IF
1477 : END IF
1478 : END IF
1479 :
1480 384 : IF (localized_wfn_control%set_of_states == 0) THEN
1481 364 : CALL section_vals_val_get(loc_section, "ENERGY_RANGE", r_vals=ene)
1482 364 : IF (ene(1) /= ene(2)) THEN
1483 10 : localized_wfn_control%set_of_states = energy_loc_range
1484 10 : localized_wfn_control%lu_ene_bound(1) = ene(1)
1485 10 : localized_wfn_control%lu_ene_bound(2) = ene(2)
1486 : END IF
1487 : END IF
1488 :
1489 : ! All States or XAS specific states
1490 384 : IF (localized_wfn_control%set_of_states == 0) THEN
1491 354 : IF (my_do_xas) THEN
1492 72 : localized_wfn_control%set_of_states = state_loc_range
1493 216 : localized_wfn_control%nloc_states(:) = 0
1494 216 : localized_wfn_control%lu_bound_states(1, :) = 0
1495 216 : localized_wfn_control%lu_bound_states(2, :) = 0
1496 72 : localized_wfn_control%nloc_states(spin_xas) = nloc_xas
1497 72 : localized_wfn_control%lu_bound_states(1, spin_xas) = 1
1498 72 : localized_wfn_control%lu_bound_states(2, spin_xas) = nloc_xas
1499 282 : ELSE IF (my_do_mixed) THEN
1500 2 : localized_wfn_control%set_of_states = state_loc_mixed
1501 2 : nextra = localized_wfn_control%nextra
1502 : ELSE
1503 280 : localized_wfn_control%set_of_states = state_loc_all
1504 : END IF
1505 : END IF
1506 :
1507 : localized_wfn_control%print_centers = &
1508 : BTEST(cp_print_key_should_output(logger%iter_info, loc_print_section, &
1509 384 : "WANNIER_CENTERS"), cp_p_file)
1510 : localized_wfn_control%print_spreads = &
1511 : BTEST(cp_print_key_should_output(logger%iter_info, loc_print_section, &
1512 384 : "WANNIER_SPREADS"), cp_p_file)
1513 : localized_wfn_control%print_cubes = &
1514 : BTEST(cp_print_key_should_output(logger%iter_info, loc_print_section, &
1515 384 : "WANNIER_CUBES"), cp_p_file)
1516 :
1517 : output_unit = cp_print_key_unit_nr(logger, loc_print_section, "PROGRAM_RUN_INFO", &
1518 384 : extension=".Log")
1519 :
1520 384 : IF (output_unit > 0) THEN
1521 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
1522 192 : "LOCALIZE| The spread relative to a set of orbitals is computed"
1523 :
1524 332 : SELECT CASE (localized_wfn_control%set_of_states)
1525 : CASE (state_loc_all)
1526 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
1527 140 : "LOCALIZE| Orbitals to be localized: All orbitals"
1528 : WRITE (UNIT=output_unit, FMT="(T2,A,/,T12,A,F16.8)") &
1529 140 : "LOCALIZE| If fractional occupation, fully occupied MOs are those ", &
1530 280 : "within occupation tolerance of ", localized_wfn_control%eps_occ
1531 : CASE (state_loc_range)
1532 : WRITE (UNIT=output_unit, FMT="(T2,A,T65,I8,A,I8)") &
1533 36 : "LOCALIZE| Orbitals to be localized: Those with index between ", &
1534 36 : localized_wfn_control%lu_bound_states(1, spin_xas), " and ", &
1535 72 : localized_wfn_control%lu_bound_states(2, spin_xas)
1536 : CASE (state_loc_list)
1537 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
1538 10 : "LOCALIZE| Orbitals to be localized: Those with index in the following list"
1539 10 : nline = localized_wfn_control%nloc_states(1)/10 + 1
1540 10 : ind = 0
1541 21 : DO i = 1, nline
1542 21 : IF (ind + 10 < localized_wfn_control%nloc_states(1)) THEN
1543 11 : WRITE (UNIT=output_unit, FMT="(T8,10I7)") localized_wfn_control%loc_states(ind + 1:ind + 10, 1)
1544 1 : ind = ind + 10
1545 : ELSE
1546 : WRITE (UNIT=output_unit, FMT="(T8,10I7)") &
1547 58 : localized_wfn_control%loc_states(ind + 1:localized_wfn_control%nloc_states(1), 1)
1548 10 : ind = localized_wfn_control%nloc_states(1)
1549 : END IF
1550 : END DO
1551 : CASE (energy_loc_range)
1552 : WRITE (UNIT=output_unit, FMT="(T2,A,T65,/,f16.6,A,f16.6,A)") &
1553 5 : "LOCALIZE| Orbitals to be localized: Those with energy in the range between ", &
1554 10 : localized_wfn_control%lu_ene_bound(1), " and ", localized_wfn_control%lu_ene_bound(2), " a.u."
1555 : CASE (state_loc_mixed)
1556 : WRITE (UNIT=output_unit, FMT="(T2,A,I4,A)") &
1557 1 : "LOCALIZE| Orbitals to be localized: Occupied orbitals + ", nextra, " orbitals"
1558 : CASE DEFAULT
1559 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
1560 192 : "LOCALIZE| Orbitals to be localized: None "
1561 : END SELECT
1562 :
1563 381 : SELECT CASE (localized_wfn_control%operator_type)
1564 : CASE (op_loc_berry)
1565 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
1566 189 : "LOCALIZE| Spread defined by the Berry phase operator "
1567 : CASE (op_loc_boys)
1568 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
1569 0 : "LOCALIZE| Spread defined by the Boys phase operator "
1570 : CASE DEFAULT
1571 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
1572 192 : "LOCALIZE| Spread defined by the Pipek phase operator "
1573 : END SELECT
1574 :
1575 328 : SELECT CASE (localized_wfn_control%localization_method)
1576 : CASE (do_loc_jacobi)
1577 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
1578 136 : "LOCALIZE| Optimal unitary transformation generated by Jacobi algorithm"
1579 : CASE (do_loc_crazy)
1580 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
1581 38 : "LOCALIZE| Optimal unitary transformation generated by Crazy angle algorithm"
1582 : WRITE (UNIT=output_unit, FMT="(T2,A,F16.8)") &
1583 38 : "LOCALIZE| maximum angle: ", localized_wfn_control%max_crazy_angle
1584 : WRITE (UNIT=output_unit, FMT="(T2,A,F16.8)") &
1585 38 : "LOCALIZE| scaling: ", localized_wfn_control%crazy_scale
1586 : WRITE (UNIT=output_unit, FMT="(T2,A,L1)") &
1587 38 : "LOCALIZE| use diag:", localized_wfn_control%crazy_use_diag
1588 : CASE (do_loc_gapo)
1589 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
1590 1 : "LOCALIZE| Optimal unitary transformation generated by gradient ascent algorithm "
1591 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
1592 1 : "LOCALIZE| for partially occupied wannier functions"
1593 : CASE (do_loc_direct)
1594 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
1595 1 : "LOCALIZE| Optimal unitary transformation generated by direct algorithm"
1596 : CASE (do_loc_l1_norm_sd)
1597 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
1598 9 : "LOCALIZE| Optimal unitary transformation generated by "
1599 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
1600 9 : "LOCALIZE| steepest descent algorithm applied on an approximate l1 norm"
1601 : CASE (do_loc_none)
1602 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
1603 0 : "LOCALIZE| No unitary transformation is applied"
1604 : CASE (do_loc_scdm)
1605 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
1606 192 : "LOCALIZE| Pivoted QR decomposition is used to transform coefficients"
1607 : END SELECT
1608 :
1609 : END IF ! process has output_unit
1610 :
1611 384 : CALL cp_print_key_finished_output(output_unit, logger, loc_print_section, "PROGRAM_RUN_INFO")
1612 :
1613 : ELSE
1614 48 : localized_wfn_control%localization_method = do_loc_none
1615 48 : localized_wfn_control%localization_method = state_loc_none
1616 48 : localized_wfn_control%print_centers = .FALSE.
1617 48 : localized_wfn_control%print_spreads = .FALSE.
1618 48 : localized_wfn_control%print_cubes = .FALSE.
1619 : END IF
1620 :
1621 432 : END SUBROUTINE read_loc_section
1622 :
1623 : ! **************************************************************************************************
1624 : !> \brief create the center and spread array and the file names for the output
1625 : !> \param localized_wfn_control ...
1626 : !> \param nmoloc ...
1627 : !> \param nspins ...
1628 : !> \par History
1629 : !> 04.2005 created [MI]
1630 : ! **************************************************************************************************
1631 484 : SUBROUTINE set_loc_centers(localized_wfn_control, nmoloc, nspins)
1632 :
1633 : TYPE(localized_wfn_control_type) :: localized_wfn_control
1634 : INTEGER, DIMENSION(2), INTENT(IN) :: nmoloc
1635 : INTEGER, INTENT(IN) :: nspins
1636 :
1637 : INTEGER :: ispin
1638 :
1639 1144 : DO ispin = 1, nspins
1640 1938 : ALLOCATE (localized_wfn_control%centers_set(ispin)%array(6, nmoloc(ispin)))
1641 32182 : localized_wfn_control%centers_set(ispin)%array = 0.0_dp
1642 : END DO
1643 :
1644 484 : END SUBROUTINE set_loc_centers
1645 :
1646 : ! **************************************************************************************************
1647 : !> \brief create the lists of mos that are taken into account
1648 : !> \param localized_wfn_control ...
1649 : !> \param nmoloc ...
1650 : !> \param nmo ...
1651 : !> \param nspins ...
1652 : !> \param my_spin ...
1653 : !> \par History
1654 : !> 04.2005 created [MI]
1655 : ! **************************************************************************************************
1656 346 : SUBROUTINE set_loc_wfn_lists(localized_wfn_control, nmoloc, nmo, nspins, my_spin)
1657 :
1658 : TYPE(localized_wfn_control_type) :: localized_wfn_control
1659 : INTEGER, DIMENSION(2), INTENT(IN) :: nmoloc, nmo
1660 : INTEGER, INTENT(IN) :: nspins
1661 : INTEGER, INTENT(IN), OPTIONAL :: my_spin
1662 :
1663 : CHARACTER(len=*), PARAMETER :: routineN = 'set_loc_wfn_lists'
1664 :
1665 : INTEGER :: i, ispin, max_iloc, max_nmoloc, state
1666 :
1667 346 : CALL timeset(routineN, state)
1668 :
1669 1038 : localized_wfn_control%nloc_states(1:2) = nmoloc(1:2)
1670 346 : max_nmoloc = MAX(nmoloc(1), nmoloc(2))
1671 :
1672 364 : SELECT CASE (localized_wfn_control%set_of_states)
1673 : CASE (state_loc_list)
1674 : ! List
1675 18 : CPASSERT(ASSOCIATED(localized_wfn_control%loc_states))
1676 52 : DO ispin = 1, nspins
1677 34 : localized_wfn_control%lu_bound_states(1, ispin) = 1
1678 34 : localized_wfn_control%lu_bound_states(2, ispin) = nmoloc(ispin)
1679 52 : IF (nmoloc(ispin) < 1) THEN
1680 4 : localized_wfn_control%lu_bound_states(1, ispin) = 0
1681 22 : localized_wfn_control%loc_states(:, ispin) = 0
1682 : END IF
1683 : END DO
1684 : CASE (state_loc_range)
1685 : ! Range
1686 114 : ALLOCATE (localized_wfn_control%loc_states(max_nmoloc, 2))
1687 446 : localized_wfn_control%loc_states = 0
1688 114 : DO ispin = 1, nspins
1689 : localized_wfn_control%lu_bound_states(1, ispin) = &
1690 76 : localized_wfn_control%lu_bound_states(1, my_spin)
1691 : localized_wfn_control%lu_bound_states(2, ispin) = &
1692 76 : localized_wfn_control%lu_bound_states(1, my_spin) + nmoloc(ispin) - 1
1693 76 : max_iloc = localized_wfn_control%lu_bound_states(2, ispin)
1694 242 : DO i = 1, nmoloc(ispin)
1695 242 : localized_wfn_control%loc_states(i, ispin) = localized_wfn_control%lu_bound_states(1, ispin) + i - 1
1696 : END DO
1697 76 : CPASSERT(max_iloc <= nmo(ispin))
1698 38 : MARK_USED(nmo)
1699 : END DO
1700 : CASE (energy_loc_range)
1701 : ! Energy
1702 30 : ALLOCATE (localized_wfn_control%loc_states(max_nmoloc, 2))
1703 202 : localized_wfn_control%loc_states = 0
1704 22 : DO ispin = 1, nspins
1705 128 : DO i = 1, nmoloc(ispin)
1706 118 : localized_wfn_control%loc_states(i, ispin) = localized_wfn_control%lu_bound_states(1, ispin) + i - 1
1707 : END DO
1708 : END DO
1709 : CASE (state_loc_all)
1710 : ! All
1711 834 : ALLOCATE (localized_wfn_control%loc_states(max_nmoloc, 2))
1712 5678 : localized_wfn_control%loc_states = 0
1713 :
1714 278 : IF (localized_wfn_control%lu_bound_states(1, 1) == 1) THEN
1715 660 : DO ispin = 1, nspins
1716 382 : localized_wfn_control%lu_bound_states(1, ispin) = 1
1717 382 : localized_wfn_control%lu_bound_states(2, ispin) = nmoloc(ispin)
1718 382 : IF (nmoloc(ispin) < 1) localized_wfn_control%lu_bound_states(1, ispin) = 0
1719 3808 : DO i = 1, nmoloc(ispin)
1720 3530 : localized_wfn_control%loc_states(i, ispin) = i
1721 : END DO
1722 : END DO
1723 : ELSE
1724 0 : DO ispin = 1, nspins
1725 0 : IF (nmoloc(ispin) < 1) localized_wfn_control%lu_bound_states(1, ispin) = 0
1726 0 : DO i = 1, nmoloc(ispin)
1727 : localized_wfn_control%loc_states(i, ispin) = &
1728 0 : localized_wfn_control%lu_bound_states(1, ispin) + i - 1
1729 : END DO
1730 : END DO
1731 : END IF
1732 : CASE (state_loc_mixed)
1733 : ! Mixed
1734 6 : ALLOCATE (localized_wfn_control%loc_states(max_nmoloc, 2))
1735 178 : localized_wfn_control%loc_states = 0
1736 350 : DO ispin = 1, nspins
1737 90 : DO i = 1, nmoloc(ispin)
1738 88 : localized_wfn_control%loc_states(i, ispin) = i
1739 : END DO
1740 : END DO
1741 : END SELECT
1742 :
1743 346 : CALL timestop(state)
1744 :
1745 346 : END SUBROUTINE set_loc_wfn_lists
1746 :
1747 : END MODULE qs_loc_utils
1748 :
|
