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 spin orbit integral matrix.
10 : !> \par History
11 : !> \author VW (27.02.2009)
12 : ! **************************************************************************************************
13 : MODULE qs_spin_orbit
14 : USE ai_spin_orbit, ONLY: pso
15 : USE basis_set_types, ONLY: gto_basis_set_p_type,&
16 : gto_basis_set_type
17 : USE block_p_types, ONLY: block_p_type
18 : USE cell_types, ONLY: cell_type,&
19 : pbc
20 : USE cp_dbcsr_api, ONLY: dbcsr_get_block_p,&
21 : dbcsr_p_type
22 : USE cp_dbcsr_output, ONLY: cp_dbcsr_write_sparse_matrix
23 : USE cp_log_handling, ONLY: cp_get_default_logger,&
24 : cp_logger_type
25 : USE cp_output_handling, ONLY: cp_p_file,&
26 : cp_print_key_finished_output,&
27 : cp_print_key_should_output,&
28 : cp_print_key_unit_nr
29 : USE input_section_types, ONLY: section_vals_val_get
30 : USE kinds, ONLY: dp
31 : USE message_passing, ONLY: mp_para_env_type
32 : USE orbital_pointers, ONLY: init_orbital_pointers,&
33 : ncoset
34 : USE particle_types, ONLY: particle_type
35 : USE qs_environment_types, ONLY: get_qs_env,&
36 : qs_environment_type
37 : USE qs_kind_types, ONLY: get_qs_kind,&
38 : get_qs_kind_set,&
39 : qs_kind_type
40 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
41 : neighbor_list_iterate,&
42 : neighbor_list_iterator_create,&
43 : neighbor_list_iterator_p_type,&
44 : neighbor_list_iterator_release,&
45 : neighbor_list_set_p_type
46 : #include "./base/base_uses.f90"
47 :
48 : IMPLICIT NONE
49 :
50 : PRIVATE
51 :
52 : ! *** Global parameters ***
53 :
54 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_spin_orbit'
55 :
56 : ! *** Public subroutines ***
57 :
58 : PUBLIC :: build_pso_matrix
59 :
60 : CONTAINS
61 :
62 : ! **************************************************************************************************
63 : !> \brief Calculation of the paramagnetic spin orbit matrix over
64 : !> Cartesian Gaussian functions.
65 : !> \param qs_env ...
66 : !> \param matrix_so ...
67 : !> \param rc ...
68 : !> \date 27.02.2009
69 : !> \author VW
70 : !> \version 1.0
71 : ! **************************************************************************************************
72 :
73 172 : SUBROUTINE build_pso_matrix(qs_env, matrix_so, rc)
74 :
75 : TYPE(qs_environment_type), POINTER :: qs_env
76 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_so
77 : REAL(dp), DIMENSION(3), INTENT(IN) :: rc
78 :
79 : CHARACTER(len=*), PARAMETER :: routineN = 'build_pso_matrix'
80 :
81 : INTEGER :: after, handle, i, iatom, icol, ikind, inode, irow, iset, iw, jatom, jkind, jset, &
82 : last_jatom, ldai, maxco, maxlgto, maxsgf, natom, ncoa, ncob, neighbor_list_id, nkind, &
83 : nseta, nsetb, sgfa, sgfb
84 172 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
85 172 : npgfb, nsgfa, nsgfb
86 172 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
87 : LOGICAL :: found, new_atom_b, omit_headers
88 : REAL(KIND=dp) :: dab, rab2
89 172 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: work
90 172 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: rr_work, soab
91 : REAL(KIND=dp), DIMENSION(3) :: ra, rab, rac, rb, rbc
92 172 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
93 172 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb
94 172 : TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: soint
95 : TYPE(cell_type), POINTER :: cell
96 : TYPE(cp_logger_type), POINTER :: logger
97 172 : 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 172 : DIMENSION(:), POINTER :: nl_iterator
102 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
103 172 : POINTER :: sab_orb
104 172 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
105 172 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
106 : TYPE(qs_kind_type), POINTER :: qs_kind
107 :
108 172 : CALL timeset(routineN, handle)
109 :
110 172 : NULLIFY (cell, sab_orb, qs_kind_set, particle_set, para_env)
111 172 : NULLIFY (logger)
112 :
113 172 : 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 : neighbor_list_id=neighbor_list_id, &
119 : para_env=para_env, &
120 : sab_orb=sab_orb, &
121 172 : cell=cell)
122 :
123 172 : nkind = SIZE(qs_kind_set)
124 172 : natom = SIZE(particle_set)
125 :
126 : ! *** Allocate work storage ***
127 :
128 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
129 : maxco=maxco, &
130 : maxlgto=maxlgto, &
131 172 : maxsgf=maxsgf)
132 :
133 172 : ldai = ncoset(maxlgto + 1)
134 172 : CALL init_orbital_pointers(ldai)
135 :
136 860 : ALLOCATE (rr_work(0:2*maxlgto + 2, ldai, ldai))
137 860 : ALLOCATE (soab(maxco, maxco, 3))
138 688 : ALLOCATE (work(maxco, maxsgf))
139 688 : ALLOCATE (soint(3))
140 :
141 105092 : rr_work(:, :, :) = 0.0_dp
142 37528 : soab(:, :, :) = 0.0_dp
143 13992 : work(:, :) = 0.0_dp
144 :
145 858 : ALLOCATE (basis_set_list(nkind))
146 514 : DO ikind = 1, nkind
147 342 : qs_kind => qs_kind_set(ikind)
148 342 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a)
149 514 : IF (ASSOCIATED(basis_set_a)) THEN
150 342 : 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 172 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
156 21886 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
157 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
158 21714 : iatom=iatom, jatom=jatom, r=rab)
159 21714 : basis_set_a => basis_set_list(ikind)%gto_basis_set
160 21714 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
161 21714 : basis_set_b => basis_set_list(jkind)%gto_basis_set
162 21714 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
163 21714 : ra = pbc(particle_set(iatom)%r, cell)
164 : ! basis ikind
165 21714 : first_sgfa => basis_set_a%first_sgf
166 21714 : la_max => basis_set_a%lmax
167 21714 : la_min => basis_set_a%lmin
168 21714 : npgfa => basis_set_a%npgf
169 21714 : nseta = basis_set_a%nset
170 21714 : nsgfa => basis_set_a%nsgf_set
171 21714 : rpgfa => basis_set_a%pgf_radius
172 21714 : set_radius_a => basis_set_a%set_radius
173 21714 : sphi_a => basis_set_a%sphi
174 21714 : zeta => basis_set_a%zet
175 : ! basis jkind
176 21714 : first_sgfb => basis_set_b%first_sgf
177 21714 : lb_max => basis_set_b%lmax
178 21714 : lb_min => basis_set_b%lmin
179 21714 : npgfb => basis_set_b%npgf
180 21714 : nsetb = basis_set_b%nset
181 21714 : nsgfb => basis_set_b%nsgf_set
182 21714 : rpgfb => basis_set_b%pgf_radius
183 21714 : set_radius_b => basis_set_b%set_radius
184 21714 : sphi_b => basis_set_b%sphi
185 21714 : zetb => basis_set_b%zet
186 :
187 21714 : IF (inode == 1) last_jatom = 0
188 :
189 86856 : rb = rab + ra
190 21714 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
191 21714 : dab = SQRT(rab2)
192 21714 : rac = pbc(ra, rc, cell)
193 86856 : rbc = rac - rab
194 :
195 21714 : 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 1184 : IF (iatom <= jatom) THEN
204 712 : irow = iatom
205 712 : icol = jatom
206 : ELSE
207 472 : irow = jatom
208 472 : icol = iatom
209 : END IF
210 :
211 4736 : DO i = 1, 3
212 3552 : NULLIFY (soint(i)%block)
213 : CALL dbcsr_get_block_p(matrix=matrix_so(i)%matrix, &
214 4736 : row=irow, col=icol, BLOCK=soint(i)%block, found=found)
215 : END DO
216 : END IF
217 :
218 69988 : DO iset = 1, nseta
219 :
220 48102 : ncoa = npgfa(iset)*ncoset(la_max(iset))
221 48102 : sgfa = first_sgfa(1, iset)
222 :
223 177118 : DO jset = 1, nsetb
224 :
225 107302 : IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
226 :
227 40738 : ncob = npgfb(jset)*ncoset(lb_max(jset))
228 40738 : sgfb = first_sgfb(1, jset)
229 :
230 : ! *** Calculate the primitive fermi contact integrals ***
231 :
232 : CALL pso(la_max(iset), la_min(iset), npgfa(iset), &
233 : rpgfa(:, iset), zeta(:, iset), &
234 : lb_max(jset), lb_min(jset), npgfb(jset), &
235 : rpgfb(:, jset), zetb(:, jset), &
236 40738 : rac, rbc, rab, soab, SIZE(rr_work, 1), SIZE(rr_work, 2), rr_work)
237 :
238 : ! *** Contraction step ***
239 :
240 211054 : DO i = 1, 3
241 :
242 : CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
243 : 1.0_dp, soab(1, 1, i), SIZE(soab, 1), &
244 : sphi_b(1, sgfb), SIZE(sphi_b, 1), &
245 122214 : 0.0_dp, work(1, 1), SIZE(work, 1))
246 :
247 229516 : IF (iatom <= jatom) THEN
248 :
249 : CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
250 : 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
251 : work(1, 1), SIZE(work, 1), &
252 : 1.0_dp, soint(i)%block(sgfa, sgfb), &
253 73815 : SIZE(soint(i)%block, 1))
254 :
255 : ELSE
256 :
257 : CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
258 : -1.0_dp, work(1, 1), SIZE(work, 1), &
259 : sphi_a(1, sgfa), SIZE(sphi_a, 1), &
260 : 1.0_dp, soint(i)%block(sgfb, sgfa), &
261 48399 : SIZE(soint(i)%block, 1))
262 : END IF
263 :
264 : END DO
265 :
266 : END DO
267 :
268 : END DO
269 :
270 : END DO
271 172 : CALL neighbor_list_iterator_release(nl_iterator)
272 :
273 : ! *** Release work storage ***
274 :
275 172 : DEALLOCATE (basis_set_list)
276 :
277 172 : DEALLOCATE (soab)
278 :
279 172 : DEALLOCATE (work)
280 :
281 172 : DEALLOCATE (soint)
282 :
283 : ! *** Print the spin orbit matrix, if requested ***
284 :
285 172 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
286 : qs_env%input, "DFT%PRINT%AO_MATRICES/PSO"), cp_p_file)) THEN
287 : iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/PSO", &
288 8 : extension=".Log")
289 8 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
290 8 : after = MIN(MAX(after, 1), 16)
291 8 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers)
292 : CALL cp_dbcsr_write_sparse_matrix(matrix_so(1)%matrix, 4, after, qs_env, &
293 8 : para_env, output_unit=iw, omit_headers=omit_headers)
294 : CALL cp_dbcsr_write_sparse_matrix(matrix_so(2)%matrix, 4, after, qs_env, &
295 8 : para_env, output_unit=iw, omit_headers=omit_headers)
296 : CALL cp_dbcsr_write_sparse_matrix(matrix_so(3)%matrix, 4, after, qs_env, &
297 8 : para_env, output_unit=iw, omit_headers=omit_headers)
298 : CALL cp_print_key_finished_output(iw, logger, qs_env%input, &
299 8 : "DFT%PRINT%AO_MATRICES/PSO")
300 : END IF
301 :
302 172 : CALL timestop(handle)
303 :
304 516 : END SUBROUTINE build_pso_matrix
305 :
306 : ! **************************************************************************************************
307 :
308 : END MODULE qs_spin_orbit
309 :
|
