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 Fermi contact integral matrix.
10 : !> \par History
11 : !> \author VW (27.02.2009)
12 : ! **************************************************************************************************
13 : MODULE 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 :
61 : CONTAINS
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 44 : 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 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(:, :) :: fcab, work
91 : REAL(KIND=dp), DIMENSION(3) :: ra, rab, rac, rb, rbc
92 44 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
93 44 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb
94 44 : TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: fcint
95 : TYPE(cell_type), POINTER :: cell
96 : TYPE(cp_logger_type), POINTER :: logger
97 44 : 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 44 : DIMENSION(:), POINTER :: nl_iterator
102 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
103 44 : POINTER :: sab_orb
104 44 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
105 44 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
106 : TYPE(qs_kind_type), POINTER :: qs_kind
107 :
108 44 : CALL timeset(routineN, handle)
109 :
110 44 : NULLIFY (cell, sab_orb, qs_kind_set, particle_set, para_env)
111 44 : NULLIFY (logger)
112 :
113 44 : 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 44 : cell=cell)
121 :
122 44 : nkind = SIZE(qs_kind_set)
123 44 : 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 44 : maxsgf=maxsgf)
131 :
132 44 : ldai = ncoset(maxlgto)
133 44 : CALL init_orbital_pointers(ldai)
134 :
135 44 : ldfc = maxco
136 176 : ALLOCATE (fcab(ldfc, ldfc))
137 3108 : fcab(:, :) = 0.0_dp
138 :
139 176 : ALLOCATE (work(maxco, maxsgf))
140 3496 : work(:, :) = 0.0_dp
141 :
142 88 : ALLOCATE (fcint(1))
143 44 : NULLIFY (fcint(1)%block)
144 :
145 218 : ALLOCATE (basis_set_list(nkind))
146 130 : DO ikind = 1, nkind
147 86 : qs_kind => qs_kind_set(ikind)
148 86 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a)
149 130 : IF (ASSOCIATED(basis_set_a)) THEN
150 86 : basis_set_list(ikind)%gto_basis_set => basis_set_a
151 : ELSE
152 0 : NULLIFY (basis_set_list(ikind)%gto_basis_set)
153 : END IF
154 : END DO
155 44 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
156 4918 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
157 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
158 4874 : iatom=iatom, jatom=jatom, r=rab)
159 4874 : basis_set_a => basis_set_list(ikind)%gto_basis_set
160 4874 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
161 4874 : basis_set_b => basis_set_list(jkind)%gto_basis_set
162 4874 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
163 4874 : ra = pbc(particle_set(iatom)%r, cell)
164 : ! basis ikind
165 4874 : first_sgfa => basis_set_a%first_sgf
166 4874 : la_max => basis_set_a%lmax
167 4874 : la_min => basis_set_a%lmin
168 4874 : npgfa => basis_set_a%npgf
169 4874 : nseta = basis_set_a%nset
170 4874 : nsgfa => basis_set_a%nsgf_set
171 4874 : rpgfa => basis_set_a%pgf_radius
172 4874 : set_radius_a => basis_set_a%set_radius
173 4874 : sphi_a => basis_set_a%sphi
174 4874 : zeta => basis_set_a%zet
175 : ! basis jkind
176 4874 : first_sgfb => basis_set_b%first_sgf
177 4874 : lb_max => basis_set_b%lmax
178 4874 : lb_min => basis_set_b%lmin
179 4874 : npgfb => basis_set_b%npgf
180 4874 : nsetb = basis_set_b%nset
181 4874 : nsgfb => basis_set_b%nsgf_set
182 4874 : rpgfb => basis_set_b%pgf_radius
183 4874 : set_radius_b => basis_set_b%set_radius
184 4874 : sphi_b => basis_set_b%sphi
185 4874 : zetb => basis_set_b%zet
186 :
187 4874 : IF (inode == 1) last_jatom = 0
188 :
189 19496 : rb = rab + ra
190 4874 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
191 4874 : dab = SQRT(rab2)
192 4874 : rac = pbc(ra, rc, cell)
193 19496 : rbc = rac - rab
194 :
195 4874 : 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 272 : IF (new_atom_b) THEN
203 272 : IF (iatom <= jatom) THEN
204 164 : irow = iatom
205 164 : icol = jatom
206 : ELSE
207 108 : irow = jatom
208 108 : icol = iatom
209 : END IF
210 :
211 272 : NULLIFY (fcint(1)%block)
212 : CALL dbcsr_get_block_p(matrix=matrix_fc(1)%matrix, &
213 272 : row=irow, col=icol, BLOCK=fcint(1)%block, found=found)
214 : END IF
215 :
216 15700 : DO iset = 1, nseta
217 :
218 10782 : ncoa = npgfa(iset)*ncoset(la_max(iset))
219 10782 : sgfa = first_sgfa(1, iset)
220 :
221 39798 : DO jset = 1, nsetb
222 :
223 24142 : IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
224 :
225 9222 : ncob = npgfb(jset)*ncoset(lb_max(jset))
226 9222 : 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 9222 : 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 9222 : 0.0_dp, work(1, 1), SIZE(work, 1))
242 :
243 20004 : 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 5597 : 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 3625 : SIZE(fcint(1)%block, 1))
258 :
259 : END IF
260 :
261 : END DO
262 :
263 : END DO
264 :
265 : END DO
266 44 : CALL neighbor_list_iterator_release(nl_iterator)
267 :
268 : ! *** Release work storage ***
269 44 : DEALLOCATE (basis_set_list)
270 :
271 44 : DEALLOCATE (fcab)
272 :
273 44 : DEALLOCATE (work)
274 :
275 44 : NULLIFY (fcint(1)%block)
276 44 : DEALLOCATE (fcint)
277 :
278 : ! *** Print the Fermi contact matrix, if requested ***
279 :
280 44 : 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 8 : extension=".Log")
284 8 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
285 8 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers)
286 8 : after = MIN(MAX(after, 1), 16)
287 : CALL cp_dbcsr_write_sparse_matrix(matrix_fc(1)%matrix, 4, after, qs_env, &
288 8 : para_env, output_unit=iw, omit_headers=omit_headers)
289 : CALL cp_print_key_finished_output(iw, logger, qs_env%input, &
290 8 : "DFT%PRINT%AO_MATRICES/FERMI_CONTACT")
291 : END IF
292 :
293 44 : CALL timestop(handle)
294 :
295 88 : END SUBROUTINE build_fermi_contact_matrix
296 :
297 : ! **************************************************************************************************
298 :
299 : END MODULE qs_fermi_contact
300 :
|
