1!--------------------------------------------------------------------------------------------------!
2!   CP2K: A general program to perform molecular dynamics simulations                              !
3!   Copyright (C) 2000 - 2020  CP2K developers group                                               !
4!--------------------------------------------------------------------------------------------------!
5
6! **************************************************************************************************
7!> \brief given the response wavefunctions obtained by the application
8!>      of the (rxp), p, and ((dk-dl)xp) operators,
9!>      here the current density vector (jx, jy, jz)
10!>      is computed for the 3 directions of the magnetic field (Bx, By, Bz)
11!> \par History
12!>      created 02-2006 [MI]
13!> \author MI
14! **************************************************************************************************
15MODULE qs_linres_nmr_epr_common_utils
16   USE cell_types,                      ONLY: cell_type
17   USE kinds,                           ONLY: dp
18   USE mathconstants,                   ONLY: gaussi
19   USE pw_grid_types,                   ONLY: pw_grid_type
20   USE pw_methods,                      ONLY: pw_transfer
21   USE pw_pool_types,                   ONLY: pw_pool_create_pw,&
22                                              pw_pool_give_back_pw,&
23                                              pw_pool_type
24   USE pw_types,                        ONLY: COMPLEXDATA1D,&
25                                              RECIPROCALSPACE,&
26                                              pw_p_type,&
27                                              pw_type
28#include "./base/base_uses.f90"
29
30   IMPLICIT NONE
31
32   PRIVATE
33
34   ! *** Public subroutines ***
35   PUBLIC :: mult_G_ov_G2_grid
36
37   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_linres_nmr_epr_common_utils'
38
39CONTAINS
40
41! **************************************************************************************************
42!> \brief Given the current density on the PW grid in reciprcal space
43!>       (obtained by FFT), calculate the integral
44!>         \int_{r}[ ((r-r') x j(r))/|r-r'|^3 ] = Bind(r')
45!>       which in reciprcal space reads  (for G/=0)
46!>          i G/|G|^2 x J(G)
47!> \param cell ...
48!> \param pw_pool ...
49!> \param rho_gspace ...
50!> \param funcG_times_rho ...
51!> \param idir ...
52!> \param my_chi ...
53!> \author MI
54!> \note
55!>      The G=0 component is not comnputed here, but can be evaluated
56!>      through the susceptibility and added to the shift in a second time
57!>
58!>      This method would not work for a non periodic system
59!>      It should be generalized like the calculation of Hartree
60! **************************************************************************************************
61   SUBROUTINE mult_G_ov_G2_grid(cell, pw_pool, rho_gspace, funcG_times_rho, idir, my_chi)
62
63      TYPE(cell_type), POINTER                           :: cell
64      TYPE(pw_pool_type), POINTER                        :: pw_pool
65      TYPE(pw_p_type), POINTER                           :: rho_gspace
66      TYPE(pw_p_type)                                    :: funcG_times_rho
67      INTEGER, INTENT(IN)                                :: idir
68      REAL(dp), INTENT(IN)                               :: my_chi
69
70      INTEGER                                            :: handle, ig, ng
71      REAL(dp)                                           :: g2
72      TYPE(pw_grid_type), POINTER                        :: grid
73      CHARACTER(len=*), PARAMETER :: routineN = 'mult_G_ov_G2_grid', &
74         routineP = moduleN//':'//routineN
75
76      TYPE(pw_type), POINTER                             :: frho, influence_fn
77
78      CALL timeset(routineN, handle)
79
80      CPASSERT(ASSOCIATED(cell))
81
82      CALL pw_pool_create_pw(pw_pool, influence_fn, &
83                             use_data=COMPLEXDATA1D, in_space=RECIPROCALSPACE)
84
85      grid => influence_fn%pw_grid
86      DO ig = grid%first_gne0, grid%ngpts_cut_local
87         g2 = grid%gsq(ig)
88         influence_fn%cc(ig) = gaussi*grid%g(idir, ig)/g2
89      END DO ! ig
90      IF (grid%have_g0) influence_fn%cc(1) = 0.0_dp
91
92      frho => funcG_times_rho%pw
93      CALL pw_transfer(rho_gspace%pw, frho)
94
95      ng = SIZE(grid%gsq)
96      frho%cc(1:ng) = frho%cc(1:ng)*influence_fn%cc(1:ng)
97      IF (grid%have_g0) frho%cc(1) = my_chi
98
99      CALL pw_pool_give_back_pw(pw_pool, influence_fn, &
100                                accept_non_compatible=.TRUE.)
101
102      CALL timestop(handle)
103
104   END SUBROUTINE mult_G_ov_G2_grid
105
106END MODULE qs_linres_nmr_epr_common_utils
107