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 : !> none
13 : ! **************************************************************************************************
14 : MODULE gopt_f_methods
15 :
16 : USE atomic_kind_list_types, ONLY: atomic_kind_list_type
17 : USE atomic_kind_types, ONLY: atomic_kind_type, &
18 : get_atomic_kind_set
19 : USE cell_methods, ONLY: cell_create, &
20 : init_cell
21 : USE cell_types, ONLY: cell_copy, &
22 : cell_release, &
23 : cell_type, &
24 : real_to_scaled, &
25 : scaled_to_real
26 : USE cp_log_handling, ONLY: cp_to_string
27 : USE cp_subsys_types, ONLY: cp_subsys_get, &
28 : cp_subsys_set, &
29 : cp_subsys_type, &
30 : pack_subsys_particles
31 : USE cp_units, ONLY: cp_unit_from_cp2k
32 : USE dimer_types, ONLY: dimer_env_type
33 : USE dimer_utils, ONLY: update_dimer_vec
34 : USE distribution_1d_types, ONLY: distribution_1d_type
35 : USE force_env_types, ONLY: force_env_get, &
36 : force_env_get_natom, &
37 : force_env_get_nparticle, &
38 : force_env_type, &
39 : use_qmmm, &
40 : use_qmmmx
41 : USE gopt_f_types, ONLY: gopt_f_type
42 : USE gopt_param_types, ONLY: gopt_param_type
43 : USE input_constants, ONLY: default_cell_direct_id, &
44 : default_cell_geo_opt_id, &
45 : default_cell_md_id, &
46 : default_cell_method_id, &
47 : default_minimization_method_id, &
48 : default_shellcore_method_id, &
49 : default_ts_method_id, &
50 : fix_none, &
51 : fix_x, &
52 : fix_xy, &
53 : fix_xz, &
54 : fix_y, &
55 : fix_yz, &
56 : fix_z
57 : USE input_cp2k_restarts, ONLY: write_restart
58 : USE input_section_types, ONLY: section_vals_type, &
59 : section_vals_val_get
60 : USE kinds, ONLY: default_string_length, &
61 : dp, &
62 : int_8
63 : USE machine, ONLY: m_flush
64 : USE md_energies, ONLY: sample_memory
65 : USE message_passing, ONLY: mp_para_env_type
66 : USE motion_utils, ONLY: write_simulation_cell, &
67 : write_stress_tensor_to_file, &
68 : write_trajectory
69 : USE particle_list_types, ONLY: particle_list_type
70 : USE particle_methods, ONLY: write_final_cif, &
71 : write_structure_data
72 : USE particle_types, ONLY: particle_type
73 : USE qmmm_util, ONLY: apply_qmmm_translate
74 : USE qmmmx_util, ONLY: apply_qmmmx_translate
75 : USE virial_methods, ONLY: virial_evaluate
76 : USE virial_types, ONLY: virial_type
77 : #include "../base/base_uses.f90"
78 :
79 : IMPLICIT NONE
80 : PRIVATE
81 :
82 : #:include "gopt_f77_methods.fypp"
83 :
84 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
85 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = "gopt_f_methods"
86 :
87 : PUBLIC :: gopt_f_create_x0, &
88 : print_geo_opt_header, print_geo_opt_nc, &
89 : gopt_f_io_init, gopt_f_io, gopt_f_io_finalize, gopt_f_ii, &
90 : apply_cell_change
91 :
92 : CONTAINS
93 :
94 : ! **************************************************************************************************
95 : !> \brief returns the value of the parameters for the actual configuration
96 : !> \param gopt_env the geometry optimization environment you want the info about
97 : !> x0: the parameter vector (is allocated by this routine)
98 : !> \param x0 ...
99 : !> \par History
100 : !> - Cell optimization revised (06.11.2012,MK)
101 : ! **************************************************************************************************
102 1995 : SUBROUTINE gopt_f_create_x0(gopt_env, x0)
103 :
104 : TYPE(gopt_f_type), POINTER :: gopt_env
105 : REAL(KIND=dp), DIMENSION(:), POINTER :: x0
106 :
107 : INTEGER :: i, idg, j, nparticle
108 : TYPE(cell_type), POINTER :: cell
109 : TYPE(cp_subsys_type), POINTER :: subsys
110 :
111 1995 : NULLIFY (cell)
112 1995 : NULLIFY (subsys)
113 :
114 3766 : SELECT CASE (gopt_env%type_id)
115 : CASE (default_minimization_method_id, default_ts_method_id)
116 1771 : CALL force_env_get(gopt_env%force_env, subsys=subsys)
117 : ! before starting we handle the case of translating coordinates (QM/MM)
118 1771 : IF (gopt_env%force_env%in_use == use_qmmm) &
119 36 : CALL apply_qmmm_translate(gopt_env%force_env%qmmm_env)
120 1771 : IF (gopt_env%force_env%in_use == use_qmmmx) &
121 0 : CALL apply_qmmmx_translate(gopt_env%force_env%qmmmx_env)
122 1771 : nparticle = force_env_get_nparticle(gopt_env%force_env)
123 5313 : ALLOCATE (x0(3*nparticle))
124 1771 : CALL pack_subsys_particles(subsys=subsys, r=x0)
125 : CASE (default_cell_method_id)
126 406 : SELECT CASE (gopt_env%cell_method_id)
127 : CASE (default_cell_direct_id)
128 182 : CALL force_env_get(gopt_env%force_env, subsys=subsys, cell=cell)
129 : ! Store reference cell
130 4732 : gopt_env%h_ref = cell%hmat
131 : ! before starting we handle the case of translating coordinates (QM/MM)
132 182 : IF (gopt_env%force_env%in_use == use_qmmm) &
133 0 : CALL apply_qmmm_translate(gopt_env%force_env%qmmm_env)
134 182 : IF (gopt_env%force_env%in_use == use_qmmmx) &
135 0 : CALL apply_qmmmx_translate(gopt_env%force_env%qmmmx_env)
136 182 : nparticle = force_env_get_nparticle(gopt_env%force_env)
137 546 : ALLOCATE (x0(3*nparticle + 6))
138 182 : CALL pack_subsys_particles(subsys=subsys, r=x0)
139 182 : idg = 3*nparticle
140 728 : DO i = 1, 3
141 1820 : DO j = 1, i
142 1092 : idg = idg + 1
143 1638 : x0(idg) = cell%hmat(j, i)
144 : END DO
145 : END DO
146 : CASE (default_cell_geo_opt_id, default_cell_md_id)
147 42 : CALL force_env_get(gopt_env%force_env, cell=cell)
148 42 : ALLOCATE (x0(6))
149 42 : idg = 0
150 168 : DO i = 1, 3
151 420 : DO j = 1, i
152 252 : idg = idg + 1
153 378 : x0(idg) = cell%hmat(j, i)
154 : END DO
155 : END DO
156 : CASE DEFAULT
157 224 : CPABORT("Invalid or not yet implemented type of cell optimization")
158 : END SELECT
159 : CASE DEFAULT
160 1995 : CPABORT("Invalid or not yet implemented type of optimization")
161 : END SELECT
162 :
163 1995 : END SUBROUTINE gopt_f_create_x0
164 :
165 : ! **************************************************************************************************
166 : !> \brief Prints iteration step of the optimization procedure on screen
167 : !> \param its ...
168 : !> \param output_unit ...
169 : !> \author Teodoro Laino [tlaino] - University of Zurich - 03.2008
170 : ! **************************************************************************************************
171 13448 : SUBROUTINE gopt_f_ii(its, output_unit)
172 :
173 : INTEGER, INTENT(IN) :: its, output_unit
174 :
175 13448 : IF (output_unit > 0) THEN
176 6762 : WRITE (UNIT=output_unit, FMT="(/,T2,26('-'))")
177 6762 : WRITE (UNIT=output_unit, FMT="(T2,A,I6)") "OPTIMIZATION STEP: ", its
178 6762 : WRITE (UNIT=output_unit, FMT="(T2,26('-'))")
179 6762 : CALL m_flush(output_unit)
180 : END IF
181 :
182 13448 : END SUBROUTINE gopt_f_ii
183 :
184 : ! **************************************************************************************************
185 : !> \brief Handles the Output during an optimization run
186 : !> \param gopt_env ...
187 : !> \param output_unit ...
188 : !> \param opt_energy ...
189 : !> \param wildcard ...
190 : !> \param its ...
191 : !> \param used_time ...
192 : !> \author Teodoro Laino [tlaino] - University of Zurich - 03.2008
193 : ! **************************************************************************************************
194 1729 : SUBROUTINE gopt_f_io_init(gopt_env, output_unit, opt_energy, wildcard, its, used_time)
195 :
196 : TYPE(gopt_f_type), POINTER :: gopt_env
197 : INTEGER, INTENT(IN) :: output_unit
198 : REAL(KIND=dp) :: opt_energy
199 : CHARACTER(LEN=5) :: wildcard
200 : INTEGER, INTENT(IN) :: its
201 : REAL(KIND=dp) :: used_time
202 :
203 : TYPE(mp_para_env_type), POINTER :: para_env
204 : CHARACTER(LEN=default_string_length) :: energy_unit, stress_unit
205 : REAL(KIND=dp) :: pres_int
206 : INTEGER(KIND=int_8) :: max_memory
207 : LOGICAL :: print_memory
208 :
209 1729 : NULLIFY (para_env)
210 1729 : CALL section_vals_val_get(gopt_env%motion_section, "PRINT%MEMORY_INFO", l_val=print_memory)
211 1729 : max_memory = 0
212 1729 : IF (print_memory) THEN
213 1729 : CALL force_env_get(gopt_env%force_env, para_env=para_env)
214 1729 : max_memory = sample_memory(para_env)
215 : END IF
216 :
217 : CALL section_vals_val_get(gopt_env%force_env%force_env_section, &
218 : "PRINT%PROGRAM_RUN_INFO%ENERGY_UNIT", &
219 1729 : c_val=energy_unit)
220 : CALL section_vals_val_get(gopt_env%force_env%force_env_section, &
221 : "PRINT%STRESS_TENSOR%STRESS_UNIT", &
222 1729 : c_val=stress_unit)
223 :
224 3262 : SELECT CASE (gopt_env%type_id)
225 : CASE (default_ts_method_id, default_minimization_method_id)
226 : ! Geometry Optimization (Minimization and Transition State Search)
227 1533 : IF (.NOT. gopt_env%dimer_rotation) THEN
228 : CALL write_cycle_infos(output_unit, &
229 : it=its, &
230 : etot=opt_energy, &
231 : wildcard=wildcard, &
232 : used_time=used_time, &
233 : max_memory=max_memory, &
234 : energy_unit=energy_unit, &
235 1401 : stress_unit=stress_unit)
236 : ELSE
237 : CALL write_rot_cycle_infos(output_unit, &
238 : it=its, &
239 : etot=opt_energy, &
240 : dimer_env=gopt_env%dimer_env, &
241 : wildcard=wildcard, &
242 : used_time=used_time, &
243 132 : max_memory=max_memory)
244 : END IF
245 : CASE (default_cell_method_id)
246 : ! Cell Optimization
247 176 : pres_int = gopt_env%cell_env%pres_int
248 : CALL write_cycle_infos(output_unit, &
249 : it=its, &
250 : etot=opt_energy, &
251 : pres_int=pres_int, &
252 : wildcard=wildcard, &
253 : used_time=used_time, &
254 : max_memory=max_memory, &
255 : energy_unit=energy_unit, &
256 176 : stress_unit=stress_unit)
257 : CASE (default_shellcore_method_id)
258 : CALL write_cycle_infos(output_unit, &
259 : it=its, &
260 : etot=opt_energy, &
261 : wildcard=wildcard, &
262 : used_time=used_time, &
263 : max_memory=max_memory, &
264 : energy_unit=energy_unit, &
265 1729 : stress_unit=stress_unit)
266 : END SELECT
267 :
268 1729 : END SUBROUTINE gopt_f_io_init
269 :
270 : ! **************************************************************************************************
271 : !> \brief Handles the Output during an optimization run
272 : !> \param gopt_env ...
273 : !> \param force_env ...
274 : !> \param root_section ...
275 : !> \param its ...
276 : !> \param opt_energy ...
277 : !> \param output_unit ...
278 : !> \param eold ...
279 : !> \param emin ...
280 : !> \param wildcard ...
281 : !> \param gopt_param ...
282 : !> \param ndf ...
283 : !> \param dx ...
284 : !> \param xi ...
285 : !> \param conv ...
286 : !> \param pred ...
287 : !> \param rat ...
288 : !> \param step ...
289 : !> \param rad ...
290 : !> \param used_time ...
291 : !> \author Teodoro Laino [tlaino] - University of Zurich - 03.2008
292 : ! **************************************************************************************************
293 26896 : SUBROUTINE gopt_f_io(gopt_env, force_env, root_section, its, opt_energy, &
294 13448 : output_unit, eold, emin, wildcard, gopt_param, ndf, dx, xi, conv, pred, rat, &
295 : step, rad, used_time)
296 :
297 : TYPE(gopt_f_type), POINTER :: gopt_env
298 : TYPE(force_env_type), POINTER :: force_env
299 : TYPE(section_vals_type), POINTER :: root_section
300 : INTEGER, INTENT(IN) :: its
301 : REAL(KIND=dp), INTENT(IN) :: opt_energy
302 : INTEGER, INTENT(IN) :: output_unit
303 : REAL(KIND=dp) :: eold, emin
304 : CHARACTER(LEN=5) :: wildcard
305 : TYPE(gopt_param_type), POINTER :: gopt_param
306 : INTEGER, INTENT(IN), OPTIONAL :: ndf
307 : REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: dx
308 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: xi
309 : LOGICAL, OPTIONAL :: conv
310 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: pred, rat, step, rad
311 : REAL(KIND=dp) :: used_time
312 :
313 : CHARACTER(LEN=default_string_length) :: energy_unit, stress_unit
314 : INTEGER(KIND=int_8) :: max_memory
315 : LOGICAL :: print_memory
316 : REAL(KIND=dp) :: pres_diff, pres_diff_constr, pres_int, &
317 : pres_tol
318 : TYPE(mp_para_env_type), POINTER :: para_env
319 :
320 13448 : NULLIFY (para_env)
321 13448 : CALL section_vals_val_get(gopt_env%motion_section, "PRINT%MEMORY_INFO", l_val=print_memory)
322 13448 : max_memory = 0
323 13448 : IF (print_memory) THEN
324 13448 : CALL force_env_get(force_env, para_env=para_env)
325 13448 : max_memory = sample_memory(para_env)
326 : END IF
327 :
328 : CALL section_vals_val_get(gopt_env%force_env%force_env_section, &
329 : "PRINT%PROGRAM_RUN_INFO%ENERGY_UNIT", &
330 13448 : c_val=energy_unit)
331 : CALL section_vals_val_get(gopt_env%force_env%force_env_section, &
332 : "PRINT%STRESS_TENSOR%STRESS_UNIT", &
333 13448 : c_val=stress_unit)
334 :
335 22382 : SELECT CASE (gopt_env%type_id)
336 : CASE (default_ts_method_id, default_minimization_method_id)
337 : ! Geometry Optimization (Minimization and Transition State Search)
338 8934 : IF (.NOT. gopt_env%dimer_rotation) THEN
339 : CALL geo_opt_io(force_env=force_env, root_section=root_section, &
340 8208 : motion_section=gopt_env%motion_section, its=its, opt_energy=opt_energy)
341 : CALL write_cycle_infos(output_unit, &
342 : it=its, &
343 : etot=opt_energy, &
344 : ediff=(opt_energy - eold), &
345 : pred=pred, &
346 : rat=rat, &
347 : step=step, &
348 : rad=rad, &
349 : emin=emin, &
350 : wildcard=wildcard, &
351 : used_time=used_time, &
352 : max_memory=max_memory, &
353 : energy_unit=energy_unit, &
354 8208 : stress_unit=stress_unit)
355 : ! Possibly check convergence
356 8208 : IF (PRESENT(conv)) THEN
357 8208 : CPASSERT(PRESENT(ndf))
358 8208 : CPASSERT(PRESENT(dx))
359 8208 : CPASSERT(PRESENT(xi))
360 8208 : CALL check_converg(ndf, dx, xi, output_unit, conv, gopt_param, max_memory, stress_unit)
361 : END IF
362 : ELSE
363 726 : CALL update_dimer_vec(gopt_env%dimer_env, gopt_env%motion_section)
364 726 : CALL write_restart(force_env=force_env, root_section=root_section)
365 : CALL write_rot_cycle_infos(output_unit, its, opt_energy, opt_energy - eold, emin, gopt_env%dimer_env, &
366 726 : wildcard=wildcard, used_time=used_time, max_memory=max_memory)
367 : ! Possibly check convergence
368 726 : IF (PRESENT(conv)) THEN
369 726 : CPASSERT(ASSOCIATED(gopt_env%dimer_env))
370 726 : CALL check_rot_conv(gopt_env%dimer_env, output_unit, conv)
371 : END IF
372 : END IF
373 : CASE (default_cell_method_id)
374 : ! Cell Optimization
375 4344 : pres_diff = gopt_env%cell_env%pres_int - gopt_env%cell_env%pres_ext
376 4344 : pres_int = gopt_env%cell_env%pres_int
377 4344 : pres_tol = gopt_env%cell_env%pres_tol
378 : CALL geo_opt_io(force_env=force_env, root_section=root_section, &
379 4344 : motion_section=gopt_env%motion_section, its=its, opt_energy=opt_energy)
380 : CALL write_cycle_infos(output_unit, &
381 : it=its, &
382 : etot=opt_energy, &
383 : ediff=(opt_energy - eold), &
384 : pred=pred, &
385 : rat=rat, &
386 : step=step, &
387 : rad=rad, &
388 : emin=emin, &
389 : pres_int=pres_int, &
390 : wildcard=wildcard, &
391 : used_time=used_time, &
392 : max_memory=max_memory, &
393 : energy_unit=energy_unit, &
394 4344 : stress_unit=stress_unit)
395 : ! Possibly check convergence
396 4344 : IF (PRESENT(conv)) THEN
397 4344 : CPASSERT(PRESENT(ndf))
398 4344 : CPASSERT(PRESENT(dx))
399 4344 : CPASSERT(PRESENT(xi))
400 4344 : IF (gopt_env%cell_env%constraint_id == fix_none) THEN
401 : CALL check_converg(ndf, dx, xi, output_unit, conv, gopt_param, max_memory, stress_unit, &
402 4326 : pres_diff, pres_tol)
403 : ELSE
404 18 : pres_diff_constr = gopt_env%cell_env%pres_constr - gopt_env%cell_env%pres_ext
405 : CALL check_converg(ndf, dx, xi, output_unit, conv, gopt_param, max_memory, stress_unit, &
406 18 : pres_diff, pres_tol, pres_diff_constr)
407 : END IF
408 : END IF
409 : CASE (default_shellcore_method_id)
410 : CALL write_cycle_infos(output_unit, &
411 : it=its, &
412 : etot=opt_energy, &
413 : ediff=(opt_energy - eold), &
414 : pred=pred, &
415 : rat=rat, &
416 : step=step, &
417 : rad=rad, &
418 : emin=emin, &
419 : wildcard=wildcard, &
420 : used_time=used_time, &
421 : max_memory=max_memory, &
422 : energy_unit=energy_unit, &
423 170 : stress_unit=stress_unit)
424 : ! Possibly check convergence
425 13618 : IF (PRESENT(conv)) THEN
426 170 : CPASSERT(PRESENT(ndf))
427 170 : CPASSERT(PRESENT(dx))
428 170 : CPASSERT(PRESENT(xi))
429 170 : CALL check_converg(ndf, dx, xi, output_unit, conv, gopt_param, max_memory, stress_unit)
430 : END IF
431 : END SELECT
432 :
433 13448 : END SUBROUTINE gopt_f_io
434 :
435 : ! **************************************************************************************************
436 : !> \brief Handles the Output at the end of an optimization run
437 : !> \param gopt_env ...
438 : !> \param force_env ...
439 : !> \param x0 ...
440 : !> \param conv ...
441 : !> \param its ...
442 : !> \param root_section ...
443 : !> \param para_env ...
444 : !> \param master ...
445 : !> \param output_unit ...
446 : !> \author Teodoro Laino [tlaino] - University of Zurich - 03.2008
447 : ! **************************************************************************************************
448 2147 : RECURSIVE SUBROUTINE gopt_f_io_finalize(gopt_env, force_env, x0, conv, its, root_section, &
449 : para_env, master, output_unit)
450 : TYPE(gopt_f_type), POINTER :: gopt_env
451 : TYPE(force_env_type), POINTER :: force_env
452 : REAL(KIND=dp), DIMENSION(:), POINTER :: x0
453 : LOGICAL :: conv
454 : INTEGER :: its
455 : TYPE(section_vals_type), POINTER :: root_section
456 : TYPE(mp_para_env_type), POINTER :: para_env
457 : INTEGER, INTENT(IN) :: master, output_unit
458 :
459 2147 : IF (gopt_env%eval_opt_geo) THEN
460 1217 : IF (.NOT. gopt_env%dimer_rotation) THEN
461 : CALL write_final_info(output_unit, conv, its, gopt_env, x0, master, &
462 1085 : para_env, force_env, gopt_env%motion_section, root_section)
463 : ELSE
464 132 : CALL update_dimer_vec(gopt_env%dimer_env, gopt_env%motion_section)
465 132 : CALL write_restart(force_env=force_env, root_section=root_section)
466 : END IF
467 : END IF
468 :
469 2147 : END SUBROUTINE gopt_f_io_finalize
470 :
471 : ! **************************************************************************************************
472 : !> \brief ...
473 : !> \param output_unit ...
474 : !> \param it ...
475 : !> \param etot ...
476 : !> \param ediff ...
477 : !> \param pred ...
478 : !> \param rat ...
479 : !> \param step ...
480 : !> \param rad ...
481 : !> \param emin ...
482 : !> \param pres_int ...
483 : !> \param wildcard ...
484 : !> \param used_time ...
485 : ! **************************************************************************************************
486 14319 : SUBROUTINE write_cycle_infos(output_unit, it, etot, ediff, pred, rat, step, rad, emin, &
487 : pres_int, wildcard, used_time, max_memory, energy_unit, stress_unit)
488 :
489 : INTEGER, INTENT(IN) :: output_unit, it
490 : REAL(KIND=dp), INTENT(IN) :: etot
491 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: ediff, pred, rat, step, rad, emin, &
492 : pres_int
493 : CHARACTER(LEN=5), INTENT(IN) :: wildcard
494 : REAL(KIND=dp), INTENT(IN) :: used_time
495 : INTEGER(KIND=int_8), INTENT(IN) :: max_memory
496 : CHARACTER(LEN=default_string_length), INTENT(IN) :: energy_unit, stress_unit
497 :
498 : CHARACTER(LEN=5) :: tag
499 :
500 14319 : IF (output_unit > 0) THEN
501 7215 : tag = "OPT| "
502 7215 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") tag//REPEAT("*", 74)
503 : WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,I25)") &
504 7215 : tag//"Step number", it
505 : WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,A25)") &
506 7215 : tag//"Optimization method", wildcard
507 : WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
508 7215 : tag//"Total energy ["//TRIM(ADJUSTL(energy_unit))//"]", &
509 14430 : cp_unit_from_cp2k(etot, TRIM(energy_unit))
510 7215 : IF (PRESENT(pres_int)) THEN
511 : WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
512 2260 : tag//"Internal pressure ["//TRIM(ADJUSTL(stress_unit))//"]", &
513 4520 : cp_unit_from_cp2k(pres_int, TRIM(stress_unit))
514 : END IF
515 7215 : IF (PRESENT(ediff)) THEN
516 : WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
517 6399 : tag//"Effective energy change ["//TRIM(ADJUSTL(energy_unit))//"]", &
518 12798 : cp_unit_from_cp2k(ediff, TRIM(energy_unit))
519 : END IF
520 7215 : IF (PRESENT(pred)) THEN
521 : WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
522 3312 : tag//"Predicted energy change ["//TRIM(ADJUSTL(energy_unit))//"]", &
523 6624 : cp_unit_from_cp2k(pred, TRIM(energy_unit))
524 : END IF
525 7215 : IF (PRESENT(rat)) THEN
526 : WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
527 3312 : tag//"Scaling factor", rat
528 : END IF
529 7215 : IF (PRESENT(step)) THEN
530 : WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
531 3312 : tag//"Step size", step
532 : END IF
533 7215 : IF (PRESENT(rad)) THEN
534 : WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
535 3312 : tag//"Trust radius", rad
536 : END IF
537 7215 : IF (PRESENT(emin)) THEN
538 6399 : IF (etot < emin) THEN
539 : WRITE (UNIT=output_unit, FMT="(T2,A,T77,A4)") &
540 5835 : tag//"Decrease in energy", " YES"
541 : ELSE
542 : WRITE (UNIT=output_unit, FMT="(T2,A,T77,A4)") &
543 564 : tag//"Decrease in energy", " NO"
544 : END IF
545 : END IF
546 : WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.3)") &
547 7215 : tag//"Used time [s]", used_time
548 7215 : IF (it == 0) THEN
549 810 : WRITE (UNIT=output_unit, FMT="(T2,A)") tag//REPEAT("*", 74)
550 810 : IF (max_memory /= 0) THEN
551 : WRITE (UNIT=output_unit, FMT="(T2,A,T60,1X,I20)") &
552 810 : tag//"Estimated peak process memory [MiB]", &
553 1620 : (max_memory + (1024*1024) - 1)/(1024*1024)
554 : END IF
555 : END IF
556 : END IF
557 :
558 14319 : END SUBROUTINE write_cycle_infos
559 :
560 : ! **************************************************************************************************
561 : !> \brief ...
562 : !> \param output_unit ...
563 : !> \param it ...
564 : !> \param etot ...
565 : !> \param ediff ...
566 : !> \param emin ...
567 : !> \param dimer_env ...
568 : !> \param used_time ...
569 : !> \param wildcard ...
570 : !> \date 01.2008
571 : !> \author Luca Bellucci and Teodoro Laino - created [tlaino]
572 : ! **************************************************************************************************
573 858 : SUBROUTINE write_rot_cycle_infos(output_unit, it, etot, ediff, emin, dimer_env, used_time, &
574 : wildcard, max_memory)
575 :
576 : INTEGER, INTENT(IN) :: output_unit, it
577 : REAL(KIND=dp), INTENT(IN) :: etot
578 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: ediff, emin
579 : TYPE(dimer_env_type), POINTER :: dimer_env
580 : REAL(KIND=dp), INTENT(IN) :: used_time
581 : CHARACTER(LEN=5), INTENT(IN) :: wildcard
582 : INTEGER(KIND=int_8), INTENT(IN) :: max_memory
583 :
584 : CHARACTER(LEN=5) :: tag
585 :
586 858 : IF (output_unit > 0) THEN
587 429 : tag = "OPT| "
588 429 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") tag//REPEAT("*", 74)
589 : WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,I25)") &
590 429 : tag//"Rotational step number", it
591 : WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,A25)") &
592 429 : tag//"Optimization method", wildcard
593 : WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
594 429 : tag//"Local curvature", dimer_env%rot%curvature, &
595 858 : tag//"Total rotational force", etot
596 429 : IF (PRESENT(ediff)) THEN
597 : WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
598 363 : tag//"Rotational force change", ediff
599 : END IF
600 429 : IF (PRESENT(emin)) THEN
601 363 : IF (etot < emin) THEN
602 : WRITE (UNIT=output_unit, FMT="(T2,A,T77,A4)") &
603 157 : tag//"Decrease in rotational force", " YES"
604 : ELSE
605 : WRITE (UNIT=output_unit, FMT="(T2,A,T77,A4)") &
606 206 : tag//"Decrease in rotational force", " NO"
607 : END IF
608 : END IF
609 : WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.3)") &
610 429 : tag//"Used time [s]", used_time
611 429 : IF (it == 0) THEN
612 66 : WRITE (UNIT=output_unit, FMT="(T2,A)") tag//REPEAT("*", 74)
613 66 : IF (max_memory /= 0) THEN
614 : WRITE (UNIT=output_unit, FMT="(T2,A,T60,1X,I20)") &
615 66 : tag//"Estimated peak process memory [MiB]", &
616 132 : (max_memory + (1024*1024) - 1)/(1024*1024)
617 : END IF
618 : END IF
619 : END IF
620 :
621 858 : END SUBROUTINE write_rot_cycle_infos
622 :
623 : ! **************************************************************************************************
624 : !> \brief ...
625 : !> \param ndf ...
626 : !> \param dr ...
627 : !> \param g ...
628 : !> \param output_unit ...
629 : !> \param conv ...
630 : !> \param gopt_param ...
631 : !> \param max_memory ...
632 : !> \param pres_diff ...
633 : !> \param pres_tol ...
634 : !> \param pres_diff_constr ...
635 : ! **************************************************************************************************
636 12722 : SUBROUTINE check_converg(ndf, dr, g, output_unit, conv, gopt_param, max_memory, stress_unit, &
637 : pres_diff, pres_tol, pres_diff_constr)
638 :
639 : INTEGER, INTENT(IN) :: ndf
640 : REAL(KIND=dp), INTENT(IN) :: dr(ndf), g(ndf)
641 : INTEGER, INTENT(IN) :: output_unit
642 : LOGICAL, INTENT(OUT) :: conv
643 : TYPE(gopt_param_type), POINTER :: gopt_param
644 : INTEGER(KIND=int_8), INTENT(IN) :: max_memory
645 : CHARACTER(LEN=default_string_length), INTENT(IN) :: stress_unit
646 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: pres_diff, pres_tol, pres_diff_constr
647 :
648 : CHARACTER(LEN=5) :: tag
649 : INTEGER :: indf
650 : LOGICAL :: conv_dx, conv_g, conv_p, conv_rdx, &
651 : conv_rg
652 : REAL(KIND=dp) :: dumm, dxcon, gcon, maxdum(4), rmsgcon, &
653 : rmsxcon
654 :
655 12722 : dxcon = gopt_param%max_dr
656 12722 : gcon = gopt_param%max_force
657 12722 : rmsgcon = gopt_param%rms_force
658 12722 : rmsxcon = gopt_param%rms_dr
659 :
660 12722 : conv = .FALSE.
661 12722 : conv_dx = .TRUE.
662 12722 : conv_rdx = .TRUE.
663 12722 : conv_g = .TRUE.
664 12722 : conv_rg = .TRUE.
665 12722 : conv_p = .TRUE.
666 :
667 12722 : dumm = 0.0_dp
668 2954411 : DO indf = 1, ndf
669 2941689 : IF (indf == 1) maxdum(1) = ABS(dr(indf))
670 2941689 : dumm = dumm + dr(indf)**2
671 2941689 : IF (ABS(dr(indf)) > dxcon) conv_dx = .FALSE.
672 2954411 : IF (ABS(dr(indf)) > maxdum(1)) maxdum(1) = ABS(dr(indf))
673 : END DO
674 : ! SQRT(dumm/ndf) > rmsxcon
675 12722 : IF (dumm > (rmsxcon*rmsxcon*ndf)) conv_rdx = .FALSE.
676 12722 : maxdum(2) = SQRT(dumm/ndf)
677 :
678 12722 : dumm = 0.0_dp
679 2954411 : DO indf = 1, ndf
680 2941689 : IF (indf == 1) maxdum(3) = ABS(g(indf))
681 2941689 : dumm = dumm + g(indf)**2
682 2941689 : IF (ABS(g(indf)) > gcon) conv_g = .FALSE.
683 2954411 : IF (ABS(g(indf)) > maxdum(3)) maxdum(3) = ABS(g(indf))
684 : END DO
685 : ! SQRT(dumm/ndf) > rmsgcon
686 12722 : IF (dumm > (rmsgcon*rmsgcon*ndf)) conv_rg = .FALSE.
687 12722 : maxdum(4) = SQRT(dumm/ndf)
688 :
689 12722 : IF (PRESENT(pres_diff_constr) .AND. PRESENT(pres_tol)) THEN
690 18 : conv_p = ABS(pres_diff_constr) < ABS(pres_tol)
691 12704 : ELSEIF (PRESENT(pres_diff) .AND. PRESENT(pres_tol)) THEN
692 4326 : conv_p = ABS(pres_diff) < ABS(pres_tol)
693 : END IF
694 :
695 12722 : IF (output_unit > 0) THEN
696 :
697 6399 : tag = "OPT| "
698 :
699 6399 : WRITE (UNIT=output_unit, FMT="(T2,A)") TRIM(tag)
700 : WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
701 6399 : tag//"Maximum step size", maxdum(1), &
702 12798 : tag//"Convergence limit for maximum step size", dxcon
703 6399 : IF (conv_dx) THEN
704 : WRITE (UNIT=output_unit, FMT="(T2,A,T77,A4)") &
705 1458 : tag//"Maximum step size is converged", " YES"
706 : ELSE
707 : WRITE (UNIT=output_unit, FMT="(T2,A,T77,A4)") &
708 4941 : tag//"Maximum step size is converged", " NO"
709 : END IF
710 :
711 6399 : WRITE (UNIT=output_unit, FMT="(T2,A)") TRIM(tag)
712 : WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
713 6399 : tag//"RMS step size", maxdum(2), &
714 12798 : tag//"Convergence limit for RMS step size", rmsxcon
715 6399 : IF (conv_rdx) THEN
716 : WRITE (UNIT=output_unit, FMT="(T2,A,T77,A4)") &
717 1588 : tag//"RMS step size is converged", " YES"
718 : ELSE
719 : WRITE (UNIT=output_unit, FMT="(T2,A,T77,A4)") &
720 4811 : tag//"RMS step size is converged", " NO"
721 : END IF
722 :
723 6399 : WRITE (UNIT=output_unit, FMT="(T2,A)") TRIM(tag)
724 : WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
725 6399 : tag//"Maximum gradient", maxdum(3), &
726 12798 : tag//"Convergence limit for maximum gradient", gcon
727 6399 : IF (conv_g) THEN
728 : WRITE (UNIT=output_unit, FMT="(T2,A,T77,A4)") &
729 1516 : tag//"Maximum gradient is converged", " YES"
730 : ELSE
731 : WRITE (UNIT=output_unit, FMT="(T2,A,T77,A4)") &
732 4883 : tag//"Maximum gradient is converged", " NO"
733 : END IF
734 :
735 6399 : WRITE (UNIT=output_unit, FMT="(T2,A)") TRIM(tag)
736 : WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
737 6399 : tag//"RMS gradient", maxdum(4), &
738 12798 : tag//"Convergence limit for RMS gradient", rmsgcon
739 6399 : IF (conv_rg) THEN
740 : WRITE (UNIT=output_unit, FMT="(T2,A,T77,A4)") &
741 1799 : tag//"RMS gradient is converged", " YES"
742 : ELSE
743 : WRITE (UNIT=output_unit, FMT="(T2,A,T77,A4)") &
744 4600 : tag//"RMS gradient is converged", " NO"
745 : END IF
746 :
747 6399 : IF (PRESENT(pres_diff) .AND. PRESENT(pres_tol)) THEN
748 2172 : WRITE (UNIT=output_unit, FMT="(T2,A)") TRIM(tag)
749 2172 : IF (PRESENT(pres_diff_constr)) THEN
750 : WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
751 : tag//"Pressure deviation without constraint ["// &
752 9 : TRIM(ADJUSTL(stress_unit))//"]", &
753 18 : cp_unit_from_cp2k(pres_diff, TRIM(stress_unit))
754 : WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
755 : tag//"Pressure deviation with constraint ["// &
756 9 : TRIM(ADJUSTL(stress_unit))//"]", &
757 18 : cp_unit_from_cp2k(pres_diff_constr, TRIM(stress_unit))
758 : ELSE
759 : WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
760 2163 : tag//"Pressure deviation ["//TRIM(ADJUSTL(stress_unit))//"]", &
761 4326 : cp_unit_from_cp2k(pres_diff, TRIM(stress_unit))
762 : END IF
763 : WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
764 2172 : tag//"Pressure tolerance ["//TRIM(ADJUSTL(stress_unit))//"]", &
765 4344 : cp_unit_from_cp2k(pres_tol, TRIM(stress_unit))
766 2172 : IF (conv_p) THEN
767 : WRITE (UNIT=output_unit, FMT="(T2,A,T77,A4)") &
768 309 : tag//"Pressure is converged", " YES"
769 : ELSE
770 : WRITE (UNIT=output_unit, FMT="(T2,A,T77,A4)") &
771 1863 : tag//"Pressure is converged", " NO"
772 : END IF
773 : END IF
774 :
775 6399 : WRITE (UNIT=output_unit, FMT="(T2,A)") tag//REPEAT("*", 74)
776 :
777 6399 : IF (max_memory /= 0) THEN
778 : WRITE (UNIT=output_unit, FMT="(T2,A,T60,1X,I20)") &
779 6399 : tag//"Estimated peak process memory after this step [MiB]", &
780 12798 : (max_memory + (1024*1024) - 1)/(1024*1024)
781 : END IF
782 :
783 : END IF
784 :
785 12722 : IF (conv_dx .AND. conv_rdx .AND. conv_g .AND. conv_rg .AND. conv_p) conv = .TRUE.
786 :
787 12722 : IF ((conv) .AND. (output_unit > 0)) THEN
788 662 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") REPEAT("*", 79)
789 : WRITE (UNIT=output_unit, FMT="(T2,A,T25,A,T78,A)") &
790 662 : "***", "GEOMETRY OPTIMIZATION COMPLETED", "***"
791 662 : WRITE (UNIT=output_unit, FMT="(T2,A)") REPEAT("*", 79)
792 : END IF
793 :
794 12722 : END SUBROUTINE check_converg
795 :
796 : ! **************************************************************************************************
797 : !> \brief ...
798 : !> \param dimer_env ...
799 : !> \param output_unit ...
800 : !> \param conv ...
801 : !> \date 01.2008
802 : !> \author Luca Bellucci and Teodoro Laino - created [tlaino]
803 : ! **************************************************************************************************
804 726 : SUBROUTINE check_rot_conv(dimer_env, output_unit, conv)
805 :
806 : TYPE(dimer_env_type), POINTER :: dimer_env
807 : INTEGER, INTENT(IN) :: output_unit
808 : LOGICAL, INTENT(OUT) :: conv
809 :
810 : CHARACTER(LEN=5) :: tag
811 :
812 726 : conv = (ABS(dimer_env%rot%angle2) < dimer_env%rot%angle_tol)
813 :
814 726 : IF (output_unit > 0) THEN
815 363 : tag = "OPT| "
816 363 : WRITE (UNIT=output_unit, FMT="(T2,A)") TRIM(tag)
817 : WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
818 363 : tag//"Predicted angle step size", dimer_env%rot%angle1, &
819 363 : tag//"Effective angle step size", dimer_env%rot%angle2, &
820 726 : tag//"Convergence limit for angle step size", dimer_env%rot%angle_tol
821 363 : IF (conv) THEN
822 : WRITE (UNIT=output_unit, FMT="(T2,A,T77,A4)") &
823 55 : tag//"Angle step size is converged", " YES"
824 : ELSE
825 : WRITE (UNIT=output_unit, FMT="(T2,A,T77,A4)") &
826 308 : tag//"Angle step size is converged", " NO"
827 : END IF
828 363 : WRITE (UNIT=output_unit, FMT="(T2,A)") tag//REPEAT("*", 74)
829 : END IF
830 :
831 726 : IF ((conv) .AND. (output_unit > 0)) THEN
832 55 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") REPEAT("*", 79)
833 : WRITE (UNIT=output_unit, FMT="(T2,A,T25,A,T78,A)") &
834 55 : "***", "ROTATION OPTIMIZATION COMPLETED", "***"
835 55 : WRITE (UNIT=output_unit, FMT="(T2,A)") REPEAT("*", 79)
836 : END IF
837 :
838 726 : END SUBROUTINE check_rot_conv
839 :
840 : ! **************************************************************************************************
841 : !> \brief ...
842 : !> \param output_unit ...
843 : !> \param conv ...
844 : !> \param it ...
845 : !> \param gopt_env ...
846 : !> \param x0 ...
847 : !> \param master ...
848 : !> \param para_env ...
849 : !> \param force_env ...
850 : !> \param motion_section ...
851 : !> \param root_section ...
852 : !> \date 11.2007
853 : !> \author Teodoro Laino [tlaino] - University of Zurich
854 : ! **************************************************************************************************
855 1085 : RECURSIVE SUBROUTINE write_final_info(output_unit, conv, it, gopt_env, x0, master, para_env, force_env, &
856 : motion_section, root_section)
857 : INTEGER, INTENT(IN) :: output_unit
858 : LOGICAL, INTENT(IN) :: conv
859 : INTEGER, INTENT(INOUT) :: it
860 : TYPE(gopt_f_type), POINTER :: gopt_env
861 : REAL(KIND=dp), DIMENSION(:), POINTER :: x0
862 : INTEGER, INTENT(IN) :: master
863 : TYPE(mp_para_env_type), POINTER :: para_env
864 : TYPE(force_env_type), POINTER :: force_env
865 : TYPE(section_vals_type), POINTER :: motion_section, root_section
866 :
867 : CHARACTER(LEN=4) :: constraint_label
868 : LOGICAL :: keep_angles, keep_symmetry, &
869 : keep_volume
870 : REAL(KIND=dp) :: etot
871 : TYPE(cell_type), POINTER :: cell
872 : TYPE(cp_subsys_type), POINTER :: subsys
873 : TYPE(particle_list_type), POINTER :: particles
874 1085 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
875 :
876 1085 : CALL force_env_get(force_env, cell=cell, subsys=subsys)
877 1085 : CALL cp_subsys_get(subsys=subsys, particles=particles)
878 1085 : particle_set => particles%els
879 :
880 : ! Passing gopt_f_type pointer gopt_env to particle_methods where
881 : ! write_final_cif is defined causes a circular dependency, so it
882 : ! is necessary to get some flags by preprocessing...
883 1085 : keep_angles = .TRUE.
884 1085 : keep_symmetry = .TRUE.
885 1085 : keep_volume = .TRUE.
886 1085 : constraint_label = "NONE"
887 1085 : IF (gopt_env%type_id == default_cell_method_id) THEN
888 224 : keep_angles = gopt_env%cell_env%keep_angles
889 224 : keep_symmetry = gopt_env%cell_env%keep_symmetry
890 224 : keep_volume = gopt_env%cell_env%keep_volume
891 224 : SELECT CASE (gopt_env%cell_env%constraint_id)
892 : CASE (fix_x)
893 0 : constraint_label = " X"
894 : CASE (fix_y)
895 0 : constraint_label = " Y"
896 : CASE (fix_z)
897 2 : constraint_label = " Z"
898 : CASE (fix_xy)
899 2 : constraint_label = " XY"
900 : CASE (fix_xz)
901 0 : constraint_label = " XZ"
902 : CASE (fix_yz)
903 0 : constraint_label = " YZ"
904 : CASE (fix_none)
905 224 : constraint_label = "NONE"
906 : END SELECT
907 : END IF
908 : CALL write_final_cif(particle_set, cell, motion_section, conv, &
909 : keep_angles, keep_symmetry, keep_volume, &
910 1085 : gopt_env%label, constraint_label)
911 :
912 1085 : IF (conv) THEN
913 372 : it = it + 1
914 372 : CALL write_structure_data(particle_set, cell, motion_section)
915 372 : CALL write_restart(force_env=force_env, root_section=root_section)
916 :
917 372 : IF (output_unit > 0) &
918 203 : WRITE (UNIT=output_unit, FMT="(/,T20,' Reevaluating energy at the minimum')")
919 :
920 : CALL cp_eval_at(gopt_env, x0, f=etot, master=master, final_evaluation=.TRUE., &
921 372 : para_env=para_env)
922 372 : CALL write_geo_traj(force_env, root_section, it, etot)
923 : END IF
924 :
925 1085 : END SUBROUTINE write_final_info
926 :
927 : ! **************************************************************************************************
928 : !> \brief Specific driver for dumping trajectory during a GEO_OPT
929 : !> \param force_env ...
930 : !> \param root_section ...
931 : !> \param it ...
932 : !> \param etot ...
933 : !> \date 11.2007
934 : !> \par History
935 : !> 09.2010: Output of core and shell positions and forces (MK)
936 : !> \author Teodoro Laino [tlaino] - University of Zurich
937 : ! **************************************************************************************************
938 25848 : SUBROUTINE write_geo_traj(force_env, root_section, it, etot)
939 :
940 : TYPE(force_env_type), POINTER :: force_env
941 : TYPE(section_vals_type), POINTER :: root_section
942 : INTEGER, INTENT(IN) :: it
943 : REAL(KIND=dp), INTENT(IN) :: etot
944 :
945 : LOGICAL :: shell_adiabatic, shell_present
946 : TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
947 12924 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
948 : TYPE(cp_subsys_type), POINTER :: subsys
949 : TYPE(particle_list_type), POINTER :: core_particles, shell_particles
950 :
951 12924 : NULLIFY (atomic_kinds)
952 12924 : NULLIFY (atomic_kind_set)
953 12924 : NULLIFY (core_particles)
954 12924 : NULLIFY (shell_particles)
955 12924 : NULLIFY (subsys)
956 :
957 12924 : CALL write_trajectory(force_env, root_section, it, 0.0_dp, 0.0_dp, etot)
958 : ! Print Force
959 12924 : CALL write_trajectory(force_env, root_section, it, 0.0_dp, 0.0_dp, etot, "FORCES", middle_name="frc")
960 12924 : CALL force_env_get(force_env, subsys=subsys)
961 12924 : CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds)
962 12924 : atomic_kind_set => atomic_kinds%els
963 : CALL get_atomic_kind_set(atomic_kind_set, &
964 : shell_present=shell_present, &
965 12924 : shell_adiabatic=shell_adiabatic)
966 12924 : IF (shell_present) THEN
967 : CALL cp_subsys_get(subsys, &
968 : core_particles=core_particles, &
969 4320 : shell_particles=shell_particles)
970 : CALL write_trajectory(force_env, root_section, it=it, time=0.0_dp, dtime=0.0_dp, &
971 : etot=etot, pk_name="SHELL_TRAJECTORY", middle_name="shpos", &
972 4320 : particles=shell_particles)
973 4320 : IF (shell_adiabatic) THEN
974 : CALL write_trajectory(force_env, root_section, it=it, time=0.0_dp, dtime=0.0_dp, &
975 : etot=etot, pk_name="SHELL_FORCES", middle_name="shfrc", &
976 4320 : particles=shell_particles)
977 : CALL write_trajectory(force_env, root_section, it=it, time=0.0_dp, dtime=0.0_dp, &
978 : etot=etot, pk_name="CORE_TRAJECTORY", middle_name="copos", &
979 4320 : particles=core_particles)
980 : CALL write_trajectory(force_env, root_section, it=it, time=0.0_dp, dtime=0.0_dp, &
981 : etot=etot, pk_name="CORE_FORCES", middle_name="cofrc", &
982 4320 : particles=core_particles)
983 : END IF
984 : END IF
985 :
986 12924 : END SUBROUTINE write_geo_traj
987 :
988 : ! **************************************************************************************************
989 : !> \brief ...
990 : !> \param gopt_env ...
991 : !> \param output_unit ...
992 : !> \param label ...
993 : !> \date 01.2008
994 : !> \author Teodoro Laino [tlaino] - University of Zurich
995 : ! **************************************************************************************************
996 2147 : SUBROUTINE print_geo_opt_header(gopt_env, output_unit, label)
997 :
998 : TYPE(gopt_f_type), POINTER :: gopt_env
999 : INTEGER, INTENT(IN) :: output_unit
1000 : CHARACTER(LEN=*), INTENT(IN) :: label
1001 :
1002 : CHARACTER(LEN=default_string_length) :: my_format, my_label
1003 : INTEGER :: ix
1004 :
1005 2147 : IF (output_unit > 0) THEN
1006 1091 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") REPEAT("*", 79)
1007 1091 : IF (gopt_env%dimer_rotation) THEN
1008 66 : my_label = "OPTIMIZING DIMER ROTATION"
1009 : ELSE
1010 1025 : my_label = "STARTING "//gopt_env%tag(1:8)//" OPTIMIZATION"
1011 : END IF
1012 :
1013 1091 : ix = (80 - 7 - LEN_TRIM(my_label))/2
1014 1091 : ix = ix + 5
1015 1091 : my_format = "(T2,A,T"//cp_to_string(ix)//",A,T78,A)"
1016 1091 : WRITE (UNIT=output_unit, FMT=TRIM(my_format)) "***", TRIM(my_label), "***"
1017 :
1018 1091 : ix = (80 - 7 - LEN_TRIM(label))/2
1019 1091 : ix = ix + 5
1020 1091 : my_format = "(T2,A,T"//cp_to_string(ix)//",A,T78,A)"
1021 1091 : WRITE (UNIT=output_unit, FMT=TRIM(my_format)) "***", TRIM(label), "***"
1022 :
1023 1091 : WRITE (UNIT=output_unit, FMT="(T2,A)") REPEAT("*", 79)
1024 1091 : CALL m_flush(output_unit)
1025 : END IF
1026 2147 : END SUBROUTINE print_geo_opt_header
1027 :
1028 : ! **************************************************************************************************
1029 : !> \brief ...
1030 : !> \param gopt_env ...
1031 : !> \param output_unit ...
1032 : !> \date 01.2008
1033 : !> \author Teodoro Laino [tlaino] - University of Zurich
1034 : ! **************************************************************************************************
1035 735 : SUBROUTINE print_geo_opt_nc(gopt_env, output_unit)
1036 :
1037 : TYPE(gopt_f_type), POINTER :: gopt_env
1038 : INTEGER, INTENT(IN) :: output_unit
1039 :
1040 735 : IF (output_unit > 0) THEN
1041 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
1042 368 : "*** MAXIMUM NUMBER OF OPTIMIZATION STEPS REACHED ***"
1043 368 : IF (.NOT. gopt_env%dimer_rotation) THEN
1044 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
1045 357 : "*** EXITING GEOMETRY OPTIMIZATION ***"
1046 : ELSE
1047 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
1048 11 : "*** EXITING ROTATION OPTIMIZATION ***"
1049 : END IF
1050 368 : CALL m_flush(output_unit)
1051 : END IF
1052 :
1053 735 : END SUBROUTINE print_geo_opt_nc
1054 :
1055 : ! **************************************************************************************************
1056 : !> \brief Prints information during GEO_OPT common to all optimizers
1057 : !> \param force_env ...
1058 : !> \param root_section ...
1059 : !> \param motion_section ...
1060 : !> \param its ...
1061 : !> \param opt_energy ...
1062 : !> \date 02.2008
1063 : !> \author Teodoro Laino [tlaino] - University of Zurich
1064 : !> \version 1.0
1065 : ! **************************************************************************************************
1066 12552 : SUBROUTINE geo_opt_io(force_env, root_section, motion_section, its, opt_energy)
1067 :
1068 : TYPE(force_env_type), POINTER :: force_env
1069 : TYPE(section_vals_type), POINTER :: root_section, motion_section
1070 : INTEGER, INTENT(IN) :: its
1071 : REAL(KIND=dp), INTENT(IN) :: opt_energy
1072 :
1073 : TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
1074 12552 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1075 : TYPE(cell_type), POINTER :: cell
1076 : TYPE(cp_subsys_type), POINTER :: subsys
1077 : TYPE(distribution_1d_type), POINTER :: local_particles
1078 : TYPE(mp_para_env_type), POINTER :: para_env
1079 : TYPE(particle_list_type), POINTER :: particles
1080 12552 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1081 : TYPE(virial_type), POINTER :: virial
1082 :
1083 12552 : NULLIFY (para_env, atomic_kind_set, subsys, particle_set, &
1084 12552 : local_particles, atomic_kinds, particles)
1085 :
1086 : ! Write Restart File
1087 12552 : CALL write_restart(force_env=force_env, root_section=root_section)
1088 :
1089 : ! Write Trajectory
1090 12552 : CALL write_geo_traj(force_env, root_section, its, opt_energy)
1091 :
1092 : ! Write the stress Tensor
1093 : CALL force_env_get(force_env, cell=cell, para_env=para_env, &
1094 12552 : subsys=subsys)
1095 : CALL cp_subsys_get(subsys=subsys, atomic_kinds=atomic_kinds, local_particles=local_particles, &
1096 12552 : particles=particles, virial=virial)
1097 12552 : atomic_kind_set => atomic_kinds%els
1098 12552 : particle_set => particles%els
1099 : CALL virial_evaluate(atomic_kind_set, particle_set, local_particles, &
1100 12552 : virial, para_env)
1101 12552 : CALL write_stress_tensor_to_file(virial, cell, motion_section, its, 0.0_dp)
1102 :
1103 : ! Write the cell
1104 12552 : CALL write_simulation_cell(cell, motion_section, its, 0.0_dp)
1105 :
1106 12552 : END SUBROUTINE geo_opt_io
1107 :
1108 : ! **************************************************************************************************
1109 : !> \brief Apply coordinate transformations after cell (shape) change
1110 : !> \param gopt_env ...
1111 : !> \param cell ...
1112 : !> \param x ...
1113 : !> \param update_forces ...
1114 : !> \date 05.11.2012 (revised version of unbiase_coordinates moved here, MK)
1115 : !> \author Matthias Krack
1116 : !> \version 1.0
1117 : ! **************************************************************************************************
1118 13150 : SUBROUTINE apply_cell_change(gopt_env, cell, x, update_forces)
1119 :
1120 : TYPE(gopt_f_type), POINTER :: gopt_env
1121 : TYPE(cell_type), POINTER :: cell
1122 : REAL(KIND=dp), DIMENSION(:), POINTER :: x
1123 : LOGICAL, INTENT(IN) :: update_forces
1124 :
1125 : INTEGER :: i, iatom, idg, j, natom, nparticle, &
1126 : shell_index
1127 : REAL(KIND=dp) :: fc, fs, mass
1128 : REAL(KIND=dp), DIMENSION(3) :: s
1129 : TYPE(cell_type), POINTER :: cell_ref
1130 : TYPE(cp_subsys_type), POINTER :: subsys
1131 : TYPE(particle_list_type), POINTER :: core_particles, particles, &
1132 : shell_particles
1133 :
1134 13150 : NULLIFY (cell_ref)
1135 13150 : NULLIFY (core_particles)
1136 13150 : NULLIFY (particles)
1137 13150 : NULLIFY (shell_particles)
1138 13150 : NULLIFY (subsys)
1139 :
1140 13150 : natom = force_env_get_natom(gopt_env%force_env)
1141 13150 : nparticle = force_env_get_nparticle(gopt_env%force_env)
1142 : CALL force_env_get(gopt_env%force_env, &
1143 13150 : subsys=subsys)
1144 : CALL cp_subsys_get(subsys=subsys, &
1145 : core_particles=core_particles, &
1146 : particles=particles, &
1147 13150 : shell_particles=shell_particles)
1148 :
1149 : ! Retrieve the reference cell
1150 13150 : CALL cell_create(cell_ref)
1151 13150 : CALL cell_copy(cell, cell_ref, tag="CELL_OPT_REF")
1152 :
1153 : ! Load the updated cell information
1154 25346 : SELECT CASE (gopt_env%cell_method_id)
1155 : CASE (default_cell_direct_id)
1156 12196 : idg = 3*nparticle
1157 12196 : CALL init_cell(cell_ref, hmat=gopt_env%h_ref)
1158 : CASE (default_cell_geo_opt_id, default_cell_md_id)
1159 13150 : idg = 0
1160 : END SELECT
1161 13150 : CPASSERT((SIZE(x) == idg + 6))
1162 :
1163 13150 : IF (update_forces) THEN
1164 :
1165 : ! Transform particle forces back to reference cell
1166 : idg = 1
1167 256636 : DO iatom = 1, natom
1168 251892 : CALL real_to_scaled(s, x(idg:idg + 2), cell)
1169 251892 : CALL scaled_to_real(x(idg:idg + 2), s, cell_ref)
1170 256636 : idg = idg + 3
1171 : END DO
1172 :
1173 : ELSE
1174 :
1175 : ! Update cell
1176 33624 : DO i = 1, 3
1177 84060 : DO j = 1, i
1178 50436 : idg = idg + 1
1179 75654 : cell%hmat(j, i) = x(idg)
1180 : END DO
1181 : END DO
1182 8406 : CALL init_cell(cell)
1183 8406 : CALL cp_subsys_set(subsys, cell=cell)
1184 :
1185 : ! Retrieve particle coordinates for the current cell
1186 8406 : SELECT CASE (gopt_env%cell_method_id)
1187 : CASE (default_cell_direct_id)
1188 : idg = 1
1189 471282 : DO iatom = 1, natom
1190 463830 : CALL real_to_scaled(s, x(idg:idg + 2), cell_ref)
1191 463830 : shell_index = particles%els(iatom)%shell_index
1192 463830 : IF (shell_index == 0) THEN
1193 143874 : CALL scaled_to_real(particles%els(iatom)%r, s, cell)
1194 : ELSE
1195 319956 : CALL scaled_to_real(core_particles%els(shell_index)%r, s, cell)
1196 319956 : i = 3*(natom + shell_index - 1) + 1
1197 319956 : CALL real_to_scaled(s, x(i:i + 2), cell_ref)
1198 319956 : CALL scaled_to_real(shell_particles%els(shell_index)%r, s, cell)
1199 : ! Update atomic position due to core and shell motion
1200 319956 : mass = particles%els(iatom)%atomic_kind%mass
1201 319956 : fc = core_particles%els(shell_index)%atomic_kind%shell%mass_core/mass
1202 319956 : fs = shell_particles%els(shell_index)%atomic_kind%shell%mass_shell/mass
1203 : particles%els(iatom)%r(1:3) = fc*core_particles%els(shell_index)%r(1:3) + &
1204 2559648 : fs*shell_particles%els(shell_index)%r(1:3)
1205 : END IF
1206 471282 : idg = idg + 3
1207 : END DO
1208 : CASE (default_cell_geo_opt_id, default_cell_md_id)
1209 49056 : DO iatom = 1, natom
1210 39696 : shell_index = particles%els(iatom)%shell_index
1211 40650 : IF (shell_index == 0) THEN
1212 35448 : CALL real_to_scaled(s, particles%els(iatom)%r, cell_ref)
1213 35448 : CALL scaled_to_real(particles%els(iatom)%r, s, cell)
1214 : ELSE
1215 4248 : CALL real_to_scaled(s, core_particles%els(shell_index)%r, cell_ref)
1216 4248 : CALL scaled_to_real(core_particles%els(shell_index)%r, s, cell)
1217 4248 : i = 3*(natom + shell_index - 1) + 1
1218 4248 : CALL real_to_scaled(s, shell_particles%els(shell_index)%r, cell_ref)
1219 4248 : CALL scaled_to_real(shell_particles%els(shell_index)%r, s, cell)
1220 : ! Update atomic position due to core and shell motion
1221 4248 : mass = particles%els(iatom)%atomic_kind%mass
1222 4248 : fc = core_particles%els(shell_index)%atomic_kind%shell%mass_core/mass
1223 4248 : fs = shell_particles%els(shell_index)%atomic_kind%shell%mass_shell/mass
1224 : particles%els(iatom)%r(1:3) = fc*core_particles%els(shell_index)%r(1:3) + &
1225 33984 : fs*shell_particles%els(shell_index)%r(1:3)
1226 : END IF
1227 : END DO
1228 : END SELECT
1229 :
1230 : END IF
1231 :
1232 13150 : CALL cell_release(cell_ref)
1233 :
1234 13150 : END SUBROUTINE apply_cell_change
1235 :
1236 : END MODULE gopt_f_methods
|
