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 : !> \par History
10 : !> \author
11 : ! **************************************************************************************************
12 : MODULE gle_system_dynamics
13 :
14 : USE distribution_1d_types, ONLY: distribution_1d_type
15 : USE extended_system_types, ONLY: map_info_type
16 : USE gle_system_types, ONLY: gle_thermo_create,&
17 : gle_type
18 : USE input_constants, ONLY: &
19 : do_thermo_communication, do_thermo_no_communication, isokin_ensemble, langevin_ensemble, &
20 : npe_f_ensemble, npe_i_ensemble, nph_uniaxial_damped_ensemble, nph_uniaxial_ensemble, &
21 : npt_f_ensemble, npt_i_ensemble, npt_ia_ensemble, nve_ensemble, nvt_ensemble, &
22 : reftraj_ensemble
23 : USE input_section_types, ONLY: section_vals_get,&
24 : section_vals_get_subs_vals,&
25 : section_vals_remove_values,&
26 : section_vals_type,&
27 : section_vals_val_get
28 : USE kinds, ONLY: dp
29 : USE message_passing, ONLY: mp_comm_type,&
30 : mp_para_env_type
31 : USE molecule_kind_types, ONLY: molecule_kind_type
32 : USE molecule_types, ONLY: global_constraint_type,&
33 : molecule_type
34 : USE parallel_rng_types, ONLY: rng_record_length,&
35 : rng_stream_type_from_record
36 : USE particle_types, ONLY: particle_type
37 : USE simpar_types, ONLY: simpar_type
38 : USE thermostat_mapping, ONLY: thermostat_mapping_region
39 : USE thermostat_types, ONLY: thermostat_info_type
40 : USE thermostat_utils, ONLY: ke_region_particles,&
41 : momentum_region_particles,&
42 : vel_rescale_particles
43 : #include "../../base/base_uses.f90"
44 :
45 : IMPLICIT NONE
46 : PRIVATE
47 :
48 : PUBLIC :: gle_particles, &
49 : initialize_gle_part, &
50 : gle_cholesky_stab, &
51 : gle_matrix_exp, &
52 : restart_gle
53 :
54 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gle_system_dynamics'
55 :
56 : CONTAINS
57 :
58 : ! **************************************************************************************************
59 : !> \brief ...
60 : !> \param gle ...
61 : !> \param molecule_kind_set ...
62 : !> \param molecule_set ...
63 : !> \param particle_set ...
64 : !> \param local_molecules ...
65 : !> \param group ...
66 : !> \param shell_adiabatic ...
67 : !> \param shell_particle_set ...
68 : !> \param core_particle_set ...
69 : !> \param vel ...
70 : !> \param shell_vel ...
71 : !> \param core_vel ...
72 : !> \date
73 : !> \par History
74 : ! **************************************************************************************************
75 800 : SUBROUTINE gle_particles(gle, molecule_kind_set, molecule_set, particle_set, local_molecules, &
76 800 : group, shell_adiabatic, shell_particle_set, core_particle_set, vel, shell_vel, core_vel)
77 :
78 : TYPE(gle_type), POINTER :: gle
79 : TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
80 : TYPE(molecule_type), POINTER :: molecule_set(:)
81 : TYPE(particle_type), POINTER :: particle_set(:)
82 : TYPE(distribution_1d_type), POINTER :: local_molecules
83 : TYPE(mp_comm_type), INTENT(IN) :: group
84 : LOGICAL, INTENT(IN), OPTIONAL :: shell_adiabatic
85 : TYPE(particle_type), OPTIONAL, POINTER :: shell_particle_set(:), &
86 : core_particle_set(:)
87 : REAL(KIND=dp), INTENT(INOUT), OPTIONAL :: vel(:, :), shell_vel(:, :), &
88 : core_vel(:, :)
89 :
90 : CHARACTER(len=*), PARAMETER :: routineN = 'gle_particles'
91 :
92 : INTEGER :: handle, iadd, ideg, imap, ndim, num
93 : LOGICAL :: my_shell_adiabatic, present_vel
94 : REAL(dp) :: alpha, beta, rr
95 800 : REAL(dp), DIMENSION(:, :), POINTER :: a_mat, e_tmp, h_tmp, s_tmp
96 : TYPE(map_info_type), POINTER :: map_info
97 :
98 800 : CALL timeset(routineN, handle)
99 800 : my_shell_adiabatic = .FALSE.
100 800 : IF (PRESENT(shell_adiabatic)) my_shell_adiabatic = shell_adiabatic
101 800 : present_vel = PRESENT(vel)
102 800 : ndim = gle%ndim
103 3200 : ALLOCATE (s_tmp(ndim, gle%loc_num_gle))
104 778400 : s_tmp = 0.0_dp
105 3200 : ALLOCATE (e_tmp(ndim, gle%loc_num_gle))
106 3200 : ALLOCATE (h_tmp(ndim, gle%loc_num_gle))
107 :
108 800 : map_info => gle%map_info
109 : CALL ke_region_particles(map_info, particle_set, molecule_kind_set, &
110 1200 : local_molecules, molecule_set, group, vel)
111 130400 : DO ideg = 1, gle%loc_num_gle
112 129600 : imap = gle%map_info%map_index(ideg)
113 130400 : gle%nvt(ideg)%kin_energy = map_info%s_kin(imap)
114 : END DO
115 :
116 : CALL momentum_region_particles(map_info, particle_set, molecule_kind_set, &
117 1200 : local_molecules, molecule_set, group, vel)
118 :
119 130400 : DO ideg = 1, gle%loc_num_gle
120 129600 : imap = gle%map_info%map_index(ideg)
121 129600 : IF (gle%nvt(ideg)%nkt == 0.0_dp) CYCLE
122 129600 : gle%nvt(ideg)%s(1) = map_info%s_kin(imap)
123 129600 : s_tmp(1, imap) = map_info%s_kin(imap)
124 129600 : rr = gle%nvt(ideg)%gaussian_rng_stream%next()
125 129600 : e_tmp(1, imap) = rr
126 648800 : DO iadd = 2, ndim
127 518400 : s_tmp(iadd, imap) = gle%nvt(ideg)%s(iadd)
128 518400 : rr = gle%nvt(ideg)%gaussian_rng_stream%next()
129 648000 : e_tmp(iadd, imap) = rr
130 : END DO
131 : END DO
132 800 : num = gle%loc_num_gle
133 800 : a_mat => gle%gle_s
134 800 : alpha = 1.0_dp
135 800 : beta = 0.0_dp
136 800 : CALL DGEMM('N', 'N', ndim, num, ndim, alpha, a_mat(1, 1), ndim, e_tmp(1, 1), ndim, beta, h_tmp(1, 1), ndim)
137 : !
138 800 : a_mat => gle%gle_t
139 800 : beta = 1.0_dp
140 800 : CALL dgemm("N", "N", ndim, num, ndim, alpha, a_mat(1, 1), ndim, s_tmp(1, 1), ndim, beta, h_tmp(1, 1), ndim)
141 :
142 130400 : DO ideg = 1, gle%loc_num_gle
143 129600 : imap = gle%map_info%map_index(ideg)
144 129600 : IF (gle%nvt(ideg)%nkt == 0.0_dp) CYCLE
145 :
146 129600 : map_info%v_scale(imap) = h_tmp(1, imap)/s_tmp(1, imap)
147 648800 : DO iadd = 2, ndim
148 648000 : gle%nvt(ideg)%s(iadd) = h_tmp(iadd, ideg)
149 : END DO
150 : END DO
151 :
152 : CALL vel_rescale_particles(map_info, molecule_kind_set, molecule_set, particle_set, &
153 : local_molecules, my_shell_adiabatic, shell_particle_set, core_particle_set, &
154 4400 : vel, shell_vel, core_vel)
155 :
156 : CALL ke_region_particles(map_info, particle_set, molecule_kind_set, &
157 1200 : local_molecules, molecule_set, group, vel)
158 130400 : DO ideg = 1, gle%loc_num_gle
159 129600 : imap = gle%map_info%map_index(ideg)
160 : gle%nvt(ideg)%thermostat_energy = gle%nvt(ideg)%thermostat_energy + &
161 130400 : 0.5_dp*(gle%nvt(ideg)%kin_energy - map_info%s_kin(imap))
162 : END DO
163 :
164 800 : DEALLOCATE (e_tmp, s_tmp, h_tmp)
165 :
166 800 : CALL timestop(handle)
167 800 : END SUBROUTINE gle_particles
168 :
169 : ! **************************************************************************************************
170 : !> \brief ...
171 : !> \param thermostat_info ...
172 : !> \param simpar ...
173 : !> \param local_molecules ...
174 : !> \param molecule ...
175 : !> \param molecule_kind_set ...
176 : !> \param para_env ...
177 : !> \param gle ...
178 : !> \param gle_section ...
179 : !> \param gci ...
180 : !> \param save_mem ...
181 : !> \author
182 : ! **************************************************************************************************
183 4 : SUBROUTINE initialize_gle_part(thermostat_info, simpar, local_molecules, &
184 : molecule, molecule_kind_set, para_env, gle, gle_section, &
185 : gci, save_mem)
186 :
187 : TYPE(thermostat_info_type), POINTER :: thermostat_info
188 : TYPE(simpar_type), POINTER :: simpar
189 : TYPE(distribution_1d_type), POINTER :: local_molecules
190 : TYPE(molecule_type), POINTER :: molecule(:)
191 : TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
192 : TYPE(mp_para_env_type), POINTER :: para_env
193 : TYPE(gle_type), POINTER :: gle
194 : TYPE(section_vals_type), POINTER :: gle_section
195 : TYPE(global_constraint_type), POINTER :: gci
196 : LOGICAL, INTENT(IN) :: save_mem
197 :
198 : LOGICAL :: restart
199 8 : REAL(dp) :: Mtmp(gle%ndim, gle%ndim)
200 :
201 : restart = .FALSE.
202 :
203 : CALL gle_to_particle_mapping(thermostat_info, simpar, local_molecules, &
204 4 : molecule, molecule_kind_set, gle, para_env, gci)
205 :
206 4 : IF (gle%ndim /= 0) THEN
207 4 : CALL init_gle_variables(gle)
208 : END IF
209 4 : CALL restart_gle(gle, gle_section, save_mem, restart)
210 :
211 : ! here we should have read a_mat and c_mat; whe can therefore compute S and T
212 : ! deterministic part of the propagator
213 124 : CALL gle_matrix_exp((-simpar%dt*0.5_dp)*gle%a_mat, gle%ndim, 15, 15, gle%gle_t)
214 : ! stochastic part
215 4624 : Mtmp = gle%c_mat - MATMUL(gle%gle_t, MATMUL(gle%c_mat, TRANSPOSE(gle%gle_t)))
216 4 : CALL gle_cholesky_stab(Mtmp, gle%gle_s, gle%ndim)
217 :
218 4 : END SUBROUTINE initialize_gle_part
219 :
220 : ! **************************************************************************************************
221 : !> \brief ...
222 : !> \param M ...
223 : !> \param n ...
224 : !> \param j ...
225 : !> \param k ...
226 : !> \param EM ...
227 : !> \author Michele
228 : !> routine to compute matrix exponential via scale & square
229 : ! **************************************************************************************************
230 30 : SUBROUTINE gle_matrix_exp(M, n, j, k, EM)
231 :
232 : INTEGER, INTENT(in) :: n
233 : REAL(dp), INTENT(in) :: M(n, n)
234 : INTEGER, INTENT(in) :: j, k
235 : REAL(dp), INTENT(out) :: EM(n, n)
236 :
237 : INTEGER :: i, p
238 60 : REAL(dp) :: SM(n, n), tc(j + 1)
239 :
240 30 : tc(1) = 1._dp
241 480 : DO i = 1, j
242 480 : tc(i + 1) = tc(i)/REAL(i, KIND=dp)
243 : END DO
244 :
245 : !scale
246 2370 : SM = M*(1._dp/2._dp**k)
247 2370 : EM = 0._dp
248 276 : DO i = 1, n
249 276 : EM(i, i) = tc(j + 1)
250 : END DO
251 :
252 : !taylor exp of scaled matrix
253 480 : DO p = j, 1, -1
254 653130 : EM = MATMUL(SM, EM)
255 4170 : DO i = 1, n
256 4140 : EM(i, i) = EM(i, i) + tc(p)
257 : END DO
258 : END DO
259 :
260 : !square
261 480 : DO p = 1, k
262 1235640 : EM = MATMUL(EM, EM)
263 : END DO
264 30 : END SUBROUTINE gle_matrix_exp
265 :
266 : ! **************************************************************************************************
267 : !> \brief ...
268 : !> \param SST ...
269 : !> \param S ...
270 : !> \param n ...
271 : !> \author Michele
272 : !> "stabilized" cholesky to deal with small & negative diagonal elements,
273 : !> in practice a LDL^T decomposition is computed, and negative els. are zeroed
274 : !> \par History
275 : !> 05.11.2015: Bug fix: In rare cases D(j) is negative due to numerics [Felix Uhl]
276 : ! **************************************************************************************************
277 12 : SUBROUTINE gle_cholesky_stab(SST, S, n)
278 : INTEGER, INTENT(in) :: n
279 : REAL(dp), INTENT(out) :: S(n, n)
280 : REAL(dp), INTENT(in) :: SST(n, n)
281 :
282 : INTEGER :: i, j, k
283 12 : REAL(dp) :: D(n), L(n, n)
284 :
285 72 : D = 0._dp
286 372 : L = 0._dp
287 372 : S = 0._dp
288 72 : DO i = 1, n
289 60 : L(i, i) = 1.0_dp
290 60 : D(i) = SST(i, i)
291 180 : DO j = 1, i - 1
292 120 : L(i, j) = SST(i, j);
293 240 : DO k = 1, j - 1
294 240 : L(i, j) = L(i, j) - L(i, k)*L(j, k)*D(k)
295 : END DO
296 180 : IF (ABS(D(j)) > EPSILON(1.0_dp)) L(i, j) = L(i, j)/D(j)
297 : END DO
298 192 : DO k = 1, i - 1
299 180 : D(i) = D(i) - L(i, k)*L(i, k)*D(k)
300 : END DO
301 : END DO
302 72 : DO i = 1, n
303 252 : DO j = 1, i
304 240 : IF ((ABS(D(j)) > EPSILON(1.0_dp)) .AND. (D(j) > 0.0_dp)) THEN
305 180 : S(i, j) = S(i, j) + L(i, j)*SQRT(D(j))
306 : END IF
307 : END DO
308 : END DO
309 :
310 12 : END SUBROUTINE gle_cholesky_stab
311 :
312 : ! **************************************************************************************************
313 : !> \brief ...
314 : !> \param thermostat_info ...
315 : !> \param simpar ...
316 : !> \param local_molecules ...
317 : !> \param molecule_set ...
318 : !> \param molecule_kind_set ...
319 : !> \param gle ...
320 : !> \param para_env ...
321 : !> \param gci ...
322 : !> \author
323 : ! **************************************************************************************************
324 4 : SUBROUTINE gle_to_particle_mapping(thermostat_info, simpar, local_molecules, &
325 : molecule_set, molecule_kind_set, gle, para_env, gci)
326 :
327 : TYPE(thermostat_info_type), POINTER :: thermostat_info
328 : TYPE(simpar_type), POINTER :: simpar
329 : TYPE(distribution_1d_type), POINTER :: local_molecules
330 : TYPE(molecule_type), POINTER :: molecule_set(:)
331 : TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
332 : TYPE(gle_type), POINTER :: gle
333 : TYPE(mp_para_env_type), POINTER :: para_env
334 : TYPE(global_constraint_type), POINTER :: gci
335 :
336 : INTEGER :: i, imap, j, mal_size, natoms_local, &
337 : nkind, number, region, &
338 : sum_of_thermostats
339 4 : INTEGER, DIMENSION(:), POINTER :: deg_of_freedom, massive_atom_list
340 : LOGICAL :: do_shell
341 : REAL(KIND=dp) :: fac
342 : TYPE(map_info_type), POINTER :: map_info
343 :
344 4 : do_shell = .FALSE.
345 4 : NULLIFY (massive_atom_list, deg_of_freedom)
346 4 : SELECT CASE (simpar%ensemble)
347 : CASE DEFAULT
348 0 : CPABORT("Unknown ensemble!")
349 : CASE (nve_ensemble, isokin_ensemble, npe_f_ensemble, npe_i_ensemble, nph_uniaxial_ensemble, &
350 : nph_uniaxial_damped_ensemble, reftraj_ensemble, langevin_ensemble)
351 0 : CPABORT("Never reach this point!")
352 : CASE (nvt_ensemble, npt_i_ensemble, npt_f_ensemble, npt_ia_ensemble)
353 :
354 4 : map_info => gle%map_info
355 4 : nkind = SIZE(molecule_kind_set)
356 4 : sum_of_thermostats = thermostat_info%sum_of_thermostats
357 4 : map_info%dis_type = thermostat_info%dis_type
358 4 : number = thermostat_info%number_of_thermostats
359 4 : region = gle%region
360 :
361 : CALL thermostat_mapping_region(map_info, deg_of_freedom, massive_atom_list, &
362 : molecule_kind_set, local_molecules, molecule_set, para_env, natoms_local, &
363 : simpar, number, region, gci, do_shell, thermostat_info%map_loc_thermo_gen, &
364 4 : sum_of_thermostats)
365 :
366 : ! This is the local number of available thermostats
367 4 : gle%loc_num_gle = number
368 4 : gle%glob_num_gle = sum_of_thermostats
369 4 : mal_size = SIZE(massive_atom_list)
370 4 : CPASSERT(mal_size /= 0)
371 4 : CALL gle_thermo_create(gle, mal_size)
372 872 : gle%mal(1:mal_size) = massive_atom_list(1:mal_size)
373 :
374 : ! Sum up the number of degrees of freedom on each thermostat.
375 : ! first: initialize the target
376 652 : map_info%s_kin = 0.0_dp
377 16 : DO i = 1, 3
378 664 : DO j = 1, natoms_local
379 660 : map_info%p_kin(i, j)%point = map_info%p_kin(i, j)%point + 1
380 : END DO
381 : END DO
382 :
383 : ! If thermostats are replicated but molecules distributed, we have to
384 : ! sum s_kin over all processors
385 4 : IF (map_info%dis_type == do_thermo_communication) CALL para_env%sum(map_info%s_kin)
386 :
387 : ! We know the total number of system thermostats.
388 4 : IF ((sum_of_thermostats == 1) .AND. (map_info%dis_type /= do_thermo_no_communication)) THEN
389 0 : fac = map_info%s_kin(1) - deg_of_freedom(1) - simpar%nfree_rot_transl
390 0 : IF (fac == 0.0_dp) THEN
391 0 : CPABORT("Zero degrees of freedom. Nothing to thermalize!")
392 : END IF
393 0 : gle%nvt(1)%nkt = simpar%temp_ext*fac
394 0 : gle%nvt(1)%degrees_of_freedom = FLOOR(fac)
395 : ELSE
396 652 : DO i = 1, gle%loc_num_gle
397 648 : imap = map_info%map_index(i)
398 648 : fac = (map_info%s_kin(imap) - deg_of_freedom(i))
399 648 : gle%nvt(i)%nkt = simpar%temp_ext*fac
400 652 : gle%nvt(i)%degrees_of_freedom = FLOOR(fac)
401 : END DO
402 : END IF
403 4 : DEALLOCATE (deg_of_freedom)
404 8 : DEALLOCATE (massive_atom_list)
405 : END SELECT
406 :
407 4 : END SUBROUTINE gle_to_particle_mapping
408 :
409 : ! **************************************************************************************************
410 : !> \brief ...
411 : !> \param gle ...
412 : !> \param gle_section ...
413 : !> \param save_mem ...
414 : !> \param restart ...
415 : ! **************************************************************************************************
416 6 : SUBROUTINE restart_gle(gle, gle_section, save_mem, restart)
417 :
418 : TYPE(gle_type), POINTER :: gle
419 : TYPE(section_vals_type), POINTER :: gle_section
420 : LOGICAL, INTENT(IN) :: save_mem
421 : LOGICAL, INTENT(OUT) :: restart
422 :
423 : CHARACTER(LEN=rng_record_length) :: rng_record
424 : INTEGER :: glob_num, i, ind, j, loc_num, n_rep
425 : LOGICAL :: explicit
426 6 : REAL(KIND=dp), DIMENSION(:), POINTER :: buffer
427 : TYPE(map_info_type), POINTER :: map_info
428 : TYPE(section_vals_type), POINTER :: work_section
429 :
430 6 : NULLIFY (buffer, work_section)
431 :
432 6 : explicit = .FALSE.
433 6 : restart = .FALSE.
434 :
435 12 : IF (ASSOCIATED(gle_section)) THEN
436 6 : work_section => section_vals_get_subs_vals(gle_section, "S")
437 6 : CALL section_vals_get(work_section, explicit=explicit)
438 6 : restart = explicit
439 : END IF
440 :
441 6 : IF (restart) THEN
442 2 : map_info => gle%map_info
443 2 : CALL section_vals_val_get(gle_section, "S%_DEFAULT_KEYWORD_", r_vals=buffer)
444 326 : DO i = 1, SIZE(gle%nvt, 1)
445 324 : ind = map_info%index(i)
446 324 : ind = (ind - 1)*(gle%ndim)
447 1946 : DO j = 1, SIZE(gle%nvt(i)%s, 1)
448 1620 : ind = ind + 1
449 1944 : gle%nvt(i)%s(j) = buffer(ind)
450 : END DO
451 : END DO
452 :
453 2 : IF (save_mem) THEN
454 0 : NULLIFY (work_section)
455 0 : work_section => section_vals_get_subs_vals(gle_section, "S")
456 0 : CALL section_vals_remove_values(work_section)
457 : END IF
458 :
459 : ! Possibly restart the initial thermostat energy
460 : work_section => section_vals_get_subs_vals(section_vals=gle_section, &
461 2 : subsection_name="THERMOSTAT_ENERGY")
462 2 : CALL section_vals_get(work_section, explicit=explicit)
463 2 : IF (explicit) THEN
464 : CALL section_vals_val_get(section_vals=work_section, keyword_name="_DEFAULT_KEYWORD_", &
465 2 : n_rep_val=n_rep)
466 2 : IF (n_rep == gle%glob_num_gle) THEN
467 326 : DO i = 1, gle%loc_num_gle
468 324 : ind = map_info%index(i)
469 : CALL section_vals_val_get(section_vals=work_section, keyword_name="_DEFAULT_KEYWORD_", &
470 326 : i_rep_val=ind, r_val=gle%nvt(i)%thermostat_energy)
471 : END DO
472 : ELSE
473 : CALL cp_abort(__LOCATION__, &
474 : "Number of thermostat energies not equal to the number of "// &
475 0 : "total thermostats!")
476 : END IF
477 : END IF
478 :
479 : ! Possibly restart the random number generators for the different thermostats
480 : work_section => section_vals_get_subs_vals(section_vals=gle_section, &
481 2 : subsection_name="RNG_INIT")
482 :
483 2 : CALL section_vals_get(work_section, explicit=explicit)
484 2 : IF (explicit) THEN
485 : CALL section_vals_val_get(section_vals=work_section, keyword_name="_DEFAULT_KEYWORD_", &
486 2 : n_rep_val=n_rep)
487 :
488 2 : glob_num = gle%glob_num_gle
489 2 : loc_num = gle%loc_num_gle
490 2 : IF (n_rep == glob_num) THEN
491 326 : DO i = 1, loc_num
492 324 : ind = map_info%index(i)
493 : CALL section_vals_val_get(section_vals=work_section, keyword_name="_DEFAULT_KEYWORD_", &
494 324 : i_rep_val=ind, c_val=rng_record)
495 326 : gle%nvt(i)%gaussian_rng_stream = rng_stream_type_from_record(rng_record)
496 : END DO
497 : ELSE
498 : CALL cp_abort(__LOCATION__, &
499 : "Number pf restartable stream not equal to the number of "// &
500 0 : "total thermostats!")
501 : END IF
502 : END IF
503 : END IF
504 :
505 6 : END SUBROUTINE restart_gle
506 :
507 : ! **************************************************************************************************
508 : !> \brief ...
509 : !> \param gle ...
510 : ! **************************************************************************************************
511 4 : SUBROUTINE init_gle_variables(gle)
512 :
513 : TYPE(gle_type), POINTER :: gle
514 :
515 : INTEGER :: i, j
516 8 : REAL(dp) :: rr(gle%ndim), cc(gle%ndim, gle%ndim)
517 :
518 4 : CALL gle_cholesky_stab(gle%c_mat, cc, gle%ndim)
519 652 : DO i = 1, gle%loc_num_gle
520 3888 : DO j = 1, gle%ndim
521 : ! here s should be properly initialized, when it is not read from restart
522 3888 : rr(j) = gle%nvt(i)%gaussian_rng_stream%next()
523 : END DO
524 23980 : gle%nvt(i)%s = MATMUL(cc, rr)
525 : END DO
526 :
527 4 : END SUBROUTINE init_gle_variables
528 :
529 454 : END MODULE gle_system_dynamics
|
