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