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