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