Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 : USE message_passing, ONLY: mp_para_env_type
53 : #include "../base/base_uses.f90"
54 :
55 : IMPLICIT NONE
56 : PRIVATE
57 :
58 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
59 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cell_opt_utils'
60 :
61 : PUBLIC :: get_dg_dh, gopt_new_logger_create, &
62 : gopt_new_logger_release, read_external_press_tensor, &
63 : apply_cell_constraints
64 :
65 : CONTAINS
66 :
67 : ! **************************************************************************************************
68 : !> \brief creates a new logger used for cell optimization algorithm
69 : !> \param new_logger ...
70 : !> \param root_section ...
71 : !> \param para_env ...
72 : !> \param project_name ...
73 : !> \param id_run ...
74 : !> \author Teodoro Laino [tlaino] - University of Zurich - 03.2008
75 : ! **************************************************************************************************
76 930 : SUBROUTINE gopt_new_logger_create(new_logger, root_section, para_env, project_name, &
77 : id_run)
78 : TYPE(cp_logger_type), POINTER :: new_logger
79 : TYPE(section_vals_type), POINTER :: root_section
80 : TYPE(mp_para_env_type), POINTER :: para_env
81 : CHARACTER(len=default_string_length), INTENT(OUT) :: project_name
82 : INTEGER, INTENT(IN) :: id_run
83 :
84 : CHARACTER(len=default_path_length) :: c_val, input_file_path, output_file_path
85 : CHARACTER(len=default_string_length) :: label
86 : INTEGER :: i, lp, unit_nr
87 : TYPE(cp_logger_type), POINTER :: logger
88 : TYPE(enumeration_type), POINTER :: enum
89 : TYPE(keyword_type), POINTER :: keyword
90 : TYPE(section_type), POINTER :: section
91 :
92 930 : NULLIFY (new_logger, logger, enum, keyword, section)
93 930 : logger => cp_get_default_logger()
94 :
95 930 : CALL create_global_section(section)
96 930 : keyword => section_get_keyword(section, "RUN_TYPE")
97 930 : CALL keyword_get(keyword, enum=enum)
98 930 : label = TRIM(enum_i2c(enum, id_run))
99 930 : CALL section_release(section)
100 :
101 : ! Redirecting output of the sub_calculation to a different file
102 930 : CALL section_vals_val_get(root_section, "GLOBAL%PROJECT_NAME", c_val=project_name)
103 930 : input_file_path = TRIM(project_name)
104 930 : lp = LEN_TRIM(input_file_path)
105 930 : i = logger%iter_info%iteration(logger%iter_info%n_rlevel)
106 930 : input_file_path(lp + 1:LEN(input_file_path)) = "-"//TRIM(label)//"-"//ADJUSTL(cp_to_string(i))
107 930 : lp = LEN_TRIM(input_file_path)
108 930 : CALL section_vals_val_set(root_section, "GLOBAL%PROJECT_NAME", c_val=input_file_path(1:lp))
109 930 : CALL section_vals_val_set(root_section, "GLOBAL%RUN_TYPE", i_val=id_run)
110 :
111 : ! Redirecting output into a new file
112 930 : output_file_path = input_file_path(1:lp)//".out"
113 930 : IF (para_env%is_source()) THEN
114 : CALL open_file(file_name=output_file_path, file_status="UNKNOWN", &
115 465 : file_action="WRITE", file_position="APPEND", unit_number=unit_nr)
116 : ELSE
117 465 : unit_nr = -1
118 : END IF
119 : CALL cp_logger_create(new_logger, para_env=para_env, default_global_unit_nr=unit_nr, &
120 930 : close_global_unit_on_dealloc=.FALSE.)
121 930 : CALL section_vals_val_get(root_section, "GLOBAL%PROJECT", c_val=c_val)
122 930 : IF (c_val /= "") THEN
123 930 : CALL cp_logger_set(new_logger, local_filename=TRIM(c_val)//"_localLog")
124 : END IF
125 930 : new_logger%iter_info%project_name = TRIM(c_val)
126 : CALL section_vals_val_get(root_section, "GLOBAL%PRINT_LEVEL", &
127 930 : i_val=new_logger%iter_info%print_level)
128 :
129 930 : END SUBROUTINE gopt_new_logger_create
130 :
131 : ! **************************************************************************************************
132 : !> \brief releases a new logger used for cell optimization algorithm
133 : !> \param new_logger ...
134 : !> \param root_section ...
135 : !> \param para_env ...
136 : !> \param project_name ...
137 : !> \param id_run ...
138 : !> \author Teodoro Laino [tlaino] - University of Zurich - 03.2008
139 : ! **************************************************************************************************
140 930 : SUBROUTINE gopt_new_logger_release(new_logger, root_section, para_env, project_name, id_run)
141 : TYPE(cp_logger_type), POINTER :: new_logger
142 : TYPE(section_vals_type), POINTER :: root_section
143 : TYPE(mp_para_env_type), POINTER :: para_env
144 : CHARACTER(len=default_string_length), INTENT(IN) :: project_name
145 : INTEGER, INTENT(IN) :: id_run
146 :
147 : INTEGER :: unit_nr
148 :
149 930 : IF (para_env%is_source()) THEN
150 465 : unit_nr = cp_logger_get_default_unit_nr(new_logger)
151 465 : CALL close_file(unit_number=unit_nr)
152 : END IF
153 930 : CALL cp_logger_release(new_logger)
154 930 : CALL section_vals_val_set(root_section, "GLOBAL%RUN_TYPE", i_val=id_run)
155 930 : CALL section_vals_val_set(root_section, "GLOBAL%PROJECT_NAME", c_val=project_name)
156 :
157 930 : END SUBROUTINE gopt_new_logger_release
158 :
159 : ! **************************************************************************************************
160 : !> \brief Reads the external pressure tensor
161 : !> \param geo_section ...
162 : !> \param cell ...
163 : !> \param pres_ext ...
164 : !> \param mtrx ...
165 : !> \param rot ...
166 : !> \author Teodoro Laino [tlaino] - University of Zurich - 03.2008
167 : ! **************************************************************************************************
168 218 : SUBROUTINE read_external_press_tensor(geo_section, cell, pres_ext, mtrx, rot)
169 : TYPE(section_vals_type), POINTER :: geo_section
170 : TYPE(cell_type), POINTER :: cell
171 : REAL(KIND=dp), INTENT(OUT) :: pres_ext
172 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT) :: mtrx
173 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: rot
174 :
175 : INTEGER :: i, ind, j
176 : LOGICAL :: check
177 : REAL(KIND=dp), DIMENSION(3, 3) :: pres_ext_tens
178 218 : REAL(KIND=dp), DIMENSION(:), POINTER :: pvals
179 :
180 218 : NULLIFY (pvals)
181 218 : mtrx = 0.0_dp
182 218 : pres_ext_tens = 0.0_dp
183 218 : pres_ext = 0.0_dp
184 218 : CALL section_vals_val_get(geo_section, "EXTERNAL_PRESSURE", r_vals=pvals)
185 218 : check = (SIZE(pvals) == 1) .OR. (SIZE(pvals) == 9)
186 218 : IF (.NOT. check) &
187 0 : CPABORT("EXTERNAL_PRESSURE can have 1 or 9 components only!")
188 :
189 218 : IF (SIZE(pvals) == 9) THEN
190 : ind = 0
191 424 : DO i = 1, 3
192 1378 : DO j = 1, 3
193 954 : ind = ind + 1
194 1272 : pres_ext_tens(j, i) = pvals(ind)
195 : END DO
196 : END DO
197 : ! Also the pressure tensor must be oriented in the same canonical directions
198 : ! of the simulation cell
199 4240 : pres_ext_tens = MATMUL(TRANSPOSE(rot), pres_ext_tens)
200 424 : DO i = 1, 3
201 424 : pres_ext = pres_ext + pres_ext_tens(i, i)
202 : END DO
203 106 : pres_ext = pres_ext/3.0_dp
204 424 : DO i = 1, 3
205 424 : pres_ext_tens(i, i) = pres_ext_tens(i, i) - pres_ext
206 : END DO
207 : ELSE
208 112 : pres_ext = pvals(1)
209 : END IF
210 :
211 2834 : IF (ANY(pres_ext_tens > 1.0E-5_dp)) THEN
212 0 : mtrx = cell%deth*MATMUL(cell%h_inv, MATMUL(pres_ext_tens, TRANSPOSE(cell%h_inv)))
213 : END IF
214 :
215 218 : END SUBROUTINE read_external_press_tensor
216 :
217 : ! **************************************************************************************************
218 : !> \brief Computes the derivatives for the cell
219 : !> \param gradient ...
220 : !> \param av_ptens ...
221 : !> \param pres_ext ...
222 : !> \param cell ...
223 : !> \param mtrx ...
224 : !> \param keep_angles ...
225 : !> \param keep_symmetry ...
226 : !> \param pres_int ...
227 : !> \param pres_constr ...
228 : !> \param constraint_id ...
229 : !> \author Teodoro Laino [tlaino] - University of Zurich - 03.2008
230 : ! **************************************************************************************************
231 4070 : SUBROUTINE get_dg_dh(gradient, av_ptens, pres_ext, cell, mtrx, keep_angles, &
232 : keep_symmetry, pres_int, pres_constr, constraint_id)
233 :
234 : REAL(KIND=dp), DIMENSION(:), POINTER :: gradient
235 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: av_ptens
236 : REAL(KIND=dp), INTENT(IN) :: pres_ext
237 : TYPE(cell_type), POINTER :: cell
238 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: mtrx
239 : LOGICAL, INTENT(IN), OPTIONAL :: keep_angles, keep_symmetry
240 : REAL(KIND=dp), INTENT(OUT) :: pres_int, pres_constr
241 : INTEGER, INTENT(IN), OPTIONAL :: constraint_id
242 :
243 : INTEGER :: i, my_constraint_id
244 : LOGICAL :: my_keep_angles, my_keep_symmetry
245 : REAL(KIND=dp), DIMENSION(3, 3) :: correction, pten_hinv_old, ptens
246 :
247 4070 : my_keep_angles = .FALSE.
248 4070 : IF (PRESENT(keep_angles)) my_keep_angles = keep_angles
249 4070 : my_keep_symmetry = .FALSE.
250 4070 : IF (PRESENT(keep_symmetry)) my_keep_symmetry = keep_symmetry
251 28490 : gradient = 0.0_dp
252 4070 : IF (PRESENT(constraint_id)) THEN
253 4070 : my_constraint_id = constraint_id
254 : ELSE
255 0 : my_constraint_id = fix_none
256 : END IF
257 :
258 28490 : gradient = 0.0_dp
259 :
260 4070 : ptens = av_ptens
261 :
262 : ! Evaluating the internal pressure
263 4070 : pres_int = 0.0_dp
264 16280 : DO i = 1, 3
265 16280 : pres_int = pres_int + ptens(i, i)
266 : END DO
267 4070 : pres_int = pres_int/3.0_dp
268 :
269 0 : SELECT CASE (my_constraint_id)
270 : CASE (fix_x)
271 0 : pres_constr = ptens(2, 2) + ptens(3, 3)
272 : CASE (fix_y)
273 0 : pres_constr = ptens(1, 1) + ptens(3, 3)
274 : CASE (fix_z)
275 8 : pres_constr = ptens(1, 1) + ptens(2, 2)
276 : CASE (fix_xy)
277 6 : pres_constr = ptens(3, 3)
278 : CASE (fix_xz)
279 0 : pres_constr = ptens(2, 2)
280 : CASE (fix_yz)
281 0 : pres_constr = ptens(1, 1)
282 : CASE (fix_none)
283 4070 : pres_constr = ptens(1, 1) + ptens(2, 2) + ptens(3, 3)
284 : END SELECT
285 4070 : pres_constr = pres_constr/3.0_dp
286 :
287 4070 : ptens(1, 1) = av_ptens(1, 1) - pres_ext
288 4070 : ptens(2, 2) = av_ptens(2, 2) - pres_ext
289 4070 : ptens(3, 3) = av_ptens(3, 3) - pres_ext
290 :
291 211640 : pten_hinv_old = cell%deth*MATMUL(cell%h_inv, ptens)
292 162800 : correction = MATMUL(mtrx, cell%hmat)
293 :
294 4070 : gradient(1) = pten_hinv_old(1, 1) - correction(1, 1)
295 4070 : gradient(2) = pten_hinv_old(2, 1) - correction(2, 1)
296 4070 : gradient(3) = pten_hinv_old(2, 2) - correction(2, 2)
297 4070 : gradient(4) = pten_hinv_old(3, 1) - correction(3, 1)
298 4070 : gradient(5) = pten_hinv_old(3, 2) - correction(3, 2)
299 4070 : gradient(6) = pten_hinv_old(3, 3) - correction(3, 3)
300 :
301 4070 : CALL apply_cell_constraints(gradient, cell, my_keep_angles, my_keep_symmetry, my_constraint_id)
302 :
303 28490 : gradient = -gradient
304 :
305 4070 : 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 4070 : 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 4070 : 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 624 : norm_b = DOT_PRODUCT(cell%hmat(:, 2), cell%hmat(:, 2))
331 468 : norm = DOT_PRODUCT(cell%hmat(1:2, 2), gradient(2:3))
332 780 : gradient(2:3) = cell%hmat(1:2, 2)/norm_b*norm
333 624 : norm_c = DOT_PRODUCT(cell%hmat(:, 3), cell%hmat(:, 3))
334 624 : norm = DOT_PRODUCT(cell%hmat(1:3, 3), gradient(4:6))
335 1092 : 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 156 : IF (cell%orthorhombic) THEN
339 98 : gradient(2) = 0.0_dp
340 98 : gradient(4) = 0.0_dp
341 98 : gradient(5) = 0.0_dp
342 : END IF
343 : END IF
344 :
345 4070 : IF (keep_symmetry) THEN
346 2028 : 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 125 : gradient(2) = 0.0_dp
415 125 : 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 1464 : gradient(5) = 0.0_dp
436 : CASE (cell_sym_triclinic)
437 : ! Nothing to do
438 : END SELECT
439 : END IF
440 :
441 4070 : 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 32 : 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 4070 : gradient(2:6) = 0.0_dp
457 : CASE (fix_none)
458 : ! Nothing to do
459 : END SELECT
460 :
461 4070 : END SUBROUTINE apply_cell_constraints
462 :
463 : END MODULE cell_opt_utils
|
