LCOV - code coverage report
Current view: top level - src/motion/thermostat - gle_system_dynamics.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 93.7 % 189 177
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 7 7

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 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         2400 :       ALLOCATE (e_tmp(ndim, gle%loc_num_gle))
     106         2400 :       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         2800 :                                  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       379440 :          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       414570 :          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        23332 :          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
        

Generated by: LCOV version 2.0-1