Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief routines for DFT+NEGF calculations (coupling with the quantum transport code OMEN)
10 : !> \par History
11 : !> 12.2012 created external_scf_method [Hossein Bani-Hashemian]
12 : !> 05.2013 created rotines to work with C-interoperable matrices [Hossein Bani-Hashemian]
13 : !> 07.2013 created transport_env routines [Hossein Bani-Hashemian]
14 : !> 11.2014 switch to CSR matrices [Hossein Bani-Hashemian]
15 : !> 12.2014 merged [Hossein Bani-Hashemian]
16 : !> 04.2016 added imaginary DM and current density cube [Hossein Bani-Hashemian]
17 : !> \author Mohammad Hossein Bani-Hashemian
18 : ! **************************************************************************************************
19 : MODULE transport
20 : USE ISO_C_BINDING, ONLY: C_ASSOCIATED,&
21 : C_BOOL,&
22 : C_DOUBLE,&
23 : C_F_PROCPOINTER,&
24 : C_INT,&
25 : C_LOC,&
26 : C_NULL_PTR,&
27 : C_PTR
28 : USE atomic_kind_types, ONLY: atomic_kind_type,&
29 : get_atomic_kind
30 : USE bibliography, ONLY: Bruck2014,&
31 : cite_reference
32 : USE cp_control_types, ONLY: dft_control_type
33 : USE cp_dbcsr_api, ONLY: &
34 : dbcsr_convert_csr_to_dbcsr, dbcsr_convert_dbcsr_to_csr, dbcsr_copy, dbcsr_create, &
35 : dbcsr_csr_create, dbcsr_csr_create_from_dbcsr, dbcsr_csr_dbcsr_blkrow_dist, &
36 : dbcsr_csr_print_sparsity, dbcsr_csr_type, dbcsr_deallocate_matrix, dbcsr_desymmetrize, &
37 : dbcsr_has_symmetry, dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry, dbcsr_type_real_8
38 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_to_csr_screening
39 : USE cp_log_handling, ONLY: cp_get_default_logger,&
40 : cp_logger_get_default_unit_nr,&
41 : cp_logger_type
42 : USE cp_output_handling, ONLY: cp_p_file,&
43 : cp_print_key_finished_output,&
44 : cp_print_key_should_output,&
45 : cp_print_key_unit_nr
46 : USE cp_realspace_grid_cube, ONLY: cp_pw_to_cube
47 : USE input_section_types, ONLY: section_get_ivals,&
48 : section_vals_get,&
49 : section_vals_get_subs_vals,&
50 : section_vals_type,&
51 : section_vals_val_get
52 : USE kinds, ONLY: dp
53 : USE particle_list_types, ONLY: particle_list_type
54 : USE particle_methods, ONLY: get_particle_set
55 : USE particle_types, ONLY: particle_type
56 : USE physcon, ONLY: boltzmann,&
57 : e_charge,&
58 : evolt,&
59 : h_bar
60 : USE pw_env_types, ONLY: pw_env_get,&
61 : pw_env_type
62 : USE pw_methods, ONLY: pw_zero
63 : USE pw_pool_types, ONLY: pw_pool_type
64 : USE pw_types, ONLY: pw_c1d_gs_type,&
65 : pw_r3d_rs_type
66 : USE qs_environment_types, ONLY: get_qs_env,&
67 : qs_environment_type,&
68 : set_qs_env
69 : USE qs_kind_types, ONLY: get_qs_kind,&
70 : qs_kind_type
71 : USE qs_ks_types, ONLY: qs_ks_env_type
72 : USE qs_linres_current, ONLY: calculate_jrho_resp
73 : USE qs_linres_types, ONLY: current_env_type
74 : USE qs_rho_types, ONLY: qs_rho_get,&
75 : qs_rho_type
76 : USE qs_subsys_types, ONLY: qs_subsys_get,&
77 : qs_subsys_type
78 : USE transport_env_types, ONLY: cp2k_csr_interop_type,&
79 : cp2k_transport_parameters,&
80 : csr_interop_matrix_get_info,&
81 : csr_interop_nullify,&
82 : transport_env_type
83 : #include "./base/base_uses.f90"
84 :
85 : IMPLICIT NONE
86 :
87 : PRIVATE
88 :
89 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'transport'
90 :
91 : PUBLIC :: transport_env_create, transport_initialize, external_scf_method
92 : PUBLIC :: qs_scf_post_transport
93 :
94 : !> interface between C/C++ and FORTRAN
95 : INTERFACE c_func_interface
96 : ! **************************************************************************************************
97 : !> \brief C routine that takes the S and H matrices as input and outputs a P matrix
98 : !> \param cp2k_transport_params transport parameters read form a CP2K input file
99 : !> \param s_mat C-interoperable overlap matrix
100 : !> \param ks_mat C-interoperable Kohn-Sham matrix
101 : !> \param p_mat C-interoperable density matrix
102 : !> \param imagp_mat C-interoperable imaginary density matrix
103 : !> \author Mohammad Hossein Bani-Hashemian
104 : ! **************************************************************************************************
105 : SUBROUTINE c_scf_routine(cp2k_transport_params, s_mat, ks_mat, p_mat, imagp_mat) BIND(C)
106 : IMPORT :: C_INT, C_PTR, cp2k_csr_interop_type, cp2k_transport_parameters
107 : IMPLICIT NONE
108 : TYPE(cp2k_transport_parameters), VALUE, INTENT(IN) :: cp2k_transport_params
109 : TYPE(cp2k_csr_interop_type), VALUE, INTENT(IN) :: s_mat
110 : TYPE(cp2k_csr_interop_type), VALUE, INTENT(IN) :: ks_mat
111 : TYPE(cp2k_csr_interop_type), INTENT(INOUT) :: p_mat
112 : TYPE(cp2k_csr_interop_type), INTENT(INOUT) :: imagp_mat
113 : END SUBROUTINE c_scf_routine
114 : END INTERFACE c_func_interface
115 :
116 : CONTAINS
117 :
118 : ! **************************************************************************************************
119 : !> \brief creates the transport environment
120 : !> \param[inout] qs_env the qs_env containing the transport_env
121 : !> \author Mohammad Hossein Bani-Hashemian
122 : ! **************************************************************************************************
123 0 : SUBROUTINE transport_env_create(qs_env)
124 : TYPE(qs_environment_type), POINTER :: qs_env
125 :
126 : CHARACTER(LEN=*), PARAMETER :: routineN = 'transport_env_create'
127 :
128 : INTEGER :: handle
129 : TYPE(dft_control_type), POINTER :: dft_control
130 : TYPE(section_vals_type), POINTER :: input
131 : TYPE(transport_env_type), POINTER :: transport_env
132 :
133 0 : CALL timeset(routineN, handle)
134 :
135 : CALL get_qs_env(qs_env, &
136 : transport_env=transport_env, &
137 : dft_control=dft_control, &
138 0 : input=input)
139 :
140 0 : CPASSERT(.NOT. ASSOCIATED(transport_env))
141 :
142 0 : ALLOCATE (transport_env)
143 :
144 0 : CALL transport_init_read_input(input, transport_env)
145 0 : CALL transport_set_contact_params(qs_env, transport_env)
146 :
147 0 : CALL set_qs_env(qs_env, transport_env=transport_env)
148 :
149 0 : CALL timestop(handle)
150 :
151 0 : END SUBROUTINE transport_env_create
152 :
153 : ! **************************************************************************************************
154 : !> \brief intitializes all fields of transport_env using the parameters read from
155 : !> the corresponding input section
156 : !> \param[inout] input the input file
157 : !> \param[inout] transport_env the transport_env to be initialized
158 : !> \author Mohammad Hossein Bani-Hashemian
159 : ! **************************************************************************************************
160 0 : SUBROUTINE transport_init_read_input(input, transport_env)
161 : TYPE(section_vals_type), POINTER :: input
162 : TYPE(transport_env_type), INTENT(INOUT) :: transport_env
163 :
164 : CHARACTER(len=*), PARAMETER :: routineN = 'transport_init_read_input'
165 :
166 : INTEGER :: contact_bandwidth, contact_injsign, &
167 : contact_natoms, contact_start, handle, &
168 : i, n_contacts, stride_contacts
169 0 : INTEGER, DIMENSION(:), POINTER :: i_vals
170 : LOGICAL :: contact_explicit, injecting_contact, &
171 : obc_equilibrium, one_circle
172 : TYPE(section_vals_type), POINTER :: beyn_section, contact_section, &
173 : pexsi_section, transport_section
174 :
175 0 : CALL timeset(routineN, handle)
176 :
177 0 : transport_section => section_vals_get_subs_vals(input, "DFT%TRANSPORT")
178 0 : contact_section => section_vals_get_subs_vals(transport_section, "CONTACT")
179 0 : beyn_section => section_vals_get_subs_vals(transport_section, "BEYN")
180 0 : pexsi_section => section_vals_get_subs_vals(transport_section, "PEXSI")
181 0 : CALL section_vals_get(contact_section, explicit=contact_explicit, n_repetition=n_contacts)
182 :
183 0 : NULLIFY (i_vals)
184 : ! read from input
185 0 : CALL section_vals_val_get(transport_section, "TRANSPORT_METHOD", i_val=transport_env%params%method)
186 0 : CALL section_vals_val_get(transport_section, "INJECTION_METHOD", i_val=transport_env%params%injection_method)
187 : CALL section_vals_val_get(transport_section, "REAL_AXIS_INTEGRATION_METHOD", &
188 0 : i_val=transport_env%params%rlaxis_integration_method)
189 0 : CALL section_vals_val_get(transport_section, "QT_FORMALISM", i_val=transport_env%params%qt_formalism)
190 0 : CALL section_vals_val_get(transport_section, "LINEAR_SOLVER", i_val=transport_env%params%linear_solver)
191 : CALL section_vals_val_get(transport_section, "MATRIX_INVERSION_METHOD", &
192 0 : i_val=transport_env%params%matrixinv_method)
193 0 : CALL section_vals_val_get(transport_section, "CONTACT_FILLING", i_val=transport_env%params%transport_neutral)
194 0 : CALL section_vals_val_get(transport_section, "N_KPOINTS", i_val=transport_env%params%n_kpoint)
195 0 : CALL section_vals_val_get(transport_section, "NUM_INTERVAL", i_val=transport_env%params%num_interval)
196 : CALL section_vals_val_get(transport_section, "TASKS_PER_ENERGY_POINT", &
197 0 : i_val=transport_env%params%tasks_per_energy_point)
198 0 : CALL section_vals_val_get(transport_section, "TASKS_PER_POLE", i_val=transport_env%params%tasks_per_pole)
199 0 : CALL section_vals_val_get(transport_section, "NUM_POLE", i_val=transport_env%params%num_pole)
200 0 : CALL section_vals_val_get(transport_section, "GPUS_PER_POINT", i_val=transport_env%params%gpus_per_point)
201 0 : CALL section_vals_val_get(transport_section, "N_POINTS_INV", i_val=transport_env%params%n_points_inv)
202 0 : CALL section_vals_val_get(transport_section, "COLZERO_THRESHOLD", r_val=transport_env%params%colzero_threshold)
203 0 : CALL section_vals_val_get(transport_section, "EPS_LIMIT", r_val=transport_env%params%eps_limit)
204 0 : CALL section_vals_val_get(transport_section, "EPS_LIMIT_CC", r_val=transport_env%params%eps_limit_cc)
205 0 : CALL section_vals_val_get(transport_section, "EPS_DECAY", r_val=transport_env%params%eps_decay)
206 : CALL section_vals_val_get(transport_section, "EPS_SINGULARITY_CURVATURES", &
207 0 : r_val=transport_env%params%eps_singularity_curvatures)
208 0 : CALL section_vals_val_get(transport_section, "EPS_MU", r_val=transport_env%params%eps_mu)
209 0 : CALL section_vals_val_get(transport_section, "EPS_EIGVAL_DEGEN", r_val=transport_env%params%eps_eigval_degen)
210 0 : CALL section_vals_val_get(transport_section, "EPS_FERMI", r_val=transport_env%params%eps_fermi)
211 0 : CALL section_vals_val_get(transport_section, "ENERGY_INTERVAL", r_val=transport_env%params%energy_interval)
212 0 : CALL section_vals_val_get(transport_section, "MIN_INTERVAL", r_val=transport_env%params%min_interval)
213 0 : CALL section_vals_val_get(transport_section, "TEMPERATURE", r_val=transport_env%params%temperature)
214 0 : CALL section_vals_val_get(transport_section, "DENSITY_MIXING", r_val=transport_env%params%dens_mixing)
215 0 : CALL section_vals_val_get(transport_section, "CSR_SCREENING", l_val=transport_env%csr_screening)
216 :
217 : ! logical*1 to logical*4 , l_val is logical*1 and c_bool is equivalent to logical*4
218 0 : CALL section_vals_val_get(transport_section, "OBC_EQUILIBRIUM", l_val=obc_equilibrium)
219 0 : IF (obc_equilibrium) THEN
220 0 : transport_env%params%obc_equilibrium = .TRUE.
221 : ELSE
222 0 : transport_env%params%obc_equilibrium = .FALSE.
223 : END IF
224 :
225 0 : CALL section_vals_val_get(transport_section, "CUTOUT", i_vals=i_vals)
226 0 : transport_env%params%cutout = i_vals
227 :
228 : CALL section_vals_val_get(beyn_section, "TASKS_PER_INTEGRATION_POINT", &
229 0 : i_val=transport_env%params%tasks_per_integration_point)
230 0 : CALL section_vals_val_get(beyn_section, "N_POINTS_BEYN", i_val=transport_env%params%n_points_beyn)
231 0 : CALL section_vals_val_get(beyn_section, "N_RAND", r_val=transport_env%params%n_rand_beyn)
232 0 : CALL section_vals_val_get(beyn_section, "N_RAND_CC", r_val=transport_env%params%n_rand_cc_beyn)
233 0 : CALL section_vals_val_get(beyn_section, "SVD_CUTOFF", r_val=transport_env%params%svd_cutoff)
234 0 : CALL section_vals_val_get(beyn_section, "ONE_CIRCLE", l_val=one_circle)
235 0 : IF (one_circle) THEN
236 0 : transport_env%params%ncrc_beyn = 1
237 : ELSE
238 0 : transport_env%params%ncrc_beyn = 2
239 : END IF
240 :
241 0 : CALL section_vals_val_get(pexsi_section, "ORDERING", i_val=transport_env%params%ordering)
242 0 : CALL section_vals_val_get(pexsi_section, "ROW_ORDERING", i_val=transport_env%params%row_ordering)
243 0 : CALL section_vals_val_get(pexsi_section, "VERBOSITY", i_val=transport_env%params%verbosity)
244 0 : CALL section_vals_val_get(pexsi_section, "NP_SYMB_FACT", i_val=transport_env%params%pexsi_np_symb_fact)
245 :
246 0 : IF (contact_explicit) THEN
247 0 : transport_env%params%num_contacts = n_contacts
248 0 : stride_contacts = 5
249 0 : transport_env%params%stride_contacts = stride_contacts
250 0 : ALLOCATE (transport_env%contacts_data(stride_contacts*n_contacts))
251 :
252 0 : DO i = 1, n_contacts
253 0 : CALL section_vals_val_get(contact_section, "BANDWIDTH", i_rep_section=i, i_val=contact_bandwidth)
254 0 : CALL section_vals_val_get(contact_section, "START", i_rep_section=i, i_val=contact_start)
255 0 : CALL section_vals_val_get(contact_section, "N_ATOMS", i_rep_section=i, i_val=contact_natoms)
256 0 : CALL section_vals_val_get(contact_section, "INJECTION_SIGN", i_rep_section=i, i_val=contact_injsign)
257 0 : CALL section_vals_val_get(contact_section, "INJECTING_CONTACT", i_rep_section=i, l_val=injecting_contact)
258 :
259 0 : IF (contact_natoms .LE. 0) CPABORT("Number of atoms in contact region needs to be defined.")
260 :
261 0 : transport_env%contacts_data((i - 1)*stride_contacts + 1) = contact_bandwidth
262 0 : transport_env%contacts_data((i - 1)*stride_contacts + 2) = contact_start - 1 ! C indexing
263 0 : transport_env%contacts_data((i - 1)*stride_contacts + 3) = contact_natoms
264 0 : transport_env%contacts_data((i - 1)*stride_contacts + 4) = contact_injsign
265 :
266 0 : IF (injecting_contact) THEN
267 0 : transport_env%contacts_data((i - 1)*stride_contacts + 5) = 1
268 : ELSE
269 0 : transport_env%contacts_data((i - 1)*stride_contacts + 5) = 0
270 : END IF
271 : END DO
272 0 : transport_env%params%contacts_data = C_LOC(transport_env%contacts_data(1))
273 : ELSE
274 0 : CPABORT("No contact region is defined.")
275 : END IF
276 :
277 0 : CALL timestop(handle)
278 :
279 0 : END SUBROUTINE transport_init_read_input
280 :
281 : ! **************************************************************************************************
282 : !> \brief initializes the transport environment
283 : !> \param ks_env ...
284 : !> \param[inout] transport_env the transport env to be initialized
285 : !> \param[in] template_matrix template matrix to keep the sparsity of matrices fixed
286 : !> \author Mohammad Hossein Bani-Hashemian
287 : ! **************************************************************************************************
288 0 : SUBROUTINE transport_initialize(ks_env, transport_env, template_matrix)
289 : TYPE(qs_ks_env_type), POINTER :: ks_env
290 : TYPE(transport_env_type), INTENT(INOUT) :: transport_env
291 : TYPE(dbcsr_type), INTENT(IN) :: template_matrix
292 :
293 : CHARACTER(len=*), PARAMETER :: routineN = 'transport_initialize'
294 :
295 : INTEGER :: handle, numnodes, unit_nr
296 : TYPE(cp_logger_type), POINTER :: logger
297 :
298 0 : CALL timeset(routineN, handle)
299 :
300 0 : CALL cite_reference(Bruck2014)
301 :
302 0 : logger => cp_get_default_logger()
303 0 : IF (logger%para_env%is_source()) THEN
304 0 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
305 : ELSE
306 0 : unit_nr = -1
307 : END IF
308 :
309 0 : numnodes = logger%para_env%num_pe
310 :
311 0 : IF (dbcsr_has_symmetry(template_matrix)) THEN
312 0 : CALL dbcsr_copy(transport_env%template_matrix_sym, template_matrix)
313 0 : CALL dbcsr_desymmetrize(transport_env%template_matrix_sym, transport_env%template_matrix_nosym)
314 : ELSE
315 0 : CALL dbcsr_copy(transport_env%template_matrix_nosym, template_matrix)
316 0 : CALL dbcsr_copy(transport_env%template_matrix_sym, template_matrix)
317 : END IF
318 :
319 0 : ALLOCATE (transport_env%dm_imag)
320 : CALL dbcsr_create(transport_env%dm_imag, "imaginary DM", &
321 0 : template=template_matrix, matrix_type=dbcsr_type_no_symmetry)
322 0 : CALL dbcsr_set(transport_env%dm_imag, 0.0_dp)
323 :
324 : CALL dbcsr_create(transport_env%csr_sparsity, "CSR sparsity", &
325 : template=transport_env%template_matrix_sym, &
326 0 : data_type=dbcsr_type_real_8)
327 0 : CALL dbcsr_copy(transport_env%csr_sparsity, transport_env%template_matrix_sym)
328 :
329 0 : CALL cp_dbcsr_to_csr_screening(ks_env, transport_env%csr_sparsity)
330 :
331 0 : IF (.NOT. transport_env%csr_screening) CALL dbcsr_set(transport_env%csr_sparsity, 1.0_dp)
332 : CALL dbcsr_csr_create_from_dbcsr(transport_env%template_matrix_nosym, &
333 : transport_env%s_matrix, &
334 : dbcsr_csr_dbcsr_blkrow_dist, &
335 : csr_sparsity=transport_env%csr_sparsity, &
336 0 : numnodes=numnodes)
337 :
338 0 : CALL dbcsr_csr_print_sparsity(transport_env%s_matrix, unit_nr)
339 :
340 0 : CALL dbcsr_convert_dbcsr_to_csr(transport_env%template_matrix_nosym, transport_env%s_matrix)
341 :
342 0 : CALL dbcsr_csr_create(transport_env%ks_matrix, transport_env%s_matrix)
343 0 : CALL dbcsr_csr_create(transport_env%p_matrix, transport_env%s_matrix)
344 0 : CALL dbcsr_csr_create(transport_env%imagp_matrix, transport_env%s_matrix)
345 :
346 0 : CALL timestop(handle)
347 :
348 0 : END SUBROUTINE transport_initialize
349 :
350 : ! **************************************************************************************************
351 : !> \brief SCF calcualtion with an externally evaluated density matrix
352 : !> \param[inout] transport_env transport environment
353 : !> \param[in] matrix_s DBCSR overlap matrix
354 : !> \param[in] matrix_ks DBCSR Kohn-Sham matrix
355 : !> \param[inout] matrix_p DBCSR density matrix
356 : !> \param[in] nelectron_spin number of electrons
357 : !> \param[in] natoms number of atoms
358 : !> \param[in] energy_diff scf energy difference
359 : !> \param[in] iscf the current scf iteration
360 : !> \param[in] extra_scf whether or not an extra scf step will be performed
361 : !> \par History
362 : !> 12.2012 created [Hossein Bani-Hashemian]
363 : !> 12.2014 revised [Hossein Bani-Hashemian]
364 : !> \author Mohammad Hossein Bani-Hashemian
365 : ! **************************************************************************************************
366 0 : SUBROUTINE external_scf_method(transport_env, matrix_s, matrix_ks, matrix_p, &
367 : nelectron_spin, natoms, energy_diff, iscf, extra_scf)
368 :
369 : TYPE(transport_env_type), INTENT(INOUT) :: transport_env
370 : TYPE(dbcsr_type), INTENT(IN) :: matrix_s, matrix_ks
371 : TYPE(dbcsr_type), INTENT(INOUT) :: matrix_p
372 : INTEGER, INTENT(IN) :: nelectron_spin, natoms
373 : REAL(dp), INTENT(IN) :: energy_diff
374 : INTEGER, INTENT(IN) :: iscf
375 : LOGICAL, INTENT(IN) :: extra_scf
376 :
377 : CHARACTER(len=*), PARAMETER :: routineN = 'external_scf_method'
378 :
379 : TYPE(cp2k_csr_interop_type) :: imagp_mat, ks_mat, p_mat, s_mat
380 :
381 : PROCEDURE(c_scf_routine), POINTER :: c_method
382 : INTEGER :: handle
383 :
384 0 : CALL timeset(routineN, handle)
385 :
386 0 : CALL C_F_PROCPOINTER(transport_env%ext_c_method_ptr, c_method)
387 0 : IF (.NOT. C_ASSOCIATED(transport_env%ext_c_method_ptr)) &
388 : CALL cp_abort(__LOCATION__, &
389 : "MISSING C/C++ ROUTINE: The TRANSPORT section is meant to be used together with an external "// &
390 0 : "program, e.g. the quantum transport code OMEN, that provides CP2K with a density matrix.")
391 :
392 0 : transport_env%params%n_occ = nelectron_spin
393 0 : transport_env%params%n_atoms = natoms
394 0 : transport_env%params%energy_diff = energy_diff
395 0 : transport_env%params%evoltfactor = evolt
396 0 : transport_env%params%e_charge = e_charge
397 0 : transport_env%params%boltzmann = boltzmann
398 0 : transport_env%params%h_bar = h_bar
399 0 : transport_env%params%iscf = iscf
400 0 : transport_env%params%extra_scf = LOGICAL(extra_scf, C_BOOL) ! C_BOOL and Fortran logical don't have the same size
401 :
402 0 : CALL csr_interop_nullify(s_mat)
403 0 : CALL csr_interop_nullify(ks_mat)
404 0 : CALL csr_interop_nullify(p_mat)
405 0 : CALL csr_interop_nullify(imagp_mat)
406 :
407 0 : CALL dbcsr_copy(transport_env%template_matrix_sym, matrix_s, keep_sparsity=.TRUE.)
408 0 : CALL convert_dbcsr_to_csr_interop(transport_env%template_matrix_sym, transport_env%s_matrix, s_mat)
409 :
410 0 : CALL dbcsr_copy(transport_env%template_matrix_sym, matrix_ks, keep_sparsity=.TRUE.)
411 0 : CALL convert_dbcsr_to_csr_interop(transport_env%template_matrix_sym, transport_env%ks_matrix, ks_mat)
412 :
413 0 : CALL dbcsr_copy(transport_env%template_matrix_sym, matrix_p, keep_sparsity=.TRUE.)
414 0 : CALL convert_dbcsr_to_csr_interop(transport_env%template_matrix_sym, transport_env%p_matrix, p_mat)
415 :
416 0 : CALL dbcsr_copy(transport_env%template_matrix_sym, matrix_s, keep_sparsity=.TRUE.)
417 0 : CALL convert_dbcsr_to_csr_interop(transport_env%template_matrix_sym, transport_env%imagp_matrix, imagp_mat)
418 :
419 0 : CALL c_method(transport_env%params, s_mat, ks_mat, p_mat, imagp_mat)
420 :
421 0 : CALL convert_csr_interop_to_dbcsr(p_mat, transport_env%p_matrix, transport_env%template_matrix_nosym)
422 0 : CALL dbcsr_copy(matrix_p, transport_env%template_matrix_nosym)
423 :
424 0 : CALL convert_csr_interop_to_dbcsr(imagp_mat, transport_env%imagp_matrix, transport_env%template_matrix_nosym)
425 0 : CALL dbcsr_copy(transport_env%dm_imag, transport_env%template_matrix_nosym)
426 :
427 0 : CALL timestop(handle)
428 :
429 0 : END SUBROUTINE external_scf_method
430 :
431 : ! **************************************************************************************************
432 : !> \brief converts a DBCSR matrix to a C-interoperable CSR matrix
433 : !> \param[in] dbcsr_mat DBCSR matrix to be converted
434 : !> \param[inout] csr_mat auxiliary CSR matrix
435 : !> \param[inout] csr_interop_mat C-interoperable CSR matrix
436 : !> \author Mohammad Hossein Bani-Hashemian
437 : ! **************************************************************************************************
438 0 : SUBROUTINE convert_dbcsr_to_csr_interop(dbcsr_mat, csr_mat, csr_interop_mat)
439 :
440 : TYPE(dbcsr_type), INTENT(IN) :: dbcsr_mat
441 : TYPE(dbcsr_csr_type), INTENT(INOUT) :: csr_mat
442 : TYPE(cp2k_csr_interop_type), INTENT(INOUT) :: csr_interop_mat
443 :
444 : CHARACTER(LEN=*), PARAMETER :: routineN = 'convert_dbcsr_to_csr_interop'
445 :
446 : INTEGER :: handle
447 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: nrows_local_all, first_row_all
448 0 : INTEGER(C_INT), DIMENSION(:), POINTER :: colind_local, rowptr_local, nzerow_local
449 0 : REAL(C_DOUBLE), DIMENSION(:), POINTER :: nzvals_local
450 : TYPE(cp_logger_type), POINTER :: logger
451 :
452 0 : CALL timeset(routineN, handle)
453 :
454 0 : logger => cp_get_default_logger()
455 :
456 : ! dbcsr to csr
457 0 : CALL dbcsr_convert_dbcsr_to_csr(dbcsr_mat, csr_mat)
458 :
459 : ! csr to csr_interop
460 0 : rowptr_local => csr_mat%rowptr_local
461 0 : colind_local => csr_mat%colind_local
462 0 : nzerow_local => csr_mat%nzerow_local
463 0 : nzvals_local => csr_mat%nzval_local%r_dp ! support real double percision for now
464 :
465 0 : IF (SIZE(rowptr_local) .EQ. 0) THEN
466 0 : csr_interop_mat%rowptr_local = C_NULL_PTR
467 : ELSE
468 0 : csr_interop_mat%rowptr_local = C_LOC(rowptr_local(1))
469 : END IF
470 :
471 0 : IF (SIZE(colind_local) .EQ. 0) THEN
472 0 : csr_interop_mat%colind_local = C_NULL_PTR
473 : ELSE
474 0 : csr_interop_mat%colind_local = C_LOC(colind_local(1))
475 : END IF
476 :
477 0 : IF (SIZE(nzerow_local) .EQ. 0) THEN
478 0 : csr_interop_mat%nzerow_local = C_NULL_PTR
479 : ELSE
480 0 : csr_interop_mat%nzerow_local = C_LOC(nzerow_local(1))
481 : END IF
482 :
483 0 : IF (SIZE(nzvals_local) .EQ. 0) THEN
484 0 : csr_interop_mat%nzvals_local = C_NULL_PTR
485 : ELSE
486 0 : csr_interop_mat%nzvals_local = C_LOC(nzvals_local(1))
487 : END IF
488 :
489 : ASSOCIATE (mp_group => logger%para_env, mepos => logger%para_env%mepos, num_pe => logger%para_env%num_pe)
490 0 : ALLOCATE (nrows_local_all(0:num_pe - 1), first_row_all(0:num_pe - 1))
491 0 : CALL mp_group%allgather(csr_mat%nrows_local, nrows_local_all)
492 0 : CALL cumsum_i(nrows_local_all, first_row_all)
493 :
494 0 : IF (mepos .EQ. 0) THEN
495 0 : csr_interop_mat%first_row = 0
496 : ELSE
497 0 : csr_interop_mat%first_row = first_row_all(mepos - 1)
498 : END IF
499 : END ASSOCIATE
500 0 : csr_interop_mat%nrows_total = csr_mat%nrows_total
501 0 : csr_interop_mat%ncols_total = csr_mat%ncols_total
502 0 : csr_interop_mat%nze_local = csr_mat%nze_local
503 0 : IF (csr_mat%nze_total > HUGE(csr_interop_mat%nze_total)) THEN
504 0 : CPABORT("overflow in nze")
505 : END IF
506 0 : csr_interop_mat%nze_total = INT(csr_mat%nze_total, KIND=KIND(csr_interop_mat%nze_total))
507 0 : csr_interop_mat%nrows_local = csr_mat%nrows_local
508 0 : csr_interop_mat%data_type = csr_mat%nzval_local%data_type
509 :
510 0 : CALL timestop(handle)
511 :
512 : CONTAINS
513 : ! **************************************************************************************************
514 : !> \brief cumulative sum of a 1d array of integers
515 : !> \param[in] arr input array
516 : !> \param[out] cumsum cumulative sum of the input array
517 : ! **************************************************************************************************
518 0 : SUBROUTINE cumsum_i(arr, cumsum)
519 : INTEGER, DIMENSION(:), INTENT(IN) :: arr
520 : INTEGER, DIMENSION(SIZE(arr)), INTENT(OUT) :: cumsum
521 :
522 : INTEGER :: i
523 :
524 0 : cumsum(1) = arr(1)
525 0 : DO i = 2, SIZE(arr)
526 0 : cumsum(i) = cumsum(i - 1) + arr(i)
527 : END DO
528 0 : END SUBROUTINE cumsum_i
529 :
530 : END SUBROUTINE convert_dbcsr_to_csr_interop
531 :
532 : ! **************************************************************************************************
533 : !> \brief converts a C-interoperable CSR matrix to a DBCSR matrix
534 : !> \param[in] csr_interop_mat C-interoperable CSR matrix
535 : !> \param[inout] csr_mat auxiliary CSR matrix
536 : !> \param[inout] dbcsr_mat DBCSR matrix
537 : !> \author Mohammad Hossein Bani-Hashemian
538 : ! **************************************************************************************************
539 0 : SUBROUTINE convert_csr_interop_to_dbcsr(csr_interop_mat, csr_mat, dbcsr_mat)
540 :
541 : TYPE(cp2k_csr_interop_type), INTENT(IN) :: csr_interop_mat
542 : TYPE(dbcsr_csr_type), INTENT(INOUT) :: csr_mat
543 : TYPE(dbcsr_type), INTENT(INOUT) :: dbcsr_mat
544 :
545 : CHARACTER(LEN=*), PARAMETER :: routineN = 'convert_csr_interop_to_dbcsr'
546 :
547 : INTEGER :: data_type, handle, ncols_total, &
548 : nrows_local, nrows_total, nze_local, &
549 : nze_total
550 0 : INTEGER, DIMENSION(:), POINTER :: colind_local, nzerow_local, rowptr_local
551 0 : REAL(dp), DIMENSION(:), POINTER :: nzvals_local
552 :
553 0 : CALL timeset(routineN, handle)
554 :
555 : ! csr_interop to csr
556 : CALL csr_interop_matrix_get_info(csr_interop_mat, &
557 : nrows_total=nrows_total, ncols_total=ncols_total, nze_local=nze_local, &
558 : nze_total=nze_total, nrows_local=nrows_local, data_type=data_type, &
559 : rowptr_local=rowptr_local, colind_local=colind_local, &
560 0 : nzerow_local=nzerow_local, nzvals_local=nzvals_local)
561 :
562 0 : csr_mat%nrows_total = nrows_total
563 0 : csr_mat%ncols_total = ncols_total
564 0 : csr_mat%nze_local = nze_local
565 0 : csr_mat%nze_total = nze_total
566 0 : csr_mat%nrows_local = nrows_local
567 0 : csr_mat%nzval_local%data_type = data_type
568 :
569 0 : csr_mat%rowptr_local = rowptr_local
570 0 : csr_mat%colind_local = colind_local
571 0 : csr_mat%nzerow_local = nzerow_local
572 0 : csr_mat%nzval_local%r_dp = nzvals_local
573 :
574 : ! csr to dbcsr
575 0 : CALL dbcsr_convert_csr_to_dbcsr(dbcsr_mat, csr_mat)
576 :
577 0 : CALL timestop(handle)
578 :
579 0 : END SUBROUTINE convert_csr_interop_to_dbcsr
580 :
581 : ! **************************************************************************************************
582 : !> \brief extraxts zeff (effective nuclear charges per atom) and nsgf (the size
583 : !> of spherical Gaussian basis functions per atom) from qs_env and initializes
584 : !> the corresponding arrays in transport_env%params
585 : !> \param[in] qs_env qs environment
586 : !> \param[inout] transport_env transport environment
587 : !> \author Mohammad Hossein Bani-Hashemian
588 : ! **************************************************************************************************
589 0 : SUBROUTINE transport_set_contact_params(qs_env, transport_env)
590 : TYPE(qs_environment_type), POINTER :: qs_env
591 : TYPE(transport_env_type), INTENT(INOUT) :: transport_env
592 :
593 : INTEGER :: i, iat, ikind, natom, nkind
594 0 : INTEGER, DIMENSION(:), POINTER :: atom_list
595 : REAL(KIND=dp) :: zeff
596 0 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
597 : TYPE(atomic_kind_type), POINTER :: atomic_kind
598 0 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
599 0 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
600 :
601 0 : CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
602 : CALL get_qs_env(qs_env, particle_set=particle_set, &
603 : qs_kind_set=qs_kind_set, &
604 0 : atomic_kind_set=atomic_kind_set)
605 :
606 0 : ALLOCATE (transport_env%nsgf(natom))
607 0 : ALLOCATE (transport_env%zeff(natom))
608 0 : CALL get_particle_set(particle_set, qs_kind_set, nsgf=transport_env%nsgf)
609 :
610 : ! reference charges
611 0 : DO ikind = 1, nkind
612 0 : CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff)
613 0 : atomic_kind => atomic_kind_set(ikind)
614 0 : CALL get_atomic_kind(atomic_kind, atom_list=atom_list)
615 0 : DO iat = 1, SIZE(atom_list)
616 0 : i = atom_list(iat)
617 0 : transport_env%zeff(i) = zeff
618 : END DO
619 : END DO
620 :
621 0 : IF (natom .EQ. 0) THEN
622 0 : transport_env%params%nsgf = C_NULL_PTR
623 0 : transport_env%params%zeff = C_NULL_PTR
624 : ELSE
625 0 : transport_env%params%nsgf = C_LOC(transport_env%nsgf(1))
626 0 : transport_env%params%zeff = C_LOC(transport_env%zeff(1))
627 : END IF
628 :
629 0 : END SUBROUTINE transport_set_contact_params
630 :
631 : !**************************************************************************************************
632 : !> \brief evaluates current density using the imaginary part of the density matrix and prints it
633 : !> into cube files
634 : !> \param[in] input the input
635 : !> \param[in] qs_env qs environment
636 : !> \author Mohammad Hossein Bani-Hashemian
637 : ! **************************************************************************************************
638 10821 : SUBROUTINE transport_current(input, qs_env)
639 : TYPE(section_vals_type), POINTER :: input
640 : TYPE(qs_environment_type), POINTER :: qs_env
641 :
642 : CHARACTER(len=*), PARAMETER :: routineN = 'transport_current'
643 :
644 : CHARACTER(len=14) :: ext
645 : CHARACTER(len=2) :: sdir
646 : INTEGER :: dir, handle, unit_nr
647 : LOGICAL :: do_current_cube, do_transport, mpi_io
648 : TYPE(cp_logger_type), POINTER :: logger
649 : TYPE(current_env_type) :: current_env
650 : TYPE(dbcsr_type), POINTER :: zero
651 : TYPE(dft_control_type), POINTER :: dft_control
652 : TYPE(particle_list_type), POINTER :: particles
653 : TYPE(pw_c1d_gs_type) :: gs
654 10821 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
655 : TYPE(pw_env_type), POINTER :: pw_env
656 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
657 : TYPE(pw_r3d_rs_type) :: rs
658 10821 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
659 : TYPE(qs_rho_type), POINTER :: rho
660 : TYPE(qs_subsys_type), POINTER :: subsys
661 : TYPE(section_vals_type), POINTER :: dft_section
662 : TYPE(transport_env_type), POINTER :: transport_env
663 :
664 10821 : CALL timeset(routineN, handle)
665 :
666 10821 : logger => cp_get_default_logger()
667 10821 : dft_section => section_vals_get_subs_vals(input, "DFT")
668 : CALL get_qs_env(qs_env=qs_env, &
669 : subsys=subsys, &
670 : pw_env=pw_env, &
671 : transport_env=transport_env, &
672 : do_transport=do_transport, &
673 : dft_control=dft_control, &
674 10821 : rho=rho)
675 10821 : CALL qs_subsys_get(subsys, particles=particles)
676 10821 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
677 10821 : CALL qs_rho_get(rho, rho_r=rho_r, rho_g=rho_g)
678 :
679 : do_current_cube = BTEST(cp_print_key_should_output(logger%iter_info, input, &
680 10821 : "DFT%TRANSPORT%PRINT%CURRENT"), cp_p_file)
681 :
682 : ! check if transport is active i.e. the imaginary dm has been actually passed to cp2k by the external code
683 10821 : IF (do_transport) THEN
684 0 : IF (C_ASSOCIATED(transport_env%ext_c_method_ptr)) THEN
685 :
686 : ! calculate_jrho_resp uses sab_all which is not associated in DFTB environment
687 0 : IF (dft_control%qs_control%dftb) CPABORT("Current is not available for DFTB.")
688 0 : IF (dft_control%qs_control%xtb) CPABORT("Current is not available for xTB.")
689 :
690 : ! now do the current
691 0 : IF (do_current_cube) THEN
692 0 : current_env%gauge = -1
693 0 : current_env%gauge_init = .FALSE.
694 :
695 0 : CALL auxbas_pw_pool%create_pw(rs)
696 0 : CALL auxbas_pw_pool%create_pw(gs)
697 :
698 : NULLIFY (zero)
699 0 : ALLOCATE (zero)
700 0 : CALL dbcsr_create(zero, template=transport_env%dm_imag)
701 0 : CALL dbcsr_copy(zero, transport_env%dm_imag)
702 0 : CALL dbcsr_set(zero, 0.0_dp)
703 :
704 0 : DO dir = 1, 3
705 0 : CALL pw_zero(rs)
706 0 : CALL pw_zero(gs)
707 : CALL calculate_jrho_resp(zero, transport_env%dm_imag, &
708 : zero, zero, dir, dir, rs, gs, qs_env, current_env, &
709 0 : retain_rsgrid=.TRUE.)
710 0 : SELECT CASE (dir)
711 : CASE (1)
712 0 : sdir = "-x"
713 : CASE (2)
714 0 : sdir = "-y"
715 : CASE (3)
716 0 : sdir = "-z"
717 : END SELECT
718 0 : ext = sdir//".cube"
719 0 : mpi_io = .TRUE.
720 : unit_nr = cp_print_key_unit_nr(logger, dft_section, "TRANSPORT%PRINT%CURRENT", &
721 : extension=ext, file_status="REPLACE", file_action="WRITE", &
722 0 : log_filename=.FALSE., mpi_io=mpi_io)
723 : CALL cp_pw_to_cube(rs, unit_nr, "Transport current", particles=particles, &
724 : stride=section_get_ivals(dft_section, "TRANSPORT%PRINT%CURRENT%STRIDE"), &
725 0 : mpi_io=mpi_io)
726 : CALL cp_print_key_finished_output(unit_nr, logger, dft_section, "TRANSPORT%PRINT%CURRENT", &
727 0 : mpi_io=mpi_io)
728 : END DO
729 :
730 0 : CALL dbcsr_deallocate_matrix(zero)
731 0 : CALL auxbas_pw_pool%give_back_pw(rs)
732 0 : CALL auxbas_pw_pool%give_back_pw(gs)
733 : END IF
734 : END IF
735 :
736 : END IF
737 :
738 10821 : CALL timestop(handle)
739 :
740 789933 : END SUBROUTINE transport_current
741 :
742 : !**************************************************************************************************
743 : !> \brief post scf calculations for transport
744 : !> \param[in] qs_env qs environment
745 : !> \author Mohammad Hossein Bani-Hashemian
746 : ! **************************************************************************************************
747 10821 : SUBROUTINE qs_scf_post_transport(qs_env)
748 : TYPE(qs_environment_type), POINTER :: qs_env
749 :
750 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_transport'
751 :
752 : INTEGER :: handle
753 : TYPE(section_vals_type), POINTER :: input
754 :
755 10821 : CALL timeset(routineN, handle)
756 :
757 10821 : NULLIFY (input)
758 :
759 : CALL get_qs_env(qs_env=qs_env, &
760 10821 : input=input)
761 :
762 10821 : CALL transport_current(input, qs_env)
763 :
764 10821 : CALL timestop(handle)
765 :
766 10821 : END SUBROUTINE qs_scf_post_transport
767 :
768 : END MODULE transport
|
