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 : !> \brief
10 : !> \par History
11 : !> \author MI 07.2009
12 : ! **************************************************************************************************
13 : MODULE gle_system_types
14 : USE bibliography, ONLY: Ceriotti2009,&
15 : Ceriotti2009b,&
16 : cite_reference
17 : USE extended_system_types, ONLY: create_map_info_type,&
18 : map_info_type,&
19 : release_map_info_type
20 : USE input_section_types, ONLY: section_vals_type,&
21 : section_vals_val_get
22 : USE kinds, ONLY: dp
23 : USE parallel_rng_types, ONLY: GAUSSIAN,&
24 : next_rng_seed,&
25 : rng_stream_type
26 : USE string_utilities, ONLY: compress
27 : #include "./base/base_uses.f90"
28 :
29 : IMPLICIT NONE
30 : PRIVATE
31 :
32 : PUBLIC :: gle_dealloc, &
33 : gle_init, gle_thermo_create, &
34 : gle_type
35 :
36 : !
37 : TYPE gle_thermo_type
38 : INTEGER :: degrees_of_freedom = -1
39 : REAL(KIND=dp) :: nkt = 0.0_dp, kin_energy = 0.0_dp, thermostat_energy = 0.0_dp
40 : REAL(KIND=dp), DIMENSION(:), POINTER :: s => NULL()
41 : TYPE(rng_stream_type) :: gaussian_rng_stream = rng_stream_type()
42 : END TYPE gle_thermo_type
43 :
44 : ! **************************************************************************************************
45 : TYPE gle_type
46 : INTEGER :: ndim = -1
47 : INTEGER :: glob_num_gle = -1, loc_num_gle = -1, region = -1
48 : INTEGER, DIMENSION(:), POINTER :: mal => NULL()
49 : REAL(dp) :: temp = 0.0_dp, dt = 0.0_dp, dt_fact = 0.0_dp
50 : REAL(dp), POINTER :: gle_s(:, :) => NULL(), gle_t(:, :) => NULL()
51 : REAL(dp), POINTER :: a_mat(:, :) => NULL(), c_mat(:, :) => NULL()
52 : TYPE(gle_thermo_type), POINTER :: nvt(:) => NULL()
53 : TYPE(map_info_type), POINTER :: map_info => NULL()
54 : END TYPE gle_type
55 :
56 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gle_system_types'
57 :
58 : CONTAINS
59 :
60 : ! **************************************************************************************************
61 : !> \brief ...
62 : !> \param gle ...
63 : !> \param dt ...
64 : !> \param temp ...
65 : !> \param section ...
66 : !> \param
67 : ! **************************************************************************************************
68 6 : SUBROUTINE gle_init(gle, dt, temp, section)
69 : TYPE(gle_type), POINTER :: gle
70 : REAL(dp), INTENT(IN) :: dt, temp
71 : TYPE(section_vals_type), POINTER :: section
72 :
73 : INTEGER :: i, ir, j, k, n_rep
74 6 : REAL(dp), DIMENSION(:), POINTER :: list
75 : REAL(KIND=dp) :: a_scale
76 :
77 6 : NULLIFY (gle%nvt)
78 6 : NULLIFY (gle%gle_s)
79 6 : NULLIFY (gle%gle_t)
80 6 : NULLIFY (gle%map_info)
81 6 : gle%loc_num_gle = 0
82 6 : gle%glob_num_gle = 0
83 6 : gle%temp = temp
84 6 : gle%dt = dt*0.5_dp
85 :
86 6 : CALL cite_reference(Ceriotti2009)
87 6 : CALL cite_reference(Ceriotti2009b)
88 6 : CALL section_vals_val_get(section, "NDIM", i_val=gle%ndim)
89 6 : CALL section_vals_val_get(section, "A_SCALE", r_val=a_scale)
90 :
91 24 : ALLOCATE (gle%a_mat(gle%ndim, gle%ndim))
92 24 : ALLOCATE (gle%c_mat(gle%ndim, gle%ndim))
93 24 : ALLOCATE (gle%gle_s(gle%ndim, gle%ndim))
94 24 : ALLOCATE (gle%gle_t(gle%ndim, gle%ndim))
95 :
96 6 : CALL section_vals_val_get(section, "A_LIST", n_rep_val=n_rep)
97 :
98 6 : j = 1
99 6 : k = 1
100 36 : DO ir = 1, n_rep
101 30 : NULLIFY (list)
102 : CALL section_vals_val_get(section, "A_LIST", &
103 30 : i_rep_val=ir, r_vals=list)
104 :
105 36 : IF (ASSOCIATED(list)) THEN
106 180 : DO i = 1, SIZE(list)
107 150 : IF (j > gle%ndim) THEN
108 0 : CPABORT("GLE: Too many elements in A_LIST")
109 : END IF
110 150 : gle%a_mat(j, k) = list(i)
111 150 : k = k + 1
112 180 : IF (k > gle%ndim) THEN
113 30 : k = 1
114 30 : j = j + 1
115 : END IF
116 : END DO
117 : END IF
118 : END DO ! ir
119 6 : IF (j < gle%ndim + 1) THEN
120 0 : CPABORT("GLE: Too few elements in A_LIST")
121 : END IF
122 186 : gle%a_mat = gle%a_mat*a_scale
123 :
124 6 : CALL section_vals_val_get(section, "C_LIST", n_rep_val=n_rep)
125 6 : IF (n_rep > 0) THEN
126 0 : j = 1
127 0 : k = 1
128 0 : DO ir = 1, n_rep
129 0 : NULLIFY (list)
130 : CALL section_vals_val_get(section, "C_LIST", &
131 0 : i_rep_val=ir, r_vals=list)
132 :
133 0 : IF (ASSOCIATED(list)) THEN
134 0 : DO i = 1, SIZE(list)
135 0 : IF (j > gle%ndim) THEN
136 0 : CPABORT("GLE: Too many elements in C_LIST")
137 : END IF
138 0 : gle%c_mat(j, k) = list(i)
139 0 : k = k + 1
140 0 : IF (k > gle%ndim) THEN
141 0 : k = 1
142 0 : j = j + 1
143 : END IF
144 : END DO
145 : END IF
146 : END DO ! ir
147 0 : IF (j < gle%ndim + 1) THEN
148 0 : CPABORT("GLE: Too few elements in C_LIST")
149 : END IF
150 : ELSE
151 186 : gle%c_mat = 0.0_dp
152 36 : DO i = 1, gle%ndim
153 36 : gle%c_mat(i, i) = gle%temp
154 : END DO
155 : END IF
156 6 : CALL create_map_info_type(gle%map_info)
157 24 : END SUBROUTINE gle_init
158 :
159 : ! **************************************************************************************************
160 : !> \brief ...
161 : !> \param gle ...
162 : !> \param mal_size ...
163 : !> \param
164 : ! **************************************************************************************************
165 6 : SUBROUTINE gle_thermo_create(gle, mal_size)
166 : TYPE(gle_type), POINTER :: gle
167 : INTEGER, INTENT(IN) :: mal_size
168 :
169 : CHARACTER(LEN=40) :: name
170 : INTEGER :: i, ithermo, my_index
171 6 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: seed
172 : REAL(KIND=dp), DIMENSION(3, 2) :: initial_seed, my_seed
173 :
174 6 : CPASSERT(ASSOCIATED(gle))
175 6 : CPASSERT(.NOT. ASSOCIATED(gle%nvt))
176 :
177 19254 : ALLOCATE (gle%nvt(gle%loc_num_gle))
178 19086 : DO i = 1, gle%loc_num_gle
179 19080 : NULLIFY (gle%nvt(i)%s)
180 57240 : ALLOCATE (gle%nvt(i)%s(gle%ndim))
181 19080 : gle%nvt(i)%kin_energy = 0.0_dp
182 19086 : gle%nvt(i)%thermostat_energy = 0.0_dp
183 : END DO
184 :
185 18 : ALLOCATE (gle%mal(mal_size))
186 18870 : gle%mal(:) = 0
187 :
188 : ! Initialize the gaussian stream random number
189 6 : initial_seed = next_rng_seed()
190 18 : ALLOCATE (seed(3, 2, gle%glob_num_gle))
191 :
192 54 : seed(:, :, 1) = initial_seed
193 19728 : DO ithermo = 2, gle%glob_num_gle
194 19728 : seed(:, :, ithermo) = next_rng_seed(seed(:, :, ithermo - 1))
195 : END DO
196 :
197 : ! Update initial seed
198 6 : initial_seed = next_rng_seed(seed(:, :, gle%glob_num_gle))
199 19086 : DO ithermo = 1, gle%loc_num_gle
200 19080 : my_index = gle%map_info%index(ithermo)
201 171720 : my_seed = seed(:, :, my_index)
202 19080 : WRITE (UNIT=name, FMT="(A,I8)") "Wiener process for Thermostat #", my_index
203 19080 : CALL compress(name)
204 : gle%nvt(ithermo)%gaussian_rng_stream = rng_stream_type( &
205 19086 : name=name, distribution_type=GAUSSIAN, extended_precision=.TRUE., seed=my_seed)
206 : END DO
207 :
208 6 : DEALLOCATE (seed)
209 :
210 6 : END SUBROUTINE gle_thermo_create
211 :
212 : ! **************************************************************************************************
213 : !> \brief Deallocate type for GLE thermostat
214 : !> \param gle ...
215 : ! **************************************************************************************************
216 6 : SUBROUTINE gle_dealloc(gle)
217 : TYPE(gle_type), POINTER :: gle
218 :
219 : INTEGER :: i
220 :
221 6 : IF (ASSOCIATED(gle)) THEN
222 6 : IF (ASSOCIATED(gle%a_mat)) THEN
223 6 : DEALLOCATE (gle%a_mat)
224 : END IF
225 6 : IF (ASSOCIATED(gle%c_mat)) THEN
226 6 : DEALLOCATE (gle%c_mat)
227 : END IF
228 6 : IF (ASSOCIATED(gle%gle_t)) THEN
229 6 : DEALLOCATE (gle%gle_t)
230 : END IF
231 6 : IF (ASSOCIATED(gle%gle_s)) THEN
232 6 : DEALLOCATE (gle%gle_s)
233 : END IF
234 6 : IF (ASSOCIATED(gle%nvt)) THEN
235 19086 : DO i = 1, SIZE(gle%nvt)
236 19086 : DEALLOCATE (gle%nvt(i)%s)
237 : END DO
238 6 : DEALLOCATE (gle%nvt)
239 : END IF
240 6 : IF (ASSOCIATED(gle%mal)) THEN
241 6 : DEALLOCATE (gle%mal)
242 : END IF
243 :
244 6 : CALL release_map_info_type(gle%map_info)
245 6 : DEALLOCATE (gle)
246 : END IF
247 :
248 6 : END SUBROUTINE gle_dealloc
249 :
250 0 : END MODULE gle_system_types
|
