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 input_constants, ONLY: fix_none
21 : USE input_section_types, ONLY: section_vals_type,&
22 : section_vals_val_get
23 : USE kinds, ONLY: default_string_length,&
24 : dp
25 : USE mathconstants, ONLY: sqrt3
26 : USE mathlib, ONLY: angle,&
27 : det_3x3
28 : USE string_utilities, ONLY: uppercase
29 : #include "../base/base_uses.f90"
30 :
31 : IMPLICIT NONE
32 : PRIVATE
33 :
34 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cell_opt_utils'
35 :
36 : PUBLIC :: get_dg_dh, keep_angles_conflict, &
37 : keep_symmetry_coupled_set, keep_symmetry_conflict, &
38 : pres_crit_shear_pinned, parse_token, &
39 : read_external_press_tensor, &
40 : apply_cell_constraints, cell_opt_constraint_label, &
41 : parse_cell_constraint, &
42 : rescale_new_cell_volume
43 :
44 : ! Masks of the six independent cell variables in the canonical frame
45 : ! (a along x, b in the xy plane): AX, BX, BY, CX, CY, CZ.
46 : ! A set bit means that a CONSTRAINT fixes the component.
47 : ! constraint_id = 0 (fix_none) leaves every component relaxable.
48 : INTEGER, PARAMETER, PUBLIC :: comp_ax = 2**0, & ! a(x)
49 : comp_bx = 2**1, & ! b(x)
50 : comp_by = 2**2, & ! b(y)
51 : comp_cx = 2**3, & ! c(x)
52 : comp_cy = 2**4, & ! c(y)
53 : comp_cz = 2**5 ! c(z)
54 : ! Shorthands for complete cell vectors, Cartesian directions and coupled
55 : ! component sets (unions of disjoint bits, hence representable as sums)
56 : INTEGER, PARAMETER, PUBLIC :: mask_a = comp_ax, &
57 : mask_b = comp_bx + comp_by, &
58 : mask_c = comp_cx + comp_cy + comp_cz, &
59 : mask_x = comp_ax + comp_bx + comp_cx, &
60 : mask_y = comp_bx + comp_by + comp_cy, &
61 : mask_z = mask_c, &
62 : mask_xy = comp_ax + comp_bx + comp_by + comp_cx + comp_cy, &
63 : mask_xz = comp_ax + comp_bx + comp_cx + comp_cy + comp_cz, &
64 : mask_yz = comp_bx + comp_by + comp_cx + comp_cy + comp_cz, &
65 : mask_diag = comp_ax + comp_by + comp_cz, &
66 : mask_diag_ab = comp_ax + comp_by, &
67 : mask_diag_ac = comp_ax + comp_cz, &
68 : mask_diag_bc = comp_by + comp_cz, &
69 : mask_ab_plane = comp_ax + comp_bx + comp_by, &
70 : mask_shear = comp_bx + comp_cx + comp_cy, &
71 : mask_all = comp_ax + comp_bx + comp_by + comp_cx + comp_cy + comp_cz
72 :
73 : CONTAINS
74 :
75 : ! **************************************************************************************************
76 : !> \brief Reads the external pressure tensor
77 : !> \param geo_section ...
78 : !> \param cell ...
79 : !> \param pres_ext ...
80 : !> \param mtrx ...
81 : !> \param input_to_opt rotation from the input Cartesian frame to the optimization frame
82 : !> \author Teodoro Laino [tlaino] - University of Zurich - 03.2008
83 : ! **************************************************************************************************
84 234 : SUBROUTINE read_external_press_tensor(geo_section, cell, pres_ext, mtrx, input_to_opt)
85 : TYPE(section_vals_type), POINTER :: geo_section
86 : TYPE(cell_type), POINTER :: cell
87 : REAL(KIND=dp), INTENT(OUT) :: pres_ext
88 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT) :: mtrx
89 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: input_to_opt
90 :
91 : INTEGER :: i, ind, j
92 : LOGICAL :: check
93 : REAL(KIND=dp), DIMENSION(3, 3) :: pres_ext_tens
94 234 : REAL(KIND=dp), DIMENSION(:), POINTER :: pvals
95 :
96 234 : NULLIFY (pvals)
97 234 : mtrx = 0.0_dp
98 234 : pres_ext_tens = 0.0_dp
99 234 : pres_ext = 0.0_dp
100 234 : CALL section_vals_val_get(geo_section, "EXTERNAL_PRESSURE", r_vals=pvals)
101 234 : check = (SIZE(pvals) == 1) .OR. (SIZE(pvals) == 9)
102 234 : IF (.NOT. check) THEN
103 0 : CPABORT("EXTERNAL_PRESSURE can have 1 or 9 components only!")
104 : END IF
105 :
106 234 : IF (SIZE(pvals) == 9) THEN
107 : ind = 0
108 472 : DO i = 1, 3
109 1534 : DO j = 1, 3
110 1062 : ind = ind + 1
111 1416 : pres_ext_tens(j, i) = pvals(ind)
112 : END DO
113 : END DO
114 : ! Express the pressure tensor in the canonical optimization frame.
115 : pres_ext_tens = MATMUL(input_to_opt, &
116 9322 : MATMUL(pres_ext_tens, TRANSPOSE(input_to_opt)))
117 472 : DO i = 1, 3
118 472 : pres_ext = pres_ext + pres_ext_tens(i, i)
119 : END DO
120 118 : pres_ext = pres_ext/3.0_dp
121 472 : DO i = 1, 3
122 472 : pres_ext_tens(i, i) = pres_ext_tens(i, i) - pres_ext
123 : END DO
124 : ELSE
125 116 : pres_ext = pvals(1)
126 : END IF
127 :
128 3042 : IF (ANY(ABS(pres_ext_tens) > 1.0E-5_dp)) THEN
129 0 : mtrx = cell%deth*MATMUL(cell%h_inv, MATMUL(pres_ext_tens, TRANSPOSE(cell%h_inv)))
130 : END IF
131 :
132 234 : END SUBROUTINE read_external_press_tensor
133 :
134 : ! **************************************************************************************************
135 : !> \brief Parses the tokens of the CONSTRAINT keyword into a cell component mask
136 : !> \param constraint_tokens value tokens of the CONSTRAINT keyword (char_t list)
137 : !> \param cell simulation cell (only the periodicity is used)
138 : !> \param explicit whether the CONSTRAINT keyword was given explicitly
139 : !> \param mask union of the fixed cell components (comp_* and mask_* bits)
140 : ! **************************************************************************************************
141 422 : SUBROUTINE parse_cell_constraint(constraint_tokens, cell, explicit, mask)
142 : CHARACTER(LEN=*), DIMENSION(:), INTENT(IN) :: constraint_tokens
143 : TYPE(cell_type), POINTER :: cell
144 : LOGICAL, INTENT(IN) :: explicit
145 : INTEGER, INTENT(OUT) :: mask
146 :
147 : CHARACTER(LEN=default_string_length) :: reject_reason, token
148 : INTEGER :: axis_index, nonperiodic_axis, overlap, &
149 : required_mask, tok_index, tok_mask
150 : LOGICAL :: accepted
151 :
152 : ! The input tokenizer delivers a bare keyword as one blank token.
153 422 : IF (SIZE(constraint_tokens) == 0 .OR. &
154 422 : (SIZE(constraint_tokens) == 1 .AND. LEN_TRIM(constraint_tokens(1)) == 0)) THEN
155 : CALL cp_abort(__LOCATION__, &
156 0 : "CELL_OPT CONSTRAINT requires a value.")
157 : END IF
158 :
159 422 : mask = 0
160 1094 : DO tok_index = 1, SIZE(constraint_tokens)
161 672 : token = ADJUSTL(constraint_tokens(tok_index))
162 672 : CALL uppercase(token)
163 672 : IF (token == "NONE") THEN
164 : ! NONE must be the first and only token: leftover tokens on the
165 : ! line or an already-accumulated mask mean it is being combined.
166 224 : IF (tok_index < SIZE(constraint_tokens) .OR. mask /= 0) THEN
167 : CALL cp_abort(__LOCATION__, &
168 : "CELL_OPT CONSTRAINT NONE cannot be combined "// &
169 0 : "with other tokens.")
170 : END IF
171 224 : mask = 0
172 224 : CYCLE
173 : END IF
174 448 : CALL parse_token(token, accepted, tok_mask, reject_reason)
175 448 : IF (.NOT. accepted) THEN
176 0 : CALL cp_abort(__LOCATION__, TRIM(reject_reason))
177 : END IF
178 448 : overlap = IAND(tok_mask, mask)
179 448 : IF (overlap == tok_mask) THEN
180 : CALL cp_warn(__LOCATION__, &
181 : "Redundant CELL_OPT CONSTRAINT token: "//TRIM(token)// &
182 : " (already constrained: "// &
183 2 : TRIM(ADJUSTL(cell_opt_constraint_label(overlap)))//").")
184 446 : ELSE IF (overlap /= 0) THEN
185 : CALL cp_warn(__LOCATION__, &
186 : "CELL_OPT CONSTRAINT token "//TRIM(token)//" partially overlaps: "// &
187 : "already constrained "// &
188 : TRIM(ADJUSTL(cell_opt_constraint_label(overlap)))// &
189 : ", newly fixed "// &
190 2 : TRIM(ADJUSTL(cell_opt_constraint_label(IAND(NOT(mask), tok_mask))))//".")
191 : END IF
192 1318 : mask = IOR(mask, tok_mask)
193 : END DO
194 :
195 : ! For 2D periodicity the non-periodic cell direction must be constrained.
196 : ! Complete the constraint automatically if the keyword was not given.
197 1688 : IF (COUNT(cell%perd /= 0) == 2) THEN
198 : nonperiodic_axis = 0
199 80 : DO axis_index = 1, 3
200 80 : IF (cell%perd(axis_index) == 0) nonperiodic_axis = axis_index
201 : END DO
202 20 : required_mask = 0
203 6 : SELECT CASE (nonperiodic_axis)
204 : CASE (1)
205 6 : required_mask = mask_x
206 : CASE (2)
207 6 : required_mask = mask_y
208 : CASE (3)
209 20 : required_mask = mask_z
210 : END SELECT
211 20 : IF (.NOT. explicit .AND. mask == 0) THEN
212 8 : mask = required_mask
213 : CALL cp_warn(__LOCATION__, &
214 8 : "2D CELL_OPT: constraining non-periodic cell direction.")
215 12 : ELSE IF (IAND(mask, required_mask) /= required_mask) THEN
216 : CALL cp_abort(__LOCATION__, &
217 0 : "2D CELL_OPT needs constrained non-periodic cell direction.")
218 : END IF
219 : END IF
220 :
221 422 : END SUBROUTINE parse_cell_constraint
222 :
223 : ! **************************************************************************************************
224 : !> \brief Returns the label of the consolidated cell constraint mask
225 : !> \param mask bit mask of fixed cell components (comp_* / mask_* bits)
226 : !> \return the historical four-character labels (" X", " XY", "NONE") for
227 : !> the legacy constraint values (padded for backward compatible
228 : !> printing), the fixed components spelled out ("AX BX BY") otherwise
229 : ! **************************************************************************************************
230 503 : PURE FUNCTION cell_opt_constraint_label(mask) RESULT(label)
231 : INTEGER, INTENT(IN) :: mask
232 : CHARACTER(LEN=default_string_length) :: label
233 :
234 822 : SELECT CASE (mask)
235 : CASE (fix_none)
236 319 : label = "NONE"
237 : CASE (mask_x)
238 4 : label = " X"
239 : CASE (mask_y)
240 4 : label = " Y"
241 : CASE (mask_z)
242 18 : label = " Z"
243 : CASE (mask_xy)
244 7 : label = " XY"
245 : CASE (mask_xz)
246 4 : label = " XZ"
247 : CASE (mask_yz)
248 4 : label = " YZ"
249 : CASE DEFAULT
250 : ! Composite constraint: spell out the fixed components
251 143 : label = ""
252 143 : IF (IAND(mask, comp_ax) /= 0) label = TRIM(label)//" AX"
253 143 : IF (IAND(mask, comp_bx) /= 0) label = TRIM(label)//" BX"
254 143 : IF (IAND(mask, comp_by) /= 0) label = TRIM(label)//" BY"
255 143 : IF (IAND(mask, comp_cx) /= 0) label = TRIM(label)//" CX"
256 143 : IF (IAND(mask, comp_cy) /= 0) label = TRIM(label)//" CY"
257 143 : IF (IAND(mask, comp_cz) /= 0) label = TRIM(label)//" CZ"
258 646 : label = ADJUSTL(label)
259 : END SELECT
260 :
261 503 : END FUNCTION cell_opt_constraint_label
262 :
263 : ! **************************************************************************************************
264 : !> \brief Whether the mask fixes some but not all components of a coupled slot set
265 : !> \param mask bit mask of fixed cell components (comp_* / mask_* bits)
266 : !> \param slots the coupled set of components (a mask_* union)
267 : !> \return ...
268 : !> \note How it compares: IAND(mask, slots) is a set intersection. It holds
269 : !> the components that the mask fixes inside the coupled set. The
270 : !> predicate is true when the intersection is non-empty and is not the
271 : !> full set. In words: the constraint fixes at least one component of
272 : !> the set, and it leaves at least one component free.
273 : !> Why the guards need it: a complete freeze and no freeze both keep
274 : !> the KEEP_ANGLES and KEEP_SYMMETRY projections consistent. A partial
275 : !> fixation breaks them, so the guards reject exactly this case.
276 : ! **************************************************************************************************
277 110 : PURE FUNCTION fixes_partially(mask, slots) RESULT(partial)
278 : INTEGER, INTENT(IN) :: mask, slots
279 : LOGICAL :: partial
280 :
281 110 : partial = IAND(mask, slots) /= 0 .AND. IAND(mask, slots) /= slots
282 110 : END FUNCTION fixes_partially
283 :
284 : ! **************************************************************************************************
285 : !> \brief Whether the constraint mask fixes the given component slot
286 : !> \param mask bit mask of fixed cell components
287 : !> \param slot component index in the canonical order AX..CZ (1..6)
288 : !> \return ...
289 : !> \note How it compares: BTEST(mask, slot - 1) reads one bit of the mask.
290 : !> Each bit position corresponds to one component slot, in the order
291 : !> comp_ax ... comp_cz = 2**0 ... 2**5. The slot argument counts from
292 : !> one so that it indexes the gradient components directly:
293 : !> gradient(1) = AX ... gradient(6) = CZ. BTEST counts bits from zero.
294 : ! **************************************************************************************************
295 36030 : PURE FUNCTION slot_is_fixed(mask, slot) RESULT(fixed)
296 : INTEGER, INTENT(IN) :: mask, slot
297 : LOGICAL :: fixed
298 :
299 36030 : fixed = BTEST(mask, slot - 1)
300 36030 : END FUNCTION slot_is_fixed
301 :
302 : ! **************************************************************************************************
303 : !> \brief Whether a CONSTRAINT mask conflicts with the KEEP_ANGLES projection
304 : !> \param mask bit mask of fixed cell components
305 : !> \param orthorhombic whether the cell is orthorhombic
306 : !> \return conflict flag; the caller aborts when it is true
307 : !> \note KEEP_ANGLES projects the gradient slots {BX,BY} onto b and {CX,CY,CZ}
308 : !> onto c. A conflict arises only on a non-orthorhombic cell, and only
309 : !> when the mask fixes part of the b vector or part of the c vector.
310 : ! **************************************************************************************************
311 74 : PURE FUNCTION keep_angles_conflict(mask, orthorhombic) RESULT(conflict)
312 : INTEGER, INTENT(IN) :: mask
313 : LOGICAL, INTENT(IN) :: orthorhombic
314 : LOGICAL :: conflict
315 :
316 : conflict = .NOT. orthorhombic .AND. &
317 74 : (fixes_partially(mask, mask_b) .OR. fixes_partially(mask, mask_c))
318 74 : END FUNCTION keep_angles_conflict
319 :
320 : ! **************************************************************************************************
321 : !> \brief Whether a CONSTRAINT mask conflicts with the KEEP_SYMMETRY projection
322 : !> \param mask bit mask of fixed cell components
323 : !> \param symmetry_id crystal system of the cell (cell_sym_* constants)
324 : !> \return conflict flag; the caller aborts when it is true
325 : !> \note Each crystal system couples a set of cell components. A conflict
326 : !> arises when the mask fixes part of such a set. Systems without
327 : !> coupling never conflict.
328 : ! **************************************************************************************************
329 74 : PURE FUNCTION keep_symmetry_conflict(mask, symmetry_id) RESULT(conflict)
330 : INTEGER, INTENT(IN) :: mask, symmetry_id
331 : LOGICAL :: conflict
332 :
333 : INTEGER :: coupled_set
334 :
335 74 : coupled_set = keep_symmetry_coupled_set(symmetry_id)
336 74 : conflict = coupled_set /= 0 .AND. fixes_partially(mask, coupled_set)
337 74 : END FUNCTION keep_symmetry_conflict
338 :
339 : ! **************************************************************************************************
340 : !> \brief The set of cell components coupled by a crystal system (for diagnostics)
341 : !> \param symmetry_id crystal system of the cell (cell_sym_* constants)
342 : !> \return mask of the coupled set; zero for systems without coupling
343 : !> \note Mirrors the SELECT CASE in keep_symmetry_conflict. The abort messages
344 : !> print this label so the user can see which components the symmetry
345 : !> couples against their constraint.
346 : ! **************************************************************************************************
347 74 : PURE FUNCTION keep_symmetry_coupled_set(symmetry_id) RESULT(coupled_set)
348 : INTEGER, INTENT(IN) :: symmetry_id
349 : INTEGER :: coupled_set
350 :
351 74 : SELECT CASE (symmetry_id)
352 : CASE (cell_sym_cubic)
353 6 : coupled_set = mask_diag
354 : CASE (cell_sym_rhombohedral)
355 6 : coupled_set = mask_all
356 : CASE (cell_sym_hexagonal_gamma_60, &
357 : cell_sym_hexagonal_gamma_120, &
358 : cell_sym_monoclinic_gamma_ab)
359 18 : coupled_set = mask_ab_plane
360 : CASE (cell_sym_tetragonal_ab)
361 8 : coupled_set = mask_diag_ab
362 : CASE (cell_sym_tetragonal_ac)
363 6 : coupled_set = mask_diag_ac
364 : CASE (cell_sym_tetragonal_bc)
365 6 : coupled_set = mask_diag_bc
366 : CASE DEFAULT
367 74 : coupled_set = 0
368 : END SELECT
369 74 : END FUNCTION keep_symmetry_coupled_set
370 :
371 : ! **************************************************************************************************
372 : !> \brief Whether the masked pressure criterion is structurally off target
373 : !> \param mask bit mask of fixed cell components (comp_* / mask_* bits)
374 : !> \param orthorhombic whether the cell is orthorhombic
375 : !> \param keep_angles whether the KEEP_ANGLES projection is active
376 : !> \return pinned flag; the caller warns, the run stays legal
377 : !> \note On a non-orthorhombic cell the masked pressure criterion (the
378 : !> average of the diagonal stresses over the free diagonal slots)
379 : !> equals the external pressure at the constrained minimum only when
380 : !> the shear stresses conjugate to the pinned shear components vanish.
381 : !> That needs the shear components BX, CX, CY to stay free. The
382 : !> criterion is therefore not guaranteed to converge when the mask
383 : !> pins a shear component, or when KEEP_ANGLES pins the cell angles
384 : !> while a diagonal slot is fixed. Zero or all diagonal slots fixed
385 : !> stay fine (the hydrostatic trace, resp. the external-pressure
386 : !> bypass). Measured on a UO2 cell: the deviation stays at -117 bar
387 : !> for CONSTRAINT Z on a triclinic cell without KEEP_ANGLES and at
388 : !> +110 bar for CONSTRAINT A with KEEP_ANGLES, static under forced
389 : !> continuation, while a symmetric-angle setup converges to -0.003
390 : !> bar. The projected-stress criterion planned in #6031 removes the
391 : !> limitation.
392 : !> The orthorhombic branch is evaluated at setup: a run that starts
393 : !> orthorhombic can grow oblique cell components when a shear slot is
394 : !> pinned while its partners stay free. That residual is second order
395 : !> (two evolved shears must couple into a checked diagonal) and the
396 : !> legacy direction masks X, Y, Z are immune to it because their
397 : !> pinned slots stay zero, where the coupling vanishes.
398 : ! **************************************************************************************************
399 270 : PURE LOGICAL FUNCTION pres_crit_shear_pinned(mask, orthorhombic, keep_angles)
400 : INTEGER, INTENT(IN) :: mask
401 : LOGICAL, INTENT(IN) :: orthorhombic, keep_angles
402 :
403 270 : pres_crit_shear_pinned = .FALSE.
404 270 : IF (.NOT. orthorhombic) THEN
405 94 : IF (IAND(mask, mask_diag) /= mask_diag) THEN
406 90 : pres_crit_shear_pinned = IAND(mask, mask_shear) /= 0
407 90 : IF (.NOT. pres_crit_shear_pinned) THEN
408 82 : pres_crit_shear_pinned = keep_angles .AND. IAND(mask, mask_diag) /= 0
409 : END IF
410 : END IF
411 : END IF
412 270 : END FUNCTION pres_crit_shear_pinned
413 :
414 : ! **************************************************************************************************
415 : !> \brief Classify one CONSTRAINT token without aborting (decision half of the parser)
416 : !> \param token one upper-case token of the CONSTRAINT value line
417 : !> \param accepted .TRUE. when the token contributes a mask
418 : !> \param token_mask the component mask contributed by an accepted token
419 : !> \param reject_reason empty for accepted tokens; the abort message otherwise
420 : !> \note parse_cell_constraint applies this classification and turns a rejection
421 : !> into a cp_abort. Keeping the decision here keeps it unit-testable:
422 : !> the unit test checks the classification table without triggering aborts.
423 : ! **************************************************************************************************
424 458 : PURE SUBROUTINE parse_token(token, accepted, token_mask, reject_reason)
425 : CHARACTER(LEN=*), INTENT(IN) :: token
426 : LOGICAL, INTENT(OUT) :: accepted
427 : INTEGER, INTENT(OUT) :: token_mask
428 : CHARACTER(LEN=default_string_length), INTENT(OUT) :: reject_reason
429 :
430 458 : accepted = .TRUE.
431 458 : token_mask = 0
432 458 : reject_reason = ""
433 82 : SELECT CASE (token)
434 : CASE ("A", "AX")
435 82 : token_mask = mask_a
436 : CASE ("B")
437 20 : token_mask = mask_b
438 : CASE ("BX")
439 56 : token_mask = comp_bx
440 : CASE ("BY")
441 64 : token_mask = comp_by
442 : CASE ("C", "Z")
443 22 : token_mask = mask_c
444 : CASE ("CX")
445 56 : token_mask = comp_cx
446 : CASE ("CY")
447 56 : token_mask = comp_cy
448 : CASE ("CZ")
449 64 : token_mask = comp_cz
450 : CASE ("X")
451 8 : token_mask = mask_x
452 : CASE ("Y")
453 6 : token_mask = mask_y
454 : CASE ("XY")
455 6 : token_mask = mask_xy
456 : CASE ("XZ")
457 4 : token_mask = mask_xz
458 : CASE ("YZ")
459 4 : token_mask = mask_yz
460 : CASE ("ALPHA", "BETA", "GAMMA")
461 2 : accepted = .FALSE.
462 2 : reject_reason = "CELL_OPT CONSTRAINT on angles (ALPHA, BETA, GAMMA) is not supported yet."
463 : CASE ("AY", "AZ", "BZ")
464 6 : accepted = .FALSE.
465 : reject_reason = "CELL_OPT CONSTRAINT token "//TRIM(token)//" is always zero in "// &
466 : "the cell convention (a lies along x, b lies in the xy plane) "// &
467 6 : "and cannot be constrained."
468 : CASE DEFAULT
469 2 : accepted = .FALSE.
470 : reject_reason = "Unknown CELL_OPT CONSTRAINT token: "//TRIM(token)// &
471 : ". Accepted tokens: NONE, the directions X, Y, Z,"// &
472 : " XY, XZ, YZ, the vectors A, B, C, and the components"// &
473 458 : " AX, BX, BY, CX, CY, CZ. Several tokens may be combined."
474 : END SELECT
475 458 : END SUBROUTINE parse_token
476 :
477 : ! **************************************************************************************************
478 : !> \brief Computes the derivatives for the cell
479 : !> \param gradient ...
480 : !> \param av_ptens ...
481 : !> \param pres_ext ...
482 : !> \param cell ...
483 : !> \param mtrx ...
484 : !> \param keep_angles ...
485 : !> \param keep_symmetry ...
486 : !> \param pres_int ...
487 : !> \param pres_constr ...
488 : !> \param constraint_id ...
489 : !> \author Teodoro Laino [tlaino] - University of Zurich - 03.2008
490 : ! **************************************************************************************************
491 5333 : SUBROUTINE get_dg_dh(gradient, av_ptens, pres_ext, cell, mtrx, keep_angles, &
492 : keep_symmetry, pres_int, pres_constr, constraint_id)
493 :
494 : REAL(KIND=dp), DIMENSION(:), POINTER :: gradient
495 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: av_ptens
496 : REAL(KIND=dp), INTENT(IN) :: pres_ext
497 : TYPE(cell_type), POINTER :: cell
498 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: mtrx
499 : LOGICAL, INTENT(IN), OPTIONAL :: keep_angles, keep_symmetry
500 : REAL(KIND=dp), INTENT(OUT) :: pres_int, pres_constr
501 : INTEGER, INTENT(IN), OPTIONAL :: constraint_id
502 :
503 : INTEGER :: i, my_constraint_id
504 : LOGICAL :: my_keep_angles, my_keep_symmetry
505 : LOGICAL, DIMENSION(3) :: diag_free
506 : REAL(KIND=dp), DIMENSION(3) :: diag_ptens
507 : REAL(KIND=dp), DIMENSION(3, 3) :: correction, pten_hinv_old, ptens
508 :
509 5333 : my_keep_angles = .FALSE.
510 5333 : IF (PRESENT(keep_angles)) my_keep_angles = keep_angles
511 5333 : my_keep_symmetry = .FALSE.
512 5333 : IF (PRESENT(keep_symmetry)) my_keep_symmetry = keep_symmetry
513 37331 : gradient = 0.0_dp
514 5333 : IF (PRESENT(constraint_id)) THEN
515 5333 : my_constraint_id = constraint_id
516 : ELSE
517 0 : my_constraint_id = fix_none
518 : END IF
519 :
520 37331 : gradient = 0.0_dp
521 :
522 5333 : ptens = av_ptens
523 :
524 : ! Evaluating the internal pressure
525 5333 : pres_int = 0.0_dp
526 21332 : DO i = 1, 3
527 21332 : pres_int = pres_int + ptens(i, i)
528 : END DO
529 5333 : pres_int = pres_int/3.0_dp
530 :
531 : ! Nominal pressure averaged over the relaxed diagonal cell components
532 : ! (AX, BY, CZ couple to the diagonal of the pressure tensor). If all
533 : ! diagonal components are fixed, no pressure criterion can be relaxed
534 : ! and the external pressure trivially satisfies the tolerance.
535 : diag_free = [IAND(my_constraint_id, comp_ax) == 0, &
536 : IAND(my_constraint_id, comp_by) == 0, &
537 21332 : IAND(my_constraint_id, comp_cz) == 0]
538 21332 : diag_ptens = [ptens(1, 1), ptens(2, 2), ptens(3, 3)]
539 5395 : IF (ANY(diag_free)) THEN
540 42640 : pres_constr = SUM(diag_ptens, MASK=diag_free)/COUNT(diag_free)
541 : ELSE
542 3 : pres_constr = pres_ext
543 : END IF
544 :
545 5333 : ptens(1, 1) = av_ptens(1, 1) - pres_ext
546 5333 : ptens(2, 2) = av_ptens(2, 2) - pres_ext
547 5333 : ptens(3, 3) = av_ptens(3, 3) - pres_ext
548 :
549 277316 : pten_hinv_old = cell%deth*MATMUL(cell%h_inv, ptens)
550 213320 : correction = MATMUL(mtrx, cell%hmat)
551 :
552 5333 : gradient(1) = pten_hinv_old(1, 1) - correction(1, 1)
553 5333 : gradient(2) = pten_hinv_old(2, 1) - correction(2, 1)
554 5333 : gradient(3) = pten_hinv_old(2, 2) - correction(2, 2)
555 5333 : gradient(4) = pten_hinv_old(3, 1) - correction(3, 1)
556 5333 : gradient(5) = pten_hinv_old(3, 2) - correction(3, 2)
557 5333 : gradient(6) = pten_hinv_old(3, 3) - correction(3, 3)
558 :
559 5333 : CALL apply_cell_constraints(gradient, cell, my_keep_angles, my_keep_symmetry, my_constraint_id)
560 :
561 37331 : gradient = -gradient
562 :
563 5333 : END SUBROUTINE get_dg_dh
564 :
565 : ! **************************************************************************************************
566 : !> \brief Apply cell constraints
567 : !> \param gradient ...
568 : !> \param cell ...
569 : !> \param keep_angles ...
570 : !> \param keep_symmetry ...
571 : !> \param constraint_id ...
572 : !> \author Matthias Krack (October 26, 2017, MK)
573 : ! **************************************************************************************************
574 6005 : SUBROUTINE apply_cell_constraints(gradient, cell, keep_angles, keep_symmetry, constraint_id)
575 :
576 : REAL(KIND=dp), DIMENSION(:), POINTER :: gradient
577 : TYPE(cell_type), POINTER :: cell
578 : LOGICAL, INTENT(IN) :: keep_angles, keep_symmetry
579 : INTEGER, INTENT(IN) :: constraint_id
580 :
581 : INTEGER :: slot
582 : REAL(KIND=dp) :: a, a_length, ab_length, b_length, cosa, &
583 : cosah, cosg, deriv_gamma, g, gamma, &
584 : norm, norm_b, norm_c, sina, sinah, sing
585 :
586 6005 : IF (keep_angles) THEN
587 : ! If we want to keep the angles constant we have to project out the
588 : ! components of the cell angles
589 3568 : norm_b = DOT_PRODUCT(cell%hmat(:, 2), cell%hmat(:, 2))
590 2676 : norm = DOT_PRODUCT(cell%hmat(1:2, 2), gradient(2:3))
591 4460 : gradient(2:3) = cell%hmat(1:2, 2)/norm_b*norm
592 3568 : norm_c = DOT_PRODUCT(cell%hmat(:, 3), cell%hmat(:, 3))
593 3568 : norm = DOT_PRODUCT(cell%hmat(1:3, 3), gradient(4:6))
594 6244 : gradient(4:6) = cell%hmat(1:3, 3)/norm_c*norm
595 : ! Retain an exact orthorhombic cell
596 : ! (off-diagonal elements must remain zero identically to keep QS fast)
597 892 : IF (cell%orthorhombic) THEN
598 238 : gradient(2) = 0.0_dp
599 238 : gradient(4) = 0.0_dp
600 238 : gradient(5) = 0.0_dp
601 : END IF
602 : END IF
603 :
604 6005 : IF (keep_symmetry) THEN
605 3184 : SELECT CASE (cell%symmetry_id)
606 : CASE (cell_sym_cubic, &
607 : cell_sym_tetragonal_ab, &
608 : cell_sym_tetragonal_ac, &
609 : cell_sym_tetragonal_bc, &
610 : cell_sym_orthorhombic)
611 104 : SELECT CASE (cell%symmetry_id)
612 : CASE (cell_sym_cubic)
613 104 : g = (gradient(1) + gradient(3) + gradient(6))/3.0_dp
614 104 : gradient(1) = g
615 104 : gradient(3) = g
616 104 : gradient(6) = g
617 : CASE (cell_sym_tetragonal_ab, &
618 : cell_sym_tetragonal_ac, &
619 : cell_sym_tetragonal_bc)
620 588 : SELECT CASE (cell%symmetry_id)
621 : CASE (cell_sym_tetragonal_ab)
622 190 : g = 0.5_dp*(gradient(1) + gradient(3))
623 190 : gradient(1) = g
624 190 : gradient(3) = g
625 : CASE (cell_sym_tetragonal_ac)
626 97 : g = 0.5_dp*(gradient(1) + gradient(6))
627 97 : gradient(1) = g
628 97 : gradient(6) = g
629 : CASE (cell_sym_tetragonal_bc)
630 90 : g = 0.5_dp*(gradient(3) + gradient(6))
631 90 : gradient(3) = g
632 377 : gradient(6) = g
633 : END SELECT
634 : CASE (cell_sym_orthorhombic)
635 : ! Nothing else to do
636 : END SELECT
637 588 : gradient(2) = 0.0_dp
638 588 : gradient(4) = 0.0_dp
639 588 : gradient(5) = 0.0_dp
640 : CASE (cell_sym_hexagonal_gamma_60)
641 240 : g = 0.5_dp*(gradient(1) + 0.5_dp*(gradient(2) + sqrt3*gradient(3)))
642 240 : gradient(1) = g
643 240 : gradient(2) = 0.5_dp*g
644 240 : gradient(3) = sqrt3*gradient(2)
645 240 : gradient(4) = 0.0_dp
646 240 : gradient(5) = 0.0_dp
647 : CASE (cell_sym_hexagonal_gamma_120)
648 122 : g = 0.5_dp*(gradient(1) - 0.5_dp*(gradient(2) - sqrt3*gradient(3)))
649 122 : gradient(1) = g
650 122 : gradient(2) = -0.5_dp*g
651 122 : gradient(3) = -sqrt3*gradient(2)
652 122 : gradient(4) = 0.0_dp
653 122 : gradient(5) = 0.0_dp
654 : CASE (cell_sym_rhombohedral)
655 : a = (angle(cell%hmat(:, 3), cell%hmat(:, 2)) + &
656 : angle(cell%hmat(:, 1), cell%hmat(:, 3)) + &
657 105 : angle(cell%hmat(:, 1), cell%hmat(:, 2)))/3.0_dp
658 105 : cosa = COS(a)
659 105 : sina = SIN(a)
660 105 : cosah = COS(0.5_dp*a)
661 105 : sinah = SIN(0.5_dp*a)
662 105 : norm = cosa/cosah
663 105 : norm_c = SQRT(1.0_dp - norm*norm)
664 : g = (gradient(1) + gradient(2)*cosa + gradient(3)*sina + &
665 105 : gradient(4)*cosah*norm + gradient(5)*sinah*norm + gradient(6)*norm_c)/3.0_dp
666 105 : gradient(1) = g
667 105 : gradient(2) = g*cosa
668 105 : gradient(3) = g*sina
669 105 : gradient(4) = g*cosah*norm
670 105 : gradient(5) = g*sinah*norm
671 105 : gradient(6) = g*norm_c
672 : CASE (cell_sym_monoclinic)
673 1191 : gradient(2) = 0.0_dp
674 1191 : gradient(5) = 0.0_dp
675 : CASE (cell_sym_monoclinic_gamma_ab)
676 : ! Cell symmetry with a = b, alpha = beta = 90 degree and gammma not equal 90 degree
677 : a_length = SQRT(cell%hmat(1, 1)*cell%hmat(1, 1) + &
678 : cell%hmat(2, 1)*cell%hmat(2, 1) + &
679 112 : cell%hmat(3, 1)*cell%hmat(3, 1))
680 : b_length = SQRT(cell%hmat(1, 2)*cell%hmat(1, 2) + &
681 : cell%hmat(2, 2)*cell%hmat(2, 2) + &
682 112 : cell%hmat(3, 2)*cell%hmat(3, 2))
683 112 : ab_length = 0.5_dp*(a_length + b_length)
684 112 : gamma = angle(cell%hmat(:, 1), cell%hmat(:, 2))
685 112 : cosg = COS(gamma)
686 112 : sing = SIN(gamma)
687 : ! Here, g is the average derivative of the cell vector length ab_length,
688 : ! and deriv_gamma is the derivative of the angle gamma
689 112 : g = 0.5_dp*(gradient(1) + cosg*gradient(2) + sing*gradient(3))
690 112 : deriv_gamma = (gradient(3)*cosg - gradient(2)*sing)/b_length
691 112 : gradient(1) = g
692 112 : gradient(2) = g*cosg - ab_length*sing*deriv_gamma
693 112 : gradient(3) = g*sing + ab_length*cosg*deriv_gamma
694 112 : gradient(4) = 0.0_dp
695 2596 : gradient(5) = 0.0_dp
696 : CASE (cell_sym_triclinic)
697 : ! Nothing to do
698 : END SELECT
699 : END IF
700 :
701 : ! Invariant: the projections above run first, the hard constraint
702 : ! zeroing runs last. KEEP_VOLUME is handled separately by
703 : ! rescale_new_cell_volume. No projected gradient survives into a fixed
704 : ! slot. The cell block is always the last six components of the
705 : ! optimization vector.
706 42035 : DO slot = 1, 6
707 42035 : IF (slot_is_fixed(constraint_id, slot)) gradient(slot) = 0.0_dp
708 : END DO
709 :
710 6005 : END SUBROUTINE apply_cell_constraints
711 :
712 : ! **************************************************************************************************
713 : !> \brief Rescale x(idg+1:idg+6) according to vol
714 : !> \param vol ...
715 : !> \param x ...
716 : !> \param idg ...
717 : ! **************************************************************************************************
718 2284 : SUBROUTINE rescale_new_cell_volume(vol, x, idg)
719 :
720 : REAL(KIND=dp), INTENT(IN) :: vol
721 : REAL(KIND=dp), DIMENSION(:), POINTER :: x
722 : INTEGER, INTENT(IN) :: idg
723 :
724 : INTEGER :: i, j, my_idg
725 : REAL(KIND=dp) :: ratio
726 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat_tmp
727 :
728 2284 : CPASSERT((SIZE(x) == idg + 6))
729 2284 : my_idg = idg
730 :
731 2284 : hmat_tmp(:, :) = 0.0_dp
732 9136 : DO i = 1, 3
733 22840 : DO j = 1, i
734 13704 : my_idg = my_idg + 1
735 20556 : hmat_tmp(j, i) = x(my_idg)
736 : END DO
737 : END DO
738 :
739 2284 : ratio = vol/ABS(det_3x3(hmat_tmp))
740 2284 : ratio = ratio**(1.0_dp/3.0_dp) ! no need to use EXP(LOG(ratio)/3) for cube root
741 :
742 2284 : my_idg = idg
743 9136 : DO i = 1, 3
744 22840 : DO j = 1, i
745 13704 : my_idg = my_idg + 1
746 20556 : x(my_idg) = x(my_idg)*ratio
747 : END DO
748 : END DO
749 :
750 2284 : END SUBROUTINE rescale_new_cell_volume
751 :
752 : END MODULE cell_opt_utils
|
