Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 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 normal mode coords
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 : !> 2015-10 added alternative normal mode transformation needed by RPMD
16 : !> [Felix Uhl
17 : ! **************************************************************************************************
18 : MODULE pint_normalmode
19 : USE input_constants, ONLY: propagator_bcmd,&
20 : propagator_cmd,&
21 : propagator_pimd,&
22 : propagator_rpmd
23 : USE input_section_types, ONLY: section_vals_type,&
24 : section_vals_val_get
25 : USE kinds, ONLY: dp
26 : USE mathconstants, ONLY: pi,&
27 : twopi
28 : USE pint_types, ONLY: normalmode_env_type
29 : #include "../base/base_uses.f90"
30 :
31 : IMPLICIT NONE
32 : PRIVATE
33 :
34 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
35 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pint_normalmode'
36 :
37 : PUBLIC :: normalmode_env_create
38 : PUBLIC :: normalmode_release
39 : PUBLIC :: normalmode_init_masses
40 : PUBLIC :: normalmode_x2u
41 : PUBLIC :: normalmode_u2x
42 : PUBLIC :: normalmode_f2uf
43 : PUBLIC :: normalmode_calc_uf_h
44 :
45 : CONTAINS
46 :
47 : ! ***************************************************************************
48 : !> \brief creates the data needed for a normal mode transformation
49 : !> \param normalmode_env ...
50 : !> \param normalmode_section ...
51 : !> \param p ...
52 : !> \param kT ...
53 : !> \param propagator ...
54 : !> \author Harald Forbert
55 : ! **************************************************************************************************
56 58 : SUBROUTINE normalmode_env_create(normalmode_env, normalmode_section, p, kT, propagator)
57 : TYPE(normalmode_env_type), INTENT(OUT) :: normalmode_env
58 : TYPE(section_vals_type), POINTER :: normalmode_section
59 : INTEGER, INTENT(in) :: p
60 : REAL(kind=dp), INTENT(in) :: kT
61 : INTEGER, INTENT(in) :: propagator
62 :
63 : INTEGER :: i, j, k, li
64 : LOGICAL :: explicit_gamma, explicit_modefactor
65 : REAL(kind=dp) :: gamma_parameter, invsqrtp, pip, sqrt2p, &
66 : twopip
67 :
68 232 : ALLOCATE (normalmode_env%x2u(p, p))
69 174 : ALLOCATE (normalmode_env%u2x(p, p))
70 174 : ALLOCATE (normalmode_env%lambda(p))
71 :
72 58 : normalmode_env%p = p
73 :
74 : CALL section_vals_val_get(normalmode_section, "Q_CENTROID", &
75 58 : r_val=normalmode_env%Q_centroid)
76 : CALL section_vals_val_get(normalmode_section, "Q_BEAD", &
77 58 : r_val=normalmode_env%Q_bead)
78 : CALL section_vals_val_get(normalmode_section, "MODEFACTOR", &
79 : explicit=explicit_modefactor, &
80 58 : r_val=normalmode_env%modefactor)
81 : CALL section_vals_val_get(normalmode_section, "GAMMA", &
82 : r_val=gamma_parameter, &
83 58 : explicit=explicit_gamma)
84 :
85 58 : IF (explicit_modefactor .AND. explicit_gamma) THEN
86 0 : CPABORT("Both GAMMA and MODEFACTOR have been declared. Please use only one.")
87 : END IF
88 58 : IF (explicit_gamma) THEN
89 4 : normalmode_env%modefactor = 1.0_dp/gamma_parameter**2
90 : END IF
91 :
92 58 : IF (propagator == propagator_cmd) THEN
93 4 : IF (.NOT. explicit_gamma) THEN
94 0 : CPABORT("GAMMA needs to be specified with CMD PROPAGATOR")
95 : END IF
96 4 : IF (gamma_parameter <= 1.0_dp) THEN
97 0 : CPWARN("GAMMA should be larger than 1.0 for CMD PROPAGATOR")
98 : END IF
99 : END IF
100 :
101 58 : IF (normalmode_env%Q_centroid < 0.0_dp) THEN
102 58 : normalmode_env%Q_centroid = -normalmode_env%Q_centroid/(kT*p)
103 : END IF
104 58 : IF (normalmode_env%Q_bead < 0.0_dp) THEN
105 58 : normalmode_env%Q_bead = -normalmode_env%Q_bead/(kT*p)
106 : END IF
107 :
108 : !Use different normal mode transformations depending on the propagator
109 58 : IF (propagator == propagator_pimd .OR. propagator == propagator_cmd .OR. &
110 : propagator == propagator_bcmd) THEN
111 :
112 40 : IF (propagator == propagator_pimd .OR. propagator == propagator_bcmd) THEN
113 36 : normalmode_env%harm = p*kT*kT/normalmode_env%modefactor
114 4 : ELSE IF (propagator == propagator_cmd) THEN
115 4 : normalmode_env%harm = p*kT*kT*gamma_parameter*gamma_parameter
116 4 : normalmode_env%modefactor = 1.0_dp/(gamma_parameter*gamma_parameter)
117 : END IF
118 :
119 : ! set up the transformation matrices
120 224 : DO i = 1, p
121 184 : normalmode_env%lambda(i) = 2.0_dp*(1.0_dp - COS(pi*(i/2)*2.0_dp/p))
122 1248 : DO j = 1, p
123 1024 : k = ((i/2)*(j - 1))/p
124 1024 : k = (i/2)*(j - 1) - k*p
125 1024 : li = 2*(i - 2*(i/2))*p - p
126 1208 : normalmode_env%u2x(j, i) = SQRT(2.0_dp/p)*SIN(twopi*(k + 0.125_dp*li)/p)
127 : END DO
128 : END DO
129 40 : normalmode_env%lambda(1) = 1.0_dp/(p*normalmode_env%modefactor)
130 224 : DO i = 1, p
131 1248 : DO j = 1, p
132 : normalmode_env%x2u(i, j) = SQRT(normalmode_env%lambda(i)* &
133 : normalmode_env%modefactor)* &
134 1208 : normalmode_env%u2x(j, i)
135 : END DO
136 : END DO
137 224 : DO i = 1, p
138 1248 : DO j = 1, p
139 : normalmode_env%u2x(i, j) = normalmode_env%u2x(i, j)/ &
140 : SQRT(normalmode_env%lambda(j)* &
141 1208 : normalmode_env%modefactor)
142 : END DO
143 : END DO
144 224 : normalmode_env%lambda(:) = normalmode_env%harm
145 :
146 18 : ELSE IF (propagator == propagator_rpmd) THEN
147 18 : normalmode_env%harm = kT/normalmode_env%modefactor
148 18 : sqrt2p = SQRT(2.0_dp/REAL(p, dp))
149 18 : twopip = twopi/REAL(p, dp)
150 18 : pip = pi/REAL(p, dp)
151 18 : invsqrtp = 1.0_dp/SQRT(REAL(p, dp))
152 3226 : normalmode_env%x2u(:, :) = 0.0_dp
153 186 : normalmode_env%x2u(1, :) = invsqrtp
154 186 : DO j = 1, p
155 1688 : DO i = 2, p/2 + 1
156 1688 : normalmode_env%x2u(i, j) = sqrt2p*COS(twopip*(i - 1)*(j - 1))
157 : END DO
158 1538 : DO i = p/2 + 2, p
159 1520 : normalmode_env%x2u(i, j) = sqrt2p*SIN(twopip*(i - 1)*(j - 1))
160 : END DO
161 : END DO
162 18 : IF (MOD(p, 2) == 0) THEN
163 18 : DO i = 1, p - 1, 2
164 84 : normalmode_env%x2u(p/2 + 1, i) = invsqrtp
165 84 : normalmode_env%x2u(p/2 + 1, i + 1) = -1.0_dp*invsqrtp
166 : END DO
167 : END IF
168 :
169 6434 : normalmode_env%u2x = TRANSPOSE(normalmode_env%x2u)
170 :
171 : ! Setting up propagator frequencies for rpmd
172 18 : normalmode_env%lambda(1) = 0.0_dp
173 168 : DO i = 2, p
174 150 : normalmode_env%lambda(i) = 2.0_dp*normalmode_env%harm*SIN((i - 1)*pip)
175 168 : normalmode_env%lambda(i) = normalmode_env%lambda(i)*normalmode_env%lambda(i)
176 : END DO
177 18 : normalmode_env%harm = kT*kT
178 : ELSE
179 0 : CPABORT("UNKNOWN PROPAGATOR FOR PINT SELECTED")
180 : END IF
181 :
182 58 : END SUBROUTINE normalmode_env_create
183 :
184 : ! ***************************************************************************
185 : !> \brief releases the normalmode environment
186 : !> \param normalmode_env the normalmode_env to release
187 : !> \author Harald Forbert
188 : ! **************************************************************************************************
189 58 : PURE SUBROUTINE normalmode_release(normalmode_env)
190 :
191 : TYPE(normalmode_env_type), INTENT(INOUT) :: normalmode_env
192 :
193 58 : DEALLOCATE (normalmode_env%x2u)
194 58 : DEALLOCATE (normalmode_env%u2x)
195 58 : DEALLOCATE (normalmode_env%lambda)
196 :
197 58 : END SUBROUTINE normalmode_release
198 :
199 : ! ***************************************************************************
200 : !> \brief initializes the masses and fictitious masses compatible with the
201 : !> normal mode information
202 : !> \param normalmode_env the definition of the normal mode transformation
203 : !> \param mass *input* the masses of the particles
204 : !> \param mass_beads masses of the beads
205 : !> \param mass_fict the fictitious masses
206 : !> \param Q masses of the nose thermostats
207 : !> \author Harald Forbert
208 : ! **************************************************************************************************
209 58 : PURE SUBROUTINE normalmode_init_masses(normalmode_env, mass, mass_beads, mass_fict, &
210 58 : Q)
211 :
212 : TYPE(normalmode_env_type), INTENT(IN) :: normalmode_env
213 : REAL(kind=dp), DIMENSION(:), INTENT(in) :: mass
214 : REAL(kind=dp), DIMENSION(:, :), INTENT(out), &
215 : OPTIONAL :: mass_beads, mass_fict
216 : REAL(kind=dp), DIMENSION(:), INTENT(out), OPTIONAL :: Q
217 :
218 : INTEGER :: iat, ib
219 :
220 58 : IF (PRESENT(Q)) THEN
221 410 : Q = normalmode_env%Q_bead
222 58 : Q(1) = normalmode_env%Q_centroid
223 : END IF
224 58 : IF (PRESENT(mass_beads) .OR. PRESENT(mass_fict)) THEN
225 58 : IF (PRESENT(mass_beads)) THEN
226 64984 : DO iat = 1, SIZE(mass)
227 64926 : mass_beads(1, iat) = 0.0_dp
228 297706 : DO ib = 2, normalmode_env%p
229 297648 : mass_beads(ib, iat) = mass(iat)
230 : END DO
231 : END DO
232 : END IF
233 58 : IF (PRESENT(mass_fict)) THEN
234 64984 : DO iat = 1, SIZE(mass)
235 362632 : DO ib = 1, normalmode_env%p
236 362574 : mass_fict(ib, iat) = mass(iat)
237 : END DO
238 : END DO
239 : END IF
240 : END IF
241 :
242 58 : END SUBROUTINE normalmode_init_masses
243 :
244 : ! ***************************************************************************
245 : !> \brief Transforms from the x into the u variables using a normal mode
246 : !> transformation for the positions
247 : !> \param normalmode_env the environment for the normal mode transformation
248 : !> \param ux will contain the u variable
249 : !> \param x the positions to transform
250 : !> \author Harald Forbert
251 : ! **************************************************************************************************
252 97 : SUBROUTINE normalmode_x2u(normalmode_env, ux, x)
253 : TYPE(normalmode_env_type), INTENT(INOUT) :: normalmode_env
254 : REAL(kind=dp), DIMENSION(:, :), INTENT(out) :: ux
255 : REAL(kind=dp), DIMENSION(:, :), INTENT(in) :: x
256 :
257 : CALL DGEMM('N', 'N', normalmode_env%p, SIZE(x, 2), normalmode_env%p, 1.0_dp, &
258 : normalmode_env%x2u(1, 1), SIZE(normalmode_env%x2u, 1), x(1, 1), SIZE(x, 1), &
259 97 : 0.0_dp, ux, SIZE(ux, 1))
260 97 : END SUBROUTINE normalmode_x2u
261 :
262 : ! ***************************************************************************
263 : !> \brief transform from the u variable to the x (back normal mode
264 : !> transformation for the positions)
265 : !> \param normalmode_env the environment for the normal mode transformation
266 : !> \param ux the u variable (positions to be backtransformed)
267 : !> \param x will contain the positions
268 : !> \author Harald Forbert
269 : ! **************************************************************************************************
270 3247 : SUBROUTINE normalmode_u2x(normalmode_env, ux, x)
271 : TYPE(normalmode_env_type), INTENT(INOUT) :: normalmode_env
272 : REAL(kind=dp), DIMENSION(:, :), INTENT(in) :: ux
273 : REAL(kind=dp), DIMENSION(:, :), INTENT(out) :: x
274 :
275 : CALL DGEMM('N', 'N', normalmode_env%p, SIZE(ux, 2), normalmode_env%p, 1.0_dp, &
276 : normalmode_env%u2x(1, 1), SIZE(normalmode_env%u2x, 1), ux(1, 1), SIZE(ux, 1), &
277 3247 : 0.0_dp, x, SIZE(x, 1))
278 3247 : END SUBROUTINE normalmode_u2x
279 :
280 : ! ***************************************************************************
281 : !> \brief normalmode transformation for the forces
282 : !> \param normalmode_env the environment for the normal mode transformation
283 : !> \param uf will contain the forces for the transformed variables afterwards
284 : !> \param f the forces to transform
285 : !> \author Harald Forbert
286 : ! **************************************************************************************************
287 712 : SUBROUTINE normalmode_f2uf(normalmode_env, uf, f)
288 : TYPE(normalmode_env_type), INTENT(INOUT) :: normalmode_env
289 : REAL(kind=dp), DIMENSION(:, :), INTENT(out) :: uf
290 : REAL(kind=dp), DIMENSION(:, :), INTENT(in) :: f
291 :
292 : CALL DGEMM('T', 'N', normalmode_env%p, SIZE(f, 2), normalmode_env%p, 1.0_dp, &
293 : normalmode_env%u2x(1, 1), SIZE(normalmode_env%u2x, 1), f(1, 1), SIZE(f, 1), &
294 712 : 0.0_dp, uf, SIZE(uf, 1))
295 712 : END SUBROUTINE normalmode_f2uf
296 :
297 : ! ***************************************************************************
298 : !> \brief calculates the harmonic force in the normal mode basis
299 : !> \param normalmode_env the normal mode environment
300 : !> \param mass_beads the masses of the beads
301 : !> \param ux the positions of the beads in the staging basis
302 : !> \param uf_h the harmonic forces (not accelerations)
303 : !> \param e_h ...
304 : !> \author Harald Forbert
305 : ! **************************************************************************************************
306 1446 : PURE SUBROUTINE normalmode_calc_uf_h(normalmode_env, mass_beads, ux, uf_h, e_h)
307 : TYPE(normalmode_env_type), INTENT(IN) :: normalmode_env
308 : REAL(kind=dp), DIMENSION(:, :), POINTER :: mass_beads, ux, uf_h
309 : REAL(KIND=dp), INTENT(OUT) :: e_h
310 :
311 : INTEGER :: ibead, idim
312 : REAL(kind=dp) :: f
313 :
314 1446 : e_h = 0.0_dp
315 448782 : DO idim = 1, SIZE(mass_beads, 2)
316 :
317 : ! starting at 2 since the centroid is at 1 and it's mass_beads
318 : ! SHOULD be zero anyways:
319 :
320 447336 : uf_h(1, idim) = 0.0_dp
321 2124294 : DO ibead = 2, normalmode_env%p
322 1675512 : f = -mass_beads(ibead, idim)*normalmode_env%lambda(ibead)*ux(ibead, idim)
323 1675512 : uf_h(ibead, idim) = f
324 : ! - to cancel the - in the force f.
325 2122848 : e_h = e_h - 0.5_dp*ux(ibead, idim)*f
326 : END DO
327 :
328 : END DO
329 1446 : END SUBROUTINE normalmode_calc_uf_h
330 :
331 : END MODULE pint_normalmode
|
