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