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