(git:b195825)
negf_env_types.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 Environment for NEGF based quantum transport calculations
10 ! **************************************************************************************************
12  USE cell_types, ONLY: cell_type,&
14  USE cp_blacs_env, ONLY: cp_blacs_env_type
15  USE cp_control_types, ONLY: dft_control_type
18  cp_fm_struct_type
19  USE cp_fm_types, ONLY: cp_fm_create,&
20  cp_fm_release,&
21  cp_fm_type
22  USE dbcsr_api, ONLY: dbcsr_copy,&
23  dbcsr_deallocate_matrix,&
24  dbcsr_init_p,&
25  dbcsr_p_type,&
26  dbcsr_set,&
27  dbcsr_type
28  USE force_env_types, ONLY: force_env_get,&
29  force_env_p_type,&
30  force_env_type,&
32  USE input_section_types, ONLY: section_vals_type,&
34  USE kinds, ONLY: default_string_length,&
35  dp
36  USE kpoint_types, ONLY: get_kpoint_env,&
38  kpoint_env_p_type,&
39  kpoint_type
40  USE message_passing, ONLY: mp_para_env_type
41  USE negf_atom_map, ONLY: negf_atom_map_type,&
43  USE negf_control_types, ONLY: negf_control_contact_type,&
44  negf_control_type
50  USE negf_subgroup_types, ONLY: negf_subgroup_env_type
53  USE particle_types, ONLY: particle_type
54  USE pw_env_types, ONLY: pw_env_get,&
55  pw_env_type
56  USE pw_pool_types, ONLY: pw_pool_type
57  USE pw_types, ONLY: pw_r3d_rs_type
60  mixing_storage_type
61  USE qs_energy, ONLY: qs_energies
63  USE qs_environment_types, ONLY: get_qs_env,&
64  qs_environment_type
65  USE qs_integrate_potential, ONLY: integrate_v_rspace
66  USE qs_mo_types, ONLY: get_mo_set,&
67  mo_set_type
68  USE qs_rho_types, ONLY: qs_rho_get,&
69  qs_rho_type
70  USE qs_subsys_types, ONLY: qs_subsys_get,&
71  qs_subsys_type
72 #include "./base/base_uses.f90"
73 
74  IMPLICIT NONE
75  PRIVATE
76 
77  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'negf_env_types'
78  LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .true.
79 
80  PUBLIC :: negf_env_type, negf_env_contact_type
82 
83 ! **************************************************************************************************
84 !> \brief Contact-specific NEGF environment.
85 !> \author Sergey Chulkov
86 ! **************************************************************************************************
87  TYPE negf_env_contact_type
88  REAL(kind=dp), DIMENSION(3) :: direction_vector, origin
89  REAL(kind=dp), DIMENSION(3) :: direction_vector_bias, origin_bias
90  !> an axis towards the secondary contact unit cell which coincides with the transport direction
91  !> 0 (undefined), 1 (+x), 2 (+y), 3 (+z), -1 (-x), -2 (-y), -3 (-z)
92  INTEGER :: direction_axis
93  !> atoms belonging to a primary contact unit cell
94  INTEGER, ALLOCATABLE, DIMENSION(:) :: atomlist_cell0
95  !> atoms belonging to a secondary contact unit cell (will be removed one day ...)
96  INTEGER, ALLOCATABLE, DIMENSION(:) :: atomlist_cell1
97  !> list of equivalent atoms in an appropriate contact force environment
98  TYPE(negf_atom_map_type), ALLOCATABLE, &
99  DIMENSION(:) :: atom_map_cell0, atom_map_cell1
100  !> energy of the HOMO
101  REAL(kind=dp) :: homo_energy
102  !> diagonal (h_00) and off-diagonal (h_01) blocks of the contact Kohn-Sham matrix ([number_of_spins]).
103  !> The matrix h_01 is of the shape [nao_cell0 x nao_cell1]
104  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: h_00, h_01
105  !> diagonal and off-diagonal blocks of the density matrix
106  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: rho_00, rho_01
107  !> diagonal and off-diagonal blocks of the overlap matrix
108  TYPE(cp_fm_type), POINTER :: s_00 => null(), s_01 => null()
109  END TYPE negf_env_contact_type
110 
111 ! **************************************************************************************************
112 !> \brief NEGF environment.
113 !> \author Sergey Chulkov
114 ! **************************************************************************************************
115  TYPE negf_env_type
116  !> contact-specific NEGF environments
117  TYPE(negf_env_contact_type), ALLOCATABLE, &
118  DIMENSION(:) :: contacts
119  !> Kohn-Sham matrix of the scattering region
120  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: h_s
121  !> Kohn-Sham matrix of the scattering region -- contact interface ([nspins, ncontacts])
122  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: h_sc
123  !> overlap matrix of the scattering region
124  TYPE(cp_fm_type), POINTER :: s_s => null()
125  !> an external Hartree potential in atomic basis set representation
126  TYPE(cp_fm_type), POINTER :: v_hartree_s => null()
127  !> overlap matrix of the scattering region -- contact interface for every contact ([ncontacts])
128  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: s_sc
129  !> structure needed for density mixing
130  TYPE(mixing_storage_type), POINTER :: mixing_storage
131  !> density mixing method
132  INTEGER :: mixing_method
133  END TYPE negf_env_type
134 
135 ! **************************************************************************************************
136 !> \brief Allocatable list of the type 'negf_atom_map_type'.
137 !> \author Sergey Chulkov
138 ! **************************************************************************************************
139  TYPE negf_atom_map_contact_type
140  TYPE(negf_atom_map_type), ALLOCATABLE, DIMENSION(:) :: atom_map
141  END TYPE negf_atom_map_contact_type
142 
143 CONTAINS
144 
145 ! **************************************************************************************************
146 !> \brief Create a new NEGF environment and compute the relevant Kohn-Sham matrices.
147 !> \param negf_env NEGF environment to create
148 !> \param sub_env NEGF parallel (sub)group environment
149 !> \param negf_control NEGF control
150 !> \param force_env the primary force environment
151 !> \param negf_mixing_section pointer to a mixing section within the NEGF input section
152 !> \param log_unit output unit number
153 !> \par History
154 !> * 01.2017 created [Sergey Chulkov]
155 ! **************************************************************************************************
156  SUBROUTINE negf_env_create(negf_env, sub_env, negf_control, force_env, negf_mixing_section, log_unit)
157  TYPE(negf_env_type), INTENT(inout) :: negf_env
158  TYPE(negf_subgroup_env_type), INTENT(in) :: sub_env
159  TYPE(negf_control_type), POINTER :: negf_control
160  TYPE(force_env_type), POINTER :: force_env
161  TYPE(section_vals_type), POINTER :: negf_mixing_section
162  INTEGER, INTENT(in) :: log_unit
163 
164  CHARACTER(len=*), PARAMETER :: routinen = 'negf_env_create'
165 
166  CHARACTER(len=default_string_length) :: contact_str, force_env_str, &
167  n_force_env_str
168  INTEGER :: handle, icontact, in_use, n_force_env, &
169  ncontacts
170  LOGICAL :: do_kpoints
171  TYPE(cp_blacs_env_type), POINTER :: blacs_env
172  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp, matrix_s_kp
173  TYPE(dft_control_type), POINTER :: dft_control
174  TYPE(force_env_p_type), DIMENSION(:), POINTER :: sub_force_env
175  TYPE(negf_atom_map_contact_type), ALLOCATABLE, &
176  DIMENSION(:) :: map_contact
177  TYPE(pw_r3d_rs_type), POINTER :: v_hartree_rspace
178  TYPE(qs_environment_type), POINTER :: qs_env, qs_env_contact
179  TYPE(qs_subsys_type), POINTER :: subsys, subsys_contact
180 
181  CALL timeset(routinen, handle)
182 
183  ! ensure we have Quickstep enabled for all force_env
184  NULLIFY (sub_force_env)
185  CALL force_env_get(force_env, in_use=in_use, qs_env=qs_env, sub_force_env=sub_force_env)
186 
187  IF (ASSOCIATED(sub_force_env)) THEN
188  n_force_env = SIZE(sub_force_env)
189  ELSE
190  n_force_env = 0
191  END IF
192 
193  IF (in_use == use_qs_force) THEN
194  DO icontact = 1, n_force_env
195  CALL force_env_get(sub_force_env(icontact)%force_env, in_use=in_use)
196  IF (in_use /= use_qs_force) EXIT
197  END DO
198  END IF
199 
200  IF (in_use /= use_qs_force) THEN
201  cpabort("Quickstep is required for NEGF run.")
202  END IF
203 
204  ! check that all mentioned FORCE_EVAL sections are actually present
205  ncontacts = SIZE(negf_control%contacts)
206 
207  DO icontact = 1, ncontacts
208  IF (negf_control%contacts(icontact)%force_env_index > n_force_env) THEN
209  WRITE (contact_str, '(I11)') icontact
210  WRITE (force_env_str, '(I11)') negf_control%contacts(icontact)%force_env_index
211  WRITE (n_force_env_str, '(I11)') n_force_env
212 
213  CALL cp_abort(__location__, &
214  "Contact number "//trim(adjustl(contact_str))//" is linked with the FORCE_EVAL section number "// &
215  trim(adjustl(force_env_str))//", however only "//trim(adjustl(n_force_env_str))// &
216  " FORCE_EVAL sections have been found. Note that FORCE_EVAL sections are enumerated from 0"// &
217  " and that the primary (0-th) section must contain all the atoms.")
218  END IF
219  END DO
220 
221  ! create basic matrices and neighbour lists for the primary force_env,
222  ! so we know how matrix elements are actually distributed across CPUs.
223  CALL qs_energies_init(qs_env, calc_forces=.false.)
224  CALL get_qs_env(qs_env, blacs_env=blacs_env, do_kpoints=do_kpoints, &
225  matrix_s_kp=matrix_s_kp, matrix_ks_kp=matrix_ks_kp, &
226  subsys=subsys, v_hartree_rspace=v_hartree_rspace)
227 
228  IF (do_kpoints) THEN
229  cpabort("k-points are currently not supported for device FORCE_EVAL")
230  END IF
231 
232  ! stage 1: map the atoms between the device force_env and all contact force_env-s
233  ALLOCATE (negf_env%contacts(ncontacts))
234  ALLOCATE (map_contact(ncontacts))
235 
236  DO icontact = 1, ncontacts
237  IF (negf_control%contacts(icontact)%force_env_index > 0) THEN
238  CALL force_env_get(sub_force_env(negf_control%contacts(icontact)%force_env_index)%force_env, qs_env=qs_env_contact)
239  CALL get_qs_env(qs_env_contact, subsys=subsys_contact)
240 
241  CALL negf_env_contact_init_maps(contact_env=negf_env%contacts(icontact), &
242  contact_control=negf_control%contacts(icontact), &
243  atom_map=map_contact(icontact)%atom_map, &
244  eps_geometry=negf_control%eps_geometry, &
245  subsys_device=subsys, &
246  subsys_contact=subsys_contact)
247 
248  IF (negf_env%contacts(icontact)%direction_axis == 0) THEN
249  WRITE (contact_str, '(I11)') icontact
250  WRITE (force_env_str, '(I11)') negf_control%contacts(icontact)%force_env_index
251  CALL cp_abort(__location__, &
252  "One lattice vector of the contact unit cell (FORCE_EVAL section "// &
253  trim(adjustl(force_env_str))//") must be parallel to the direction of the contact "// &
254  trim(adjustl(contact_str))//".")
255  END IF
256  END IF
257  END DO
258 
259  ! stage 2: obtain relevant Kohn-Sham matrix blocks for each contact (separate bulk DFT calculation)
260  DO icontact = 1, ncontacts
261  IF (negf_control%contacts(icontact)%force_env_index > 0) THEN
262  IF (log_unit > 0) &
263  WRITE (log_unit, '(/,T2,A,T70,I11)') "NEGF| Construct the Kohn-Sham matrix for the contact ", icontact
264 
265  CALL force_env_get(sub_force_env(negf_control%contacts(icontact)%force_env_index)%force_env, qs_env=qs_env_contact)
266 
267  CALL qs_energies(qs_env_contact, consistent_energies=.false., calc_forces=.false.)
268 
269  CALL negf_env_contact_init_matrices(contact_env=negf_env%contacts(icontact), sub_env=sub_env, &
270  qs_env_contact=qs_env_contact, matrix_s_device=matrix_s_kp(1, 1)%matrix)
271  END IF
272  END DO
273 
274  ! obtain an initial KS-matrix for the scattering region
275  CALL qs_energies(qs_env, consistent_energies=.false., calc_forces=.false.)
276 
277  ! *** obtain relevant Kohn-Sham matrix blocks for each contact with no separate FORCE_ENV ***
278  DO icontact = 1, ncontacts
279  IF (negf_control%contacts(icontact)%force_env_index <= 0) THEN
280  CALL negf_env_contact_init_matrices_gamma(contact_env=negf_env%contacts(icontact), &
281  contact_control=negf_control%contacts(icontact), &
282  sub_env=sub_env, qs_env=qs_env, &
283  eps_geometry=negf_control%eps_geometry)
284  END IF
285  END DO
286 
287  ! extract device-related matrix blocks
288  CALL negf_env_device_init_matrices(negf_env, negf_control, sub_env, qs_env)
289 
290  ! electron density mixing;
291  ! the input section below should be consistent with the subroutine create_negf_section()
292  NULLIFY (negf_env%mixing_storage)
293  CALL section_vals_val_get(negf_mixing_section, "METHOD", i_val=negf_env%mixing_method)
294 
295  CALL get_qs_env(qs_env, dft_control=dft_control)
296  ALLOCATE (negf_env%mixing_storage)
297  CALL mixing_storage_create(negf_env%mixing_storage, negf_mixing_section, &
298  negf_env%mixing_method, dft_control%qs_control%cutoff)
299 
300  CALL timestop(handle)
301  END SUBROUTINE negf_env_create
302 
303 ! **************************************************************************************************
304 !> \brief Establish mapping between the primary and the contact force environments
305 !> \param contact_env NEGF environment for the given contact (modified on exit)
306 !> \param contact_control NEGF control
307 !> \param atom_map atomic map
308 !> \param eps_geometry accuracy in mapping atoms between different force environments
309 !> \param subsys_device QuickStep subsystem of the device force environment
310 !> \param subsys_contact QuickStep subsystem of the contact force environment
311 !> \author Sergey Chulkov
312 ! **************************************************************************************************
313  SUBROUTINE negf_env_contact_init_maps(contact_env, contact_control, atom_map, &
314  eps_geometry, subsys_device, subsys_contact)
315  TYPE(negf_env_contact_type), INTENT(inout) :: contact_env
316  TYPE(negf_control_contact_type), INTENT(in) :: contact_control
317  TYPE(negf_atom_map_type), ALLOCATABLE, &
318  DIMENSION(:), INTENT(inout) :: atom_map
319  REAL(kind=dp), INTENT(in) :: eps_geometry
320  TYPE(qs_subsys_type), POINTER :: subsys_device, subsys_contact
321 
322  CHARACTER(LEN=*), PARAMETER :: routinen = 'negf_env_contact_init_maps'
323 
324  INTEGER :: handle, natoms
325 
326  CALL timeset(routinen, handle)
327 
328  CALL contact_direction_vector(contact_env%origin, &
329  contact_env%direction_vector, &
330  contact_env%origin_bias, &
331  contact_env%direction_vector_bias, &
332  contact_control%atomlist_screening, &
333  contact_control%atomlist_bulk, &
334  subsys_device)
335 
336  contact_env%direction_axis = contact_direction_axis(contact_env%direction_vector, subsys_contact, eps_geometry)
337 
338  IF (contact_env%direction_axis /= 0) THEN
339  natoms = SIZE(contact_control%atomlist_bulk)
340  ALLOCATE (atom_map(natoms))
341 
342  ! map atom listed in 'contact_control%atomlist_bulk' to the corresponding atom/cell replica from the contact force_env
343  CALL negf_map_atomic_indices(atom_map=atom_map, &
344  atom_list=contact_control%atomlist_bulk, &
345  subsys_device=subsys_device, &
346  subsys_contact=subsys_contact, &
347  eps_geometry=eps_geometry)
348 
349  ! list atoms from 'contact_control%atomlist_bulk' which belong to
350  ! the primary unit cell of the bulk region for the given contact
351  CALL list_atoms_in_bulk_primary_unit_cell(atomlist_cell0=contact_env%atomlist_cell0, &
352  atom_map_cell0=contact_env%atom_map_cell0, &
353  atomlist_bulk=contact_control%atomlist_bulk, &
354  atom_map=atom_map, &
355  origin=contact_env%origin, &
356  direction_vector=contact_env%direction_vector, &
357  direction_axis=contact_env%direction_axis, &
358  subsys_device=subsys_device)
359 
360  ! secondary unit cell
361  CALL list_atoms_in_bulk_secondary_unit_cell(atomlist_cell1=contact_env%atomlist_cell1, &
362  atom_map_cell1=contact_env%atom_map_cell1, &
363  atomlist_bulk=contact_control%atomlist_bulk, &
364  atom_map=atom_map, &
365  origin=contact_env%origin, &
366  direction_vector=contact_env%direction_vector, &
367  direction_axis=contact_env%direction_axis, &
368  subsys_device=subsys_device)
369  END IF
370 
371  CALL timestop(handle)
372  END SUBROUTINE negf_env_contact_init_maps
373 
374 ! **************************************************************************************************
375 !> \brief Extract relevant matrix blocks for the given contact.
376 !> \param contact_env NEGF environment for the contact (modified on exit)
377 !> \param sub_env NEGF parallel (sub)group environment
378 !> \param qs_env_contact QuickStep environment for the contact force environment
379 !> \param matrix_s_device overlap matrix from device force environment
380 !> \author Sergey Chulkov
381 ! **************************************************************************************************
382  SUBROUTINE negf_env_contact_init_matrices(contact_env, sub_env, qs_env_contact, matrix_s_device)
383  TYPE(negf_env_contact_type), INTENT(inout) :: contact_env
384  TYPE(negf_subgroup_env_type), INTENT(in) :: sub_env
385  TYPE(qs_environment_type), POINTER :: qs_env_contact
386  TYPE(dbcsr_type), POINTER :: matrix_s_device
387 
388  CHARACTER(LEN=*), PARAMETER :: routinen = 'negf_env_contact_init_matrices'
389 
390  INTEGER :: handle, iatom, ispin, nao, natoms, &
391  nimages, nspins
392  INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_list0, atom_list1
393  INTEGER, ALLOCATABLE, DIMENSION(:, :) :: index_to_cell, is_same_cell
394  INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
395  LOGICAL :: do_kpoints
396  TYPE(cp_fm_struct_type), POINTER :: fm_struct
397  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp, matrix_s_kp, rho_ao_kp
398  TYPE(dbcsr_type), POINTER :: matrix_s_ref
399  TYPE(dft_control_type), POINTER :: dft_control
400  TYPE(kpoint_type), POINTER :: kpoints
401  TYPE(mp_para_env_type), POINTER :: para_env
402  TYPE(qs_rho_type), POINTER :: rho_struct
403  TYPE(qs_subsys_type), POINTER :: subsys
404 
405  CALL timeset(routinen, handle)
406 
407  CALL get_qs_env(qs_env_contact, &
408  dft_control=dft_control, &
409  do_kpoints=do_kpoints, &
410  kpoints=kpoints, &
411  matrix_ks_kp=matrix_ks_kp, &
412  matrix_s_kp=matrix_s_kp, &
413  para_env=para_env, &
414  rho=rho_struct, &
415  subsys=subsys)
416  CALL qs_rho_get(rho_struct, rho_ao_kp=rho_ao_kp)
417 
418  CALL negf_homo_energy_estimate(contact_env%homo_energy, qs_env_contact)
419 
420  natoms = SIZE(contact_env%atomlist_cell0)
421  ALLOCATE (atom_list0(natoms))
422  DO iatom = 1, natoms
423  atom_list0(iatom) = contact_env%atom_map_cell0(iatom)%iatom
424 
425  ! with no k-points there is one-to-one correspondence between the primary unit cell
426  ! of the contact force_env and the first contact unit cell of the device force_env
427  IF (sum(abs(contact_env%atom_map_cell0(iatom)%cell(:))) > 0) THEN
428  cpabort("NEGF K-points are not currently supported")
429  END IF
430  END DO
431 
432  cpassert(SIZE(contact_env%atomlist_cell1) == natoms)
433  ALLOCATE (atom_list1(natoms))
434  DO iatom = 1, natoms
435  atom_list1(iatom) = contact_env%atom_map_cell1(iatom)%iatom
436  END DO
437 
438  nspins = dft_control%nspins
439  nimages = dft_control%nimages
440 
441  IF (do_kpoints) THEN
442  CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index)
443  ELSE
444  ALLOCATE (cell_to_index(0:0, 0:0, 0:0))
445  cell_to_index(0, 0, 0) = 1
446  END IF
447 
448  ALLOCATE (index_to_cell(3, nimages))
449  CALL invert_cell_to_index(cell_to_index, nimages, index_to_cell)
450  IF (.NOT. do_kpoints) DEALLOCATE (cell_to_index)
451 
452  NULLIFY (fm_struct)
453  nao = number_of_atomic_orbitals(subsys, atom_list0)
454  CALL cp_fm_struct_create(fm_struct, nrow_global=nao, ncol_global=nao, context=sub_env%blacs_env)
455 
456  ! ++ create matrices: s_00, s_01
457  ALLOCATE (contact_env%s_00, contact_env%s_01)
458  CALL cp_fm_create(contact_env%s_00, fm_struct)
459  CALL cp_fm_create(contact_env%s_01, fm_struct)
460 
461  ! ++ create matrices: h_00, h_01, rho_00, rho_01
462  ALLOCATE (contact_env%h_00(nspins), contact_env%h_01(nspins))
463  ALLOCATE (contact_env%rho_00(nspins), contact_env%rho_01(nspins))
464  DO ispin = 1, nspins
465  CALL cp_fm_create(contact_env%h_00(ispin), fm_struct)
466  CALL cp_fm_create(contact_env%h_01(ispin), fm_struct)
467  CALL cp_fm_create(contact_env%rho_00(ispin), fm_struct)
468  CALL cp_fm_create(contact_env%rho_01(ispin), fm_struct)
469  END DO
470 
471  CALL cp_fm_struct_release(fm_struct)
472 
473  NULLIFY (matrix_s_ref)
474  CALL dbcsr_init_p(matrix_s_ref)
475  CALL dbcsr_copy(matrix_s_ref, matrix_s_kp(1, 1)%matrix)
476  CALL dbcsr_set(matrix_s_ref, 0.0_dp)
477 
478  ALLOCATE (is_same_cell(natoms, natoms))
479 
480  CALL negf_reference_contact_matrix(matrix_contact=matrix_s_ref, &
481  matrix_device=matrix_s_device, &
482  atom_list=contact_env%atomlist_cell0, &
483  atom_map=contact_env%atom_map_cell0, &
484  para_env=para_env)
485 
486  ! extract matrices: s_00, s_01
487  CALL negf_copy_contact_matrix(fm_cell0=contact_env%s_00, &
488  fm_cell1=contact_env%s_01, &
489  direction_axis=contact_env%direction_axis, &
490  matrix_kp=matrix_s_kp(1, :), &
491  index_to_cell=index_to_cell, &
492  atom_list0=atom_list0, atom_list1=atom_list1, &
493  subsys=subsys, mpi_comm_global=para_env, &
494  is_same_cell=is_same_cell, matrix_ref=matrix_s_ref)
495 
496  ! extract matrices: h_00, h_01, rho_00, rho_01
497  DO ispin = 1, nspins
498  CALL negf_copy_contact_matrix(fm_cell0=contact_env%h_00(ispin), &
499  fm_cell1=contact_env%h_01(ispin), &
500  direction_axis=contact_env%direction_axis, &
501  matrix_kp=matrix_ks_kp(ispin, :), &
502  index_to_cell=index_to_cell, &
503  atom_list0=atom_list0, atom_list1=atom_list1, &
504  subsys=subsys, mpi_comm_global=para_env, &
505  is_same_cell=is_same_cell)
506 
507  CALL negf_copy_contact_matrix(fm_cell0=contact_env%rho_00(ispin), &
508  fm_cell1=contact_env%rho_01(ispin), &
509  direction_axis=contact_env%direction_axis, &
510  matrix_kp=rho_ao_kp(ispin, :), &
511  index_to_cell=index_to_cell, &
512  atom_list0=atom_list0, atom_list1=atom_list1, &
513  subsys=subsys, mpi_comm_global=para_env, &
514  is_same_cell=is_same_cell)
515  END DO
516 
517  DEALLOCATE (is_same_cell)
518  CALL dbcsr_deallocate_matrix(matrix_s_ref)
519  DEALLOCATE (index_to_cell)
520  DEALLOCATE (atom_list0, atom_list1)
521 
522  CALL timestop(handle)
523  END SUBROUTINE negf_env_contact_init_matrices
524 
525 ! **************************************************************************************************
526 !> \brief Extract relevant matrix blocks for the given contact using the device's force environment.
527 !> \param contact_env NEGF environment for the contact (modified on exit)
528 !> \param contact_control NEGF control for the contact
529 !> \param sub_env NEGF parallel (sub)group environment
530 !> \param qs_env QuickStep environment for the device force environment
531 !> \param eps_geometry accuracy in Cartesian coordinates
532 !> \author Sergey Chulkov
533 ! **************************************************************************************************
534  SUBROUTINE negf_env_contact_init_matrices_gamma(contact_env, contact_control, sub_env, qs_env, eps_geometry)
535  TYPE(negf_env_contact_type), INTENT(inout) :: contact_env
536  TYPE(negf_control_contact_type), INTENT(in) :: contact_control
537  TYPE(negf_subgroup_env_type), INTENT(in) :: sub_env
538  TYPE(qs_environment_type), POINTER :: qs_env
539  REAL(kind=dp), INTENT(in) :: eps_geometry
540 
541  CHARACTER(LEN=*), PARAMETER :: routinen = 'negf_env_contact_init_matrices_gamma'
542 
543  INTEGER :: handle, iatom, icell, ispin, nao_c, &
544  nspins
545  LOGICAL :: do_kpoints
546  REAL(kind=dp), DIMENSION(2) :: r2_origin_cell
547  REAL(kind=dp), DIMENSION(3) :: direction_vector, origin
548  TYPE(cp_fm_struct_type), POINTER :: fm_struct
549  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp, matrix_s_kp, rho_ao_kp
550  TYPE(dft_control_type), POINTER :: dft_control
551  TYPE(mp_para_env_type), POINTER :: para_env
552  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
553  TYPE(qs_rho_type), POINTER :: rho_struct
554  TYPE(qs_subsys_type), POINTER :: subsys
555 
556  CALL timeset(routinen, handle)
557 
558  CALL get_qs_env(qs_env, &
559  dft_control=dft_control, &
560  do_kpoints=do_kpoints, &
561  matrix_ks_kp=matrix_ks_kp, &
562  matrix_s_kp=matrix_s_kp, &
563  para_env=para_env, &
564  rho=rho_struct, &
565  subsys=subsys)
566  CALL qs_rho_get(rho_struct, rho_ao_kp=rho_ao_kp)
567 
568  IF (do_kpoints) THEN
569  CALL cp_abort(__location__, &
570  "K-points in device region have not been implemented yet.")
571  END IF
572 
573  nspins = dft_control%nspins
574 
575  nao_c = number_of_atomic_orbitals(subsys, contact_control%atomlist_cell(1)%vector)
576  IF (number_of_atomic_orbitals(subsys, contact_control%atomlist_cell(2)%vector) /= nao_c) THEN
577  CALL cp_abort(__location__, &
578  "Primary and secondary bulk contact cells should be identical "// &
579  "in terms of the number of atoms of each kind, and their basis sets. "// &
580  "No single atom, however, can be shared between these two cells.")
581  END IF
582 
583  contact_env%homo_energy = 0.0_dp
584 
585  CALL contact_direction_vector(contact_env%origin, &
586  contact_env%direction_vector, &
587  contact_env%origin_bias, &
588  contact_env%direction_vector_bias, &
589  contact_control%atomlist_screening, &
590  contact_control%atomlist_bulk, &
591  subsys)
592 
593  contact_env%direction_axis = contact_direction_axis(contact_env%direction_vector, subsys, eps_geometry)
594 
595  ! choose the primary and secondary contact unit cells
596  CALL qs_subsys_get(subsys, particle_set=particle_set)
597 
598  origin = particle_set(contact_control%atomlist_screening(1))%r
599  DO iatom = 2, SIZE(contact_control%atomlist_screening)
600  origin = origin + particle_set(contact_control%atomlist_screening(iatom))%r
601  END DO
602  origin = origin/real(SIZE(contact_control%atomlist_screening), kind=dp)
603 
604  DO icell = 1, 2
605  direction_vector = particle_set(contact_control%atomlist_cell(icell)%vector(1))%r
606  DO iatom = 2, SIZE(contact_control%atomlist_cell(icell)%vector)
607  direction_vector = direction_vector + particle_set(contact_control%atomlist_cell(icell)%vector(iatom))%r
608  END DO
609  direction_vector = direction_vector/real(SIZE(contact_control%atomlist_cell(icell)%vector), kind=dp)
610  direction_vector = direction_vector - origin
611  r2_origin_cell(icell) = dot_product(direction_vector, direction_vector)
612  END DO
613 
614  IF (sqrt(abs(r2_origin_cell(1) - r2_origin_cell(2))) < eps_geometry) THEN
615  ! primary and secondary bulk unit cells should not overlap;
616  ! currently we check that they are different by at least one atom that is, indeed, not sufficient.
617  CALL cp_abort(__location__, &
618  "Primary and secondary bulk contact cells should not overlap ")
619  ELSE IF (r2_origin_cell(1) < r2_origin_cell(2)) THEN
620  ALLOCATE (contact_env%atomlist_cell0(SIZE(contact_control%atomlist_cell(1)%vector)))
621  contact_env%atomlist_cell0(:) = contact_control%atomlist_cell(1)%vector(:)
622  ALLOCATE (contact_env%atomlist_cell1(SIZE(contact_control%atomlist_cell(2)%vector)))
623  contact_env%atomlist_cell1(:) = contact_control%atomlist_cell(2)%vector(:)
624  ELSE
625  ALLOCATE (contact_env%atomlist_cell0(SIZE(contact_control%atomlist_cell(2)%vector)))
626  contact_env%atomlist_cell0(:) = contact_control%atomlist_cell(2)%vector(:)
627  ALLOCATE (contact_env%atomlist_cell1(SIZE(contact_control%atomlist_cell(1)%vector)))
628  contact_env%atomlist_cell1(:) = contact_control%atomlist_cell(1)%vector(:)
629  END IF
630 
631  NULLIFY (fm_struct)
632  CALL cp_fm_struct_create(fm_struct, nrow_global=nao_c, ncol_global=nao_c, context=sub_env%blacs_env)
633  ALLOCATE (contact_env%h_00(nspins), contact_env%h_01(nspins))
634  ALLOCATE (contact_env%rho_00(nspins), contact_env%rho_01(nspins))
635  DO ispin = 1, nspins
636  CALL cp_fm_create(contact_env%h_00(ispin), fm_struct)
637  CALL cp_fm_create(contact_env%h_01(ispin), fm_struct)
638  CALL cp_fm_create(contact_env%rho_00(ispin), fm_struct)
639  CALL cp_fm_create(contact_env%rho_01(ispin), fm_struct)
640  END DO
641  ALLOCATE (contact_env%s_00, contact_env%s_01)
642  CALL cp_fm_create(contact_env%s_00, fm_struct)
643  CALL cp_fm_create(contact_env%s_01, fm_struct)
644  CALL cp_fm_struct_release(fm_struct)
645 
646  DO ispin = 1, nspins
647  CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=matrix_ks_kp(ispin, 1)%matrix, &
648  fm=contact_env%h_00(ispin), &
649  atomlist_row=contact_env%atomlist_cell0, &
650  atomlist_col=contact_env%atomlist_cell0, &
651  subsys=subsys, mpi_comm_global=para_env, &
652  do_upper_diag=.true., do_lower=.true.)
653  CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=matrix_ks_kp(ispin, 1)%matrix, &
654  fm=contact_env%h_01(ispin), &
655  atomlist_row=contact_env%atomlist_cell0, &
656  atomlist_col=contact_env%atomlist_cell1, &
657  subsys=subsys, mpi_comm_global=para_env, &
658  do_upper_diag=.true., do_lower=.true.)
659 
660  CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=rho_ao_kp(ispin, 1)%matrix, &
661  fm=contact_env%rho_00(ispin), &
662  atomlist_row=contact_env%atomlist_cell0, &
663  atomlist_col=contact_env%atomlist_cell0, &
664  subsys=subsys, mpi_comm_global=para_env, &
665  do_upper_diag=.true., do_lower=.true.)
666  CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=rho_ao_kp(ispin, 1)%matrix, &
667  fm=contact_env%rho_01(ispin), &
668  atomlist_row=contact_env%atomlist_cell0, &
669  atomlist_col=contact_env%atomlist_cell1, &
670  subsys=subsys, mpi_comm_global=para_env, &
671  do_upper_diag=.true., do_lower=.true.)
672  END DO
673 
674  CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=matrix_s_kp(1, 1)%matrix, &
675  fm=contact_env%s_00, &
676  atomlist_row=contact_env%atomlist_cell0, &
677  atomlist_col=contact_env%atomlist_cell0, &
678  subsys=subsys, mpi_comm_global=para_env, &
679  do_upper_diag=.true., do_lower=.true.)
680  CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=matrix_s_kp(1, 1)%matrix, &
681  fm=contact_env%s_01, &
682  atomlist_row=contact_env%atomlist_cell0, &
683  atomlist_col=contact_env%atomlist_cell1, &
684  subsys=subsys, mpi_comm_global=para_env, &
685  do_upper_diag=.true., do_lower=.true.)
686 
687  CALL timestop(handle)
688  END SUBROUTINE negf_env_contact_init_matrices_gamma
689 
690 ! **************************************************************************************************
691 !> \brief Extract relevant matrix blocks for the scattering region as well as
692 !> all the scattering -- contact interface regions.
693 !> \param negf_env NEGF environment (modified on exit)
694 !> \param negf_control NEGF control
695 !> \param sub_env NEGF parallel (sub)group environment
696 !> \param qs_env Primary QuickStep environment
697 !> \author Sergey Chulkov
698 ! **************************************************************************************************
699  SUBROUTINE negf_env_device_init_matrices(negf_env, negf_control, sub_env, qs_env)
700  TYPE(negf_env_type), INTENT(inout) :: negf_env
701  TYPE(negf_control_type), POINTER :: negf_control
702  TYPE(negf_subgroup_env_type), INTENT(in) :: sub_env
703  TYPE(qs_environment_type), POINTER :: qs_env
704 
705  CHARACTER(LEN=*), PARAMETER :: routinen = 'negf_env_device_init_matrices'
706 
707  INTEGER :: handle, icontact, ispin, nao_c, nao_s, &
708  ncontacts, nspins
709  LOGICAL :: do_kpoints
710  TYPE(cp_fm_struct_type), POINTER :: fm_struct
711  TYPE(dbcsr_p_type) :: hmat
712  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp, matrix_s_kp
713  TYPE(dft_control_type), POINTER :: dft_control
714  TYPE(mp_para_env_type), POINTER :: para_env
715  TYPE(pw_env_type), POINTER :: pw_env
716  TYPE(pw_pool_type), POINTER :: pw_pool
717  TYPE(pw_r3d_rs_type) :: v_hartree
718  TYPE(qs_subsys_type), POINTER :: subsys
719 
720  CALL timeset(routinen, handle)
721 
722  IF (ALLOCATED(negf_control%atomlist_S_screening)) THEN
723  CALL get_qs_env(qs_env, &
724  dft_control=dft_control, &
725  do_kpoints=do_kpoints, &
726  matrix_ks_kp=matrix_ks_kp, &
727  matrix_s_kp=matrix_s_kp, &
728  para_env=para_env, &
729  pw_env=pw_env, &
730  subsys=subsys)
731  CALL pw_env_get(pw_env, auxbas_pw_pool=pw_pool)
732 
733  IF (do_kpoints) THEN
734  CALL cp_abort(__location__, &
735  "K-points in device region have not been implemented yet.")
736  END IF
737 
738  ncontacts = SIZE(negf_control%contacts)
739  nspins = dft_control%nspins
740 
741  NULLIFY (fm_struct)
742  nao_s = number_of_atomic_orbitals(subsys, negf_control%atomlist_S_screening)
743 
744  ! ++ create matrices: h_s, s_s
745  NULLIFY (negf_env%s_s, negf_env%v_hartree_s, fm_struct)
746  ALLOCATE (negf_env%h_s(nspins))
747 
748  CALL cp_fm_struct_create(fm_struct, nrow_global=nao_s, ncol_global=nao_s, context=sub_env%blacs_env)
749  ALLOCATE (negf_env%s_s)
750  CALL cp_fm_create(negf_env%s_s, fm_struct)
751  DO ispin = 1, nspins
752  CALL cp_fm_create(negf_env%h_s(ispin), fm_struct)
753  END DO
754  ALLOCATE (negf_env%v_hartree_s)
755  CALL cp_fm_create(negf_env%v_hartree_s, fm_struct)
756  CALL cp_fm_struct_release(fm_struct)
757 
758  ! ++ create matrices: h_sc, s_sc
759  ALLOCATE (negf_env%h_sc(nspins, ncontacts), negf_env%s_sc(ncontacts))
760  DO icontact = 1, ncontacts
761  nao_c = number_of_atomic_orbitals(subsys, negf_env%contacts(icontact)%atomlist_cell0)
762  CALL cp_fm_struct_create(fm_struct, nrow_global=nao_s, ncol_global=nao_c, context=sub_env%blacs_env)
763 
764  CALL cp_fm_create(negf_env%s_sc(icontact), fm_struct)
765 
766  DO ispin = 1, nspins
767  CALL cp_fm_create(negf_env%h_sc(ispin, icontact), fm_struct)
768  END DO
769 
770  CALL cp_fm_struct_release(fm_struct)
771  END DO
772 
773  ! extract matrices: h_s, s_s
774  DO ispin = 1, nspins
775  CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=matrix_ks_kp(ispin, 1)%matrix, &
776  fm=negf_env%h_s(ispin), &
777  atomlist_row=negf_control%atomlist_S_screening, &
778  atomlist_col=negf_control%atomlist_S_screening, &
779  subsys=subsys, mpi_comm_global=para_env, &
780  do_upper_diag=.true., do_lower=.true.)
781  END DO
782 
783  CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=matrix_s_kp(1, 1)%matrix, &
784  fm=negf_env%s_s, &
785  atomlist_row=negf_control%atomlist_S_screening, &
786  atomlist_col=negf_control%atomlist_S_screening, &
787  subsys=subsys, mpi_comm_global=para_env, &
788  do_upper_diag=.true., do_lower=.true.)
789 
790  ! v_hartree_s
791  NULLIFY (hmat%matrix)
792  CALL dbcsr_init_p(hmat%matrix)
793  CALL dbcsr_copy(matrix_b=hmat%matrix, matrix_a=matrix_s_kp(1, 1)%matrix)
794  CALL dbcsr_set(hmat%matrix, 0.0_dp)
795 
796  CALL pw_pool%create_pw(v_hartree)
797  CALL negf_env_init_v_hartree(v_hartree, negf_env%contacts, negf_control%contacts)
798 
799  CALL integrate_v_rspace(v_rspace=v_hartree, hmat=hmat, qs_env=qs_env, &
800  calculate_forces=.false., compute_tau=.false., gapw=.false.)
801 
802  CALL pw_pool%give_back_pw(v_hartree)
803 
804  CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=hmat%matrix, &
805  fm=negf_env%v_hartree_s, &
806  atomlist_row=negf_control%atomlist_S_screening, &
807  atomlist_col=negf_control%atomlist_S_screening, &
808  subsys=subsys, mpi_comm_global=para_env, &
809  do_upper_diag=.true., do_lower=.true.)
810 
811  CALL dbcsr_deallocate_matrix(hmat%matrix)
812 
813  ! extract matrices: h_sc, s_sc
814  DO icontact = 1, ncontacts
815  DO ispin = 1, nspins
816  CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=matrix_ks_kp(ispin, 1)%matrix, &
817  fm=negf_env%h_sc(ispin, icontact), &
818  atomlist_row=negf_control%atomlist_S_screening, &
819  atomlist_col=negf_env%contacts(icontact)%atomlist_cell0, &
820  subsys=subsys, mpi_comm_global=para_env, &
821  do_upper_diag=.true., do_lower=.true.)
822  END DO
823 
824  CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=matrix_s_kp(1, 1)%matrix, &
825  fm=negf_env%s_sc(icontact), &
826  atomlist_row=negf_control%atomlist_S_screening, &
827  atomlist_col=negf_env%contacts(icontact)%atomlist_cell0, &
828  subsys=subsys, mpi_comm_global=para_env, &
829  do_upper_diag=.true., do_lower=.true.)
830  END DO
831  END IF
832 
833  CALL timestop(handle)
834  END SUBROUTINE negf_env_device_init_matrices
835 
836 ! **************************************************************************************************
837 !> \brief Contribution to the Hartree potential related to the external bias voltage.
838 !> \param v_hartree Hartree potential (modified on exit)
839 !> \param contact_env NEGF environment for every contact
840 !> \param contact_control NEGF control for every contact
841 !> \author Sergey Chulkov
842 ! **************************************************************************************************
843  SUBROUTINE negf_env_init_v_hartree(v_hartree, contact_env, contact_control)
844  TYPE(pw_r3d_rs_type), INTENT(IN) :: v_hartree
845  TYPE(negf_env_contact_type), DIMENSION(:), &
846  INTENT(in) :: contact_env
847  TYPE(negf_control_contact_type), DIMENSION(:), &
848  INTENT(in) :: contact_control
849 
850  CHARACTER(len=*), PARAMETER :: routinen = 'negf_env_init_v_hartree'
851  REAL(kind=dp), PARAMETER :: threshold = 16.0_dp*epsilon(0.0_dp)
852 
853  INTEGER :: dx, dy, dz, handle, icontact, ix, iy, &
854  iz, lx, ly, lz, ncontacts, ux, uy, uz
855  REAL(kind=dp) :: dvol, pot, proj, v1, v2
856  REAL(kind=dp), DIMENSION(3) :: dirvector_bias, point_coord, &
857  point_indices, vector
858 
859  CALL timeset(routinen, handle)
860 
861  ncontacts = SIZE(contact_env)
862  cpassert(SIZE(contact_control) == ncontacts)
863  cpassert(ncontacts == 2)
864 
865  dirvector_bias = contact_env(2)%origin_bias - contact_env(1)%origin_bias
866  v1 = contact_control(1)%v_external
867  v2 = contact_control(2)%v_external
868 
869  lx = v_hartree%pw_grid%bounds_local(1, 1)
870  ux = v_hartree%pw_grid%bounds_local(2, 1)
871  ly = v_hartree%pw_grid%bounds_local(1, 2)
872  uy = v_hartree%pw_grid%bounds_local(2, 2)
873  lz = v_hartree%pw_grid%bounds_local(1, 3)
874  uz = v_hartree%pw_grid%bounds_local(2, 3)
875 
876  dx = v_hartree%pw_grid%npts(1)/2
877  dy = v_hartree%pw_grid%npts(2)/2
878  dz = v_hartree%pw_grid%npts(3)/2
879 
880  dvol = v_hartree%pw_grid%dvol
881 
882  DO iz = lz, uz
883  point_indices(3) = real(iz + dz, kind=dp)
884  DO iy = ly, uy
885  point_indices(2) = real(iy + dy, kind=dp)
886 
887  DO ix = lx, ux
888  point_indices(1) = real(ix + dx, kind=dp)
889  point_coord(:) = matmul(v_hartree%pw_grid%dh, point_indices)
890 
891  vector = point_coord - contact_env(1)%origin_bias
892  proj = projection_on_direction_vector(vector, dirvector_bias)
893  IF (proj + threshold >= 0.0_dp .AND. proj - threshold <= 1.0_dp) THEN
894  ! scattering region
895  ! proj == 0 we are at the first contact boundary
896  ! proj == 1 we are at the second contact boundary
897  IF (proj < 0.0_dp) THEN
898  proj = 0.0_dp
899  ELSE IF (proj > 1.0_dp) THEN
900  proj = 1.0_dp
901  END IF
902  pot = v1 + (v2 - v1)*proj
903  ELSE
904  pot = 0.0_dp
905  DO icontact = 1, ncontacts
906  vector = point_coord - contact_env(icontact)%origin_bias
907  proj = projection_on_direction_vector(vector, contact_env(icontact)%direction_vector_bias)
908 
909  IF (proj + threshold >= 0.0_dp .AND. proj - threshold <= 1.0_dp) THEN
910  pot = contact_control(icontact)%v_external
911  EXIT
912  END IF
913  END DO
914  END IF
915 
916  v_hartree%array(ix, iy, iz) = pot*dvol
917  END DO
918  END DO
919  END DO
920 
921  CALL timestop(handle)
922  END SUBROUTINE negf_env_init_v_hartree
923 
924 ! **************************************************************************************************
925 !> \brief Detect the axis towards secondary unit cell.
926 !> \param direction_vector direction vector
927 !> \param subsys_contact QuickStep subsystem of the contact force environment
928 !> \param eps_geometry accuracy in mapping atoms between different force environments
929 !> \return direction axis: 0 (undefined), 1 (x), 2(y), 3 (z)
930 !> \par History
931 !> * 08.2017 created [Sergey Chulkov]
932 ! **************************************************************************************************
933  FUNCTION contact_direction_axis(direction_vector, subsys_contact, eps_geometry) RESULT(direction_axis)
934  REAL(kind=dp), DIMENSION(3), INTENT(in) :: direction_vector
935  TYPE(qs_subsys_type), POINTER :: subsys_contact
936  REAL(kind=dp), INTENT(in) :: eps_geometry
937  INTEGER :: direction_axis
938 
939  INTEGER :: i, naxes
940  REAL(kind=dp), DIMENSION(3) :: scaled
941  TYPE(cell_type), POINTER :: cell
942 
943  CALL qs_subsys_get(subsys_contact, cell=cell)
944  CALL real_to_scaled(scaled, direction_vector, cell)
945 
946  naxes = 0
947  direction_axis = 0 ! initialize to make GCC<=6 happy
948 
949  DO i = 1, 3
950  IF (abs(scaled(i)) > eps_geometry) THEN
951  IF (scaled(i) > 0.0_dp) THEN
952  direction_axis = i
953  ELSE
954  direction_axis = -i
955  END IF
956  naxes = naxes + 1
957  END IF
958  END DO
959 
960  ! direction_vector is not parallel to one of the unit cell's axis
961  IF (naxes /= 1) direction_axis = 0
962  END FUNCTION contact_direction_axis
963 
964 ! **************************************************************************************************
965 !> \brief Estimate energy of the highest spin-alpha occupied molecular orbital.
966 !> \param homo_energy HOMO energy (initialised on exit)
967 !> \param qs_env QuickStep environment
968 !> \par History
969 !> * 01.2017 created [Sergey Chulkov]
970 ! **************************************************************************************************
971  SUBROUTINE negf_homo_energy_estimate(homo_energy, qs_env)
972  REAL(kind=dp), INTENT(out) :: homo_energy
973  TYPE(qs_environment_type), POINTER :: qs_env
974 
975  CHARACTER(LEN=*), PARAMETER :: routinen = 'negf_homo_energy_estimate'
976  INTEGER, PARAMETER :: gamma_point = 1
977 
978  INTEGER :: handle, homo, ikpgr, ikpoint, imo, &
979  ispin, kplocal, nmo, nspins
980  INTEGER, DIMENSION(2) :: kp_range
981  LOGICAL :: do_kpoints
982  REAL(kind=dp) :: my_homo_energy
983  REAL(kind=dp), DIMENSION(:), POINTER :: eigenvalues
984  TYPE(kpoint_env_p_type), DIMENSION(:), POINTER :: kp_env
985  TYPE(kpoint_type), POINTER :: kpoints
986  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
987  TYPE(mo_set_type), DIMENSION(:, :), POINTER :: mos_kp
988  TYPE(mp_para_env_type), POINTER :: para_env, para_env_kp
989 
990  CALL timeset(routinen, handle)
991  my_homo_energy = 0.0_dp
992 
993  CALL get_qs_env(qs_env, para_env=para_env, mos=mos, kpoints=kpoints, do_kpoints=do_kpoints)
994 
995  IF (do_kpoints) THEN
996  CALL get_kpoint_info(kpoints, kp_env=kp_env, kp_range=kp_range, para_env_kp=para_env_kp)
997 
998  ! looking for a processor that holds the gamma point
999  IF (para_env_kp%mepos == 0 .AND. kp_range(1) <= gamma_point .AND. kp_range(2) >= gamma_point) THEN
1000  kplocal = kp_range(2) - kp_range(1) + 1
1001 
1002  DO ikpgr = 1, kplocal
1003  CALL get_kpoint_env(kp_env(ikpgr)%kpoint_env, nkpoint=ikpoint, mos=mos_kp)
1004 
1005  IF (ikpoint == gamma_point) THEN
1006  ! mos_kp(component, spin), where component = 1 (real), or 2 (imaginary)
1007  CALL get_mo_set(mos_kp(1, 1), homo=homo, eigenvalues=eigenvalues) ! mu=fermi_level
1008 
1009  my_homo_energy = eigenvalues(homo)
1010  EXIT
1011  END IF
1012  END DO
1013  END IF
1014 
1015  CALL para_env%sum(my_homo_energy)
1016  ELSE
1017  ! Hamiltonian of the bulk contact region has been computed without k-points.
1018  ! Try to obtain the HOMO energy assuming there is no OT. We probably should abort here
1019  ! as we do need a second replica of the bulk contact unit cell along transport
1020  ! direction anyway which is not available without k-points.
1021 
1022  cpwarn("It is strongly advised to use k-points within all contact FORCE_EVAL-s")
1023 
1024  nspins = SIZE(mos)
1025 
1026  spin_loop: DO ispin = 1, nspins
1027  CALL get_mo_set(mos(ispin), homo=homo, nmo=nmo, eigenvalues=eigenvalues)
1028 
1029  DO imo = nmo, 1, -1
1030  IF (eigenvalues(imo) /= 0.0_dp) EXIT spin_loop
1031  END DO
1032  END DO spin_loop
1033 
1034  IF (imo == 0) THEN
1035  cpabort("Orbital transformation (OT) for contact FORCE_EVAL-s is not supported")
1036  END IF
1037 
1038  my_homo_energy = eigenvalues(homo)
1039  END IF
1040 
1041  homo_energy = my_homo_energy
1042  CALL timestop(handle)
1043  END SUBROUTINE negf_homo_energy_estimate
1044 
1045 ! **************************************************************************************************
1046 !> \brief List atoms from the contact's primary unit cell.
1047 !> \param atomlist_cell0 list of atoms belonging to the contact's primary unit cell
1048 !> (allocate and initialised on exit)
1049 !> \param atom_map_cell0 atomic map of atoms from 'atomlist_cell0' (allocate and initialised on exit)
1050 !> \param atomlist_bulk list of atoms belonging to the bulk contact region
1051 !> \param atom_map atomic map of atoms from 'atomlist_bulk'
1052 !> \param origin origin of the contact
1053 !> \param direction_vector direction vector of the contact
1054 !> \param direction_axis axis towards secondary unit cell
1055 !> \param subsys_device QuickStep subsystem of the device force environment
1056 !> \par History
1057 !> * 08.2017 created [Sergey Chulkov]
1058 ! **************************************************************************************************
1059  SUBROUTINE list_atoms_in_bulk_primary_unit_cell(atomlist_cell0, atom_map_cell0, atomlist_bulk, atom_map, &
1060  origin, direction_vector, direction_axis, subsys_device)
1061  INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(inout) :: atomlist_cell0
1062  TYPE(negf_atom_map_type), ALLOCATABLE, &
1063  DIMENSION(:), INTENT(inout) :: atom_map_cell0
1064  INTEGER, DIMENSION(:), INTENT(in) :: atomlist_bulk
1065  TYPE(negf_atom_map_type), DIMENSION(:), INTENT(in) :: atom_map
1066  REAL(kind=dp), DIMENSION(3), INTENT(in) :: origin, direction_vector
1067  INTEGER, INTENT(in) :: direction_axis
1068  TYPE(qs_subsys_type), POINTER :: subsys_device
1069 
1070  CHARACTER(LEN=*), PARAMETER :: routinen = 'list_atoms_in_bulk_primary_unit_cell'
1071 
1072  INTEGER :: atom_min, dir_axis_min, &
1073  direction_axis_abs, handle, iatom, &
1074  natoms_bulk, natoms_cell0
1075  REAL(kind=dp) :: proj, proj_min
1076  REAL(kind=dp), DIMENSION(3) :: vector
1077  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1078 
1079  CALL timeset(routinen, handle)
1080  CALL qs_subsys_get(subsys_device, particle_set=particle_set)
1081 
1082  natoms_bulk = SIZE(atomlist_bulk)
1083  cpassert(SIZE(atom_map, 1) == natoms_bulk)
1084  direction_axis_abs = abs(direction_axis)
1085 
1086  ! looking for the nearest atom from the scattering region
1087  proj_min = 1.0_dp
1088  atom_min = 1
1089  DO iatom = 1, natoms_bulk
1090  vector = particle_set(atomlist_bulk(iatom))%r - origin
1091  proj = projection_on_direction_vector(vector, direction_vector)
1092 
1093  IF (proj < proj_min) THEN
1094  proj_min = proj
1095  atom_min = iatom
1096  END IF
1097  END DO
1098 
1099  dir_axis_min = atom_map(atom_min)%cell(direction_axis_abs)
1100 
1101  natoms_cell0 = 0
1102  DO iatom = 1, natoms_bulk
1103  IF (atom_map(iatom)%cell(direction_axis_abs) == dir_axis_min) &
1104  natoms_cell0 = natoms_cell0 + 1
1105  END DO
1106 
1107  ALLOCATE (atomlist_cell0(natoms_cell0))
1108  ALLOCATE (atom_map_cell0(natoms_cell0))
1109 
1110  natoms_cell0 = 0
1111  DO iatom = 1, natoms_bulk
1112  IF (atom_map(iatom)%cell(direction_axis_abs) == dir_axis_min) THEN
1113  natoms_cell0 = natoms_cell0 + 1
1114  atomlist_cell0(natoms_cell0) = atomlist_bulk(iatom)
1115  atom_map_cell0(natoms_cell0) = atom_map(iatom)
1116  END IF
1117  END DO
1118 
1119  CALL timestop(handle)
1120  END SUBROUTINE list_atoms_in_bulk_primary_unit_cell
1121 
1122 ! **************************************************************************************************
1123 !> \brief List atoms from the contact's secondary unit cell.
1124 !> \param atomlist_cell1 list of atoms belonging to the contact's secondary unit cell
1125 !> (allocate and initialised on exit)
1126 !> \param atom_map_cell1 atomic map of atoms from 'atomlist_cell1'
1127 !> (allocate and initialised on exit)
1128 !> \param atomlist_bulk list of atoms belonging to the bulk contact region
1129 !> \param atom_map atomic map of atoms from 'atomlist_bulk'
1130 !> \param origin origin of the contact
1131 !> \param direction_vector direction vector of the contact
1132 !> \param direction_axis axis towards the secondary unit cell
1133 !> \param subsys_device QuickStep subsystem of the device force environment
1134 !> \par History
1135 !> * 11.2017 created [Sergey Chulkov]
1136 !> \note Cloned from list_atoms_in_bulk_primary_unit_cell. Will be removed once we can managed to
1137 !> maintain consistency between real-space matrices from different force_eval sections.
1138 ! **************************************************************************************************
1139  SUBROUTINE list_atoms_in_bulk_secondary_unit_cell(atomlist_cell1, atom_map_cell1, atomlist_bulk, atom_map, &
1140  origin, direction_vector, direction_axis, subsys_device)
1141  INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(inout) :: atomlist_cell1
1142  TYPE(negf_atom_map_type), ALLOCATABLE, &
1143  DIMENSION(:), INTENT(inout) :: atom_map_cell1
1144  INTEGER, DIMENSION(:), INTENT(in) :: atomlist_bulk
1145  TYPE(negf_atom_map_type), DIMENSION(:), INTENT(in) :: atom_map
1146  REAL(kind=dp), DIMENSION(3), INTENT(in) :: origin, direction_vector
1147  INTEGER, INTENT(in) :: direction_axis
1148  TYPE(qs_subsys_type), POINTER :: subsys_device
1149 
1150  CHARACTER(LEN=*), PARAMETER :: routinen = 'list_atoms_in_bulk_secondary_unit_cell'
1151 
1152  INTEGER :: atom_min, dir_axis_min, &
1153  direction_axis_abs, handle, iatom, &
1154  natoms_bulk, natoms_cell1, offset
1155  REAL(kind=dp) :: proj, proj_min
1156  REAL(kind=dp), DIMENSION(3) :: vector
1157  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1158 
1159  CALL timeset(routinen, handle)
1160  CALL qs_subsys_get(subsys_device, particle_set=particle_set)
1161 
1162  natoms_bulk = SIZE(atomlist_bulk)
1163  cpassert(SIZE(atom_map, 1) == natoms_bulk)
1164  direction_axis_abs = abs(direction_axis)
1165  offset = sign(1, direction_axis)
1166 
1167  ! looking for the nearest atom from the scattering region
1168  proj_min = 1.0_dp
1169  atom_min = 1
1170  DO iatom = 1, natoms_bulk
1171  vector = particle_set(atomlist_bulk(iatom))%r - origin
1172  proj = projection_on_direction_vector(vector, direction_vector)
1173 
1174  IF (proj < proj_min) THEN
1175  proj_min = proj
1176  atom_min = iatom
1177  END IF
1178  END DO
1179 
1180  dir_axis_min = atom_map(atom_min)%cell(direction_axis_abs)
1181 
1182  natoms_cell1 = 0
1183  DO iatom = 1, natoms_bulk
1184  IF (atom_map(iatom)%cell(direction_axis_abs) == dir_axis_min + offset) &
1185  natoms_cell1 = natoms_cell1 + 1
1186  END DO
1187 
1188  ALLOCATE (atomlist_cell1(natoms_cell1))
1189  ALLOCATE (atom_map_cell1(natoms_cell1))
1190 
1191  natoms_cell1 = 0
1192  DO iatom = 1, natoms_bulk
1193  IF (atom_map(iatom)%cell(direction_axis_abs) == dir_axis_min + offset) THEN
1194  natoms_cell1 = natoms_cell1 + 1
1195  atomlist_cell1(natoms_cell1) = atomlist_bulk(iatom)
1196  atom_map_cell1(natoms_cell1) = atom_map(iatom)
1197  atom_map_cell1(natoms_cell1)%cell(direction_axis_abs) = dir_axis_min
1198  END IF
1199  END DO
1200 
1201  CALL timestop(handle)
1202  END SUBROUTINE list_atoms_in_bulk_secondary_unit_cell
1203 
1204 ! **************************************************************************************************
1205 !> \brief Release a NEGF environment variable.
1206 !> \param negf_env NEGF environment to release
1207 !> \par History
1208 !> * 01.2017 created [Sergey Chulkov]
1209 ! **************************************************************************************************
1210  SUBROUTINE negf_env_release(negf_env)
1211  TYPE(negf_env_type), INTENT(inout) :: negf_env
1212 
1213  CHARACTER(len=*), PARAMETER :: routinen = 'negf_env_release'
1214 
1215  INTEGER :: handle, icontact
1216 
1217  CALL timeset(routinen, handle)
1218 
1219  IF (ALLOCATED(negf_env%contacts)) THEN
1220  DO icontact = SIZE(negf_env%contacts), 1, -1
1221  CALL negf_env_contact_release(negf_env%contacts(icontact))
1222  END DO
1223 
1224  DEALLOCATE (negf_env%contacts)
1225  END IF
1226 
1227  ! h_s
1228  CALL cp_fm_release(negf_env%h_s)
1229 
1230  ! h_sc
1231  CALL cp_fm_release(negf_env%h_sc)
1232 
1233  ! s_s
1234  IF (ASSOCIATED(negf_env%s_s)) THEN
1235  CALL cp_fm_release(negf_env%s_s)
1236  DEALLOCATE (negf_env%s_s)
1237  NULLIFY (negf_env%s_s)
1238  END IF
1239 
1240  ! s_sc
1241  CALL cp_fm_release(negf_env%s_sc)
1242 
1243  ! v_hartree_s
1244  IF (ASSOCIATED(negf_env%v_hartree_s)) THEN
1245  CALL cp_fm_release(negf_env%v_hartree_s)
1246  DEALLOCATE (negf_env%v_hartree_s)
1247  NULLIFY (negf_env%v_hartree_s)
1248  END IF
1249 
1250  ! density mixing
1251  IF (ASSOCIATED(negf_env%mixing_storage)) THEN
1252  CALL mixing_storage_release(negf_env%mixing_storage)
1253  DEALLOCATE (negf_env%mixing_storage)
1254  END IF
1255 
1256  CALL timestop(handle)
1257  END SUBROUTINE negf_env_release
1258 
1259 ! **************************************************************************************************
1260 !> \brief Release a NEGF contact environment variable.
1261 !> \param contact_env NEGF contact environment to release
1262 ! **************************************************************************************************
1263  SUBROUTINE negf_env_contact_release(contact_env)
1264  TYPE(negf_env_contact_type), INTENT(inout) :: contact_env
1265 
1266  CHARACTER(len=*), PARAMETER :: routinen = 'negf_env_contact_release'
1267 
1268  INTEGER :: handle
1269 
1270  CALL timeset(routinen, handle)
1271 
1272  ! h_00
1273  CALL cp_fm_release(contact_env%h_00)
1274 
1275  ! h_01
1276  CALL cp_fm_release(contact_env%h_01)
1277 
1278  ! rho_00
1279  CALL cp_fm_release(contact_env%rho_00)
1280 
1281  ! rho_01
1282  CALL cp_fm_release(contact_env%rho_01)
1283 
1284  ! s_00
1285  IF (ASSOCIATED(contact_env%s_00)) THEN
1286  CALL cp_fm_release(contact_env%s_00)
1287  DEALLOCATE (contact_env%s_00)
1288  NULLIFY (contact_env%s_00)
1289  END IF
1290 
1291  ! s_01
1292  IF (ASSOCIATED(contact_env%s_01)) THEN
1293  CALL cp_fm_release(contact_env%s_01)
1294  DEALLOCATE (contact_env%s_01)
1295  NULLIFY (contact_env%s_01)
1296  END IF
1297 
1298  IF (ALLOCATED(contact_env%atomlist_cell0)) DEALLOCATE (contact_env%atomlist_cell0)
1299  IF (ALLOCATED(contact_env%atomlist_cell1)) DEALLOCATE (contact_env%atomlist_cell1)
1300  IF (ALLOCATED(contact_env%atom_map_cell0)) DEALLOCATE (contact_env%atom_map_cell0)
1301 
1302  CALL timestop(handle)
1303  END SUBROUTINE negf_env_contact_release
1304 END MODULE negf_env_types
Handles all functions related to the CELL.
Definition: cell_types.F:15
subroutine, public real_to_scaled(s, r, cell)
Transform real to scaled cell coordinates. s=h_inv*r.
Definition: cell_types.F:486
methods related to the blacs parallel environment
Definition: cp_blacs_env.F:15
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
represent the structure of a full matrix
Definition: cp_fm_struct.F:14
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
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition: cp_fm_struct.F:320
represent a full matrix distributed on many processors
Definition: cp_fm_types.F:15
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
Interface for the force calculations.
recursive subroutine, public force_env_get(force_env, in_use, fist_env, qs_env, meta_env, fp_env, subsys, para_env, potential_energy, additional_potential, kinetic_energy, harmonic_shell, kinetic_shell, cell, sub_force_env, qmmm_env, qmmmx_env, eip_env, pwdft_env, globenv, input, force_env_section, method_name_id, root_section, mixed_env, nnp_env, embed_env)
returns various attributes about the force environment
integer, parameter, public use_qs_force
objects that represent the structure of input sections and the data contained in an input section
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
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
integer, parameter, public default_string_length
Definition: kinds.F:57
Types and basic routines needed for a kpoint calculation.
Definition: kpoint_types.F:15
subroutine, public get_kpoint_env(kpoint_env, nkpoint, wkp, xkp, is_local, mos)
Get information from a single kpoint environment.
Definition: kpoint_types.F:747
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym)
Retrieve information from a kpoint environment.
Definition: kpoint_types.F:333
Interface to the message passing library MPI.
Map atoms between various force environments.
Definition: negf_atom_map.F:13
subroutine, public negf_map_atomic_indices(atom_map, atom_list, subsys_device, subsys_contact, eps_geometry)
Map atoms in the cell 'subsys_device' listed in 'atom_list' with the corresponding atoms in the cell ...
Definition: negf_atom_map.F:76
Input control types for NEGF based quantum transport calculations.
Environment for NEGF based quantum transport calculations.
subroutine, public negf_env_release(negf_env)
Release a NEGF environment variable.
subroutine, public negf_env_create(negf_env, sub_env, negf_control, force_env, negf_mixing_section, log_unit)
Create a new NEGF environment and compute the relevant Kohn-Sham matrices.
Helper routines to manipulate with matrices.
subroutine, public negf_reference_contact_matrix(matrix_contact, matrix_device, atom_list, atom_map, para_env)
Extract part of the DBCSR matrix based on selected atoms and copy it into another DBCSR matrix.
integer function, public number_of_atomic_orbitals(subsys, atom_list)
Compute the number of atomic orbitals of the given set of atoms.
subroutine, public invert_cell_to_index(cell_to_index, nimages, index_to_cell)
Invert cell_to_index mapping between unit cells and DBCSR matrix images.
subroutine, public negf_copy_contact_matrix(fm_cell0, fm_cell1, direction_axis, matrix_kp, index_to_cell, atom_list0, atom_list1, subsys, mpi_comm_global, is_same_cell, matrix_ref)
Driver routine to extract diagonal and off-diagonal blocks from a symmetric DBCSR matrix.
subroutine, public negf_copy_sym_dbcsr_to_fm_submat(matrix, fm, atomlist_row, atomlist_col, subsys, mpi_comm_global, do_upper_diag, do_lower)
Extract part of the DBCSR matrix based on selected atoms and copy it into a dense matrix.
Environment for NEGF based quantum transport calculations.
Routines to deal with vectors in 3-D real space.
Definition: negf_vectors.F:12
pure real(kind=dp) function, public projection_on_direction_vector(vector, vector0)
project the 'vector' onto the direction 'vector0'. Both vectors should have the same origin.
Definition: negf_vectors.F:134
subroutine, public contact_direction_vector(origin, direction_vector, origin_bias, direction_vector_bias, atomlist_screening, atomlist_bulk, subsys)
compute direction vector of the given contact
Definition: negf_vectors.F:44
Define the data structure for the particle information.
container for various plainwaves related things
Definition: pw_env_types.F:14
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
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
module that contains the definitions of the scf types
subroutine, public mixing_storage_release(mixing_store)
releases a mixing_storage
subroutine, public mixing_storage_create(mixing_store, mixing_section, mixing_method, ecut)
creates a mixing_storage
Utility subroutine for qs energy calculation.
subroutine, public qs_energies_init(qs_env, calc_forces)
Refactoring of qs_energies_scf. Driver routine for the initial setup and calculations for a qs energy...
Perform a QUICKSTEP wavefunction optimization (single point)
Definition: qs_energy.F:14
subroutine, public qs_energies(qs_env, consistent_energies, calc_forces)
Driver routine for QUICKSTEP single point wavefunction optimization.
Definition: qs_energy.F:65
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.
Integrate single or product functions over a potential on a RS grid.
Definition and initialisation of the mo data type.
Definition: qs_mo_types.F:22
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kTS, mu, flexible_electron_count)
Get the components of a MO set data structure.
Definition: qs_mo_types.F:397
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition: qs_rho_types.F:18
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition: qs_rho_types.F:229
types that represent a quickstep subsys
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)
...