Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Data type and methods dealing with PI calcs in staging coordinates
10 : !> \author fawzi
11 : !> \par History
12 : !> 2006-02 created
13 : !> 2006-11 modified so it might actually work [hforbert]
14 : !> 2009-04-07 moved from pint_types module to a separate file [lwalewski]
15 : ! **************************************************************************************************
16 : MODULE pint_staging
17 : USE input_section_types, ONLY: section_vals_type,&
18 : section_vals_val_get
19 : USE kinds, ONLY: dp
20 : USE pint_types, ONLY: staging_env_type
21 : #include "../base/base_uses.f90"
22 :
23 : IMPLICIT NONE
24 : PRIVATE
25 :
26 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
27 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pint_staging'
28 :
29 : PUBLIC :: staging_env_create, staging_release, &
30 : staging_init_masses, &
31 : staging_x2u, staging_u2x, staging_f2uf, &
32 : staging_calc_uf_h
33 :
34 : CONTAINS
35 :
36 : ! ***************************************************************************
37 : !> \brief creates the data needed for a staging transformation
38 : !> \param staging_env ...
39 : !> \param staging_section ...
40 : !> \param p ...
41 : !> \param kT ...
42 : !> \author fawzi
43 : ! **************************************************************************************************
44 0 : SUBROUTINE staging_env_create(staging_env, staging_section, p, kT)
45 : TYPE(staging_env_type), INTENT(OUT) :: staging_env
46 : TYPE(section_vals_type), POINTER :: staging_section
47 : INTEGER, INTENT(in) :: p
48 : REAL(kind=dp), INTENT(in) :: kT
49 :
50 0 : CALL section_vals_val_get(staging_section, "j", i_val=staging_env%j)
51 0 : CALL section_vals_val_get(staging_section, "Q_end", i_val=staging_env%j)
52 0 : staging_env%p = p
53 0 : staging_env%nseg = staging_env%p/staging_env%j
54 :
55 0 : staging_env%w_p = SQRT(REAL(staging_env%p, dp))*kT
56 0 : staging_env%w_j = kT*SQRT(REAL(staging_env%nseg, dp))
57 0 : staging_env%Q_stage = kT/staging_env%w_p**2
58 0 : IF (staging_env%Q_end <= 0._dp) THEN
59 0 : staging_env%Q_end = staging_env%j*staging_env%Q_stage
60 : END IF
61 0 : END SUBROUTINE staging_env_create
62 :
63 : ! ***************************************************************************
64 : !> \brief releases the staging environment, kept for symmetry reasons with staging_env_create
65 : !> \param staging_env the staging_env to release
66 : !> \author Fawzi Mohamed
67 : ! **************************************************************************************************
68 0 : ELEMENTAL SUBROUTINE staging_release(staging_env)
69 : TYPE(staging_env_type), INTENT(IN) :: staging_env
70 :
71 : MARK_USED(staging_env)
72 0 : END SUBROUTINE staging_release
73 :
74 : ! ***************************************************************************
75 : !> \brief initializes the masses and fictitious masses compatibly with the
76 : !> staging information
77 : !> \param staging_env the definition of the staging transformation
78 : !> \param mass *input* the masses of the particles
79 : !> \param mass_beads masses of the beads
80 : !> \param mass_fict the fictitious masses
81 : !> \param Q masses of the nose thermostats
82 : !> \par History
83 : !> 11.2003 created [fawzi]
84 : !> \author Fawzi Mohamed
85 : ! **************************************************************************************************
86 0 : PURE SUBROUTINE staging_init_masses(staging_env, mass, mass_beads, mass_fict, &
87 0 : Q)
88 : TYPE(staging_env_type), INTENT(IN) :: staging_env
89 : REAL(kind=dp), DIMENSION(:), INTENT(in) :: mass
90 : REAL(kind=dp), DIMENSION(:, :), INTENT(out), &
91 : OPTIONAL :: mass_beads, mass_fict
92 : REAL(kind=dp), DIMENSION(:), INTENT(out), OPTIONAL :: Q
93 :
94 : INTEGER :: i, iat, ib, iseg
95 0 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: scal
96 :
97 0 : IF (PRESENT(Q)) THEN
98 0 : Q = staging_env%Q_stage
99 0 : DO i = 1, staging_env%p, staging_env%j
100 0 : Q(i) = staging_env%Q_end
101 : END DO
102 : END IF
103 0 : IF (PRESENT(mass_beads) .OR. PRESENT(mass_fict)) THEN
104 :
105 0 : ALLOCATE (scal(staging_env%p))
106 0 : DO iseg = 1, staging_env%nseg
107 0 : DO i = 1, staging_env%j ! check order!!!
108 0 : scal(staging_env%j*(iseg - 1) + i) = REAL(i, dp)/REAL(MAX(1, i - 1), dp)
109 : END DO
110 : END DO
111 : ! scal=zeros(staging_env%j,Float64)
112 : ! divide(arange(2,staging_env%j+1,typecode=Float64),
113 : ! arange(1,staging_env%j,typecode=Float64),scal[1:])
114 : ! scal[0]=1.
115 : ! scal=outerproduct(ones(staging_env%nseg),scal)
116 :
117 0 : IF (PRESENT(mass_beads)) THEN
118 0 : DO iat = 1, SIZE(mass)
119 0 : DO ib = 1, staging_env%p
120 0 : mass_beads(ib, iat) = scal(ib)*mass(iat)
121 : END DO
122 : END DO
123 : END IF
124 0 : IF (PRESENT(mass_fict)) THEN
125 0 : DO iat = 1, SIZE(mass)
126 0 : DO ib = 1, staging_env%p
127 0 : mass_fict(ib, iat) = scal(ib)*mass(iat)
128 : END DO
129 : END DO
130 : END IF
131 : END IF
132 0 : END SUBROUTINE staging_init_masses
133 :
134 : ! ***************************************************************************
135 : !> \brief Transforms from the x into the u variables using a staging
136 : !> transformation for the positions
137 : !> \param staging_env the environment for the staging transformation
138 : !> \param ux will contain the u variable
139 : !> \param x the positions to transform
140 : !> \author fawzi
141 : ! **************************************************************************************************
142 0 : PURE SUBROUTINE staging_x2u(staging_env, ux, x)
143 : TYPE(staging_env_type), INTENT(IN) :: staging_env
144 : REAL(kind=dp), DIMENSION(:, :), INTENT(out) :: ux
145 : REAL(kind=dp), DIMENSION(:, :), INTENT(in) :: x
146 :
147 : INTEGER :: k, s
148 :
149 0 : ux = x
150 0 : DO s = 0, staging_env%nseg - 1
151 0 : DO k = 2, staging_env%j
152 : ux(staging_env%j*s + k, :) = ux(staging_env%j*s + k, :) &
153 : - ((REAL(k - 1, dp)/REAL(k, dp) &
154 : *x(MODULO((staging_env%j*s + k + 1), staging_env%p), :) + &
155 0 : x(staging_env%j*s + 1, :)/REAL(k, dp)))
156 : END DO
157 : END DO
158 0 : END SUBROUTINE staging_x2u
159 :
160 : ! ***************************************************************************
161 : !> \brief transform from the u variable to the x (back staging transformation
162 : !> for the positions)
163 : !> \param staging_env the environment for the staging transformation
164 : !> \param ux the u variable (positions to be backtransformed)
165 : !> \param x will contain the positions
166 : !> \author fawzi
167 : ! **************************************************************************************************
168 0 : PURE SUBROUTINE staging_u2x(staging_env, ux, x)
169 : TYPE(staging_env_type), INTENT(IN) :: staging_env
170 : REAL(kind=dp), DIMENSION(:, :), INTENT(in) :: ux
171 : REAL(kind=dp), DIMENSION(:, :), INTENT(out) :: x
172 :
173 : INTEGER :: i, ist, j
174 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: iii, jjj
175 : REAL(kind=dp) :: const, const2
176 :
177 0 : j = staging_env%j
178 0 : const = REAL(j - 1, dp)/REAL(j, dp)
179 0 : const2 = 1._dp/REAL(j, dp)
180 0 : ALLOCATE (iii(staging_env%nseg), jjj(staging_env%nseg))
181 0 : DO i = 1, staging_env%nseg
182 0 : iii(i) = staging_env%j*(i - 1) + 1 !first el
183 : END DO
184 0 : DO i = 1, staging_env%nseg - 1
185 0 : jjj(i) = iii(i) + j ! next first el (pbc)
186 : END DO
187 0 : jjj(staging_env%nseg) = 1
188 :
189 0 : x = ux
190 0 : DO i = 1, staging_env%nseg
191 : x(j - 1 + iii(i), :) = x(j - 1 + iii(i), :) + &
192 0 : const*ux(jjj(i), :) + ux(iii(i), :)*const2
193 : END DO
194 0 : DO ist = 1, staging_env%nseg
195 0 : DO i = staging_env%j - 2, 2, -1
196 : x(i + iii(ist), :) = x(i + iii(ist), :) + &
197 : REAL(i - 1, dp)/REAL(i, dp)*x(i + iii(ist) + 1, :) &
198 0 : + ux(iii(ist), :)/REAL(i, dp)
199 : END DO
200 : END DO
201 0 : END SUBROUTINE staging_u2x
202 :
203 : ! ***************************************************************************
204 : !> \brief staging transformation for the forces
205 : !> \param staging_env the environment for the staging transformation
206 : !> \param uf will contain the forces after for the transformed variable
207 : !> \param f the forces to transform
208 : !> \author fawzi
209 : ! **************************************************************************************************
210 0 : PURE SUBROUTINE staging_f2uf(staging_env, uf, f)
211 : TYPE(staging_env_type), INTENT(IN) :: staging_env
212 : REAL(kind=dp), DIMENSION(:, :), INTENT(out) :: uf
213 : REAL(kind=dp), DIMENSION(:, :), INTENT(in) :: f
214 :
215 : INTEGER :: i, idim, ij, ist, k
216 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: iii, jjj, kkk
217 : REAL(kind=dp) :: const, sum_f
218 :
219 0 : const = REAL(staging_env%j - 1, dp)/REAL(staging_env%j, dp)
220 : ALLOCATE (iii(staging_env%j), jjj(staging_env%j), &
221 0 : kkk(staging_env%j))
222 0 : DO ist = 1, staging_env%j - 1
223 0 : iii(ist) = (ist - 1)*staging_env%j + 1 ! first el
224 0 : jjj(ist) = iii(ist) + staging_env%j - 1 ! last el
225 0 : kkk(ist) = iii(ist) - 1 ! prev el
226 : END DO
227 0 : kkk(1) = staging_env%p
228 :
229 0 : uf = f
230 : ! staging beads
231 0 : DO k = 1, staging_env%nseg
232 0 : DO i = 2, staging_env%j
233 : uf(i + iii(k), :) = uf(i + iii(k), :) &
234 0 : + REAL(i - 1, dp)/REAL(i, dp)*uf(i + iii(k) - 1, :)
235 : END DO
236 : END DO
237 : ! end point beads
238 0 : DO idim = 1, SIZE(uf, 2)
239 0 : DO k = 1, staging_env%nseg
240 : sum_f = 0._dp
241 0 : DO ij = 2, staging_env%j - 1
242 0 : sum_f = sum_f + uf((k - 1)*staging_env%j + ij, idim)
243 : END DO
244 : uf(iii(k), idim) = uf(iii(k), idim) + &
245 0 : sum_f - const*(uf(jjj(k), idim) - uf(kkk(k), idim))
246 : END DO
247 : END DO
248 0 : END SUBROUTINE staging_f2uf
249 :
250 : ! ***************************************************************************
251 : !> \brief calculates the harmonic force in the staging basis
252 : !> \param staging_env the staging environment
253 : !> \param mass_beads the masses of the beads
254 : !> \param ux the positions of the beads in the staging basis
255 : !> \param uf_h the harmonic forces (not accelerations)
256 : !> \param e_h ...
257 : !> \author fawzi
258 : ! **************************************************************************************************
259 0 : PURE SUBROUTINE staging_calc_uf_h(staging_env, mass_beads, ux, uf_h, e_h)
260 : TYPE(staging_env_type), INTENT(IN) :: staging_env
261 : REAL(kind=dp), DIMENSION(:, :), POINTER :: mass_beads, ux, uf_h
262 : REAL(KIND=dp), INTENT(OUT) :: e_h
263 :
264 : INTEGER :: idim, isg, ist
265 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: iii, jjj, kkk
266 : REAL(kind=dp) :: d, f
267 :
268 0 : e_h = 0.0_dp
269 :
270 : ALLOCATE (iii(staging_env%nseg), jjj(staging_env%nseg), &
271 0 : kkk(staging_env%nseg))
272 :
273 0 : DO ist = 1, staging_env%nseg
274 0 : iii(ist) = (ist - 1)*staging_env%j + 1 ! first el
275 0 : jjj(ist) = iii(ist) + staging_env%j ! next first (pbc)
276 0 : kkk(ist) = iii(ist) - staging_env%j ! prev first el (pbc)
277 : END DO
278 0 : jjj(staging_env%nseg) = 1
279 0 : kkk(1) = staging_env%p - staging_env%j
280 :
281 0 : DO idim = 1, SIZE(mass_beads, 2)
282 0 : DO ist = 1, staging_env%nseg
283 : e_h = e_h + 0.5*mass_beads(1, idim)*staging_env%w_j**2* &
284 0 : (ux(iii(ist), idim) - ux(jjj(ist), idim))**2
285 : uf_h(iii(ist), idim) = mass_beads(1, idim)*staging_env%w_j**2*( &
286 : 2._dp*ux(iii(ist), idim) &
287 : - ux(jjj(ist), idim) &
288 : - ux(kkk(ist), idim) &
289 0 : )
290 0 : DO isg = 2, staging_env%j ! use 3 as start?
291 0 : d = ux((ist - 1)*staging_env%j + isg, idim)
292 0 : f = mass_beads((ist - 1)*staging_env%j + isg, idim)*staging_env%w_j**2*d
293 0 : e_h = e_h + 0.5_dp*f*d
294 0 : uf_h((ist - 1)*staging_env%j + isg, idim) = f
295 : END DO
296 : END DO
297 : END DO
298 0 : END SUBROUTINE staging_calc_uf_h
299 :
300 : END MODULE pint_staging
|
