1!--------------------------------------------------------------------------------------------------!
2!   CP2K: A general program to perform molecular dynamics simulations                              !
3!   Copyright (C) 2000 - 2020  CP2K developers group                                               !
4!--------------------------------------------------------------------------------------------------!
5! **************************************************************************************************
6MODULE qs_local_rho_types
7
8   USE kinds,                           ONLY: dp
9   USE mathconstants,                   ONLY: fourpi,&
10                                              pi
11   USE memory_utilities,                ONLY: reallocate
12   USE qs_grid_atom,                    ONLY: grid_atom_type
13   USE qs_harmonics_atom,               ONLY: harmonics_atom_type
14   USE qs_rho0_types,                   ONLY: deallocate_rho0_atom,&
15                                              deallocate_rho0_mpole,&
16                                              rho0_atom_type,&
17                                              rho0_mpole_type
18   USE qs_rho_atom_types,               ONLY: deallocate_rho_atom_set,&
19                                              rho_atom_type
20#include "./base/base_uses.f90"
21
22   IMPLICIT NONE
23
24   PRIVATE
25
26! *** Global parameters (only in this module)
27
28   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_local_rho_types'
29
30! *** Define rhoz and local_rho types ***
31
32! **************************************************************************************************
33   TYPE rhoz_type
34      REAL(dp)                             ::  one_atom
35      REAL(dp), DIMENSION(:), POINTER      ::  r_coef
36      REAL(dp), DIMENSION(:), POINTER      ::  dr_coef
37      REAL(dp), DIMENSION(:), POINTER      ::  vr_coef
38   END TYPE rhoz_type
39
40! **************************************************************************************************
41   TYPE local_rho_type
42      TYPE(rho_atom_type), DIMENSION(:), POINTER            :: rho_atom_set
43      TYPE(rho0_mpole_type), POINTER                        :: rho0_mpole
44      TYPE(rho0_atom_type), DIMENSION(:), POINTER           :: rho0_atom_set
45      TYPE(rhoz_type), DIMENSION(:), POINTER               :: rhoz_set
46      REAL(dp)                                              :: rhoz_tot
47   END TYPE local_rho_type
48
49! Public Types
50   PUBLIC :: local_rho_type, rhoz_type
51
52! Public Subroutine
53   PUBLIC :: allocate_rhoz, calculate_rhoz, &
54             get_local_rho, local_rho_set_create, &
55             local_rho_set_release, set_local_rho
56
57CONTAINS
58
59! **************************************************************************************************
60!> \brief ...
61!> \param rhoz_set ...
62!> \param nkind ...
63! **************************************************************************************************
64   SUBROUTINE allocate_rhoz(rhoz_set, nkind)
65
66      TYPE(rhoz_type), DIMENSION(:), POINTER             :: rhoz_set
67      INTEGER                                            :: nkind
68
69      CHARACTER(len=*), PARAMETER :: routineN = 'allocate_rhoz', routineP = moduleN//':'//routineN
70
71      INTEGER                                            :: ikind
72
73      IF (ASSOCIATED(rhoz_set)) THEN
74         CALL deallocate_rhoz(rhoz_set)
75      END IF
76
77      ALLOCATE (rhoz_set(nkind))
78
79      DO ikind = 1, nkind
80         NULLIFY (rhoz_set(ikind)%r_coef)
81         NULLIFY (rhoz_set(ikind)%dr_coef)
82         NULLIFY (rhoz_set(ikind)%vr_coef)
83      ENDDO
84
85   END SUBROUTINE allocate_rhoz
86
87! **************************************************************************************************
88!> \brief ...
89!> \param rhoz ...
90!> \param grid_atom ...
91!> \param alpha ...
92!> \param zeff ...
93!> \param natom ...
94!> \param rhoz_tot ...
95!> \param harmonics ...
96! **************************************************************************************************
97   SUBROUTINE calculate_rhoz(rhoz, grid_atom, alpha, zeff, natom, rhoz_tot, harmonics)
98
99      TYPE(rhoz_type)                                    :: rhoz
100      TYPE(grid_atom_type)                               :: grid_atom
101      REAL(dp), INTENT(IN)                               :: alpha
102      REAL(dp)                                           :: zeff
103      INTEGER                                            :: natom
104      REAL(dp), INTENT(INOUT)                            :: rhoz_tot
105      TYPE(harmonics_atom_type)                          :: harmonics
106
107      CHARACTER(len=*), PARAMETER :: routineN = 'calculate_rhoz', routineP = moduleN//':'//routineN
108
109      INTEGER                                            :: ir, na, nr
110      REAL(dp)                                           :: c1, c2, c3, prefactor1, prefactor2, &
111                                                            prefactor3, sum
112
113      nr = grid_atom%nr
114      na = grid_atom%ng_sphere
115      CALL reallocate(rhoz%r_coef, 1, nr)
116      CALL reallocate(rhoz%dr_coef, 1, nr)
117      CALL reallocate(rhoz%vr_coef, 1, nr)
118
119      c1 = alpha/pi
120      c2 = c1*c1*c1*fourpi
121      c3 = SQRT(alpha)
122      prefactor1 = zeff*SQRT(c2)
123      prefactor2 = -2.0_dp*alpha
124      prefactor3 = -zeff*SQRT(fourpi)
125
126      sum = 0.0_dp
127      DO ir = 1, nr
128         c1 = -alpha*grid_atom%rad2(ir)
129         rhoz%r_coef(ir) = -EXP(c1)*prefactor1
130         IF (ABS(rhoz%r_coef(ir)) < 1.0E-30_dp) THEN
131            rhoz%r_coef(ir) = 0.0_dp
132            rhoz%dr_coef(ir) = 0.0_dp
133         ELSE
134            rhoz%dr_coef(ir) = prefactor2*rhoz%r_coef(ir)
135         END IF
136         rhoz%vr_coef(ir) = prefactor3*erf(grid_atom%rad(ir)*c3)/grid_atom%rad(ir)
137         sum = sum + rhoz%r_coef(ir)*grid_atom%wr(ir)
138      END DO
139      rhoz%one_atom = sum*harmonics%slm_int(1)
140      rhoz_tot = rhoz_tot + natom*rhoz%one_atom
141
142   END SUBROUTINE calculate_rhoz
143
144! **************************************************************************************************
145!> \brief ...
146!> \param rhoz_set ...
147! **************************************************************************************************
148   SUBROUTINE deallocate_rhoz(rhoz_set)
149
150      TYPE(rhoz_type), DIMENSION(:), POINTER             :: rhoz_set
151
152      CHARACTER(len=*), PARAMETER :: routineN = 'deallocate_rhoz', &
153         routineP = moduleN//':'//routineN
154
155      INTEGER                                            :: ikind, nkind
156
157      nkind = SIZE(rhoz_set)
158
159      DO ikind = 1, nkind
160         DEALLOCATE (rhoz_set(ikind)%r_coef)
161         DEALLOCATE (rhoz_set(ikind)%dr_coef)
162         DEALLOCATE (rhoz_set(ikind)%vr_coef)
163      END DO
164
165      DEALLOCATE (rhoz_set)
166
167   END SUBROUTINE deallocate_rhoz
168
169! **************************************************************************************************
170!> \brief ...
171!> \param local_rho_set ...
172!> \param rho_atom_set ...
173!> \param rho0_atom_set ...
174!> \param rho0_mpole ...
175!> \param rhoz_set ...
176! **************************************************************************************************
177   SUBROUTINE get_local_rho(local_rho_set, rho_atom_set, rho0_atom_set, rho0_mpole, rhoz_set)
178
179      TYPE(local_rho_type), POINTER                      :: local_rho_set
180      TYPE(rho_atom_type), DIMENSION(:), OPTIONAL, &
181         POINTER                                         :: rho_atom_set
182      TYPE(rho0_atom_type), DIMENSION(:), OPTIONAL, &
183         POINTER                                         :: rho0_atom_set
184      TYPE(rho0_mpole_type), OPTIONAL, POINTER           :: rho0_mpole
185      TYPE(rhoz_type), DIMENSION(:), OPTIONAL, POINTER   :: rhoz_set
186
187      CHARACTER(len=*), PARAMETER :: routineN = 'get_local_rho', routineP = moduleN//':'//routineN
188
189      IF (PRESENT(rho_atom_set)) rho_atom_set => local_rho_set%rho_atom_set
190      IF (PRESENT(rho0_atom_set)) rho0_atom_set => local_rho_set%rho0_atom_set
191      IF (PRESENT(rho0_mpole)) rho0_mpole => local_rho_set%rho0_mpole
192      IF (PRESENT(rhoz_set)) rhoz_set => local_rho_set%rhoz_set
193
194   END SUBROUTINE get_local_rho
195
196! **************************************************************************************************
197!> \brief ...
198!> \param local_rho_set ...
199! **************************************************************************************************
200   SUBROUTINE local_rho_set_create(local_rho_set)
201
202      TYPE(local_rho_type), POINTER                      :: local_rho_set
203
204      CHARACTER(len=*), PARAMETER :: routineN = 'local_rho_set_create', &
205         routineP = moduleN//':'//routineN
206
207      ALLOCATE (local_rho_set)
208
209      NULLIFY (local_rho_set%rho_atom_set)
210      NULLIFY (local_rho_set%rho0_atom_set)
211      NULLIFY (local_rho_set%rho0_mpole)
212      NULLIFY (local_rho_set%rhoz_set)
213
214   END SUBROUTINE local_rho_set_create
215
216! **************************************************************************************************
217!> \brief ...
218!> \param local_rho_set ...
219! **************************************************************************************************
220   SUBROUTINE local_rho_set_release(local_rho_set)
221
222      TYPE(local_rho_type), POINTER                      :: local_rho_set
223
224      CHARACTER(len=*), PARAMETER :: routineN = 'local_rho_set_release', &
225         routineP = moduleN//':'//routineN
226
227      IF (ASSOCIATED(local_rho_set)) THEN
228         IF (ASSOCIATED(local_rho_set%rho_atom_set)) THEN
229            CALL deallocate_rho_atom_set(local_rho_set%rho_atom_set)
230         END IF
231
232         IF (ASSOCIATED(local_rho_set%rho0_atom_set)) THEN
233            CALL deallocate_rho0_atom(local_rho_set%rho0_atom_set)
234         END IF
235
236         IF (ASSOCIATED(local_rho_set%rho0_mpole)) THEN
237            CALL deallocate_rho0_mpole(local_rho_set%rho0_mpole)
238         END IF
239
240         IF (ASSOCIATED(local_rho_set%rhoz_set)) THEN
241            CALL deallocate_rhoz(local_rho_set%rhoz_set)
242         ENDIF
243
244         DEALLOCATE (local_rho_set)
245      END IF
246
247   END SUBROUTINE local_rho_set_release
248
249! **************************************************************************************************
250!> \brief ...
251!> \param local_rho_set ...
252!> \param rho_atom_set ...
253!> \param rho0_atom_set ...
254!> \param rho0_mpole ...
255!> \param rhoz_set ...
256! **************************************************************************************************
257   SUBROUTINE set_local_rho(local_rho_set, rho_atom_set, rho0_atom_set, rho0_mpole, &
258                            rhoz_set)
259
260      TYPE(local_rho_type), POINTER                      :: local_rho_set
261      TYPE(rho_atom_type), DIMENSION(:), OPTIONAL, &
262         POINTER                                         :: rho_atom_set
263      TYPE(rho0_atom_type), DIMENSION(:), OPTIONAL, &
264         POINTER                                         :: rho0_atom_set
265      TYPE(rho0_mpole_type), OPTIONAL, POINTER           :: rho0_mpole
266      TYPE(rhoz_type), DIMENSION(:), OPTIONAL, POINTER   :: rhoz_set
267
268      CHARACTER(len=*), PARAMETER :: routineN = 'set_local_rho', routineP = moduleN//':'//routineN
269
270      IF (PRESENT(rho_atom_set)) THEN
271         IF (ASSOCIATED(local_rho_set%rho_atom_set)) THEN
272            CALL deallocate_rho_atom_set(local_rho_set%rho_atom_set)
273         ENDIF
274         local_rho_set%rho_atom_set => rho_atom_set
275      END IF
276
277      IF (PRESENT(rho0_atom_set)) THEN
278         IF (ASSOCIATED(local_rho_set%rho0_atom_set)) THEN
279            CALL deallocate_rho0_atom(local_rho_set%rho0_atom_set)
280         ENDIF
281         local_rho_set%rho0_atom_set => rho0_atom_set
282      END IF
283
284      IF (PRESENT(rho0_mpole)) THEN
285         IF (ASSOCIATED(local_rho_set%rho0_mpole)) THEN
286            CALL deallocate_rho0_mpole(local_rho_set%rho0_mpole)
287         ENDIF
288         local_rho_set%rho0_mpole => rho0_mpole
289      END IF
290
291      IF (PRESENT(rhoz_set)) THEN
292         IF (ASSOCIATED(local_rho_set%rhoz_set)) THEN
293            CALL deallocate_rhoz(local_rho_set%rhoz_set)
294         ENDIF
295         local_rho_set%rhoz_set => rhoz_set
296      END IF
297
298   END SUBROUTINE set_local_rho
299
300END MODULE qs_local_rho_types
301
302