1!--------------------------------------------------------------------------------------------------!
2!   CP2K: A general program to perform molecular dynamics simulations                              !
3!   Copyright (C) 2000 - 2019  CP2K developers group                                               !
4!--------------------------------------------------------------------------------------------------!
5
6! **************************************************************************************************
7!> \brief Type for the canonical sampling through velocity rescaling
8!> \author Teodoro Laino - 09.2007 University of Zurich [tlaino]
9! **************************************************************************************************
10MODULE csvr_system_types
11   USE bibliography,                    ONLY: Bussi2007,&
12                                              cite_reference
13   USE extended_system_types,           ONLY: create_map_info_type,&
14                                              map_info_type,&
15                                              release_map_info_type
16   USE input_section_types,             ONLY: section_vals_type,&
17                                              section_vals_val_get
18   USE kinds,                           ONLY: dp
19   USE parallel_rng_types,              ONLY: GAUSSIAN,&
20                                              create_rng_stream,&
21                                              delete_rng_stream,&
22                                              next_rng_seed,&
23                                              rng_stream_type
24   USE simpar_types,                    ONLY: simpar_type
25   USE string_utilities,                ONLY: compress
26#include "./base/base_uses.f90"
27
28   IMPLICIT NONE
29
30   PRIVATE
31   PUBLIC :: csvr_system_type, &
32             csvr_init, &
33             csvr_dealloc, &
34             csvr_thermo_create
35
36! **************************************************************************************************
37   TYPE csvr_thermo_type
38      INTEGER                                 :: degrees_of_freedom
39      REAL(KIND=dp)                           :: nkt
40      REAL(KIND=dp)                           :: thermostat_energy
41      REAL(KIND=dp)                           :: region_kin_energy
42      TYPE(rng_stream_type), POINTER          :: gaussian_rng_stream
43   END TYPE csvr_thermo_type
44
45! **************************************************************************************************
46   TYPE csvr_system_type
47      INTEGER                                 :: region, glob_num_csvr, loc_num_csvr
48      REAL(KIND=dp)                           :: tau_csvr, dt_fact
49      TYPE(csvr_thermo_type), POINTER         :: nvt(:)
50      TYPE(map_info_type), POINTER            :: map_info
51   END TYPE csvr_system_type
52
53! *** Global parameters ***
54   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'csvr_system_types'
55
56CONTAINS
57
58! **************************************************************************************************
59!> \brief Initialize type for Canonical Sampling through Velocity Rescaling (CSVR)
60!> \param csvr ...
61!> \param simpar ...
62!> \param section ...
63!> \author Teodoro Laino [tlaino] 10.2007- University of Zurich
64! **************************************************************************************************
65   SUBROUTINE csvr_init(csvr, simpar, section)
66      TYPE(csvr_system_type), POINTER                    :: csvr
67      TYPE(simpar_type), POINTER                         :: simpar
68      TYPE(section_vals_type), POINTER                   :: section
69
70      CHARACTER(LEN=*), PARAMETER :: routineN = 'csvr_init', routineP = moduleN//':'//routineN
71
72      NULLIFY (csvr%nvt)
73      NULLIFY (csvr%map_info)
74      csvr%loc_num_csvr = 0
75      csvr%glob_num_csvr = 0
76      csvr%dt_fact = 1.0_dp
77      CALL cite_reference(Bussi2007)
78      CALL section_vals_val_get(section, "TIMECON", r_val=csvr%tau_csvr)
79      ! The CSVR library expects the tau_csv to be in unit of integration timestep
80      ! if applied once.. divided by two if the process is applied both to the first
81      ! and the second verlet step
82      csvr%tau_csvr = csvr%tau_csvr/(0.5_dp*simpar%dt)
83      CALL create_map_info_type(csvr%map_info)
84
85   END SUBROUTINE csvr_init
86
87! **************************************************************************************************
88!> \brief Initialize NVT type for CSVR thermostat
89!> \param csvr ...
90!> \author Teodoro Laino [tlaino] 10.2007- University of Zurich
91! **************************************************************************************************
92   SUBROUTINE csvr_thermo_create(csvr)
93      TYPE(csvr_system_type), POINTER                    :: csvr
94
95      CHARACTER(LEN=*), PARAMETER :: routineN = 'csvr_thermo_create', &
96         routineP = moduleN//':'//routineN
97
98      CHARACTER(LEN=40)                                  :: name
99      INTEGER                                            :: i, ithermo, my_index
100      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: seed
101      REAL(KIND=dp), DIMENSION(3, 2)                     :: initial_seed, my_seed
102
103      CPASSERT(ASSOCIATED(csvr))
104      CPASSERT(.NOT. ASSOCIATED(csvr%nvt))
105
106      ALLOCATE (csvr%nvt(csvr%loc_num_csvr))
107      DO i = 1, csvr%loc_num_csvr
108         csvr%nvt(i)%thermostat_energy = 0.0_dp
109         NULLIFY (csvr%nvt(i)%gaussian_rng_stream)
110      END DO
111      ! Initialize the gaussian stream random number
112      ALLOCATE (seed(3, 2, csvr%glob_num_csvr))
113      initial_seed = next_rng_seed()
114
115      seed(:, :, 1) = initial_seed
116      DO ithermo = 2, csvr%glob_num_csvr
117         seed(:, :, ithermo) = next_rng_seed(seed(:, :, ithermo - 1))
118      END DO
119      ! Update initial seed
120      initial_seed = next_rng_seed(seed(:, :, csvr%glob_num_csvr))
121      DO ithermo = 1, csvr%loc_num_csvr
122         my_index = csvr%map_info%index(ithermo)
123         my_seed = seed(:, :, my_index)
124         WRITE (UNIT=name, FMT="(A,I8)") "Wiener process for Thermostat #", my_index
125         CALL compress(name)
126         CALL create_rng_stream(rng_stream=csvr%nvt(ithermo)%gaussian_rng_stream, &
127                                name=name, distribution_type=GAUSSIAN, extended_precision=.TRUE., &
128                                seed=my_seed)
129      END DO
130      DEALLOCATE (seed)
131
132   END SUBROUTINE csvr_thermo_create
133
134! **************************************************************************************************
135!> \brief Deallocate type for CSVR thermostat
136!> \param csvr ...
137!> \author Teodoro Laino [tlaino] 10.2007- University of Zurich
138! **************************************************************************************************
139   SUBROUTINE csvr_dealloc(csvr)
140      TYPE(csvr_system_type), POINTER                    :: csvr
141
142      CHARACTER(LEN=*), PARAMETER :: routineN = 'csvr_dealloc', routineP = moduleN//':'//routineN
143
144      IF (ASSOCIATED(csvr)) THEN
145         CALL csvr_thermo_dealloc(csvr%nvt)
146         CALL release_map_info_type(csvr%map_info)
147         DEALLOCATE (csvr)
148      ENDIF
149
150   END SUBROUTINE csvr_dealloc
151
152! **************************************************************************************************
153!> \brief Deallocate NVT type for CSVR thermostat
154!> \param nvt ...
155!> \author Teodoro Laino [tlaino] 10.2007- University of Zurich
156! **************************************************************************************************
157   SUBROUTINE csvr_thermo_dealloc(nvt)
158      TYPE(csvr_thermo_type), DIMENSION(:), POINTER      :: nvt
159
160      CHARACTER(LEN=*), PARAMETER :: routineN = 'csvr_thermo_dealloc', &
161         routineP = moduleN//':'//routineN
162
163      INTEGER                                            :: i
164
165      IF (ASSOCIATED(nvt)) THEN
166         DO i = 1, SIZE(nvt)
167            IF (ASSOCIATED(nvt(i)%gaussian_rng_stream)) THEN
168               CALL delete_rng_stream(nvt(i)%gaussian_rng_stream)
169            ENDIF
170         END DO
171         DEALLOCATE (nvt)
172      ENDIF
173   END SUBROUTINE csvr_thermo_dealloc
174
175END MODULE csvr_system_types
176
177