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