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