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 : !> \par History
10 : !> created 06-2006 [RD]
11 : !> \author RD
12 : ! **************************************************************************************************
13 : MODULE qs_linres_epr_ownutils
14 : USE atomic_kind_types, ONLY: atomic_kind_type,&
15 : get_atomic_kind
16 : USE cell_types, ONLY: cell_type
17 : USE cp_control_types, ONLY: dft_control_type
18 : USE cp_dbcsr_api, ONLY: dbcsr_p_type
19 : USE cp_log_handling, ONLY: cp_get_default_logger,&
20 : cp_logger_type
21 : USE cp_output_handling, ONLY: cp_p_file,&
22 : cp_print_key_finished_output,&
23 : cp_print_key_should_output,&
24 : cp_print_key_unit_nr
25 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
26 : section_vals_type,&
27 : section_vals_val_get
28 : USE kinds, ONLY: default_string_length,&
29 : dp
30 : USE mathlib, ONLY: diamat_all
31 : USE message_passing, ONLY: mp_para_env_type
32 : USE particle_types, ONLY: particle_type
33 : USE pw_env_types, ONLY: pw_env_get,&
34 : pw_env_type
35 : USE pw_methods, ONLY: pw_axpy,&
36 : pw_integral_ab,&
37 : pw_scale,&
38 : pw_transfer,&
39 : pw_zero
40 : USE pw_pool_types, ONLY: pw_pool_p_type,&
41 : pw_pool_type
42 : USE pw_spline_utils, ONLY: Eval_Interp_Spl3_pbc,&
43 : find_coeffs,&
44 : pw_spline_do_precond,&
45 : pw_spline_precond_create,&
46 : pw_spline_precond_release,&
47 : pw_spline_precond_set_kind,&
48 : pw_spline_precond_type,&
49 : spl3_pbc
50 : USE pw_types, ONLY: pw_c1d_gs_type,&
51 : pw_r3d_rs_type
52 : USE qs_core_energies, ONLY: calculate_ptrace
53 : USE qs_environment_types, ONLY: get_qs_env,&
54 : qs_environment_type
55 : USE qs_grid_atom, ONLY: grid_atom_type
56 : USE qs_harmonics_atom, ONLY: harmonics_atom_type
57 : USE qs_kind_types, ONLY: get_qs_kind,&
58 : qs_kind_type
59 : USE qs_linres_nmr_epr_common_utils, ONLY: mult_G_ov_G2_grid
60 : USE qs_linres_op, ONLY: fac_vecp,&
61 : set_vecp,&
62 : set_vecp_rev
63 : USE qs_linres_types, ONLY: current_env_type,&
64 : epr_env_type,&
65 : get_current_env,&
66 : get_epr_env,&
67 : jrho_atom_type,&
68 : nablavks_atom_type
69 : USE qs_rho_atom_types, ONLY: get_rho_atom,&
70 : rho_atom_coeff,&
71 : rho_atom_type
72 : USE qs_rho_types, ONLY: qs_rho_get,&
73 : qs_rho_p_type,&
74 : qs_rho_type
75 : USE realspace_grid_types, ONLY: realspace_grid_desc_type
76 : USE util, ONLY: get_limit
77 : #include "./base/base_uses.f90"
78 :
79 : IMPLICIT NONE
80 :
81 : PRIVATE
82 : PUBLIC :: epr_g_print, epr_g_zke, epr_g_so, epr_g_soo, epr_ind_magnetic_field
83 :
84 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_linres_epr_ownutils'
85 :
86 : CONTAINS
87 :
88 : ! **************************************************************************************************
89 : !> \brief Prints the g tensor
90 : !> \param epr_env ...
91 : !> \param qs_env ...
92 : !> \par History
93 : !> 06.2006 created [RD]
94 : !> \author RD
95 : ! **************************************************************************************************
96 14 : SUBROUTINE epr_g_print(epr_env, qs_env)
97 :
98 : TYPE(epr_env_type) :: epr_env
99 : TYPE(qs_environment_type), POINTER :: qs_env
100 :
101 : CHARACTER(LEN=default_string_length) :: title
102 : INTEGER :: idir1, idir2, output_unit, unit_nr
103 : REAL(KIND=dp) :: eigenv_g(3), g_sym(3, 3), gsum
104 : TYPE(cp_logger_type), POINTER :: logger
105 : TYPE(section_vals_type), POINTER :: lr_section
106 :
107 14 : NULLIFY (logger, lr_section)
108 :
109 14 : logger => cp_get_default_logger()
110 14 : lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES")
111 :
112 : output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
113 14 : extension=".linresLog")
114 :
115 14 : gsum = 0.0_dp
116 :
117 56 : DO idir1 = 1, 3
118 182 : DO idir2 = 1, 3
119 168 : gsum = gsum + epr_env%g_total(idir1, idir2)
120 : END DO
121 : END DO
122 :
123 14 : IF (output_unit > 0) THEN
124 : WRITE (UNIT=output_unit, FMT="(T2,A,T56,E14.6)") &
125 7 : "epr|TOT:checksum", gsum
126 : END IF
127 :
128 : CALL cp_print_key_finished_output(output_unit, logger, lr_section, &
129 14 : "PRINT%PROGRAM_RUN_INFO")
130 :
131 14 : IF (BTEST(cp_print_key_should_output(logger%iter_info, lr_section, &
132 : "EPR%PRINT%G_TENSOR"), cp_p_file)) THEN
133 :
134 : unit_nr = cp_print_key_unit_nr(logger, lr_section, "EPR%PRINT%G_TENSOR", &
135 : extension=".data", middle_name="GTENSOR", &
136 14 : log_filename=.FALSE.)
137 :
138 14 : IF (unit_nr > 0) THEN
139 :
140 7 : WRITE (title, "(A)") "G tensor "
141 7 : WRITE (unit_nr, "(T2,A)") title
142 :
143 7 : WRITE (unit_nr, "(T2,A)") "gmatrix_zke"
144 7 : WRITE (unit_nr, "(3(A,f15.10))") " XX=", epr_env%g_zke, &
145 14 : " XY=", 0.0_dp, " XZ=", 0.0_dp
146 7 : WRITE (unit_nr, "(3(A,f15.10))") " YX=", 0.0_dp, &
147 14 : " YY=", epr_env%g_zke, " YZ=", 0.0_dp
148 7 : WRITE (unit_nr, "(3(A,f15.10))") " ZX=", 0.0_dp, &
149 14 : " ZY=", 0.0_dp, " ZZ=", epr_env%g_zke
150 :
151 7 : WRITE (unit_nr, "(T2,A)") "gmatrix_so"
152 7 : WRITE (unit_nr, "(3(A,f15.10))") " XX=", epr_env%g_so(1, 1), &
153 14 : " XY=", epr_env%g_so(1, 2), " XZ=", epr_env%g_so(1, 3)
154 7 : WRITE (unit_nr, "(3(A,f15.10))") " YX=", epr_env%g_so(2, 1), &
155 14 : " YY=", epr_env%g_so(2, 2), " YZ=", epr_env%g_so(2, 3)
156 7 : WRITE (unit_nr, "(3(A,f15.10))") " ZX=", epr_env%g_so(3, 1), &
157 14 : " ZY=", epr_env%g_so(3, 2), " ZZ=", epr_env%g_so(3, 3)
158 :
159 7 : WRITE (unit_nr, "(T2,A)") "gmatrix_soo"
160 7 : WRITE (unit_nr, "(3(A,f15.10))") " XX=", epr_env%g_soo(1, 1), &
161 14 : " XY=", epr_env%g_soo(1, 2), " XZ=", epr_env%g_soo(1, 3)
162 7 : WRITE (unit_nr, "(3(A,f15.10))") " YX=", epr_env%g_soo(2, 1), &
163 14 : " YY=", epr_env%g_soo(2, 2), " YZ=", epr_env%g_soo(2, 3)
164 7 : WRITE (unit_nr, "(3(A,f15.10))") " ZX=", epr_env%g_soo(3, 1), &
165 14 : " ZY=", epr_env%g_soo(3, 2), " ZZ=", epr_env%g_soo(3, 3)
166 :
167 7 : WRITE (unit_nr, "(T2,A)") "gmatrix_total"
168 7 : WRITE (unit_nr, "(3(A,f15.10))") " XX=", epr_env%g_total(1, 1) + epr_env%g_free_factor, &
169 14 : " XY=", epr_env%g_total(1, 2), " XZ=", epr_env%g_total(1, 3)
170 7 : WRITE (unit_nr, "(3(A,f15.10))") " YX=", epr_env%g_total(2, 1), &
171 14 : " YY=", epr_env%g_total(2, 2) + epr_env%g_free_factor, " YZ=", epr_env%g_total(2, 3)
172 7 : WRITE (unit_nr, "(3(A,f15.10))") " ZX=", epr_env%g_total(3, 1), &
173 14 : " ZY=", epr_env%g_total(3, 2), " ZZ=", epr_env%g_total(3, 3) + epr_env%g_free_factor
174 :
175 28 : DO idir1 = 1, 3
176 91 : DO idir2 = 1, 3
177 : g_sym(idir1, idir2) = (epr_env%g_total(idir1, idir2) + &
178 84 : epr_env%g_total(idir2, idir1))/2.0_dp
179 : END DO
180 : END DO
181 :
182 7 : WRITE (unit_nr, "(T2,A)") "gtensor_total"
183 7 : WRITE (unit_nr, "(3(A,f15.10))") " XX=", g_sym(1, 1) + epr_env%g_free_factor, &
184 14 : " XY=", g_sym(1, 2), " XZ=", g_sym(1, 3)
185 7 : WRITE (unit_nr, "(3(A,f15.10))") " YX=", g_sym(2, 1), &
186 14 : " YY=", g_sym(2, 2) + epr_env%g_free_factor, " YZ=", g_sym(2, 3)
187 7 : WRITE (unit_nr, "(3(A,f15.10))") " ZX=", g_sym(3, 1), &
188 14 : " ZY=", g_sym(3, 2), " ZZ=", g_sym(3, 3) + epr_env%g_free_factor
189 :
190 7 : CALL diamat_all(g_sym, eigenv_g)
191 28 : eigenv_g(:) = eigenv_g(:)*1.0e6_dp
192 :
193 7 : WRITE (unit_nr, "(T2,A)") "delta_g principal values in ppm"
194 7 : WRITE (unit_nr, "(f15.3,3(A,f15.10))") eigenv_g(1), " X=", g_sym(1, 1), &
195 14 : " Y=", g_sym(2, 1), " Z=", g_sym(3, 1)
196 7 : WRITE (unit_nr, "(f15.3,3(A,f15.10))") eigenv_g(2), " X=", g_sym(1, 2), &
197 14 : " Y=", g_sym(2, 2), " Z=", g_sym(3, 2)
198 7 : WRITE (unit_nr, "(f15.3,3(A,f15.10))") eigenv_g(3), " X=", g_sym(1, 3), &
199 14 : " Y=", g_sym(2, 3), " Z=", g_sym(3, 3)
200 :
201 : END IF
202 :
203 : CALL cp_print_key_finished_output(unit_nr, logger, lr_section,&
204 14 : & "EPR%PRINT%G_TENSOR")
205 :
206 : END IF
207 :
208 14 : END SUBROUTINE epr_g_print
209 :
210 : ! **************************************************************************************************
211 : !> \brief Calculate zke part of the g tensor
212 : !> \param epr_env ...
213 : !> \param qs_env ...
214 : !> \par History
215 : !> 06.2006 created [RD]
216 : !> \author RD
217 : ! **************************************************************************************************
218 14 : SUBROUTINE epr_g_zke(epr_env, qs_env)
219 :
220 : TYPE(epr_env_type) :: epr_env
221 : TYPE(qs_environment_type), POINTER :: qs_env
222 :
223 : INTEGER :: i1, ispin, output_unit
224 : REAL(KIND=dp) :: epr_g_zke_temp(2)
225 : TYPE(cp_logger_type), POINTER :: logger
226 14 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: kinetic, rho_ao
227 : TYPE(dft_control_type), POINTER :: dft_control
228 : TYPE(mp_para_env_type), POINTER :: para_env
229 : TYPE(qs_rho_type), POINTER :: rho
230 : TYPE(section_vals_type), POINTER :: lr_section
231 :
232 14 : NULLIFY (dft_control, logger, lr_section, rho, kinetic, para_env, rho_ao)
233 :
234 28 : logger => cp_get_default_logger()
235 14 : lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES")
236 :
237 : output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
238 14 : extension=".linresLog")
239 :
240 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, &
241 14 : kinetic=kinetic, rho=rho, para_env=para_env)
242 :
243 14 : CALL qs_rho_get(rho, rho_ao=rho_ao)
244 :
245 42 : DO ispin = 1, dft_control%nspins
246 : CALL calculate_ptrace(kinetic(1)%matrix, rho_ao(ispin)%matrix, &
247 42 : ecore=epr_g_zke_temp(ispin))
248 : END DO
249 :
250 14 : epr_env%g_zke = epr_env%g_zke_factor*(epr_g_zke_temp(1) - epr_g_zke_temp(2))
251 56 : DO i1 = 1, 3
252 56 : epr_env%g_total(i1, i1) = epr_env%g_total(i1, i1) + epr_env%g_zke
253 : END DO
254 :
255 14 : IF (output_unit > 0) THEN
256 : WRITE (UNIT=output_unit, FMT="(T2,A,T56,E24.16)") &
257 7 : "epr|ZKE:g_zke", epr_env%g_zke
258 : END IF
259 :
260 : CALL cp_print_key_finished_output(output_unit, logger, lr_section, &
261 14 : "PRINT%PROGRAM_RUN_INFO")
262 :
263 14 : END SUBROUTINE epr_g_zke
264 :
265 : ! **************************************************************************************************
266 : !> \brief Calculates g_so
267 : !> \param epr_env ...
268 : !> \param current_env ...
269 : !> \param qs_env ...
270 : !> \param iB ...
271 : !> \par History
272 : !> 06.2006 created [RD]
273 : !> \author RD
274 : ! **************************************************************************************************
275 42 : SUBROUTINE epr_g_so(epr_env, current_env, qs_env, iB)
276 :
277 : TYPE(epr_env_type) :: epr_env
278 : TYPE(current_env_type) :: current_env
279 : TYPE(qs_environment_type), POINTER :: qs_env
280 : INTEGER, INTENT(IN) :: iB
281 :
282 : INTEGER :: aint_precond, ia, iat, iatom, idir1, &
283 : idir2, idir3, ikind, ir, ispin, &
284 : max_iter, natom, nkind, nspins, &
285 : output_unit, precond_kind
286 : INTEGER, DIMENSION(2) :: bo
287 42 : INTEGER, DIMENSION(:), POINTER :: atom_list
288 : LOGICAL :: gapw, paw_atom, success
289 : REAL(dp) :: eps_r, eps_x, hard_radius, temp_so_soft, &
290 : vks_ra_idir2, vks_ra_idir3
291 : REAL(dp), DIMENSION(3, 3) :: temp_so_gapw
292 42 : REAL(dp), DIMENSION(:, :), POINTER :: g_so, g_total
293 : REAL(KIND=dp), DIMENSION(3) :: ra
294 42 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
295 : TYPE(cp_logger_type), POINTER :: logger
296 : TYPE(dft_control_type), POINTER :: dft_control
297 : TYPE(grid_atom_type), POINTER :: grid_atom
298 : TYPE(harmonics_atom_type), POINTER :: harmonics
299 42 : TYPE(jrho_atom_type), DIMENSION(:), POINTER :: jrho1_atom_set
300 : TYPE(jrho_atom_type), POINTER :: jrho1_atom
301 : TYPE(mp_para_env_type), POINTER :: para_env
302 42 : TYPE(nablavks_atom_type), DIMENSION(:), POINTER :: nablavks_atom_set
303 : TYPE(nablavks_atom_type), POINTER :: nablavks_atom
304 42 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
305 : TYPE(pw_env_type), POINTER :: pw_env
306 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
307 42 : TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:, :) :: vks_pw_spline
308 42 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: jrho2_r, jrho3_r, nrho1_r, nrho2_r, &
309 42 : nrho3_r
310 : TYPE(pw_spline_precond_type) :: precond
311 42 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
312 42 : TYPE(qs_rho_p_type), DIMENSION(:), POINTER :: jrho1_set
313 42 : TYPE(qs_rho_p_type), DIMENSION(:, :), POINTER :: nablavks_set
314 : TYPE(section_vals_type), POINTER :: interp_section, lr_section
315 :
316 42 : NULLIFY (atomic_kind_set, qs_kind_set, atom_list, dft_control, &
317 42 : grid_atom, g_so, g_total, harmonics, interp_section, jrho1_atom_set, &
318 42 : jrho1_set, logger, lr_section, nablavks_atom, nablavks_atom_set, &
319 42 : nablavks_set, para_env, particle_set, jrho2_r, jrho3_r, nrho1_r, nrho2_r, nrho3_r)
320 :
321 84 : logger => cp_get_default_logger()
322 42 : lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES")
323 :
324 : output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
325 42 : extension=".linresLog")
326 :
327 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, &
328 : atomic_kind_set=atomic_kind_set, &
329 : qs_kind_set=qs_kind_set, &
330 : para_env=para_env, pw_env=pw_env, &
331 42 : particle_set=particle_set)
332 :
333 : CALL get_epr_env(epr_env=epr_env, &
334 : nablavks_set=nablavks_set, &
335 : nablavks_atom_set=nablavks_atom_set, &
336 42 : g_total=g_total, g_so=g_so)
337 :
338 : CALL get_current_env(current_env=current_env, &
339 42 : jrho1_set=jrho1_set, jrho1_atom_set=jrho1_atom_set)
340 :
341 42 : gapw = dft_control%qs_control%gapw
342 42 : nkind = SIZE(qs_kind_set, 1)
343 42 : nspins = dft_control%nspins
344 :
345 168 : DO idir1 = 1, 3
346 126 : CALL set_vecp(idir1, idir2, idir3)
347 : ! j_pw x nabla_vks_pw
348 126 : temp_so_soft = 0.0_dp
349 378 : DO ispin = 1, nspins
350 252 : CALL qs_rho_get(jrho1_set(idir2)%rho, rho_r=jrho2_r)
351 252 : CALL qs_rho_get(jrho1_set(idir3)%rho, rho_r=jrho3_r)
352 252 : CALL qs_rho_get(nablavks_set(idir2, ispin)%rho, rho_r=nrho2_r)
353 252 : CALL qs_rho_get(nablavks_set(idir3, ispin)%rho, rho_r=nrho3_r)
354 : temp_so_soft = temp_so_soft + (-1.0_dp)**(1 + ispin)*( &
355 : pw_integral_ab(jrho2_r(ispin), nrho3_r(1)) - &
356 378 : pw_integral_ab(jrho3_r(ispin), nrho2_r(1)))
357 : END DO
358 126 : temp_so_soft = -1.0_dp*epr_env%g_so_factor*temp_so_soft
359 126 : IF (output_unit > 0) THEN
360 : WRITE (UNIT=output_unit, FMT="(T2,A,T18,I1,I1,T56,E24.16)") &
361 63 : "epr|SOX:soft", iB, idir1, temp_so_soft
362 : END IF
363 294 : g_so(iB, idir1) = temp_so_soft
364 : END DO !idir1
365 :
366 42 : IF (gapw) THEN
367 :
368 30 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
369 330 : ALLOCATE (vks_pw_spline(3, nspins))
370 :
371 : interp_section => section_vals_get_subs_vals(lr_section, &
372 30 : "EPR%INTERPOLATOR")
373 : CALL section_vals_val_get(interp_section, "aint_precond", &
374 30 : i_val=aint_precond)
375 30 : CALL section_vals_val_get(interp_section, "precond", i_val=precond_kind)
376 30 : CALL section_vals_val_get(interp_section, "max_iter", i_val=max_iter)
377 30 : CALL section_vals_val_get(interp_section, "eps_r", r_val=eps_r)
378 30 : CALL section_vals_val_get(interp_section, "eps_x", r_val=eps_x)
379 :
380 90 : DO ispin = 1, nspins
381 270 : DO idir1 = 1, 3
382 180 : CALL auxbas_pw_pool%create_pw(vks_pw_spline(idir1, ispin))
383 : ! calculate spline coefficients
384 : CALL pw_spline_precond_create(precond, precond_kind=aint_precond, &
385 180 : pool=auxbas_pw_pool, pbc=.TRUE., transpose=.FALSE.)
386 180 : CALL qs_rho_get(nablavks_set(idir1, ispin)%rho, rho_r=nrho1_r)
387 : CALL pw_spline_do_precond(precond, nrho1_r(1), &
388 180 : vks_pw_spline(idir1, ispin))
389 180 : CALL pw_spline_precond_set_kind(precond, precond_kind)
390 : success = find_coeffs(values=nrho1_r(1), &
391 : coeffs=vks_pw_spline(idir1, ispin), linOp=spl3_pbc, &
392 : preconditioner=precond, pool=auxbas_pw_pool, &
393 180 : eps_r=eps_r, eps_x=eps_x, max_iter=max_iter)
394 180 : CPASSERT(success)
395 240 : CALL pw_spline_precond_release(precond)
396 : END DO ! idir1
397 : END DO ! ispin
398 :
399 30 : temp_so_gapw = 0.0_dp
400 :
401 90 : DO ikind = 1, nkind
402 60 : NULLIFY (atom_list, grid_atom, harmonics)
403 60 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
404 : CALL get_qs_kind(qs_kind_set(ikind), &
405 : hard_radius=hard_radius, &
406 : grid_atom=grid_atom, &
407 : harmonics=harmonics, &
408 60 : paw_atom=paw_atom)
409 :
410 60 : IF (.NOT. paw_atom) CYCLE
411 :
412 : ! Distribute the atoms of this kind
413 :
414 60 : bo = get_limit(natom, para_env%num_pe, para_env%mepos)
415 :
416 240 : DO iat = 1, natom !bo(1),bo(2)! this partitioning blocks the interpolation
417 : ! ! routines (i.e. waiting for parallel sum)
418 90 : iatom = atom_list(iat)
419 90 : NULLIFY (jrho1_atom, nablavks_atom)
420 90 : jrho1_atom => jrho1_atom_set(iatom)
421 90 : nablavks_atom => nablavks_atom_set(iatom)
422 420 : DO idir1 = 1, 3
423 270 : CALL set_vecp(idir1, idir2, idir3)
424 900 : DO ispin = 1, nspins
425 27810 : DO ir = 1, grid_atom%nr
426 :
427 27000 : IF (grid_atom%rad(ir) >= hard_radius) CYCLE
428 :
429 771660 : DO ia = 1, grid_atom%ng_sphere
430 :
431 3024000 : ra = particle_set(iatom)%r
432 3024000 : ra(:) = ra(:) + grid_atom%rad(ir)*harmonics%a(:, ia)
433 : vks_ra_idir2 = Eval_Interp_Spl3_pbc(ra, &
434 756000 : vks_pw_spline(idir2, ispin))
435 : vks_ra_idir3 = Eval_Interp_Spl3_pbc(ra, &
436 756000 : vks_pw_spline(idir3, ispin))
437 :
438 756000 : IF (iat .LT. bo(1) .OR. iat .GT. bo(2)) CYCLE !quick and dirty:
439 : ! !here take care of the partition
440 :
441 : ! + sum_A j_loc_h_A x nabla_vks_s_A
442 : temp_so_gapw(iB, idir1) = temp_so_gapw(iB, idir1) + &
443 : (-1.0_dp)**(1 + ispin)*( &
444 : jrho1_atom%jrho_vec_rad_h(idir2, ispin)%r_coef(ir, ia)* &
445 : vks_ra_idir3 - &
446 : jrho1_atom%jrho_vec_rad_h(idir3, ispin)%r_coef(ir, ia)* &
447 : vks_ra_idir2 &
448 378000 : )*grid_atom%wr(ir)*grid_atom%wa(ia)
449 :
450 : ! - sum_A j_loc_s_A x nabla_vks_s_A
451 : temp_so_gapw(iB, idir1) = temp_so_gapw(iB, idir1) - &
452 : (-1.0_dp)**(1 + ispin)*( &
453 : jrho1_atom%jrho_vec_rad_s(idir2, ispin)%r_coef(ir, ia)* &
454 : vks_ra_idir3 - &
455 : jrho1_atom%jrho_vec_rad_s(idir3, ispin)%r_coef(ir, ia)* &
456 : vks_ra_idir2 &
457 378000 : )*grid_atom%wr(ir)*grid_atom%wa(ia)
458 :
459 : ! + sum_A j_loc_h_A x nabla_vks_loc_h_A
460 : temp_so_gapw(iB, idir1) = temp_so_gapw(iB, idir1) + &
461 : (-1.0_dp)**(1 + ispin)*( &
462 : jrho1_atom%jrho_vec_rad_h(idir2, ispin)%r_coef(ir, ia)* &
463 : nablavks_atom%nablavks_vec_rad_h(idir3, ispin)%r_coef(ir, ia) - &
464 : jrho1_atom%jrho_vec_rad_h(idir3, ispin)%r_coef(ir, ia)* &
465 : nablavks_atom%nablavks_vec_rad_h(idir2, ispin)%r_coef(ir, ia) &
466 378000 : )*grid_atom%wr(ir)*grid_atom%wa(ia)
467 :
468 : ! - sum_A j_loc_h_A x nabla_vks_loc_s_A
469 : temp_so_gapw(iB, idir1) = temp_so_gapw(iB, idir1) - &
470 : (-1.0_dp)**(1 + ispin)*( &
471 : jrho1_atom%jrho_vec_rad_h(idir2, ispin)%r_coef(ir, ia)* &
472 : nablavks_atom%nablavks_vec_rad_s(idir3, ispin)%r_coef(ir, ia) - &
473 : jrho1_atom%jrho_vec_rad_h(idir3, ispin)%r_coef(ir, ia)* &
474 : nablavks_atom%nablavks_vec_rad_s(idir2, ispin)%r_coef(ir, ia) &
475 783000 : )*grid_atom%wr(ir)*grid_atom%wa(ia)
476 :
477 : ! ORIGINAL
478 : ! ! + sum_A j_loc_h_A x nabla_vks_loc_h_A
479 : ! temp_so_gapw(iB,idir1) = temp_so_gapw(iB,idir1) + &
480 : ! (-1.0_dp)**(1.0_dp + REAL(ispin,KIND=dp)) * ( &
481 : ! jrho1_atom%jrho_vec_rad_h(idir2,iB,ispin)%r_coef(ir,ia) * &
482 : ! nablavks_atom%nablavks_vec_rad_h(idir3,ispin)%r_coef(ir,ia) - &
483 : ! jrho1_atom%jrho_vec_rad_h(idir3,iB,ispin)%r_coef(ir,ia) * &
484 : ! nablavks_atom%nablavks_vec_rad_h(idir2,ispin)%r_coef(ir,ia) &
485 : ! ) * grid_atom%wr(ir)*grid_atom%wa(ia)
486 : ! ! - sum_A j_loc_s_A x nabla_vks_loc_s_A
487 : ! temp_so_gapw(iB,idir1) = temp_so_gapw(iB,idir1) - &
488 : ! (-1.0_dp)**(1.0_dp + REAL(ispin,KIND=dp)) * ( &
489 : ! jrho1_atom%jrho_vec_rad_s(idir2,iB,ispin)%r_coef(ir,ia) * &
490 : ! nablavks_atom%nablavks_vec_rad_s(idir3,ispin)%r_coef(ir,ia) - &
491 : ! jrho1_atom%jrho_vec_rad_s(idir3,iB,ispin)%r_coef(ir,ia) * &
492 : ! nablavks_atom%nablavks_vec_rad_s(idir2,ispin)%r_coef(ir,ia) &
493 : ! ) * grid_atom%wr(ir)*grid_atom%wa(ia)
494 : END DO !ia
495 : END DO !ir
496 : END DO !ispin
497 : END DO !idir1
498 : END DO !iat
499 : END DO !ikind
500 :
501 30 : CALL para_env%sum(temp_so_gapw)
502 390 : temp_so_gapw(:, :) = -1.0_dp*epr_env%g_so_factor_gapw*temp_so_gapw(:, :)
503 :
504 30 : IF (output_unit > 0) THEN
505 60 : DO idir1 = 1, 3
506 : WRITE (UNIT=output_unit, FMT="(T2,A,T18,I1,I1,T56,E24.16)") &
507 60 : "epr|SOX:gapw", iB, idir1, temp_so_gapw(iB, idir1)
508 : END DO
509 : END IF
510 :
511 120 : g_so(iB, :) = g_so(iB, :) + temp_so_gapw(iB, :)
512 :
513 90 : DO ispin = 1, nspins
514 270 : DO idir1 = 1, 3
515 240 : CALL auxbas_pw_pool%give_back_pw(vks_pw_spline(idir1, ispin))
516 : END DO
517 : END DO
518 60 : DEALLOCATE (vks_pw_spline)
519 :
520 : END IF ! gapw
521 :
522 294 : g_total(iB, :) = g_total(iB, :) + g_so(iB, :)
523 :
524 : CALL cp_print_key_finished_output(output_unit, logger, lr_section, &
525 42 : "PRINT%PROGRAM_RUN_INFO")
526 :
527 336 : END SUBROUTINE epr_g_so
528 :
529 : ! **************************************************************************************************
530 : !> \brief Calculates g_soo (soft part only for now)
531 : !> \param epr_env ...
532 : !> \param current_env ...
533 : !> \param qs_env ...
534 : !> \param iB ...
535 : !> \par History
536 : !> 06.2006 created [RD]
537 : !> \author RD
538 : ! **************************************************************************************************
539 42 : SUBROUTINE epr_g_soo(epr_env, current_env, qs_env, iB)
540 :
541 : TYPE(epr_env_type) :: epr_env
542 : TYPE(current_env_type) :: current_env
543 : TYPE(qs_environment_type), POINTER :: qs_env
544 : INTEGER, INTENT(IN) :: iB
545 :
546 : INTEGER :: aint_precond, ia, iat, iatom, idir1, &
547 : ikind, ir, iso, ispin, max_iter, &
548 : natom, nkind, nspins, output_unit, &
549 : precond_kind
550 : INTEGER, DIMENSION(2) :: bo
551 42 : INTEGER, DIMENSION(:), POINTER :: atom_list
552 : LOGICAL :: gapw, paw_atom, soo_rho_hard, success
553 : REAL(dp) :: bind_ra_idir1, chi_tensor(3, 3, 2), &
554 : eps_r, eps_x, hard_radius, rho_spin, &
555 : temp_soo_soft
556 : REAL(dp), DIMENSION(3, 3) :: temp_soo_gapw
557 42 : REAL(dp), DIMENSION(:, :), POINTER :: g_soo, g_total
558 : REAL(KIND=dp), DIMENSION(3) :: ra
559 42 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
560 : TYPE(cp_logger_type), POINTER :: logger
561 : TYPE(dft_control_type), POINTER :: dft_control
562 : TYPE(grid_atom_type), POINTER :: grid_atom
563 : TYPE(harmonics_atom_type), POINTER :: harmonics
564 : TYPE(mp_para_env_type), POINTER :: para_env
565 42 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
566 : TYPE(pw_env_type), POINTER :: pw_env
567 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
568 42 : TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:, :) :: bind_pw_spline
569 84 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: brho1_r, rho_r
570 : TYPE(pw_spline_precond_type) :: precond
571 42 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
572 42 : TYPE(qs_rho_p_type), DIMENSION(:, :), POINTER :: bind_set
573 : TYPE(qs_rho_type), POINTER :: rho
574 42 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: rho_rad_h, rho_rad_s
575 42 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
576 : TYPE(rho_atom_type), POINTER :: rho_atom
577 : TYPE(section_vals_type), POINTER :: g_section, interp_section, lr_section
578 :
579 42 : NULLIFY (atomic_kind_set, qs_kind_set, atom_list, bind_set, dft_control, &
580 42 : grid_atom, g_section, g_soo, g_total, harmonics, interp_section, &
581 42 : logger, lr_section, para_env, particle_set, rho, rho_atom, &
582 42 : rho_atom_set, rho_r, brho1_r)
583 :
584 84 : logger => cp_get_default_logger()
585 42 : lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES")
586 :
587 : output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
588 42 : extension=".linresLog")
589 :
590 : g_section => section_vals_get_subs_vals(lr_section, &
591 42 : "EPR%PRINT%G_TENSOR")
592 :
593 42 : CALL section_vals_val_get(g_section, "soo_rho_hard", l_val=soo_rho_hard)
594 :
595 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
596 : dft_control=dft_control, para_env=para_env, particle_set=particle_set, &
597 42 : pw_env=pw_env, rho=rho, rho_atom_set=rho_atom_set)
598 :
599 : CALL get_epr_env(epr_env=epr_env, bind_set=bind_set, &
600 42 : g_soo=g_soo, g_total=g_total)
601 :
602 : CALL get_current_env(current_env=current_env, &
603 42 : chi_tensor=chi_tensor)
604 42 : CALL qs_rho_get(rho, rho_r=rho_r)
605 :
606 42 : gapw = dft_control%qs_control%gapw
607 42 : nkind = SIZE(qs_kind_set, 1)
608 42 : nspins = dft_control%nspins
609 :
610 168 : DO idir1 = 1, 3
611 126 : temp_soo_soft = 0.0_dp
612 378 : DO ispin = 1, nspins
613 252 : CALL qs_rho_get(bind_set(idir1, iB)%rho, rho_r=brho1_r)
614 : temp_soo_soft = temp_soo_soft + (-1.0_dp)**(1 + ispin)* &
615 378 : pw_integral_ab(brho1_r(1), rho_r(ispin))
616 : END DO
617 126 : temp_soo_soft = 1.0_dp*epr_env%g_soo_factor*temp_soo_soft
618 126 : IF (output_unit > 0) THEN
619 : WRITE (UNIT=output_unit, FMT="(T2,A,T18,i1,i1,T56,E24.16)") &
620 63 : "epr|SOO:soft", iB, idir1, temp_soo_soft
621 : END IF
622 168 : g_soo(iB, idir1) = temp_soo_soft
623 : END DO
624 :
625 168 : DO idir1 = 1, 3
626 : temp_soo_soft = 1.0_dp*epr_env%g_soo_chicorr_factor*chi_tensor(idir1, iB, 2)* &
627 126 : (REAL(dft_control%multiplicity, KIND=dp) - 1.0_dp)
628 126 : IF (output_unit > 0) THEN
629 : WRITE (UNIT=output_unit, FMT="(T2,A,T18,i1,i1,T56,E24.16)") &
630 63 : "epr|SOO:soft_g0", iB, idir1, temp_soo_soft
631 : END IF
632 168 : g_soo(iB, idir1) = g_soo(iB, idir1) + temp_soo_soft
633 : END DO
634 :
635 42 : IF (gapw .AND. soo_rho_hard) THEN
636 :
637 12 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
638 156 : ALLOCATE (bind_pw_spline(3, 3))
639 :
640 : interp_section => section_vals_get_subs_vals(lr_section, &
641 12 : "EPR%INTERPOLATOR")
642 : CALL section_vals_val_get(interp_section, "aint_precond", &
643 12 : i_val=aint_precond)
644 12 : CALL section_vals_val_get(interp_section, "precond", i_val=precond_kind)
645 12 : CALL section_vals_val_get(interp_section, "max_iter", i_val=max_iter)
646 12 : CALL section_vals_val_get(interp_section, "eps_r", r_val=eps_r)
647 12 : CALL section_vals_val_get(interp_section, "eps_x", r_val=eps_x)
648 :
649 48 : DO idir1 = 1, 3
650 36 : CALL auxbas_pw_pool%create_pw(bind_pw_spline(idir1, iB))
651 : ! calculate spline coefficients
652 : CALL pw_spline_precond_create(precond, precond_kind=aint_precond, &
653 36 : pool=auxbas_pw_pool, pbc=.TRUE., transpose=.FALSE.)
654 36 : CALL qs_rho_get(bind_set(idir1, iB)%rho, rho_r=brho1_r)
655 : CALL pw_spline_do_precond(precond, brho1_r(1), &
656 36 : bind_pw_spline(idir1, iB))
657 36 : CALL pw_spline_precond_set_kind(precond, precond_kind)
658 : success = find_coeffs(values=brho1_r(1), &
659 : coeffs=bind_pw_spline(idir1, iB), linOp=spl3_pbc, &
660 : preconditioner=precond, pool=auxbas_pw_pool, &
661 36 : eps_r=eps_r, eps_x=eps_x, max_iter=max_iter)
662 36 : CPASSERT(success)
663 48 : CALL pw_spline_precond_release(precond)
664 : END DO ! idir1
665 :
666 12 : temp_soo_gapw = 0.0_dp
667 :
668 36 : DO ikind = 1, nkind
669 24 : NULLIFY (atom_list, grid_atom, harmonics)
670 24 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
671 : CALL get_qs_kind(qs_kind_set(ikind), &
672 : hard_radius=hard_radius, &
673 : grid_atom=grid_atom, &
674 : harmonics=harmonics, &
675 24 : paw_atom=paw_atom)
676 :
677 24 : IF (.NOT. paw_atom) CYCLE
678 :
679 : ! Distribute the atoms of this kind
680 :
681 24 : bo = get_limit(natom, para_env%num_pe, para_env%mepos)
682 :
683 96 : DO iat = 1, natom !bo(1),bo(2)! this partitioning blocks the interpolation
684 : ! ! routines (i.e. waiting for parallel sum)
685 36 : iatom = atom_list(iat)
686 36 : rho_atom => rho_atom_set(iatom)
687 36 : NULLIFY (rho_rad_h, rho_rad_s)
688 : CALL get_rho_atom(rho_atom=rho_atom, rho_rad_h=rho_rad_h, &
689 36 : rho_rad_s=rho_rad_s)
690 168 : DO idir1 = 1, 3
691 360 : DO ispin = 1, nspins
692 11124 : DO ir = 1, grid_atom%nr
693 :
694 10800 : IF (grid_atom%rad(ir) >= hard_radius) CYCLE
695 :
696 308664 : DO ia = 1, grid_atom%ng_sphere
697 :
698 1209600 : ra = particle_set(iatom)%r
699 1209600 : ra(:) = ra(:) + grid_atom%rad(ir)*harmonics%a(:, ia)
700 : bind_ra_idir1 = Eval_Interp_Spl3_pbc(ra, &
701 302400 : bind_pw_spline(idir1, iB))
702 :
703 302400 : IF (iat .LT. bo(1) .OR. iat .GT. bo(2)) CYCLE !quick and dirty:
704 : ! !here take care of the partition
705 :
706 151200 : rho_spin = 0.0_dp
707 :
708 2721600 : DO iso = 1, harmonics%max_iso_not0
709 : rho_spin = rho_spin + &
710 : (rho_rad_h(ispin)%r_coef(ir, iso) - &
711 : rho_rad_s(ispin)%r_coef(ir, iso))* &
712 2721600 : harmonics%slm(ia, iso)
713 : END DO
714 :
715 : temp_soo_gapw(iB, idir1) = temp_soo_gapw(iB, idir1) + &
716 : (-1.0_dp)**(1 + ispin)*( &
717 : bind_ra_idir1*rho_spin &
718 313200 : )*grid_atom%wr(ir)*grid_atom%wa(ia)
719 :
720 : END DO !ia
721 : END DO !ir
722 : END DO ! ispin
723 : END DO !idir1
724 : END DO !iat
725 : END DO !ikind
726 :
727 12 : CALL para_env%sum(temp_soo_gapw)
728 156 : temp_soo_gapw(:, :) = 1.0_dp*epr_env%g_soo_factor*temp_soo_gapw(:, :)
729 :
730 12 : IF (output_unit > 0) THEN
731 24 : DO idir1 = 1, 3
732 : WRITE (UNIT=output_unit, FMT="(T2,A,T18,I1,I1,T56,E24.16)") &
733 24 : "epr|SOO:gapw", iB, idir1, temp_soo_gapw(iB, idir1)
734 : END DO
735 : END IF
736 :
737 48 : g_soo(iB, :) = g_soo(iB, :) + temp_soo_gapw(iB, :)
738 :
739 48 : DO idir1 = 1, 3
740 48 : CALL auxbas_pw_pool%give_back_pw(bind_pw_spline(idir1, iB))
741 : END DO
742 24 : DEALLOCATE (bind_pw_spline)
743 :
744 : END IF ! gapw
745 :
746 294 : g_total(iB, :) = g_total(iB, :) + g_soo(iB, :)
747 :
748 : CALL cp_print_key_finished_output(output_unit, logger, lr_section, &
749 42 : "PRINT%PROGRAM_RUN_INFO")
750 :
751 336 : END SUBROUTINE epr_g_soo
752 :
753 : ! **************************************************************************************************
754 : !> \brief ...
755 : !> \param epr_env ...
756 : !> \param current_env ...
757 : !> \param qs_env ...
758 : !> \param iB ...
759 : ! **************************************************************************************************
760 42 : SUBROUTINE epr_ind_magnetic_field(epr_env, current_env, qs_env, iB)
761 :
762 : TYPE(epr_env_type) :: epr_env
763 : TYPE(current_env_type) :: current_env
764 : TYPE(qs_environment_type), POINTER :: qs_env
765 : INTEGER, INTENT(IN) :: iB
766 :
767 : CHARACTER(LEN=*), PARAMETER :: routineN = 'epr_ind_magnetic_field'
768 :
769 : INTEGER :: handle, idir, idir2, idir3, iiB, iiiB, &
770 : ispin, natom, nspins
771 : LOGICAL :: gapw
772 : REAL(dp) :: scale_fac
773 : TYPE(cell_type), POINTER :: cell
774 : TYPE(dft_control_type), POINTER :: dft_control
775 42 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
776 : TYPE(pw_c1d_gs_type) :: pw_gspace_work
777 42 : TYPE(pw_c1d_gs_type), ALLOCATABLE, DIMENSION(:, :) :: shift_pw_gspace
778 42 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: jrho1_g
779 : TYPE(pw_env_type), POINTER :: pw_env
780 42 : TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
781 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
782 : TYPE(pw_r3d_rs_type) :: shift_pw_rspace
783 42 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: epr_rho_r
784 : TYPE(realspace_grid_desc_type), POINTER :: auxbas_rs_desc
785 :
786 42 : CALL timeset(routineN, handle)
787 :
788 42 : NULLIFY (cell, dft_control, pw_env, auxbas_rs_desc, auxbas_pw_pool, &
789 42 : pw_pools, particle_set, jrho1_g, epr_rho_r)
790 :
791 : CALL get_qs_env(qs_env=qs_env, cell=cell, dft_control=dft_control, &
792 42 : particle_set=particle_set)
793 :
794 42 : gapw = dft_control%qs_control%gapw
795 42 : natom = SIZE(particle_set, 1)
796 42 : nspins = dft_control%nspins
797 :
798 42 : CALL get_epr_env(epr_env=epr_env)
799 :
800 42 : CALL get_current_env(current_env=current_env)
801 :
802 42 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
803 : CALL pw_env_get(pw_env, auxbas_rs_desc=auxbas_rs_desc, &
804 42 : auxbas_pw_pool=auxbas_pw_pool, pw_pools=pw_pools)
805 : !
806 : ! Initialize
807 : ! Allocate grids for the calculation of jrho and the shift
808 462 : ALLOCATE (shift_pw_gspace(3, nspins))
809 126 : DO ispin = 1, nspins
810 378 : DO idir = 1, 3
811 252 : CALL auxbas_pw_pool%create_pw(shift_pw_gspace(idir, ispin))
812 336 : CALL pw_zero(shift_pw_gspace(idir, ispin))
813 : END DO
814 : END DO
815 42 : CALL auxbas_pw_pool%create_pw(shift_pw_rspace)
816 42 : CALL pw_zero(shift_pw_rspace)
817 42 : CALL auxbas_pw_pool%create_pw(pw_gspace_work)
818 42 : CALL pw_zero(pw_gspace_work)
819 : !
820 42 : CALL set_vecp(iB, iiB, iiiB)
821 : !
822 126 : DO ispin = 1, nspins
823 : !
824 336 : DO idir3 = 1, 3
825 : ! set to zero for the calculation of the shift
826 336 : CALL pw_zero(shift_pw_gspace(idir3, ispin))
827 : END DO
828 378 : DO idir = 1, 3
829 252 : CALL qs_rho_get(current_env%jrho1_set(idir)%rho, rho_g=jrho1_g)
830 : ! Field gradient
831 : ! loop over the Gvec components: x,y,z
832 1092 : DO idir2 = 1, 3
833 1008 : IF (idir /= idir2) THEN
834 : ! in reciprocal space multiply (G_idir2(i)/G(i)^2)J_(idir)(G(i))
835 : CALL mult_G_ov_G2_grid(auxbas_pw_pool, jrho1_g(ispin), &
836 504 : pw_gspace_work, idir2, 0.0_dp)
837 : !
838 : ! scale and add to the correct component of the shift column
839 504 : CALL set_vecp_rev(idir, idir2, idir3)
840 504 : scale_fac = fac_vecp(idir3, idir2, idir)
841 504 : CALL pw_scale(pw_gspace_work, scale_fac)
842 504 : CALL pw_axpy(pw_gspace_work, shift_pw_gspace(idir3, ispin))
843 : END IF
844 : END DO
845 : !
846 : END DO ! idir
847 : END DO ! ispin
848 :
849 : ! Store the total soft induced magnetic field (corrected for sic)
850 42 : IF (dft_control%nspins == 2) THEN
851 168 : DO idir = 1, 3
852 126 : CALL qs_rho_get(epr_env%bind_set(idir, iB)%rho, rho_r=epr_rho_r)
853 168 : CALL pw_transfer(shift_pw_gspace(idir, 2), epr_rho_r(1))
854 : END DO
855 : END IF
856 : !
857 : ! Dellocate grids for the calculation of jrho and the shift
858 42 : CALL auxbas_pw_pool%give_back_pw(pw_gspace_work)
859 126 : DO ispin = 1, dft_control%nspins
860 378 : DO idir = 1, 3
861 336 : CALL auxbas_pw_pool%give_back_pw(shift_pw_gspace(idir, ispin))
862 : END DO
863 : END DO
864 42 : DEALLOCATE (shift_pw_gspace)
865 42 : CALL auxbas_pw_pool%give_back_pw(shift_pw_rspace)
866 : !
867 : ! Finalize
868 42 : CALL timestop(handle)
869 : !
870 84 : END SUBROUTINE epr_ind_magnetic_field
871 :
872 : END MODULE qs_linres_epr_ownutils
873 :
|
