(git:0de0cc2)
qs_linres_nmr_shift.F
Go to the documentation of this file.
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 ! **************************************************************************************************
16  USE atomic_kind_types, ONLY: atomic_kind_type,&
18  USE cell_types, ONLY: cell_type,&
19  pbc
20  USE cp_control_types, ONLY: dft_control_type
23  cp_logger_type
24  USE cp_output_handling, ONLY: cp_p_file,&
29  section_vals_type,&
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
47  find_coeffs,&
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
63  USE qs_linres_op, ONLY: fac_vecp,&
64  set_vecp,&
66  USE qs_linres_types, ONLY: current_env_type,&
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  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  REAL(dp), DIMENSION(:, :, :), POINTER :: chemical_shift, chemical_shift_loc, &
111  chemical_shift_nics, &
112  chemical_shift_nics_loc
113  TYPE(cell_type), POINTER :: cell
114  TYPE(dft_control_type), POINTER :: dft_control
115  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
116  TYPE(pw_c1d_gs_type) :: pw_gspace_work
117  TYPE(pw_c1d_gs_type), ALLOCATABLE, DIMENSION(:, :) :: shift_pw_gspace
118  TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: jrho1_g
119  TYPE(pw_env_type), POINTER :: pw_env
120  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  CALL timeset(routinen, handle)
127 
128  NULLIFY (chemical_shift, chemical_shift_loc, chemical_shift_nics, &
129  chemical_shift_nics_loc, cell, dft_control, pw_env, &
130  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  particle_set=particle_set)
134 
135  gapw = dft_control%qs_control%gapw
136  natom = SIZE(particle_set, 1)
137  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  interpolate_shift=interpolate_shift)
144 
145  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  auxbas_pw_pool=auxbas_pw_pool, pw_pools=pw_pools)
148  !
149  !
150  nmr_section => section_vals_get_subs_vals(qs_env%input, &
151  & "PROPERTIES%LINRES%NMR")
152  !
153  ! Initialize
154  ! Allocate grids for the calculation of jrho and the shift
155  ALLOCATE (shift_pw_gspace(3, nspins))
156  DO ispin = 1, nspins
157  DO idir = 1, 3
158  CALL auxbas_pw_pool%create_pw(shift_pw_gspace(idir, ispin))
159  CALL pw_zero(shift_pw_gspace(idir, ispin))
160  END DO
161  END DO
162  !
163  !
164  CALL set_vecp(ib, iib, iiib)
165  !
166  CALL auxbas_pw_pool%create_pw(pw_gspace_work)
167  CALL pw_zero(pw_gspace_work)
168  DO ispin = 1, nspins
169  !
170  DO idir = 1, 3
171  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  DO idir2 = 1, 3
175  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  pw_gspace_work, idir2, 0.0_dp)
179  !
180  ! scale and add to the correct component of the shift column
181  CALL set_vecp_rev(idir, idir2, idir3)
182  scale_fac = fac_vecp(idir3, idir2, idir)
183  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  CALL auxbas_pw_pool%give_back_pw(pw_gspace_work)
191  !
192  ! compute shildings
193  IF (interpolate_shift) THEN
194  CALL auxbas_pw_pool%create_pw(shift_pw_rspace)
195  DO ispin = 1, nspins
196  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  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  ib, idir, nmr_section)
203  END DO
204  END DO
205  CALL auxbas_pw_pool%give_back_pw(shift_pw_rspace)
206  ELSE
207  DO ispin = 1, nspins
208  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  shift_pw_gspace(idir, ispin), ib, idir)
212  END DO
213  END DO
214  END IF
215  !
216  IF (gapw) THEN
217  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  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  DO ispin = 1, nspins
227  DO idir = 1, 3
228  CALL auxbas_pw_pool%give_back_pw(shift_pw_gspace(idir, ispin))
229  END DO
230  END DO
231  DEALLOCATE (shift_pw_gspace)
232  !
233  ! Finalize
234  CALL timestop(handle)
235  !
236  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  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  INTEGER, ALLOCATABLE, DIMENSION(:) :: list_j, list_nics_j
259  INTEGER, DIMENSION(2) :: bo
260  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  REAL(dp), ALLOCATABLE, DIMENSION(:) :: j_grid
266  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: cs_loc_tmp, cs_nics_loc_tmp, dist_ij, &
267  dist_nics_ij, r_grid
268  REAL(dp), DIMENSION(:, :), POINTER :: jrho_h_grid, jrho_s_grid, r_nics
269  REAL(dp), DIMENSION(:, :, :), POINTER :: chemical_shift_loc, &
270  chemical_shift_nics_loc
271  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  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  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
281  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
282 
283  CALL timeset(routinen, handle)
284  !
285  NULLIFY (atomic_kind_set, qs_kind_set, cell, dft_control, para_env, particle_set, &
286  chemical_shift_loc, chemical_shift_nics_loc, jrho1_atom_set, &
287  jrho1_atom, r_nics, jrho_h_grid, jrho_s_grid, &
288  atom_list, grid_atom, harmonics, logger)
289  !
290  logger => cp_get_default_logger()
291  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  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  do_nics=do_nics)
308 
309  CALL get_current_env(current_env=current_env, &
310  jrho1_atom_set=jrho1_atom_set)
311  !
312  nkind = SIZE(atomic_kind_set, 1)
313  natom_tot = SIZE(particle_set, 1)
314  nspins = dft_control%nspins
315  itegrated_jrho = 0.0_dp
316  !
317  idir2_1 = modulo(idir, 3) + 1
318  idir3_1 = modulo(idir + 1, 3) + 1
319  scale_fac_1 = fac_vecp(idir3_1, idir2_1, idir)
320  !
321  ALLOCATE (cs_loc_tmp(3, natom_tot), list_j(natom_tot), &
322  dist_ij(3, natom_tot))
323  cs_loc_tmp = 0.0_dp
324  IF (do_nics) THEN
325  ALLOCATE (cs_nics_loc_tmp(3, n_nics), list_nics_j(n_nics), &
326  dist_nics_ij(3, n_nics))
327  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  DO ikind = 1, nkind
333 
334  NULLIFY (atom_list, grid_atom, harmonics)
335  CALL get_atomic_kind(atomic_kind_set(ikind), &
336  atom_list=atom_list, &
337  natom=natom)
338 
339  CALL get_qs_kind(qs_kind_set(ikind), &
340  paw_atom=paw_atom, &
341  harmonics=harmonics, &
342  grid_atom=grid_atom)
343  !
344  na = grid_atom%ng_sphere
345  nr = grid_atom%nr
346  nra = nr*na
347  ALLOCATE (r_grid(3, nra), j_grid(nra))
348  ira = 1
349  DO ia = 1, na
350  DO ir = 1, nr
351  r_grid(:, ira) = grid_atom%rad(ir)*harmonics%a(:, ia)
352  ira = ira + 1
353  END DO
354  END DO
355  !
356  ! Quick cycle if needed
357  IF (paw_atom) THEN
358  !
359  ! Distribute the atoms of this kind
360  bo = get_limit(natom, para_env%num_pe, para_env%mepos)
361  !
362  DO iat = bo(1), bo(2)
363  iatom = atom_list(iat)
364  !
365  ! find all the atoms within the radius
366  natom_local = 0
367  DO jatom = 1, natom_tot
368  rij(:) = pbc(particle_set(iatom)%r, particle_set(jatom)%r, cell)
369  dist = sqrt(rij(1)*rij(1) + rij(2)*rij(2) + rij(3)*rij(3))
370  IF (dist .LE. shift_gapw_radius) THEN
371  natom_local = natom_local + 1
372  list_j(natom_local) = jatom
373  dist_ij(:, natom_local) = rij(:)
374  END IF
375  END DO
376  !
377  ! ... also for nics
378  IF (do_nics) THEN
379  nnics_local = 0
380  DO jatom = 1, n_nics
381  r_jatom(:) = r_nics(:, jatom)
382  rij(:) = pbc(particle_set(iatom)%r, r_jatom, cell)
383  dist = sqrt(rij(1)*rij(1) + rij(2)*rij(2) + rij(3)*rij(3))
384  IF (dist .LE. shift_gapw_radius) THEN
385  nnics_local = nnics_local + 1
386  list_nics_j(nnics_local) = jatom
387  dist_nics_ij(:, nnics_local) = rij(:)
388  END IF
389  END DO
390  END IF
391  !
392  NULLIFY (jrho1_atom, jrho_h_grid, jrho_s_grid)
393  jrho1_atom => jrho1_atom_set(iatom)
394  !
395  DO ispin = 1, nspins
396  jrho_h_grid => jrho1_atom%jrho_vec_rad_h(idir, ispin)%r_coef
397  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  ira = 1
403  DO ia = 1, na
404  DO ir = 1, nr
405  j_grid(ira) = (jrho_h_grid(ir, ia) - jrho_s_grid(ir, ia))*grid_atom%weight(ia, ir)
406  itegrated_jrho = itegrated_jrho + j_grid(ira)
407  ira = ira + 1
408  END DO
409  END DO
410  !
411  DO j = 1, natom_local
412  jatom = list_j(j)
413  rij(:) = dist_ij(:, j)
414  !
415  s_1 = 0.0_dp
416  s_2 = 0.0_dp
417  DO ira = 1, nra
418  !
419  rdiff(:) = rij(:) - r_grid(:, ira)
420  ddiff = sqrt(rdiff(1)*rdiff(1) + rdiff(2)*rdiff(2) + rdiff(3)*rdiff(3))
421  IF (ddiff .GT. 1.0e-12_dp) THEN
422  dum = scale_fac_1*j_grid(ira)/(ddiff*ddiff*ddiff)
423  s_1 = s_1 + rdiff(idir2_1)*dum
424  s_2 = s_2 + rdiff(idir3_1)*dum
425  END IF ! ddiff
426  END DO ! ira
427  cs_loc_tmp(idir3_1, jatom) = cs_loc_tmp(idir3_1, jatom) + s_1
428  cs_loc_tmp(idir2_1, jatom) = cs_loc_tmp(idir2_1, jatom) - s_2
429  END DO ! j
430  !
431  IF (do_nics) THEN
432  DO j = 1, nnics_local
433  jatom = list_nics_j(j)
434  rij(:) = dist_nics_ij(:, j)
435  !
436  s_1 = 0.0_dp
437  s_2 = 0.0_dp
438  DO ira = 1, nra
439  !
440  rdiff(:) = rij(:) - r_grid(:, ira)
441  ddiff = sqrt(rdiff(1)*rdiff(1) + rdiff(2)*rdiff(2) + rdiff(3)*rdiff(3))
442  IF (ddiff .GT. 1.0e-12_dp) THEN
443  dum = scale_fac_1*j_grid(ira)/(ddiff*ddiff*ddiff)
444  s_1 = s_1 + rdiff(idir2_1)*dum
445  s_2 = s_2 + rdiff(idir3_1)*dum
446  END IF ! ddiff
447  END DO ! ira
448  cs_nics_loc_tmp(idir3_1, jatom) = cs_nics_loc_tmp(idir3_1, jatom) + s_1
449  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  DEALLOCATE (r_grid, j_grid)
456  END DO ! ikind
457  !
458  !
459  CALL para_env%sum(itegrated_jrho)
460  IF (output_unit > 0) THEN
461  WRITE (output_unit, '(T2,A,E24.16)') 'Integrated local j_'&
462  &//achar(idir + 119)//achar(ib + 119)//'(r)=', itegrated_jrho
463  END IF
464  !
465  CALL para_env%sum(cs_loc_tmp)
466  chemical_shift_loc(:, ib, :) = chemical_shift_loc(:, ib, :) &
467  & - nmr_env%shift_factor_gapw*cs_loc_tmp(:, :)/2.0_dp
468  !
469  DEALLOCATE (cs_loc_tmp, list_j, dist_ij)
470  !
471  IF (do_nics) THEN
472  CALL para_env%sum(cs_nics_loc_tmp)
473  chemical_shift_nics_loc(:, ib, :) = chemical_shift_nics_loc(:, ib, :) &
474  & - nmr_env%shift_factor_gapw*cs_nics_loc_tmp(:, :)/2.0_dp
475  !
476  DEALLOCATE (cs_nics_loc_tmp, list_nics_j, dist_nics_ij)
477  END IF
478  !
479  CALL timestop(handle)
480  !
481  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  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  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  REAL(dp), DIMENSION(:, :), POINTER :: r_nics
516  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  CALL timeset(routinen, handle)
523  !
524  NULLIFY (interp_section, auxbas_pw_pool)
525  NULLIFY (cs_atom_list, chemical_shift, chemical_shift_nics, r_nics)
526 
527  interp_section => section_vals_get_subs_vals(nmr_section, &
528  "INTERPOLATOR")
529  CALL section_vals_val_get(interp_section, "aint_precond", &
530  i_val=aint_precond)
531  CALL section_vals_val_get(interp_section, "precond", i_val=precond_kind)
532  CALL section_vals_val_get(interp_section, "max_iter", i_val=max_iter)
533  CALL section_vals_val_get(interp_section, "eps_r", r_val=eps_r)
534  CALL section_vals_val_get(interp_section, "eps_x", r_val=eps_x)
535 
536  ! calculate spline coefficients
537  CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
538  CALL auxbas_pw_pool%create_pw(shiftspl)
539 
540  CALL pw_spline_precond_create(precond, precond_kind=aint_precond, &
541  pool=auxbas_pw_pool, pbc=.true., transpose=.false.)
542  CALL pw_spline_do_precond(precond, shift_pw_rspace, shiftspl)
543  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  eps_r=eps_r, eps_x=eps_x, max_iter=max_iter)
547  cpassert(success)
548  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  do_nics=do_nics)
555 
556  IF (ASSOCIATED(cs_atom_list)) THEN
557  natom = SIZE(cs_atom_list, 1)
558  ELSE
559  natom = -1
560  END IF
561 
562  DO iat = 1, natom
563  iatom = cs_atom_list(iat)
564  r_iatom = pbc(particle_set(iatom)%r, cell)
565  shift_val = eval_interp_spl3_pbc(r_iatom, shiftspl)
566  chemical_shift(idir, i_b, iatom) = chemical_shift(idir, i_b, iatom) + &
567  nmr_env%shift_factor*twopi**2*shift_val
568  END DO
569 
570  IF (do_nics) THEN
571  DO iatom = 1, n_nics
572  ra(:) = r_nics(:, iatom)
573  r_iatom = pbc(ra, cell)
574  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  nmr_env%shift_factor*twopi**2*shift_val
577  END DO
578  END IF
579 
580  CALL auxbas_pw_pool%give_back_pw(shiftspl)
581 
582  !
583  CALL timestop(handle)
584  !
585  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  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  INTEGER, DIMENSION(:), POINTER :: cs_atom_list
609  LOGICAL :: do_nics
610  REAL(dp) :: r_iatom(3), ra(3)
611  REAL(dp), DIMENSION(:, :), POINTER :: r_nics
612  REAL(dp), DIMENSION(:, :, :), POINTER :: chemical_shift, chemical_shift_nics
613 
614  CALL timeset(routinen, handle)
615  !
616  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  n_nics=n_nics, r_nics=r_nics, do_nics=do_nics)
622  !
623  IF (ASSOCIATED(cs_atom_list)) THEN
624  natom = SIZE(cs_atom_list, 1)
625  ELSE
626  natom = -1
627  END IF
628  !
629  ! compute the chemical shift
630  DO iat = 1, natom
631  iatom = cs_atom_list(iat)
632  r_iatom = pbc(particle_set(iatom)%r, cell)
633  CALL gsumr(r_iatom, shift_pw_gspace, cplx)
634  chemical_shift(idir, i_b, iatom) = chemical_shift(idir, i_b, iatom) + &
635  nmr_env%shift_factor*twopi**2*real(cplx, dp)
636  END DO
637  !
638  ! compute nics
639  IF (do_nics) THEN
640  DO iat = 1, n_nics
641  ra = pbc(r_nics(:, iat), cell)
642  CALL gsumr(ra, shift_pw_gspace, cplx)
643  chemical_shift_nics(idir, i_b, iat) = chemical_shift_nics(idir, i_b, iat) + &
644  nmr_env%shift_factor*twopi**2*real(cplx, dp)
645  END DO
646  END IF
647 
648  CALL timestop(handle)
649 
650  END SUBROUTINE gsum_shift_pwgrid
651 
652 ! **************************************************************************************************
653 !> \brief ...
654 !> \param r ...
655 !> \param pw ...
656 !> \param cplx ...
657 ! **************************************************************************************************
658  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  grid => pw%pw_grid
668  cplx = cmplx(0.0_dp, 0.0_dp, kind=dp)
669  DO ig = grid%first_gne0, grid%ngpts_cut_local
670  rg = (grid%g(1, ig)*r(1) + grid%g(2, ig)*r(2) + grid%g(3, ig)*r(3))*gaussi
671  cplx = cplx + pw%array(ig)*exp(rg)
672  END DO
673  IF (grid%have_g0) cplx = cplx + pw%array(1)
674  CALL grid%para%group%sum(cplx)
675  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  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  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  REAL(dp), DIMENSION(:, :), POINTER :: r_nics
702  REAL(dp), DIMENSION(:, :, :), POINTER :: cs, cs_loc, cs_nics, cs_nics_loc, &
703  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  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
709  TYPE(section_vals_type), POINTER :: nmr_section
710 
711  NULLIFY (cs, cs_nics, r_nics, cs_loc, cs_nics_loc, logger, particle_set, atom_kind, dft_control)
712 
713  logger => cp_get_default_logger()
714  output_unit = cp_logger_get_default_io_unit(logger)
715 
716  nmr_section => section_vals_get_subs_vals(qs_env%input, &
717  "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  do_nics=do_nics)
728  !
729  CALL get_current_env(current_env=current_env, &
730  chi_tensor=chi_tensor, &
731  chi_tensor_loc=chi_tensor_loc)
732  !
733  ! multiply by the appropriate factor
734  chi_tensor_tmp(:, :) = 0.0_dp
735  chi_tensor_loc_tmp(:, :) = 0.0_dp
736  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  particle_set=particle_set)
742 
743  natom = SIZE(particle_set, 1)
744  gapw = dft_control%qs_control%gapw
745  nat_print = SIZE(cs_atom_list, 1)
746 
747  ALLOCATE (cs_tot(3, 3, nat_print))
748  IF (do_nics) THEN
749  ALLOCATE (cs_nics_tot(3, 3, n_nics))
750  END IF
751  ! Finalize Chi calculation
752  ! Symmetrize
753  chi_sym_tot(:, :) = (chi_tensor_tmp(:, :) + transpose(chi_tensor_tmp(:, :)))/2.0_dp
754  IF (gapw) THEN
755  chi_sym_tot(:, :) = chi_sym_tot(:, :) &
756  & + (chi_tensor_loc_tmp(:, :) + transpose(chi_tensor_loc_tmp(:, :)))/2.0_dp
757  END IF
758  chi_tmp(:, :) = chi_sym_tot(:, :)
759  CALL diamat_all(chi_tmp, eig)
760  chi_iso = (eig(1) + eig(2) + eig(3))/3.0_dp
761  chi_aniso = eig(3) - (eig(2) + eig(1))/2.0_dp
762  !
763  IF (output_unit > 0) THEN
764  WRITE (output_unit, '(T2,A,E14.6)') 'CheckSum Chi =', &
765  sqrt(ddot(9, chi_tensor_tmp(1, 1), 1, chi_tensor_tmp(1, 1), 1))
766  END IF
767  !
768  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  extension=".data", middle_name="CHI", log_filename=.false.)
773 
774  WRITE (title, '(A)') "Magnetic Susceptibility Tensor "
775  IF (unit_atoms > 0) THEN
776  WRITE (unit_atoms, '(T2,A)') title
777  WRITE (unit_atoms, '(T1,A)') " CHI from SOFT J in 10^-30 J/T^2 units"
778  WRITE (unit_atoms, '(3(A,f10.4))') ' XX = ', chi_tensor_tmp(1, 1),&
779  & ' XY = ', chi_tensor_tmp(1, 2),&
780  & ' XZ = ', chi_tensor_tmp(1, 3)
781  WRITE (unit_atoms, '(3(A,f10.4))') ' YX = ', chi_tensor_tmp(2, 1),&
782  & ' YY = ', chi_tensor_tmp(2, 2),&
783  & ' YZ = ', chi_tensor_tmp(2, 3)
784  WRITE (unit_atoms, '(3(A,f10.4))') ' ZX = ', chi_tensor_tmp(3, 1),&
785  & ' ZY = ', chi_tensor_tmp(3, 2),&
786  & ' ZZ = ', chi_tensor_tmp(3, 3)
787  IF (gapw) THEN
788  WRITE (unit_atoms, '(T1,A)') " CHI from LOCAL J in 10^-30 J/T^2 units"
789  WRITE (unit_atoms, '(3(A,f10.4))') ' XX = ', chi_tensor_loc_tmp(1, 1),&
790  & ' XY = ', chi_tensor_loc_tmp(1, 2),&
791  & ' XZ = ', chi_tensor_loc_tmp(1, 3)
792  WRITE (unit_atoms, '(3(A,f10.4))') ' YX = ', chi_tensor_loc_tmp(2, 1),&
793  & ' YY = ', chi_tensor_loc_tmp(2, 2),&
794  & ' YZ = ', chi_tensor_loc_tmp(2, 3)
795  WRITE (unit_atoms, '(3(A,f10.4))') ' ZX = ', chi_tensor_loc_tmp(3, 1),&
796  & ' ZY = ', chi_tensor_loc_tmp(3, 2),&
797  & ' ZZ = ', chi_tensor_loc_tmp(3, 3)
798  END IF
799  WRITE (unit_atoms, '(T1,A)') " Total CHI in 10^-30 J/T^2 units"
800  WRITE (unit_atoms, '(3(A,f10.4))') ' XX = ', chi_sym_tot(1, 1),&
801  & ' XY = ', chi_sym_tot(1, 2),&
802  & ' XZ = ', chi_sym_tot(1, 3)
803  WRITE (unit_atoms, '(3(A,f10.4))') ' YX = ', chi_sym_tot(2, 1),&
804  & ' YY = ', chi_sym_tot(2, 2),&
805  & ' YZ = ', chi_sym_tot(2, 3)
806  WRITE (unit_atoms, '(3(A,f10.4))') ' ZX = ', chi_sym_tot(3, 1),&
807  & ' ZY = ', chi_sym_tot(3, 2),&
808  & ' ZZ = ', chi_sym_tot(3, 3)
809  chi_sym_tot(:, :) = chi_sym_tot(:, :)*nmr_env%chi_SI2ppmcgs
810  WRITE (unit_atoms, '(T1,A)') " Total CHI in ppm-cgs units"
811  WRITE (unit_atoms, '(3(A,f10.4))') ' XX = ', chi_sym_tot(1, 1),&
812  & ' XY = ', chi_sym_tot(1, 2),&
813  & ' XZ = ', chi_sym_tot(1, 3)
814  WRITE (unit_atoms, '(3(A,f10.4))') ' YX = ', chi_sym_tot(2, 1),&
815  & ' YY = ', chi_sym_tot(2, 2),&
816  & ' YZ = ', chi_sym_tot(2, 3)
817  WRITE (unit_atoms, '(3(A,f10.4))') ' ZX = ', chi_sym_tot(3, 1),&
818  & ' ZY = ', chi_sym_tot(3, 2),&
819  & ' ZZ = ', chi_sym_tot(3, 3)
820  WRITE (unit_atoms, '(/T1,3(A,f10.4))') &
821  ' PV1=', nmr_env%chi_SI2ppmcgs*eig(1), &
822  ' PV2=', nmr_env%chi_SI2ppmcgs*eig(2), &
823  ' PV3=', nmr_env%chi_SI2ppmcgs*eig(3)
824  WRITE (unit_atoms, '(T1,A,F10.4,10X,A,F10.4)') &
825  ' ISO=', nmr_env%chi_SI2ppmcgs*chi_iso, &
826  'ANISO=', nmr_env%chi_SI2ppmcgs*chi_aniso
827  END IF
828 
829  CALL cp_print_key_finished_output(unit_atoms, logger, nmr_section,&
830  & "PRINT%CHI_TENSOR")
831  END IF ! print chi
832  !
833  ! Add the chi part to the shifts
834  cs_tot(:, :, :) = 0.0_dp
835  DO ir = 1, nat_print
836  iatom = cs_atom_list(ir)
837  rpos(1:3) = particle_set(iatom)%r(1:3)
838  atom_kind => particle_set(iatom)%atomic_kind
839  CALL get_atomic_kind(atom_kind, name=name, element_symbol=element_symbol)
840  cs_tot(:, :, ir) = chi_tensor_tmp(:, :)*nmr_env%chi_SI2shiftppm + cs(:, :, iatom)
841  IF (gapw) cs_tot(:, :, ir) = cs_tot(:, :, ir) + cs_loc(:, :, iatom)
842  END DO ! ir
843  IF (output_unit > 0) THEN
844  WRITE (output_unit, '(T2,A,E14.6)') 'CheckSum Shifts =', &
845  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  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  log_filename=.false.)
855 
856  nat_print = SIZE(cs_atom_list, 1)
857  IF (unit_atoms > 0) THEN
858  WRITE (title, '(A,1X,I5)') "Shielding atom at atomic positions. # tensors printed ", nat_print
859  WRITE (unit_atoms, '(T2,A)') title
860  DO ir = 1, nat_print
861  iatom = cs_atom_list(ir)
862  rpos(1:3) = particle_set(iatom)%r(1:3)
863  atom_kind => particle_set(iatom)%atomic_kind
864  CALL get_atomic_kind(atom_kind, name=name, element_symbol=element_symbol)
865  shift_sym_tot(:, :) = 0.5_dp*(cs_tot(:, :, ir) + transpose(cs_tot(:, :, ir)))
866  CALL diamat_all(shift_sym_tot, eig)
867  shift_iso = (eig(1) + eig(2) + eig(3))/3.0_dp
868  shift_aniso = eig(3) - (eig(2) + eig(1))/2.0_dp
869  !
870  WRITE (unit_atoms, '(T2,I5,A,2X,A2,2X,3f15.6)') iatom, trim(name), element_symbol, rpos(1:3)
871  !
872  IF (gapw) THEN
873  WRITE (unit_atoms, '(T1,A)') " SIGMA from SOFT J"
874  WRITE (unit_atoms, '(3(A,f10.4))') ' XX = ', cs(1, 1, iatom),&
875  & ' XY = ', cs(1, 2, iatom),&
876  & ' XZ = ', cs(1, 3, iatom)
877  WRITE (unit_atoms, '(3(A,f10.4))') ' YX = ', cs(2, 1, iatom),&
878  & ' YY = ', cs(2, 2, iatom),&
879  & ' YZ = ', cs(2, 3, iatom)
880  WRITE (unit_atoms, '(3(A,f10.4))') ' ZX = ', cs(3, 1, iatom),&
881  & ' ZY = ', cs(3, 2, iatom),&
882  & ' ZZ = ', cs(3, 3, iatom)
883  WRITE (unit_atoms, '(T1,A)') " SIGMA from LOCAL J"
884  WRITE (unit_atoms, '(3(A,f10.4))') ' XX = ', cs_loc(1, 1, iatom),&
885  & ' XY = ', cs_loc(1, 2, iatom),&
886  & ' XZ = ', cs_loc(1, 3, iatom)
887  WRITE (unit_atoms, '(3(A,f10.4))') ' YX = ', cs_loc(2, 1, iatom),&
888  & ' YY = ', cs_loc(2, 2, iatom),&
889  & ' YZ = ', cs_loc(2, 3, iatom)
890  WRITE (unit_atoms, '(3(A,f10.4))') ' ZX = ', cs_loc(3, 1, iatom),&
891  & ' ZY = ', cs_loc(3, 2, iatom),&
892  & ' ZZ = ', cs_loc(3, 3, iatom)
893  END IF
894  WRITE (unit_atoms, '(T1,A)') " SIGMA TOTAL"
895  WRITE (unit_atoms, '(3(A,f10.4))') ' XX = ', cs_tot(1, 1, ir),&
896  & ' XY = ', cs_tot(1, 2, ir),&
897  & ' XZ = ', cs_tot(1, 3, ir)
898  WRITE (unit_atoms, '(3(A,f10.4))') ' YX = ', cs_tot(2, 1, ir),&
899  & ' YY = ', cs_tot(2, 2, ir),&
900  & ' YZ = ', cs_tot(2, 3, ir)
901  WRITE (unit_atoms, '(3(A,f10.4))') ' ZX = ', cs_tot(3, 1, ir),&
902  & ' ZY = ', cs_tot(3, 2, ir),&
903  & ' ZZ = ', cs_tot(3, 3, ir)
904  WRITE (unit_atoms, '(T1,2(A,f12.4))') ' ISOTROPY = ', shift_iso,&
905  & ' ANISOTROPY = ', shift_aniso
906  END DO ! ir
907  END IF
908  CALL cp_print_key_finished_output(unit_atoms, logger, nmr_section,&
909  & "PRINT%SHIELDING_TENSOR")
910 
911  IF (do_nics) THEN
912  !
913  ! Add the chi part to the nics
914  cs_nics_tot(:, :, :) = 0.0_dp
915  DO ir = 1, n_nics
916  cs_nics_tot(:, :, ir) = chi_tensor_tmp(:, :)*nmr_env%chi_SI2shiftppm + cs_nics(:, :, ir)
917  IF (gapw) cs_nics_tot(:, :, ir) = cs_nics_tot(:, :, ir) + cs_nics_loc(:, :, ir)
918  END DO ! ir
919  IF (output_unit > 0) THEN
920  WRITE (output_unit, '(T2,A,E14.6)') 'CheckSum NICS =', &
921  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  log_filename=.false.)
927  IF (unit_nics > 0) THEN
928  WRITE (title, '(A,1X,I5)') "Shielding at nics positions. # tensors printed ", n_nics
929  WRITE (unit_nics, '(T2,A)') title
930  DO ir = 1, n_nics
931  shift_sym_tot(:, :) = 0.5_dp*(cs_nics_tot(:, :, ir) + transpose(cs_nics_tot(:, :, ir)))
932  CALL diamat_all(shift_sym_tot, eig)
933  shift_iso = (eig(1) + eig(2) + eig(3))/3.0_dp
934  shift_aniso = eig(3) - (eig(2) + eig(1))/2.0_dp
935  !
936  WRITE (unit_nics, '(T2,I5,2X,3f15.6)') ir, r_nics(1:3, ir)
937  !
938  IF (gapw) THEN
939  WRITE (unit_nics, '(T1,A)') " SIGMA from SOFT J"
940  WRITE (unit_nics, '(3(A,f10.4))') ' XX = ', cs_nics(1, 1, ir),&
941  & ' XY = ', cs_nics(1, 2, ir),&
942  & ' XZ = ', cs_nics(1, 3, ir)
943  WRITE (unit_nics, '(3(A,f10.4))') ' YX = ', cs_nics(2, 1, ir),&
944  & ' YY = ', cs_nics(2, 2, ir),&
945  & ' YZ = ', cs_nics(2, 3, ir)
946  WRITE (unit_nics, '(3(A,f10.4))') ' ZX = ', cs_nics(3, 1, ir),&
947  & ' ZY = ', cs_nics(3, 2, ir),&
948  & ' ZZ = ', cs_nics(3, 3, ir)
949  WRITE (unit_nics, '(T1,A)') " SIGMA from LOCAL J"
950  WRITE (unit_nics, '(3(A,f10.4))') ' XX = ', cs_nics_loc(1, 1, ir),&
951  & ' XY = ', cs_nics_loc(1, 2, ir),&
952  & ' XZ = ', cs_nics_loc(1, 3, ir)
953  WRITE (unit_nics, '(3(A,f10.4))') ' YX = ', cs_nics_loc(2, 1, ir),&
954  & ' YY = ', cs_nics_loc(2, 2, ir),&
955  & ' YZ = ', cs_nics_loc(2, 3, ir)
956  WRITE (unit_nics, '(3(A,f10.4))') ' ZX = ', cs_nics_loc(3, 1, ir),&
957  & ' ZY = ', cs_nics_loc(3, 2, ir),&
958  & ' ZZ = ', cs_nics_loc(3, 3, ir)
959  END IF
960  WRITE (unit_nics, '(T1,A)') " SIGMA TOTAL"
961  WRITE (unit_nics, '(3(A,f10.4))') ' XX = ', cs_nics_tot(1, 1, ir),&
962  & ' XY = ', cs_nics_tot(1, 2, ir),&
963  & ' XZ = ', cs_nics_tot(1, 3, ir)
964  WRITE (unit_nics, '(3(A,f10.4))') ' YX = ', cs_nics_tot(2, 1, ir),&
965  & ' YY = ', cs_nics_tot(2, 2, ir),&
966  & ' YZ = ', cs_nics_tot(2, 3, ir)
967  WRITE (unit_nics, '(3(A,f10.4))') ' ZX = ', cs_nics_tot(3, 1, ir),&
968  & ' ZY = ', cs_nics_tot(3, 2, ir),&
969  & ' ZZ = ', cs_nics_tot(3, 3, ir)
970  WRITE (unit_nics, '(T1,2(A,f12.4))') ' ISOTROPY = ', shift_iso,&
971  & ' ANISOTROPY = ', shift_aniso
972  END DO
973  END IF
974  CALL cp_print_key_finished_output(unit_nics, logger, nmr_section,&
975  & "PRINT%SHIELDING_TENSOR")
976  END IF
977  END IF ! print shift
978  !
979  ! clean up
980  DEALLOCATE (cs_tot)
981  IF (do_nics) THEN
982  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  END SUBROUTINE nmr_shift_print
993 
994 END MODULE qs_linres_nmr_shift
995 
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Definition: dumpdcd.F:1203
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition: grid_common.h:117
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Handles all functions related to the CELL.
Definition: cell_types.F:15
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
various routines to log and control the output. The idea is that decisions about where to log should ...
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
integer, parameter, public cp_p_file
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
objects that represent the structure of input sections and the data contained in an input section
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
integer, parameter, public default_string_length
Definition: kinds.F:57
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
complex(kind=dp), parameter, public gaussi
real(kind=dp), parameter, public twopi
Collection of simple mathematical functions and subroutines.
Definition: mathlib.F:15
subroutine, public diamat_all(a, eigval, dac)
Diagonalize the symmetric n by n matrix a using the LAPACK library. Only the upper triangle of matrix...
Definition: mathlib.F:376
Interface to the message passing library MPI.
Define the data structure for the particle information.
computes preconditioners, and implements methods to apply them currently used in qs_ot
container for various plainwaves related things
Definition: pw_env_types.F:14
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Definition: pw_env_types.F:113
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition: pw_pool_types.F:24
different utils that are useful to manipulate splines on the regular grid of a pw
logical function, public find_coeffs(values, coeffs, linOp, preconditioner, pool, eps_r, eps_x, max_iter, sumtype)
solves iteratively (CG) a systmes of linear equations linOp(coeffs)=values (for example those needed ...
subroutine, public pw_spline_precond_release(preconditioner)
releases the preconditioner
subroutine, public pw_spline_precond_create(preconditioner, precond_kind, pool, pbc, transpose)
...
subroutine, public pw_spline_do_precond(preconditioner, in_v, out_v)
applies the preconditioner to the system of equations to find the coefficients of the spline
subroutine, public pw_spline_precond_set_kind(preconditioner, precond_kind, pbc, transpose)
switches the types of precoditioner to use
real(kind=dp) function, public eval_interp_spl3_pbc(vec, pw)
Evaluates the PBC interpolated Spline (pw) function on the generic input vector (vec)
subroutine, public spl3_pbc(pw_in, pw_out)
...
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_r3d_rs_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, U_of_dft_plus_u, J_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, J0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
given the response wavefunctions obtained by the application of the (rxp), p, and ((dk-dl)xp) operato...
subroutine, public mult_g_ov_g2_grid(pw_pool, rho_gspace, funcG_times_rho, idir, my_chi)
Given the current density on the PW grid in reciprcal space (obtained by FFT), calculate the integral...
from the response current density calculates the shift tensor and the susceptibility
subroutine, public nmr_shift(nmr_env, current_env, qs_env, iB)
...
subroutine, public nmr_shift_print(nmr_env, current_env, qs_env)
Shielding tensor and Chi are printed into a file if required from input It is possible to print only ...
Calculate the operators p rxp and D needed in the optimization of the different contribution of the f...
Definition: qs_linres_op.F:18
subroutine, public set_vecp(i1, i2, i3)
...
real(dp) function, public fac_vecp(a, b, c)
...
subroutine, public set_vecp_rev(i1, i2, i3)
...
Type definitiona for linear response calculations.
subroutine, public get_nmr_env(nmr_env, n_nics, cs_atom_list, do_calc_cs_atom, r_nics, chemical_shift, chemical_shift_loc, chemical_shift_nics_loc, chemical_shift_nics, shift_gapw_radius, do_nics, interpolate_shift)
...
subroutine, public get_current_env(current_env, simple_done, simple_converged, full_done, nao, nstates, gauge, list_cubes, statetrueindex, gauge_name, basisfun_center, nbr_center, center_list, centers_set, psi1_p, psi1_rxp, psi1_D, p_psi0, rxp_psi0, jrho1_atom_set, jrho1_set, chi_tensor, chi_tensor_loc, gauge_atom_radius, rs_gauge, use_old_gauge_atom, chi_pbc, psi0_order)
...
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition: qs_rho_types.F:18
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition: qs_rho_types.F:229
All kind of helpful little routines.
Definition: util.F:14
pure integer function, dimension(2), public get_limit(m, n, me)
divide m entries into n parts, return size of part me
Definition: util.F:333