43#include "../../base/base_uses.f90"
54 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'gle_system_dynamics'
75 SUBROUTINE gle_particles(gle, molecule_kind_set, molecule_set, particle_set, local_molecules, &
76 group, shell_adiabatic, shell_particle_set, core_particle_set, vel, shell_vel, core_vel)
84 LOGICAL,
INTENT(IN),
OPTIONAL :: shell_adiabatic
85 TYPE(
particle_type),
OPTIONAL,
POINTER :: shell_particle_set(:), &
87 REAL(kind=
dp),
INTENT(INOUT),
OPTIONAL :: vel(:, :), shell_vel(:, :), &
90 CHARACTER(len=*),
PARAMETER :: routinen =
'gle_particles'
92 INTEGER :: handle, iadd, ideg, imap, ndim, num
93 LOGICAL :: my_shell_adiabatic, present_vel
94 REAL(
dp) :: alpha, beta, rr
95 REAL(
dp),
DIMENSION(:, :),
POINTER :: a_mat, e_tmp, h_tmp, s_tmp
98 CALL timeset(routinen, handle)
99 my_shell_adiabatic = .false.
100 IF (
PRESENT(shell_adiabatic)) my_shell_adiabatic = shell_adiabatic
101 present_vel =
PRESENT(vel)
103 ALLOCATE (s_tmp(ndim, gle%loc_num_gle))
105 ALLOCATE (e_tmp(ndim, gle%loc_num_gle))
106 ALLOCATE (h_tmp(ndim, gle%loc_num_gle))
108 map_info => gle%map_info
110 local_molecules, molecule_set, group, vel)
111 DO ideg = 1, gle%loc_num_gle
112 imap = gle%map_info%map_index(ideg)
113 gle%nvt(ideg)%kin_energy = map_info%s_kin(imap)
117 local_molecules, molecule_set, group, vel)
119 DO ideg = 1, gle%loc_num_gle
120 imap = gle%map_info%map_index(ideg)
121 IF (gle%nvt(ideg)%nkt == 0.0_dp) cycle
122 gle%nvt(ideg)%s(1) = map_info%s_kin(imap)
123 s_tmp(1, imap) = map_info%s_kin(imap)
124 rr = gle%nvt(ideg)%gaussian_rng_stream%next()
127 s_tmp(iadd, imap) = gle%nvt(ideg)%s(iadd)
128 rr = gle%nvt(ideg)%gaussian_rng_stream%next()
129 e_tmp(iadd, imap) = rr
132 num = gle%loc_num_gle
136 CALL dgemm(
'N',
'N', ndim, num, ndim, alpha, a_mat(1, 1), ndim, e_tmp(1, 1), ndim, beta, h_tmp(1, 1), ndim)
140 CALL dgemm(
"N",
"N", ndim, num, ndim, alpha, a_mat(1, 1), ndim, s_tmp(1, 1), ndim, beta, h_tmp(1, 1), ndim)
142 DO ideg = 1, gle%loc_num_gle
143 imap = gle%map_info%map_index(ideg)
144 IF (gle%nvt(ideg)%nkt == 0.0_dp) cycle
146 map_info%v_scale(imap) = h_tmp(1, imap)/s_tmp(1, imap)
148 gle%nvt(ideg)%s(iadd) = h_tmp(iadd, ideg)
153 local_molecules, my_shell_adiabatic, shell_particle_set, core_particle_set, &
154 vel, shell_vel, core_vel)
157 local_molecules, molecule_set, group, vel)
158 DO ideg = 1, gle%loc_num_gle
159 imap = gle%map_info%map_index(ideg)
160 gle%nvt(ideg)%thermostat_energy = gle%nvt(ideg)%thermostat_energy + &
161 0.5_dp*(gle%nvt(ideg)%kin_energy - map_info%s_kin(imap))
164 DEALLOCATE (e_tmp, s_tmp, h_tmp)
166 CALL timestop(handle)
75 SUBROUTINE gle_particles(gle, molecule_kind_set, molecule_set, particle_set, local_molecules, &
…
184 molecule, molecule_kind_set, para_env, gle, gle_section, &
196 LOGICAL,
INTENT(IN) :: save_mem
199 REAL(
dp) :: mtmp(gle%ndim, gle%ndim)
203 CALL gle_to_particle_mapping(thermostat_info, simpar, local_molecules, &
204 molecule, molecule_kind_set, gle, para_env, gci)
206 IF (gle%ndim /= 0)
THEN
207 CALL init_gle_variables(gle)
209 CALL restart_gle(gle, gle_section, save_mem, restart)
213 CALL gle_matrix_exp((-simpar%dt*0.5_dp)*gle%a_mat, gle%ndim, 15, 15, gle%gle_t)
215 mtmp = gle%c_mat - matmul(gle%gle_t, matmul(gle%c_mat, transpose(gle%gle_t)))
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)
238 REAL(
dp) :: sm(n, n), tc(j + 1)
242 tc(i + 1) = tc(i)/real(i, kind=
dp)
246 sm = m*(1._dp/2._dp**k)
256 em(i, i) = em(i, i) + tc(p)
278 INTEGER,
INTENT(in) :: n
279 REAL(
dp),
INTENT(out) :: s(n, n)
280 REAL(
dp),
INTENT(in) :: sst(n, n)
283 REAL(
dp) :: d(n), l(n, n)
294 l(i, j) = l(i, j) - l(i, k)*l(j, k)*d(k)
296 IF (abs(d(j)) > epsilon(1.0_dp)) l(i, j) = l(i, j)/d(j)
299 d(i) = d(i) - l(i, k)*l(i, k)*d(k)
304 IF ((abs(d(j)) > epsilon(1.0_dp)) .AND. (d(j) > 0.0_dp))
THEN
305 s(i, j) = s(i, j) + l(i, j)*sqrt(d(j))
324 SUBROUTINE gle_to_particle_mapping(thermostat_info, simpar, local_molecules, &
325 molecule_set, molecule_kind_set, gle, para_env, gci)
336 INTEGER :: i, imap, j, mal_size, natoms_local, &
337 nkind, number, region, &
339 INTEGER,
DIMENSION(:),
POINTER :: deg_of_freedom, massive_atom_list
345 NULLIFY (massive_atom_list, deg_of_freedom)
346 SELECT CASE (simpar%ensemble)
348 cpabort(
"Unknown ensemble!")
351 cpabort(
"Never reach this point!")
354 map_info => gle%map_info
355 nkind =
SIZE(molecule_kind_set)
356 sum_of_thermostats = thermostat_info%sum_of_thermostats
357 map_info%dis_type = thermostat_info%dis_type
358 number = thermostat_info%number_of_thermostats
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, &
367 gle%loc_num_gle = number
368 gle%glob_num_gle = sum_of_thermostats
369 mal_size =
SIZE(massive_atom_list)
370 cpassert(mal_size /= 0)
372 gle%mal(1:mal_size) = massive_atom_list(1:mal_size)
376 map_info%s_kin = 0.0_dp
378 DO j = 1, natoms_local
379 map_info%p_kin(i, j)%point = map_info%p_kin(i, j)%point + 1
389 fac = map_info%s_kin(1) - deg_of_freedom(1) - simpar%nfree_rot_transl
390 IF (
fac == 0.0_dp)
THEN
391 cpabort(
"Zero degrees of freedom. Nothing to thermalize!")
393 gle%nvt(1)%nkt = simpar%temp_ext*
fac
394 gle%nvt(1)%degrees_of_freedom = floor(
fac)
396 DO i = 1, gle%loc_num_gle
397 imap = map_info%map_index(i)
398 fac = (map_info%s_kin(imap) - deg_of_freedom(i))
399 gle%nvt(i)%nkt = simpar%temp_ext*
fac
400 gle%nvt(i)%degrees_of_freedom = floor(
fac)
403 DEALLOCATE (deg_of_freedom)
404 DEALLOCATE (massive_atom_list)
407 END SUBROUTINE gle_to_particle_mapping
420 LOGICAL,
INTENT(IN) :: save_mem
421 LOGICAL,
INTENT(OUT) :: restart
423 CHARACTER(LEN=rng_record_length) :: rng_record
424 INTEGER :: glob_num, i, ind, j, loc_num, n_rep
426 REAL(kind=
dp),
DIMENSION(:),
POINTER :: buffer
430 NULLIFY (buffer, work_section)
435 IF (
ASSOCIATED(gle_section))
THEN
442 map_info => gle%map_info
444 DO i = 1,
SIZE(gle%nvt, 1)
445 ind = map_info%index(i)
446 ind = (ind - 1)*(gle%ndim)
447 DO j = 1,
SIZE(gle%nvt(i)%s, 1)
449 gle%nvt(i)%s(j) = buffer(ind)
454 NULLIFY (work_section)
461 subsection_name=
"THERMOSTAT_ENERGY")
466 IF (n_rep == gle%glob_num_gle)
THEN
467 DO i = 1, gle%loc_num_gle
468 ind = map_info%index(i)
470 i_rep_val=ind, r_val=gle%nvt(i)%thermostat_energy)
473 CALL cp_abort(__location__, &
474 "Number of thermostat energies not equal to the number of "// &
475 "total thermostats!")
481 subsection_name=
"RNG_INIT")
488 glob_num = gle%glob_num_gle
489 loc_num = gle%loc_num_gle
490 IF (n_rep == glob_num)
THEN
492 ind = map_info%index(i)
494 i_rep_val=ind, c_val=rng_record)
498 CALL cp_abort(__location__, &
499 "Number pf restartable stream not equal to the number of "// &
500 "total thermostats!")
511 SUBROUTINE init_gle_variables(gle)
516 REAL(
dp) :: rr(gle%ndim), cc(gle%ndim, gle%ndim)
519 DO i = 1, gle%loc_num_gle
522 rr(j) = gle%nvt(i)%gaussian_rng_stream%next()
524 gle%nvt(i)%s = matmul(cc, rr)
527 END SUBROUTINE init_gle_variables
static GRID_HOST_DEVICE double fac(const int i)
Factorial function, e.g. fac(5) = 5! = 120.
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
Lumps all possible extended system variables into one type for easy access and passing.
subroutine, public gle_particles(gle, molecule_kind_set, molecule_set, particle_set, local_molecules, group, shell_adiabatic, shell_particle_set, core_particle_set, vel, shell_vel, core_vel)
...
subroutine, public restart_gle(gle, gle_section, save_mem, restart)
...
subroutine, public gle_matrix_exp(m, n, j, k, em)
...
subroutine, public gle_cholesky_stab(sst, s, n)
...
subroutine, public initialize_gle_part(thermostat_info, simpar, local_molecules, molecule, molecule_kind_set, para_env, gle, gle_section, gci, save_mem)
...
subroutine, public gle_thermo_create(gle, mal_size)
...
Defines the basic variable types.
integer, parameter, public dp
Interface to the message passing library MPI.
Define the molecule kind structure types and the corresponding functionality.
Define the data structure for the molecule information.
Parallel (pseudo)random number generator (RNG) for multiple streams and substreams of random numbers.
type(rng_stream_type) function, public rng_stream_type_from_record(rng_record)
Create a RNG stream from a record given as an internal file (string).
integer, parameter, public rng_record_length
Define the data structure for the particle information.
Type for storing MD parameters.
subroutine, public thermostat_mapping_region(map_info, deg_of_freedom, massive_atom_list, molecule_kind_set, local_molecules, molecule_set, para_env, natoms_local, simpar, number, region, gci, shell, map_loc_thermo_gen, sum_of_thermostats)
Main general setup thermostat regions (thermostat independent)
Thermostat structure: module containing thermostat available for MD.
Utilities for thermostats.
subroutine, public momentum_region_particles(map_info, particle_set, molecule_kind_set, local_molecules, molecule_set, group, vel)
...
subroutine, public vel_rescale_particles(map_info, molecule_kind_set, molecule_set, particle_set, local_molecules, shell_adiabatic, shell_particle_set, core_particle_set, vel, shell_vel, core_vel)
...
subroutine, public ke_region_particles(map_info, particle_set, molecule_kind_set, local_molecules, molecule_set, group, vel)
...
structure to store local (to a processor) ordered lists of integers.
stores all the informations relevant to an mpi environment
Simulation parameter type for molecular dynamics.