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 Distribution of the electric field gradient integral matrix.
10 : !> \par History
11 : !> \author VW (27.02.2009)
12 : ! **************************************************************************************************
13 : MODULE qs_elec_field
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
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_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 44 : 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 44 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
86 44 : npgfb, nsgfa, nsgfb
87 44 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
88 : LOGICAL :: found, new_atom_b, omit_headers
89 : REAL(KIND=dp) :: dab, rab2
90 44 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: work
91 44 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: efgab, rr_work
92 : REAL(KIND=dp), DIMENSION(3) :: ra, rab, rac, rb, rbc
93 44 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
94 44 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb
95 44 : TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: efgint
96 : TYPE(cell_type), POINTER :: cell
97 : TYPE(cp_logger_type), POINTER :: logger
98 44 : 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 44 : DIMENSION(:), POINTER :: nl_iterator
103 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
104 44 : POINTER :: sab_orb
105 44 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
106 44 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
107 : TYPE(qs_kind_type), POINTER :: qs_kind
108 :
109 44 : CALL timeset(routineN, handle)
110 :
111 44 : NULLIFY (cell, sab_orb, qs_kind_set, particle_set, para_env)
112 44 : NULLIFY (logger)
113 :
114 44 : 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 44 : cell=cell)
123 :
124 44 : nkind = SIZE(qs_kind_set)
125 44 : 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 44 : maxsgf=maxsgf)
132 :
133 44 : ldai = ncoset(maxlgto + 2)
134 44 : CALL init_orbital_pointers(ldai)
135 :
136 220 : ALLOCATE (rr_work(0:2*maxlgto + 4, ldai, ldai))
137 :
138 220 : ALLOCATE (efgab(maxco, maxco, 6))
139 :
140 176 : ALLOCATE (work(maxco, maxsgf))
141 :
142 308 : ALLOCATE (efgint(6))
143 :
144 141724 : rr_work(:, :, :) = 0.0_dp
145 18692 : efgab(:, :, :) = 0.0_dp
146 3496 : work(:, :) = 0.0_dp
147 :
148 218 : ALLOCATE (basis_set_list(nkind))
149 130 : DO ikind = 1, nkind
150 86 : qs_kind => qs_kind_set(ikind)
151 86 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a)
152 130 : IF (ASSOCIATED(basis_set_a)) THEN
153 86 : basis_set_list(ikind)%gto_basis_set => basis_set_a
154 : ELSE
155 0 : NULLIFY (basis_set_list(ikind)%gto_basis_set)
156 : END IF
157 : END DO
158 44 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
159 4918 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
160 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
161 4874 : iatom=iatom, jatom=jatom, r=rab)
162 4874 : basis_set_a => basis_set_list(ikind)%gto_basis_set
163 4874 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
164 4874 : basis_set_b => basis_set_list(jkind)%gto_basis_set
165 4874 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
166 4874 : ra = pbc(particle_set(iatom)%r, cell)
167 : ! basis ikind
168 4874 : first_sgfa => basis_set_a%first_sgf
169 4874 : la_max => basis_set_a%lmax
170 4874 : la_min => basis_set_a%lmin
171 4874 : npgfa => basis_set_a%npgf
172 4874 : nseta = basis_set_a%nset
173 4874 : nsgfa => basis_set_a%nsgf_set
174 4874 : rpgfa => basis_set_a%pgf_radius
175 4874 : set_radius_a => basis_set_a%set_radius
176 4874 : sphi_a => basis_set_a%sphi
177 4874 : zeta => basis_set_a%zet
178 : ! basis jkind
179 4874 : first_sgfb => basis_set_b%first_sgf
180 4874 : lb_max => basis_set_b%lmax
181 4874 : lb_min => basis_set_b%lmin
182 4874 : npgfb => basis_set_b%npgf
183 4874 : nsetb = basis_set_b%nset
184 4874 : nsgfb => basis_set_b%nsgf_set
185 4874 : rpgfb => basis_set_b%pgf_radius
186 4874 : set_radius_b => basis_set_b%set_radius
187 4874 : sphi_b => basis_set_b%sphi
188 4874 : zetb => basis_set_b%zet
189 :
190 4874 : IF (inode == 1) last_jatom = 0
191 :
192 19496 : rb = rab + ra
193 4874 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
194 4874 : dab = SQRT(rab2)
195 4874 : rac = pbc(ra, rc, cell)
196 19496 : rbc = rac - rab
197 :
198 4874 : 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 272 : IF (iatom <= jatom) THEN
207 164 : irow = iatom
208 164 : icol = jatom
209 : ELSE
210 108 : irow = jatom
211 108 : icol = iatom
212 : END IF
213 :
214 1904 : DO i = 1, 6
215 1632 : NULLIFY (efgint(i)%block)
216 : CALL dbcsr_get_block_p(matrix=matrix_efg(i)%matrix, &
217 1904 : row=irow, col=icol, BLOCK=efgint(i)%block, found=found)
218 : END DO
219 : END IF
220 :
221 15700 : DO iset = 1, nseta
222 :
223 10782 : ncoa = npgfa(iset)*ncoset(la_max(iset))
224 10782 : sgfa = first_sgfa(1, iset)
225 :
226 39798 : DO jset = 1, nsetb
227 :
228 24142 : IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
229 :
230 9222 : ncob = npgfb(jset)*ncoset(lb_max(jset))
231 9222 : 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 9222 : rac, rbc, rab, efgab, SIZE(rr_work, 1), SIZE(rr_work, 2), rr_work)
240 :
241 : ! *** Contraction step ***
242 :
243 75336 : 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 55332 : 0.0_dp, work(1, 1), SIZE(work, 1))
249 :
250 79474 : 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 33582 : 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 21750 : 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 44 : CALL neighbor_list_iterator_release(nl_iterator)
274 :
275 : ! *** Release work storage ***
276 44 : DEALLOCATE (basis_set_list)
277 :
278 44 : DEALLOCATE (efgab)
279 :
280 44 : DEALLOCATE (work)
281 :
282 44 : DEALLOCATE (efgint)
283 :
284 : ! Print the electric field gradient matrix, if requested
285 :
286 44 : 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 8 : extension=".Log")
290 8 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
291 8 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers)
292 8 : after = MIN(MAX(after, 1), 16)
293 : CALL cp_dbcsr_write_sparse_matrix(matrix_efg(1)%matrix, 4, after, qs_env, &
294 8 : 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 8 : 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 8 : 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 8 : 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 8 : 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 8 : para_env, output_unit=iw, omit_headers=omit_headers)
305 : CALL cp_print_key_finished_output(iw, logger, qs_env%input, &
306 8 : "DFT%PRINT%AO_MATRICES/EFG")
307 : END IF
308 :
309 44 : CALL timestop(handle)
310 :
311 132 : END SUBROUTINE build_efg_matrix
312 :
313 : ! **************************************************************************************************
314 :
315 : END MODULE qs_elec_field
316 :
|
