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