1!--------------------------------------------------------------------------------------------------!
2!   CP2K: A general program to perform molecular dynamics simulations                              !
3!   Copyright (C) 2000 - 2020  CP2K developers group                                               !
4!--------------------------------------------------------------------------------------------------!
5
6! **************************************************************************************************
7!> \brief contains information regarding the decoupling/recoupling method of Bloechl
8!> \author Teodoro Laino
9! **************************************************************************************************
10MODULE cp_ddapc_types
11   USE cell_methods,                    ONLY: read_cell
12   USE cell_types,                      ONLY: cell_release,&
13                                              cell_type
14   USE cp_ddapc_methods,                ONLY: ddapc_eval_AmI,&
15                                              ddapc_eval_gfunc,&
16                                              ewald_ddapc_pot,&
17                                              solvation_ddapc_pot
18   USE cp_log_handling,                 ONLY: cp_get_default_logger,&
19                                              cp_logger_get_default_io_unit,&
20                                              cp_logger_type
21   USE cp_output_handling,              ONLY: cp_printkey_is_on
22   USE cp_para_types,                   ONLY: cp_para_env_type
23   USE ewald_spline_util,               ONLY: Setup_Ewald_Spline
24   USE input_section_types,             ONLY: section_vals_get,&
25                                              section_vals_get_subs_vals,&
26                                              section_vals_type,&
27                                              section_vals_val_get
28   USE kinds,                           ONLY: dp
29   USE mathconstants,                   ONLY: pi
30   USE particle_types,                  ONLY: particle_type
31   USE pw_grid_types,                   ONLY: pw_grid_type
32   USE pw_grids,                        ONLY: pw_grid_release
33   USE pw_poisson_types,                ONLY: pw_poisson_multipole
34   USE pw_pool_types,                   ONLY: pw_pool_give_back_pw,&
35                                              pw_pool_release,&
36                                              pw_pool_type
37   USE pw_types,                        ONLY: pw_type
38#include "./base/base_uses.f90"
39
40   IMPLICIT NONE
41   PRIVATE
42   LOGICAL, PRIVATE, PARAMETER          :: debug_this_module = .TRUE.
43   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_ddapc_types'
44   INTEGER, PRIVATE, SAVE               :: last_cp_ddapc_id = 0
45
46   PUBLIC :: cp_ddapc_type, cp_ddapc_create, cp_ddapc_retain, cp_ddapc_release
47   PUBLIC :: cp_ddapc_ewald_type, cp_ddapc_ewald_create, cp_ddapc_ewald_release
48
49! **************************************************************************************************
50!> \author Teodoro Laino
51! **************************************************************************************************
52   TYPE cp_ddapc_type
53      INTEGER :: ref_count, id_nr
54      REAL(KIND=dp) :: c0
55      REAL(KIND=dp), DIMENSION(:, :), POINTER     :: AmI
56      REAL(KIND=dp), DIMENSION(:, :), POINTER     :: Md ! decoupling
57      REAL(KIND=dp), DIMENSION(:, :), POINTER     :: Mr ! recoupling
58      REAL(KIND=dp), DIMENSION(:, :), POINTER     :: Mt ! decoupling+recoupling
59      REAL(KIND=dp), DIMENSION(:, :), POINTER     :: Ms ! solvation
60      REAL(KIND=dp), POINTER, DIMENSION(:, :)  :: gfunc
61      REAL(KIND=dp), POINTER, DIMENSION(:)    :: w
62   END TYPE cp_ddapc_type
63
64! **************************************************************************************************
65   TYPE cp_ddapc_ewald_type
66      LOGICAL                                    :: do_decoupling
67      LOGICAL                                    :: do_qmmm_periodic_decpl
68      LOGICAL                                    :: do_solvation
69      LOGICAL                                    :: do_property
70      LOGICAL                                    :: do_restraint
71      TYPE(section_vals_type), POINTER :: ewald_section
72      TYPE(pw_pool_type), POINTER :: pw_pool_qm, pw_pool_mm
73      TYPE(pw_grid_type), POINTER :: pw_grid_qm, pw_grid_mm
74      TYPE(pw_type), POINTER :: coeff_qm, coeff_mm
75   END TYPE cp_ddapc_ewald_type
76
77CONTAINS
78
79! **************************************************************************************************
80!> \brief ...
81!> \param cp_para_env ...
82!> \param cp_ddapc_env ...
83!> \param cp_ddapc_ewald ...
84!> \param particle_set ...
85!> \param radii ...
86!> \param cell ...
87!> \param super_cell ...
88!> \param rho_tot_g ...
89!> \param gcut ...
90!> \param iw2 ...
91!> \param Vol ...
92!> \param force_env_section ...
93!> \author Tedoro Laino
94!> \note NB receive cp_para_env to pass down to parallelized ewald_ddapc_pot()
95! **************************************************************************************************
96   SUBROUTINE cp_ddapc_create(cp_para_env, cp_ddapc_env, cp_ddapc_ewald, &
97                              particle_set, radii, cell, super_cell, rho_tot_g, gcut, iw2, Vol, &
98                              force_env_section)
99      TYPE(cp_para_env_type), POINTER                    :: cp_para_env
100      TYPE(cp_ddapc_type), POINTER                       :: cp_ddapc_env
101      TYPE(cp_ddapc_ewald_type), POINTER                 :: cp_ddapc_ewald
102      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
103      REAL(kind=dp), DIMENSION(:), POINTER               :: radii
104      TYPE(cell_type), POINTER                           :: cell, super_cell
105      TYPE(pw_type), POINTER                             :: rho_tot_g
106      REAL(KIND=dp), INTENT(IN)                          :: gcut
107      INTEGER, INTENT(IN)                                :: iw2
108      REAL(KIND=dp), INTENT(IN)                          :: Vol
109      TYPE(section_vals_type), POINTER                   :: force_env_section
110
111      CHARACTER(len=*), PARAMETER :: routineN = 'cp_ddapc_create', &
112         routineP = moduleN//':'//routineN
113
114      INTEGER                                            :: handle
115      TYPE(section_vals_type), POINTER                   :: param_section, solvation_section
116
117      CALL timeset(routineN, handle)
118      ALLOCATE (cp_ddapc_env)
119      cp_ddapc_env%ref_count = 1
120      last_cp_ddapc_id = last_cp_ddapc_id + 1
121      cp_ddapc_env%id_nr = last_cp_ddapc_id
122      NULLIFY (cp_ddapc_env%AmI, &
123               cp_ddapc_env%Md, &
124               cp_ddapc_env%Mt, &
125               cp_ddapc_env%Mr, &
126               cp_ddapc_env%Ms, &
127               cp_ddapc_env%gfunc, &
128               cp_ddapc_env%w)
129      ! Evaluates gfunc and AmI
130      CALL ddapc_eval_gfunc(cp_ddapc_env%gfunc, cp_ddapc_env%w, gcut, rho_tot_g, radii)
131      CALL ddapc_eval_AmI(cp_ddapc_env%AmI, &
132                          cp_ddapc_env%c0, &
133                          cp_ddapc_env%gfunc, &
134                          cp_ddapc_env%w, &
135                          particle_set, &
136                          gcut, &
137                          rho_tot_g, &
138                          radii, &
139                          iw2, &
140                          Vol)
141      IF (cp_ddapc_ewald%do_qmmm_periodic_decpl .OR. &
142          cp_ddapc_ewald%do_decoupling) THEN
143         !
144         ! Evaluate the matrix for the Classical contribution to the coupling/decoupling scheme
145         !
146         param_section => cp_ddapc_ewald%ewald_section
147         !NB parallelized ewald_ddapc_pot() needs cp_para_env
148         CALL ewald_ddapc_pot(cp_para_env, cp_ddapc_ewald%coeff_qm, &
149                              1.0_dp, &
150                              cell, &
151                              param_section, &
152                              particle_set, &
153                              cp_ddapc_env%Md, &
154                              radii)
155         IF (cp_ddapc_ewald%do_qmmm_periodic_decpl .OR. cp_ddapc_ewald%do_decoupling) THEN
156            ALLOCATE (cp_ddapc_env%Mt(SIZE(cp_ddapc_env%Md, 1), SIZE(cp_ddapc_env%Md, 2)))
157            IF (cp_ddapc_ewald%do_decoupling) THEN
158               ! Just decoupling
159               cp_ddapc_env%Mt = cp_ddapc_env%Md
160            ELSE
161               ! QMMM periodic calculation
162               !NB parallelized ewald_ddapc_pot() needs cp_para_env
163               CALL ewald_ddapc_pot(cp_para_env, cp_ddapc_ewald%coeff_mm, -1.0_dp, super_cell, param_section, &
164                                    particle_set, cp_ddapc_env%Mr, radii)
165               cp_ddapc_env%Mt = cp_ddapc_env%Md + cp_ddapc_env%Mr
166            END IF
167         END IF
168      END IF
169      IF (cp_ddapc_ewald%do_solvation) THEN
170         ! Spherical Solvation model
171         solvation_section => section_vals_get_subs_vals(force_env_section, "DFT%SCRF")
172         CALL solvation_ddapc_pot(solvation_section, &
173                                  particle_set, cp_ddapc_env%Ms, radii)
174      END IF
175      CALL timestop(handle)
176   END SUBROUTINE cp_ddapc_create
177
178! **************************************************************************************************
179!> \brief ...
180!> \param cp_ddapc_env ...
181!> \par History
182!>      none
183!> \author Teodoro Laino - [tlaino]
184! **************************************************************************************************
185   SUBROUTINE cp_ddapc_retain(cp_ddapc_env)
186      TYPE(cp_ddapc_type), POINTER                       :: cp_ddapc_env
187
188      CHARACTER(len=*), PARAMETER :: routineN = 'cp_ddapc_retain', &
189         routineP = moduleN//':'//routineN
190
191      CPASSERT(ASSOCIATED(cp_ddapc_env))
192      CPASSERT(cp_ddapc_env%ref_count > 0)
193      cp_ddapc_env%ref_count = cp_ddapc_env%ref_count + 1
194   END SUBROUTINE cp_ddapc_retain
195
196! **************************************************************************************************
197!> \brief ...
198!> \param cp_ddapc_env ...
199!> \par History
200!>      none
201!> \author Teodoro Laino - [tlaino]
202! **************************************************************************************************
203   SUBROUTINE cp_ddapc_release(cp_ddapc_env)
204      TYPE(cp_ddapc_type), POINTER                       :: cp_ddapc_env
205
206      CHARACTER(len=*), PARAMETER :: routineN = 'cp_ddapc_release', &
207         routineP = moduleN//':'//routineN
208
209      IF (ASSOCIATED(cp_ddapc_env)) THEN
210         CPASSERT(cp_ddapc_env%ref_count > 0)
211         cp_ddapc_env%ref_count = cp_ddapc_env%ref_count - 1
212         IF (cp_ddapc_env%ref_count == 0) THEN
213            IF (ASSOCIATED(cp_ddapc_env%AmI)) THEN
214               DEALLOCATE (cp_ddapc_env%AmI)
215            END IF
216            IF (ASSOCIATED(cp_ddapc_env%Mt)) THEN
217               DEALLOCATE (cp_ddapc_env%Mt)
218            END IF
219            IF (ASSOCIATED(cp_ddapc_env%Md)) THEN
220               DEALLOCATE (cp_ddapc_env%Md)
221            END IF
222            IF (ASSOCIATED(cp_ddapc_env%Mr)) THEN
223               DEALLOCATE (cp_ddapc_env%Mr)
224            END IF
225            IF (ASSOCIATED(cp_ddapc_env%Ms)) THEN
226               DEALLOCATE (cp_ddapc_env%Ms)
227            END IF
228            IF (ASSOCIATED(cp_ddapc_env%gfunc)) THEN
229               DEALLOCATE (cp_ddapc_env%gfunc)
230            END IF
231            IF (ASSOCIATED(cp_ddapc_env%w)) THEN
232               DEALLOCATE (cp_ddapc_env%w)
233            END IF
234            DEALLOCATE (cp_ddapc_env)
235         END IF
236      END IF
237   END SUBROUTINE cp_ddapc_release
238
239! **************************************************************************************************
240!> \brief ...
241!> \param cp_ddapc_ewald ...
242!> \param qmmm_decoupl ...
243!> \param qm_cell ...
244!> \param force_env_section ...
245!> \param subsys_section ...
246!> \param para_env ...
247!> \par History
248!>      none
249!> \author Teodoro Laino - [tlaino]
250! **************************************************************************************************
251   SUBROUTINE cp_ddapc_ewald_create(cp_ddapc_ewald, qmmm_decoupl, qm_cell, &
252                                    force_env_section, subsys_section, para_env)
253      TYPE(cp_ddapc_ewald_type), POINTER                 :: cp_ddapc_ewald
254      LOGICAL, INTENT(IN)                                :: qmmm_decoupl
255      TYPE(cell_type), POINTER                           :: qm_cell
256      TYPE(section_vals_type), POINTER                   :: force_env_section, subsys_section
257      TYPE(cp_para_env_type), POINTER                    :: para_env
258
259      CHARACTER(len=*), PARAMETER :: routineN = 'cp_ddapc_ewald_create', &
260         routineP = moduleN//':'//routineN
261
262      INTEGER                                            :: my_val, npts(3), output_unit
263      INTEGER, DIMENSION(:), POINTER                     :: ngrids
264      LOGICAL                                            :: analyt, decoupling, &
265                                                            do_qmmm_periodic_decpl, do_restraint, &
266                                                            do_restraintB, do_solvation
267      REAL(KIND=dp)                                      :: hmat(3, 3)
268      REAL(KIND=dp), DIMENSION(:), POINTER               :: gx, gy, gz, LG
269      TYPE(cell_type), POINTER                           :: dummy_cell, mm_cell
270      TYPE(cp_logger_type), POINTER                      :: logger
271      TYPE(section_vals_type), POINTER :: cell_section, grid_print_section, multipole_section, &
272         poisson_section, printC_section, qmmm_per_section, restraint_section, restraint_sectionB, &
273         solvation_section
274
275      logger => cp_get_default_logger()
276      output_unit = cp_logger_get_default_io_unit(logger)
277      CPASSERT(.NOT. ASSOCIATED(cp_ddapc_ewald))
278      ALLOCATE (cp_ddapc_ewald)
279      NULLIFY (cp_ddapc_ewald%pw_grid_mm, &
280               cp_ddapc_ewald%pw_grid_qm, &
281               cp_ddapc_ewald%ewald_section, &
282               cp_ddapc_ewald%pw_pool_mm, &
283               cp_ddapc_ewald%pw_pool_qm, &
284               cp_ddapc_ewald%coeff_mm, &
285               cp_ddapc_ewald%coeff_qm)
286      NULLIFY (multipole_section)
287
288      poisson_section => section_vals_get_subs_vals(force_env_section, "DFT%POISSON")
289      solvation_section => section_vals_get_subs_vals(force_env_section, "DFT%SCRF")
290      qmmm_per_section => section_vals_get_subs_vals(force_env_section, "QMMM%PERIODIC")
291      printC_section => section_vals_get_subs_vals(force_env_section, "PROPERTIES%FIT_CHARGE")
292      restraint_section => section_vals_get_subs_vals(force_env_section, "DFT%QS%DDAPC_RESTRAINT")
293      restraint_sectionB => section_vals_get_subs_vals(force_env_section, &
294                                                       "PROPERTIES%ET_COUPLING%DDAPC_RESTRAINT_A")
295      CALL section_vals_get(solvation_section, explicit=do_solvation)
296      CALL section_vals_get(poisson_section, explicit=decoupling)
297      CALL section_vals_get(restraint_section, explicit=do_restraint)
298      CALL section_vals_get(restraint_sectionB, explicit=do_restraintB)
299      do_qmmm_periodic_decpl = qmmm_decoupl
300      cp_ddapc_ewald%do_solvation = do_solvation
301      cp_ddapc_ewald%do_qmmm_periodic_decpl = do_qmmm_periodic_decpl
302      cp_ddapc_ewald%do_property = cp_printkey_is_on(logger%iter_info, printC_section)
303      cp_ddapc_ewald%do_restraint = do_restraint .OR. do_restraintB
304      ! Determining the tasks and further check
305      IF (do_qmmm_periodic_decpl .AND. decoupling) THEN
306         ! Check than an additional POISSON section has not been defined. In case write a warning
307         IF (output_unit > 0) &
308            WRITE (output_unit, '(T2,"WARNING",A)') &
309            "A calculation with the QMMM periodic model has been requested.", &
310            "The explicit POISSON section in DFT section will be IGNORED.", &
311            "QM Electrostatic controlled only by the PERIODIC section in QMMM section"
312         decoupling = .FALSE.
313      END IF
314      IF (decoupling) THEN
315         ! Simple decoupling technique
316         CALL section_vals_val_get(poisson_section, "POISSON_SOLVER", i_val=my_val)
317         SELECT CASE (my_val)
318         CASE (pw_poisson_multipole)
319            multipole_section => section_vals_get_subs_vals(poisson_section, "MULTIPOLE")
320         CASE DEFAULT
321            decoupling = .FALSE.
322         END SELECT
323      END IF
324      cp_ddapc_ewald%do_decoupling = decoupling
325      IF (cp_ddapc_ewald%do_qmmm_periodic_decpl) THEN
326         ! QMMM periodic
327         multipole_section => section_vals_get_subs_vals(qmmm_per_section, "MULTIPOLE")
328      END IF
329      cp_ddapc_ewald%ewald_section => multipole_section
330      IF (cp_ddapc_ewald%do_decoupling .OR. cp_ddapc_ewald%do_qmmm_periodic_decpl) THEN
331         ! Do we do the calculation analytically or interpolating the g-space factor?
332         CALL section_vals_val_get(multipole_section, "ANALYTICAL_GTERM", l_val=analyt)
333         IF (.NOT. analyt) THEN
334            CALL section_vals_val_get(multipole_section, "ngrids", i_vals=ngrids)
335            npts = ngrids
336
337            NULLIFY (LG, gx, gy, gz)
338            hmat = qm_cell%hmat
339            CALL eval_lg(multipole_section, hmat, qm_cell%deth, LG, gx, gy, gz)
340            grid_print_section => section_vals_get_subs_vals(force_env_section, "PRINT%GRID_INFORMATION")
341            CALL Setup_Ewald_Spline(pw_grid=cp_ddapc_ewald%pw_grid_qm, pw_pool=cp_ddapc_ewald%pw_pool_qm, &
342                                    coeff=cp_ddapc_ewald%coeff_qm, LG=LG, gx=gx, gy=gy, gz=gz, hmat=hmat, npts=npts, &
343                                    param_section=multipole_section, tag="ddapc", &
344                                    para_env=para_env, print_section=grid_print_section)
345            DEALLOCATE (LG)
346            DEALLOCATE (gx)
347            DEALLOCATE (gy)
348            DEALLOCATE (gz)
349            IF (cp_ddapc_ewald%do_qmmm_periodic_decpl) THEN
350               NULLIFY (mm_cell, dummy_cell)
351               cell_section => section_vals_get_subs_vals(subsys_section, "CELL")
352               CALL read_cell(mm_cell, dummy_cell, cell_section=cell_section, para_env=para_env)
353               hmat = mm_cell%hmat
354               CALL eval_lg(multipole_section, hmat, mm_cell%deth, LG, gx, gy, gz)
355               grid_print_section => section_vals_get_subs_vals(force_env_section, "PRINT%GRID_INFORMATION")
356               CALL Setup_Ewald_Spline(pw_grid=cp_ddapc_ewald%pw_grid_mm, pw_pool=cp_ddapc_ewald%pw_pool_mm, &
357                                       coeff=cp_ddapc_ewald%coeff_mm, LG=LG, gx=gx, gy=gy, gz=gz, hmat=hmat, npts=npts, &
358                                       param_section=multipole_section, tag="ddapc", para_env=para_env, &
359                                       print_section=grid_print_section)
360               DEALLOCATE (LG)
361               DEALLOCATE (gx)
362               DEALLOCATE (gy)
363               DEALLOCATE (gz)
364               CALL cell_release(dummy_cell)
365               CALL cell_release(mm_cell)
366            END IF
367         END IF
368      END IF
369   END SUBROUTINE cp_ddapc_ewald_create
370
371! **************************************************************************************************
372!> \brief ...
373!> \param multipole_section ...
374!> \param hmat ...
375!> \param deth ...
376!> \param LG ...
377!> \param gx ...
378!> \param gy ...
379!> \param gz ...
380!> \par History
381!>      none
382!> \author Teodoro Laino - [tlaino]
383! **************************************************************************************************
384   SUBROUTINE eval_lg(multipole_section, hmat, deth, LG, gx, gy, gz)
385      TYPE(section_vals_type), POINTER                   :: multipole_section
386      REAL(KIND=dp), INTENT(IN)                          :: hmat(3, 3), deth
387      REAL(KIND=dp), DIMENSION(:), POINTER               :: LG, gx, gy, gz
388
389      CHARACTER(len=*), PARAMETER :: routineN = 'eval_lg', routineP = moduleN//':'//routineN
390
391      INTEGER                                            :: i, k1, k2, k3, n_rep, ndim, nmax1, &
392                                                            nmax2, nmax3
393      REAL(KIND=dp)                                      :: alpha, eps, fac, fs, fvec(3), galpha, &
394                                                            gsq, gsqi, rcut, tol, tol1
395
396      rcut = MIN(hmat(1, 1), hmat(2, 2), hmat(3, 3))/2.0_dp
397      CALL section_vals_val_get(multipole_section, "RCUT", n_rep_val=n_rep)
398      IF (n_rep == 1) CALL section_vals_val_get(multipole_section, "RCUT", r_val=rcut)
399      CALL section_vals_val_get(multipole_section, "EWALD_PRECISION", r_val=eps)
400      eps = MIN(ABS(eps), 0.5_dp)
401      tol = SQRT(ABS(LOG(eps*rcut)))
402      alpha = SQRT(ABS(LOG(eps*rcut*tol)))/rcut
403      galpha = 1.0_dp/(4.0_dp*alpha*alpha)
404      tol1 = SQRT(-LOG(eps*rcut*(2.0_dp*tol*alpha)**2))
405      nmax1 = NINT(0.25_dp + hmat(1, 1)*alpha*tol1/pi)
406      nmax2 = NINT(0.25_dp + hmat(2, 2)*alpha*tol1/pi)
407      nmax3 = NINT(0.25_dp + hmat(3, 3)*alpha*tol1/pi)
408      fac = 1.e0_dp/deth
409      fvec = 2.0_dp*pi/(/hmat(1, 1), hmat(2, 2), hmat(3, 3)/)
410      ndim = (nmax1 + 1)*(2*nmax2 + 1)*(2*nmax3 + 1) - 1
411      ALLOCATE (LG(ndim))
412      ALLOCATE (gx(ndim))
413      ALLOCATE (gy(ndim))
414      ALLOCATE (gz(ndim))
415
416      i = 0
417      DO k1 = 0, nmax1
418         DO k2 = -nmax2, nmax2
419            DO k3 = -nmax3, nmax3
420               IF (k1 == 0 .AND. k2 == 0 .AND. k3 == 0) CYCLE
421               i = i + 1
422               fs = 2.0_dp; IF (k1 == 0) fs = 1.0_dp
423               gx(i) = fvec(1)*REAL(k1, KIND=dp)
424               gy(i) = fvec(2)*REAL(k2, KIND=dp)
425               gz(i) = fvec(3)*REAL(k3, KIND=dp)
426               gsq = gx(i)*gx(i) + gy(i)*gy(i) + gz(i)*gz(i)
427               gsqi = fs/gsq
428               LG(i) = fac*gsqi*EXP(-galpha*gsq)
429            END DO
430         END DO
431      END DO
432
433   END SUBROUTINE eval_lg
434
435! **************************************************************************************************
436!> \brief ...
437!> \param cp_ddapc_ewald ...
438!> \par History
439!>      none
440!> \author Teodoro Laino - [tlaino]
441! **************************************************************************************************
442   SUBROUTINE cp_ddapc_ewald_release(cp_ddapc_ewald)
443      TYPE(cp_ddapc_ewald_type), POINTER                 :: cp_ddapc_ewald
444
445      CHARACTER(len=*), PARAMETER :: routineN = 'cp_ddapc_ewald_release', &
446         routineP = moduleN//':'//routineN
447
448      IF (ASSOCIATED(cp_ddapc_ewald)) THEN
449         IF (ASSOCIATED(cp_ddapc_ewald%coeff_qm)) THEN
450            CALL pw_pool_give_back_pw(cp_ddapc_ewald%pw_pool_qm, cp_ddapc_ewald%coeff_qm)
451         END IF
452         IF (ASSOCIATED(cp_ddapc_ewald%coeff_mm)) THEN
453            CALL pw_pool_give_back_pw(cp_ddapc_ewald%pw_pool_mm, cp_ddapc_ewald%coeff_mm)
454         END IF
455         IF (ASSOCIATED(cp_ddapc_ewald%pw_pool_qm)) THEN
456            CALL pw_pool_release(cp_ddapc_ewald%pw_pool_qm)
457            CPASSERT(.NOT. ASSOCIATED(cp_ddapc_ewald%pw_pool_qm))
458         END IF
459         IF (ASSOCIATED(cp_ddapc_ewald%pw_pool_mm)) THEN
460            CALL pw_pool_release(cp_ddapc_ewald%pw_pool_mm)
461            CPASSERT(.NOT. ASSOCIATED(cp_ddapc_ewald%pw_pool_mm))
462         END IF
463         IF (ASSOCIATED(cp_ddapc_ewald%pw_grid_qm)) THEN
464            CALL pw_grid_release(cp_ddapc_ewald%pw_grid_qm)
465            CPASSERT(.NOT. ASSOCIATED(cp_ddapc_ewald%pw_grid_qm))
466         END IF
467         IF (ASSOCIATED(cp_ddapc_ewald%pw_grid_mm)) THEN
468            CALL pw_grid_release(cp_ddapc_ewald%pw_grid_mm)
469            CPASSERT(.NOT. ASSOCIATED(cp_ddapc_ewald%pw_grid_mm))
470         END IF
471         DEALLOCATE (cp_ddapc_ewald)
472      END IF
473
474   END SUBROUTINE cp_ddapc_ewald_release
475
476END MODULE cp_ddapc_types
477