1!--------------------------------------------------------------------------------------------------!
2!   CP2K: A general program to perform molecular dynamics simulations                              !
3!   Copyright (C) 2000 - 2020  CP2K developers group                                               !
4!--------------------------------------------------------------------------------------------------!
5
6! **************************************************************************************************
7!> \brief pw_methods
8!> \author CJM
9! **************************************************************************************************
10MODULE ewald_pw_methods
11   USE dg_rho0_types,                   ONLY: dg_rho0_get,&
12                                              dg_rho0_init,&
13                                              dg_rho0_set,&
14                                              dg_rho0_type
15   USE dg_types,                        ONLY: dg_get,&
16                                              dg_set,&
17                                              dg_type
18   USE dgs,                             ONLY: dg_grid_change
19   USE ewald_environment_types,         ONLY: ewald_env_get,&
20                                              ewald_env_set,&
21                                              ewald_environment_type
22   USE ewald_pw_types,                  ONLY: ewald_pw_get,&
23                                              ewald_pw_set,&
24                                              ewald_pw_type
25   USE input_section_types,             ONLY: section_vals_type
26   USE kinds,                           ONLY: dp
27   USE pw_grid_types,                   ONLY: pw_grid_type
28   USE pw_grids,                        ONLY: pw_grid_change
29   USE pw_poisson_methods,              ONLY: pw_poisson_set
30! CJM, pw_poisson_rebuild
31   USE pw_poisson_read_input,           ONLY: pw_poisson_read_parameters
32   USE pw_poisson_types,                ONLY: do_ewald_ewald,&
33                                              do_ewald_none,&
34                                              do_ewald_pme,&
35                                              do_ewald_spme,&
36                                              pw_poisson_create,&
37                                              pw_poisson_parameter_type,&
38                                              pw_poisson_release,&
39                                              pw_poisson_type
40   USE pw_pool_types,                   ONLY: pw_pool_p_type,&
41                                              pw_pool_type
42#include "./base/base_uses.f90"
43
44   IMPLICIT NONE
45
46   PRIVATE
47
48   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ewald_pw_methods'
49
50   PUBLIC :: ewald_pw_grid_update
51
52CONTAINS
53
54! **************************************************************************************************
55!> \brief Rescales pw_grids for given box, if necessary
56!> \param ewald_pw ...
57!> \param ewald_env ...
58!> \param cell_hmat ...
59!> \par History
60!>      none
61!> \author JGH (15-Mar-2001)
62! **************************************************************************************************
63   SUBROUTINE ewald_pw_grid_update(ewald_pw, ewald_env, cell_hmat)
64      TYPE(ewald_pw_type), POINTER                       :: ewald_pw
65      TYPE(ewald_environment_type), POINTER              :: ewald_env
66      REAL(KIND=dp), DIMENSION(3, 3)                     :: cell_hmat
67
68      CHARACTER(len=*), PARAMETER :: routineN = 'ewald_pw_grid_update', &
69         routineP = moduleN//':'//routineN
70
71      INTEGER                                            :: ewald_type, o_spline
72      REAL(dp)                                           :: alpha
73      REAL(KIND=dp), DIMENSION(3, 3)                     :: old_cell_hmat
74      TYPE(dg_type), POINTER                             :: dg
75      TYPE(pw_poisson_parameter_type)                    :: poisson_params
76      TYPE(pw_poisson_type), POINTER                     :: poisson_env
77      TYPE(pw_pool_p_type), DIMENSION(:), POINTER        :: pw_pools
78      TYPE(pw_pool_type), POINTER                        :: pw_big_pool, pw_small_pool
79      TYPE(section_vals_type), POINTER                   :: poisson_section
80
81      NULLIFY (pw_big_pool)
82      NULLIFY (pw_small_pool)
83      NULLIFY (dg, poisson_env, poisson_section)
84
85      CALL ewald_env_get(ewald_env, ewald_type=ewald_type, &
86                         alpha=alpha, o_spline=o_spline, &
87                         poisson_section=poisson_section, &
88                         cell_hmat=old_cell_hmat)
89
90      IF (ALL(cell_hmat == old_cell_hmat)) RETURN ! rebuild not needed
91
92      CALL ewald_env_set(ewald_env, cell_hmat=cell_hmat)
93
94      SELECT CASE (ewald_type)
95      CASE (do_ewald_ewald)
96         CALL ewald_pw_get(ewald_pw, pw_big_pool=pw_big_pool, &
97                           dg=dg, poisson_env=poisson_env)
98         CALL pw_grid_change(cell_hmat, pw_big_pool%pw_grid)
99         CALL ewald_pw_rho0_setup(ewald_env, pw_big_pool%pw_grid, dg)
100         CALL pw_poisson_release(poisson_env)
101         CALL ewald_pw_set(ewald_pw, pw_big_pool=pw_big_pool, dg=dg, &
102                           poisson_env=poisson_env)
103      CASE (do_ewald_pme)
104         CALL ewald_pw_get(ewald_pw, pw_big_pool=pw_big_pool, &
105                           pw_small_pool=pw_small_pool, dg=dg, &
106                           poisson_env=poisson_env)
107         IF (.NOT. ASSOCIATED(poisson_env)) THEN
108            CALL pw_poisson_create(poisson_env)
109            CALL ewald_pw_set(ewald_pw, poisson_env=poisson_env)
110         END IF
111         CALL pw_grid_change(cell_hmat, pw_big_pool%pw_grid)
112         CALL dg_grid_change(cell_hmat, pw_big_pool%pw_grid, pw_small_pool%pw_grid)
113         CALL ewald_pw_rho0_setup(ewald_env, pw_small_pool%pw_grid, dg)
114         CALL ewald_pw_set(ewald_pw, pw_big_pool=pw_big_pool, &
115                           pw_small_pool=pw_small_pool, dg=dg, &
116                           poisson_env=poisson_env)
117      CASE (do_ewald_spme)
118         CALL ewald_pw_get(ewald_pw, pw_big_pool=pw_big_pool, &
119                           poisson_env=poisson_env)
120         IF (.NOT. ASSOCIATED(poisson_env)) THEN
121            CALL pw_poisson_create(poisson_env)
122         END IF
123         CALL pw_grid_change(cell_hmat, pw_big_pool%pw_grid)
124         CALL ewald_pw_set(ewald_pw, pw_big_pool=pw_big_pool, &
125                           poisson_env=poisson_env)
126      CASE (do_ewald_none)
127      CASE default
128         CPABORT("")
129      END SELECT
130      IF (ASSOCIATED(poisson_env)) THEN
131         ALLOCATE (pw_pools(1))
132         pw_pools(1)%pool => pw_big_pool
133         CALL pw_poisson_read_parameters(poisson_section, poisson_params)
134         poisson_params%ewald_type = ewald_type
135         poisson_params%ewald_o_spline = o_spline
136         poisson_params%ewald_alpha = alpha
137         CALL pw_poisson_set(poisson_env, cell_hmat=cell_hmat, parameters=poisson_params, &
138                             use_level=1, pw_pools=pw_pools)
139         DEALLOCATE (pw_pools)
140      END IF
141
142   END SUBROUTINE ewald_pw_grid_update
143
144! **************************************************************************************************
145!> \brief Calculates the Fourier transform of the "Ewald function"
146!> \param ewald_env ...
147!> \param pw_grid ...
148!> \param dg ...
149!> \par History
150!>      none
151!> \author JGH (15-Mar-2001)
152! **************************************************************************************************
153   SUBROUTINE ewald_pw_rho0_setup(ewald_env, pw_grid, dg)
154      TYPE(ewald_environment_type), POINTER              :: ewald_env
155      TYPE(pw_grid_type), POINTER                        :: pw_grid
156      TYPE(dg_type), POINTER                             :: dg
157
158      CHARACTER(len=*), PARAMETER :: routineN = 'ewald_pw_rho0_setup', &
159         routineP = moduleN//':'//routineN
160
161      INTEGER                                            :: ewald_type, grid_index
162      REAL(dp)                                           :: alpha
163      REAL(dp), POINTER                                  :: gcc(:), zet(:)
164      TYPE(dg_rho0_type), POINTER                        :: dg_rho0
165
166      CALL ewald_env_get(ewald_env, alpha=alpha, ewald_type=ewald_type)
167      CALL dg_get(dg, dg_rho0=dg_rho0)
168      CALL dg_rho0_get(dg_rho0, gcc=gcc, zet=zet)
169
170! This is the first ( and only ) double grid
171      grid_index = 1
172
173      IF (.NOT. ASSOCIATED(zet)) THEN
174         ALLOCATE (zet(1))
175      ENDIF
176
177! No contracted Gaussians are used here
178      NULLIFY (gcc)
179
180      zet(1) = alpha
181      CALL dg_rho0_set(dg_rho0, TYPE=ewald_type, zet=zet)
182
183      CALL dg_rho0_init(dg_rho0, pw_grid)
184
185      CALL dg_set(dg, dg_rho0=dg_rho0, grid_index=grid_index)
186
187   END SUBROUTINE ewald_pw_rho0_setup
188
189END MODULE ewald_pw_methods
190
191