(git:374b731)
Loading...
Searching...
No Matches
qs_fermi_contact.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 Distribution of the Fermi contact integral matrix.
10!> \par History
11!> \author VW (27.02.2009)
12! **************************************************************************************************
13MODULE qs_fermi_contact
14
15 USE ai_fermi_contact, ONLY: fermi_contact
16 USE basis_set_types, ONLY: gto_basis_set_p_type,&
17 gto_basis_set_type
18 USE block_p_types, ONLY: block_p_type
19 USE cell_types, ONLY: cell_type,&
20 pbc
21 USE cp_dbcsr_output, ONLY: cp_dbcsr_write_sparse_matrix
22 USE cp_log_handling, ONLY: cp_get_default_logger,&
23 cp_logger_type
24 USE cp_output_handling, ONLY: cp_p_file,&
25 cp_print_key_finished_output,&
26 cp_print_key_should_output,&
27 cp_print_key_unit_nr
28 USE dbcsr_api, ONLY: dbcsr_get_block_p,&
29 dbcsr_p_type
30 USE input_section_types, ONLY: section_vals_val_get
31 USE kinds, ONLY: dp
32 USE message_passing, ONLY: mp_para_env_type
33 USE orbital_pointers, ONLY: init_orbital_pointers,&
34 ncoset
35 USE particle_types, ONLY: particle_type
36 USE qs_environment_types, ONLY: get_qs_env,&
37 qs_environment_type
38 USE qs_kind_types, ONLY: get_qs_kind,&
39 get_qs_kind_set,&
40 qs_kind_type
41 USE qs_neighbor_list_types, ONLY: get_iterator_info,&
42 neighbor_list_iterate,&
43 neighbor_list_iterator_create,&
44 neighbor_list_iterator_p_type,&
45 neighbor_list_iterator_release,&
46 neighbor_list_set_p_type
47#include "./base/base_uses.f90"
48
49 IMPLICIT NONE
50
51 PRIVATE
52
53! *** Global parameters ***
54
55 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_fermi_contact'
56
57! *** Public subroutines ***
58
59 PUBLIC :: build_fermi_contact_matrix
60
61CONTAINS
62
63! **************************************************************************************************
64!> \brief Calculation of the Fermi contact matrix over Cartesian
65!> Gaussian functions.
66!> \param qs_env ...
67!> \param matrix_fc ...
68!> \param rc ...
69!> \date 27.02.2009
70!> \author VW
71!> \version 1.0
72! **************************************************************************************************
73
74 SUBROUTINE build_fermi_contact_matrix(qs_env, matrix_fc, rc)
75
76 TYPE(qs_environment_type), POINTER :: qs_env
77 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_fc
78 REAL(dp), DIMENSION(3), INTENT(IN) :: rc
79
80 CHARACTER(len=*), PARAMETER :: routinen = 'build_fermi_contact_matrix'
81
82 INTEGER :: after, handle, iatom, icol, ikind, inode, irow, iset, iw, jatom, jkind, jset, &
83 last_jatom, ldai, ldfc, maxco, maxlgto, maxsgf, natom, ncoa, ncob, nkind, nseta, nsetb, &
84 sgfa, sgfb
85 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
86 npgfb, nsgfa, nsgfb
87 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
88 LOGICAL :: found, new_atom_b, omit_headers
89 REAL(kind=dp) :: dab, rab2
90 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: fcab, work
91 REAL(kind=dp), DIMENSION(3) :: ra, rab, rac, rb, rbc
92 REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
93 REAL(kind=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb
94 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: fcint
95 TYPE(cell_type), POINTER :: cell
96 TYPE(cp_logger_type), POINTER :: logger
97 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
98 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
99 TYPE(mp_para_env_type), POINTER :: para_env
100 TYPE(neighbor_list_iterator_p_type), &
101 DIMENSION(:), POINTER :: nl_iterator
102 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
103 POINTER :: sab_orb
104 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
105 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
106 TYPE(qs_kind_type), POINTER :: qs_kind
107
108 CALL timeset(routinen, handle)
109
110 NULLIFY (cell, sab_orb, qs_kind_set, particle_set, para_env)
111 NULLIFY (logger)
112
113 logger => cp_get_default_logger()
114
115 CALL get_qs_env(qs_env=qs_env, &
116 qs_kind_set=qs_kind_set, &
117 particle_set=particle_set, &
118 para_env=para_env, &
119 sab_orb=sab_orb, &
120 cell=cell)
121
122 nkind = SIZE(qs_kind_set)
123 natom = SIZE(particle_set)
124
125 ! *** Allocate work storage ***
126
127 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
128 maxco=maxco, &
129 maxlgto=maxlgto, &
130 maxsgf=maxsgf)
131
132 ldai = ncoset(maxlgto)
133 CALL init_orbital_pointers(ldai)
134
135 ldfc = maxco
136 ALLOCATE (fcab(ldfc, ldfc))
137 fcab(:, :) = 0.0_dp
138
139 ALLOCATE (work(maxco, maxsgf))
140 work(:, :) = 0.0_dp
141
142 ALLOCATE (fcint(1))
143 NULLIFY (fcint(1)%block)
144
145 ALLOCATE (basis_set_list(nkind))
146 DO ikind = 1, nkind
147 qs_kind => qs_kind_set(ikind)
148 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a)
149 IF (ASSOCIATED(basis_set_a)) THEN
150 basis_set_list(ikind)%gto_basis_set => basis_set_a
151 ELSE
152 NULLIFY (basis_set_list(ikind)%gto_basis_set)
153 END IF
154 END DO
155 CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
156 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
157 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
158 iatom=iatom, jatom=jatom, r=rab)
159 basis_set_a => basis_set_list(ikind)%gto_basis_set
160 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
161 basis_set_b => basis_set_list(jkind)%gto_basis_set
162 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
163 ra = pbc(particle_set(iatom)%r, cell)
164 ! basis ikind
165 first_sgfa => basis_set_a%first_sgf
166 la_max => basis_set_a%lmax
167 la_min => basis_set_a%lmin
168 npgfa => basis_set_a%npgf
169 nseta = basis_set_a%nset
170 nsgfa => basis_set_a%nsgf_set
171 rpgfa => basis_set_a%pgf_radius
172 set_radius_a => basis_set_a%set_radius
173 sphi_a => basis_set_a%sphi
174 zeta => basis_set_a%zet
175 ! basis jkind
176 first_sgfb => basis_set_b%first_sgf
177 lb_max => basis_set_b%lmax
178 lb_min => basis_set_b%lmin
179 npgfb => basis_set_b%npgf
180 nsetb = basis_set_b%nset
181 nsgfb => basis_set_b%nsgf_set
182 rpgfb => basis_set_b%pgf_radius
183 set_radius_b => basis_set_b%set_radius
184 sphi_b => basis_set_b%sphi
185 zetb => basis_set_b%zet
186
187 IF (inode == 1) last_jatom = 0
188
189 rb = rab + ra
190 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
191 dab = sqrt(rab2)
192 rac = pbc(ra, rc, cell)
193 rbc = rac - rab
194
195 IF (jatom /= last_jatom) THEN
196 new_atom_b = .true.
197 last_jatom = jatom
198 ELSE
199 new_atom_b = .false.
200 END IF
201
202 IF (new_atom_b) THEN
203 IF (iatom <= jatom) THEN
204 irow = iatom
205 icol = jatom
206 ELSE
207 irow = jatom
208 icol = iatom
209 END IF
210
211 NULLIFY (fcint(1)%block)
212 CALL dbcsr_get_block_p(matrix=matrix_fc(1)%matrix, &
213 row=irow, col=icol, block=fcint(1)%block, found=found)
214 END IF
215
216 DO iset = 1, nseta
217
218 ncoa = npgfa(iset)*ncoset(la_max(iset))
219 sgfa = first_sgfa(1, iset)
220
221 DO jset = 1, nsetb
222
223 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
224
225 ncob = npgfb(jset)*ncoset(lb_max(jset))
226 sgfb = first_sgfb(1, jset)
227
228 ! *** Calculate the primitive fermi contact integrals ***
229
230 CALL fermi_contact(la_max(iset), la_min(iset), npgfa(iset), &
231 rpgfa(:, iset), zeta(:, iset), &
232 lb_max(jset), lb_min(jset), npgfb(jset), &
233 rpgfb(:, jset), zetb(:, jset), &
234 rac, rbc, dab, fcab, SIZE(fcab, 1))
235
236 ! *** Contraction step ***
237
238 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
239 1.0_dp, fcab(1, 1), SIZE(fcab, 1), &
240 sphi_b(1, sgfb), SIZE(sphi_b, 1), &
241 0.0_dp, work(1, 1), SIZE(work, 1))
242
243 IF (iatom <= jatom) THEN
244
245 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
246 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
247 work(1, 1), SIZE(work, 1), &
248 1.0_dp, fcint(1)%block(sgfa, sgfb), &
249 SIZE(fcint(1)%block, 1))
250
251 ELSE
252
253 CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
254 1.0_dp, work(1, 1), SIZE(work, 1), &
255 sphi_a(1, sgfa), SIZE(sphi_a, 1), &
256 1.0_dp, fcint(1)%block(sgfb, sgfa), &
257 SIZE(fcint(1)%block, 1))
258
259 END IF
260
261 END DO
262
263 END DO
264
265 END DO
266 CALL neighbor_list_iterator_release(nl_iterator)
267
268 ! *** Release work storage ***
269 DEALLOCATE (basis_set_list)
270
271 DEALLOCATE (fcab)
272
273 DEALLOCATE (work)
274
275 NULLIFY (fcint(1)%block)
276 DEALLOCATE (fcint)
277
278! *** Print the Fermi contact matrix, if requested ***
279
280 IF (btest(cp_print_key_should_output(logger%iter_info, &
281 qs_env%input, "DFT%PRINT%AO_MATRICES/FERMI_CONTACT"), cp_p_file)) THEN
282 iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/FERMI_CONTACT", &
283 extension=".Log")
284 CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
285 CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers)
286 after = min(max(after, 1), 16)
287 CALL cp_dbcsr_write_sparse_matrix(matrix_fc(1)%matrix, 4, after, qs_env, &
288 para_env, output_unit=iw, omit_headers=omit_headers)
289 CALL cp_print_key_finished_output(iw, logger, qs_env%input, &
290 "DFT%PRINT%AO_MATRICES/FERMI_CONTACT")
291 END IF
292
293 CALL timestop(handle)
294
295 END SUBROUTINE build_fermi_contact_matrix
296
297! **************************************************************************************************
298
299END MODULE qs_fermi_contact
300
dgemm
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Definition grid_cpu_task_list.c:195
cell_types::pbc
Definition cell_types.F:92
ai_fermi_contact
Calculation of the Fermi contact integrals over Cartesian Gaussian-type functions.
Definition ai_fermi_contact.F:29
ai_fermi_contact::fermi_contact
subroutine, public fermi_contact(la_max, la_min, npgfa, rpgfa, zeta, lb_max, lb_min, npgfb, rpgfb, zetb, rac, rbc, dab, fcab, ldfc)
Purpose: Calculation of the two-center Fermi contact integrals 4/3*pi*[a|delta(r-c)|b] over Cartesian...
Definition ai_fermi_contact.F:73
basis_set_types
Definition basis_set_types.F:15
block_p_types
collect pointers to a block of reals
Definition block_p_types.F:14
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cp_dbcsr_output
DBCSR output in CP2K.
Definition cp_dbcsr_output.F:17
cp_dbcsr_output::cp_dbcsr_write_sparse_matrix
subroutine, public cp_dbcsr_write_sparse_matrix(sparse_matrix, before, after, qs_env, para_env, first_row, last_row, first_col, last_col, scale, output_unit, omit_headers)
...
Definition cp_dbcsr_output.F:163
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:803
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1013
cp_output_handling::cp_p_file
integer, parameter, public cp_p_file
Definition cp_output_handling.F:103
cp_output_handling::cp_print_key_should_output
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Definition cp_output_handling.F:345
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1040
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition orbital_pointers.F:15
orbital_pointers::init_orbital_pointers
subroutine, public init_orbital_pointers(maxl)
Initialize or update the orbital pointers.
Definition orbital_pointers.F:253
orbital_pointers::ncoset
integer, dimension(:), allocatable, public ncoset
Definition orbital_pointers.F:42
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:502
qs_fermi_contact
Distribution of the Fermi contact integral matrix.
Definition qs_fermi_contact.F:13
qs_fermi_contact::build_fermi_contact_matrix
subroutine, public build_fermi_contact_matrix(qs_env, matrix_fc, rc)
Calculation of the Fermi contact matrix over Cartesian Gaussian functions.
Definition qs_fermi_contact.F:75
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_r3d_rs_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:451
qs_kind_types::get_qs_kind_set
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, basis_type, total_zeff_corr)
Get attributes of an atomic kind set.
Definition qs_kind_types.F:864
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition qs_neighbor_list_types.F:17
qs_neighbor_list_types::neighbor_list_iterator_create
subroutine, public neighbor_list_iterator_create(iterator_set, nl, search, nthread)
Neighbor list iterator functions.
Definition qs_neighbor_list_types.F:165
qs_neighbor_list_types::neighbor_list_iterator_release
subroutine, public neighbor_list_iterator_release(iterator_set)
...
Definition qs_neighbor_list_types.F:251
qs_neighbor_list_types::neighbor_list_iterate
integer function, public neighbor_list_iterate(iterator_set, mepos)
...
Definition qs_neighbor_list_types.F:351
qs_neighbor_list_types::get_iterator_info
subroutine, public get_iterator_info(iterator_set, mepos, ikind, jkind, nkind, ilist, nlist, inode, nnode, iatom, jatom, r, cell)
...
Definition qs_neighbor_list_types.F:549
basis_set_types::gto_basis_set_p_type
Definition basis_set_types.F:93
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:71
block_p_types::block_p_type
Definition block_p_types.F:23
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
particle_types::particle_type
Definition particle_types.F:35
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:215
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:171
qs_neighbor_list_types::neighbor_list_iterator_p_type
Definition qs_neighbor_list_types.F:117
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67