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