Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief from the response current density calculates the shift tensor
10 : !> and the susceptibility
11 : !> \par History
12 : !> created 02-2006 [MI]
13 : !> \author MI
14 : ! **************************************************************************************************
15 : MODULE qs_linres_nmr_shift
16 : USE atomic_kind_types, ONLY: atomic_kind_type,&
17 : get_atomic_kind
18 : USE cell_types, ONLY: cell_type,&
19 : pbc
20 : USE cp_control_types, ONLY: dft_control_type
21 : USE cp_log_handling, ONLY: cp_get_default_logger,&
22 : cp_logger_get_default_io_unit,&
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 input_section_types, ONLY: section_vals_get_subs_vals,&
29 : section_vals_type,&
30 : section_vals_val_get
31 : USE kinds, ONLY: default_string_length,&
32 : dp
33 : USE mathconstants, ONLY: gaussi,&
34 : twopi,&
35 : z_zero
36 : USE mathlib, ONLY: diamat_all
37 : USE message_passing, ONLY: mp_para_env_type
38 : USE particle_types, ONLY: particle_type
39 : USE pw_env_types, ONLY: pw_env_get,&
40 : pw_env_type
41 : USE pw_grid_types, ONLY: pw_grid_type
42 : USE pw_methods, ONLY: pw_axpy,&
43 : pw_transfer,&
44 : pw_zero
45 : USE pw_pool_types, ONLY: pw_pool_p_type,&
46 : pw_pool_type
47 : USE pw_spline_utils, ONLY: Eval_Interp_Spl3_pbc,&
48 : find_coeffs,&
49 : pw_spline_do_precond,&
50 : pw_spline_precond_create,&
51 : pw_spline_precond_release,&
52 : pw_spline_precond_set_kind,&
53 : pw_spline_precond_type,&
54 : spl3_pbc
55 : USE pw_types, ONLY: pw_c1d_gs_type,&
56 : pw_r3d_rs_type
57 : USE qs_environment_types, ONLY: get_qs_env,&
58 : qs_environment_type
59 : USE qs_grid_atom, ONLY: grid_atom_type
60 : USE qs_harmonics_atom, ONLY: harmonics_atom_type
61 : USE qs_kind_types, ONLY: get_qs_kind,&
62 : qs_kind_type
63 : USE qs_linres_nmr_epr_common_utils, ONLY: mult_G_ov_G2_grid
64 : USE qs_linres_op, ONLY: fac_vecp,&
65 : set_vecp,&
66 : set_vecp_rev
67 : USE qs_linres_types, ONLY: current_env_type,&
68 : get_current_env,&
69 : get_nmr_env,&
70 : jrho_atom_type,&
71 : nmr_env_type
72 : USE qs_rho_types, ONLY: qs_rho_get
73 : USE realspace_grid_types, ONLY: realspace_grid_desc_type
74 : USE util, ONLY: get_limit
75 : #include "./base/base_uses.f90"
76 :
77 : IMPLICIT NONE
78 :
79 : PRIVATE
80 :
81 : ! *** Public subroutines ***
82 : PUBLIC :: nmr_shift_print, &
83 : nmr_shift
84 :
85 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_linres_nmr_shift'
86 :
87 : ! **************
88 :
89 : CONTAINS
90 :
91 : ! **************************************************************************************************
92 : !> \brief ...
93 : !> \param nmr_env ...
94 : !> \param current_env ...
95 : !> \param qs_env ...
96 : !> \param iB ...
97 : ! **************************************************************************************************
98 480 : SUBROUTINE nmr_shift(nmr_env, current_env, qs_env, iB)
99 :
100 : TYPE(nmr_env_type) :: nmr_env
101 : TYPE(current_env_type) :: current_env
102 : TYPE(qs_environment_type), POINTER :: qs_env
103 : INTEGER, INTENT(IN) :: iB
104 :
105 : CHARACTER(LEN=*), PARAMETER :: routineN = 'nmr_shift'
106 :
107 : INTEGER :: handle, idir, idir2, idir3, iiB, iiiB, &
108 : ispin, natom, nspins
109 : LOGICAL :: gapw, interpolate_shift
110 : REAL(dp) :: scale_fac
111 480 : REAL(dp), DIMENSION(:, :, :), POINTER :: chemical_shift, chemical_shift_loc, &
112 480 : chemical_shift_nics, &
113 480 : chemical_shift_nics_loc
114 : TYPE(cell_type), POINTER :: cell
115 : TYPE(dft_control_type), POINTER :: dft_control
116 480 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
117 : TYPE(pw_c1d_gs_type) :: pw_gspace_work
118 480 : TYPE(pw_c1d_gs_type), ALLOCATABLE, DIMENSION(:, :) :: shift_pw_gspace
119 480 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: jrho1_g
120 : TYPE(pw_env_type), POINTER :: pw_env
121 480 : TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
122 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
123 : TYPE(pw_r3d_rs_type) :: shift_pw_rspace
124 : TYPE(realspace_grid_desc_type), POINTER :: auxbas_rs_desc
125 : TYPE(section_vals_type), POINTER :: nmr_section
126 :
127 480 : CALL timeset(routineN, handle)
128 :
129 480 : NULLIFY (chemical_shift, chemical_shift_loc, chemical_shift_nics, &
130 480 : chemical_shift_nics_loc, cell, dft_control, pw_env, &
131 480 : auxbas_rs_desc, auxbas_pw_pool, pw_pools, particle_set, jrho1_g)
132 :
133 : CALL get_qs_env(qs_env=qs_env, cell=cell, dft_control=dft_control, &
134 480 : particle_set=particle_set)
135 :
136 480 : gapw = dft_control%qs_control%gapw
137 480 : natom = SIZE(particle_set, 1)
138 480 : nspins = dft_control%nspins
139 :
140 : CALL get_nmr_env(nmr_env=nmr_env, chemical_shift=chemical_shift, &
141 : chemical_shift_loc=chemical_shift_loc, &
142 : chemical_shift_nics=chemical_shift_nics, &
143 : chemical_shift_nics_loc=chemical_shift_nics_loc, &
144 480 : interpolate_shift=interpolate_shift)
145 :
146 480 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
147 : CALL pw_env_get(pw_env, auxbas_rs_desc=auxbas_rs_desc, &
148 480 : auxbas_pw_pool=auxbas_pw_pool, pw_pools=pw_pools)
149 : !
150 : !
151 : nmr_section => section_vals_get_subs_vals(qs_env%input, &
152 480 : & "PROPERTIES%LINRES%NMR")
153 : !
154 : ! Initialize
155 : ! Allocate grids for the calculation of jrho and the shift
156 4104 : ALLOCATE (shift_pw_gspace(3, nspins))
157 1146 : DO ispin = 1, nspins
158 3144 : DO idir = 1, 3
159 1998 : CALL auxbas_pw_pool%create_pw(shift_pw_gspace(idir, ispin))
160 2664 : CALL pw_zero(shift_pw_gspace(idir, ispin))
161 : END DO
162 : END DO
163 : !
164 : !
165 480 : CALL set_vecp(iB, iiB, iiiB)
166 : !
167 480 : CALL auxbas_pw_pool%create_pw(pw_gspace_work)
168 480 : CALL pw_zero(pw_gspace_work)
169 1146 : DO ispin = 1, nspins
170 : !
171 3144 : DO idir = 1, 3
172 1998 : CALL qs_rho_get(current_env%jrho1_set(idir)%rho, rho_g=jrho1_g)
173 : ! Field gradient
174 : ! loop over the Gvec components: x,y,z
175 8658 : DO idir2 = 1, 3
176 7992 : IF (idir /= idir2) THEN
177 : ! in reciprocal space multiply (G_idir2(i)/G(i)^2)J_(idir)(G(i))
178 : CALL mult_G_ov_G2_grid(auxbas_pw_pool, jrho1_g(ispin), &
179 3996 : pw_gspace_work, idir2, 0.0_dp)
180 : !
181 : ! scale and add to the correct component of the shift column
182 3996 : CALL set_vecp_rev(idir, idir2, idir3)
183 3996 : scale_fac = fac_vecp(idir3, idir2, idir)
184 3996 : CALL pw_axpy(pw_gspace_work, shift_pw_gspace(idir3, ispin), scale_fac)
185 : END IF
186 : END DO
187 : !
188 : END DO ! idir
189 : END DO ! ispin
190 : !
191 480 : CALL auxbas_pw_pool%give_back_pw(pw_gspace_work)
192 : !
193 : ! compute shildings
194 480 : IF (interpolate_shift) THEN
195 24 : CALL auxbas_pw_pool%create_pw(shift_pw_rspace)
196 60 : DO ispin = 1, nspins
197 168 : DO idir = 1, 3
198 : ! Here first G->R and then interpolation to get the shifts.
199 : ! The interpolation doesnt work in parallel yet.
200 : ! The choice between both methods should be left to the user.
201 108 : CALL pw_transfer(shift_pw_gspace(idir, ispin), shift_pw_rspace)
202 : CALL interpolate_shift_pwgrid(nmr_env, pw_env, particle_set, cell, shift_pw_rspace, &
203 144 : iB, idir, nmr_section)
204 : END DO
205 : END DO
206 24 : CALL auxbas_pw_pool%give_back_pw(shift_pw_rspace)
207 : ELSE
208 1086 : DO ispin = 1, nspins
209 2976 : DO idir = 1, 3
210 : ! Here the shifts are computed from summation of the coeff on the G-grip .
211 : CALL gsum_shift_pwgrid(nmr_env, particle_set, cell, &
212 2520 : shift_pw_gspace(idir, ispin), iB, idir)
213 : END DO
214 : END DO
215 : END IF
216 : !
217 480 : IF (gapw) THEN
218 1032 : DO idir = 1, 3
219 : ! Finally the radial functions are multiplied by the YLM and properly summed
220 : ! The resulting array is J on the local grid. One array per atom.
221 : ! Local contributions by numerical integration over the spherical grids
222 1032 : CALL nmr_shift_gapw(nmr_env, current_env, qs_env, iB, idir)
223 : END DO ! idir
224 : END IF
225 : !
226 : ! Dellocate grids for the calculation of jrho and the shift
227 1146 : DO ispin = 1, nspins
228 3144 : DO idir = 1, 3
229 2664 : CALL auxbas_pw_pool%give_back_pw(shift_pw_gspace(idir, ispin))
230 : END DO
231 : END DO
232 480 : DEALLOCATE (shift_pw_gspace)
233 : !
234 : ! Finalize
235 480 : CALL timestop(handle)
236 : !
237 1440 : END SUBROUTINE nmr_shift
238 :
239 : ! **************************************************************************************************
240 : !> \brief ...
241 : !> \param nmr_env ...
242 : !> \param current_env ...
243 : !> \param qs_env ...
244 : !> \param iB ...
245 : !> \param idir ...
246 : ! **************************************************************************************************
247 774 : SUBROUTINE nmr_shift_gapw(nmr_env, current_env, qs_env, iB, idir)
248 :
249 : TYPE(nmr_env_type) :: nmr_env
250 : TYPE(current_env_type) :: current_env
251 : TYPE(qs_environment_type), POINTER :: qs_env
252 : INTEGER, INTENT(IN) :: IB, idir
253 :
254 : CHARACTER(len=*), PARAMETER :: routineN = 'nmr_shift_gapw'
255 :
256 : INTEGER :: handle, ia, iat, iatom, idir2_1, idir3_1, ikind, ir, ira, ispin, j, jatom, &
257 : n_nics, na, natom, natom_local, natom_tot, nkind, nnics_local, nr, nra, nspins, &
258 : output_unit
259 774 : INTEGER, ALLOCATABLE, DIMENSION(:) :: list_j, list_nics_j
260 : INTEGER, DIMENSION(2) :: bo
261 774 : INTEGER, DIMENSION(:), POINTER :: atom_list
262 : LOGICAL :: do_nics, paw_atom
263 : REAL(dp) :: ddiff, dist, dum, itegrated_jrho, &
264 : r_jatom(3), rdiff(3), rij(3), s_1, &
265 : s_2, scale_fac_1, shift_gapw_radius
266 774 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: j_grid
267 774 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: cs_loc_tmp, cs_nics_loc_tmp, dist_ij, &
268 774 : dist_nics_ij, r_grid
269 774 : REAL(dp), DIMENSION(:, :), POINTER :: jrho_h_grid, jrho_s_grid, r_nics
270 774 : REAL(dp), DIMENSION(:, :, :), POINTER :: chemical_shift_loc, &
271 774 : chemical_shift_nics_loc
272 774 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
273 : TYPE(cell_type), POINTER :: cell
274 : TYPE(cp_logger_type), POINTER :: logger
275 : TYPE(dft_control_type), POINTER :: dft_control
276 : TYPE(grid_atom_type), POINTER :: grid_atom
277 : TYPE(harmonics_atom_type), POINTER :: harmonics
278 774 : TYPE(jrho_atom_type), DIMENSION(:), POINTER :: jrho1_atom_set
279 : TYPE(jrho_atom_type), POINTER :: jrho1_atom
280 : TYPE(mp_para_env_type), POINTER :: para_env
281 774 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
282 774 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
283 :
284 774 : CALL timeset(routineN, handle)
285 : !
286 774 : NULLIFY (atomic_kind_set, qs_kind_set, cell, dft_control, para_env, particle_set, &
287 774 : chemical_shift_loc, chemical_shift_nics_loc, jrho1_atom_set, &
288 774 : jrho1_atom, r_nics, jrho_h_grid, jrho_s_grid, &
289 774 : atom_list, grid_atom, harmonics, logger)
290 : !
291 774 : logger => cp_get_default_logger()
292 774 : output_unit = cp_logger_get_default_io_unit(logger)
293 : !
294 : CALL get_qs_env(qs_env=qs_env, &
295 : atomic_kind_set=atomic_kind_set, &
296 : qs_kind_set=qs_kind_set, &
297 : cell=cell, &
298 : dft_control=dft_control, &
299 : para_env=para_env, &
300 774 : particle_set=particle_set)
301 :
302 : CALL get_nmr_env(nmr_env=nmr_env, &
303 : chemical_shift_loc=chemical_shift_loc, &
304 : chemical_shift_nics_loc=chemical_shift_nics_loc, &
305 : shift_gapw_radius=shift_gapw_radius, &
306 : n_nics=n_nics, &
307 : r_nics=r_nics, &
308 774 : do_nics=do_nics)
309 :
310 : CALL get_current_env(current_env=current_env, &
311 774 : jrho1_atom_set=jrho1_atom_set)
312 : !
313 774 : nkind = SIZE(atomic_kind_set, 1)
314 774 : natom_tot = SIZE(particle_set, 1)
315 774 : nspins = dft_control%nspins
316 774 : itegrated_jrho = 0.0_dp
317 : !
318 774 : idir2_1 = MODULO(idir, 3) + 1
319 774 : idir3_1 = MODULO(idir + 1, 3) + 1
320 774 : scale_fac_1 = fac_vecp(idir3_1, idir2_1, idir)
321 : !
322 : ALLOCATE (cs_loc_tmp(3, natom_tot), list_j(natom_tot), &
323 4644 : dist_ij(3, natom_tot))
324 11070 : cs_loc_tmp = 0.0_dp
325 774 : IF (do_nics) THEN
326 : ALLOCATE (cs_nics_loc_tmp(3, n_nics), list_nics_j(n_nics), &
327 108 : dist_nics_ij(3, n_nics))
328 234 : cs_nics_loc_tmp = 0.0_dp
329 : END IF
330 : !
331 : ! Loop over atoms to collocate the current on each atomic grid, JA
332 : ! Per each JA, loop over the points where the shift needs to be computed
333 2196 : DO ikind = 1, nkind
334 :
335 1422 : NULLIFY (atom_list, grid_atom, harmonics)
336 : CALL get_atomic_kind(atomic_kind_set(ikind), &
337 : atom_list=atom_list, &
338 1422 : natom=natom)
339 :
340 : CALL get_qs_kind(qs_kind_set(ikind), &
341 : paw_atom=paw_atom, &
342 : harmonics=harmonics, &
343 1422 : grid_atom=grid_atom)
344 : !
345 1422 : na = grid_atom%ng_sphere
346 1422 : nr = grid_atom%nr
347 1422 : nra = nr*na
348 7110 : ALLOCATE (r_grid(3, nra), j_grid(nra))
349 69282 : ira = 1
350 69282 : DO ia = 1, na
351 3469482 : DO ir = 1, nr
352 13600800 : r_grid(:, ira) = grid_atom%rad(ir)*harmonics%a(:, ia)
353 3468060 : ira = ira + 1
354 : END DO
355 : END DO
356 : !
357 : ! Quick cycle if needed
358 1422 : IF (paw_atom) THEN
359 : !
360 : ! Distribute the atoms of this kind
361 1062 : bo = get_limit(natom, para_env%num_pe, para_env%mepos)
362 : !
363 1773 : DO iat = bo(1), bo(2)
364 711 : iatom = atom_list(iat)
365 : !
366 : ! find all the atoms within the radius
367 711 : natom_local = 0
368 3150 : DO jatom = 1, natom_tot
369 2439 : rij(:) = pbc(particle_set(iatom)%r, particle_set(jatom)%r, cell)
370 2439 : dist = SQRT(rij(1)*rij(1) + rij(2)*rij(2) + rij(3)*rij(3))
371 3150 : IF (dist <= shift_gapw_radius) THEN
372 2403 : natom_local = natom_local + 1
373 2403 : list_j(natom_local) = jatom
374 9612 : dist_ij(:, natom_local) = rij(:)
375 : END IF
376 : END DO
377 : !
378 : ! ... also for nics
379 711 : IF (do_nics) THEN
380 18 : nnics_local = 0
381 72 : DO jatom = 1, n_nics
382 216 : r_jatom(:) = r_nics(:, jatom)
383 54 : rij(:) = pbc(particle_set(iatom)%r, r_jatom, cell)
384 54 : dist = SQRT(rij(1)*rij(1) + rij(2)*rij(2) + rij(3)*rij(3))
385 72 : IF (dist <= shift_gapw_radius) THEN
386 54 : nnics_local = nnics_local + 1
387 54 : list_nics_j(nnics_local) = jatom
388 216 : dist_nics_ij(:, nnics_local) = rij(:)
389 : END IF
390 : END DO
391 : END IF
392 : !
393 711 : NULLIFY (jrho1_atom, jrho_h_grid, jrho_s_grid)
394 711 : jrho1_atom => jrho1_atom_set(iatom)
395 : !
396 2826 : DO ispin = 1, nspins
397 1053 : jrho_h_grid => jrho1_atom%jrho_vec_rad_h(idir, ispin)%r_coef
398 1053 : jrho_s_grid => jrho1_atom%jrho_vec_rad_s(idir, ispin)%r_coef
399 : !
400 : ! loop over the atoms neighbors of iatom in terms of the current density
401 : ! for each compute the contribution to the shift coming from the
402 : ! local current density at iatom
403 1053 : ira = 1
404 50787 : DO ia = 1, na
405 2548467 : DO ir = 1, nr
406 2497680 : j_grid(ira) = (jrho_h_grid(ir, ia) - jrho_s_grid(ir, ia))*grid_atom%weight(ia, ir)
407 2497680 : itegrated_jrho = itegrated_jrho + j_grid(ira)
408 2547414 : ira = ira + 1
409 : END DO
410 : END DO
411 : !
412 4410 : DO j = 1, natom_local
413 3357 : jatom = list_j(j)
414 13428 : rij(:) = dist_ij(:, j)
415 : !
416 : s_1 = 0.0_dp
417 : s_2 = 0.0_dp
418 7826517 : DO ira = 1, nra
419 : !
420 31292640 : rdiff(:) = rij(:) - r_grid(:, ira)
421 7823160 : ddiff = SQRT(rdiff(1)*rdiff(1) + rdiff(2)*rdiff(2) + rdiff(3)*rdiff(3))
422 7826517 : IF (ddiff > 1.0E-12_dp) THEN
423 7823160 : dum = scale_fac_1*j_grid(ira)/(ddiff*ddiff*ddiff)
424 7823160 : s_1 = s_1 + rdiff(idir2_1)*dum
425 7823160 : s_2 = s_2 + rdiff(idir3_1)*dum
426 : END IF ! ddiff
427 : END DO ! ira
428 3357 : cs_loc_tmp(idir3_1, jatom) = cs_loc_tmp(idir3_1, jatom) + s_1
429 4410 : cs_loc_tmp(idir2_1, jatom) = cs_loc_tmp(idir2_1, jatom) - s_2
430 : END DO ! j
431 : !
432 1764 : IF (do_nics) THEN
433 144 : DO j = 1, nnics_local
434 108 : jatom = list_nics_j(j)
435 432 : rij(:) = dist_nics_ij(:, j)
436 : !
437 : s_1 = 0.0_dp
438 : s_2 = 0.0_dp
439 270108 : DO ira = 1, nra
440 : !
441 1080000 : rdiff(:) = rij(:) - r_grid(:, ira)
442 270000 : ddiff = SQRT(rdiff(1)*rdiff(1) + rdiff(2)*rdiff(2) + rdiff(3)*rdiff(3))
443 270108 : IF (ddiff > 1.0E-12_dp) THEN
444 270000 : dum = scale_fac_1*j_grid(ira)/(ddiff*ddiff*ddiff)
445 270000 : s_1 = s_1 + rdiff(idir2_1)*dum
446 270000 : s_2 = s_2 + rdiff(idir3_1)*dum
447 : END IF ! ddiff
448 : END DO ! ira
449 108 : cs_nics_loc_tmp(idir3_1, jatom) = cs_nics_loc_tmp(idir3_1, jatom) + s_1
450 144 : cs_nics_loc_tmp(idir2_1, jatom) = cs_nics_loc_tmp(idir2_1, jatom) - s_2
451 : END DO ! j
452 : END IF ! do_nics
453 : END DO ! ispin
454 : END DO ! iat
455 : END IF
456 3618 : DEALLOCATE (r_grid, j_grid)
457 : END DO ! ikind
458 : !
459 : !
460 774 : CALL para_env%sum(itegrated_jrho)
461 774 : IF (output_unit > 0) THEN
462 : WRITE (output_unit, '(T2,A,E24.16)') 'Integrated local j_'&
463 387 : &//ACHAR(idir + 119)//ACHAR(iB + 119)//'(r)=', itegrated_jrho
464 : END IF
465 : !
466 774 : CALL para_env%sum(cs_loc_tmp)
467 : chemical_shift_loc(:, iB, :) = chemical_shift_loc(:, iB, :) &
468 11070 : & - nmr_env%shift_factor_gapw*cs_loc_tmp(:, :)/2.0_dp
469 : !
470 774 : DEALLOCATE (cs_loc_tmp, list_j, dist_ij)
471 : !
472 774 : IF (do_nics) THEN
473 18 : CALL para_env%sum(cs_nics_loc_tmp)
474 : chemical_shift_nics_loc(:, iB, :) = chemical_shift_nics_loc(:, iB, :) &
475 234 : & - nmr_env%shift_factor_gapw*cs_nics_loc_tmp(:, :)/2.0_dp
476 : !
477 18 : DEALLOCATE (cs_nics_loc_tmp, list_nics_j, dist_nics_ij)
478 : END IF
479 : !
480 774 : CALL timestop(handle)
481 : !
482 2322 : END SUBROUTINE nmr_shift_gapw
483 :
484 : ! **************************************************************************************************
485 : !> \brief interpolate the shift calculated on the PW grid in order to ger
486 : !> the value on arbitrary points in real space
487 : !> \param nmr_env to get the shift tensor and the list of additional points
488 : !> \param pw_env ...
489 : !> \param particle_set for the atomic position
490 : !> \param cell to take into account the pbs, and to have the volume
491 : !> \param shift_pw_rspace specific component of the shift tensor on the pw grid
492 : !> \param i_B component of the magnetic field for which the shift is calculated (row)
493 : !> \param idir component of the vector \int_{r}[ ((r-r') x j(r))/|r-r'|^3 ] = Bind(r')
494 : !> \param nmr_section ...
495 : !> \author MI
496 : ! **************************************************************************************************
497 864 : SUBROUTINE interpolate_shift_pwgrid(nmr_env, pw_env, particle_set, cell, shift_pw_rspace, &
498 : i_B, idir, nmr_section)
499 :
500 : TYPE(nmr_env_type) :: nmr_env
501 : TYPE(pw_env_type), POINTER :: pw_env
502 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
503 : TYPE(cell_type), POINTER :: cell
504 : TYPE(pw_r3d_rs_type), INTENT(IN) :: shift_pw_rspace
505 : INTEGER, INTENT(IN) :: i_B, idir
506 : TYPE(section_vals_type), POINTER :: nmr_section
507 :
508 : CHARACTER(LEN=*), PARAMETER :: routineN = 'interpolate_shift_pwgrid'
509 :
510 : INTEGER :: aint_precond, handle, iat, iatom, &
511 : max_iter, n_nics, natom, precond_kind
512 108 : INTEGER, DIMENSION(:), POINTER :: cs_atom_list
513 : LOGICAL :: do_nics, success
514 : REAL(dp) :: eps_r, eps_x, R_iatom(3), ra(3), &
515 : shift_val
516 108 : REAL(dp), DIMENSION(:, :), POINTER :: r_nics
517 108 : REAL(dp), DIMENSION(:, :, :), POINTER :: chemical_shift, chemical_shift_nics
518 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
519 : TYPE(pw_r3d_rs_type) :: shiftspl
520 : TYPE(pw_spline_precond_type) :: precond
521 : TYPE(section_vals_type), POINTER :: interp_section
522 :
523 108 : CALL timeset(routineN, handle)
524 : !
525 108 : NULLIFY (interp_section, auxbas_pw_pool)
526 108 : NULLIFY (cs_atom_list, chemical_shift, chemical_shift_nics, r_nics)
527 :
528 : interp_section => section_vals_get_subs_vals(nmr_section, &
529 108 : "INTERPOLATOR")
530 : CALL section_vals_val_get(interp_section, "aint_precond", &
531 108 : i_val=aint_precond)
532 108 : CALL section_vals_val_get(interp_section, "precond", i_val=precond_kind)
533 108 : CALL section_vals_val_get(interp_section, "max_iter", i_val=max_iter)
534 108 : CALL section_vals_val_get(interp_section, "eps_r", r_val=eps_r)
535 108 : CALL section_vals_val_get(interp_section, "eps_x", r_val=eps_x)
536 :
537 : ! calculate spline coefficients
538 108 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
539 108 : CALL auxbas_pw_pool%create_pw(shiftspl)
540 :
541 : CALL pw_spline_precond_create(precond, precond_kind=aint_precond, &
542 108 : pool=auxbas_pw_pool, pbc=.TRUE., transpose=.FALSE.)
543 108 : CALL pw_spline_do_precond(precond, shift_pw_rspace, shiftspl)
544 108 : CALL pw_spline_precond_set_kind(precond, precond_kind)
545 : success = find_coeffs(values=shift_pw_rspace, coeffs=shiftspl, &
546 : linOp=spl3_pbc, preconditioner=precond, pool=auxbas_pw_pool, &
547 108 : eps_r=eps_r, eps_x=eps_x, max_iter=max_iter)
548 108 : CPASSERT(success)
549 108 : CALL pw_spline_precond_release(precond)
550 :
551 : CALL get_nmr_env(nmr_env=nmr_env, cs_atom_list=cs_atom_list, &
552 : chemical_shift=chemical_shift, &
553 : chemical_shift_nics=chemical_shift_nics, &
554 : n_nics=n_nics, r_nics=r_nics, &
555 108 : do_nics=do_nics)
556 :
557 108 : IF (ASSOCIATED(cs_atom_list)) THEN
558 108 : natom = SIZE(cs_atom_list, 1)
559 : ELSE
560 : natom = -1
561 : END IF
562 :
563 378 : DO iat = 1, natom
564 270 : iatom = cs_atom_list(iat)
565 270 : R_iatom = pbc(particle_set(iatom)%r, cell)
566 270 : shift_val = Eval_Interp_Spl3_pbc(R_iatom, shiftspl)
567 : chemical_shift(idir, i_B, iatom) = chemical_shift(idir, i_B, iatom) + &
568 378 : nmr_env%shift_factor*twopi**2*shift_val
569 : END DO
570 :
571 108 : IF (do_nics) THEN
572 144 : DO iatom = 1, n_nics
573 432 : ra(:) = r_nics(:, iatom)
574 108 : R_iatom = pbc(ra, cell)
575 108 : shift_val = Eval_Interp_Spl3_pbc(R_iatom, shiftspl)
576 : chemical_shift_nics(idir, i_B, iatom) = chemical_shift_nics(idir, i_B, iatom) + &
577 144 : nmr_env%shift_factor*twopi**2*shift_val
578 : END DO
579 : END IF
580 :
581 108 : CALL auxbas_pw_pool%give_back_pw(shiftspl)
582 :
583 : !
584 108 : CALL timestop(handle)
585 : !
586 972 : END SUBROUTINE interpolate_shift_pwgrid
587 :
588 : ! **************************************************************************************************
589 : !> \brief ...
590 : !> \param nmr_env ...
591 : !> \param particle_set ...
592 : !> \param cell ...
593 : !> \param shift_pw_gspace ...
594 : !> \param i_B ...
595 : !> \param idir ...
596 : ! **************************************************************************************************
597 3780 : SUBROUTINE gsum_shift_pwgrid(nmr_env, particle_set, cell, shift_pw_gspace, &
598 : i_B, idir)
599 : TYPE(nmr_env_type) :: nmr_env
600 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
601 : TYPE(cell_type), POINTER :: cell
602 : TYPE(pw_c1d_gs_type), INTENT(IN) :: shift_pw_gspace
603 : INTEGER, INTENT(IN) :: i_B, idir
604 :
605 : CHARACTER(LEN=*), PARAMETER :: routineN = 'gsum_shift_pwgrid'
606 :
607 : COMPLEX(dp) :: cplx
608 : INTEGER :: handle, iat, iatom, n_nics, natom
609 1890 : INTEGER, DIMENSION(:), POINTER :: cs_atom_list
610 : LOGICAL :: do_nics
611 : REAL(dp) :: R_iatom(3), ra(3)
612 1890 : REAL(dp), DIMENSION(:, :), POINTER :: r_nics
613 1890 : REAL(dp), DIMENSION(:, :, :), POINTER :: chemical_shift, chemical_shift_nics
614 :
615 1890 : CALL timeset(routineN, handle)
616 : !
617 1890 : NULLIFY (cs_atom_list, chemical_shift, chemical_shift_nics, r_nics)
618 : !
619 : CALL get_nmr_env(nmr_env=nmr_env, cs_atom_list=cs_atom_list, &
620 : chemical_shift=chemical_shift, &
621 : chemical_shift_nics=chemical_shift_nics, &
622 1890 : n_nics=n_nics, r_nics=r_nics, do_nics=do_nics)
623 : !
624 1890 : IF (ASSOCIATED(cs_atom_list)) THEN
625 1890 : natom = SIZE(cs_atom_list, 1)
626 : ELSE
627 : natom = -1
628 : END IF
629 : !
630 : ! compute the chemical shift
631 7254 : DO iat = 1, natom
632 5364 : iatom = cs_atom_list(iat)
633 5364 : R_iatom = pbc(particle_set(iatom)%r, cell)
634 5364 : CALL gsumr(R_iatom, shift_pw_gspace, cplx)
635 : chemical_shift(idir, i_B, iatom) = chemical_shift(idir, i_B, iatom) + &
636 7254 : nmr_env%shift_factor*twopi**2*REAL(cplx, dp)
637 : END DO
638 : !
639 : ! compute nics
640 1890 : IF (do_nics) THEN
641 198 : DO iat = 1, n_nics
642 162 : ra = pbc(r_nics(:, iat), cell)
643 162 : CALL gsumr(ra, shift_pw_gspace, cplx)
644 : chemical_shift_nics(idir, i_B, iat) = chemical_shift_nics(idir, i_B, iat) + &
645 198 : nmr_env%shift_factor*twopi**2*REAL(cplx, dp)
646 : END DO
647 : END IF
648 :
649 1890 : CALL timestop(handle)
650 :
651 1890 : END SUBROUTINE gsum_shift_pwgrid
652 :
653 : ! **************************************************************************************************
654 : !> \brief ...
655 : !> \param r ...
656 : !> \param pw ...
657 : !> \param cplx ...
658 : ! **************************************************************************************************
659 5526 : SUBROUTINE gsumr(r, pw, cplx)
660 : REAL(dp), INTENT(IN) :: r(3)
661 : TYPE(pw_c1d_gs_type), INTENT(IN) :: pw
662 : COMPLEX(dp) :: cplx
663 :
664 : COMPLEX(dp) :: rg
665 : INTEGER :: ig
666 : TYPE(pw_grid_type), POINTER :: grid
667 :
668 5526 : grid => pw%pw_grid
669 5526 : cplx = z_zero
670 102609198 : DO ig = grid%first_gne0, grid%ngpts_cut_local
671 102603672 : rg = (grid%g(1, ig)*r(1) + grid%g(2, ig)*r(2) + grid%g(3, ig)*r(3))*gaussi
672 102609198 : cplx = cplx + pw%array(ig)*EXP(rg)
673 : END DO
674 5526 : IF (grid%have_g0) cplx = cplx + pw%array(1)
675 5526 : CALL grid%para%group%sum(cplx)
676 5526 : END SUBROUTINE gsumr
677 :
678 : ! **************************************************************************************************
679 : !> \brief Shielding tensor and Chi are printed into a file
680 : !> if required from input
681 : !> It is possible to print only for a subset of atoms or
682 : !> or points in non-ionic positions
683 : !> \param nmr_env ...
684 : !> \param current_env ...
685 : !> \param qs_env ...
686 : !> \author MI
687 : ! **************************************************************************************************
688 160 : SUBROUTINE nmr_shift_print(nmr_env, current_env, qs_env)
689 : TYPE(nmr_env_type) :: nmr_env
690 : TYPE(current_env_type) :: current_env
691 : TYPE(qs_environment_type), POINTER :: qs_env
692 :
693 : CHARACTER(LEN=2) :: element_symbol
694 : CHARACTER(LEN=default_string_length) :: name, title
695 : INTEGER :: iatom, ir, n_nics, nat_print, natom, &
696 : output_unit, unit_atoms, unit_nics
697 160 : INTEGER, DIMENSION(:), POINTER :: cs_atom_list
698 : LOGICAL :: do_nics, gapw
699 : REAL(dp) :: chi_aniso, chi_iso, chi_sym_tot(3, 3), chi_tensor(3, 3, 2), &
700 : chi_tensor_loc(3, 3, 2), chi_tensor_loc_tmp(3, 3), chi_tensor_tmp(3, 3), chi_tmp(3, 3), &
701 : eig(3), rpos(3), shift_aniso, shift_iso, shift_sym_tot(3, 3)
702 160 : REAL(dp), DIMENSION(:, :), POINTER :: r_nics
703 160 : REAL(dp), DIMENSION(:, :, :), POINTER :: cs, cs_loc, cs_nics, cs_nics_loc, &
704 160 : cs_nics_tot, cs_tot
705 : REAL(dp), EXTERNAL :: DDOT
706 : TYPE(atomic_kind_type), POINTER :: atom_kind
707 : TYPE(cp_logger_type), POINTER :: logger
708 : TYPE(dft_control_type), POINTER :: dft_control
709 160 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
710 : TYPE(section_vals_type), POINTER :: nmr_section
711 :
712 160 : NULLIFY (cs, cs_nics, r_nics, cs_loc, cs_nics_loc, logger, particle_set, atom_kind, dft_control)
713 :
714 320 : logger => cp_get_default_logger()
715 160 : output_unit = cp_logger_get_default_io_unit(logger)
716 :
717 : nmr_section => section_vals_get_subs_vals(qs_env%input, &
718 160 : "PROPERTIES%LINRES%NMR")
719 :
720 : CALL get_nmr_env(nmr_env=nmr_env, &
721 : chemical_shift=cs, &
722 : chemical_shift_nics=cs_nics, &
723 : chemical_shift_loc=cs_loc, &
724 : chemical_shift_nics_loc=cs_nics_loc, &
725 : cs_atom_list=cs_atom_list, &
726 : n_nics=n_nics, &
727 : r_nics=r_nics, &
728 160 : do_nics=do_nics)
729 : !
730 : CALL get_current_env(current_env=current_env, &
731 : chi_tensor=chi_tensor, &
732 160 : chi_tensor_loc=chi_tensor_loc)
733 : !
734 : ! multiply by the appropriate factor
735 160 : chi_tensor_tmp(:, :) = 0.0_dp
736 160 : chi_tensor_loc_tmp(:, :) = 0.0_dp
737 2080 : chi_tensor_tmp(:, :) = (chi_tensor(:, :, 1) + chi_tensor(:, :, 2))*nmr_env%chi_factor
738 : !chi_tensor_loc_tmp(:,:) = (chi_tensor_loc(:,:,1) + chi_tensor_loc(:,:,2)) * here there is another factor
739 : !
740 : CALL get_qs_env(qs_env=qs_env, &
741 : dft_control=dft_control, &
742 160 : particle_set=particle_set)
743 :
744 160 : natom = SIZE(particle_set, 1)
745 160 : gapw = dft_control%qs_control%gapw
746 160 : nat_print = SIZE(cs_atom_list, 1)
747 :
748 480 : ALLOCATE (cs_tot(3, 3, nat_print))
749 160 : IF (do_nics) THEN
750 18 : ALLOCATE (cs_nics_tot(3, 3, n_nics))
751 : END IF
752 : ! Finalize Chi calculation
753 : ! Symmetrize
754 2080 : chi_sym_tot(:, :) = (chi_tensor_tmp(:, :) + TRANSPOSE(chi_tensor_tmp(:, :)))/2.0_dp
755 160 : IF (gapw) THEN
756 : chi_sym_tot(:, :) = chi_sym_tot(:, :) &
757 1118 : & + (chi_tensor_loc_tmp(:, :) + TRANSPOSE(chi_tensor_loc_tmp(:, :)))/2.0_dp
758 : END IF
759 160 : chi_tmp(:, :) = chi_sym_tot(:, :)
760 160 : CALL diamat_all(chi_tmp, eig)
761 160 : chi_iso = (eig(1) + eig(2) + eig(3))/3.0_dp
762 160 : chi_aniso = eig(3) - (eig(2) + eig(1))/2.0_dp
763 : !
764 160 : IF (output_unit > 0) THEN
765 80 : WRITE (output_unit, '(T2,A,E14.6)') 'CheckSum Chi =', &
766 160 : SQRT(DDOT(9, chi_tensor_tmp(1, 1), 1, chi_tensor_tmp(1, 1), 1))
767 : END IF
768 : !
769 160 : IF (BTEST(cp_print_key_should_output(logger%iter_info, nmr_section, &
770 : "PRINT%CHI_TENSOR"), cp_p_file)) THEN
771 :
772 : unit_atoms = cp_print_key_unit_nr(logger, nmr_section, "PRINT%CHI_TENSOR", &
773 22 : extension=".data", middle_name="CHI", log_filename=.FALSE.)
774 :
775 22 : WRITE (title, '(A)') "Magnetic Susceptibility Tensor "
776 22 : IF (unit_atoms > 0) THEN
777 11 : WRITE (unit_atoms, '(T2,A)') title
778 11 : WRITE (unit_atoms, '(T1,A)') " CHI from SOFT J in 10^-30 J/T^2 units"
779 11 : WRITE (unit_atoms, '(3(A,f10.4))') ' XX = ', chi_tensor_tmp(1, 1),&
780 11 : & ' XY = ', chi_tensor_tmp(1, 2),&
781 22 : & ' XZ = ', chi_tensor_tmp(1, 3)
782 11 : WRITE (unit_atoms, '(3(A,f10.4))') ' YX = ', chi_tensor_tmp(2, 1),&
783 11 : & ' YY = ', chi_tensor_tmp(2, 2),&
784 22 : & ' YZ = ', chi_tensor_tmp(2, 3)
785 11 : WRITE (unit_atoms, '(3(A,f10.4))') ' ZX = ', chi_tensor_tmp(3, 1),&
786 11 : & ' ZY = ', chi_tensor_tmp(3, 2),&
787 22 : & ' ZZ = ', chi_tensor_tmp(3, 3)
788 11 : IF (gapw) THEN
789 11 : WRITE (unit_atoms, '(T1,A)') " CHI from LOCAL J in 10^-30 J/T^2 units"
790 11 : WRITE (unit_atoms, '(3(A,f10.4))') ' XX = ', chi_tensor_loc_tmp(1, 1),&
791 11 : & ' XY = ', chi_tensor_loc_tmp(1, 2),&
792 22 : & ' XZ = ', chi_tensor_loc_tmp(1, 3)
793 11 : WRITE (unit_atoms, '(3(A,f10.4))') ' YX = ', chi_tensor_loc_tmp(2, 1),&
794 11 : & ' YY = ', chi_tensor_loc_tmp(2, 2),&
795 22 : & ' YZ = ', chi_tensor_loc_tmp(2, 3)
796 11 : WRITE (unit_atoms, '(3(A,f10.4))') ' ZX = ', chi_tensor_loc_tmp(3, 1),&
797 11 : & ' ZY = ', chi_tensor_loc_tmp(3, 2),&
798 22 : & ' ZZ = ', chi_tensor_loc_tmp(3, 3)
799 : END IF
800 11 : WRITE (unit_atoms, '(T1,A)') " Total CHI in 10^-30 J/T^2 units"
801 11 : WRITE (unit_atoms, '(3(A,f10.4))') ' XX = ', chi_sym_tot(1, 1),&
802 11 : & ' XY = ', chi_sym_tot(1, 2),&
803 22 : & ' XZ = ', chi_sym_tot(1, 3)
804 11 : WRITE (unit_atoms, '(3(A,f10.4))') ' YX = ', chi_sym_tot(2, 1),&
805 11 : & ' YY = ', chi_sym_tot(2, 2),&
806 22 : & ' YZ = ', chi_sym_tot(2, 3)
807 11 : WRITE (unit_atoms, '(3(A,f10.4))') ' ZX = ', chi_sym_tot(3, 1),&
808 11 : & ' ZY = ', chi_sym_tot(3, 2),&
809 22 : & ' ZZ = ', chi_sym_tot(3, 3)
810 143 : chi_sym_tot(:, :) = chi_sym_tot(:, :)*nmr_env%chi_SI2ppmcgs
811 11 : WRITE (unit_atoms, '(T1,A)') " Total CHI in ppm-cgs units"
812 11 : WRITE (unit_atoms, '(3(A,f10.4))') ' XX = ', chi_sym_tot(1, 1),&
813 11 : & ' XY = ', chi_sym_tot(1, 2),&
814 22 : & ' XZ = ', chi_sym_tot(1, 3)
815 11 : WRITE (unit_atoms, '(3(A,f10.4))') ' YX = ', chi_sym_tot(2, 1),&
816 11 : & ' YY = ', chi_sym_tot(2, 2),&
817 22 : & ' YZ = ', chi_sym_tot(2, 3)
818 11 : WRITE (unit_atoms, '(3(A,f10.4))') ' ZX = ', chi_sym_tot(3, 1),&
819 11 : & ' ZY = ', chi_sym_tot(3, 2),&
820 22 : & ' ZZ = ', chi_sym_tot(3, 3)
821 : WRITE (unit_atoms, '(/T1,3(A,f10.4))') &
822 11 : ' PV1=', nmr_env%chi_SI2ppmcgs*eig(1), &
823 11 : ' PV2=', nmr_env%chi_SI2ppmcgs*eig(2), &
824 22 : ' PV3=', nmr_env%chi_SI2ppmcgs*eig(3)
825 : WRITE (unit_atoms, '(T1,A,F10.4,10X,A,F10.4)') &
826 11 : ' ISO=', nmr_env%chi_SI2ppmcgs*chi_iso, &
827 22 : 'ANISO=', nmr_env%chi_SI2ppmcgs*chi_aniso
828 : END IF
829 :
830 : CALL cp_print_key_finished_output(unit_atoms, logger, nmr_section,&
831 22 : & "PRINT%CHI_TENSOR")
832 : END IF ! print chi
833 : !
834 : ! Add the chi part to the shifts
835 6504 : cs_tot(:, :, :) = 0.0_dp
836 648 : DO ir = 1, nat_print
837 488 : iatom = cs_atom_list(ir)
838 1952 : rpos(1:3) = particle_set(iatom)%r(1:3)
839 488 : atom_kind => particle_set(iatom)%atomic_kind
840 488 : CALL get_atomic_kind(atom_kind, name=name, element_symbol=element_symbol)
841 12200 : cs_tot(:, :, ir) = chi_tensor_tmp(:, :)*nmr_env%chi_SI2shiftppm + cs(:, :, iatom)
842 7368 : IF (gapw) cs_tot(:, :, ir) = cs_tot(:, :, ir) + cs_loc(:, :, iatom)
843 : END DO ! ir
844 160 : IF (output_unit > 0) THEN
845 80 : WRITE (output_unit, '(T2,A,E14.6)') 'CheckSum Shifts =', &
846 160 : SQRT(DDOT(9*SIZE(cs_tot, 3), cs_tot(1, 1, 1), 1, cs_tot(1, 1, 1), 1))
847 : END IF
848 : !
849 : ! print shifts
850 160 : IF (BTEST(cp_print_key_should_output(logger%iter_info, nmr_section, &
851 : "PRINT%SHIELDING_TENSOR"), cp_p_file)) THEN
852 :
853 : unit_atoms = cp_print_key_unit_nr(logger, nmr_section, "PRINT%SHIELDING_TENSOR", &
854 : extension=".data", middle_name="SHIFT", &
855 160 : log_filename=.FALSE.)
856 :
857 160 : nat_print = SIZE(cs_atom_list, 1)
858 160 : IF (unit_atoms > 0) THEN
859 80 : WRITE (title, '(A,1X,I5)') "Shielding atom at atomic positions. # tensors printed ", nat_print
860 80 : WRITE (unit_atoms, '(T2,A)') title
861 324 : DO ir = 1, nat_print
862 244 : iatom = cs_atom_list(ir)
863 976 : rpos(1:3) = particle_set(iatom)%r(1:3)
864 244 : atom_kind => particle_set(iatom)%atomic_kind
865 244 : CALL get_atomic_kind(atom_kind, name=name, element_symbol=element_symbol)
866 3172 : shift_sym_tot(:, :) = 0.5_dp*(cs_tot(:, :, ir) + TRANSPOSE(cs_tot(:, :, ir)))
867 244 : CALL diamat_all(shift_sym_tot, eig)
868 244 : shift_iso = (eig(1) + eig(2) + eig(3))/3.0_dp
869 244 : shift_aniso = eig(3) - (eig(2) + eig(1))/2.0_dp
870 : !
871 244 : WRITE (unit_atoms, '(T2,I5,A,2X,A2,2X,3f15.6)') iatom, TRIM(name), element_symbol, rpos(1:3)
872 : !
873 244 : IF (gapw) THEN
874 140 : WRITE (unit_atoms, '(T1,A)') " SIGMA from SOFT J"
875 140 : WRITE (unit_atoms, '(3(A,f10.4))') ' XX = ', cs(1, 1, iatom),&
876 140 : & ' XY = ', cs(1, 2, iatom),&
877 280 : & ' XZ = ', cs(1, 3, iatom)
878 140 : WRITE (unit_atoms, '(3(A,f10.4))') ' YX = ', cs(2, 1, iatom),&
879 140 : & ' YY = ', cs(2, 2, iatom),&
880 280 : & ' YZ = ', cs(2, 3, iatom)
881 140 : WRITE (unit_atoms, '(3(A,f10.4))') ' ZX = ', cs(3, 1, iatom),&
882 140 : & ' ZY = ', cs(3, 2, iatom),&
883 280 : & ' ZZ = ', cs(3, 3, iatom)
884 140 : WRITE (unit_atoms, '(T1,A)') " SIGMA from LOCAL J"
885 140 : WRITE (unit_atoms, '(3(A,f10.4))') ' XX = ', cs_loc(1, 1, iatom),&
886 140 : & ' XY = ', cs_loc(1, 2, iatom),&
887 280 : & ' XZ = ', cs_loc(1, 3, iatom)
888 140 : WRITE (unit_atoms, '(3(A,f10.4))') ' YX = ', cs_loc(2, 1, iatom),&
889 140 : & ' YY = ', cs_loc(2, 2, iatom),&
890 280 : & ' YZ = ', cs_loc(2, 3, iatom)
891 140 : WRITE (unit_atoms, '(3(A,f10.4))') ' ZX = ', cs_loc(3, 1, iatom),&
892 140 : & ' ZY = ', cs_loc(3, 2, iatom),&
893 280 : & ' ZZ = ', cs_loc(3, 3, iatom)
894 : END IF
895 244 : WRITE (unit_atoms, '(T1,A)') " SIGMA TOTAL"
896 244 : WRITE (unit_atoms, '(3(A,f10.4))') ' XX = ', cs_tot(1, 1, ir),&
897 244 : & ' XY = ', cs_tot(1, 2, ir),&
898 488 : & ' XZ = ', cs_tot(1, 3, ir)
899 244 : WRITE (unit_atoms, '(3(A,f10.4))') ' YX = ', cs_tot(2, 1, ir),&
900 244 : & ' YY = ', cs_tot(2, 2, ir),&
901 488 : & ' YZ = ', cs_tot(2, 3, ir)
902 244 : WRITE (unit_atoms, '(3(A,f10.4))') ' ZX = ', cs_tot(3, 1, ir),&
903 244 : & ' ZY = ', cs_tot(3, 2, ir),&
904 488 : & ' ZZ = ', cs_tot(3, 3, ir)
905 244 : WRITE (unit_atoms, '(T1,2(A,f12.4))') ' ISOTROPY = ', shift_iso,&
906 568 : & ' ANISOTROPY = ', shift_aniso
907 : END DO ! ir
908 : END IF
909 : CALL cp_print_key_finished_output(unit_atoms, logger, nmr_section,&
910 160 : & "PRINT%SHIELDING_TENSOR")
911 :
912 160 : IF (do_nics) THEN
913 : !
914 : ! Add the chi part to the nics
915 318 : cs_nics_tot(:, :, :) = 0.0_dp
916 30 : DO ir = 1, n_nics
917 600 : cs_nics_tot(:, :, ir) = chi_tensor_tmp(:, :)*nmr_env%chi_SI2shiftppm + cs_nics(:, :, ir)
918 174 : IF (gapw) cs_nics_tot(:, :, ir) = cs_nics_tot(:, :, ir) + cs_nics_loc(:, :, ir)
919 : END DO ! ir
920 6 : IF (output_unit > 0) THEN
921 3 : WRITE (output_unit, '(T2,A,E14.6)') 'CheckSum NICS =', &
922 6 : SQRT(DDOT(9*SIZE(cs_nics_tot, 3), cs_nics_tot(1, 1, 1), 1, cs_nics_tot(1, 1, 1), 1))
923 : END IF
924 : !
925 : unit_nics = cp_print_key_unit_nr(logger, nmr_section, "PRINT%SHIELDING_TENSOR", &
926 : extension=".data", middle_name="NICS", &
927 6 : log_filename=.FALSE.)
928 6 : IF (unit_nics > 0) THEN
929 3 : WRITE (title, '(A,1X,I5)') "Shielding at nics positions. # tensors printed ", n_nics
930 3 : WRITE (unit_nics, '(T2,A)') title
931 15 : DO ir = 1, n_nics
932 156 : shift_sym_tot(:, :) = 0.5_dp*(cs_nics_tot(:, :, ir) + TRANSPOSE(cs_nics_tot(:, :, ir)))
933 12 : CALL diamat_all(shift_sym_tot, eig)
934 12 : shift_iso = (eig(1) + eig(2) + eig(3))/3.0_dp
935 12 : shift_aniso = eig(3) - (eig(2) + eig(1))/2.0_dp
936 : !
937 12 : WRITE (unit_nics, '(T2,I5,2X,3f15.6)') ir, r_nics(1:3, ir)
938 : !
939 12 : IF (gapw) THEN
940 3 : WRITE (unit_nics, '(T1,A)') " SIGMA from SOFT J"
941 3 : WRITE (unit_nics, '(3(A,f10.4))') ' XX = ', cs_nics(1, 1, ir),&
942 3 : & ' XY = ', cs_nics(1, 2, ir),&
943 6 : & ' XZ = ', cs_nics(1, 3, ir)
944 3 : WRITE (unit_nics, '(3(A,f10.4))') ' YX = ', cs_nics(2, 1, ir),&
945 3 : & ' YY = ', cs_nics(2, 2, ir),&
946 6 : & ' YZ = ', cs_nics(2, 3, ir)
947 3 : WRITE (unit_nics, '(3(A,f10.4))') ' ZX = ', cs_nics(3, 1, ir),&
948 3 : & ' ZY = ', cs_nics(3, 2, ir),&
949 6 : & ' ZZ = ', cs_nics(3, 3, ir)
950 3 : WRITE (unit_nics, '(T1,A)') " SIGMA from LOCAL J"
951 3 : WRITE (unit_nics, '(3(A,f10.4))') ' XX = ', cs_nics_loc(1, 1, ir),&
952 3 : & ' XY = ', cs_nics_loc(1, 2, ir),&
953 6 : & ' XZ = ', cs_nics_loc(1, 3, ir)
954 3 : WRITE (unit_nics, '(3(A,f10.4))') ' YX = ', cs_nics_loc(2, 1, ir),&
955 3 : & ' YY = ', cs_nics_loc(2, 2, ir),&
956 6 : & ' YZ = ', cs_nics_loc(2, 3, ir)
957 3 : WRITE (unit_nics, '(3(A,f10.4))') ' ZX = ', cs_nics_loc(3, 1, ir),&
958 3 : & ' ZY = ', cs_nics_loc(3, 2, ir),&
959 6 : & ' ZZ = ', cs_nics_loc(3, 3, ir)
960 : END IF
961 12 : WRITE (unit_nics, '(T1,A)') " SIGMA TOTAL"
962 12 : WRITE (unit_nics, '(3(A,f10.4))') ' XX = ', cs_nics_tot(1, 1, ir),&
963 12 : & ' XY = ', cs_nics_tot(1, 2, ir),&
964 24 : & ' XZ = ', cs_nics_tot(1, 3, ir)
965 12 : WRITE (unit_nics, '(3(A,f10.4))') ' YX = ', cs_nics_tot(2, 1, ir),&
966 12 : & ' YY = ', cs_nics_tot(2, 2, ir),&
967 24 : & ' YZ = ', cs_nics_tot(2, 3, ir)
968 12 : WRITE (unit_nics, '(3(A,f10.4))') ' ZX = ', cs_nics_tot(3, 1, ir),&
969 12 : & ' ZY = ', cs_nics_tot(3, 2, ir),&
970 24 : & ' ZZ = ', cs_nics_tot(3, 3, ir)
971 12 : WRITE (unit_nics, '(T1,2(A,f12.4))') ' ISOTROPY = ', shift_iso,&
972 27 : & ' ANISOTROPY = ', shift_aniso
973 : END DO
974 : END IF
975 : CALL cp_print_key_finished_output(unit_nics, logger, nmr_section,&
976 6 : & "PRINT%SHIELDING_TENSOR")
977 : END IF
978 : END IF ! print shift
979 : !
980 : ! clean up
981 160 : DEALLOCATE (cs_tot)
982 160 : IF (do_nics) THEN
983 6 : DEALLOCATE (cs_nics_tot)
984 : END IF
985 : !
986 : !100 FORMAT(A,1X,I5)
987 : !101 FORMAT(T2,A)
988 : !102 FORMAT(T2,I5,A,2X,A2,2X,3f15.6)
989 : !103 FORMAT(T2,I5,2X,3f15.6)
990 : !104 FORMAT(T1,A)
991 : !105 FORMAT(3(A,f10.4))
992 : !106 FORMAT(T1,2(A,f12.4))
993 320 : END SUBROUTINE nmr_shift_print
994 :
995 : END MODULE qs_linres_nmr_shift
996 :
|
