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