Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2023 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief contains a functional that calculates the energy and its derivatives
10 : !> for the geometry optimizer
11 : !> \par History
12 : !> 03.2008 - Teodoro Laino [tlaino] - University of Zurich - Cell Optimization
13 : ! **************************************************************************************************
14 : MODULE cell_opt_utils
15 : USE cell_types, ONLY: &
16 : cell_sym_cubic, cell_sym_hexagonal, cell_sym_monoclinic, cell_sym_monoclinic_gamma_ab, &
17 : cell_sym_orthorhombic, cell_sym_rhombohedral, cell_sym_tetragonal_ab, &
18 : cell_sym_tetragonal_ac, cell_sym_tetragonal_bc, cell_sym_triclinic, cell_type, &
19 : get_cell_param, set_cell_param
20 : USE cp_files, ONLY: close_file,&
21 : open_file
22 : USE cp_log_handling, ONLY: cp_get_default_logger,&
23 : cp_logger_create,&
24 : cp_logger_get_default_unit_nr,&
25 : cp_logger_release,&
26 : cp_logger_set,&
27 : cp_logger_type,&
28 : cp_to_string
29 : USE cp_units, ONLY: cp_units_rad
30 : USE input_constants, ONLY: fix_none,&
31 : fix_x,&
32 : fix_xy,&
33 : fix_xz,&
34 : fix_y,&
35 : fix_yz,&
36 : fix_z
37 : USE input_cp2k_global, ONLY: create_global_section
38 : USE input_enumeration_types, ONLY: enum_i2c,&
39 : enumeration_type
40 : USE input_keyword_types, ONLY: keyword_get,&
41 : keyword_type
42 : USE input_section_types, ONLY: section_get_keyword,&
43 : section_release,&
44 : section_type,&
45 : section_vals_type,&
46 : section_vals_val_get,&
47 : section_vals_val_set
48 : USE kinds, ONLY: default_path_length,&
49 : default_string_length,&
50 : dp
51 : USE mathconstants, ONLY: sqrt3
52 : USE mathlib, ONLY: angle
53 : USE message_passing, ONLY: mp_para_env_type
54 : #include "../base/base_uses.f90"
55 :
56 : IMPLICIT NONE
57 : PRIVATE
58 :
59 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
60 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cell_opt_utils'
61 :
62 : PUBLIC :: get_ut_cell_matrix, get_dg_dh, gopt_new_logger_create, &
63 : gopt_new_logger_release, read_external_press_tensor, &
64 : apply_cell_constraints
65 :
66 : CONTAINS
67 : ! **************************************************************************************************
68 : !> \brief Transform a general cell matrix into an upper triangular one
69 : !> \param cell ...
70 : !> \author Teodoro Laino [tlaino] - University of Zurich - 03.2008
71 : ! **************************************************************************************************
72 192 : SUBROUTINE get_ut_cell_matrix(cell)
73 : TYPE(cell_type), POINTER :: cell
74 :
75 : INTEGER, DIMENSION(3) :: periodic
76 : REAL(KIND=dp), DIMENSION(3) :: cell_angle, cell_length
77 :
78 : ! orthorhombic cells are diagonal, and should stay exactly like this
79 :
80 192 : IF (.NOT. cell%orthorhombic) THEN
81 56 : CALL get_cell_param(cell, cell_length, cell_angle, cp_units_rad, periodic=periodic)
82 : CALL set_cell_param(cell, cell_length, cell_angle, periodic=periodic, &
83 56 : do_init_cell=.TRUE.)
84 : END IF
85 :
86 192 : END SUBROUTINE get_ut_cell_matrix
87 :
88 : ! **************************************************************************************************
89 : !> \brief creates a new logger used for cell optimization algorithm
90 : !> \param new_logger ...
91 : !> \param root_section ...
92 : !> \param para_env ...
93 : !> \param project_name ...
94 : !> \param id_run ...
95 : !> \author Teodoro Laino [tlaino] - University of Zurich - 03.2008
96 : ! **************************************************************************************************
97 850 : SUBROUTINE gopt_new_logger_create(new_logger, root_section, para_env, project_name, &
98 : id_run)
99 : TYPE(cp_logger_type), POINTER :: new_logger
100 : TYPE(section_vals_type), POINTER :: root_section
101 : TYPE(mp_para_env_type), POINTER :: para_env
102 : CHARACTER(len=default_string_length), INTENT(OUT) :: project_name
103 : INTEGER, INTENT(IN) :: id_run
104 :
105 : CHARACTER(len=default_path_length) :: c_val, input_file_path, output_file_path
106 : CHARACTER(len=default_string_length) :: label
107 : INTEGER :: i, lp, unit_nr
108 : TYPE(cp_logger_type), POINTER :: logger
109 : TYPE(enumeration_type), POINTER :: enum
110 : TYPE(keyword_type), POINTER :: keyword
111 : TYPE(section_type), POINTER :: section
112 :
113 850 : NULLIFY (new_logger, logger, enum, keyword, section)
114 850 : logger => cp_get_default_logger()
115 :
116 850 : CALL create_global_section(section)
117 850 : keyword => section_get_keyword(section, "RUN_TYPE")
118 850 : CALL keyword_get(keyword, enum=enum)
119 850 : label = TRIM(enum_i2c(enum, id_run))
120 850 : CALL section_release(section)
121 :
122 : ! Redirecting output of the sub_calculation to a different file
123 850 : CALL section_vals_val_get(root_section, "GLOBAL%PROJECT_NAME", c_val=project_name)
124 850 : input_file_path = TRIM(project_name)
125 850 : lp = LEN_TRIM(input_file_path)
126 850 : i = logger%iter_info%iteration(logger%iter_info%n_rlevel)
127 850 : input_file_path(lp + 1:LEN(input_file_path)) = "-"//TRIM(label)//"-"//ADJUSTL(cp_to_string(i))
128 850 : lp = LEN_TRIM(input_file_path)
129 850 : CALL section_vals_val_set(root_section, "GLOBAL%PROJECT_NAME", c_val=input_file_path(1:lp))
130 850 : CALL section_vals_val_set(root_section, "GLOBAL%RUN_TYPE", i_val=id_run)
131 :
132 : ! Redirecting output into a new file
133 850 : output_file_path = input_file_path(1:lp)//".out"
134 850 : IF (para_env%is_source()) THEN
135 : CALL open_file(file_name=output_file_path, file_status="UNKNOWN", &
136 425 : file_action="WRITE", file_position="APPEND", unit_number=unit_nr)
137 : ELSE
138 425 : unit_nr = -1
139 : END IF
140 : CALL cp_logger_create(new_logger, para_env=para_env, default_global_unit_nr=unit_nr, &
141 850 : close_global_unit_on_dealloc=.FALSE.)
142 850 : CALL section_vals_val_get(root_section, "GLOBAL%PROJECT", c_val=c_val)
143 850 : IF (c_val /= "") THEN
144 850 : CALL cp_logger_set(new_logger, local_filename=TRIM(c_val)//"_localLog")
145 : END IF
146 850 : new_logger%iter_info%project_name = c_val
147 : CALL section_vals_val_get(root_section, "GLOBAL%PRINT_LEVEL", &
148 850 : i_val=new_logger%iter_info%print_level)
149 :
150 850 : END SUBROUTINE gopt_new_logger_create
151 :
152 : ! **************************************************************************************************
153 : !> \brief releases a new logger used for cell optimization algorithm
154 : !> \param new_logger ...
155 : !> \param root_section ...
156 : !> \param para_env ...
157 : !> \param project_name ...
158 : !> \param id_run ...
159 : !> \author Teodoro Laino [tlaino] - University of Zurich - 03.2008
160 : ! **************************************************************************************************
161 850 : SUBROUTINE gopt_new_logger_release(new_logger, root_section, para_env, project_name, id_run)
162 : TYPE(cp_logger_type), POINTER :: new_logger
163 : TYPE(section_vals_type), POINTER :: root_section
164 : TYPE(mp_para_env_type), POINTER :: para_env
165 : CHARACTER(len=default_string_length), INTENT(IN) :: project_name
166 : INTEGER, INTENT(IN) :: id_run
167 :
168 : INTEGER :: unit_nr
169 :
170 850 : IF (para_env%is_source()) THEN
171 425 : unit_nr = cp_logger_get_default_unit_nr(new_logger)
172 425 : CALL close_file(unit_number=unit_nr)
173 : END IF
174 850 : CALL cp_logger_release(new_logger)
175 850 : CALL section_vals_val_set(root_section, "GLOBAL%RUN_TYPE", i_val=id_run)
176 850 : CALL section_vals_val_set(root_section, "GLOBAL%PROJECT_NAME", c_val=project_name)
177 :
178 850 : END SUBROUTINE gopt_new_logger_release
179 :
180 : ! **************************************************************************************************
181 : !> \brief Reads the external pressure tensor
182 : !> \param geo_section ...
183 : !> \param cell ...
184 : !> \param pres_ext ...
185 : !> \param mtrx ...
186 : !> \param rot ...
187 : !> \author Teodoro Laino [tlaino] - University of Zurich - 03.2008
188 : ! **************************************************************************************************
189 192 : SUBROUTINE read_external_press_tensor(geo_section, cell, pres_ext, mtrx, rot)
190 : TYPE(section_vals_type), POINTER :: geo_section
191 : TYPE(cell_type), POINTER :: cell
192 : REAL(KIND=dp), INTENT(OUT) :: pres_ext
193 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT) :: mtrx
194 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: rot
195 :
196 : INTEGER :: i, ind, j
197 : LOGICAL :: check
198 : REAL(KIND=dp), DIMENSION(3, 3) :: pres_ext_tens
199 192 : REAL(KIND=dp), DIMENSION(:), POINTER :: pvals
200 :
201 192 : NULLIFY (pvals)
202 192 : mtrx = 0.0_dp
203 192 : pres_ext_tens = 0.0_dp
204 192 : pres_ext = 0.0_dp
205 192 : CALL section_vals_val_get(geo_section, "EXTERNAL_PRESSURE", r_vals=pvals)
206 192 : check = (SIZE(pvals) == 1) .OR. (SIZE(pvals) == 9)
207 192 : IF (.NOT. check) &
208 0 : CPABORT("EXTERNAL_PRESSURE can have 1 or 9 components only!")
209 :
210 192 : IF (SIZE(pvals) == 9) THEN
211 : ind = 0
212 392 : DO i = 1, 3
213 1274 : DO j = 1, 3
214 882 : ind = ind + 1
215 1176 : pres_ext_tens(j, i) = pvals(ind)
216 : END DO
217 : END DO
218 : ! Also the pressure tensor must be oriented in the same canonical directions
219 : ! of the simulation cell
220 3920 : pres_ext_tens = MATMUL(TRANSPOSE(rot), pres_ext_tens)
221 392 : DO i = 1, 3
222 392 : pres_ext = pres_ext + pres_ext_tens(i, i)
223 : END DO
224 98 : pres_ext = pres_ext/3.0_dp
225 392 : DO i = 1, 3
226 392 : pres_ext_tens(i, i) = pres_ext_tens(i, i) - pres_ext
227 : END DO
228 : ELSE
229 94 : pres_ext = pvals(1)
230 : END IF
231 :
232 2496 : IF (ANY(pres_ext_tens > 1.0E-5_dp)) THEN
233 0 : mtrx = cell%deth*MATMUL(cell%h_inv, MATMUL(pres_ext_tens, TRANSPOSE(cell%h_inv)))
234 : END IF
235 :
236 192 : END SUBROUTINE read_external_press_tensor
237 :
238 : ! **************************************************************************************************
239 : !> \brief Computes the derivatives for the cell
240 : !> \param gradient ...
241 : !> \param av_ptens ...
242 : !> \param pres_ext ...
243 : !> \param cell ...
244 : !> \param mtrx ...
245 : !> \param keep_angles ...
246 : !> \param keep_symmetry ...
247 : !> \param pres_int ...
248 : !> \param pres_constr ...
249 : !> \param constraint_id ...
250 : !> \author Teodoro Laino [tlaino] - University of Zurich - 03.2008
251 : ! **************************************************************************************************
252 3377 : SUBROUTINE get_dg_dh(gradient, av_ptens, pres_ext, cell, mtrx, keep_angles, &
253 : keep_symmetry, pres_int, pres_constr, constraint_id)
254 :
255 : REAL(KIND=dp), DIMENSION(:), POINTER :: gradient
256 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: av_ptens
257 : REAL(KIND=dp), INTENT(IN) :: pres_ext
258 : TYPE(cell_type), POINTER :: cell
259 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: mtrx
260 : LOGICAL, INTENT(IN), OPTIONAL :: keep_angles, keep_symmetry
261 : REAL(KIND=dp), INTENT(OUT) :: pres_int, pres_constr
262 : INTEGER, INTENT(IN), OPTIONAL :: constraint_id
263 :
264 : INTEGER :: i, my_constraint_id
265 : LOGICAL :: my_keep_angles, my_keep_symmetry
266 : REAL(KIND=dp), DIMENSION(3, 3) :: correction, pten_hinv_old, ptens
267 :
268 3377 : my_keep_angles = .FALSE.
269 3377 : IF (PRESENT(keep_angles)) my_keep_angles = keep_angles
270 3377 : my_keep_symmetry = .FALSE.
271 3377 : IF (PRESENT(keep_symmetry)) my_keep_symmetry = keep_symmetry
272 23639 : gradient = 0.0_dp
273 3377 : IF (PRESENT(constraint_id)) THEN
274 3377 : my_constraint_id = constraint_id
275 : ELSE
276 0 : my_constraint_id = fix_none
277 : END IF
278 :
279 23639 : gradient = 0.0_dp
280 :
281 3377 : ptens = av_ptens
282 :
283 : ! Evaluating the internal pressure
284 3377 : pres_int = 0.0_dp
285 13508 : DO i = 1, 3
286 13508 : pres_int = pres_int + ptens(i, i)
287 : END DO
288 3377 : pres_int = pres_int/3.0_dp
289 :
290 0 : SELECT CASE (my_constraint_id)
291 : CASE (fix_x)
292 0 : pres_constr = ptens(2, 2) + ptens(3, 3)
293 : CASE (fix_y)
294 0 : pres_constr = ptens(1, 1) + ptens(3, 3)
295 : CASE (fix_z)
296 5 : pres_constr = ptens(1, 1) + ptens(2, 2)
297 : CASE (fix_xy)
298 10 : pres_constr = ptens(3, 3)
299 : CASE (fix_xz)
300 0 : pres_constr = ptens(2, 2)
301 : CASE (fix_yz)
302 0 : pres_constr = ptens(1, 1)
303 : CASE (fix_none)
304 3377 : pres_constr = ptens(1, 1) + ptens(2, 2) + ptens(3, 3)
305 : END SELECT
306 3377 : pres_constr = pres_constr/3.0_dp
307 :
308 3377 : ptens(1, 1) = av_ptens(1, 1) - pres_ext
309 3377 : ptens(2, 2) = av_ptens(2, 2) - pres_ext
310 3377 : ptens(3, 3) = av_ptens(3, 3) - pres_ext
311 :
312 175604 : pten_hinv_old = cell%deth*MATMUL(cell%h_inv, ptens)
313 135080 : correction = MATMUL(mtrx, cell%hmat)
314 :
315 3377 : gradient(1) = pten_hinv_old(1, 1) - correction(1, 1)
316 3377 : gradient(2) = pten_hinv_old(2, 1) - correction(2, 1)
317 3377 : gradient(3) = pten_hinv_old(2, 2) - correction(2, 2)
318 3377 : gradient(4) = pten_hinv_old(3, 1) - correction(3, 1)
319 3377 : gradient(5) = pten_hinv_old(3, 2) - correction(3, 2)
320 3377 : gradient(6) = pten_hinv_old(3, 3) - correction(3, 3)
321 :
322 3377 : CALL apply_cell_constraints(gradient, cell, my_keep_angles, my_keep_symmetry, my_constraint_id)
323 :
324 23639 : gradient = -gradient
325 :
326 3377 : END SUBROUTINE get_dg_dh
327 :
328 : ! **************************************************************************************************
329 : !> \brief Apply cell constraints
330 : !> \param gradient ...
331 : !> \param cell ...
332 : !> \param keep_angles ...
333 : !> \param keep_symmetry ...
334 : !> \param constraint_id ...
335 : !> \author Matthias Krack (October 26, 2017, MK)
336 : ! **************************************************************************************************
337 3377 : SUBROUTINE apply_cell_constraints(gradient, cell, keep_angles, keep_symmetry, constraint_id)
338 :
339 : REAL(KIND=dp), DIMENSION(:), POINTER :: gradient
340 : TYPE(cell_type), POINTER :: cell
341 : LOGICAL, INTENT(IN) :: keep_angles, keep_symmetry
342 : INTEGER, INTENT(IN) :: constraint_id
343 :
344 : REAL(KIND=dp) :: a, a_length, ab_length, b_length, cosa, &
345 : cosah, cosgamma, deriv_gamma, g, &
346 : gamma, norm, norm_b, norm_c, sina, &
347 : sinah, singamma
348 :
349 3377 : IF (keep_angles) THEN
350 : ! If we want to keep the angles constant we have to project out the
351 : ! components of the cell angles
352 400 : norm_b = DOT_PRODUCT(cell%hmat(:, 2), cell%hmat(:, 2))
353 300 : norm = DOT_PRODUCT(cell%hmat(1:2, 2), gradient(2:3))
354 500 : gradient(2:3) = cell%hmat(1:2, 2)/norm_b*norm
355 400 : norm_c = DOT_PRODUCT(cell%hmat(:, 3), cell%hmat(:, 3))
356 400 : norm = DOT_PRODUCT(cell%hmat(1:3, 3), gradient(4:6))
357 700 : gradient(4:6) = cell%hmat(1:3, 3)/norm_c*norm
358 : ! Retain an exact orthorhombic cell
359 : ! (off-diagonal elements must remain zero identically to keep QS fast)
360 100 : IF (cell%orthorhombic) THEN
361 70 : gradient(2) = 0.0_dp
362 70 : gradient(4) = 0.0_dp
363 70 : gradient(5) = 0.0_dp
364 : END IF
365 : END IF
366 :
367 3377 : IF (keep_symmetry) THEN
368 568 : SELECT CASE (cell%symmetry_id)
369 : CASE (cell_sym_cubic, &
370 : cell_sym_tetragonal_ab, &
371 : cell_sym_tetragonal_ac, &
372 : cell_sym_tetragonal_bc, &
373 : cell_sym_orthorhombic)
374 100 : SELECT CASE (cell%symmetry_id)
375 : CASE (cell_sym_cubic)
376 100 : g = (gradient(1) + gradient(3) + gradient(6))/3.0_dp
377 100 : gradient(1) = g
378 100 : gradient(3) = g
379 465 : gradient(6) = g
380 : CASE (cell_sym_tetragonal_ab, &
381 : cell_sym_tetragonal_ac, &
382 : cell_sym_tetragonal_bc)
383 186 : SELECT CASE (cell%symmetry_id)
384 : CASE (cell_sym_tetragonal_ab)
385 186 : g = 0.5_dp*(gradient(1) + gradient(3))
386 186 : gradient(1) = g
387 186 : gradient(3) = g
388 : CASE (cell_sym_tetragonal_ac)
389 93 : g = 0.5_dp*(gradient(1) + gradient(6))
390 93 : gradient(1) = g
391 93 : gradient(6) = g
392 : CASE (cell_sym_tetragonal_bc)
393 86 : g = 0.5_dp*(gradient(3) + gradient(6))
394 86 : gradient(3) = g
395 86 : gradient(6) = g
396 : END SELECT
397 : CASE (cell_sym_orthorhombic)
398 : ! Nothing else to do
399 : END SELECT
400 568 : gradient(2) = 0.0_dp
401 568 : gradient(4) = 0.0_dp
402 568 : gradient(5) = 0.0_dp
403 : CASE (cell_sym_hexagonal)
404 118 : g = 0.5_dp*(gradient(1) + 0.5_dp*(gradient(2) + sqrt3*gradient(3)))
405 118 : gradient(1) = g
406 118 : gradient(2) = 0.5_dp*g
407 118 : gradient(3) = sqrt3*gradient(2)
408 118 : gradient(4) = 0.0_dp
409 118 : gradient(5) = 0.0_dp
410 : CASE (cell_sym_rhombohedral)
411 : a = (angle(cell%hmat(:, 3), cell%hmat(:, 2)) + &
412 : angle(cell%hmat(:, 1), cell%hmat(:, 3)) + &
413 83 : angle(cell%hmat(:, 1), cell%hmat(:, 2)))/3.0_dp
414 83 : cosa = COS(a)
415 83 : sina = SIN(a)
416 83 : cosah = COS(0.5_dp*a)
417 83 : sinah = SIN(0.5_dp*a)
418 83 : norm = cosa/cosah
419 83 : norm_c = SQRT(1.0_dp - norm*norm)
420 : g = (gradient(1) + gradient(2)*cosa + gradient(3)*sina + &
421 83 : gradient(4)*cosah*norm + gradient(5)*sinah*norm + gradient(6)*norm_c)/3.0_dp
422 83 : gradient(1) = g
423 83 : gradient(2) = g*cosa
424 83 : gradient(3) = g*sina
425 83 : gradient(4) = g*cosah*norm
426 83 : gradient(5) = g*sinah*norm
427 83 : gradient(6) = g*norm_c
428 : CASE (cell_sym_monoclinic)
429 126 : gradient(2) = 0.0_dp
430 126 : gradient(5) = 0.0_dp
431 : CASE (cell_sym_monoclinic_gamma_ab)
432 : ! Cell symmetry with a=b, alpha=beta=90deg and gammma unequal 90deg
433 : a_length = SQRT(cell%hmat(1, 1)*cell%hmat(1, 1) + &
434 : cell%hmat(2, 1)*cell%hmat(2, 1) + &
435 106 : cell%hmat(3, 1)*cell%hmat(3, 1))
436 : b_length = SQRT(cell%hmat(1, 2)*cell%hmat(1, 2) + &
437 : cell%hmat(2, 2)*cell%hmat(2, 2) + &
438 106 : cell%hmat(3, 2)*cell%hmat(3, 2))
439 106 : ab_length = 0.5_dp*(a_length + b_length)
440 106 : gamma = angle(cell%hmat(:, 1), cell%hmat(:, 2))
441 106 : cosgamma = COS(gamma)
442 106 : singamma = SIN(gamma)
443 : ! Here, g is the average derivative of the cell vector length ab_length, and deriv_gamma is the derivative of the angle gamma
444 106 : g = 0.5_dp*(gradient(1) + cosgamma*gradient(2) + singamma*gradient(3))
445 106 : deriv_gamma = (gradient(3)*cosgamma - gradient(2)*singamma)/b_length
446 106 : gradient(1) = g
447 106 : gradient(2) = g*cosgamma - ab_length*singamma*deriv_gamma
448 106 : gradient(3) = g*singamma + ab_length*cosgamma*deriv_gamma
449 106 : gradient(4) = 0.0_dp
450 1223 : gradient(5) = 0.0_dp
451 : CASE (cell_sym_triclinic)
452 : ! Nothing to do
453 : END SELECT
454 : END IF
455 :
456 3377 : SELECT CASE (constraint_id)
457 : CASE (fix_x)
458 0 : gradient(1:2) = 0.0_dp
459 0 : gradient(4) = 0.0_dp
460 : CASE (fix_y)
461 0 : gradient(2:3) = 0.0_dp
462 0 : gradient(5) = 0.0_dp
463 : CASE (fix_z)
464 20 : gradient(4:6) = 0.0_dp
465 : CASE (fix_xy)
466 60 : gradient(1:5) = 0.0_dp
467 : CASE (fix_xz)
468 0 : gradient(1:2) = 0.0_dp
469 0 : gradient(4:6) = 0.0_dp
470 : CASE (fix_yz)
471 3377 : gradient(2:6) = 0.0_dp
472 : CASE (fix_none)
473 : ! Nothing to do
474 : END SELECT
475 :
476 3377 : END SUBROUTINE apply_cell_constraints
477 :
478 : END MODULE cell_opt_utils
|