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 0 : 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 0 : NULLIFY (new_logger, logger, enum, keyword, section)
94 0 : logger => cp_get_default_logger()
95 :
96 0 : CALL create_global_section(section)
97 0 : keyword => section_get_keyword(section, "RUN_TYPE")
98 0 : CALL keyword_get(keyword, enum=enum)
99 0 : label = TRIM(enum_i2c(enum, id_run))
100 0 : CALL section_release(section)
101 :
102 : ! Redirecting output of the sub_calculation to a different file
103 0 : CALL section_vals_val_get(root_section, "GLOBAL%PROJECT_NAME", c_val=project_name)
104 0 : input_file_path = TRIM(project_name)
105 0 : lp = LEN_TRIM(input_file_path)
106 0 : i = logger%iter_info%iteration(logger%iter_info%n_rlevel)
107 0 : input_file_path(lp + 1:LEN(input_file_path)) = "-"//TRIM(label)//"-"//ADJUSTL(cp_to_string(i))
108 0 : lp = LEN_TRIM(input_file_path)
109 0 : CALL section_vals_val_set(root_section, "GLOBAL%PROJECT_NAME", c_val=input_file_path(1:lp))
110 0 : CALL section_vals_val_set(root_section, "GLOBAL%RUN_TYPE", i_val=id_run)
111 :
112 : ! Redirecting output into a new file
113 0 : output_file_path = input_file_path(1:lp)//".out"
114 0 : IF (para_env%is_source()) THEN
115 : CALL open_file(file_name=output_file_path, file_status="UNKNOWN", &
116 0 : file_action="WRITE", file_position="APPEND", unit_number=unit_nr)
117 : ELSE
118 0 : unit_nr = -1
119 : END IF
120 : CALL cp_logger_create(new_logger, para_env=para_env, default_global_unit_nr=unit_nr, &
121 0 : close_global_unit_on_dealloc=.FALSE.)
122 0 : CALL section_vals_val_get(root_section, "GLOBAL%PROJECT", c_val=c_val)
123 0 : IF (c_val /= "") THEN
124 0 : CALL cp_logger_set(new_logger, local_filename=TRIM(c_val)//"_localLog")
125 : END IF
126 0 : new_logger%iter_info%project_name = TRIM(c_val)
127 : CALL section_vals_val_get(root_section, "GLOBAL%PRINT_LEVEL", &
128 0 : i_val=new_logger%iter_info%print_level)
129 :
130 0 : 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 0 : 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 0 : IF (para_env%is_source()) THEN
151 0 : unit_nr = cp_logger_get_default_unit_nr(new_logger)
152 0 : CALL close_file(unit_number=unit_nr)
153 : END IF
154 0 : CALL cp_logger_release(new_logger)
155 0 : CALL section_vals_val_set(root_section, "GLOBAL%RUN_TYPE", i_val=id_run)
156 0 : CALL section_vals_val_set(root_section, "GLOBAL%PROJECT_NAME", c_val=project_name)
157 :
158 0 : 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 204 : 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 204 : REAL(KIND=dp), DIMENSION(:), POINTER :: pvals
180 :
181 204 : NULLIFY (pvals)
182 204 : mtrx = 0.0_dp
183 204 : pres_ext_tens = 0.0_dp
184 204 : pres_ext = 0.0_dp
185 204 : CALL section_vals_val_get(geo_section, "EXTERNAL_PRESSURE", r_vals=pvals)
186 204 : check = (SIZE(pvals) == 1) .OR. (SIZE(pvals) == 9)
187 204 : IF (.NOT. check) &
188 0 : CPABORT("EXTERNAL_PRESSURE can have 1 or 9 components only!")
189 :
190 204 : IF (SIZE(pvals) == 9) THEN
191 : ind = 0
192 456 : DO i = 1, 3
193 1482 : DO j = 1, 3
194 1026 : ind = ind + 1
195 1368 : 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 4560 : pres_ext_tens = MATMUL(TRANSPOSE(rot), pres_ext_tens)
201 456 : DO i = 1, 3
202 456 : pres_ext = pres_ext + pres_ext_tens(i, i)
203 : END DO
204 114 : pres_ext = pres_ext/3.0_dp
205 456 : DO i = 1, 3
206 456 : pres_ext_tens(i, i) = pres_ext_tens(i, i) - pres_ext
207 : END DO
208 : ELSE
209 90 : pres_ext = pvals(1)
210 : END IF
211 :
212 2652 : 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 204 : 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 5054 : 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 5054 : my_keep_angles = .FALSE.
249 5054 : IF (PRESENT(keep_angles)) my_keep_angles = keep_angles
250 5054 : my_keep_symmetry = .FALSE.
251 5054 : IF (PRESENT(keep_symmetry)) my_keep_symmetry = keep_symmetry
252 35378 : gradient = 0.0_dp
253 5054 : IF (PRESENT(constraint_id)) THEN
254 5054 : my_constraint_id = constraint_id
255 : ELSE
256 0 : my_constraint_id = fix_none
257 : END IF
258 :
259 35378 : gradient = 0.0_dp
260 :
261 5054 : ptens = av_ptens
262 :
263 : ! Evaluating the internal pressure
264 5054 : pres_int = 0.0_dp
265 20216 : DO i = 1, 3
266 20216 : pres_int = pres_int + ptens(i, i)
267 : END DO
268 5054 : pres_int = pres_int/3.0_dp
269 :
270 0 : SELECT CASE (my_constraint_id)
271 : CASE (fix_x)
272 0 : pres_constr = 0.5_dp*(ptens(2, 2) + ptens(3, 3))
273 : CASE (fix_y)
274 0 : pres_constr = 0.5_dp*(ptens(1, 1) + ptens(3, 3))
275 : CASE (fix_z)
276 6 : pres_constr = 0.5_dp*(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 5054 : pres_constr = (ptens(1, 1) + ptens(2, 2) + ptens(3, 3))/3.0_dp
285 : END SELECT
286 :
287 5054 : ptens(1, 1) = av_ptens(1, 1) - pres_ext
288 5054 : ptens(2, 2) = av_ptens(2, 2) - pres_ext
289 5054 : ptens(3, 3) = av_ptens(3, 3) - pres_ext
290 :
291 262808 : pten_hinv_old = cell%deth*MATMUL(cell%h_inv, ptens)
292 202160 : correction = MATMUL(mtrx, cell%hmat)
293 :
294 5054 : gradient(1) = pten_hinv_old(1, 1) - correction(1, 1)
295 5054 : gradient(2) = pten_hinv_old(2, 1) - correction(2, 1)
296 5054 : gradient(3) = pten_hinv_old(2, 2) - correction(2, 2)
297 5054 : gradient(4) = pten_hinv_old(3, 1) - correction(3, 1)
298 5054 : gradient(5) = pten_hinv_old(3, 2) - correction(3, 2)
299 5054 : gradient(6) = pten_hinv_old(3, 3) - correction(3, 3)
300 :
301 5054 : CALL apply_cell_constraints(gradient, cell, my_keep_angles, my_keep_symmetry, my_constraint_id)
302 :
303 35378 : gradient = -gradient
304 :
305 5054 : END SUBROUTINE get_dg_dh
306 :
307 : ! **************************************************************************************************
308 : !> \brief Apply cell constraints
309 : !> \param gradient ...
310 : !> \param cell ...
311 : !> \param keep_angles ...
312 : !> \param keep_symmetry ...
313 : !> \param constraint_id ...
314 : !> \author Matthias Krack (October 26, 2017, MK)
315 : ! **************************************************************************************************
316 5054 : SUBROUTINE apply_cell_constraints(gradient, cell, keep_angles, keep_symmetry, constraint_id)
317 :
318 : REAL(KIND=dp), DIMENSION(:), POINTER :: gradient
319 : TYPE(cell_type), POINTER :: cell
320 : LOGICAL, INTENT(IN) :: keep_angles, keep_symmetry
321 : INTEGER, INTENT(IN) :: constraint_id
322 :
323 : REAL(KIND=dp) :: a, a_length, ab_length, b_length, cosa, &
324 : cosah, cosg, deriv_gamma, g, gamma, &
325 : norm, norm_b, norm_c, sina, sinah, sing
326 :
327 5054 : IF (keep_angles) THEN
328 : ! If we want to keep the angles constant we have to project out the
329 : ! components of the cell angles
330 2956 : norm_b = DOT_PRODUCT(cell%hmat(:, 2), cell%hmat(:, 2))
331 2217 : norm = DOT_PRODUCT(cell%hmat(1:2, 2), gradient(2:3))
332 3695 : gradient(2:3) = cell%hmat(1:2, 2)/norm_b*norm
333 2956 : norm_c = DOT_PRODUCT(cell%hmat(:, 3), cell%hmat(:, 3))
334 2956 : norm = DOT_PRODUCT(cell%hmat(1:3, 3), gradient(4:6))
335 5173 : gradient(4:6) = cell%hmat(1:3, 3)/norm_c*norm
336 : ! Retain an exact orthorhombic cell
337 : ! (off-diagonal elements must remain zero identically to keep QS fast)
338 739 : IF (cell%orthorhombic) THEN
339 157 : gradient(2) = 0.0_dp
340 157 : gradient(4) = 0.0_dp
341 157 : gradient(5) = 0.0_dp
342 : END IF
343 : END IF
344 :
345 5054 : IF (keep_symmetry) THEN
346 3083 : SELECT CASE (cell%symmetry_id)
347 : CASE (cell_sym_cubic, &
348 : cell_sym_tetragonal_ab, &
349 : cell_sym_tetragonal_ac, &
350 : cell_sym_tetragonal_bc, &
351 : cell_sym_orthorhombic)
352 100 : SELECT CASE (cell%symmetry_id)
353 : CASE (cell_sym_cubic)
354 100 : g = (gradient(1) + gradient(3) + gradient(6))/3.0_dp
355 100 : gradient(1) = g
356 100 : gradient(3) = g
357 100 : gradient(6) = g
358 : CASE (cell_sym_tetragonal_ab, &
359 : cell_sym_tetragonal_ac, &
360 : cell_sym_tetragonal_bc)
361 564 : SELECT CASE (cell%symmetry_id)
362 : CASE (cell_sym_tetragonal_ab)
363 186 : g = 0.5_dp*(gradient(1) + gradient(3))
364 186 : gradient(1) = g
365 186 : gradient(3) = g
366 : CASE (cell_sym_tetragonal_ac)
367 89 : g = 0.5_dp*(gradient(1) + gradient(6))
368 89 : gradient(1) = g
369 89 : gradient(6) = g
370 : CASE (cell_sym_tetragonal_bc)
371 86 : g = 0.5_dp*(gradient(3) + gradient(6))
372 86 : gradient(3) = g
373 361 : gradient(6) = g
374 : END SELECT
375 : CASE (cell_sym_orthorhombic)
376 : ! Nothing else to do
377 : END SELECT
378 564 : gradient(2) = 0.0_dp
379 564 : gradient(4) = 0.0_dp
380 564 : gradient(5) = 0.0_dp
381 : CASE (cell_sym_hexagonal_gamma_60)
382 236 : g = 0.5_dp*(gradient(1) + 0.5_dp*(gradient(2) + sqrt3*gradient(3)))
383 236 : gradient(1) = g
384 236 : gradient(2) = 0.5_dp*g
385 236 : gradient(3) = sqrt3*gradient(2)
386 236 : gradient(4) = 0.0_dp
387 236 : gradient(5) = 0.0_dp
388 : CASE (cell_sym_hexagonal_gamma_120)
389 118 : g = 0.5_dp*(gradient(1) - 0.5_dp*(gradient(2) - sqrt3*gradient(3)))
390 118 : gradient(1) = g
391 118 : gradient(2) = -0.5_dp*g
392 118 : gradient(3) = -sqrt3*gradient(2)
393 118 : gradient(4) = 0.0_dp
394 118 : gradient(5) = 0.0_dp
395 : CASE (cell_sym_rhombohedral)
396 : a = (angle(cell%hmat(:, 3), cell%hmat(:, 2)) + &
397 : angle(cell%hmat(:, 1), cell%hmat(:, 3)) + &
398 93 : angle(cell%hmat(:, 1), cell%hmat(:, 2)))/3.0_dp
399 93 : cosa = COS(a)
400 93 : sina = SIN(a)
401 93 : cosah = COS(0.5_dp*a)
402 93 : sinah = SIN(0.5_dp*a)
403 93 : norm = cosa/cosah
404 93 : norm_c = SQRT(1.0_dp - norm*norm)
405 : g = (gradient(1) + gradient(2)*cosa + gradient(3)*sina + &
406 93 : gradient(4)*cosah*norm + gradient(5)*sinah*norm + gradient(6)*norm_c)/3.0_dp
407 93 : gradient(1) = g
408 93 : gradient(2) = g*cosa
409 93 : gradient(3) = g*sina
410 93 : gradient(4) = g*cosah*norm
411 93 : gradient(5) = g*sinah*norm
412 93 : gradient(6) = g*norm_c
413 : CASE (cell_sym_monoclinic)
414 1178 : gradient(2) = 0.0_dp
415 1178 : gradient(5) = 0.0_dp
416 : CASE (cell_sym_monoclinic_gamma_ab)
417 : ! Cell symmetry with a = b, alpha = beta = 90 degree and gammma not equal 90 degree
418 : a_length = SQRT(cell%hmat(1, 1)*cell%hmat(1, 1) + &
419 : cell%hmat(2, 1)*cell%hmat(2, 1) + &
420 106 : cell%hmat(3, 1)*cell%hmat(3, 1))
421 : b_length = SQRT(cell%hmat(1, 2)*cell%hmat(1, 2) + &
422 : cell%hmat(2, 2)*cell%hmat(2, 2) + &
423 106 : cell%hmat(3, 2)*cell%hmat(3, 2))
424 106 : ab_length = 0.5_dp*(a_length + b_length)
425 106 : gamma = angle(cell%hmat(:, 1), cell%hmat(:, 2))
426 106 : cosg = COS(gamma)
427 106 : sing = SIN(gamma)
428 : ! Here, g is the average derivative of the cell vector length ab_length, and deriv_gamma is the derivative of the angle gamma
429 106 : g = 0.5_dp*(gradient(1) + cosg*gradient(2) + sing*gradient(3))
430 106 : deriv_gamma = (gradient(3)*cosg - gradient(2)*sing)/b_length
431 106 : gradient(1) = g
432 106 : gradient(2) = g*cosg - ab_length*sing*deriv_gamma
433 106 : gradient(3) = g*sing + ab_length*cosg*deriv_gamma
434 106 : gradient(4) = 0.0_dp
435 2519 : gradient(5) = 0.0_dp
436 : CASE (cell_sym_triclinic)
437 : ! Nothing to do
438 : END SELECT
439 : END IF
440 :
441 5054 : SELECT CASE (constraint_id)
442 : CASE (fix_x)
443 0 : gradient(1:2) = 0.0_dp
444 0 : gradient(4) = 0.0_dp
445 : CASE (fix_y)
446 0 : gradient(2:3) = 0.0_dp
447 0 : gradient(5) = 0.0_dp
448 : CASE (fix_z)
449 24 : gradient(4:6) = 0.0_dp
450 : CASE (fix_xy)
451 36 : gradient(1:5) = 0.0_dp
452 : CASE (fix_xz)
453 0 : gradient(1:2) = 0.0_dp
454 0 : gradient(4:6) = 0.0_dp
455 : CASE (fix_yz)
456 5054 : gradient(2:6) = 0.0_dp
457 : CASE (fix_none)
458 : ! Nothing to do
459 : END SELECT
460 :
461 5054 : END SUBROUTINE apply_cell_constraints
462 :
463 : ! **************************************************************************************************
464 : !> \brief Rescale x(idg+1:idg+6) according to vol
465 : !> \param vol ...
466 : !> \param x ...
467 : !> \param idg ...
468 : ! **************************************************************************************************
469 2278 : SUBROUTINE rescale_new_cell_volume(vol, x, idg)
470 :
471 : REAL(KIND=dp), INTENT(IN) :: vol
472 : REAL(KIND=dp), DIMENSION(:), POINTER :: x
473 : INTEGER, INTENT(IN) :: idg
474 :
475 : INTEGER :: i, j, my_idg
476 : REAL(KIND=dp) :: ratio
477 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat_tmp
478 :
479 2278 : CPASSERT((SIZE(x) == idg + 6))
480 2278 : my_idg = idg
481 :
482 2278 : hmat_tmp(:, :) = 0.0_dp
483 9112 : DO i = 1, 3
484 22780 : DO j = 1, i
485 13668 : my_idg = my_idg + 1
486 20502 : hmat_tmp(j, i) = x(my_idg)
487 : END DO
488 : END DO
489 :
490 2278 : ratio = vol/ABS(det_3x3(hmat_tmp))
491 2278 : ratio = ratio**(1.0_dp/3.0_dp) ! no need to use EXP(LOG(ratio)/3) for cube root
492 :
493 2278 : my_idg = idg
494 9112 : DO i = 1, 3
495 22780 : DO j = 1, i
496 13668 : my_idg = my_idg + 1
497 20502 : x(my_idg) = x(my_idg)*ratio
498 : END DO
499 : END DO
500 :
501 2278 : END SUBROUTINE rescale_new_cell_volume
502 :
503 : END MODULE cell_opt_utils
|
