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