1!--------------------------------------------------------------------------------------------------!
2!   CP2K: A general program to perform molecular dynamics simulations                              !
3!   Copyright (C) 2000 - 2020  CP2K developers group                                               !
4!--------------------------------------------------------------------------------------------------!
5
6! **************************************************************************************************
7!> \par History
8!>      none
9! **************************************************************************************************
10MODULE dg_rho0_types
11
12   USE kinds,                           ONLY: dp
13   USE pw_grid_types,                   ONLY: pw_grid_type
14   USE pw_methods,                      ONLY: pw_zero
15   USE pw_poisson_types,                ONLY: do_ewald_ewald,&
16                                              do_ewald_none,&
17                                              do_ewald_pme,&
18                                              do_ewald_spme
19   USE pw_types,                        ONLY: REALDATA3D,&
20                                              pw_create,&
21                                              pw_p_type,&
22                                              pw_release
23#include "../base/base_uses.f90"
24
25   IMPLICIT NONE
26
27   PRIVATE
28   PUBLIC:: dg_rho0_type, dg_rho0_init, dg_rho0_set, dg_rho0_get, &
29            dg_rho0_create, dg_rho0_retain, dg_rho0_release
30
31   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'dg_rho0_types'
32   INTEGER, PRIVATE, SAVE :: last_dg_rho0_id_nr = 0
33
34! **************************************************************************************************
35!> \brief   Type for Gaussian Densities
36!>              type = type of gaussian (PME)
37!>              grid = grid number
38!>              gcc = Gaussian contraction coefficient
39!>              zet = Gaussian exponent
40! **************************************************************************************************
41   TYPE dg_rho0_type
42      INTEGER :: ref_count, id_nr
43      INTEGER :: TYPE
44      INTEGER :: grid
45      INTEGER :: kind
46      REAL(KIND=dp) :: cutoff_radius
47      REAL(KIND=dp), DIMENSION(:), POINTER :: gcc
48      REAL(KIND=dp), DIMENSION(:), POINTER :: zet
49      TYPE(pw_p_type) :: density
50   END TYPE dg_rho0_type
51
52CONTAINS
53
54! **************************************************************************************************
55!> \brief   Set the dg_rho0_type
56!> \param dg_rho0 ...
57!> \param TYPE ...
58!> \param grid ...
59!> \param kind ...
60!> \param cutoff_radius ...
61!> \param gcc ...
62!> \param zet ...
63!> \param density ...
64!> \version 1.0
65! **************************************************************************************************
66   SUBROUTINE dg_rho0_set(dg_rho0, TYPE, grid, kind, cutoff_radius, &
67                          gcc, zet, density)
68      INTEGER, OPTIONAL                                  :: TYPE
69      TYPE(dg_rho0_type), POINTER                        :: dg_rho0
70      INTEGER, OPTIONAL                                  :: grid, kind
71      REAL(KIND=dp), OPTIONAL                            :: cutoff_radius
72      REAL(KIND=dp), OPTIONAL, POINTER                   :: gcc(:), zet(:)
73      TYPE(pw_p_type), OPTIONAL                          :: density
74
75      CHARACTER(LEN=*), PARAMETER :: routineN = 'dg_rho0_set', routineP = moduleN//':'//routineN
76
77      IF (PRESENT(grid)) dg_rho0%grid = grid
78      IF (PRESENT(kind)) dg_rho0%kind = kind
79      IF (PRESENT(density)) dg_rho0%density = density
80      IF (PRESENT(gcc)) dg_rho0%gcc => gcc
81      IF (PRESENT(zet)) dg_rho0%zet => zet
82      IF (PRESENT(TYPE)) dg_rho0%type = TYPE
83      IF (PRESENT(cutoff_radius)) dg_rho0%cutoff_radius = cutoff_radius
84
85   END SUBROUTINE dg_rho0_set
86
87! **************************************************************************************************
88!> \brief  Get the dg_rho0_type
89!> \param dg_rho0 ...
90!> \param cutoff_radius ...
91!> \param TYPE ...
92!> \param grid ...
93!> \param kind ...
94!> \param gcc ...
95!> \param zet ...
96!> \param density ...
97!> \version 1.0
98! **************************************************************************************************
99   SUBROUTINE dg_rho0_get(dg_rho0, cutoff_radius, TYPE, grid, kind, gcc, zet, density)
100      INTEGER, OPTIONAL                                  :: TYPE
101      REAL(KIND=dp), OPTIONAL                            :: cutoff_radius
102      TYPE(dg_rho0_type), POINTER                        :: dg_rho0
103      INTEGER, OPTIONAL                                  :: grid, kind
104      REAL(KIND=dp), OPTIONAL, POINTER                   :: gcc(:), zet(:)
105      TYPE(pw_p_type), OPTIONAL, POINTER                 :: density
106
107      CHARACTER(LEN=*), PARAMETER :: routineN = 'dg_rho0_get', routineP = moduleN//':'//routineN
108
109      IF (PRESENT(grid)) grid = dg_rho0%grid
110      IF (PRESENT(kind)) kind = dg_rho0%kind
111      IF (PRESENT(density)) density = dg_rho0%density
112      IF (PRESENT(gcc)) gcc => dg_rho0%gcc
113      IF (PRESENT(zet)) zet => dg_rho0%zet
114      IF (PRESENT(TYPE)) TYPE = dg_rho0%type
115      IF (PRESENT(cutoff_radius)) cutoff_radius = dg_rho0%cutoff_radius
116
117   END SUBROUTINE dg_rho0_get
118
119! **************************************************************************************************
120!> \brief   create the dg_rho0 structure
121!> \param dg_rho0 ...
122!> \version 1.0
123! **************************************************************************************************
124   SUBROUTINE dg_rho0_create(dg_rho0)
125      TYPE(dg_rho0_type), POINTER                        :: dg_rho0
126
127      CHARACTER(LEN=*), PARAMETER :: routineN = 'dg_rho0_create', routineP = moduleN//':'//routineN
128
129      ALLOCATE (dg_rho0)
130      NULLIFY (dg_rho0%gcc)
131      NULLIFY (dg_rho0%zet)
132      dg_rho0%cutoff_radius = 0.0_dp
133      dg_rho0%grid = 0
134      dg_rho0%kind = 0
135      dg_rho0%type = do_ewald_none
136      last_dg_rho0_id_nr = last_dg_rho0_id_nr + 1
137      dg_rho0%id_nr = last_dg_rho0_id_nr
138      dg_rho0%ref_count = 1
139      NULLIFY (dg_rho0%density%pw)
140
141   END SUBROUTINE dg_rho0_create
142
143! **************************************************************************************************
144!> \brief retains the given dg_rho0_type
145!> \param dg_rho0 the dg_rho0_type to retain
146!> \par History
147!>      04.2003 created [fawzi]
148!> \author fawzi
149!> \note
150!>      see doc/ReferenceCounting.html
151! **************************************************************************************************
152   SUBROUTINE dg_rho0_retain(dg_rho0)
153      TYPE(dg_rho0_type), POINTER                        :: dg_rho0
154
155      CHARACTER(len=*), PARAMETER :: routineN = 'dg_rho0_retain', routineP = moduleN//':'//routineN
156
157      CPASSERT(ASSOCIATED(dg_rho0))
158      CPASSERT(dg_rho0%ref_count > 0)
159      dg_rho0%ref_count = dg_rho0%ref_count + 1
160   END SUBROUTINE dg_rho0_retain
161
162! **************************************************************************************************
163!> \brief releases the given dg_rho0_type
164!> \param dg_rho0 the dg_rho0_type to release
165!> \par History
166!>      04.2003 created [fawzi]
167!> \author fawzi
168!> \note
169!>      see doc/ReferenceCounting.html
170! **************************************************************************************************
171   SUBROUTINE dg_rho0_release(dg_rho0)
172      TYPE(dg_rho0_type), POINTER                        :: dg_rho0
173
174      CHARACTER(len=*), PARAMETER :: routineN = 'dg_rho0_release', &
175         routineP = moduleN//':'//routineN
176
177      IF (ASSOCIATED(dg_rho0)) THEN
178         CPASSERT(dg_rho0%ref_count > 0)
179         dg_rho0%ref_count = dg_rho0%ref_count - 1
180         IF (dg_rho0%ref_count == 0) THEN
181            IF (ASSOCIATED(dg_rho0%gcc)) THEN
182               DEALLOCATE (dg_rho0%gcc)
183            END IF
184            IF (ASSOCIATED(dg_rho0%zet)) THEN
185               DEALLOCATE (dg_rho0%zet)
186            END IF
187            CALL pw_release(dg_rho0%density%pw)
188            NULLIFY (dg_rho0%gcc)
189            NULLIFY (dg_rho0%zet)
190            DEALLOCATE (dg_rho0)
191         END IF
192      END IF
193      NULLIFY (dg_rho0)
194   END SUBROUTINE dg_rho0_release
195
196! **************************************************************************************************
197!> \brief ...
198!> \param dg_rho0 ...
199!> \param pw_grid ...
200! **************************************************************************************************
201   SUBROUTINE dg_rho0_init(dg_rho0, pw_grid)
202      TYPE(dg_rho0_type), POINTER                        :: dg_rho0
203      TYPE(pw_grid_type), POINTER                        :: pw_grid
204
205      CHARACTER(LEN=*), PARAMETER :: routineN = 'dg_rho0_init', routineP = moduleN//':'//routineN
206
207      CALL pw_release(dg_rho0%density%pw)
208      SELECT CASE (dg_rho0%type)
209      CASE (do_ewald_ewald)
210         CALL pw_create(dg_rho0%density%pw, pw_grid, REALDATA3D)
211         CALL dg_rho0_pme_gauss(dg_rho0%density, dg_rho0%zet(1))
212      CASE (do_ewald_pme)
213         CALL pw_create(dg_rho0%density%pw, pw_grid, REALDATA3D)
214         CALL dg_rho0_pme_gauss(dg_rho0%density, dg_rho0%zet(1))
215      CASE (do_ewald_spme)
216         CPABORT('SPME type not implemented')
217      END SELECT
218
219   END SUBROUTINE dg_rho0_init
220
221! **************************************************************************************************
222!> \brief ...
223!> \param dg_rho0 ...
224!> \param alpha ...
225! **************************************************************************************************
226   SUBROUTINE dg_rho0_pme_gauss(dg_rho0, alpha)
227
228      TYPE(pw_p_type), INTENT(INOUT)                     :: dg_rho0
229      REAL(KIND=dp), INTENT(IN)                          :: alpha
230
231      CHARACTER(LEN=*), PARAMETER :: routineN = 'dg_rho0_pme_gauss', &
232         routineP = moduleN//':'//routineN
233      INTEGER, PARAMETER                                 :: IMPOSSIBLE = 10000
234
235      INTEGER                                            :: gpt, l0, ln, lp, m0, mn, mp, n0, nn, np
236      INTEGER, DIMENSION(:), POINTER                     :: ghat
237      INTEGER, DIMENSION(:, :), POINTER                  :: bds
238      REAL(KIND=dp)                                      :: const, e_gsq
239      REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: rho0
240      TYPE(pw_grid_type), POINTER                        :: pw_grid
241
242      const = 1.0_dp/(8.0_dp*alpha**2)
243
244      pw_grid => dg_rho0%pw%pw_grid
245      bds => pw_grid%bounds
246
247      IF (-bds(1, 1) == bds(2, 1)) THEN
248         l0 = IMPOSSIBLE
249      ELSE
250         l0 = bds(1, 1)
251      END IF
252
253      IF (-bds(1, 2) == bds(2, 2)) THEN
254         m0 = IMPOSSIBLE
255      ELSE
256         m0 = bds(1, 2)
257      END IF
258
259      IF (-bds(1, 3) == bds(2, 3)) THEN
260         n0 = IMPOSSIBLE
261      ELSE
262         n0 = bds(1, 3)
263      END IF
264
265      CALL pw_zero(dg_rho0%pw)
266
267      rho0 => dg_rho0%pw%cr3d
268
269      DO gpt = 1, pw_grid%ngpts_cut_local
270         ghat => pw_grid%g_hat(:, gpt)
271
272         lp = pw_grid%mapl%pos(ghat(1))
273         ln = pw_grid%mapl%neg(ghat(1))
274         mp = pw_grid%mapm%pos(ghat(2))
275         mn = pw_grid%mapm%neg(ghat(2))
276         np = pw_grid%mapn%pos(ghat(3))
277         nn = pw_grid%mapn%neg(ghat(3))
278
279         e_gsq = EXP(-const*pw_grid%gsq(gpt))/pw_grid%vol
280
281         !*apsi
282         lp = lp + bds(1, 1)
283         mp = mp + bds(1, 2)
284         np = np + bds(1, 3)
285         ln = ln + bds(1, 1)
286         mn = mn + bds(1, 2)
287         nn = nn + bds(1, 3)
288
289         rho0(lp, mp, np) = e_gsq
290         rho0(ln, mn, nn) = e_gsq
291
292         IF (ghat(1) == l0 .OR. ghat(2) == m0 .OR. ghat(3) == n0) THEN
293            rho0(lp, mp, np) = 0.0_dp
294            rho0(ln, mn, nn) = 0.0_dp
295         END IF
296
297      END DO
298
299   END SUBROUTINE dg_rho0_pme_gauss
300
301END MODULE dg_rho0_types
302