(git:0de0cc2)
qs_elec_field.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 electric field gradient integral matrix.
10 !> \par History
11 !> \author VW (27.02.2009)
12 ! **************************************************************************************************
14 
15  USE ai_elec_field, ONLY: efg
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
23  cp_logger_type
24  USE cp_output_handling, ONLY: cp_p_file,&
28  USE dbcsr_api, ONLY: dbcsr_get_block_p,&
29  dbcsr_p_type
31  USE kinds, ONLY: dp
32  USE message_passing, ONLY: mp_para_env_type
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,&
40  qs_kind_type
44  neighbor_list_iterator_p_type,&
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_elec_field'
56 
57 ! *** Public subroutines ***
58 
59  PUBLIC :: build_efg_matrix
60 
61 CONTAINS
62 
63 ! **************************************************************************************************
64 !> \brief Calculation of the electric field gradient matrix over
65 !> Cartesian Gaussian functions.
66 !> \param qs_env ...
67 !> \param matrix_efg ...
68 !> \param rc ...
69 !> \date 27.02.2009
70 !> \author VW
71 !> \version 1.0
72 ! **************************************************************************************************
73 
74  SUBROUTINE build_efg_matrix(qs_env, matrix_efg, rc)
75 
76  TYPE(qs_environment_type), POINTER :: qs_env
77  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_efg
78  REAL(dp), DIMENSION(3), INTENT(IN) :: rc
79 
80  CHARACTER(len=*), PARAMETER :: routinen = 'build_efg_matrix'
81 
82  INTEGER :: after, handle, i, iatom, icol, ikind, inode, irow, iset, iw, jatom, jkind, jset, &
83  last_jatom, ldai, maxco, maxlgto, maxsgf, natom, ncoa, ncob, neighbor_list_id, nkind, &
84  nseta, nsetb, 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(:, :) :: work
91  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: efgab, rr_work
92  REAL(kind=dp), DIMENSION(3) :: ra, rab, rac, rb, rbc
93  REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
94  REAL(kind=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb
95  TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: efgint
96  TYPE(cell_type), POINTER :: cell
97  TYPE(cp_logger_type), POINTER :: logger
98  TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
99  TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
100  TYPE(mp_para_env_type), POINTER :: para_env
101  TYPE(neighbor_list_iterator_p_type), &
102  DIMENSION(:), POINTER :: nl_iterator
103  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
104  POINTER :: sab_orb
105  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
106  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
107  TYPE(qs_kind_type), POINTER :: qs_kind
108 
109  CALL timeset(routinen, handle)
110 
111  NULLIFY (cell, sab_orb, qs_kind_set, particle_set, para_env)
112  NULLIFY (logger)
113 
114  logger => cp_get_default_logger()
115 
116  CALL get_qs_env(qs_env=qs_env, &
117  qs_kind_set=qs_kind_set, &
118  particle_set=particle_set, &
119  neighbor_list_id=neighbor_list_id, &
120  para_env=para_env, &
121  sab_orb=sab_orb, &
122  cell=cell)
123 
124  nkind = SIZE(qs_kind_set)
125  natom = SIZE(particle_set)
126 
127  ! *** Allocate work storage ***
128  CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
129  maxco=maxco, &
130  maxlgto=maxlgto, &
131  maxsgf=maxsgf)
132 
133  ldai = ncoset(maxlgto + 2)
134  CALL init_orbital_pointers(ldai)
135 
136  ALLOCATE (rr_work(0:2*maxlgto + 4, ldai, ldai))
137 
138  ALLOCATE (efgab(maxco, maxco, 6))
139 
140  ALLOCATE (work(maxco, maxsgf))
141 
142  ALLOCATE (efgint(6))
143 
144  rr_work(:, :, :) = 0.0_dp
145  efgab(:, :, :) = 0.0_dp
146  work(:, :) = 0.0_dp
147 
148  ALLOCATE (basis_set_list(nkind))
149  DO ikind = 1, nkind
150  qs_kind => qs_kind_set(ikind)
151  CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a)
152  IF (ASSOCIATED(basis_set_a)) THEN
153  basis_set_list(ikind)%gto_basis_set => basis_set_a
154  ELSE
155  NULLIFY (basis_set_list(ikind)%gto_basis_set)
156  END IF
157  END DO
158  CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
159  DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
160  CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
161  iatom=iatom, jatom=jatom, r=rab)
162  basis_set_a => basis_set_list(ikind)%gto_basis_set
163  IF (.NOT. ASSOCIATED(basis_set_a)) cycle
164  basis_set_b => basis_set_list(jkind)%gto_basis_set
165  IF (.NOT. ASSOCIATED(basis_set_b)) cycle
166  ra = pbc(particle_set(iatom)%r, cell)
167  ! basis ikind
168  first_sgfa => basis_set_a%first_sgf
169  la_max => basis_set_a%lmax
170  la_min => basis_set_a%lmin
171  npgfa => basis_set_a%npgf
172  nseta = basis_set_a%nset
173  nsgfa => basis_set_a%nsgf_set
174  rpgfa => basis_set_a%pgf_radius
175  set_radius_a => basis_set_a%set_radius
176  sphi_a => basis_set_a%sphi
177  zeta => basis_set_a%zet
178  ! basis jkind
179  first_sgfb => basis_set_b%first_sgf
180  lb_max => basis_set_b%lmax
181  lb_min => basis_set_b%lmin
182  npgfb => basis_set_b%npgf
183  nsetb = basis_set_b%nset
184  nsgfb => basis_set_b%nsgf_set
185  rpgfb => basis_set_b%pgf_radius
186  set_radius_b => basis_set_b%set_radius
187  sphi_b => basis_set_b%sphi
188  zetb => basis_set_b%zet
189 
190  IF (inode == 1) last_jatom = 0
191 
192  rb = rab + ra
193  rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
194  dab = sqrt(rab2)
195  rac = pbc(ra, rc, cell)
196  rbc = rac - rab
197 
198  IF (jatom /= last_jatom) THEN
199  new_atom_b = .true.
200  last_jatom = jatom
201  ELSE
202  new_atom_b = .false.
203  END IF
204 
205  IF (new_atom_b) THEN
206  IF (iatom <= jatom) THEN
207  irow = iatom
208  icol = jatom
209  ELSE
210  irow = jatom
211  icol = iatom
212  END IF
213 
214  DO i = 1, 6
215  NULLIFY (efgint(i)%block)
216  CALL dbcsr_get_block_p(matrix=matrix_efg(i)%matrix, &
217  row=irow, col=icol, block=efgint(i)%block, found=found)
218  END DO
219  END IF
220 
221  DO iset = 1, nseta
222 
223  ncoa = npgfa(iset)*ncoset(la_max(iset))
224  sgfa = first_sgfa(1, iset)
225 
226  DO jset = 1, nsetb
227 
228  IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
229 
230  ncob = npgfb(jset)*ncoset(lb_max(jset))
231  sgfb = first_sgfb(1, jset)
232 
233  ! *** Calculate the primitive fermi contact integrals ***
234 
235  CALL efg(la_max(iset), la_min(iset), npgfa(iset), &
236  rpgfa(:, iset), zeta(:, iset), &
237  lb_max(jset), lb_min(jset), npgfb(jset), &
238  rpgfb(:, jset), zetb(:, jset), &
239  rac, rbc, rab, efgab, SIZE(rr_work, 1), SIZE(rr_work, 2), rr_work)
240 
241  ! *** Contraction step ***
242 
243  DO i = 1, 6
244 
245  CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
246  1.0_dp, efgab(1, 1, i), SIZE(efgab, 1), &
247  sphi_b(1, sgfb), SIZE(sphi_b, 1), &
248  0.0_dp, work(1, 1), SIZE(work, 1))
249 
250  IF (iatom <= jatom) THEN
251  CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
252  1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
253  work(1, 1), SIZE(work, 1), &
254  1.0_dp, efgint(i)%block(sgfa, sgfb), &
255  SIZE(efgint(i)%block, 1))
256 
257  ELSE
258 
259  CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
260  1.0_dp, work(1, 1), SIZE(work, 1), &
261  sphi_a(1, sgfa), SIZE(sphi_a, 1), &
262  1.0_dp, efgint(i)%block(sgfb, sgfa), &
263  SIZE(efgint(i)%block, 1))
264  END IF
265 
266  END DO
267 
268  END DO
269 
270  END DO
271 
272  END DO
273  CALL neighbor_list_iterator_release(nl_iterator)
274 
275  ! *** Release work storage ***
276  DEALLOCATE (basis_set_list)
277 
278  DEALLOCATE (efgab)
279 
280  DEALLOCATE (work)
281 
282  DEALLOCATE (efgint)
283 
284  ! Print the electric field gradient matrix, if requested
285 
286  IF (btest(cp_print_key_should_output(logger%iter_info, &
287  qs_env%input, "DFT%PRINT%AO_MATRICES/EFG"), cp_p_file)) THEN
288  iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/EFG", &
289  extension=".Log")
290  CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
291  CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers)
292  after = min(max(after, 1), 16)
293  CALL cp_dbcsr_write_sparse_matrix(matrix_efg(1)%matrix, 4, after, qs_env, &
294  para_env, output_unit=iw, omit_headers=omit_headers)
295  CALL cp_dbcsr_write_sparse_matrix(matrix_efg(2)%matrix, 4, after, qs_env, &
296  para_env, output_unit=iw, omit_headers=omit_headers)
297  CALL cp_dbcsr_write_sparse_matrix(matrix_efg(3)%matrix, 4, after, qs_env, &
298  para_env, output_unit=iw, omit_headers=omit_headers)
299  CALL cp_dbcsr_write_sparse_matrix(matrix_efg(4)%matrix, 4, after, qs_env, &
300  para_env, output_unit=iw, omit_headers=omit_headers)
301  CALL cp_dbcsr_write_sparse_matrix(matrix_efg(5)%matrix, 4, after, qs_env, &
302  para_env, output_unit=iw, omit_headers=omit_headers)
303  CALL cp_dbcsr_write_sparse_matrix(matrix_efg(6)%matrix, 4, after, qs_env, &
304  para_env, output_unit=iw, omit_headers=omit_headers)
305  CALL cp_print_key_finished_output(iw, logger, qs_env%input, &
306  "DFT%PRINT%AO_MATRICES/EFG")
307  END IF
308 
309  CALL timestop(handle)
310 
311  END SUBROUTINE build_efg_matrix
312 
313 ! **************************************************************************************************
314 
315 END MODULE qs_elec_field
316 
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Definition: dumpdcd.F:1203
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.
Calculation of Coulomb integrals over Cartesian Gaussian-type functions (electron repulsion integrals...
Definition: ai_elec_field.F:38
subroutine, public efg(la_max, la_min, npgfa, rpgfa, zeta, lb_max, lb_min, npgfb, rpgfb, zetb, rac, rbc, rab, vab, ldrr1, ldrr2, rr)
Calculation of the primitive electric field integrals over Cartesian Gaussian-type functions.
Definition: ai_elec_field.F:83
collect pointers to a block of reals
Definition: block_p_types.F:14
Handles all functions related to the CELL.
Definition: cell_types.F:15
DBCSR output in CP2K.
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)
...
various routines to log and control the output. The idea is that decisions about where to log should ...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
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)
...
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,...
integer, parameter, public cp_p_file
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...
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
Interface to the message passing library MPI.
Provides Cartesian and spherical orbital pointers and indices.
subroutine, public init_orbital_pointers(maxl)
Initialize or update the orbital pointers.
integer, dimension(:), allocatable, public ncoset
Define the data structure for the particle information.
Distribution of the electric field gradient integral matrix.
Definition: qs_elec_field.F:13
subroutine, public build_efg_matrix(qs_env, matrix_efg, rc)
Calculation of the electric field gradient matrix over Cartesian Gaussian functions.
Definition: qs_elec_field.F:75
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.
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
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.
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.
Define the neighbor list data types and the corresponding functionality.
subroutine, public neighbor_list_iterator_create(iterator_set, nl, search, nthread)
Neighbor list iterator functions.
subroutine, public neighbor_list_iterator_release(iterator_set)
...
integer function, public neighbor_list_iterate(iterator_set, mepos)
...
subroutine, public get_iterator_info(iterator_set, mepos, ikind, jkind, nkind, ilist, nlist, inode, nnode, iatom, jatom, r, cell)
...