LCOV - code coverage report
Current view: top level - src - qs_linres_nmr_shift.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 100.0 % 446 446
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 6 6

            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              :    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              : 
        

Generated by: LCOV version 2.0-1