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