Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief calculates the electron transfer coupling elements by projection-operator approach
10 : !> Kondov et al. J.Phys.Chem.C 2007, 111, 11970-11981
11 : !> \author Z. Futera (02.2017)
12 : ! **************************************************************************************************
13 : MODULE et_coupling_proj
14 :
15 : USE atomic_kind_types, ONLY: atomic_kind_type,&
16 : get_atomic_kind
17 : USE basis_set_types, ONLY: get_gto_basis_set,&
18 : gto_basis_set_type
19 : USE bibliography, ONLY: Futera2017,&
20 : cite_reference
21 : USE cell_types, ONLY: cell_type
22 : USE cp_blacs_env, ONLY: cp_blacs_env_type
23 : USE cp_control_types, ONLY: dft_control_type
24 : USE cp_dbcsr_api, ONLY: dbcsr_p_type
25 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
26 : cp_dbcsr_sm_fm_multiply
27 : USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale
28 : USE cp_fm_diag, ONLY: choose_eigv_solver,&
29 : cp_fm_power
30 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
31 : cp_fm_struct_equivalent,&
32 : cp_fm_struct_release,&
33 : cp_fm_struct_type
34 : USE cp_fm_types, ONLY: &
35 : cp_fm_create, cp_fm_get_element, cp_fm_get_submatrix, cp_fm_release, cp_fm_set_all, &
36 : cp_fm_set_element, cp_fm_to_fm, cp_fm_to_fm_submat, cp_fm_type, cp_fm_vectorssum
37 : USE cp_log_handling, ONLY: cp_get_default_logger,&
38 : cp_logger_type,&
39 : cp_to_string
40 : USE cp_output_handling, ONLY: cp_p_file,&
41 : cp_print_key_finished_output,&
42 : cp_print_key_should_output,&
43 : cp_print_key_unit_nr
44 : USE cp_realspace_grid_cube, ONLY: cp_pw_to_cube
45 : USE input_section_types, ONLY: section_get_ivals,&
46 : section_get_lval,&
47 : section_vals_get,&
48 : section_vals_get_subs_vals,&
49 : section_vals_type,&
50 : section_vals_val_get
51 : USE kinds, ONLY: default_path_length,&
52 : default_string_length,&
53 : dp
54 : USE kpoint_types, ONLY: kpoint_type
55 : USE message_passing, ONLY: mp_para_env_type
56 : USE orbital_pointers, ONLY: nso
57 : USE parallel_gemm_api, ONLY: parallel_gemm
58 : USE particle_list_types, ONLY: particle_list_type
59 : USE particle_types, ONLY: particle_type
60 : USE physcon, ONLY: evolt
61 : USE pw_env_types, ONLY: pw_env_get,&
62 : pw_env_type
63 : USE pw_pool_types, ONLY: pw_pool_p_type,&
64 : pw_pool_type
65 : USE pw_types, ONLY: pw_c1d_gs_type,&
66 : pw_r3d_rs_type
67 : USE qs_collocate_density, ONLY: calculate_wavefunction
68 : USE qs_environment_types, ONLY: get_qs_env,&
69 : qs_environment_type
70 : USE qs_kind_types, ONLY: get_qs_kind,&
71 : get_qs_kind_set,&
72 : qs_kind_type
73 : USE qs_mo_methods, ONLY: make_mo_eig
74 : USE qs_mo_occupation, ONLY: set_mo_occupation
75 : USE qs_mo_types, ONLY: allocate_mo_set,&
76 : deallocate_mo_set,&
77 : mo_set_type
78 : USE qs_subsys_types, ONLY: qs_subsys_get,&
79 : qs_subsys_type
80 : USE scf_control_types, ONLY: scf_control_type
81 : #include "./base/base_uses.f90"
82 :
83 : IMPLICIT NONE
84 :
85 : PRIVATE
86 :
87 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'et_coupling_proj'
88 :
89 : ! Electronic-coupling calculation data structure
90 : !
91 : ! n_atoms - number of atoms in the blocks
92 : ! n_blocks - number of atomic blocks (donor,acceptor,bridge,...)
93 : ! fermi - system Fermi level (alpha/beta spin component)
94 : ! m_transf - transformation matrix for basis-set orthogonalization (S^{-1/2})
95 : ! m_transf_inv - inversion transformation matrix
96 : ! block - atomic data blocks
97 : TYPE et_cpl
98 : INTEGER :: n_atoms = 0
99 : INTEGER :: n_blocks = 0
100 : REAL(KIND=dp), DIMENSION(:), POINTER :: fermi => NULL()
101 : TYPE(cp_fm_type), POINTER :: m_transf => NULL()
102 : TYPE(cp_fm_type), POINTER :: m_transf_inv => NULL()
103 : TYPE(et_cpl_block), DIMENSION(:), POINTER :: block => NULL()
104 : END TYPE et_cpl
105 :
106 : ! Electronic-coupling data block
107 : !
108 : ! n_atoms - number of atoms
109 : ! n_electrons - number of electrons
110 : ! n_ao - number of AO basis functions
111 : ! atom - list of atoms
112 : ! mo - electronic states
113 : ! hab - electronic-coupling elements
114 : TYPE et_cpl_block
115 : INTEGER :: n_atoms = 0
116 : INTEGER :: n_electrons = 0
117 : INTEGER :: n_ao = 0
118 : TYPE(et_cpl_atom), DIMENSION(:), POINTER :: atom => NULL()
119 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mo => NULL()
120 : TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: hab => NULL()
121 : END TYPE et_cpl_block
122 :
123 : ! Electronic-coupling block-atom data
124 : ! id - atom ID
125 : ! n_ao - number of AO basis functions
126 : ! ao_pos - position of atom in array of AO functions
127 : TYPE et_cpl_atom
128 : INTEGER :: id = 0
129 : INTEGER :: n_ao = 0
130 : INTEGER :: ao_pos = 0
131 : END TYPE et_cpl_atom
132 :
133 : PUBLIC :: calc_et_coupling_proj
134 :
135 : CONTAINS
136 :
137 : ! **************************************************************************************************
138 : !> \brief Release memory allocate for electronic coupling data structures
139 : !> \param ec electronic coupling data structure
140 : !> \author Z. Futera (02.2017)
141 : ! **************************************************************************************************
142 10 : SUBROUTINE release_ec_data(ec)
143 :
144 : ! Routine arguments
145 : TYPE(et_cpl), POINTER :: ec
146 :
147 : INTEGER :: i, j
148 :
149 : ! Routine name for debug purposes
150 :
151 10 : IF (ASSOCIATED(ec)) THEN
152 :
153 10 : IF (ASSOCIATED(ec%fermi)) &
154 10 : DEALLOCATE (ec%fermi)
155 10 : IF (ASSOCIATED(ec%m_transf)) THEN
156 10 : CALL cp_fm_release(matrix=ec%m_transf)
157 10 : DEALLOCATE (ec%m_transf)
158 10 : NULLIFY (ec%m_transf)
159 : END IF
160 10 : IF (ASSOCIATED(ec%m_transf_inv)) THEN
161 10 : CALL cp_fm_release(matrix=ec%m_transf_inv)
162 10 : DEALLOCATE (ec%m_transf_inv)
163 10 : NULLIFY (ec%m_transf_inv)
164 : END IF
165 :
166 10 : IF (ASSOCIATED(ec%block)) THEN
167 :
168 30 : DO i = 1, SIZE(ec%block)
169 20 : IF (ASSOCIATED(ec%block(i)%atom)) &
170 20 : DEALLOCATE (ec%block(i)%atom)
171 20 : IF (ASSOCIATED(ec%block(i)%mo)) THEN
172 60 : DO j = 1, SIZE(ec%block(i)%mo)
173 60 : CALL deallocate_mo_set(ec%block(i)%mo(j))
174 : END DO
175 20 : DEALLOCATE (ec%block(i)%mo)
176 : END IF
177 30 : CALL cp_fm_release(ec%block(i)%hab)
178 : END DO
179 :
180 10 : DEALLOCATE (ec%block)
181 :
182 : END IF
183 :
184 10 : DEALLOCATE (ec)
185 :
186 : END IF
187 :
188 10 : END SUBROUTINE release_ec_data
189 :
190 : ! **************************************************************************************************
191 : !> \brief check the electronic-coupling input section and set the atomic block data
192 : !> \param qs_env QuickStep environment containing all system data
193 : !> \param et_proj_sec the electronic-coupling input section
194 : !> \param ec electronic coupling data structure
195 : !> \author Z. Futera (02.2017)
196 : ! **************************************************************************************************
197 10 : SUBROUTINE set_block_data(qs_env, et_proj_sec, ec)
198 :
199 : ! Routine arguments
200 : TYPE(qs_environment_type), POINTER :: qs_env
201 : TYPE(section_vals_type), POINTER :: et_proj_sec
202 : TYPE(et_cpl), POINTER :: ec
203 :
204 : INTEGER :: i, j, k, l, n, n_ao, n_atoms, n_set
205 10 : INTEGER, DIMENSION(:), POINTER :: atom_id, atom_nf, atom_ps, n_shell, t
206 10 : INTEGER, DIMENSION(:, :), POINTER :: ang_mom_id
207 : LOGICAL :: found
208 : TYPE(gto_basis_set_type), POINTER :: ao_basis_set
209 10 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
210 10 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
211 : TYPE(section_vals_type), POINTER :: block_sec
212 :
213 : ! Routine name for debug purposes
214 :
215 10 : NULLIFY (ao_basis_set)
216 10 : NULLIFY (particle_set)
217 10 : NULLIFY (qs_kind_set)
218 10 : NULLIFY (n_shell)
219 10 : NULLIFY (ang_mom_id)
220 10 : NULLIFY (atom_nf)
221 10 : NULLIFY (atom_id)
222 10 : NULLIFY (block_sec)
223 :
224 : ! Initialization
225 10 : ec%n_atoms = 0
226 10 : ec%n_blocks = 0
227 10 : NULLIFY (ec%fermi)
228 10 : NULLIFY (ec%m_transf)
229 10 : NULLIFY (ec%m_transf_inv)
230 10 : NULLIFY (ec%block)
231 :
232 : ! Number of atoms / atom types
233 10 : CALL get_qs_env(qs_env, particle_set=particle_set, qs_kind_set=qs_kind_set, natom=n_atoms)
234 : ! Number of AO basis functions
235 10 : CALL get_qs_kind_set(qs_kind_set, nsgf=n_ao)
236 :
237 : ! Number of AO functions per atom
238 30 : ALLOCATE (atom_nf(n_atoms))
239 10 : CPASSERT(ASSOCIATED(atom_nf))
240 :
241 82 : atom_nf = 0
242 82 : DO i = 1, n_atoms
243 72 : CALL get_atomic_kind(particle_set(i)%atomic_kind, kind_number=j)
244 72 : CALL get_qs_kind(qs_kind_set(j), basis_set=ao_basis_set)
245 72 : IF (.NOT. ASSOCIATED(ao_basis_set)) &
246 0 : CPABORT('Unsupported basis set type. ')
247 : CALL get_gto_basis_set(gto_basis_set=ao_basis_set, &
248 72 : nset=n_set, nshell=n_shell, l=ang_mom_id)
249 238 : DO j = 1, n_set
250 280 : DO k = 1, n_shell(j)
251 208 : atom_nf(i) = atom_nf(i) + nso(ang_mom_id(k, j))
252 : END DO
253 : END DO
254 : END DO
255 :
256 : ! Sanity check
257 10 : n = 0
258 82 : DO i = 1, n_atoms
259 82 : n = n + atom_nf(i)
260 : END DO
261 10 : CPASSERT(n == n_ao)
262 :
263 : ! Atom position in AO array
264 30 : ALLOCATE (atom_ps(n_atoms))
265 10 : CPASSERT(ASSOCIATED(atom_ps))
266 82 : atom_ps = 1
267 72 : DO i = 1, n_atoms - 1
268 72 : atom_ps(i + 1) = atom_ps(i) + atom_nf(i)
269 : END DO
270 :
271 : ! Number of blocks
272 10 : block_sec => section_vals_get_subs_vals(et_proj_sec, 'BLOCK')
273 10 : CALL section_vals_get(block_sec, n_repetition=ec%n_blocks)
274 50 : ALLOCATE (ec%block(ec%n_blocks))
275 10 : CPASSERT(ASSOCIATED(ec%block))
276 :
277 : ! Block data
278 30 : ALLOCATE (t(n_atoms))
279 10 : CPASSERT(ASSOCIATED(t))
280 :
281 10 : ec%n_atoms = 0
282 30 : DO i = 1, ec%n_blocks
283 :
284 : ! Initialization
285 20 : ec%block(i)%n_atoms = 0
286 20 : ec%block(i)%n_electrons = 0
287 20 : ec%block(i)%n_ao = 0
288 20 : NULLIFY (ec%block(i)%atom)
289 20 : NULLIFY (ec%block(i)%mo)
290 20 : NULLIFY (ec%block(i)%hab)
291 :
292 : ! Number of electrons
293 : CALL section_vals_val_get(block_sec, i_rep_section=i, &
294 20 : keyword_name='NELECTRON', i_val=ec%block(i)%n_electrons)
295 :
296 : ! User-defined atom array
297 : CALL section_vals_val_get(block_sec, i_rep_section=i, &
298 20 : keyword_name='ATOMS', i_vals=atom_id)
299 :
300 : ! Count unique atoms
301 92 : DO j = 1, SIZE(atom_id)
302 : ! Check atom ID validity
303 72 : IF (atom_id(j) < 1 .OR. atom_id(j) > n_atoms) &
304 0 : CPABORT('invalid fragment atom ID ('//TRIM(ADJUSTL(cp_to_string(atom_id(j))))//')')
305 : ! Check if the atom is not in previously-defined blocks
306 : found = .FALSE.
307 108 : DO k = 1, i - 1
308 348 : DO l = 1, ec%block(k)%n_atoms
309 276 : IF (ec%block(k)%atom(l)%id == atom_id(j)) THEN
310 0 : CPWARN('multiple definition of atom'//TRIM(ADJUSTL(cp_to_string(atom_id(j)))))
311 0 : found = .TRUE.
312 0 : EXIT
313 : END IF
314 : END DO
315 : END DO
316 : ! Check if the atom is not in already defined in the present block
317 72 : IF (.NOT. found) THEN
318 276 : DO k = 1, ec%block(i)%n_atoms
319 276 : IF (t(k) == atom_id(j)) THEN
320 0 : CPWARN('multiple definition of atom'//TRIM(ADJUSTL(cp_to_string(atom_id(j)))))
321 : found = .TRUE.
322 : EXIT
323 : END IF
324 : END DO
325 : END IF
326 : ! Save the atom
327 20 : IF (.NOT. found) THEN
328 72 : ec%block(i)%n_atoms = ec%block(i)%n_atoms + 1
329 72 : t(ec%block(i)%n_atoms) = atom_id(j)
330 : END IF
331 : END DO
332 :
333 : ! Memory allocation
334 132 : ALLOCATE (ec%block(i)%atom(ec%block(i)%n_atoms))
335 20 : CPASSERT(ASSOCIATED(ec%block(i)%atom))
336 :
337 : ! Save atom IDs and number of AOs
338 92 : DO j = 1, ec%block(i)%n_atoms
339 72 : ec%block(i)%atom(j)%id = t(j)
340 72 : ec%block(i)%atom(j)%n_ao = atom_nf(ec%block(i)%atom(j)%id)
341 72 : ec%block(i)%atom(j)%ao_pos = atom_ps(ec%block(i)%atom(j)%id)
342 92 : ec%block(i)%n_ao = ec%block(i)%n_ao + ec%block(i)%atom(j)%n_ao
343 : END DO
344 :
345 30 : ec%n_atoms = ec%n_atoms + ec%block(i)%n_atoms
346 : END DO
347 :
348 : ! Clean memory
349 10 : IF (ASSOCIATED(atom_nf)) &
350 10 : DEALLOCATE (atom_nf)
351 10 : IF (ASSOCIATED(atom_ps)) &
352 10 : DEALLOCATE (atom_ps)
353 10 : IF (ASSOCIATED(t)) &
354 10 : DEALLOCATE (t)
355 :
356 10 : END SUBROUTINE set_block_data
357 :
358 : ! **************************************************************************************************
359 : !> \brief check the electronic-coupling input section and set the atomic block data
360 : !> \param ec electronic coupling data structure
361 : !> \param fa system Fermi level (alpha spin)
362 : !> \param fb system Fermi level (beta spin)
363 : !> \author Z. Futera (02.2017)
364 : ! **************************************************************************************************
365 10 : SUBROUTINE set_fermi(ec, fa, fb)
366 :
367 : ! Routine arguments
368 : TYPE(et_cpl), POINTER :: ec
369 : REAL(KIND=dp) :: fa
370 : REAL(KIND=dp), OPTIONAL :: fb
371 :
372 : ! Routine name for debug purposes
373 :
374 10 : NULLIFY (ec%fermi)
375 :
376 10 : IF (PRESENT(fb)) THEN
377 :
378 10 : ALLOCATE (ec%fermi(2))
379 10 : CPASSERT(ASSOCIATED(ec%fermi))
380 10 : ec%fermi(1) = fa
381 10 : ec%fermi(2) = fb
382 :
383 : ELSE
384 :
385 0 : ALLOCATE (ec%fermi(1))
386 0 : CPASSERT(ASSOCIATED(ec%fermi))
387 0 : ec%fermi(1) = fa
388 :
389 : END IF
390 :
391 10 : END SUBROUTINE set_fermi
392 :
393 : ! **************************************************************************************************
394 : !> \brief reorder Hamiltonian matrix according to defined atomic blocks
395 : !> \param ec electronic coupling data structure
396 : !> \param mat_h the Hamiltonian matrix
397 : !> \param mat_w working matrix of the same dimension
398 : !> \author Z. Futera (02.2017)
399 : ! **************************************************************************************************
400 20 : SUBROUTINE reorder_hamiltonian_matrix(ec, mat_h, mat_w)
401 :
402 : ! Routine arguments
403 : TYPE(et_cpl), POINTER :: ec
404 : TYPE(cp_fm_type), INTENT(IN) :: mat_h, mat_w
405 :
406 : INTEGER :: ic, ir, jc, jr, kc, kr, mc, mr, nc, nr
407 : REAL(KIND=dp) :: xh
408 :
409 : ! Routine name for debug purposes
410 : ! Local variables
411 :
412 20 : IF (.NOT. cp_fm_struct_equivalent(mat_h%matrix_struct, mat_w%matrix_struct)) &
413 0 : CPABORT('cannot reorder Hamiltonian, working-matrix structure is not equivalent')
414 :
415 : ! Matrix-element reordering
416 20 : nr = 1
417 : ! Rows
418 60 : DO ir = 1, ec%n_blocks
419 204 : DO jr = 1, ec%block(ir)%n_atoms
420 592 : DO kr = 1, ec%block(ir)%atom(jr)%n_ao
421 : ! Columns
422 408 : nc = 1
423 1224 : DO ic = 1, ec%n_blocks
424 6072 : DO jc = 1, ec%block(ic)%n_atoms
425 18384 : DO kc = 1, ec%block(ic)%atom(jc)%n_ao
426 12720 : mr = ec%block(ir)%atom(jr)%ao_pos + kr - 1
427 12720 : mc = ec%block(ic)%atom(jc)%ao_pos + kc - 1
428 12720 : CALL cp_fm_get_element(mat_h, nr, nc, xh)
429 12720 : CALL cp_fm_set_element(mat_w, nr, nc, xh)
430 30288 : nc = nc + 1
431 : END DO
432 : END DO
433 : END DO
434 552 : nr = nr + 1
435 : END DO
436 : END DO
437 : END DO
438 :
439 : ! Copy the reordered matrix to original data array
440 20 : CALL cp_fm_to_fm(mat_w, mat_h)
441 :
442 20 : END SUBROUTINE reorder_hamiltonian_matrix
443 :
444 : ! **************************************************************************************************
445 : !> \brief calculated transformation matrix for basis-set orthogonalization (S^{-1/2})
446 : !> \param qs_env QuickStep environment containing all system data
447 : !> \param mat_t storage for the transformation matrix
448 : !> \param mat_i storage for the inversion transformation matrix
449 : !> \param mat_w working matrix of the same dimension
450 : !> \author Z. Futera (02.2017)
451 : ! **************************************************************************************************
452 10 : SUBROUTINE get_s_half_inv_matrix(qs_env, mat_t, mat_i, mat_w)
453 :
454 : ! Routine arguments
455 : TYPE(qs_environment_type), POINTER :: qs_env
456 : TYPE(cp_fm_type), INTENT(INOUT) :: mat_t, mat_i
457 : TYPE(cp_fm_type), INTENT(IN) :: mat_w
458 :
459 : INTEGER :: n_deps
460 10 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_s
461 : TYPE(scf_control_type), POINTER :: scf_cntrl
462 :
463 : ! Routine name for debug purposes
464 :
465 10 : NULLIFY (mat_s)
466 10 : NULLIFY (scf_cntrl)
467 :
468 10 : CALL get_qs_env(qs_env, matrix_s=mat_s)
469 10 : CALL copy_dbcsr_to_fm(mat_s(1)%matrix, mat_t)
470 10 : CALL copy_dbcsr_to_fm(mat_s(1)%matrix, mat_i)
471 :
472 : ! Transformation S -> S^{-1/2}
473 10 : CALL get_qs_env(qs_env, scf_control=scf_cntrl)
474 10 : CALL cp_fm_power(mat_t, mat_w, -0.5_dp, scf_cntrl%eps_eigval, n_deps)
475 10 : CALL cp_fm_power(mat_i, mat_w, +0.5_dp, scf_cntrl%eps_eigval, n_deps)
476 : ! Sanity check
477 10 : IF (n_deps /= 0) THEN
478 : CALL cp_warn(__LOCATION__, &
479 : "Overlap matrix exhibits linear dependencies. At least some "// &
480 0 : "eigenvalues have been quenched.")
481 : END IF
482 :
483 10 : END SUBROUTINE get_s_half_inv_matrix
484 :
485 : ! **************************************************************************************************
486 : !> \brief transform KS hamiltonian to orthogonalized block-separated basis set
487 : !> \param qs_env QuickStep environment containing all system data
488 : !> \param ec electronic coupling data structure
489 : !> \param fm_s full-matrix structure used for allocation of KS matrices
490 : !> \param mat_t storage for pointers to the transformed KS matrices
491 : !> \param mat_w working matrix of the same dimension
492 : !> \param n_ao total number of AO basis functions
493 : !> \param n_spins number of spin components
494 : !> \author Z. Futera (02.2017)
495 : ! **************************************************************************************************
496 10 : SUBROUTINE get_block_hamiltonian(qs_env, ec, fm_s, mat_t, mat_w, n_ao, n_spins)
497 :
498 : ! Routine arguments
499 : TYPE(qs_environment_type), POINTER :: qs_env
500 : TYPE(et_cpl), POINTER :: ec
501 : TYPE(cp_fm_struct_type), POINTER :: fm_s
502 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:), &
503 : INTENT(OUT) :: mat_t
504 : TYPE(cp_fm_type), INTENT(IN) :: mat_w
505 : INTEGER :: n_ao, n_spins
506 :
507 : INTEGER :: i
508 10 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_h
509 :
510 : ! Routine name for debug purposes
511 :
512 10 : NULLIFY (mat_h)
513 :
514 : ! Memory allocation
515 50 : ALLOCATE (mat_t(n_spins))
516 :
517 : ! KS Hamiltonian
518 10 : CALL get_qs_env(qs_env, matrix_ks=mat_h)
519 : ! Transformation matrix
520 10 : ALLOCATE (ec%m_transf, ec%m_transf_inv)
521 : CALL cp_fm_create(matrix=ec%m_transf, matrix_struct=fm_s, &
522 10 : name='S^(-1/2) TRANSFORMATION MATRIX')
523 : CALL cp_fm_create(matrix=ec%m_transf_inv, matrix_struct=fm_s, &
524 10 : name='S^(+1/2) TRANSFORMATION MATRIX')
525 10 : CALL get_s_half_inv_matrix(qs_env, ec%m_transf, ec%m_transf_inv, mat_w)
526 :
527 30 : DO i = 1, n_spins
528 :
529 : ! Full-matrix format
530 : CALL cp_fm_create(matrix=mat_t(i), matrix_struct=fm_s, &
531 20 : name='KS HAMILTONIAN IN SEPARATED ORTHOGONALIZED BASIS SET')
532 20 : CALL copy_dbcsr_to_fm(mat_h(i)%matrix, mat_t(i))
533 :
534 : ! Transform KS Hamiltonian to the orthogonalized AO basis set
535 20 : CALL parallel_gemm("N", "N", n_ao, n_ao, n_ao, 1.0_dp, ec%m_transf, mat_t(i), 0.0_dp, mat_w)
536 20 : CALL parallel_gemm("N", "N", n_ao, n_ao, n_ao, 1.0_dp, mat_w, ec%m_transf, 0.0_dp, mat_t(i))
537 :
538 : ! Reorder KS Hamiltonain elements to defined block structure
539 30 : CALL reorder_hamiltonian_matrix(ec, mat_t(i), mat_w)
540 :
541 : END DO
542 :
543 10 : END SUBROUTINE get_block_hamiltonian
544 :
545 : ! **************************************************************************************************
546 : !> \brief Diagonalize diagonal blocks of the KS hamiltonian in separated orthogonalized basis set
547 : !> \param qs_env QuickStep environment containing all system data
548 : !> \param ec electronic coupling data structure
549 : !> \param mat_h Hamiltonian in separated orthogonalized basis set
550 : !> \author Z. Futera (02.2017)
551 : ! **************************************************************************************************
552 10 : SUBROUTINE hamiltonian_block_diag(qs_env, ec, mat_h)
553 :
554 : ! Routine arguments
555 : TYPE(qs_environment_type), POINTER :: qs_env
556 : TYPE(et_cpl), POINTER :: ec
557 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mat_h
558 :
559 : INTEGER :: i, j, k, l, n_spins, spin
560 10 : REAL(KIND=dp), DIMENSION(:), POINTER :: vec_e
561 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
562 : TYPE(cp_fm_struct_type), POINTER :: fm_s
563 : TYPE(cp_fm_type) :: mat_u
564 10 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: dat
565 : TYPE(mp_para_env_type), POINTER :: para_env
566 :
567 : ! Routine name for debug purposes
568 :
569 10 : NULLIFY (vec_e)
570 10 : NULLIFY (blacs_env)
571 10 : NULLIFY (para_env)
572 10 : NULLIFY (fm_s)
573 :
574 : ! Parallel environment
575 10 : CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
576 :
577 : ! Storage for block sub-matrices
578 50 : ALLOCATE (dat(ec%n_blocks))
579 10 : CPASSERT(ALLOCATED(dat))
580 :
581 : ! Storage for electronic states and couplings
582 10 : n_spins = SIZE(mat_h)
583 30 : DO i = 1, ec%n_blocks
584 100 : ALLOCATE (ec%block(i)%mo(n_spins))
585 20 : CPASSERT(ASSOCIATED(ec%block(i)%mo))
586 200 : ALLOCATE (ec%block(i)%hab(n_spins, ec%n_blocks))
587 30 : CPASSERT(ASSOCIATED(ec%block(i)%hab))
588 : END DO
589 :
590 : ! Spin components
591 30 : DO spin = 1, n_spins
592 :
593 : ! Diagonal blocks
594 20 : j = 1
595 60 : DO i = 1, ec%n_blocks
596 :
597 : ! Memory allocation
598 : CALL cp_fm_struct_create(fmstruct=fm_s, para_env=para_env, context=blacs_env, &
599 40 : nrow_global=ec%block(i)%n_ao, ncol_global=ec%block(i)%n_ao)
600 : CALL cp_fm_create(matrix=dat(i), matrix_struct=fm_s, &
601 40 : name='H_KS DIAGONAL BLOCK')
602 :
603 120 : ALLOCATE (vec_e(ec%block(i)%n_ao))
604 40 : CPASSERT(ASSOCIATED(vec_e))
605 :
606 : ! Copy block data
607 : CALL cp_fm_to_fm_submat(mat_h(spin), &
608 : dat(i), ec%block(i)%n_ao, &
609 40 : ec%block(i)%n_ao, j, j, 1, 1)
610 :
611 : ! Diagonalization
612 40 : CALL cp_fm_create(matrix=mat_u, matrix_struct=fm_s, name='UNITARY MATRIX')
613 40 : CALL choose_eigv_solver(dat(i), mat_u, vec_e)
614 40 : CALL cp_fm_to_fm(mat_u, dat(i))
615 :
616 : ! Save state energies / vectors
617 40 : CALL create_block_mo_set(qs_env, ec, i, spin, mat_u, vec_e)
618 :
619 : ! Clean memory
620 40 : CALL cp_fm_struct_release(fmstruct=fm_s)
621 40 : CALL cp_fm_release(matrix=mat_u)
622 40 : DEALLOCATE (vec_e)
623 :
624 : ! Off-set for next block
625 100 : j = j + ec%block(i)%n_ao
626 :
627 : END DO
628 :
629 : ! Off-diagonal blocks
630 20 : k = 1
631 60 : DO i = 1, ec%n_blocks
632 40 : l = 1
633 120 : DO j = 1, ec%n_blocks
634 80 : IF (i /= j) THEN
635 :
636 : ! Memory allocation
637 : CALL cp_fm_struct_create(fmstruct=fm_s, para_env=para_env, context=blacs_env, &
638 40 : nrow_global=ec%block(i)%n_ao, ncol_global=ec%block(j)%n_ao)
639 : CALL cp_fm_create(matrix=ec%block(i)%hab(spin, j), matrix_struct=fm_s, &
640 40 : name='H_KS OFF-DIAGONAL BLOCK')
641 :
642 : ! Copy block data
643 : CALL cp_fm_to_fm_submat(mat_h(spin), &
644 : ec%block(i)%hab(spin, j), ec%block(i)%n_ao, &
645 40 : ec%block(j)%n_ao, k, l, 1, 1)
646 :
647 : ! Transformation
648 40 : CALL cp_fm_create(matrix=mat_u, matrix_struct=fm_s, name='FULL WORK MATRIX')
649 : CALL parallel_gemm("T", "N", ec%block(i)%n_ao, ec%block(j)%n_ao, ec%block(i)%n_ao, &
650 40 : 1.0_dp, dat(i), ec%block(i)%hab(spin, j), 0.0_dp, mat_u)
651 : CALL parallel_gemm("N", "N", ec%block(i)%n_ao, ec%block(j)%n_ao, ec%block(j)%n_ao, &
652 40 : 1.0_dp, mat_u, dat(j), 0.0_dp, ec%block(i)%hab(spin, j))
653 :
654 : ! Clean memory
655 40 : CALL cp_fm_struct_release(fmstruct=fm_s)
656 40 : CALL cp_fm_release(matrix=mat_u)
657 :
658 : END IF
659 : ! Off-set for next block
660 120 : l = l + ec%block(j)%n_ao
661 : END DO
662 : ! Off-set for next block
663 60 : k = k + ec%block(i)%n_ao
664 : END DO
665 :
666 : ! Clean memory
667 30 : IF (ALLOCATED(dat)) THEN
668 60 : DO i = 1, SIZE(dat)
669 60 : CALL cp_fm_release(dat(i))
670 : END DO
671 : END IF
672 : END DO
673 :
674 : ! Clean memory
675 10 : IF (ALLOCATED(dat)) &
676 10 : DEALLOCATE (dat)
677 :
678 20 : END SUBROUTINE hamiltonian_block_diag
679 :
680 : ! **************************************************************************************************
681 : !> \brief Return sum of selected squared MO coefficients
682 : !> \param blk_at list of atoms in the block
683 : !> \param mo array of MO sets
684 : !> \param id state index
685 : !> \param atom list of atoms for MO coefficient summing
686 : !> \return ...
687 : !> \author Z. Futera (02.2017)
688 : ! **************************************************************************************************
689 0 : FUNCTION get_mo_c2_sum(blk_at, mo, id, atom) RESULT(c2)
690 :
691 : ! Routine arguments
692 : TYPE(et_cpl_atom), DIMENSION(:), POINTER :: blk_at
693 : TYPE(cp_fm_type), INTENT(IN) :: mo
694 : INTEGER, INTENT(IN) :: id
695 : INTEGER, DIMENSION(:), POINTER :: atom
696 : REAL(KIND=dp) :: c2
697 :
698 : INTEGER :: i, ir, j, k
699 : LOGICAL :: found
700 : REAL(KIND=dp) :: c
701 :
702 : ! Returning value
703 : ! Routine name for debug purposes
704 : ! Local variables
705 :
706 : ! initialization
707 0 : c2 = 0.0d0
708 :
709 : ! selected atoms
710 0 : DO i = 1, SIZE(atom)
711 :
712 : ! find atomic function offset
713 0 : found = .FALSE.
714 0 : DO j = 1, SIZE(blk_at)
715 0 : IF (blk_at(j)%id == atom(i)) THEN
716 : found = .TRUE.
717 : EXIT
718 : END IF
719 : END DO
720 :
721 0 : IF (.NOT. found) &
722 0 : CPABORT('MO-fraction atom ID not defined in the block')
723 :
724 : ! sum MO coefficients from the atom
725 0 : DO k = 1, blk_at(j)%n_ao
726 0 : ir = blk_at(j)%ao_pos + k - 1
727 0 : CALL cp_fm_get_element(mo, ir, id, c)
728 0 : c2 = c2 + c*c
729 : END DO
730 :
731 : END DO
732 :
733 0 : END FUNCTION get_mo_c2_sum
734 :
735 : ! **************************************************************************************************
736 : !> \brief Print out specific MO coefficients
737 : !> \param output_unit unit number of the open output stream
738 : !> \param qs_env QuickStep environment containing all system data
739 : !> \param ec electronic coupling data structure
740 : !> \param blk atomic-block ID
741 : !> \param n_spins number of spin components
742 : !> \author Z. Futera (02.2017)
743 : ! **************************************************************************************************
744 20 : SUBROUTINE print_mo_coeff(output_unit, qs_env, ec, blk, n_spins)
745 :
746 : ! Routine arguments
747 : INTEGER, INTENT(IN) :: output_unit
748 : TYPE(qs_environment_type), POINTER :: qs_env
749 : TYPE(et_cpl), POINTER :: ec
750 : INTEGER, INTENT(IN) :: blk, n_spins
751 :
752 : INTEGER :: j, k, l, m, n, n_ao, n_mo
753 20 : INTEGER, DIMENSION(:), POINTER :: list_at, list_mo
754 : REAL(KIND=dp) :: c1, c2
755 20 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: mat_w
756 20 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
757 : TYPE(section_vals_type), POINTER :: block_sec, print_sec
758 :
759 : ! Routine name for debug purposes
760 :
761 20 : NULLIFY (block_sec)
762 20 : NULLIFY (print_sec)
763 20 : NULLIFY (qs_kind_set)
764 :
765 : ! Atomic block data
766 : block_sec => section_vals_get_subs_vals(qs_env%input, &
767 40 : 'PROPERTIES%ET_COUPLING%PROJECTION%BLOCK')
768 :
769 20 : print_sec => section_vals_get_subs_vals(block_sec, 'PRINT', i_rep_section=blk)
770 :
771 : ! List of atoms
772 20 : CALL section_vals_val_get(print_sec, keyword_name='MO_COEFF_ATOM', n_rep_val=n)
773 :
774 20 : IF (n > 0) THEN
775 :
776 0 : IF (output_unit > 0) &
777 0 : WRITE (output_unit, '(/,T3,A/)') 'Block state fractions:'
778 :
779 : ! Number of AO functions
780 0 : CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set)
781 0 : CALL get_qs_kind_set(qs_kind_set, nsgf=n_ao)
782 :
783 : ! MOs in orthonormal basis set
784 0 : ALLOCATE (mat_w(n_spins))
785 0 : DO j = 1, n_spins
786 0 : n_mo = ec%block(blk)%n_ao
787 : CALL cp_fm_create(matrix=mat_w(j), &
788 : matrix_struct=ec%block(blk)%mo(j)%mo_coeff%matrix_struct, &
789 0 : name='BLOCK MOs IN ORTHONORMAL BASIS SET')
790 : CALL parallel_gemm("N", "N", n_ao, n_mo, n_ao, 1.0_dp, ec%m_transf_inv, &
791 0 : ec%block(blk)%mo(j)%mo_coeff, 0.0_dp, mat_w(j))
792 : END DO
793 :
794 0 : DO j = 1, n
795 0 : NULLIFY (list_at)
796 : CALL section_vals_val_get(print_sec, keyword_name='MO_COEFF_ATOM', &
797 0 : i_rep_val=j, i_vals=list_at)
798 0 : IF (ASSOCIATED(list_at)) THEN
799 :
800 : ! List of states
801 0 : CALL section_vals_val_get(print_sec, keyword_name='MO_COEFF_ATOM_STATE', n_rep_val=m)
802 :
803 0 : IF (m > 0) THEN
804 :
805 0 : DO k = 1, m
806 0 : NULLIFY (list_mo)
807 : CALL section_vals_val_get(print_sec, keyword_name='MO_COEFF_ATOM_STATE', &
808 0 : i_rep_val=k, i_vals=list_mo)
809 0 : IF (ASSOCIATED(list_mo)) THEN
810 :
811 0 : IF (j > 1) THEN
812 0 : IF (output_unit > 0) &
813 0 : WRITE (output_unit, *)
814 : END IF
815 :
816 0 : DO l = 1, SIZE(list_mo)
817 :
818 0 : IF (n_spins > 1) THEN
819 : c1 = get_mo_c2_sum(ec%block(blk)%atom, mat_w(1), &
820 0 : list_mo(l), list_at)
821 : c2 = get_mo_c2_sum(ec%block(blk)%atom, mat_w(2), &
822 0 : list_mo(l), list_at)
823 0 : IF (output_unit > 0) &
824 0 : WRITE (output_unit, '(I5,A,I5,2F20.10)') j, ' /', list_mo(l), c1, c2
825 : ELSE
826 : c1 = get_mo_c2_sum(ec%block(blk)%atom, mat_w(1), &
827 0 : list_mo(l), list_at)
828 0 : IF (output_unit > 0) &
829 0 : WRITE (output_unit, '(I5,A,I5,F20.10)') j, ' /', list_mo(l), c1
830 : END IF
831 :
832 : END DO
833 :
834 : END IF
835 : END DO
836 :
837 : END IF
838 :
839 : END IF
840 : END DO
841 :
842 : ! Clean memory
843 0 : CALL cp_fm_release(mat_w)
844 :
845 : END IF
846 :
847 40 : END SUBROUTINE print_mo_coeff
848 :
849 : ! **************************************************************************************************
850 : !> \brief Print out electronic states (MOs)
851 : !> \param output_unit unit number of the open output stream
852 : !> \param mo array of MO sets
853 : !> \param n_spins number of spin components
854 : !> \param label output label
855 : !> \param mx_mo_a maximum number of alpha states to print out
856 : !> \param mx_mo_b maximum number of beta states to print out
857 : !> \param fermi print out Fermi level and number of electrons
858 : !> \author Z. Futera (02.2017)
859 : ! **************************************************************************************************
860 40 : SUBROUTINE print_states(output_unit, mo, n_spins, label, mx_mo_a, mx_mo_b, fermi)
861 :
862 : ! Routine arguments
863 : INTEGER, INTENT(IN) :: output_unit
864 : TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mo
865 : INTEGER, INTENT(IN) :: n_spins
866 : CHARACTER(LEN=*), INTENT(IN) :: label
867 : INTEGER, INTENT(IN), OPTIONAL :: mx_mo_a, mx_mo_b
868 : LOGICAL, INTENT(IN), OPTIONAL :: fermi
869 :
870 : INTEGER :: i, mx_a, mx_b, n
871 : LOGICAL :: prnt_fm
872 :
873 : ! Routine name for debug purposes
874 :
875 20 : prnt_fm = .FALSE.
876 20 : IF (PRESENT(fermi)) &
877 20 : prnt_fm = fermi
878 :
879 20 : IF (output_unit > 0) THEN
880 :
881 15 : WRITE (output_unit, '(/,T3,A/)') 'State energies ('//TRIM(ADJUSTL(label))//'):'
882 :
883 : ! Spin-polarized calculation
884 15 : IF (n_spins > 1) THEN
885 :
886 15 : mx_a = mo(1)%nmo
887 15 : IF (PRESENT(mx_mo_a)) &
888 10 : mx_a = MIN(mo(1)%nmo, mx_mo_a)
889 15 : mx_b = mo(2)%nmo
890 15 : IF (PRESENT(mx_mo_b)) &
891 10 : mx_b = MIN(mo(2)%nmo, mx_mo_b)
892 15 : n = MAX(mx_a, mx_b)
893 :
894 181 : DO i = 1, n
895 166 : WRITE (output_unit, '(T3,I10)', ADVANCE='no') i
896 166 : IF (i <= mx_a) THEN
897 : WRITE (output_unit, '(2F12.4)', ADVANCE='no') &
898 166 : mo(1)%occupation_numbers(i), mo(1)%eigenvalues(i)
899 : ELSE
900 0 : WRITE (output_unit, '(A)', ADVANCE='no') ' '
901 : END IF
902 166 : WRITE (output_unit, '(A)', ADVANCE='no') ' '
903 181 : IF (i <= mx_b) THEN
904 : WRITE (output_unit, '(2F12.4)') &
905 161 : mo(2)%occupation_numbers(i), mo(2)%eigenvalues(i)
906 : ELSE
907 5 : WRITE (output_unit, *)
908 : END IF
909 : END DO
910 :
911 15 : IF (prnt_fm) THEN
912 : WRITE (output_unit, '(/,T3,I10,F24.4,I10,F19.4)') &
913 15 : mo(1)%nelectron, mo(1)%mu, &
914 30 : mo(2)%nelectron, mo(2)%mu
915 : END IF
916 :
917 : ! Spin-restricted calculation
918 : ELSE
919 :
920 0 : mx_a = mo(1)%nmo
921 0 : IF (PRESENT(mx_mo_a)) &
922 0 : mx_a = MIN(mo(1)%nmo, mx_mo_a)
923 :
924 0 : DO i = 1, mx_a
925 : WRITE (output_unit, '(T3,I10,2F12.4)') &
926 0 : i, mo(1)%occupation_numbers(i), mo(1)%eigenvalues(i)
927 : END DO
928 :
929 0 : IF (prnt_fm) THEN
930 : WRITE (output_unit, '(/,T3,I10,F24.4)') &
931 0 : mo(1)%nelectron, mo(1)%mu
932 : END IF
933 :
934 : END IF
935 :
936 : END IF
937 :
938 20 : END SUBROUTINE print_states
939 :
940 : ! **************************************************************************************************
941 : !> \brief Print out donor-acceptor state couplings
942 : !> \param ec_sec ...
943 : !> \param output_unit unit number of the open output stream
944 : !> \param logger ...
945 : !> \param ec electronic coupling data structure
946 : !> \param mo ...
947 : !> \author Z. Futera (02.2017)
948 : ! **************************************************************************************************
949 10 : SUBROUTINE print_couplings(ec_sec, output_unit, logger, ec, mo)
950 :
951 : ! Routine arguments
952 : TYPE(section_vals_type), POINTER :: ec_sec
953 : INTEGER, INTENT(IN) :: output_unit
954 : TYPE(cp_logger_type), POINTER :: logger
955 : TYPE(et_cpl), POINTER :: ec
956 : TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mo
957 :
958 : CHARACTER(LEN=default_path_length) :: filename, my_pos, title
959 : INTEGER :: i, j, k, l, n_states(2), nc, nr, nspins, &
960 : unit_nr
961 : LOGICAL :: append
962 10 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: w1, w2
963 : TYPE(section_vals_type), POINTER :: print_key
964 :
965 : ! Routine name for debug purposes
966 : ! Local variables
967 :
968 10 : n_states = 0
969 30 : DO i = 1, SIZE(mo)
970 30 : n_states(i) = mo(i)%nmo
971 : END DO
972 10 : nspins = 1
973 10 : IF (n_states(2) > 0) nspins = 2
974 :
975 : print_key => section_vals_get_subs_vals(section_vals=ec_sec, &
976 10 : subsection_name="PRINT%COUPLINGS")
977 :
978 10 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), &
979 : cp_p_file)) THEN
980 :
981 10 : my_pos = "REWIND"
982 10 : append = section_get_lval(print_key, "APPEND")
983 10 : IF (append) THEN
984 0 : my_pos = "APPEND"
985 : END IF
986 :
987 10 : IF (output_unit > 0) &
988 5 : WRITE (output_unit, '(/,T3,A/)') 'Printing coupling elements to output files'
989 :
990 30 : DO i = 1, ec%n_blocks
991 40 : DO j = i + 1, ec%n_blocks
992 :
993 10 : nr = ec%block(i)%hab(1, j)%matrix_struct%nrow_global
994 10 : nc = ec%block(i)%hab(1, j)%matrix_struct%ncol_global
995 :
996 40 : ALLOCATE (w1(nr, nc))
997 10 : CPASSERT(ASSOCIATED(w1))
998 10 : CALL cp_fm_get_submatrix(ec%block(i)%hab(1, j), w1)
999 10 : IF (nspins > 1) THEN
1000 30 : ALLOCATE (w2(nr, nc))
1001 10 : CPASSERT(ASSOCIATED(w2))
1002 10 : CALL cp_fm_get_submatrix(ec%block(i)%hab(2, j), w2)
1003 : END IF
1004 :
1005 10 : IF (output_unit > 0) THEN
1006 :
1007 5 : WRITE (filename, '(a5,I1.1,a1,I1.1)') "ET_BL_", i, "-", j
1008 : unit_nr = cp_print_key_unit_nr(logger, ec_sec, "PRINT%COUPLINGS", extension=".elcoup", &
1009 5 : middle_name=TRIM(filename), file_position=my_pos, log_filename=.FALSE.)
1010 :
1011 5 : WRITE (title, *) 'Coupling elements [meV] between blocks:', i, j
1012 :
1013 5 : WRITE (unit_nr, *) TRIM(title)
1014 5 : IF (nspins > 1) THEN
1015 5 : WRITE (unit_nr, '(T3,A8,T13,A8,T28,A,A)') 'State A', 'State B', 'Coupling spin 1', ' Coupling spin 2'
1016 : ELSE
1017 0 : WRITE (unit_nr, '(T3,A8,T13,A8,T28,A)') 'State A', 'State B', 'Coupling'
1018 : END IF
1019 :
1020 55 : DO k = 1, MIN(ec%block(i)%n_ao, n_states(1))
1021 836 : DO l = 1, MIN(ec%block(j)%n_ao, n_states(1))
1022 :
1023 836 : IF (nspins > 1) THEN
1024 :
1025 : WRITE (unit_nr, '(T3,I5,T13,I5,T22,E20.6)', ADVANCE='no') &
1026 786 : k, l, w1(k, l)*evolt*1000.0_dp
1027 786 : IF ((k <= n_states(2)) .AND. (l <= n_states(2))) THEN
1028 : WRITE (unit_nr, '(E20.6)') &
1029 779 : w2(k, l)*evolt*1000.0_dp
1030 : ELSE
1031 7 : WRITE (unit_nr, *)
1032 : END IF
1033 :
1034 : ELSE
1035 :
1036 : WRITE (unit_nr, '(T3,I5,T13,I5,T22,E20.6)') &
1037 0 : k, l, w1(k, l)*evolt*1000.0_dp
1038 : END IF
1039 :
1040 : END DO
1041 55 : WRITE (unit_nr, *)
1042 : END DO
1043 5 : CALL cp_print_key_finished_output(unit_nr, logger, ec_sec, "PRINT%COUPLINGS")
1044 :
1045 : END IF
1046 :
1047 10 : IF (ASSOCIATED(w1)) DEALLOCATE (w1)
1048 30 : IF (ASSOCIATED(w2)) DEALLOCATE (w2)
1049 :
1050 : END DO
1051 : END DO
1052 :
1053 : END IF
1054 10 : END SUBROUTINE print_couplings
1055 :
1056 : ! **************************************************************************************************
1057 : !> \brief Normalize set of MO vectors
1058 : !> \param qs_env QuickStep environment containing all system data
1059 : !> \param mo storage for the MO data set
1060 : !> \param n_ao number of AO basis functions
1061 : !> \param n_mo number of block states
1062 : !> \author Z. Futera (02.2017)
1063 : ! **************************************************************************************************
1064 40 : SUBROUTINE normalize_mo_vectors(qs_env, mo, n_ao, n_mo)
1065 :
1066 : ! Routine arguments
1067 : TYPE(qs_environment_type), POINTER :: qs_env
1068 : TYPE(mo_set_type), POINTER :: mo
1069 : INTEGER, INTENT(IN) :: n_ao, n_mo
1070 :
1071 : REAL(KIND=dp), DIMENSION(:), POINTER :: vec_t
1072 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
1073 : TYPE(cp_fm_struct_type), POINTER :: fm_s
1074 : TYPE(cp_fm_type) :: mat_sc, mat_t
1075 40 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_s
1076 : TYPE(mp_para_env_type), POINTER :: para_env
1077 :
1078 : ! Routine name for debug purposes
1079 :
1080 : ! Initialization
1081 40 : NULLIFY (blacs_env)
1082 40 : NULLIFY (para_env)
1083 40 : NULLIFY (fm_s)
1084 40 : NULLIFY (mat_s)
1085 : NULLIFY (vec_t)
1086 :
1087 : ! Overlap matrix
1088 40 : CALL get_qs_env(qs_env, matrix_s=mat_s)
1089 :
1090 : ! Calculate S*C product
1091 : CALL cp_fm_create(matrix=mat_sc, matrix_struct=mo%mo_coeff%matrix_struct, &
1092 40 : name='S*C PRODUCT MATRIX')
1093 40 : CALL cp_dbcsr_sm_fm_multiply(mat_s(1)%matrix, mo%mo_coeff, mat_sc, n_mo)
1094 :
1095 : ! Calculate C^T*S*C
1096 40 : CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
1097 : CALL cp_fm_struct_create(fmstruct=fm_s, para_env=para_env, context=blacs_env, &
1098 40 : nrow_global=n_mo, ncol_global=n_mo)
1099 : CALL cp_fm_create(matrix=mat_t, matrix_struct=fm_s, &
1100 40 : name='C^T*S*C OVERLAP PRODUCT MATRIX')
1101 40 : CALL parallel_gemm('T', 'N', n_mo, n_mo, n_ao, 1.0_dp, mo%mo_coeff, mat_sc, 0.0_dp, mat_t)
1102 :
1103 : ! Normalization
1104 120 : ALLOCATE (vec_t(n_mo))
1105 40 : CPASSERT(ASSOCIATED(vec_t))
1106 40 : CALL cp_fm_vectorssum(mat_t, vec_t)
1107 448 : vec_t = 1.0_dp/DSQRT(vec_t)
1108 40 : CALL cp_fm_column_scale(mo%mo_coeff, vec_t)
1109 :
1110 : ! Clean memory
1111 40 : CALL cp_fm_struct_release(fmstruct=fm_s)
1112 40 : CALL cp_fm_release(matrix=mat_sc)
1113 40 : CALL cp_fm_release(matrix=mat_t)
1114 40 : IF (ASSOCIATED(vec_t)) &
1115 40 : DEALLOCATE (vec_t)
1116 :
1117 80 : END SUBROUTINE normalize_mo_vectors
1118 :
1119 : ! **************************************************************************************************
1120 : !> \brief Transform block MO coefficients to original non-orthogonal basis set and save them
1121 : !> \param qs_env QuickStep environment containing all system data
1122 : !> \param ec electronic coupling data structure
1123 : !> \param id block ID
1124 : !> \param mo storage for the MO data set
1125 : !> \param mat_u matrix of the block states
1126 : !> \param n_ao number of AO basis functions
1127 : !> \param n_mo number of block states
1128 : !> \author Z. Futera (02.2017)
1129 : ! **************************************************************************************************
1130 40 : SUBROUTINE set_mo_coefficients(qs_env, ec, id, mo, mat_u, n_ao, n_mo)
1131 :
1132 : ! Routine arguments
1133 : TYPE(qs_environment_type), POINTER :: qs_env
1134 : TYPE(et_cpl), POINTER :: ec
1135 : INTEGER, INTENT(IN) :: id
1136 : TYPE(mo_set_type), POINTER :: mo
1137 : TYPE(cp_fm_type), INTENT(IN) :: mat_u
1138 : INTEGER, INTENT(IN) :: n_ao, n_mo
1139 :
1140 : INTEGER :: ic, ir, jc, jr, mr, nc, nr
1141 : REAL(KIND=dp) :: xu
1142 : TYPE(cp_fm_type) :: mat_w
1143 :
1144 : ! Routine name for debug purposes
1145 : ! Local variables
1146 :
1147 : ! Working matrix
1148 : CALL cp_fm_create(matrix=mat_w, matrix_struct=mo%mo_coeff%matrix_struct, &
1149 40 : name='BLOCK MO-TRANSFORMATION WORKING MATRIX')
1150 40 : CALL cp_fm_set_all(mat_w, 0.0_dp)
1151 :
1152 : ! Matrix-element reordering
1153 40 : nr = 1
1154 : ! Rows
1155 184 : DO ir = 1, ec%block(id)%n_atoms
1156 592 : DO jr = 1, ec%block(id)%atom(ir)%n_ao
1157 : ! Columns
1158 408 : nc = 1
1159 2832 : DO ic = 1, ec%block(id)%n_atoms
1160 9192 : DO jc = 1, ec%block(id)%atom(ic)%n_ao
1161 6360 : mr = ec%block(id)%atom(ir)%ao_pos + jr - 1
1162 6360 : CALL cp_fm_get_element(mat_u, nr, nc, xu)
1163 6360 : CALL cp_fm_set_element(mat_w, mr, nc, xu)
1164 15144 : nc = nc + 1
1165 : END DO
1166 : END DO
1167 552 : nr = nr + 1
1168 : END DO
1169 : END DO
1170 :
1171 : ! Transformation to original non-orthogonal basis set
1172 40 : CALL parallel_gemm("N", "N", n_ao, n_mo, n_ao, 1.0_dp, ec%m_transf, mat_w, 0.0_dp, mo%mo_coeff)
1173 40 : CALL normalize_mo_vectors(qs_env, mo, n_ao, n_mo)
1174 :
1175 : ! Clean memory
1176 40 : CALL cp_fm_release(matrix=mat_w)
1177 :
1178 40 : END SUBROUTINE set_mo_coefficients
1179 :
1180 : ! **************************************************************************************************
1181 : !> \brief Creates MO set corresponding to one atomic data block
1182 : !> \param qs_env QuickStep environment containing all system data
1183 : !> \param ec electronic coupling data structure
1184 : !> \param id block ID
1185 : !> \param spin spin component
1186 : !> \param mat_u matrix of the block states
1187 : !> \param vec_e array of the block eigenvalues
1188 : !> \author Z. Futera (02.2017)
1189 : ! **************************************************************************************************
1190 40 : SUBROUTINE create_block_mo_set(qs_env, ec, id, spin, mat_u, vec_e)
1191 :
1192 : ! Routine arguments
1193 : TYPE(qs_environment_type), POINTER :: qs_env
1194 : TYPE(et_cpl), POINTER :: ec
1195 : INTEGER, INTENT(IN) :: id, spin
1196 : TYPE(cp_fm_type), INTENT(IN) :: mat_u
1197 : REAL(KIND=dp), DIMENSION(:), POINTER :: vec_e
1198 :
1199 : INTEGER :: n_ao, n_el, n_mo
1200 : REAL(KIND=dp) :: mx_occ
1201 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
1202 : TYPE(cp_fm_struct_type), POINTER :: fm_s
1203 : TYPE(dft_control_type), POINTER :: dft_cntrl
1204 : TYPE(mo_set_type), POINTER :: mo
1205 : TYPE(mp_para_env_type), POINTER :: para_env
1206 40 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1207 : TYPE(scf_control_type), POINTER :: scf_cntrl
1208 :
1209 : ! Routine name for debug purposes
1210 :
1211 40 : NULLIFY (blacs_env)
1212 40 : NULLIFY (dft_cntrl)
1213 40 : NULLIFY (para_env)
1214 40 : NULLIFY (qs_kind_set)
1215 40 : NULLIFY (fm_s)
1216 40 : NULLIFY (scf_cntrl)
1217 40 : NULLIFY (mo)
1218 :
1219 : ! Number of basis functions
1220 40 : CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set)
1221 40 : CALL get_qs_kind_set(qs_kind_set, nsgf=n_ao)
1222 :
1223 : ! Number of states
1224 40 : n_mo = mat_u%matrix_struct%nrow_global
1225 40 : IF (n_mo /= mat_u%matrix_struct%ncol_global) &
1226 0 : CPABORT('block state matrix is not square')
1227 40 : IF (n_mo /= SIZE(vec_e)) &
1228 0 : CPABORT('inconsistent number of states / energies')
1229 :
1230 : ! Maximal occupancy
1231 40 : CALL get_qs_env(qs_env, dft_control=dft_cntrl)
1232 40 : mx_occ = 2.0_dp
1233 40 : IF (dft_cntrl%nspins > 1) &
1234 40 : mx_occ = 1.0_dp
1235 :
1236 : ! Number of electrons
1237 40 : n_el = ec%block(id)%n_electrons
1238 40 : IF (dft_cntrl%nspins > 1) THEN
1239 40 : n_el = n_el/2
1240 40 : IF (MOD(ec%block(id)%n_electrons, 2) == 1) THEN
1241 28 : IF (spin == 1) &
1242 14 : n_el = n_el + 1
1243 : END IF
1244 : END IF
1245 :
1246 : ! Memory allocation (Use deallocate_mo_set to prevent accidental memory leaks)
1247 40 : CALL deallocate_mo_set(ec%block(id)%mo(spin))
1248 40 : CALL allocate_mo_set(ec%block(id)%mo(spin), n_ao, n_mo, n_el, REAL(n_el, dp), mx_occ, 0.0_dp)
1249 40 : mo => ec%block(id)%mo(spin)
1250 :
1251 : ! State energies
1252 120 : ALLOCATE (mo%eigenvalues(n_mo))
1253 40 : CPASSERT(ASSOCIATED(mo%eigenvalues))
1254 896 : mo%eigenvalues = vec_e
1255 :
1256 : ! States coefficients
1257 40 : CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
1258 : CALL cp_fm_struct_create(fmstruct=fm_s, para_env=para_env, context=blacs_env, &
1259 40 : nrow_global=n_ao, ncol_global=n_mo)
1260 40 : ALLOCATE (mo%mo_coeff)
1261 40 : CALL cp_fm_create(matrix=mo%mo_coeff, matrix_struct=fm_s, name='BLOCK STATES')
1262 :
1263 : ! Transform MO coefficients to original non-orthogonal basis set
1264 40 : CALL set_mo_coefficients(qs_env, ec, id, mo, mat_u, n_ao, n_mo)
1265 :
1266 : ! Occupancies
1267 120 : ALLOCATE (mo%occupation_numbers(n_mo))
1268 40 : CPASSERT(ASSOCIATED(mo%occupation_numbers))
1269 448 : mo%occupation_numbers = 0.0_dp
1270 :
1271 40 : IF (n_el > 0) THEN
1272 28 : CALL get_qs_env(qs_env, scf_control=scf_cntrl)
1273 28 : CALL set_mo_occupation(mo_set=mo, smear=scf_cntrl%smear)
1274 : END IF
1275 :
1276 : ! Clean memory
1277 40 : CALL cp_fm_struct_release(fmstruct=fm_s)
1278 :
1279 40 : END SUBROUTINE create_block_mo_set
1280 :
1281 : ! **************************************************************************************************
1282 : !> \brief save given electronic state to cube files
1283 : !> \param qs_env QuickStep environment containing all system data
1284 : !> \param logger output logger
1285 : !> \param input input-file block print setting section
1286 : !> \param mo electronic states data
1287 : !> \param ib block ID
1288 : !> \param im state ID
1289 : !> \param is spin ID
1290 : !> \author Z. Futera (02.2017)
1291 : ! **************************************************************************************************
1292 0 : SUBROUTINE save_mo_cube(qs_env, logger, input, mo, ib, im, is)
1293 :
1294 : ! Routine arguments
1295 : TYPE(qs_environment_type), POINTER :: qs_env
1296 : TYPE(cp_logger_type), POINTER :: logger
1297 : TYPE(section_vals_type), POINTER :: input
1298 : TYPE(mo_set_type), POINTER :: mo
1299 : INTEGER, INTENT(IN) :: ib, im, is
1300 :
1301 : CHARACTER(LEN=default_path_length) :: filename
1302 : CHARACTER(LEN=default_string_length) :: title
1303 : INTEGER :: unit_nr
1304 0 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1305 : TYPE(cell_type), POINTER :: cell
1306 : TYPE(dft_control_type), POINTER :: dft_control
1307 : TYPE(particle_list_type), POINTER :: particles
1308 0 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1309 : TYPE(pw_c1d_gs_type) :: wf_g
1310 : TYPE(pw_env_type), POINTER :: pw_env
1311 0 : TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
1312 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1313 : TYPE(pw_r3d_rs_type) :: wf_r
1314 0 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1315 : TYPE(qs_subsys_type), POINTER :: subsys
1316 :
1317 : ! Routine name for debug purposes
1318 :
1319 : ! Initialization
1320 0 : NULLIFY (particles)
1321 0 : NULLIFY (subsys)
1322 :
1323 0 : NULLIFY (pw_env)
1324 0 : NULLIFY (pw_pools)
1325 0 : NULLIFY (auxbas_pw_pool)
1326 :
1327 0 : NULLIFY (atomic_kind_set)
1328 0 : NULLIFY (cell)
1329 0 : NULLIFY (dft_control)
1330 0 : NULLIFY (particle_set)
1331 0 : NULLIFY (qs_kind_set)
1332 :
1333 : ! Name of the cube file
1334 0 : WRITE (filename, '(A4,I1.1,A1,I5.5,A1,I1.1)') 'BWF_', ib, '_', im, '_', is
1335 : ! Open the file
1336 : unit_nr = cp_print_key_unit_nr(logger, input, 'MO_CUBES', extension='.cube', &
1337 0 : middle_name=TRIM(filename), file_position='REWIND', log_filename=.FALSE.)
1338 : ! Title of the file
1339 0 : WRITE (title, *) 'WAVEFUNCTION ', im, ' block ', ib, ' spin ', is
1340 :
1341 : ! List of all atoms
1342 0 : CALL get_qs_env(qs_env, subsys=subsys)
1343 0 : CALL qs_subsys_get(subsys, particles=particles)
1344 :
1345 : ! Grids for wavefunction
1346 0 : CALL get_qs_env(qs_env, pw_env=pw_env)
1347 0 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, pw_pools=pw_pools)
1348 0 : CALL auxbas_pw_pool%create_pw(wf_r)
1349 0 : CALL auxbas_pw_pool%create_pw(wf_g)
1350 :
1351 : ! Calculate the grid values
1352 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
1353 0 : cell=cell, dft_control=dft_control, particle_set=particle_set)
1354 : CALL calculate_wavefunction(mo%mo_coeff, im, wf_r, wf_g, atomic_kind_set, &
1355 0 : qs_kind_set, cell, dft_control, particle_set, pw_env)
1356 : CALL cp_pw_to_cube(wf_r, unit_nr, title, particles=particles, &
1357 0 : stride=section_get_ivals(input, 'MO_CUBES%STRIDE'))
1358 :
1359 : ! Close file
1360 0 : CALL cp_print_key_finished_output(unit_nr, logger, input, 'MO_CUBES')
1361 :
1362 : ! Clean memory
1363 0 : CALL auxbas_pw_pool%give_back_pw(wf_r)
1364 0 : CALL auxbas_pw_pool%give_back_pw(wf_g)
1365 :
1366 0 : END SUBROUTINE save_mo_cube
1367 :
1368 : ! **************************************************************************************************
1369 : !> \brief save specified electronic states to cube files
1370 : !> \param qs_env QuickStep environment containing all system data
1371 : !> \param ec electronic coupling data structure
1372 : !> \param n_spins number of spin states
1373 : !> \author Z. Futera (02.2017)
1374 : ! **************************************************************************************************
1375 10 : SUBROUTINE save_el_states(qs_env, ec, n_spins)
1376 :
1377 : ! Routine arguments
1378 : TYPE(qs_environment_type), POINTER :: qs_env
1379 : TYPE(et_cpl), POINTER :: ec
1380 : INTEGER, INTENT(IN) :: n_spins
1381 :
1382 : INTEGER :: i, j, k, l, n
1383 10 : INTEGER, DIMENSION(:), POINTER :: list
1384 : TYPE(cp_logger_type), POINTER :: logger
1385 : TYPE(mo_set_type), POINTER :: mo
1386 : TYPE(section_vals_type), POINTER :: block_sec, mo_sec, print_sec
1387 :
1388 : ! Routine name for debug purposes
1389 :
1390 10 : NULLIFY (logger)
1391 10 : NULLIFY (block_sec)
1392 10 : NULLIFY (print_sec)
1393 10 : NULLIFY (mo_sec)
1394 :
1395 : ! Output logger
1396 20 : logger => cp_get_default_logger()
1397 : block_sec => section_vals_get_subs_vals(qs_env%input, &
1398 10 : 'PROPERTIES%ET_COUPLING%PROJECTION%BLOCK')
1399 :
1400 : ! Print states of all blocks
1401 30 : DO i = 1, ec%n_blocks
1402 :
1403 20 : print_sec => section_vals_get_subs_vals(block_sec, 'PRINT', i_rep_section=i)
1404 :
1405 : ! Check if the print input section is active
1406 20 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
1407 10 : print_sec, 'MO_CUBES'), cp_p_file)) THEN
1408 :
1409 0 : mo_sec => section_vals_get_subs_vals(print_sec, 'MO_CUBES')
1410 :
1411 : ! Spin states
1412 0 : DO j = 1, n_spins
1413 :
1414 0 : mo => ec%block(i)%mo(j)
1415 :
1416 0 : CALL section_vals_val_get(mo_sec, keyword_name='MO_LIST', n_rep_val=n)
1417 :
1418 : ! List of specific MOs
1419 0 : IF (n > 0) THEN
1420 :
1421 0 : DO k = 1, n
1422 0 : NULLIFY (list)
1423 : CALL section_vals_val_get(mo_sec, keyword_name='MO_LIST', &
1424 0 : i_rep_val=k, i_vals=list)
1425 0 : IF (ASSOCIATED(list)) THEN
1426 0 : DO l = 1, SIZE(list)
1427 0 : CALL save_mo_cube(qs_env, logger, print_sec, mo, i, list(l), j)
1428 : END DO
1429 : END IF
1430 : END DO
1431 :
1432 : ! Frontier MOs
1433 : ELSE
1434 :
1435 : ! Occupied states
1436 0 : CALL section_vals_val_get(mo_sec, keyword_name='NHOMO', i_val=n)
1437 :
1438 0 : IF (n > 0) THEN
1439 0 : DO k = MAX(1, mo%homo - n + 1), mo%homo
1440 0 : CALL save_mo_cube(qs_env, logger, print_sec, mo, i, k, j)
1441 : END DO
1442 : END IF
1443 :
1444 : ! Unoccupied states
1445 0 : CALL section_vals_val_get(mo_sec, keyword_name='NLUMO', i_val=n)
1446 :
1447 0 : IF (n > 0) THEN
1448 0 : DO k = mo%lfomo, MIN(mo%lfomo + n - 1, mo%nmo)
1449 0 : CALL save_mo_cube(qs_env, logger, print_sec, mo, i, k, j)
1450 : END DO
1451 : END IF
1452 :
1453 : END IF
1454 :
1455 : END DO
1456 :
1457 : END IF
1458 :
1459 : END DO
1460 :
1461 10 : END SUBROUTINE save_el_states
1462 :
1463 : ! **************************************************************************************************
1464 : !> \brief calculates the electron transfer coupling elements by projection-operator approach
1465 : !> Kondov et al. J.Phys.Chem.C 2007, 111, 11970-11981
1466 : !> \param qs_env QuickStep environment containing all system data
1467 : !> \author Z. Futera (02.2017)
1468 : ! **************************************************************************************************
1469 10 : SUBROUTINE calc_et_coupling_proj(qs_env)
1470 :
1471 : ! Routine arguments
1472 : TYPE(qs_environment_type), POINTER :: qs_env
1473 :
1474 : INTEGER :: i, j, k, n_ao, n_atoms, output_unit
1475 : LOGICAL :: do_kp, master
1476 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
1477 : TYPE(cp_fm_struct_type), POINTER :: fm_s
1478 : TYPE(cp_fm_type) :: mat_w
1479 10 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: mat_h
1480 : TYPE(cp_logger_type), POINTER :: logger
1481 10 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks, mo_der
1482 : TYPE(dft_control_type), POINTER :: dft_cntrl
1483 : TYPE(et_cpl), POINTER :: ec
1484 : TYPE(kpoint_type), POINTER :: kpoints
1485 10 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mo
1486 : TYPE(mp_para_env_type), POINTER :: para_env
1487 10 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1488 : TYPE(scf_control_type), POINTER :: scf_control
1489 : TYPE(section_vals_type), POINTER :: et_proj_sec
1490 :
1491 : ! Routine name for debug purposes
1492 :
1493 : ! Pointer initialization
1494 10 : NULLIFY (logger)
1495 :
1496 10 : NULLIFY (blacs_env)
1497 10 : NULLIFY (para_env)
1498 10 : NULLIFY (dft_cntrl)
1499 10 : NULLIFY (kpoints)
1500 10 : NULLIFY (qs_kind_set)
1501 : NULLIFY (et_proj_sec)
1502 :
1503 10 : NULLIFY (fm_s)
1504 10 : NULLIFY (ks, mo_der)
1505 :
1506 10 : NULLIFY (ec)
1507 :
1508 : ! Reference
1509 10 : CALL cite_reference(Futera2017)
1510 :
1511 : ! Stream for output to LOG file
1512 10 : logger => cp_get_default_logger()
1513 :
1514 10 : et_proj_sec => section_vals_get_subs_vals(qs_env%input, 'PROPERTIES%ET_COUPLING%PROJECTION')
1515 :
1516 : output_unit = cp_print_key_unit_nr(logger, et_proj_sec, &
1517 10 : 'PROGRAM_RUN_INFO', extension='.log')
1518 :
1519 : ! Parallel calculation - master thread
1520 10 : master = .FALSE.
1521 10 : IF (output_unit > 0) &
1522 : master = .TRUE.
1523 :
1524 : ! Header
1525 : IF (master) THEN
1526 : WRITE (output_unit, '(/,T2,A)') &
1527 5 : '!-----------------------------------------------------------------------------!'
1528 : WRITE (output_unit, '(T17,A)') &
1529 5 : 'Electronic coupling - Projection-operator method'
1530 : END IF
1531 :
1532 : ! Main data structure
1533 10 : ALLOCATE (ec)
1534 10 : CPASSERT(ASSOCIATED(ec))
1535 10 : CALL set_block_data(qs_env, et_proj_sec, ec)
1536 :
1537 : ! Number of atoms and AO functions
1538 10 : CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, natom=n_atoms)
1539 10 : CALL get_qs_kind_set(qs_kind_set, nsgf=n_ao)
1540 :
1541 : ! Print out info about system partitioning
1542 10 : IF (master) THEN
1543 :
1544 : WRITE (output_unit, '(/,T3,A,I10)') &
1545 5 : 'Number of atoms = ', n_atoms
1546 : WRITE (output_unit, '(T3,A,I10)') &
1547 5 : 'Number of fragments = ', ec%n_blocks
1548 : WRITE (output_unit, '(T3,A,I10)') &
1549 5 : 'Number of fragment atoms = ', ec%n_atoms
1550 : WRITE (output_unit, '(T3,A,I10)') &
1551 5 : 'Number of unassigned atoms = ', n_atoms - ec%n_atoms
1552 : WRITE (output_unit, '(T3,A,I10)') &
1553 5 : 'Number of AO basis functions = ', n_ao
1554 :
1555 15 : DO i = 1, ec%n_blocks
1556 :
1557 : WRITE (output_unit, '(/,T3,A,I0,A)') &
1558 10 : 'Block ', i, ':'
1559 : WRITE (output_unit, '(T3,A,I10)') &
1560 10 : 'Number of block atoms = ', ec%block(i)%n_atoms
1561 : WRITE (output_unit, '(T3,A,I10)') &
1562 10 : 'Number of block electrons = ', ec%block(i)%n_electrons
1563 : WRITE (output_unit, '(T3,A,I10)') &
1564 10 : 'Number of block AO functions = ', ec%block(i)%n_ao
1565 :
1566 15 : IF (ec%block(i)%n_atoms < 10) THEN
1567 :
1568 : WRITE (output_unit, '(T3,A,10I6)') &
1569 10 : 'Block atom IDs = ', &
1570 56 : (ec%block(i)%atom(j)%id, j=1, ec%block(i)%n_atoms)
1571 :
1572 : ELSE
1573 :
1574 0 : WRITE (output_unit, '(T3,A)') 'Block atom IDs ='
1575 0 : DO j = 1, ec%block(i)%n_atoms/10
1576 0 : WRITE (output_unit, '(T3,A,10I6)') ' ', &
1577 0 : (ec%block(i)%atom((j - 1)*10 + k)%id, k=1, 10)
1578 : END DO
1579 0 : IF (MOD(ec%block(i)%n_atoms, 10) /= 0) THEN
1580 0 : WRITE (output_unit, '(T3,A,10I6)') ' ', &
1581 0 : (ec%block(i)%atom(k + 10*(ec%block(i)%n_atoms/10))%id, &
1582 0 : k=1, MOD(ec%block(i)%n_atoms, 10))
1583 : END IF
1584 :
1585 : END IF
1586 :
1587 : END DO
1588 :
1589 : END IF
1590 :
1591 : ! Full matrix data structure
1592 10 : CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
1593 : CALL cp_fm_struct_create(fmstruct=fm_s, para_env=para_env, context=blacs_env, &
1594 10 : nrow_global=n_ao, ncol_global=n_ao)
1595 10 : CALL cp_fm_create(matrix=mat_w, matrix_struct=fm_s, name='FULL WORK MATRIX')
1596 :
1597 : ! Spin polarization / K-point sampling
1598 10 : CALL get_qs_env(qs_env, dft_control=dft_cntrl, do_kpoints=do_kp)
1599 10 : CALL get_qs_env(qs_env, mos=mo, matrix_ks=ks, mo_derivs=mo_der, scf_control=scf_control)
1600 10 : CALL make_mo_eig(mo, dft_cntrl%nspins, ks, scf_control, mo_der)
1601 :
1602 10 : IF (do_kp) THEN
1603 0 : CPABORT('ET_COUPLING not implemented with kpoints')
1604 : ELSE
1605 : ! no K-points
1606 10 : IF (master) &
1607 5 : WRITE (output_unit, '(T3,A)') 'No K-point sampling (Gamma point only)'
1608 : END IF
1609 :
1610 10 : IF (dft_cntrl%nspins == 2) THEN
1611 :
1612 10 : IF (master) &
1613 5 : WRITE (output_unit, '(/,T3,A)') 'Spin-polarized calculation'
1614 :
1615 : !<--- Open shell / No K-points ------------------------------------------------>!
1616 :
1617 : ! State eneries of the whole system
1618 10 : IF (mo(1)%nao /= mo(2)%nao) &
1619 0 : CPABORT('different number of alpha/beta AO basis functions')
1620 10 : IF (master) THEN
1621 : WRITE (output_unit, '(/,T3,A,I10)') &
1622 5 : 'Number of AO basis functions = ', mo(1)%nao
1623 : WRITE (output_unit, '(T3,A,I10)') &
1624 5 : 'Number of alpha states = ', mo(1)%nmo
1625 : WRITE (output_unit, '(T3,A,I10)') &
1626 5 : 'Number of beta states = ', mo(2)%nmo
1627 : END IF
1628 10 : CALL print_states(output_unit, mo, dft_cntrl%nspins, 'the whole system', fermi=.TRUE.)
1629 10 : CALL set_fermi(ec, mo(1)%mu, mo(2)%mu)
1630 :
1631 : ! KS Hamiltonian
1632 10 : CALL get_block_hamiltonian(qs_env, ec, fm_s, mat_h, mat_w, n_ao, dft_cntrl%nspins)
1633 :
1634 : ! Block diagonization
1635 10 : CALL hamiltonian_block_diag(qs_env, ec, mat_h)
1636 :
1637 : ! Print out energies and couplings
1638 30 : DO i = 1, ec%n_blocks
1639 20 : IF (output_unit > 0) THEN
1640 : CALL print_states(output_unit, ec%block(i)%mo, dft_cntrl%nspins, &
1641 : 'block '//TRIM(ADJUSTL(cp_to_string(i)))//' states', &
1642 10 : mx_mo_a=mo(1)%nmo, mx_mo_b=mo(2)%nmo, fermi=.TRUE.)
1643 : END IF
1644 30 : CALL print_mo_coeff(output_unit, qs_env, ec, i, dft_cntrl%nspins)
1645 : END DO
1646 :
1647 10 : CALL print_couplings(et_proj_sec, output_unit, logger, ec, mo)
1648 :
1649 : ELSE
1650 :
1651 0 : IF (master) &
1652 0 : WRITE (output_unit, '(/,T3,A)') 'Spin-restricted calculation'
1653 :
1654 : !<--- Close shell / No K-points ----------------------------------------------->!
1655 :
1656 : ! State eneries of the whole system
1657 : IF (master) THEN
1658 : WRITE (output_unit, '(/,T3,A,I10)') &
1659 0 : 'Number of AO basis functions = ', mo(1)%nao
1660 : WRITE (output_unit, '(T3,A,I10)') &
1661 0 : 'Number of states = ', mo(1)%nmo
1662 : END IF
1663 0 : CALL print_states(output_unit, mo, dft_cntrl%nspins, 'the whole system', fermi=.TRUE.)
1664 0 : CALL set_fermi(ec, mo(1)%mu)
1665 :
1666 : ! KS Hamiltonian
1667 0 : CALL get_block_hamiltonian(qs_env, ec, fm_s, mat_h, mat_w, n_ao, dft_cntrl%nspins)
1668 :
1669 : ! Block diagonization
1670 0 : CALL hamiltonian_block_diag(qs_env, ec, mat_h)
1671 :
1672 : ! Print out energies and couplings
1673 0 : DO i = 1, ec%n_blocks
1674 0 : IF (output_unit > 0) THEN
1675 : CALL print_states(output_unit, ec%block(i)%mo, dft_cntrl%nspins, &
1676 : 'block '//TRIM(ADJUSTL(cp_to_string(i)))//' states', &
1677 0 : mx_mo_a=mo(1)%nmo, fermi=.TRUE.)
1678 : END IF
1679 0 : CALL print_mo_coeff(output_unit, qs_env, ec, i, dft_cntrl%nspins)
1680 : END DO
1681 :
1682 0 : CALL print_couplings(et_proj_sec, output_unit, logger, ec, mo)
1683 :
1684 : END IF
1685 :
1686 : ! Save electronic states
1687 10 : CALL save_el_states(qs_env, ec, dft_cntrl%nspins)
1688 :
1689 : ! Footer
1690 10 : IF (master) WRITE (output_unit, '(/,T2,A)') &
1691 5 : '!-----------------------------------------------------------------------------!'
1692 :
1693 : ! Clean memory
1694 10 : CALL cp_fm_struct_release(fmstruct=fm_s)
1695 10 : CALL cp_fm_release(matrix=mat_w)
1696 10 : IF (ALLOCATED(mat_h)) THEN
1697 30 : DO i = 1, SIZE(mat_h)
1698 30 : CALL cp_fm_release(matrix=mat_h(i))
1699 : END DO
1700 10 : DEALLOCATE (mat_h)
1701 : END IF
1702 10 : CALL release_ec_data(ec)
1703 :
1704 : ! Close output stream
1705 10 : CALL cp_print_key_finished_output(output_unit, logger, et_proj_sec, 'PROGRAM_RUN_INFO')
1706 :
1707 20 : END SUBROUTINE calc_et_coupling_proj
1708 :
1709 0 : END MODULE et_coupling_proj
|
