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