(git:34ef472)
et_coupling_proj.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 ! **************************************************************************************************
9 !> \brief 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_operations, ONLY: copy_dbcsr_to_fm,&
25  cp_dbcsr_sm_fm_multiply
26  USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale
27  USE cp_fm_diag, ONLY: choose_eigv_solver,&
28  cp_fm_power
29  USE cp_fm_struct, ONLY: cp_fm_struct_create,&
30  cp_fm_struct_equivalent,&
31  cp_fm_struct_release,&
32  cp_fm_struct_type
33  USE cp_fm_types, ONLY: &
34  cp_fm_create, cp_fm_get_element, cp_fm_get_submatrix, cp_fm_release, cp_fm_set_all, &
35  cp_fm_set_element, cp_fm_to_fm, cp_fm_to_fm_submat, cp_fm_type, cp_fm_vectorssum
36  USE cp_log_handling, ONLY: cp_get_default_logger,&
37  cp_logger_type,&
38  cp_to_string
39  USE cp_output_handling, ONLY: cp_p_file,&
40  cp_print_key_finished_output,&
41  cp_print_key_should_output,&
42  cp_print_key_unit_nr
43  USE cp_realspace_grid_cube, ONLY: cp_pw_to_cube
44  USE dbcsr_api, ONLY: dbcsr_p_type
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  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  IF (ASSOCIATED(ec)) THEN
152 
153  IF (ASSOCIATED(ec%fermi)) &
154  DEALLOCATE (ec%fermi)
155  IF (ASSOCIATED(ec%m_transf)) THEN
156  CALL cp_fm_release(matrix=ec%m_transf)
157  DEALLOCATE (ec%m_transf)
158  NULLIFY (ec%m_transf)
159  END IF
160  IF (ASSOCIATED(ec%m_transf_inv)) THEN
161  CALL cp_fm_release(matrix=ec%m_transf_inv)
162  DEALLOCATE (ec%m_transf_inv)
163  NULLIFY (ec%m_transf_inv)
164  END IF
165 
166  IF (ASSOCIATED(ec%block)) THEN
167 
168  DO i = 1, SIZE(ec%block)
169  IF (ASSOCIATED(ec%block(i)%atom)) &
170  DEALLOCATE (ec%block(i)%atom)
171  IF (ASSOCIATED(ec%block(i)%mo)) THEN
172  DO j = 1, SIZE(ec%block(i)%mo)
173  CALL deallocate_mo_set(ec%block(i)%mo(j))
174  END DO
175  DEALLOCATE (ec%block(i)%mo)
176  END IF
177  CALL cp_fm_release(ec%block(i)%hab)
178  END DO
179 
180  DEALLOCATE (ec%block)
181 
182  END IF
183 
184  DEALLOCATE (ec)
185 
186  END IF
187 
188  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  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  INTEGER, DIMENSION(:), POINTER :: atom_id, atom_nf, atom_ps, n_shell, t
206  INTEGER, DIMENSION(:, :), POINTER :: ang_mom_id
207  LOGICAL :: found
208  TYPE(gto_basis_set_type), POINTER :: ao_basis_set
209  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
210  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  NULLIFY (ao_basis_set)
216  NULLIFY (particle_set)
217  NULLIFY (qs_kind_set)
218  NULLIFY (n_shell)
219  NULLIFY (ang_mom_id)
220  NULLIFY (atom_nf)
221  NULLIFY (atom_id)
222  NULLIFY (block_sec)
223 
224  ! Initialization
225  ec%n_atoms = 0
226  ec%n_blocks = 0
227  NULLIFY (ec%fermi)
228  NULLIFY (ec%m_transf)
229  NULLIFY (ec%m_transf_inv)
230  NULLIFY (ec%block)
231 
232  ! Number of atoms / atom types
233  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  CALL get_qs_kind_set(qs_kind_set, nsgf=n_ao)
236 
237  ! Number of AO functions per atom
238  ALLOCATE (atom_nf(n_atoms))
239  cpassert(ASSOCIATED(atom_nf))
240 
241  atom_nf = 0
242  DO i = 1, n_atoms
243  CALL get_atomic_kind(particle_set(i)%atomic_kind, kind_number=j)
244  CALL get_qs_kind(qs_kind_set(j), basis_set=ao_basis_set)
245  IF (.NOT. ASSOCIATED(ao_basis_set)) &
246  cpabort('Unsupported basis set type. ')
247  CALL get_gto_basis_set(gto_basis_set=ao_basis_set, &
248  nset=n_set, nshell=n_shell, l=ang_mom_id)
249  DO j = 1, n_set
250  DO k = 1, n_shell(j)
251  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  n = 0
258  DO i = 1, n_atoms
259  n = n + atom_nf(i)
260  END DO
261  cpassert(n == n_ao)
262 
263  ! Atom position in AO array
264  ALLOCATE (atom_ps(n_atoms))
265  cpassert(ASSOCIATED(atom_ps))
266  atom_ps = 1
267  DO i = 1, n_atoms - 1
268  atom_ps(i + 1) = atom_ps(i) + atom_nf(i)
269  END DO
270 
271  ! Number of blocks
272  block_sec => section_vals_get_subs_vals(et_proj_sec, 'BLOCK')
273  CALL section_vals_get(block_sec, n_repetition=ec%n_blocks)
274  ALLOCATE (ec%block(ec%n_blocks))
275  cpassert(ASSOCIATED(ec%block))
276 
277  ! Block data
278  ALLOCATE (t(n_atoms))
279  cpassert(ASSOCIATED(t))
280 
281  ec%n_atoms = 0
282  DO i = 1, ec%n_blocks
283 
284  ! Initialization
285  ec%block(i)%n_atoms = 0
286  ec%block(i)%n_electrons = 0
287  ec%block(i)%n_ao = 0
288  NULLIFY (ec%block(i)%atom)
289  NULLIFY (ec%block(i)%mo)
290  NULLIFY (ec%block(i)%hab)
291 
292  ! Number of electrons
293  CALL section_vals_val_get(block_sec, i_rep_section=i, &
294  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  keyword_name='ATOMS', i_vals=atom_id)
299 
300  ! Count unique atoms
301  DO j = 1, SIZE(atom_id)
302  ! Check atom ID validity
303  IF (atom_id(j) < 1 .OR. atom_id(j) > n_atoms) &
304  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  DO k = 1, i - 1
308  DO l = 1, ec%block(k)%n_atoms
309  IF (ec%block(k)%atom(l)%id == atom_id(j)) THEN
310  cpwarn('multiple definition of atom'//trim(adjustl(cp_to_string(atom_id(j)))))
311  found = .true.
312  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  IF (.NOT. found) THEN
318  DO k = 1, ec%block(i)%n_atoms
319  IF (t(k) == atom_id(j)) THEN
320  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  IF (.NOT. found) THEN
328  ec%block(i)%n_atoms = ec%block(i)%n_atoms + 1
329  t(ec%block(i)%n_atoms) = atom_id(j)
330  END IF
331  END DO
332 
333  ! Memory allocation
334  ALLOCATE (ec%block(i)%atom(ec%block(i)%n_atoms))
335  cpassert(ASSOCIATED(ec%block(i)%atom))
336 
337  ! Save atom IDs and number of AOs
338  DO j = 1, ec%block(i)%n_atoms
339  ec%block(i)%atom(j)%id = t(j)
340  ec%block(i)%atom(j)%n_ao = atom_nf(ec%block(i)%atom(j)%id)
341  ec%block(i)%atom(j)%ao_pos = atom_ps(ec%block(i)%atom(j)%id)
342  ec%block(i)%n_ao = ec%block(i)%n_ao + ec%block(i)%atom(j)%n_ao
343  END DO
344 
345  ec%n_atoms = ec%n_atoms + ec%block(i)%n_atoms
346  END DO
347 
348  ! Clean memory
349  IF (ASSOCIATED(atom_nf)) &
350  DEALLOCATE (atom_nf)
351  IF (ASSOCIATED(atom_ps)) &
352  DEALLOCATE (atom_ps)
353  IF (ASSOCIATED(t)) &
354  DEALLOCATE (t)
355 
356  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  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  NULLIFY (ec%fermi)
375 
376  IF (PRESENT(fb)) THEN
377 
378  ALLOCATE (ec%fermi(2))
379  cpassert(ASSOCIATED(ec%fermi))
380  ec%fermi(1) = fa
381  ec%fermi(2) = fb
382 
383  ELSE
384 
385  ALLOCATE (ec%fermi(1))
386  cpassert(ASSOCIATED(ec%fermi))
387  ec%fermi(1) = fa
388 
389  END IF
390 
391  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  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  IF (.NOT. cp_fm_struct_equivalent(mat_h%matrix_struct, mat_w%matrix_struct)) &
413  cpabort('cannot reorder Hamiltonian, working-matrix structure is not equivalent')
414 
415  ! Matrix-element reordering
416  nr = 1
417  ! Rows
418  DO ir = 1, ec%n_blocks
419  DO jr = 1, ec%block(ir)%n_atoms
420  DO kr = 1, ec%block(ir)%atom(jr)%n_ao
421  ! Columns
422  nc = 1
423  DO ic = 1, ec%n_blocks
424  DO jc = 1, ec%block(ic)%n_atoms
425  DO kc = 1, ec%block(ic)%atom(jc)%n_ao
426  mr = ec%block(ir)%atom(jr)%ao_pos + kr - 1
427  mc = ec%block(ic)%atom(jc)%ao_pos + kc - 1
428  CALL cp_fm_get_element(mat_h, nr, nc, xh)
429  CALL cp_fm_set_element(mat_w, nr, nc, xh)
430  nc = nc + 1
431  END DO
432  END DO
433  END DO
434  nr = nr + 1
435  END DO
436  END DO
437  END DO
438 
439  ! Copy the reordered matrix to original data array
440  CALL cp_fm_to_fm(mat_w, mat_h)
441 
442  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  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(IN) :: mat_t, mat_i, mat_w
457 
458  INTEGER :: n_deps
459  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_s
460  TYPE(scf_control_type), POINTER :: scf_cntrl
461 
462 ! Routine name for debug purposes
463 
464  NULLIFY (mat_s)
465  NULLIFY (scf_cntrl)
466 
467  CALL get_qs_env(qs_env, matrix_s=mat_s)
468  CALL copy_dbcsr_to_fm(mat_s(1)%matrix, mat_t)
469  CALL copy_dbcsr_to_fm(mat_s(1)%matrix, mat_i)
470 
471  ! Transformation S -> S^{-1/2}
472  CALL get_qs_env(qs_env, scf_control=scf_cntrl)
473  CALL cp_fm_power(mat_t, mat_w, -0.5_dp, scf_cntrl%eps_eigval, n_deps)
474  CALL cp_fm_power(mat_i, mat_w, +0.5_dp, scf_cntrl%eps_eigval, n_deps)
475  ! Sanity check
476  IF (n_deps /= 0) THEN
477  CALL cp_warn(__location__, &
478  "Overlap matrix exhibits linear dependencies. At least some "// &
479  "eigenvalues have been quenched.")
480  END IF
481 
482  END SUBROUTINE get_s_half_inv_matrix
483 
484 ! **************************************************************************************************
485 !> \brief transform KS hamiltonian to orthogonalized block-separated basis set
486 !> \param qs_env QuickStep environment containing all system data
487 !> \param ec electronic coupling data structure
488 !> \param fm_s full-matrix structure used for allocation of KS matrices
489 !> \param mat_t storage for pointers to the transformed KS matrices
490 !> \param mat_w working matrix of the same dimension
491 !> \param n_ao total number of AO basis functions
492 !> \param n_spins number of spin components
493 !> \author Z. Futera (02.2017)
494 ! **************************************************************************************************
495  SUBROUTINE get_block_hamiltonian(qs_env, ec, fm_s, mat_t, mat_w, n_ao, n_spins)
496 
497  ! Routine arguments
498  TYPE(qs_environment_type), POINTER :: qs_env
499  TYPE(et_cpl), POINTER :: ec
500  TYPE(cp_fm_struct_type), POINTER :: fm_s
501  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:), &
502  INTENT(OUT) :: mat_t
503  TYPE(cp_fm_type), INTENT(IN) :: mat_w
504  INTEGER :: n_ao, n_spins
505 
506  INTEGER :: i
507  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_h
508 
509 ! Routine name for debug purposes
510 
511  NULLIFY (mat_h)
512 
513  ! Memory allocation
514  ALLOCATE (mat_t(n_spins))
515 
516  ! KS Hamiltonian
517  CALL get_qs_env(qs_env, matrix_ks=mat_h)
518  ! Transformation matrix
519  ALLOCATE (ec%m_transf, ec%m_transf_inv)
520  CALL cp_fm_create(matrix=ec%m_transf, matrix_struct=fm_s, &
521  name='S^(-1/2) TRANSFORMATION MATRIX')
522  CALL cp_fm_create(matrix=ec%m_transf_inv, matrix_struct=fm_s, &
523  name='S^(+1/2) TRANSFORMATION MATRIX')
524  CALL get_s_half_inv_matrix(qs_env, ec%m_transf, ec%m_transf_inv, mat_w)
525 
526  DO i = 1, n_spins
527 
528  ! Full-matrix format
529  CALL cp_fm_create(matrix=mat_t(i), matrix_struct=fm_s, &
530  name='KS HAMILTONIAN IN SEPARATED ORTHOGONALIZED BASIS SET')
531  CALL copy_dbcsr_to_fm(mat_h(i)%matrix, mat_t(i))
532 
533  ! Transform KS Hamiltonian to the orthogonalized AO basis set
534  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)
535  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))
536 
537  ! Reorder KS Hamiltonain elements to defined block structure
538  CALL reorder_hamiltonian_matrix(ec, mat_t(i), mat_w)
539 
540  END DO
541 
542  END SUBROUTINE get_block_hamiltonian
543 
544 ! **************************************************************************************************
545 !> \brief Diagonalize diagonal blocks of the KS hamiltonian in separated orthogonalized basis set
546 !> \param qs_env QuickStep environment containing all system data
547 !> \param ec electronic coupling data structure
548 !> \param mat_h Hamiltonian in separated orthogonalized basis set
549 !> \author Z. Futera (02.2017)
550 ! **************************************************************************************************
551  SUBROUTINE hamiltonian_block_diag(qs_env, ec, mat_h)
552 
553  ! Routine arguments
554  TYPE(qs_environment_type), POINTER :: qs_env
555  TYPE(et_cpl), POINTER :: ec
556  TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mat_h
557 
558  INTEGER :: i, j, k, l, n_spins, spin
559  REAL(KIND=dp), DIMENSION(:), POINTER :: vec_e
560  TYPE(cp_blacs_env_type), POINTER :: blacs_env
561  TYPE(cp_fm_struct_type), POINTER :: fm_s
562  TYPE(cp_fm_type) :: mat_u
563  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: dat
564  TYPE(mp_para_env_type), POINTER :: para_env
565 
566 ! Routine name for debug purposes
567 
568  NULLIFY (vec_e)
569  NULLIFY (blacs_env)
570  NULLIFY (para_env)
571  NULLIFY (fm_s)
572 
573  ! Parallel environment
574  CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
575 
576  ! Storage for block sub-matrices
577  ALLOCATE (dat(ec%n_blocks))
578  cpassert(ALLOCATED(dat))
579 
580  ! Storage for electronic states and couplings
581  n_spins = SIZE(mat_h)
582  DO i = 1, ec%n_blocks
583  ALLOCATE (ec%block(i)%mo(n_spins))
584  cpassert(ASSOCIATED(ec%block(i)%mo))
585  ALLOCATE (ec%block(i)%hab(n_spins, ec%n_blocks))
586  cpassert(ASSOCIATED(ec%block(i)%hab))
587  END DO
588 
589  ! Spin components
590  DO spin = 1, n_spins
591 
592  ! Diagonal blocks
593  j = 1
594  DO i = 1, ec%n_blocks
595 
596  ! Memory allocation
597  CALL cp_fm_struct_create(fmstruct=fm_s, para_env=para_env, context=blacs_env, &
598  nrow_global=ec%block(i)%n_ao, ncol_global=ec%block(i)%n_ao)
599  CALL cp_fm_create(matrix=dat(i), matrix_struct=fm_s, &
600  name='H_KS DIAGONAL BLOCK')
601 
602  ALLOCATE (vec_e(ec%block(i)%n_ao))
603  cpassert(ASSOCIATED(vec_e))
604 
605  ! Copy block data
606  CALL cp_fm_to_fm_submat(mat_h(spin), &
607  dat(i), ec%block(i)%n_ao, &
608  ec%block(i)%n_ao, j, j, 1, 1)
609 
610  ! Diagonalization
611  CALL cp_fm_create(matrix=mat_u, matrix_struct=fm_s, name='UNITARY MATRIX')
612  CALL choose_eigv_solver(dat(i), mat_u, vec_e)
613  CALL cp_fm_to_fm(mat_u, dat(i))
614 
615  ! Save state energies / vectors
616  CALL create_block_mo_set(qs_env, ec, i, spin, mat_u, vec_e)
617 
618  ! Clean memory
619  CALL cp_fm_struct_release(fmstruct=fm_s)
620  CALL cp_fm_release(matrix=mat_u)
621  DEALLOCATE (vec_e)
622 
623  ! Off-set for next block
624  j = j + ec%block(i)%n_ao
625 
626  END DO
627 
628  ! Off-diagonal blocks
629  k = 1
630  DO i = 1, ec%n_blocks
631  l = 1
632  DO j = 1, ec%n_blocks
633  IF (i /= j) THEN
634 
635  ! Memory allocation
636  CALL cp_fm_struct_create(fmstruct=fm_s, para_env=para_env, context=blacs_env, &
637  nrow_global=ec%block(i)%n_ao, ncol_global=ec%block(j)%n_ao)
638  CALL cp_fm_create(matrix=ec%block(i)%hab(spin, j), matrix_struct=fm_s, &
639  name='H_KS OFF-DIAGONAL BLOCK')
640 
641  ! Copy block data
642  CALL cp_fm_to_fm_submat(mat_h(spin), &
643  ec%block(i)%hab(spin, j), ec%block(i)%n_ao, &
644  ec%block(j)%n_ao, k, l, 1, 1)
645 
646  ! Transformation
647  CALL cp_fm_create(matrix=mat_u, matrix_struct=fm_s, name='FULL WORK MATRIX')
648  CALL parallel_gemm("T", "N", ec%block(i)%n_ao, ec%block(j)%n_ao, ec%block(i)%n_ao, &
649  1.0_dp, dat(i), ec%block(i)%hab(spin, j), 0.0_dp, mat_u)
650  CALL parallel_gemm("N", "N", ec%block(i)%n_ao, ec%block(j)%n_ao, ec%block(j)%n_ao, &
651  1.0_dp, mat_u, dat(j), 0.0_dp, ec%block(i)%hab(spin, j))
652 
653  ! Clean memory
654  CALL cp_fm_struct_release(fmstruct=fm_s)
655  CALL cp_fm_release(matrix=mat_u)
656 
657  END IF
658  ! Off-set for next block
659  l = l + ec%block(j)%n_ao
660  END DO
661  ! Off-set for next block
662  k = k + ec%block(i)%n_ao
663  END DO
664 
665  ! Clean memory
666  IF (ALLOCATED(dat)) THEN
667  DO i = 1, SIZE(dat)
668  CALL cp_fm_release(dat(i))
669  END DO
670  END IF
671  END DO
672 
673  ! Clean memory
674  IF (ALLOCATED(dat)) &
675  DEALLOCATE (dat)
676 
677  END SUBROUTINE hamiltonian_block_diag
678 
679 ! **************************************************************************************************
680 !> \brief Return sum of selected squared MO coefficients
681 !> \param blk_at list of atoms in the block
682 !> \param mo array of MO sets
683 !> \param id state index
684 !> \param atom list of atoms for MO coefficient summing
685 !> \return ...
686 !> \author Z. Futera (02.2017)
687 ! **************************************************************************************************
688  FUNCTION get_mo_c2_sum(blk_at, mo, id, atom) RESULT(c2)
689 
690  ! Routine arguments
691  TYPE(et_cpl_atom), DIMENSION(:), POINTER :: blk_at
692  TYPE(cp_fm_type), INTENT(IN) :: mo
693  INTEGER, INTENT(IN) :: id
694  INTEGER, DIMENSION(:), POINTER :: atom
695  REAL(KIND=dp) :: c2
696 
697  INTEGER :: i, ir, j, k
698  LOGICAL :: found
699  REAL(KIND=dp) :: c
700 
701 ! Returning value
702 ! Routine name for debug purposes
703 ! Local variables
704 
705  ! initialization
706  c2 = 0.0d0
707 
708  ! selected atoms
709  DO i = 1, SIZE(atom)
710 
711  ! find atomic function offset
712  found = .false.
713  DO j = 1, SIZE(blk_at)
714  IF (blk_at(j)%id == atom(i)) THEN
715  found = .true.
716  EXIT
717  END IF
718  END DO
719 
720  IF (.NOT. found) &
721  cpabort('MO-fraction atom ID not defined in the block')
722 
723  ! sum MO coefficients from the atom
724  DO k = 1, blk_at(j)%n_ao
725  ir = blk_at(j)%ao_pos + k - 1
726  CALL cp_fm_get_element(mo, ir, id, c)
727  c2 = c2 + c*c
728  END DO
729 
730  END DO
731 
732  END FUNCTION get_mo_c2_sum
733 
734 ! **************************************************************************************************
735 !> \brief Print out specific MO coefficients
736 !> \param output_unit unit number of the open output stream
737 !> \param qs_env QuickStep environment containing all system data
738 !> \param ec electronic coupling data structure
739 !> \param blk atomic-block ID
740 !> \param n_spins number of spin components
741 !> \author Z. Futera (02.2017)
742 ! **************************************************************************************************
743  SUBROUTINE print_mo_coeff(output_unit, qs_env, ec, blk, n_spins)
744 
745  ! Routine arguments
746  INTEGER, INTENT(IN) :: output_unit
747  TYPE(qs_environment_type), POINTER :: qs_env
748  TYPE(et_cpl), POINTER :: ec
749  INTEGER, INTENT(IN) :: blk, n_spins
750 
751  INTEGER :: j, k, l, m, n, n_ao, n_mo
752  INTEGER, DIMENSION(:), POINTER :: list_at, list_mo
753  REAL(KIND=dp) :: c1, c2
754  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: mat_w
755  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
756  TYPE(section_vals_type), POINTER :: block_sec, print_sec
757 
758 ! Routine name for debug purposes
759 
760  NULLIFY (block_sec)
761  NULLIFY (print_sec)
762  NULLIFY (qs_kind_set)
763 
764  ! Atomic block data
765  block_sec => section_vals_get_subs_vals(qs_env%input, &
766  'PROPERTIES%ET_COUPLING%PROJECTION%BLOCK')
767 
768  print_sec => section_vals_get_subs_vals(block_sec, 'PRINT', i_rep_section=blk)
769 
770  ! List of atoms
771  CALL section_vals_val_get(print_sec, keyword_name='MO_COEFF_ATOM', n_rep_val=n)
772 
773  IF (n > 0) THEN
774 
775  IF (output_unit > 0) &
776  WRITE (output_unit, '(/,T3,A/)') 'Block state fractions:'
777 
778  ! Number of AO functions
779  CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set)
780  CALL get_qs_kind_set(qs_kind_set, nsgf=n_ao)
781 
782  ! MOs in orthonormal basis set
783  ALLOCATE (mat_w(n_spins))
784  DO j = 1, n_spins
785  n_mo = ec%block(blk)%n_ao
786  CALL cp_fm_create(matrix=mat_w(j), &
787  matrix_struct=ec%block(blk)%mo(j)%mo_coeff%matrix_struct, &
788  name='BLOCK MOs IN ORTHONORMAL BASIS SET')
789  CALL parallel_gemm("N", "N", n_ao, n_mo, n_ao, 1.0_dp, ec%m_transf_inv, &
790  ec%block(blk)%mo(j)%mo_coeff, 0.0_dp, mat_w(j))
791  END DO
792 
793  DO j = 1, n
794  NULLIFY (list_at)
795  CALL section_vals_val_get(print_sec, keyword_name='MO_COEFF_ATOM', &
796  i_rep_val=j, i_vals=list_at)
797  IF (ASSOCIATED(list_at)) THEN
798 
799  ! List of states
800  CALL section_vals_val_get(print_sec, keyword_name='MO_COEFF_ATOM_STATE', n_rep_val=m)
801 
802  IF (m > 0) THEN
803 
804  DO k = 1, m
805  NULLIFY (list_mo)
806  CALL section_vals_val_get(print_sec, keyword_name='MO_COEFF_ATOM_STATE', &
807  i_rep_val=k, i_vals=list_mo)
808  IF (ASSOCIATED(list_mo)) THEN
809 
810  IF (j > 1) THEN
811  IF (output_unit > 0) &
812  WRITE (output_unit, *)
813  END IF
814 
815  DO l = 1, SIZE(list_mo)
816 
817  IF (n_spins > 1) THEN
818  c1 = get_mo_c2_sum(ec%block(blk)%atom, mat_w(1), &
819  list_mo(l), list_at)
820  c2 = get_mo_c2_sum(ec%block(blk)%atom, mat_w(2), &
821  list_mo(l), list_at)
822  IF (output_unit > 0) &
823  WRITE (output_unit, '(I5,A,I5,2F20.10)') j, ' /', list_mo(l), c1, c2
824  ELSE
825  c1 = get_mo_c2_sum(ec%block(blk)%atom, mat_w(1), &
826  list_mo(l), list_at)
827  IF (output_unit > 0) &
828  WRITE (output_unit, '(I5,A,I5,F20.10)') j, ' /', list_mo(l), c1
829  END IF
830 
831  END DO
832 
833  END IF
834  END DO
835 
836  END IF
837 
838  END IF
839  END DO
840 
841  ! Clean memory
842  CALL cp_fm_release(mat_w)
843 
844  END IF
845 
846  END SUBROUTINE print_mo_coeff
847 
848 ! **************************************************************************************************
849 !> \brief Print out electronic states (MOs)
850 !> \param output_unit unit number of the open output stream
851 !> \param mo array of MO sets
852 !> \param n_spins number of spin components
853 !> \param label output label
854 !> \param mx_mo_a maximum number of alpha states to print out
855 !> \param mx_mo_b maximum number of beta states to print out
856 !> \param fermi print out Fermi level and number of electrons
857 !> \author Z. Futera (02.2017)
858 ! **************************************************************************************************
859  SUBROUTINE print_states(output_unit, mo, n_spins, label, mx_mo_a, mx_mo_b, fermi)
860 
861  ! Routine arguments
862  INTEGER, INTENT(IN) :: output_unit
863  TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mo
864  INTEGER, INTENT(IN) :: n_spins
865  CHARACTER(LEN=*), INTENT(IN) :: label
866  INTEGER, INTENT(IN), OPTIONAL :: mx_mo_a, mx_mo_b
867  LOGICAL, INTENT(IN), OPTIONAL :: fermi
868 
869  INTEGER :: i, mx_a, mx_b, n
870  LOGICAL :: prnt_fm
871 
872 ! Routine name for debug purposes
873 
874  prnt_fm = .false.
875  IF (PRESENT(fermi)) &
876  prnt_fm = fermi
877 
878  IF (output_unit > 0) THEN
879 
880  WRITE (output_unit, '(/,T3,A/)') 'State energies ('//trim(adjustl(label))//'):'
881 
882  ! Spin-polarized calculation
883  IF (n_spins > 1) THEN
884 
885  mx_a = mo(1)%nmo
886  IF (PRESENT(mx_mo_a)) &
887  mx_a = min(mo(1)%nmo, mx_mo_a)
888  mx_b = mo(2)%nmo
889  IF (PRESENT(mx_mo_b)) &
890  mx_b = min(mo(2)%nmo, mx_mo_b)
891  n = max(mx_a, mx_b)
892 
893  DO i = 1, n
894  WRITE (output_unit, '(T3,I10)', advance='no') i
895  IF (i <= mx_a) THEN
896  WRITE (output_unit, '(2F12.4)', advance='no') &
897  mo(1)%occupation_numbers(i), mo(1)%eigenvalues(i)
898  ELSE
899  WRITE (output_unit, '(A)', advance='no') ' '
900  END IF
901  WRITE (output_unit, '(A)', advance='no') ' '
902  IF (i <= mx_b) THEN
903  WRITE (output_unit, '(2F12.4)') &
904  mo(2)%occupation_numbers(i), mo(2)%eigenvalues(i)
905  ELSE
906  WRITE (output_unit, *)
907  END IF
908  END DO
909 
910  IF (prnt_fm) THEN
911  WRITE (output_unit, '(/,T3,I10,F24.4,I10,F19.4)') &
912  mo(1)%nelectron, mo(1)%mu, &
913  mo(2)%nelectron, mo(2)%mu
914  END IF
915 
916  ! Spin-restricted calculation
917  ELSE
918 
919  mx_a = mo(1)%nmo
920  IF (PRESENT(mx_mo_a)) &
921  mx_a = min(mo(1)%nmo, mx_mo_a)
922 
923  DO i = 1, mx_a
924  WRITE (output_unit, '(T3,I10,2F12.4)') &
925  i, mo(1)%occupation_numbers(i), mo(1)%eigenvalues(i)
926  END DO
927 
928  IF (prnt_fm) THEN
929  WRITE (output_unit, '(/,T3,I10,F24.4)') &
930  mo(1)%nelectron, mo(1)%mu
931  END IF
932 
933  END IF
934 
935  END IF
936 
937  END SUBROUTINE print_states
938 
939 ! **************************************************************************************************
940 !> \brief Print out donor-acceptor state couplings
941 !> \param ec_sec ...
942 !> \param output_unit unit number of the open output stream
943 !> \param logger ...
944 !> \param ec electronic coupling data structure
945 !> \param mo ...
946 !> \author Z. Futera (02.2017)
947 ! **************************************************************************************************
948  SUBROUTINE print_couplings(ec_sec, output_unit, logger, ec, mo)
949 
950  ! Routine arguments
951  TYPE(section_vals_type), POINTER :: ec_sec
952  INTEGER, INTENT(IN) :: output_unit
953  TYPE(cp_logger_type), POINTER :: logger
954  TYPE(et_cpl), POINTER :: ec
955  TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mo
956 
957  CHARACTER(LEN=default_path_length) :: filename, my_pos, title
958  INTEGER :: i, j, k, l, n_states(2), nc, nr, nspins, &
959  unit_nr
960  LOGICAL :: append
961  REAL(KIND=dp), DIMENSION(:, :), POINTER :: w1, w2
962  TYPE(section_vals_type), POINTER :: print_key
963 
964 ! Routine name for debug purposes
965 ! Local variables
966 
967  n_states = 0
968  DO i = 1, SIZE(mo)
969  n_states(i) = mo(i)%nmo
970  END DO
971  nspins = 1
972  IF (n_states(2) > 0) nspins = 2
973 
974  print_key => section_vals_get_subs_vals(section_vals=ec_sec, &
975  subsection_name="PRINT%COUPLINGS")
976 
977  IF (btest(cp_print_key_should_output(logger%iter_info, print_key), &
978  cp_p_file)) THEN
979 
980  my_pos = "REWIND"
981  append = section_get_lval(print_key, "APPEND")
982  IF (append) THEN
983  my_pos = "APPEND"
984  END IF
985 
986  IF (output_unit > 0) &
987  WRITE (output_unit, '(/,T3,A/)') 'Printing coupling elements to output files'
988 
989  DO i = 1, ec%n_blocks
990  DO j = i + 1, ec%n_blocks
991 
992  nr = ec%block(i)%hab(1, j)%matrix_struct%nrow_global
993  nc = ec%block(i)%hab(1, j)%matrix_struct%ncol_global
994 
995  ALLOCATE (w1(nr, nc))
996  cpassert(ASSOCIATED(w1))
997  CALL cp_fm_get_submatrix(ec%block(i)%hab(1, j), w1)
998  IF (nspins > 1) THEN
999  ALLOCATE (w2(nr, nc))
1000  cpassert(ASSOCIATED(w2))
1001  CALL cp_fm_get_submatrix(ec%block(i)%hab(2, j), w2)
1002  END IF
1003 
1004  IF (output_unit > 0) THEN
1005 
1006  WRITE (filename, '(a5,I1.1,a1,I1.1)') "ET_BL_", i, "-", j
1007  unit_nr = cp_print_key_unit_nr(logger, ec_sec, "PRINT%COUPLINGS", extension=".elcoup", &
1008  middle_name=trim(filename), file_position=my_pos, log_filename=.false.)
1009 
1010  WRITE (title, *) 'Coupling elements [meV] between blocks:', i, j
1011 
1012  WRITE (unit_nr, *) trim(title)
1013  IF (nspins > 1) THEN
1014  WRITE (unit_nr, '(T3,A8,T13,A8,T28,A,A)') 'State A', 'State B', 'Coupling spin 1', ' Coupling spin 2'
1015  ELSE
1016  WRITE (unit_nr, '(T3,A8,T13,A8,T28,A)') 'State A', 'State B', 'Coupling'
1017  END IF
1018 
1019  DO k = 1, min(ec%block(i)%n_ao, n_states(1))
1020  DO l = 1, min(ec%block(j)%n_ao, n_states(1))
1021 
1022  IF (nspins > 1) THEN
1023 
1024  WRITE (unit_nr, '(T3,I5,T13,I5,T22,E20.6)', advance='no') &
1025  k, l, w1(k, l)*evolt*1000.0_dp
1026  IF ((k <= n_states(2)) .AND. (l <= n_states(2))) THEN
1027  WRITE (unit_nr, '(E20.6)') &
1028  w2(k, l)*evolt*1000.0_dp
1029  ELSE
1030  WRITE (unit_nr, *)
1031  END IF
1032 
1033  ELSE
1034 
1035  WRITE (unit_nr, '(T3,I5,T13,I5,T22,E20.6)') &
1036  k, l, w1(k, l)*evolt*1000.0_dp
1037  END IF
1038 
1039  END DO
1040  WRITE (unit_nr, *)
1041  END DO
1042  CALL cp_print_key_finished_output(unit_nr, logger, ec_sec, "PRINT%COUPLINGS")
1043 
1044  END IF
1045 
1046  IF (ASSOCIATED(w1)) DEALLOCATE (w1)
1047  IF (ASSOCIATED(w2)) DEALLOCATE (w2)
1048 
1049  END DO
1050  END DO
1051 
1052  END IF
1053  END SUBROUTINE print_couplings
1054 
1055 ! **************************************************************************************************
1056 !> \brief Normalize set of MO vectors
1057 !> \param qs_env QuickStep environment containing all system data
1058 !> \param mo storage for the MO data set
1059 !> \param n_ao number of AO basis functions
1060 !> \param n_mo number of block states
1061 !> \author Z. Futera (02.2017)
1062 ! **************************************************************************************************
1063  SUBROUTINE normalize_mo_vectors(qs_env, mo, n_ao, n_mo)
1064 
1065  ! Routine arguments
1066  TYPE(qs_environment_type), POINTER :: qs_env
1067  TYPE(mo_set_type), POINTER :: mo
1068  INTEGER, INTENT(IN) :: n_ao, n_mo
1069 
1070  REAL(KIND=dp), DIMENSION(:), POINTER :: vec_t
1071  TYPE(cp_blacs_env_type), POINTER :: blacs_env
1072  TYPE(cp_fm_struct_type), POINTER :: fm_s
1073  TYPE(cp_fm_type) :: mat_sc, mat_t
1074  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_s
1075  TYPE(mp_para_env_type), POINTER :: para_env
1076 
1077 ! Routine name for debug purposes
1078 
1079  ! Initialization
1080  NULLIFY (blacs_env)
1081  NULLIFY (para_env)
1082  NULLIFY (fm_s)
1083  NULLIFY (mat_s)
1084  NULLIFY (vec_t)
1085 
1086  ! Overlap matrix
1087  CALL get_qs_env(qs_env, matrix_s=mat_s)
1088 
1089  ! Calculate S*C product
1090  CALL cp_fm_create(matrix=mat_sc, matrix_struct=mo%mo_coeff%matrix_struct, &
1091  name='S*C PRODUCT MATRIX')
1092  CALL cp_dbcsr_sm_fm_multiply(mat_s(1)%matrix, mo%mo_coeff, mat_sc, n_mo)
1093 
1094  ! Calculate C^T*S*C
1095  CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
1096  CALL cp_fm_struct_create(fmstruct=fm_s, para_env=para_env, context=blacs_env, &
1097  nrow_global=n_mo, ncol_global=n_mo)
1098  CALL cp_fm_create(matrix=mat_t, matrix_struct=fm_s, &
1099  name='C^T*S*C OVERLAP PRODUCT MATRIX')
1100  CALL parallel_gemm('T', 'N', n_mo, n_mo, n_ao, 1.0_dp, mo%mo_coeff, mat_sc, 0.0_dp, mat_t)
1101 
1102  ! Normalization
1103  ALLOCATE (vec_t(n_mo))
1104  cpassert(ASSOCIATED(vec_t))
1105  CALL cp_fm_vectorssum(mat_t, vec_t)
1106  vec_t = 1.0_dp/dsqrt(vec_t)
1107  CALL cp_fm_column_scale(mo%mo_coeff, vec_t)
1108 
1109  ! Clean memory
1110  CALL cp_fm_struct_release(fmstruct=fm_s)
1111  CALL cp_fm_release(matrix=mat_sc)
1112  CALL cp_fm_release(matrix=mat_t)
1113  IF (ASSOCIATED(vec_t)) &
1114  DEALLOCATE (vec_t)
1115 
1116  END SUBROUTINE normalize_mo_vectors
1117 
1118 ! **************************************************************************************************
1119 !> \brief Transform block MO coefficients to original non-orthogonal basis set and save them
1120 !> \param qs_env QuickStep environment containing all system data
1121 !> \param ec electronic coupling data structure
1122 !> \param id block ID
1123 !> \param mo storage for the MO data set
1124 !> \param mat_u matrix of the block states
1125 !> \param n_ao number of AO basis functions
1126 !> \param n_mo number of block states
1127 !> \author Z. Futera (02.2017)
1128 ! **************************************************************************************************
1129  SUBROUTINE set_mo_coefficients(qs_env, ec, id, mo, mat_u, n_ao, n_mo)
1130 
1131  ! Routine arguments
1132  TYPE(qs_environment_type), POINTER :: qs_env
1133  TYPE(et_cpl), POINTER :: ec
1134  INTEGER, INTENT(IN) :: id
1135  TYPE(mo_set_type), POINTER :: mo
1136  TYPE(cp_fm_type), INTENT(IN) :: mat_u
1137  INTEGER, INTENT(IN) :: n_ao, n_mo
1138 
1139  INTEGER :: ic, ir, jc, jr, mr, nc, nr
1140  REAL(KIND=dp) :: xu
1141  TYPE(cp_fm_type) :: mat_w
1142 
1143 ! Routine name for debug purposes
1144 ! Local variables
1145 
1146  ! Working matrix
1147  CALL cp_fm_create(matrix=mat_w, matrix_struct=mo%mo_coeff%matrix_struct, &
1148  name='BLOCK MO-TRANSFORMATION WORKING MATRIX')
1149  CALL cp_fm_set_all(mat_w, 0.0_dp)
1150 
1151  ! Matrix-element reordering
1152  nr = 1
1153  ! Rows
1154  DO ir = 1, ec%block(id)%n_atoms
1155  DO jr = 1, ec%block(id)%atom(ir)%n_ao
1156  ! Columns
1157  nc = 1
1158  DO ic = 1, ec%block(id)%n_atoms
1159  DO jc = 1, ec%block(id)%atom(ic)%n_ao
1160  mr = ec%block(id)%atom(ir)%ao_pos + jr - 1
1161  CALL cp_fm_get_element(mat_u, nr, nc, xu)
1162  CALL cp_fm_set_element(mat_w, mr, nc, xu)
1163  nc = nc + 1
1164  END DO
1165  END DO
1166  nr = nr + 1
1167  END DO
1168  END DO
1169 
1170  ! Transformation to original non-orthogonal basis set
1171  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)
1172  CALL normalize_mo_vectors(qs_env, mo, n_ao, n_mo)
1173 
1174  ! Clean memory
1175  CALL cp_fm_release(matrix=mat_w)
1176 
1177  END SUBROUTINE set_mo_coefficients
1178 
1179 ! **************************************************************************************************
1180 !> \brief Creates MO set corresponding to one atomic data block
1181 !> \param qs_env QuickStep environment containing all system data
1182 !> \param ec electronic coupling data structure
1183 !> \param id block ID
1184 !> \param spin spin component
1185 !> \param mat_u matrix of the block states
1186 !> \param vec_e array of the block eigenvalues
1187 !> \author Z. Futera (02.2017)
1188 ! **************************************************************************************************
1189  SUBROUTINE create_block_mo_set(qs_env, ec, id, spin, mat_u, vec_e)
1190 
1191  ! Routine arguments
1192  TYPE(qs_environment_type), POINTER :: qs_env
1193  TYPE(et_cpl), POINTER :: ec
1194  INTEGER, INTENT(IN) :: id, spin
1195  TYPE(cp_fm_type), INTENT(IN) :: mat_u
1196  REAL(KIND=dp), DIMENSION(:), POINTER :: vec_e
1197 
1198  INTEGER :: n_ao, n_el, n_mo
1199  REAL(KIND=dp) :: mx_occ
1200  TYPE(cp_blacs_env_type), POINTER :: blacs_env
1201  TYPE(cp_fm_struct_type), POINTER :: fm_s
1202  TYPE(dft_control_type), POINTER :: dft_cntrl
1203  TYPE(mo_set_type), POINTER :: mo
1204  TYPE(mp_para_env_type), POINTER :: para_env
1205  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1206  TYPE(scf_control_type), POINTER :: scf_cntrl
1207 
1208 ! Routine name for debug purposes
1209 
1210  NULLIFY (blacs_env)
1211  NULLIFY (dft_cntrl)
1212  NULLIFY (para_env)
1213  NULLIFY (qs_kind_set)
1214  NULLIFY (fm_s)
1215  NULLIFY (scf_cntrl)
1216  NULLIFY (mo)
1217 
1218  ! Number of basis functions
1219  CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set)
1220  CALL get_qs_kind_set(qs_kind_set, nsgf=n_ao)
1221 
1222  ! Number of states
1223  n_mo = mat_u%matrix_struct%nrow_global
1224  IF (n_mo /= mat_u%matrix_struct%ncol_global) &
1225  cpabort('block state matrix is not square')
1226  IF (n_mo /= SIZE(vec_e)) &
1227  cpabort('inconsistent number of states / energies')
1228 
1229  ! Maximal occupancy
1230  CALL get_qs_env(qs_env, dft_control=dft_cntrl)
1231  mx_occ = 2.0_dp
1232  IF (dft_cntrl%nspins > 1) &
1233  mx_occ = 1.0_dp
1234 
1235  ! Number of electrons
1236  n_el = ec%block(id)%n_electrons
1237  IF (dft_cntrl%nspins > 1) THEN
1238  n_el = n_el/2
1239  IF (mod(ec%block(id)%n_electrons, 2) == 1) THEN
1240  IF (spin == 1) &
1241  n_el = n_el + 1
1242  END IF
1243  END IF
1244 
1245  ! Memory allocation (Use deallocate_mo_set to prevent accidental memory leaks)
1246  CALL deallocate_mo_set(ec%block(id)%mo(spin))
1247  CALL allocate_mo_set(ec%block(id)%mo(spin), n_ao, n_mo, n_el, real(n_el, dp), mx_occ, 0.0_dp)
1248  mo => ec%block(id)%mo(spin)
1249 
1250  ! State energies
1251  ALLOCATE (mo%eigenvalues(n_mo))
1252  cpassert(ASSOCIATED(mo%eigenvalues))
1253  mo%eigenvalues = vec_e
1254 
1255  ! States coefficients
1256  CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
1257  CALL cp_fm_struct_create(fmstruct=fm_s, para_env=para_env, context=blacs_env, &
1258  nrow_global=n_ao, ncol_global=n_mo)
1259  ALLOCATE (mo%mo_coeff)
1260  CALL cp_fm_create(matrix=mo%mo_coeff, matrix_struct=fm_s, name='BLOCK STATES')
1261 
1262  ! Transform MO coefficients to original non-orthogonal basis set
1263  CALL set_mo_coefficients(qs_env, ec, id, mo, mat_u, n_ao, n_mo)
1264 
1265  ! Occupancies
1266  ALLOCATE (mo%occupation_numbers(n_mo))
1267  cpassert(ASSOCIATED(mo%occupation_numbers))
1268  mo%occupation_numbers = 0.0_dp
1269 
1270  IF (n_el > 0) THEN
1271  CALL get_qs_env(qs_env, scf_control=scf_cntrl)
1272  CALL set_mo_occupation(mo_set=mo, smear=scf_cntrl%smear)
1273  END IF
1274 
1275  ! Clean memory
1276  CALL cp_fm_struct_release(fmstruct=fm_s)
1277 
1278  END SUBROUTINE create_block_mo_set
1279 
1280 ! **************************************************************************************************
1281 !> \brief save given electronic state to cube files
1282 !> \param qs_env QuickStep environment containing all system data
1283 !> \param logger output logger
1284 !> \param input input-file block print setting section
1285 !> \param mo electronic states data
1286 !> \param ib block ID
1287 !> \param im state ID
1288 !> \param is spin ID
1289 !> \author Z. Futera (02.2017)
1290 ! **************************************************************************************************
1291  SUBROUTINE save_mo_cube(qs_env, logger, input, mo, ib, im, is)
1292 
1293  ! Routine arguments
1294  TYPE(qs_environment_type), POINTER :: qs_env
1295  TYPE(cp_logger_type), POINTER :: logger
1296  TYPE(section_vals_type), POINTER :: input
1297  TYPE(mo_set_type), POINTER :: mo
1298  INTEGER, INTENT(IN) :: ib, im, is
1299 
1300  CHARACTER(LEN=default_path_length) :: filename
1301  CHARACTER(LEN=default_string_length) :: title
1302  INTEGER :: unit_nr
1303  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1304  TYPE(cell_type), POINTER :: cell
1305  TYPE(dft_control_type), POINTER :: dft_control
1306  TYPE(particle_list_type), POINTER :: particles
1307  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1308  TYPE(pw_c1d_gs_type) :: wf_g
1309  TYPE(pw_env_type), POINTER :: pw_env
1310  TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
1311  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1312  TYPE(pw_r3d_rs_type) :: wf_r
1313  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1314  TYPE(qs_subsys_type), POINTER :: subsys
1315 
1316 ! Routine name for debug purposes
1317 
1318  ! Initialization
1319  NULLIFY (particles)
1320  NULLIFY (subsys)
1321 
1322  NULLIFY (pw_env)
1323  NULLIFY (pw_pools)
1324  NULLIFY (auxbas_pw_pool)
1325 
1326  NULLIFY (atomic_kind_set)
1327  NULLIFY (cell)
1328  NULLIFY (dft_control)
1329  NULLIFY (particle_set)
1330  NULLIFY (qs_kind_set)
1331 
1332  ! Name of the cube file
1333  WRITE (filename, '(A4,I1.1,A1,I5.5,A1,I1.1)') 'BWF_', ib, '_', im, '_', is
1334  ! Open the file
1335  unit_nr = cp_print_key_unit_nr(logger, input, 'MO_CUBES', extension='.cube', &
1336  middle_name=trim(filename), file_position='REWIND', log_filename=.false.)
1337  ! Title of the file
1338  WRITE (title, *) 'WAVEFUNCTION ', im, ' block ', ib, ' spin ', is
1339 
1340  ! List of all atoms
1341  CALL get_qs_env(qs_env, subsys=subsys)
1342  CALL qs_subsys_get(subsys, particles=particles)
1343 
1344  ! Grids for wavefunction
1345  CALL get_qs_env(qs_env, pw_env=pw_env)
1346  CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, pw_pools=pw_pools)
1347  CALL auxbas_pw_pool%create_pw(wf_r)
1348  CALL auxbas_pw_pool%create_pw(wf_g)
1349 
1350  ! Calculate the grid values
1351  CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
1352  cell=cell, dft_control=dft_control, particle_set=particle_set)
1353  CALL calculate_wavefunction(mo%mo_coeff, im, wf_r, wf_g, atomic_kind_set, &
1354  qs_kind_set, cell, dft_control, particle_set, pw_env)
1355  CALL cp_pw_to_cube(wf_r, unit_nr, title, particles=particles, &
1356  stride=section_get_ivals(input, 'MO_CUBES%STRIDE'))
1357 
1358  ! Close file
1359  CALL cp_print_key_finished_output(unit_nr, logger, input, 'MO_CUBES')
1360 
1361  ! Clean memory
1362  CALL auxbas_pw_pool%give_back_pw(wf_r)
1363  CALL auxbas_pw_pool%give_back_pw(wf_g)
1364 
1365  END SUBROUTINE save_mo_cube
1366 
1367 ! **************************************************************************************************
1368 !> \brief save specified electronic states to cube files
1369 !> \param qs_env QuickStep environment containing all system data
1370 !> \param ec electronic coupling data structure
1371 !> \param n_spins number of spin states
1372 !> \author Z. Futera (02.2017)
1373 ! **************************************************************************************************
1374  SUBROUTINE save_el_states(qs_env, ec, n_spins)
1375 
1376  ! Routine arguments
1377  TYPE(qs_environment_type), POINTER :: qs_env
1378  TYPE(et_cpl), POINTER :: ec
1379  INTEGER, INTENT(IN) :: n_spins
1380 
1381  INTEGER :: i, j, k, l, n
1382  INTEGER, DIMENSION(:), POINTER :: list
1383  TYPE(cp_logger_type), POINTER :: logger
1384  TYPE(mo_set_type), POINTER :: mo
1385  TYPE(section_vals_type), POINTER :: block_sec, mo_sec, print_sec
1386 
1387 ! Routine name for debug purposes
1388 
1389  NULLIFY (logger)
1390  NULLIFY (block_sec)
1391  NULLIFY (print_sec)
1392  NULLIFY (mo_sec)
1393 
1394  ! Output logger
1395  logger => cp_get_default_logger()
1396  block_sec => section_vals_get_subs_vals(qs_env%input, &
1397  'PROPERTIES%ET_COUPLING%PROJECTION%BLOCK')
1398 
1399  ! Print states of all blocks
1400  DO i = 1, ec%n_blocks
1401 
1402  print_sec => section_vals_get_subs_vals(block_sec, 'PRINT', i_rep_section=i)
1403 
1404  ! Check if the print input section is active
1405  IF (btest(cp_print_key_should_output(logger%iter_info, &
1406  print_sec, 'MO_CUBES'), cp_p_file)) THEN
1407 
1408  mo_sec => section_vals_get_subs_vals(print_sec, 'MO_CUBES')
1409 
1410  ! Spin states
1411  DO j = 1, n_spins
1412 
1413  mo => ec%block(i)%mo(j)
1414 
1415  CALL section_vals_val_get(mo_sec, keyword_name='MO_LIST', n_rep_val=n)
1416 
1417  ! List of specific MOs
1418  IF (n > 0) THEN
1419 
1420  DO k = 1, n
1421  NULLIFY (list)
1422  CALL section_vals_val_get(mo_sec, keyword_name='MO_LIST', &
1423  i_rep_val=k, i_vals=list)
1424  IF (ASSOCIATED(list)) THEN
1425  DO l = 1, SIZE(list)
1426  CALL save_mo_cube(qs_env, logger, print_sec, mo, i, list(l), j)
1427  END DO
1428  END IF
1429  END DO
1430 
1431  ! Frontier MOs
1432  ELSE
1433 
1434  ! Occupied states
1435  CALL section_vals_val_get(mo_sec, keyword_name='NHOMO', i_val=n)
1436 
1437  IF (n > 0) THEN
1438  DO k = max(1, mo%homo - n + 1), mo%homo
1439  CALL save_mo_cube(qs_env, logger, print_sec, mo, i, k, j)
1440  END DO
1441  END IF
1442 
1443  ! Unoccupied states
1444  CALL section_vals_val_get(mo_sec, keyword_name='NLUMO', i_val=n)
1445 
1446  IF (n > 0) THEN
1447  DO k = mo%lfomo, min(mo%lfomo + n - 1, mo%nmo)
1448  CALL save_mo_cube(qs_env, logger, print_sec, mo, i, k, j)
1449  END DO
1450  END IF
1451 
1452  END IF
1453 
1454  END DO
1455 
1456  END IF
1457 
1458  END DO
1459 
1460  END SUBROUTINE save_el_states
1461 
1462 ! **************************************************************************************************
1463 !> \brief calculates the electron transfer coupling elements by projection-operator approach
1464 !> Kondov et al. J.Phys.Chem.C 2007, 111, 11970-11981
1465 !> \param qs_env QuickStep environment containing all system data
1466 !> \author Z. Futera (02.2017)
1467 ! **************************************************************************************************
1468  SUBROUTINE calc_et_coupling_proj(qs_env)
1469 
1470  ! Routine arguments
1471  TYPE(qs_environment_type), POINTER :: qs_env
1472 
1473  INTEGER :: i, j, k, n_ao, n_atoms, output_unit
1474  LOGICAL :: do_kp, master
1475  TYPE(cp_blacs_env_type), POINTER :: blacs_env
1476  TYPE(cp_fm_struct_type), POINTER :: fm_s
1477  TYPE(cp_fm_type) :: mat_w
1478  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: mat_h
1479  TYPE(cp_logger_type), POINTER :: logger
1480  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks, mo_der
1481  TYPE(dft_control_type), POINTER :: dft_cntrl
1482  TYPE(et_cpl), POINTER :: ec
1483  TYPE(kpoint_type), POINTER :: kpoints
1484  TYPE(mo_set_type), DIMENSION(:), POINTER :: mo
1485  TYPE(mp_para_env_type), POINTER :: para_env
1486  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1487  TYPE(scf_control_type), POINTER :: scf_control
1488  TYPE(section_vals_type), POINTER :: et_proj_sec
1489 
1490 ! Routine name for debug purposes
1491 
1492  ! Pointer initialization
1493  NULLIFY (logger)
1494 
1495  NULLIFY (blacs_env)
1496  NULLIFY (para_env)
1497  NULLIFY (dft_cntrl)
1498  NULLIFY (kpoints)
1499  NULLIFY (qs_kind_set)
1500  NULLIFY (et_proj_sec)
1501 
1502  NULLIFY (fm_s)
1503  NULLIFY (ks, mo_der)
1504 
1505  NULLIFY (ec)
1506 
1507  ! Reference
1508  CALL cite_reference(futera2017)
1509 
1510  ! Stream for output to LOG file
1511  logger => cp_get_default_logger()
1512 
1513  et_proj_sec => section_vals_get_subs_vals(qs_env%input, 'PROPERTIES%ET_COUPLING%PROJECTION')
1514 
1515  output_unit = cp_print_key_unit_nr(logger, et_proj_sec, &
1516  'PROGRAM_RUN_INFO', extension='.log')
1517 
1518  ! Parallel calculation - master thread
1519  master = .false.
1520  IF (output_unit > 0) &
1521  master = .true.
1522 
1523  ! Header
1524  IF (master) THEN
1525  WRITE (output_unit, '(/,T2,A)') &
1526  '!-----------------------------------------------------------------------------!'
1527  WRITE (output_unit, '(T17,A)') &
1528  'Electronic coupling - Projection-operator method'
1529  END IF
1530 
1531  ! Main data structure
1532  ALLOCATE (ec)
1533  cpassert(ASSOCIATED(ec))
1534  CALL set_block_data(qs_env, et_proj_sec, ec)
1535 
1536  ! Number of atoms and AO functions
1537  CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, natom=n_atoms)
1538  CALL get_qs_kind_set(qs_kind_set, nsgf=n_ao)
1539 
1540  ! Print out info about system partitioning
1541  IF (master) THEN
1542 
1543  WRITE (output_unit, '(/,T3,A,I10)') &
1544  'Number of atoms = ', n_atoms
1545  WRITE (output_unit, '(T3,A,I10)') &
1546  'Number of fragments = ', ec%n_blocks
1547  WRITE (output_unit, '(T3,A,I10)') &
1548  'Number of fragment atoms = ', ec%n_atoms
1549  WRITE (output_unit, '(T3,A,I10)') &
1550  'Number of unassigned atoms = ', n_atoms - ec%n_atoms
1551  WRITE (output_unit, '(T3,A,I10)') &
1552  'Number of AO basis functions = ', n_ao
1553 
1554  DO i = 1, ec%n_blocks
1555 
1556  WRITE (output_unit, '(/,T3,A,I0,A)') &
1557  'Block ', i, ':'
1558  WRITE (output_unit, '(T3,A,I10)') &
1559  'Number of block atoms = ', ec%block(i)%n_atoms
1560  WRITE (output_unit, '(T3,A,I10)') &
1561  'Number of block electrons = ', ec%block(i)%n_electrons
1562  WRITE (output_unit, '(T3,A,I10)') &
1563  'Number of block AO functions = ', ec%block(i)%n_ao
1564 
1565  IF (ec%block(i)%n_atoms < 10) THEN
1566 
1567  WRITE (output_unit, '(T3,A,10I6)') &
1568  'Block atom IDs = ', &
1569  (ec%block(i)%atom(j)%id, j=1, ec%block(i)%n_atoms)
1570 
1571  ELSE
1572 
1573  WRITE (output_unit, '(T3,A)') 'Block atom IDs ='
1574  DO j = 1, ec%block(i)%n_atoms/10
1575  WRITE (output_unit, '(T3,A,10I6)') ' ', &
1576  (ec%block(i)%atom((j - 1)*10 + k)%id, k=1, 10)
1577  END DO
1578  IF (mod(ec%block(i)%n_atoms, 10) /= 0) THEN
1579  WRITE (output_unit, '(T3,A,10I6)') ' ', &
1580  (ec%block(i)%atom(k + 10*(ec%block(i)%n_atoms/10))%id, &
1581  k=1, mod(ec%block(i)%n_atoms, 10))
1582  END IF
1583 
1584  END IF
1585 
1586  END DO
1587 
1588  END IF
1589 
1590  ! Full matrix data structure
1591  CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
1592  CALL cp_fm_struct_create(fmstruct=fm_s, para_env=para_env, context=blacs_env, &
1593  nrow_global=n_ao, ncol_global=n_ao)
1594  CALL cp_fm_create(matrix=mat_w, matrix_struct=fm_s, name='FULL WORK MATRIX')
1595 
1596  ! Spin polarization / K-point sampling
1597  CALL get_qs_env(qs_env, dft_control=dft_cntrl, do_kpoints=do_kp)
1598  CALL get_qs_env(qs_env, mos=mo, matrix_ks=ks, mo_derivs=mo_der, scf_control=scf_control)
1599  CALL make_mo_eig(mo, dft_cntrl%nspins, ks, scf_control, mo_der)
1600 
1601  IF (do_kp) THEN
1602  cpabort('ET_COUPLING not implemented with kpoints')
1603  ELSE
1604  ! no K-points
1605  IF (master) &
1606  WRITE (output_unit, '(T3,A)') 'No K-point sampling (Gamma point only)'
1607  END IF
1608 
1609  IF (dft_cntrl%nspins == 2) THEN
1610 
1611  IF (master) &
1612  WRITE (output_unit, '(/,T3,A)') 'Spin-polarized calculation'
1613 
1614  !<--- Open shell / No K-points ------------------------------------------------>!
1615 
1616  ! State eneries of the whole system
1617  IF (mo(1)%nao /= mo(2)%nao) &
1618  cpabort('different number of alpha/beta AO basis functions')
1619  IF (master) THEN
1620  WRITE (output_unit, '(/,T3,A,I10)') &
1621  'Number of AO basis functions = ', mo(1)%nao
1622  WRITE (output_unit, '(T3,A,I10)') &
1623  'Number of alpha states = ', mo(1)%nmo
1624  WRITE (output_unit, '(T3,A,I10)') &
1625  'Number of beta states = ', mo(2)%nmo
1626  END IF
1627  CALL print_states(output_unit, mo, dft_cntrl%nspins, 'the whole system', fermi=.true.)
1628  CALL set_fermi(ec, mo(1)%mu, mo(2)%mu)
1629 
1630  ! KS Hamiltonian
1631  CALL get_block_hamiltonian(qs_env, ec, fm_s, mat_h, mat_w, n_ao, dft_cntrl%nspins)
1632 
1633  ! Block diagonization
1634  CALL hamiltonian_block_diag(qs_env, ec, mat_h)
1635 
1636  ! Print out energies and couplings
1637  DO i = 1, ec%n_blocks
1638  IF (output_unit > 0) THEN
1639  CALL print_states(output_unit, ec%block(i)%mo, dft_cntrl%nspins, &
1640  'block '//trim(adjustl(cp_to_string(i)))//' states', &
1641  mx_mo_a=mo(1)%nmo, mx_mo_b=mo(2)%nmo, fermi=.true.)
1642  END IF
1643  CALL print_mo_coeff(output_unit, qs_env, ec, i, dft_cntrl%nspins)
1644  END DO
1645 
1646  CALL print_couplings(et_proj_sec, output_unit, logger, ec, mo)
1647 
1648  ELSE
1649 
1650  IF (master) &
1651  WRITE (output_unit, '(/,T3,A)') 'Spin-restricted calculation'
1652 
1653  !<--- Close shell / No K-points ----------------------------------------------->!
1654 
1655  ! State eneries of the whole system
1656  IF (master) THEN
1657  WRITE (output_unit, '(/,T3,A,I10)') &
1658  'Number of AO basis functions = ', mo(1)%nao
1659  WRITE (output_unit, '(T3,A,I10)') &
1660  'Number of states = ', mo(1)%nmo
1661  END IF
1662  CALL print_states(output_unit, mo, dft_cntrl%nspins, 'the whole system', fermi=.true.)
1663  CALL set_fermi(ec, mo(1)%mu)
1664 
1665  ! KS Hamiltonian
1666  CALL get_block_hamiltonian(qs_env, ec, fm_s, mat_h, mat_w, n_ao, dft_cntrl%nspins)
1667 
1668  ! Block diagonization
1669  CALL hamiltonian_block_diag(qs_env, ec, mat_h)
1670 
1671  ! Print out energies and couplings
1672  DO i = 1, ec%n_blocks
1673  IF (output_unit > 0) THEN
1674  CALL print_states(output_unit, ec%block(i)%mo, dft_cntrl%nspins, &
1675  'block '//trim(adjustl(cp_to_string(i)))//' states', &
1676  mx_mo_a=mo(1)%nmo, fermi=.true.)
1677  END IF
1678  CALL print_mo_coeff(output_unit, qs_env, ec, i, dft_cntrl%nspins)
1679  END DO
1680 
1681  CALL print_couplings(et_proj_sec, output_unit, logger, ec, mo)
1682 
1683  END IF
1684 
1685  ! Save electronic states
1686  CALL save_el_states(qs_env, ec, dft_cntrl%nspins)
1687 
1688  ! Footer
1689  IF (master) WRITE (output_unit, '(/,T2,A)') &
1690  '!-----------------------------------------------------------------------------!'
1691 
1692  ! Clean memory
1693  CALL cp_fm_struct_release(fmstruct=fm_s)
1694  CALL cp_fm_release(matrix=mat_w)
1695  IF (ALLOCATED(mat_h)) THEN
1696  DO i = 1, SIZE(mat_h)
1697  CALL cp_fm_release(matrix=mat_h(i))
1698  END DO
1699  DEALLOCATE (mat_h)
1700  END IF
1701  CALL release_ec_data(ec)
1702 
1703  ! Close output stream
1704  CALL cp_print_key_finished_output(output_unit, logger, et_proj_sec, 'PROGRAM_RUN_INFO')
1705 
1706  END SUBROUTINE calc_et_coupling_proj
1707 
1708 END MODULE et_coupling_proj
atom
Definition: atom.F:9
atomic_kind_types
Define the atomic kind types and their sub types.
Definition: atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition: atomic_kind_types.F:138
basis_set_types
Definition: basis_set_types.F:15
basis_set_types::get_gto_basis_set
subroutine, public get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, npgf_sum, nshell_sum, maxder, short_kind_radius)
...
Definition: basis_set_types.F:615
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition: bibliography.F:28
bibliography::futera2017
integer, save, public futera2017
Definition: bibliography.F:43
cell_types
Handles all functions related to the CELL.
Definition: cell_types.F:15
cp_blacs_env
methods related to the blacs parallel environment
Definition: cp_blacs_env.F:15
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition: cp_control_types.F:12
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition: cp_dbcsr_operations.F:18
cp_dbcsr_operations::cp_dbcsr_sm_fm_multiply
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
Definition: cp_dbcsr_operations.F:744
cp_dbcsr_operations::copy_dbcsr_to_fm
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
Definition: cp_dbcsr_operations.F:208
cp_fm_basic_linalg
basic linear algebra operations for full matrices
Definition: cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_column_scale
subroutine, public cp_fm_column_scale(matrixa, scaling)
scales column i of matrix a with scaling(i)
Definition: cp_fm_basic_linalg.F:1882
cp_fm_diag
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
Definition: cp_fm_diag.F:17
cp_fm_diag::cp_fm_power
subroutine, public cp_fm_power(matrix, work, exponent, threshold, n_dependent, verbose, eigvals)
...
Definition: cp_fm_diag.F:896
cp_fm_diag::choose_eigv_solver
subroutine, public choose_eigv_solver(matrix, eigenvectors, eigenvalues, info)
Choose the Eigensolver depending on which library is available ELPA seems to be unstable for small sy...
Definition: cp_fm_diag.F:208
cp_fm_struct
represent the structure of a full matrix
Definition: cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition: cp_fm_struct.F:125
cp_fm_struct::cp_fm_struct_equivalent
logical function, public cp_fm_struct_equivalent(fmstruct1, fmstruct2)
returns true if the two matrix structures are equivalent, false otherwise.
Definition: cp_fm_struct.F:357
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition: cp_fm_struct.F:320
cp_fm_types
represent a full matrix distributed on many processors
Definition: cp_fm_types.F:15
cp_fm_types::cp_fm_vectorssum
subroutine, public cp_fm_vectorssum(matrix, sum_array, dir)
summing up all the elements along the matrix's i-th index or
Definition: cp_fm_types.F:1252
cp_fm_types::cp_fm_get_element
subroutine, public cp_fm_get_element(matrix, irow_global, icol_global, alpha, local)
returns an element of a fm this value is valid on every cpu using this call is expensive
Definition: cp_fm_types.F:643
cp_fm_types::cp_fm_to_fm_submat
subroutine, public cp_fm_to_fm_submat(msource, mtarget, nrow, ncol, s_firstrow, s_firstcol, t_firstrow, t_firstcol)
copy just a part ot the matrix
Definition: cp_fm_types.F:1473
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition: cp_fm_types.F:535
cp_fm_types::cp_fm_get_submatrix
subroutine, public cp_fm_get_submatrix(fm, target_m, start_row, start_col, n_rows, n_cols, transpose)
gets a submatrix of a full matrix op(target_m)(1:n_rows,1:n_cols) =fm(start_row:start_row+n_rows,...
Definition: cp_fm_types.F:901
cp_fm_types::cp_fm_set_element
subroutine, public cp_fm_set_element(matrix, irow_global, icol_global, alpha)
sets an element of a matrix
Definition: cp_fm_types.F:700
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition: cp_fm_types.F:167
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition: cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition: cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition: cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition: cp_output_handling.F:803
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition: cp_output_handling.F:1013
cp_output_handling::cp_p_file
integer, parameter, public cp_p_file
Definition: cp_output_handling.F:103
cp_output_handling::cp_print_key_should_output
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Definition: cp_output_handling.F:345
cp_realspace_grid_cube
A wrapper around pw_to_cube() which accepts particle_list_type.
Definition: cp_realspace_grid_cube.F:12
cp_realspace_grid_cube::cp_pw_to_cube
subroutine, public cp_pw_to_cube(pw, unit_nr, title, particles, stride, zero_tails, silent, mpi_io)
...
Definition: cp_realspace_grid_cube.F:45
et_coupling_proj
calculates the electron transfer coupling elements by projection-operator approach Kondov et al....
Definition: et_coupling_proj.F:13
et_coupling_proj::calc_et_coupling_proj
subroutine, public calc_et_coupling_proj(qs_env)
calculates the electron transfer coupling elements by projection-operator approach Kondov et al....
Definition: et_coupling_proj.F:1469
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition: input_section_types.F:15
input_section_types::section_get_ivals
integer function, dimension(:), pointer, public section_get_ivals(section_vals, keyword_name)
...
Definition: input_section_types.F:989
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition: input_section_types.F:724
input_section_types::section_vals_get
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
Definition: input_section_types.F:697
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition: input_section_types.F:1040
input_section_types::section_get_lval
logical function, public section_get_lval(section_vals, keyword_name)
...
Definition: input_section_types.F:1005
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition: kinds.F:57
kinds::default_path_length
integer, parameter, public default_path_length
Definition: kinds.F:58
kpoint_types
Types and basic routines needed for a kpoint calculation.
Definition: kpoint_types.F:15
list
An array-based list which grows on demand. When the internal array is full, a new array of twice the ...
Definition: list.F:24
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition: orbital_pointers.F:15
orbital_pointers::nso
integer, dimension(:), allocatable, public nso
Definition: orbital_pointers.F:42
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition: parallel_gemm_api.F:14
particle_list_types
represent a simple array based list of the given type
Definition: particle_list_types.F:15
particle_types
Define the data structure for the particle information.
Definition: particle_types.F:19
physcon
Definition of physical constants:
Definition: physcon.F:68
physcon::evolt
real(kind=dp), parameter, public evolt
Definition: physcon.F:183
pw_env_types
container for various plainwaves related things
Definition: pw_env_types.F:14
pw_env_types::pw_env_get
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Definition: pw_env_types.F:113
pw_pool_types
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition: pw_pool_types.F:24
pw_types
Definition: pw_types.F:24
qs_collocate_density
Calculate the plane wave density by collocating the primitive Gaussian functions (pgf).
Definition: qs_collocate_density.F:38
qs_collocate_density::calculate_wavefunction
subroutine, public calculate_wavefunction(mo_vectors, ivector, rho, rho_gspace, atomic_kind_set, qs_kind_set, cell, dft_control, particle_set, pw_env, basis_type, external_vector)
maps a given wavefunction on the grid
Definition: qs_collocate_density.F:2860
qs_environment_types
Definition: qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Definition: qs_environment_types.F:502
qs_kind_types
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_r3d_rs_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, U_of_dft_plus_u, J_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, J0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition: qs_kind_types.F:451
qs_kind_types::get_qs_kind_set
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, basis_type, total_zeff_corr)
Get attributes of an atomic kind set.
Definition: qs_kind_types.F:864
qs_mo_methods
collects routines that perform operations directly related to MOs
Definition: qs_mo_methods.F:14
qs_mo_methods::make_mo_eig
subroutine, public make_mo_eig(mos, nspins, ks_rmpv, scf_control, mo_derivs, admm_env)
Calculate KS eigenvalues starting from OF MOS.
Definition: qs_mo_methods.F:817
qs_mo_occupation
Set occupation of molecular orbitals.
Definition: qs_mo_occupation.F:15
qs_mo_types
Definition and initialisation of the mo data type.
Definition: qs_mo_types.F:22
qs_mo_types::allocate_mo_set
subroutine, public allocate_mo_set(mo_set, nao, nmo, nelectron, n_el_f, maxocc, flexible_electron_count)
Allocates a mo set and partially initializes it (nao,nmo,nelectron, and flexible_electron_count are v...
Definition: qs_mo_types.F:206
qs_mo_types::deallocate_mo_set
subroutine, public deallocate_mo_set(mo_set)
Deallocate a wavefunction data structure.
Definition: qs_mo_types.F:352
qs_subsys_types
types that represent a quickstep subsys
Definition: qs_subsys_types.F:12
qs_subsys_types::qs_subsys_get
subroutine, public qs_subsys_get(subsys, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell, cell_ref, use_ref_cell, energy, force, qs_kind_set, cp_subsys, nelectron_total, nelectron_spin)
...
Definition: qs_subsys_types.F:134
scf_control_types
parameters that control an scf iteration
Definition: scf_control_types.F:17