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_types
8!> \author CJM
9! **************************************************************************************************
10MODULE ewald_pw_types
11   USE ao_util,                         ONLY: exp_radius
12   USE cell_types,                      ONLY: cell_type
13   USE cp_log_handling,                 ONLY: cp_get_default_logger,&
14                                              cp_logger_type
15   USE cp_output_handling,              ONLY: cp_print_key_finished_output,&
16                                              cp_print_key_unit_nr
17   USE cp_para_types,                   ONLY: cp_para_env_type
18   USE cp_realspace_grid_init,          ONLY: init_input_type
19   USE dg_types,                        ONLY: dg_create,&
20                                              dg_release,&
21                                              dg_retain,&
22                                              dg_type
23   USE dgs,                             ONLY: dg_pme_grid_setup
24   USE ewald_environment_types,         ONLY: ewald_env_get,&
25                                              ewald_environment_type
26   USE input_section_types,             ONLY: section_vals_get_subs_vals,&
27                                              section_vals_type
28   USE kinds,                           ONLY: dp
29   USE mathconstants,                   ONLY: pi
30   USE message_passing,                 ONLY: mp_comm_self
31   USE pw_grid_types,                   ONLY: HALFSPACE,&
32                                              pw_grid_type
33   USE pw_grids,                        ONLY: pw_grid_create,&
34                                              pw_grid_release,&
35                                              pw_grid_setup
36   USE pw_poisson_methods,              ONLY: pw_poisson_set
37   USE pw_poisson_read_input,           ONLY: pw_poisson_read_parameters
38   USE pw_poisson_types,                ONLY: &
39        do_ewald_ewald, do_ewald_none, do_ewald_pme, do_ewald_spme, pw_poisson_create, &
40        pw_poisson_parameter_type, pw_poisson_release, pw_poisson_retain, pw_poisson_type
41   USE pw_pool_types,                   ONLY: pw_pool_create,&
42                                              pw_pool_p_type,&
43                                              pw_pool_release,&
44                                              pw_pool_retain,&
45                                              pw_pool_type
46   USE realspace_grid_types,            ONLY: &
47        realspace_grid_desc_type, realspace_grid_input_type, realspace_grid_type, rs_grid_create, &
48        rs_grid_create_descriptor, rs_grid_print, rs_grid_release, rs_grid_release_descriptor, &
49        rs_grid_retain_descriptor
50#include "./base/base_uses.f90"
51
52   IMPLICIT NONE
53
54   PRIVATE
55   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ewald_pw_types'
56   INTEGER, PRIVATE, SAVE :: last_ewald_pw_id_nr = 0
57   PUBLIC :: ewald_pw_type, ewald_pw_release, &
58             ewald_pw_retain, ewald_pw_create, &
59             ewald_pw_get, ewald_pw_set
60
61! **************************************************************************************************
62   TYPE ewald_pw_type
63      PRIVATE
64      INTEGER :: ref_count, id_nr
65      TYPE(pw_pool_type), POINTER       :: pw_small_pool
66      TYPE(pw_pool_type), POINTER       :: pw_big_pool
67      TYPE(realspace_grid_desc_type), POINTER    :: rs_desc
68      TYPE(pw_poisson_type), POINTER    :: poisson_env
69      TYPE(dg_type), POINTER            :: dg
70   END TYPE ewald_pw_type
71
72CONTAINS
73
74! **************************************************************************************************
75!> \brief retains the structure ewald_pw_type
76!> \param ewald_pw ...
77! **************************************************************************************************
78   SUBROUTINE ewald_pw_retain(ewald_pw)
79      TYPE(ewald_pw_type), POINTER                       :: ewald_pw
80
81      CHARACTER(len=*), PARAMETER :: routineN = 'ewald_pw_retain', &
82         routineP = moduleN//':'//routineN
83
84      CPASSERT(ASSOCIATED(ewald_pw))
85      CPASSERT(ewald_pw%ref_count > 0)
86      ewald_pw%ref_count = ewald_pw%ref_count + 1
87   END SUBROUTINE ewald_pw_retain
88
89! **************************************************************************************************
90!> \brief creates the structure ewald_pw_type
91!> \param ewald_pw ...
92!> \param ewald_env ...
93!> \param cell ...
94!> \param cell_ref ...
95!> \param print_section ...
96! **************************************************************************************************
97   SUBROUTINE ewald_pw_create(ewald_pw, ewald_env, cell, cell_ref, print_section)
98      TYPE(ewald_pw_type), POINTER                       :: ewald_pw
99      TYPE(ewald_environment_type), POINTER              :: ewald_env
100      TYPE(cell_type), POINTER                           :: cell, cell_ref
101      TYPE(section_vals_type), POINTER                   :: print_section
102
103      CHARACTER(len=*), PARAMETER :: routineN = 'ewald_pw_create', &
104         routineP = moduleN//':'//routineN
105
106      TYPE(dg_type), POINTER                             :: dg
107
108      NULLIFY (dg)
109      ALLOCATE (ewald_pw)
110      NULLIFY (ewald_pw%pw_big_pool)
111      NULLIFY (ewald_pw%pw_small_pool)
112      NULLIFY (ewald_pw%rs_desc)
113      NULLIFY (ewald_pw%poisson_env)
114      CALL dg_create(dg)
115      ewald_pw%dg => dg
116      ewald_pw%ref_count = 1
117      last_ewald_pw_id_nr = last_ewald_pw_id_nr + 1
118      ewald_pw%id_nr = last_ewald_pw_id_nr
119      CALL ewald_pw_init(ewald_pw, ewald_env, cell, cell_ref, print_section)
120   END SUBROUTINE ewald_pw_create
121
122!****f* ewald_pw_types/ewald_pw_release [1.0] *
123
124! **************************************************************************************************
125!> \brief releases the memory used by the ewald_pw
126!> \param ewald_pw ...
127! **************************************************************************************************
128   SUBROUTINE ewald_pw_release(ewald_pw)
129      TYPE(ewald_pw_type), POINTER                       :: ewald_pw
130
131      CHARACTER(len=*), PARAMETER :: routineN = 'ewald_pw_release', &
132         routineP = moduleN//':'//routineN
133
134      INTEGER                                            :: handle
135
136      CALL timeset(routineN, handle)
137      IF (ASSOCIATED(ewald_pw)) THEN
138         CPASSERT(ewald_pw%ref_count > 0)
139         ewald_pw%ref_count = ewald_pw%ref_count - 1
140         IF (ewald_pw%ref_count < 1) THEN
141            CALL pw_pool_release(ewald_pw%pw_small_pool)
142            CALL pw_pool_release(ewald_pw%pw_big_pool)
143            CALL rs_grid_release_descriptor(ewald_pw%rs_desc)
144            CALL pw_poisson_release(ewald_pw%poisson_env)
145            CALL dg_release(ewald_pw%dg)
146            DEALLOCATE (ewald_pw)
147         END IF
148      END IF
149      NULLIFY (ewald_pw)
150      CALL timestop(handle)
151   END SUBROUTINE ewald_pw_release
152
153!****** ewald_pw_types/ewald_pw_init [1.0] *
154
155! **************************************************************************************************
156!> \brief ...
157!> \param ewald_pw ...
158!> \param ewald_env ...
159!> \param cell ...
160!> \param cell_ref ...
161!> \param print_section ...
162!> \par History
163!>      JGH (12-Jan-2001): Added SPME part
164!>      JGH (15-Mar-2001): Work newly distributed between initialize, setup,
165!>                         and force routine
166!> \author CJM
167! **************************************************************************************************
168   SUBROUTINE ewald_pw_init(ewald_pw, ewald_env, cell, cell_ref, print_section)
169      TYPE(ewald_pw_type), POINTER                       :: ewald_pw
170      TYPE(ewald_environment_type), POINTER              :: ewald_env
171      TYPE(cell_type), POINTER                           :: cell, cell_ref
172      TYPE(section_vals_type), POINTER                   :: print_section
173
174      CHARACTER(len=*), PARAMETER :: routineN = 'ewald_pw_init', routineP = moduleN//':'//routineN
175
176      INTEGER                                            :: bo(2, 3), ewald_type, gmax(3), handle, &
177                                                            npts_s(3), ns_max, o_spline, &
178                                                            output_unit
179      REAL(KIND=dp)                                      :: alpha, alphasq, cutoff_radius, epsilon, &
180                                                            norm
181      TYPE(cp_logger_type), POINTER                      :: logger
182      TYPE(cp_para_env_type), POINTER                    :: para_env
183      TYPE(pw_grid_type), POINTER                        :: pw_big_grid, pw_small_grid
184      TYPE(pw_poisson_parameter_type)                    :: poisson_params
185      TYPE(pw_pool_p_type), DIMENSION(:), POINTER        :: pw_pools
186      TYPE(pw_pool_type), POINTER                        :: pw_pool
187      TYPE(realspace_grid_desc_type), POINTER            :: rs_desc
188      TYPE(realspace_grid_input_type)                    :: input_settings
189      TYPE(realspace_grid_type), POINTER                 :: rs
190      TYPE(section_vals_type), POINTER                   :: poisson_section, rs_grid_section
191
192      CALL timeset(routineN, handle)
193
194      NULLIFY (pw_big_grid)
195      NULLIFY (pw_small_grid, poisson_section)
196
197      CPASSERT(ASSOCIATED(ewald_pw))
198      CPASSERT(ASSOCIATED(ewald_env))
199      CPASSERT(ASSOCIATED(cell))
200      CPASSERT(ewald_pw%ref_count > 0)
201      CALL ewald_env_get(ewald_env=ewald_env, &
202                         para_env=para_env, &
203                         gmax=gmax, alpha=alpha, &
204                         ns_max=ns_max, &
205                         ewald_type=ewald_type, &
206                         o_spline=o_spline, &
207                         poisson_section=poisson_section, &
208                         epsilon=epsilon)
209
210      rs_grid_section => section_vals_get_subs_vals(poisson_section, "EWALD%RS_GRID")
211
212      SELECT CASE (ewald_type)
213      CASE (do_ewald_ewald)
214         ! set up Classic EWALD sum
215         logger => cp_get_default_logger()
216         output_unit = cp_print_key_unit_nr(logger, print_section, "", extension=".Log")
217         CALL pw_grid_create(pw_big_grid, mp_comm_self)
218
219         IF (ANY(gmax == 2*(gmax/2))) THEN
220            CPABORT("gmax has to be odd.")
221         END IF
222         bo(1, :) = -gmax/2
223         bo(2, :) = +gmax/2
224         CALL pw_grid_setup(cell_ref%hmat, pw_big_grid, grid_span=HALFSPACE, bounds=bo, spherical=.TRUE., &
225                            fft_usage=.FALSE., iounit=output_unit)
226         NULLIFY (pw_pool)
227         CALL pw_pool_create(pw_pool, pw_grid=pw_big_grid)
228         ewald_pw%pw_big_pool => pw_pool
229         CALL pw_pool_retain(ewald_pw%pw_big_pool)
230         CALL pw_pool_release(pw_pool)
231         CALL pw_grid_release(pw_big_grid)
232         CALL cp_print_key_finished_output(output_unit, logger, print_section, "")
233
234      CASE (do_ewald_pme)
235         ! set up Particle-Mesh EWALD sum
236         logger => cp_get_default_logger()
237         output_unit = cp_print_key_unit_nr(logger, print_section, "", extension=".Log")
238         IF (.NOT. ASSOCIATED(ewald_pw%poisson_env)) THEN
239            CALL pw_poisson_create(ewald_pw%poisson_env)
240         END IF
241         CALL pw_grid_create(pw_small_grid, mp_comm_self)
242         CALL pw_grid_create(pw_big_grid, para_env%group)
243         IF (ns_max == 2*(ns_max/2)) THEN
244            CPABORT("ns_max has to be odd.")
245         END IF
246         npts_s(:) = ns_max
247         ! compute cut-off radius
248         alphasq = alpha**2
249         norm = (2.0_dp*alphasq/pi)**(1.5_dp)
250         cutoff_radius = exp_radius(0, 2.0_dp*alphasq, epsilon, norm)
251
252         CALL dg_pme_grid_setup(cell_ref%hmat, npts_s, cutoff_radius, &
253                                pw_small_grid, pw_big_grid, rs_dims=(/para_env%num_pe, 1/), &
254                                iounit=output_unit, fft_usage=.TRUE.)
255         ! Write some useful info
256         IF (output_unit > 0) THEN
257            WRITE (output_unit, '( A,T71,E10.4 )') &
258               ' EWALD| Gaussian tolerance (effective) ', epsilon
259            WRITE (output_unit, '( A,T63,3I6 )') &
260               ' EWALD| Small box grid ', pw_small_grid%npts
261            WRITE (output_unit, '( A,T63,3I6 )') &
262               ' EWALD| Full box grid ', pw_big_grid%npts
263         END IF
264
265         ! pw pools initialized
266         NULLIFY (pw_pool)
267         CALL pw_pool_create(pw_pool, pw_grid=pw_big_grid)
268         ewald_pw%pw_big_pool => pw_pool
269         CALL pw_pool_retain(ewald_pw%pw_big_pool)
270         CALL pw_pool_release(pw_pool)
271
272         NULLIFY (pw_pool)
273         CALL pw_pool_create(pw_pool, pw_grid=pw_small_grid)
274         ewald_pw%pw_small_pool => pw_pool
275         CALL pw_pool_retain(ewald_pw%pw_small_pool)
276         CALL pw_pool_release(pw_pool)
277
278         NULLIFY (rs_desc)
279         CALL init_input_type(input_settings, nsmax=MAXVAL(pw_small_grid%npts(1:3)), &
280                              rs_grid_section=rs_grid_section, ilevel=1, &
281                              higher_grid_layout=(/-1, -1, -1/))
282         CALL rs_grid_create_descriptor(rs_desc, pw_big_grid, input_settings)
283
284         CALL rs_grid_create(rs, rs_desc)
285         CALL rs_grid_print(rs, output_unit)
286         CALL rs_grid_release(rs)
287
288         CALL cp_print_key_finished_output(output_unit, logger, print_section, "")
289
290         ewald_pw%rs_desc => rs_desc
291
292         CALL rs_grid_retain_descriptor(ewald_pw%rs_desc)
293         CALL rs_grid_release_descriptor(rs_desc)
294
295         CALL pw_grid_release(pw_small_grid)
296         CALL pw_grid_release(pw_big_grid)
297
298      CASE (do_ewald_spme)
299         ! set up the Smooth-Particle-Mesh EWALD sum
300         logger => cp_get_default_logger()
301         output_unit = cp_print_key_unit_nr(logger, print_section, "", extension=".Log")
302         IF (.NOT. ASSOCIATED(ewald_pw%poisson_env)) THEN
303            CALL pw_poisson_create(ewald_pw%poisson_env)
304         END IF
305         CALL pw_grid_create(pw_big_grid, para_env%group)
306         npts_s = gmax
307         CALL pw_grid_setup(cell_ref%hmat, pw_big_grid, grid_span=HALFSPACE, npts=npts_s, spherical=.TRUE., &
308                            rs_dims=(/para_env%num_pe, 1/), iounit=output_unit, fft_usage=.TRUE.)
309
310         ! pw pools initialized
311         NULLIFY (pw_pool)
312         CALL pw_pool_create(pw_pool, pw_grid=pw_big_grid)
313         ewald_pw%pw_big_pool => pw_pool
314         CALL pw_pool_retain(ewald_pw%pw_big_pool)
315         CALL pw_pool_release(pw_pool)
316
317         NULLIFY (rs_desc)
318         CALL init_input_type(input_settings, nsmax=o_spline, &
319                              rs_grid_section=rs_grid_section, ilevel=1, &
320                              higher_grid_layout=(/-1, -1, -1/))
321         CALL rs_grid_create_descriptor(rs_desc, pw_big_grid, input_settings)
322
323         CALL rs_grid_create(rs, rs_desc)
324         CALL rs_grid_print(rs, output_unit)
325         CALL rs_grid_release(rs)
326         CALL cp_print_key_finished_output(output_unit, logger, print_section, "")
327
328         ewald_pw%rs_desc => rs_desc
329
330         CALL rs_grid_retain_descriptor(ewald_pw%rs_desc)
331         CALL rs_grid_release_descriptor(rs_desc)
332
333         CALL pw_grid_release(pw_big_grid)
334      CASE (do_ewald_none)
335         ! No EWALD sums..
336      CASE default
337         CPABORT("")
338      END SELECT
339      ! Poisson Environment
340      IF (ASSOCIATED(ewald_pw%poisson_env)) THEN
341         ALLOCATE (pw_pools(1))
342         pw_pools(1)%pool => ewald_pw%pw_big_pool
343         CALL pw_poisson_read_parameters(poisson_section, poisson_params)
344         poisson_params%ewald_type = ewald_type
345         poisson_params%ewald_o_spline = o_spline
346         poisson_params%ewald_alpha = alpha
347         CALL pw_poisson_set(ewald_pw%poisson_env, cell_hmat=cell%hmat, parameters=poisson_params, &
348                             use_level=1, pw_pools=pw_pools)
349         DEALLOCATE (pw_pools)
350      END IF
351      CALL timestop(handle)
352   END SUBROUTINE ewald_pw_init
353
354!****** ewald_pw_types/ewald_pw_get [1.0] *
355
356! **************************************************************************************************
357!> \brief get the ewald_pw environment to the correct program.
358!> \param ewald_pw ...
359!> \param pw_big_pool ...
360!> \param pw_small_pool ...
361!> \param rs_desc ...
362!> \param poisson_env ...
363!> \param dg ...
364!> \author CJM
365! **************************************************************************************************
366   SUBROUTINE ewald_pw_get(ewald_pw, pw_big_pool, pw_small_pool, rs_desc, poisson_env, dg)
367
368      TYPE(ewald_pw_type), POINTER                       :: ewald_pw
369      TYPE(pw_pool_type), OPTIONAL, POINTER              :: pw_big_pool, pw_small_pool
370      TYPE(realspace_grid_desc_type), OPTIONAL, POINTER  :: rs_desc
371      TYPE(pw_poisson_type), OPTIONAL, POINTER           :: poisson_env
372      TYPE(dg_type), OPTIONAL, POINTER                   :: dg
373
374      CHARACTER(LEN=*), PARAMETER :: routineN = 'ewald_pw_get', routineP = moduleN//':'//routineN
375
376      IF (PRESENT(poisson_env)) poisson_env => ewald_pw%poisson_env
377      IF (PRESENT(pw_big_pool)) pw_big_pool => ewald_pw%pw_big_pool
378      IF (PRESENT(pw_small_pool)) pw_small_pool => ewald_pw%pw_small_pool
379      IF (PRESENT(rs_desc)) rs_desc => ewald_pw%rs_desc
380      IF (PRESENT(dg)) dg => ewald_pw%dg
381
382   END SUBROUTINE ewald_pw_get
383
384!****** ewald_pw_types/ewald_pw_set [1.0] *
385
386! **************************************************************************************************
387!> \brief set the ewald_pw environment to the correct program.
388!> \param ewald_pw ...
389!> \param pw_big_pool ...
390!> \param pw_small_pool ...
391!> \param rs_desc ...
392!> \param dg ...
393!> \param poisson_env ...
394!> \author CJM
395! **************************************************************************************************
396   SUBROUTINE ewald_pw_set(ewald_pw, pw_big_pool, pw_small_pool, rs_desc, dg, &
397                           poisson_env)
398
399      TYPE(ewald_pw_type), POINTER                       :: ewald_pw
400      TYPE(pw_pool_type), OPTIONAL, POINTER              :: pw_big_pool, pw_small_pool
401      TYPE(realspace_grid_desc_type), OPTIONAL, POINTER  :: rs_desc
402      TYPE(dg_type), OPTIONAL, POINTER                   :: dg
403      TYPE(pw_poisson_type), OPTIONAL, POINTER           :: poisson_env
404
405      CHARACTER(LEN=*), PARAMETER :: routineN = 'ewald_pw_set', routineP = moduleN//':'//routineN
406
407      IF (PRESENT(pw_big_pool)) THEN
408         CALL pw_pool_retain(pw_big_pool)
409         CALL pw_pool_release(ewald_pw%pw_big_pool)
410         ewald_pw%pw_big_pool => pw_big_pool
411      ENDIF
412      IF (PRESENT(pw_small_pool)) THEN
413         CALL pw_pool_retain(pw_small_pool)
414         CALL pw_pool_release(ewald_pw%pw_small_pool)
415         ewald_pw%pw_small_pool => pw_small_pool
416      ENDIF
417      IF (PRESENT(rs_desc)) THEN
418         CALL rs_grid_retain_descriptor(rs_desc)
419         CALL rs_grid_release_descriptor(ewald_pw%rs_desc)
420         ewald_pw%rs_desc => rs_desc
421      ENDIF
422      IF (PRESENT(dg)) THEN
423         CALL dg_retain(dg)
424         CALL dg_release(ewald_pw%dg)
425         ewald_pw%dg => dg
426      ENDIF
427      IF (PRESENT(poisson_env)) THEN
428         IF (ASSOCIATED(poisson_env)) &
429            CALL pw_poisson_retain(poisson_env)
430         CALL pw_poisson_release(ewald_pw%poisson_env)
431         ewald_pw%poisson_env => poisson_env
432      END IF
433
434   END SUBROUTINE ewald_pw_set
435
436END MODULE ewald_pw_types
437