(git:a90b72f)
Loading...
Searching...
No Matches
gopt_f_methods.F
Go to the documentation of this file.
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! **************************************************************************************************
14MODULE gopt_f_methods
15
16 USE atomic_kind_list_types, ONLY: atomic_kind_list_type
17 USE atomic_kind_types, ONLY: atomic_kind_type, &
18 get_atomic_kind_set
19 USE cell_methods, ONLY: cell_create, &
20 init_cell
21 USE cell_types, ONLY: cell_copy, &
22 cell_release, &
23 cell_type, &
24 real_to_scaled, &
25 scaled_to_real
26 USE cp_log_handling, ONLY: cp_to_string
27 USE cp_subsys_types, ONLY: cp_subsys_get, &
28 cp_subsys_set, &
29 cp_subsys_type, &
30 pack_subsys_particles
31 USE cp_units, ONLY: cp_unit_from_cp2k
32 USE dimer_types, ONLY: dimer_env_type
33 USE dimer_utils, ONLY: update_dimer_vec
34 USE distribution_1d_types, ONLY: distribution_1d_type
35 USE force_env_types, ONLY: force_env_get, &
36 force_env_get_natom, &
37 force_env_get_nparticle, &
38 force_env_type, &
39 use_qmmm, &
40 use_qmmmx
41 USE gopt_f_types, ONLY: gopt_f_type
42 USE gopt_param_types, ONLY: gopt_param_type
43 USE input_constants, ONLY: default_cell_method_id, &
44 default_minimization_method_id, &
45 default_shellcore_method_id, &
46 default_ts_method_id, &
47 fix_none, &
48 fix_x, &
49 fix_xy, &
50 fix_xz, &
51 fix_y, &
52 fix_yz, &
53 fix_z
54 USE input_cp2k_restarts, ONLY: write_restart
55 USE input_section_types, ONLY: section_vals_type, &
56 section_vals_val_get
57 USE kinds, ONLY: default_string_length, &
58 dp, &
59 int_8
60 USE machine, ONLY: m_flush
61 USE md_energies, ONLY: sample_memory
62 USE message_passing, ONLY: mp_para_env_type
63 USE motion_utils, ONLY: write_simulation_cell, &
64 write_stress_tensor_to_file, &
65 write_trajectory
66 USE particle_list_types, ONLY: particle_list_type
67 USE particle_methods, ONLY: write_final_structure, &
68 write_structure_data
69 USE particle_types, ONLY: particle_type
70 USE qmmm_util, ONLY: apply_qmmm_translate
71 USE qmmmx_util, ONLY: apply_qmmmx_translate
72 USE virial_methods, ONLY: virial_evaluate
73 USE virial_types, ONLY: virial_type
74#include "../base/base_uses.f90"
75
76 IMPLICIT NONE
77 PRIVATE
78
79
80! **************************************************************************************************
81!> \brief evaluate the potential energy and its gradients using an array
82!> with same dimension as the particle_set
83!> \param gopt_env the geometry optimization environment
84!> \param x the position where the function should be evaluated
85!> \param f the function value
86!> \param gradient the value of its gradient
87!> \par History
88!> none
89!> \author Teodoro Laino [tlaino] - University of Zurich - 01.2008
90! **************************************************************************************************
91INTERFACE
92
93 SUBROUTINE cp_eval_at(gopt_env, x, f, gradient, master, &
94 final_evaluation, para_env)
95
96 USE message_passing, ONLY: mp_para_env_type
97 USE gopt_f_types, ONLY: gopt_f_type
98 USE kinds, ONLY: dp
99
100 TYPE(gopt_f_type), POINTER :: gopt_env
101 REAL(KIND=dp), DIMENSION(:), POINTER :: x
102 REAL(KIND=dp), INTENT(out), OPTIONAL :: f
103 REAL(KIND=dp), DIMENSION(:), OPTIONAL, &
104 POINTER :: gradient
105 INTEGER, INTENT(IN) :: master
106 LOGICAL, INTENT(IN), OPTIONAL :: final_evaluation
107 TYPE(mp_para_env_type), POINTER :: para_env
108
109 END SUBROUTINE cp_eval_at
110
111END INTERFACE
112
113 LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .true.
114 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = "gopt_f_methods"
115
116 PUBLIC :: gopt_f_create_x0, &
117 print_geo_opt_header, print_geo_opt_nc, &
118 gopt_f_io_init, gopt_f_io, gopt_f_io_finalize, gopt_f_ii, &
119 apply_cell_change
120
121CONTAINS
122
123! **************************************************************************************************
124!> \brief returns the value of the parameters for the actual configuration
125!> \param gopt_env the geometry optimization environment you want the info about
126!> x0: the parameter vector (is allocated by this routine)
127!> \param x0 ...
128!> \par History
129!> - Cell optimization revised (06.11.2012,MK)
130! **************************************************************************************************
131 SUBROUTINE gopt_f_create_x0(gopt_env, x0)
132
133 TYPE(gopt_f_type), POINTER :: gopt_env
134 REAL(kind=dp), DIMENSION(:), POINTER :: x0
135
136 INTEGER :: i, idg, j, nparticle
137 TYPE(cell_type), POINTER :: cell
138 TYPE(cp_subsys_type), POINTER :: subsys
139
140 NULLIFY (cell)
141 NULLIFY (subsys)
142
143 SELECT CASE (gopt_env%type_id)
144 CASE (default_minimization_method_id, default_ts_method_id)
145 CALL force_env_get(gopt_env%force_env, subsys=subsys)
146 ! before starting we handle the case of translating coordinates (QM/MM)
147 IF (gopt_env%force_env%in_use == use_qmmm) &
148 CALL apply_qmmm_translate(gopt_env%force_env%qmmm_env)
149 IF (gopt_env%force_env%in_use == use_qmmmx) &
150 CALL apply_qmmmx_translate(gopt_env%force_env%qmmmx_env)
151 nparticle = force_env_get_nparticle(gopt_env%force_env)
152 ALLOCATE (x0(3*nparticle))
153 CALL pack_subsys_particles(subsys=subsys, r=x0)
154 CASE (default_cell_method_id)
155 CALL force_env_get(gopt_env%force_env, subsys=subsys, cell=cell)
156 ! Store reference cell
157 gopt_env%h_ref = cell%hmat
158 ! before starting we handle the case of translating coordinates (QM/MM)
159 IF (gopt_env%force_env%in_use == use_qmmm) &
160 CALL apply_qmmm_translate(gopt_env%force_env%qmmm_env)
161 IF (gopt_env%force_env%in_use == use_qmmmx) &
162 CALL apply_qmmmx_translate(gopt_env%force_env%qmmmx_env)
163 nparticle = force_env_get_nparticle(gopt_env%force_env)
164 ALLOCATE (x0(3*nparticle + 6))
165 CALL pack_subsys_particles(subsys=subsys, r=x0)
166 idg = 3*nparticle
167 DO i = 1, 3
168 DO j = 1, i
169 idg = idg + 1
170 x0(idg) = cell%hmat(j, i)
171 END DO
172 END DO
173 CASE DEFAULT
174 cpabort("Invalid or not yet implemented type of optimization")
175 END SELECT
176
177 END SUBROUTINE gopt_f_create_x0
178
179! **************************************************************************************************
180!> \brief Prints iteration step of the optimization procedure on screen
181!> \param its ...
182!> \param output_unit ...
183!> \author Teodoro Laino [tlaino] - University of Zurich - 03.2008
184! **************************************************************************************************
185 SUBROUTINE gopt_f_ii(its, output_unit)
186
187 INTEGER, INTENT(IN) :: its, output_unit
188
189 IF (output_unit > 0) THEN
190 WRITE (unit=output_unit, fmt="(/,T2,26('-'))")
191 WRITE (unit=output_unit, fmt="(T2,A,I6)") "OPTIMIZATION STEP: ", its
192 WRITE (unit=output_unit, fmt="(T2,26('-'))")
193 CALL m_flush(output_unit)
194 END IF
195
196 END SUBROUTINE gopt_f_ii
197
198! **************************************************************************************************
199!> \brief Handles the Output during an optimization run
200!> \param gopt_env ...
201!> \param output_unit ...
202!> \param opt_energy ...
203!> \param wildcard ...
204!> \param its ...
205!> \param used_time ...
206!> \author Teodoro Laino [tlaino] - University of Zurich - 03.2008
207! **************************************************************************************************
208 SUBROUTINE gopt_f_io_init(gopt_env, output_unit, opt_energy, wildcard, its, used_time)
209
210 TYPE(gopt_f_type), POINTER :: gopt_env
211 INTEGER, INTENT(IN) :: output_unit
212 REAL(kind=dp) :: opt_energy
213 CHARACTER(LEN=5) :: wildcard
214 INTEGER, INTENT(IN) :: its
215 REAL(kind=dp) :: used_time
216
217 TYPE(mp_para_env_type), POINTER :: para_env
218 CHARACTER(LEN=default_string_length) :: energy_unit, stress_unit
219 REAL(kind=dp) :: pres_int
220 INTEGER(KIND=int_8) :: max_memory
221 LOGICAL :: print_memory
222
223 NULLIFY (para_env)
224 CALL section_vals_val_get(gopt_env%motion_section, "PRINT%MEMORY_INFO", l_val=print_memory)
225 max_memory = 0
226 IF (print_memory) THEN
227 CALL force_env_get(gopt_env%force_env, para_env=para_env)
228 max_memory = sample_memory(para_env)
229 END IF
230
231 CALL section_vals_val_get(gopt_env%force_env%force_env_section, &
232 "PRINT%PROGRAM_RUN_INFO%ENERGY_UNIT", &
233 c_val=energy_unit)
234 CALL section_vals_val_get(gopt_env%force_env%force_env_section, &
235 "PRINT%STRESS_TENSOR%STRESS_UNIT", &
236 c_val=stress_unit)
237
238 SELECT CASE (gopt_env%type_id)
239 CASE (default_ts_method_id, default_minimization_method_id)
240 ! Geometry Optimization (Minimization and Transition State Search)
241 IF (.NOT. gopt_env%dimer_rotation) THEN
242 CALL write_cycle_infos(output_unit, &
243 it=its, &
244 etot=opt_energy, &
245 wildcard=wildcard, &
246 used_time=used_time, &
247 max_memory=max_memory, &
248 energy_unit=energy_unit, &
249 stress_unit=stress_unit)
250 ELSE
251 CALL write_rot_cycle_infos(output_unit, &
252 it=its, &
253 etot=opt_energy, &
254 dimer_env=gopt_env%dimer_env, &
255 wildcard=wildcard, &
256 used_time=used_time, &
257 max_memory=max_memory)
258 END IF
259 CASE (default_cell_method_id)
260 ! Cell Optimization
261 pres_int = gopt_env%cell_env%pres_int
262 CALL write_cycle_infos(output_unit, &
263 it=its, &
264 etot=opt_energy, &
265 pres_int=pres_int, &
266 wildcard=wildcard, &
267 used_time=used_time, &
268 max_memory=max_memory, &
269 energy_unit=energy_unit, &
270 stress_unit=stress_unit)
271 CASE (default_shellcore_method_id)
272 CALL write_cycle_infos(output_unit, &
273 it=its, &
274 etot=opt_energy, &
275 wildcard=wildcard, &
276 used_time=used_time, &
277 max_memory=max_memory, &
278 energy_unit=energy_unit, &
279 stress_unit=stress_unit)
280 END SELECT
281
282 END SUBROUTINE gopt_f_io_init
283
284! **************************************************************************************************
285!> \brief Handles the Output during an optimization run
286!> \param gopt_env ...
287!> \param force_env ...
288!> \param root_section ...
289!> \param its ...
290!> \param opt_energy ...
291!> \param output_unit ...
292!> \param eold ...
293!> \param emin ...
294!> \param wildcard ...
295!> \param gopt_param ...
296!> \param ndf ...
297!> \param dx ...
298!> \param xi ...
299!> \param conv ...
300!> \param pred ...
301!> \param rat ...
302!> \param step ...
303!> \param rad ...
304!> \param used_time ...
305!> \author Teodoro Laino [tlaino] - University of Zurich - 03.2008
306! **************************************************************************************************
307 SUBROUTINE gopt_f_io(gopt_env, force_env, root_section, its, opt_energy, &
308 output_unit, eold, emin, wildcard, gopt_param, ndf, dx, xi, conv, pred, rat, &
309 step, rad, used_time)
310
311 TYPE(gopt_f_type), POINTER :: gopt_env
312 TYPE(force_env_type), POINTER :: force_env
313 TYPE(section_vals_type), POINTER :: root_section
314 INTEGER, INTENT(IN) :: its
315 REAL(kind=dp), INTENT(IN) :: opt_energy
316 INTEGER, INTENT(IN) :: output_unit
317 REAL(kind=dp) :: eold, emin
318 CHARACTER(LEN=5) :: wildcard
319 TYPE(gopt_param_type), POINTER :: gopt_param
320 INTEGER, INTENT(IN), OPTIONAL :: ndf
321 REAL(kind=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: dx
322 REAL(kind=dp), DIMENSION(:), OPTIONAL, POINTER :: xi
323 LOGICAL, OPTIONAL :: conv
324 REAL(kind=dp), INTENT(IN), OPTIONAL :: pred, rat, step, rad
325 REAL(kind=dp) :: used_time
326
327 CHARACTER(LEN=default_string_length) :: energy_unit, stress_unit
328 INTEGER(KIND=int_8) :: max_memory
329 LOGICAL :: print_memory
330 REAL(kind=dp) :: pres_diff, pres_diff_constr, pres_int, &
331 pres_tol
332 TYPE(mp_para_env_type), POINTER :: para_env
333
334 NULLIFY (para_env)
335 CALL section_vals_val_get(gopt_env%motion_section, "PRINT%MEMORY_INFO", l_val=print_memory)
336 max_memory = 0
337 IF (print_memory) THEN
338 CALL force_env_get(force_env, para_env=para_env)
339 max_memory = sample_memory(para_env)
340 END IF
341
342 CALL section_vals_val_get(gopt_env%force_env%force_env_section, &
343 "PRINT%PROGRAM_RUN_INFO%ENERGY_UNIT", &
344 c_val=energy_unit)
345 CALL section_vals_val_get(gopt_env%force_env%force_env_section, &
346 "PRINT%STRESS_TENSOR%STRESS_UNIT", &
347 c_val=stress_unit)
348
349 SELECT CASE (gopt_env%type_id)
350 CASE (default_ts_method_id, default_minimization_method_id)
351 ! Geometry Optimization (Minimization and Transition State Search)
352 IF (.NOT. gopt_env%dimer_rotation) THEN
353 CALL geo_opt_io(force_env=force_env, root_section=root_section, &
354 motion_section=gopt_env%motion_section, its=its, opt_energy=opt_energy)
355 CALL write_cycle_infos(output_unit, &
356 it=its, &
357 etot=opt_energy, &
358 ediff=(opt_energy - eold), &
359 pred=pred, &
360 rat=rat, &
361 step=step, &
362 rad=rad, &
363 emin=emin, &
364 wildcard=wildcard, &
365 used_time=used_time, &
366 max_memory=max_memory, &
367 energy_unit=energy_unit, &
368 stress_unit=stress_unit)
369 ! Possibly check convergence
370 IF (PRESENT(conv)) THEN
371 cpassert(PRESENT(ndf))
372 cpassert(PRESENT(dx))
373 cpassert(PRESENT(xi))
374 CALL check_converg(ndf, dx, xi, output_unit, conv, gopt_param, max_memory, stress_unit)
375 END IF
376 ELSE
377 CALL update_dimer_vec(gopt_env%dimer_env, gopt_env%motion_section)
378 CALL write_restart(force_env=force_env, root_section=root_section)
379 CALL write_rot_cycle_infos(output_unit, its, opt_energy, opt_energy - eold, emin, gopt_env%dimer_env, &
380 wildcard=wildcard, used_time=used_time, max_memory=max_memory)
381 ! Possibly check convergence
382 IF (PRESENT(conv)) THEN
383 cpassert(ASSOCIATED(gopt_env%dimer_env))
384 CALL check_rot_conv(gopt_env%dimer_env, output_unit, conv)
385 END IF
386 END IF
387 CASE (default_cell_method_id)
388 ! Cell Optimization
389 pres_diff = gopt_env%cell_env%pres_int - gopt_env%cell_env%pres_ext
390 pres_int = gopt_env%cell_env%pres_int
391 pres_tol = gopt_env%cell_env%pres_tol
392 CALL geo_opt_io(force_env=force_env, root_section=root_section, &
393 motion_section=gopt_env%motion_section, its=its, opt_energy=opt_energy)
394 CALL write_cycle_infos(output_unit, &
395 it=its, &
396 etot=opt_energy, &
397 ediff=(opt_energy - eold), &
398 pred=pred, &
399 rat=rat, &
400 step=step, &
401 rad=rad, &
402 emin=emin, &
403 pres_int=pres_int, &
404 wildcard=wildcard, &
405 used_time=used_time, &
406 max_memory=max_memory, &
407 energy_unit=energy_unit, &
408 stress_unit=stress_unit)
409 ! Possibly check convergence
410 IF (PRESENT(conv)) THEN
411 cpassert(PRESENT(ndf))
412 cpassert(PRESENT(dx))
413 cpassert(PRESENT(xi))
414 IF (gopt_env%cell_env%constraint_id == fix_none) THEN
415 CALL check_converg(ndf, dx, xi, output_unit, conv, gopt_param, max_memory, stress_unit, &
416 pres_diff, pres_tol)
417 ELSE
418 pres_diff_constr = gopt_env%cell_env%pres_constr - gopt_env%cell_env%pres_ext
419 CALL check_converg(ndf, dx, xi, output_unit, conv, gopt_param, max_memory, stress_unit, &
420 pres_diff, pres_tol, pres_diff_constr)
421 END IF
422 END IF
423 CASE (default_shellcore_method_id)
424 CALL write_cycle_infos(output_unit, &
425 it=its, &
426 etot=opt_energy, &
427 ediff=(opt_energy - eold), &
428 pred=pred, &
429 rat=rat, &
430 step=step, &
431 rad=rad, &
432 emin=emin, &
433 wildcard=wildcard, &
434 used_time=used_time, &
435 max_memory=max_memory, &
436 energy_unit=energy_unit, &
437 stress_unit=stress_unit)
438 ! Possibly check convergence
439 IF (PRESENT(conv)) THEN
440 cpassert(PRESENT(ndf))
441 cpassert(PRESENT(dx))
442 cpassert(PRESENT(xi))
443 CALL check_converg(ndf, dx, xi, output_unit, conv, gopt_param, max_memory, stress_unit)
444 END IF
445 END SELECT
446
447 END SUBROUTINE gopt_f_io
448
449! **************************************************************************************************
450!> \brief Handles the Output at the end of an optimization run
451!> \param gopt_env ...
452!> \param force_env ...
453!> \param x0 ...
454!> \param conv ...
455!> \param its ...
456!> \param root_section ...
457!> \param para_env ...
458!> \param master ...
459!> \param output_unit ...
460!> \author Teodoro Laino [tlaino] - University of Zurich - 03.2008
461! **************************************************************************************************
462 RECURSIVE SUBROUTINE gopt_f_io_finalize(gopt_env, force_env, x0, conv, its, root_section, &
463 para_env, master, output_unit)
464 TYPE(gopt_f_type), POINTER :: gopt_env
465 TYPE(force_env_type), POINTER :: force_env
466 REAL(kind=dp), DIMENSION(:), POINTER :: x0
467 LOGICAL :: conv
468 INTEGER :: its
469 TYPE(section_vals_type), POINTER :: root_section
470 TYPE(mp_para_env_type), POINTER :: para_env
471 INTEGER, INTENT(IN) :: master, output_unit
472
473 IF (gopt_env%eval_opt_geo) THEN
474 IF (.NOT. gopt_env%dimer_rotation) THEN
475 CALL write_final_info(output_unit, conv, its, gopt_env, x0, master, &
476 para_env, force_env, gopt_env%motion_section, root_section)
477 ELSE
478 CALL update_dimer_vec(gopt_env%dimer_env, gopt_env%motion_section)
479 CALL write_restart(force_env=force_env, root_section=root_section)
480 END IF
481 END IF
482
483 END SUBROUTINE gopt_f_io_finalize
484
485! **************************************************************************************************
486!> \brief ...
487!> \param output_unit ...
488!> \param it ...
489!> \param etot ...
490!> \param ediff ...
491!> \param pred ...
492!> \param rat ...
493!> \param step ...
494!> \param rad ...
495!> \param emin ...
496!> \param pres_int ...
497!> \param wildcard ...
498!> \param used_time ...
499! **************************************************************************************************
500 SUBROUTINE write_cycle_infos(output_unit, it, etot, ediff, pred, rat, step, rad, emin, &
501 pres_int, wildcard, used_time, max_memory, energy_unit, stress_unit)
502
503 INTEGER, INTENT(IN) :: output_unit, it
504 REAL(kind=dp), INTENT(IN) :: etot
505 REAL(kind=dp), INTENT(IN), OPTIONAL :: ediff, pred, rat, step, rad, emin, &
506 pres_int
507 CHARACTER(LEN=5), INTENT(IN) :: wildcard
508 REAL(kind=dp), INTENT(IN) :: used_time
509 INTEGER(KIND=int_8), INTENT(IN) :: max_memory
510 CHARACTER(LEN=default_string_length), INTENT(IN) :: energy_unit, stress_unit
511
512 CHARACTER(LEN=5) :: tag
513
514 IF (output_unit > 0) THEN
515 tag = "OPT| "
516 WRITE (unit=output_unit, fmt="(/,T2,A)") tag//repeat("*", 74)
517 WRITE (unit=output_unit, fmt="(T2,A,T55,1X,I25)") &
518 tag//"Step number", it
519 WRITE (unit=output_unit, fmt="(T2,A,T55,1X,A25)") &
520 tag//"Optimization method", wildcard
521 WRITE (unit=output_unit, fmt="(T2,A,T55,1X,F25.10)") &
522 tag//"Total energy ["//trim(adjustl(energy_unit))//"]", &
523 cp_unit_from_cp2k(etot, trim(energy_unit))
524 IF (PRESENT(pres_int)) THEN
525 WRITE (unit=output_unit, fmt="(T2,A,T55,1X,F25.10)") &
526 tag//"Internal pressure ["//trim(adjustl(stress_unit))//"]", &
527 cp_unit_from_cp2k(pres_int, trim(stress_unit))
528 END IF
529 IF (PRESENT(ediff)) THEN
530 WRITE (unit=output_unit, fmt="(T2,A,T55,1X,F25.10)") &
531 tag//"Effective energy change ["//trim(adjustl(energy_unit))//"]", &
532 cp_unit_from_cp2k(ediff, trim(energy_unit))
533 END IF
534 IF (PRESENT(pred)) THEN
535 WRITE (unit=output_unit, fmt="(T2,A,T55,1X,F25.10)") &
536 tag//"Predicted energy change ["//trim(adjustl(energy_unit))//"]", &
537 cp_unit_from_cp2k(pred, trim(energy_unit))
538 END IF
539 IF (PRESENT(rat)) THEN
540 WRITE (unit=output_unit, fmt="(T2,A,T55,1X,F25.10)") &
541 tag//"Scaling factor", rat
542 END IF
543 IF (PRESENT(step)) THEN
544 WRITE (unit=output_unit, fmt="(T2,A,T55,1X,F25.10)") &
545 tag//"Step size", step
546 END IF
547 IF (PRESENT(rad)) THEN
548 WRITE (unit=output_unit, fmt="(T2,A,T55,1X,F25.10)") &
549 tag//"Trust radius", rad
550 END IF
551 IF (PRESENT(emin)) THEN
552 IF (etot < emin) THEN
553 WRITE (unit=output_unit, fmt="(T2,A,T77,A4)") &
554 tag//"Decrease in energy", " YES"
555 ELSE
556 WRITE (unit=output_unit, fmt="(T2,A,T77,A4)") &
557 tag//"Decrease in energy", " NO"
558 END IF
559 END IF
560 WRITE (unit=output_unit, fmt="(T2,A,T55,1X,F25.3)") &
561 tag//"Used time [s]", used_time
562 IF (it == 0) THEN
563 WRITE (unit=output_unit, fmt="(T2,A)") tag//repeat("*", 74)
564 IF (max_memory /= 0) THEN
565 WRITE (unit=output_unit, fmt="(T2,A,T60,1X,I20)") &
566 tag//"Estimated peak process memory [MiB]", &
567 (max_memory + (1024*1024) - 1)/(1024*1024)
568 END IF
569 END IF
570 END IF
571
572 END SUBROUTINE write_cycle_infos
573
574! **************************************************************************************************
575!> \brief ...
576!> \param output_unit ...
577!> \param it ...
578!> \param etot ...
579!> \param ediff ...
580!> \param emin ...
581!> \param dimer_env ...
582!> \param used_time ...
583!> \param wildcard ...
584!> \date 01.2008
585!> \author Luca Bellucci and Teodoro Laino - created [tlaino]
586! **************************************************************************************************
587 SUBROUTINE write_rot_cycle_infos(output_unit, it, etot, ediff, emin, dimer_env, used_time, &
588 wildcard, max_memory)
589
590 INTEGER, INTENT(IN) :: output_unit, it
591 REAL(kind=dp), INTENT(IN) :: etot
592 REAL(kind=dp), INTENT(IN), OPTIONAL :: ediff, emin
593 TYPE(dimer_env_type), POINTER :: dimer_env
594 REAL(kind=dp), INTENT(IN) :: used_time
595 CHARACTER(LEN=5), INTENT(IN) :: wildcard
596 INTEGER(KIND=int_8), INTENT(IN) :: max_memory
597
598 CHARACTER(LEN=5) :: tag
599
600 IF (output_unit > 0) THEN
601 tag = "OPT| "
602 WRITE (unit=output_unit, fmt="(/,T2,A)") tag//repeat("*", 74)
603 WRITE (unit=output_unit, fmt="(T2,A,T55,1X,I25)") &
604 tag//"Rotational step number", it
605 WRITE (unit=output_unit, fmt="(T2,A,T55,1X,A25)") &
606 tag//"Optimization method", wildcard
607 WRITE (unit=output_unit, fmt="(T2,A,T55,1X,F25.10)") &
608 tag//"Local curvature", dimer_env%rot%curvature, &
609 tag//"Total rotational force", etot
610 IF (PRESENT(ediff)) THEN
611 WRITE (unit=output_unit, fmt="(T2,A,T55,1X,F25.10)") &
612 tag//"Rotational force change", ediff
613 END IF
614 IF (PRESENT(emin)) THEN
615 IF (etot < emin) THEN
616 WRITE (unit=output_unit, fmt="(T2,A,T77,A4)") &
617 tag//"Decrease in rotational force", " YES"
618 ELSE
619 WRITE (unit=output_unit, fmt="(T2,A,T77,A4)") &
620 tag//"Decrease in rotational force", " NO"
621 END IF
622 END IF
623 WRITE (unit=output_unit, fmt="(T2,A,T55,1X,F25.3)") &
624 tag//"Used time [s]", used_time
625 IF (it == 0) THEN
626 WRITE (unit=output_unit, fmt="(T2,A)") tag//repeat("*", 74)
627 IF (max_memory /= 0) THEN
628 WRITE (unit=output_unit, fmt="(T2,A,T60,1X,I20)") &
629 tag//"Estimated peak process memory [MiB]", &
630 (max_memory + (1024*1024) - 1)/(1024*1024)
631 END IF
632 END IF
633 END IF
634
635 END SUBROUTINE write_rot_cycle_infos
636
637! **************************************************************************************************
638!> \brief ...
639!> \param ndf ...
640!> \param dr ...
641!> \param g ...
642!> \param output_unit ...
643!> \param conv ...
644!> \param gopt_param ...
645!> \param max_memory ...
646!> \param pres_diff ...
647!> \param pres_tol ...
648!> \param pres_diff_constr ...
649! **************************************************************************************************
650 SUBROUTINE check_converg(ndf, dr, g, output_unit, conv, gopt_param, max_memory, stress_unit, &
651 pres_diff, pres_tol, pres_diff_constr)
652
653 INTEGER, INTENT(IN) :: ndf
654 REAL(kind=dp), INTENT(IN) :: dr(ndf), g(ndf)
655 INTEGER, INTENT(IN) :: output_unit
656 LOGICAL, INTENT(OUT) :: conv
657 TYPE(gopt_param_type), POINTER :: gopt_param
658 INTEGER(KIND=int_8), INTENT(IN) :: max_memory
659 CHARACTER(LEN=default_string_length), INTENT(IN) :: stress_unit
660 REAL(kind=dp), INTENT(IN), OPTIONAL :: pres_diff, pres_tol, pres_diff_constr
661
662 CHARACTER(LEN=5) :: tag
663 INTEGER :: indf
664 LOGICAL :: conv_dx, conv_g, conv_p, conv_rdx, &
665 conv_rg
666 REAL(kind=dp) :: dumm, dxcon, gcon, maxdum(4), rmsgcon, &
667 rmsxcon
668
669 dxcon = gopt_param%max_dr
670 gcon = gopt_param%max_force
671 rmsgcon = gopt_param%rms_force
672 rmsxcon = gopt_param%rms_dr
673
674 conv = .false.
675 conv_dx = .true.
676 conv_rdx = .true.
677 conv_g = .true.
678 conv_rg = .true.
679 conv_p = .true.
680
681 dumm = 0.0_dp
682 DO indf = 1, ndf
683 IF (indf == 1) maxdum(1) = abs(dr(indf))
684 dumm = dumm + dr(indf)**2
685 IF (abs(dr(indf)) > dxcon) conv_dx = .false.
686 IF (abs(dr(indf)) > maxdum(1)) maxdum(1) = abs(dr(indf))
687 END DO
688 ! SQRT(dumm/ndf) > rmsxcon
689 IF (dumm > (rmsxcon*rmsxcon*ndf)) conv_rdx = .false.
690 maxdum(2) = sqrt(dumm/ndf)
691
692 dumm = 0.0_dp
693 DO indf = 1, ndf
694 IF (indf == 1) maxdum(3) = abs(g(indf))
695 dumm = dumm + g(indf)**2
696 IF (abs(g(indf)) > gcon) conv_g = .false.
697 IF (abs(g(indf)) > maxdum(3)) maxdum(3) = abs(g(indf))
698 END DO
699 ! SQRT(dumm/ndf) > rmsgcon
700 IF (dumm > (rmsgcon*rmsgcon*ndf)) conv_rg = .false.
701 maxdum(4) = sqrt(dumm/ndf)
702
703 IF (PRESENT(pres_diff_constr) .AND. PRESENT(pres_tol)) THEN
704 conv_p = abs(pres_diff_constr) < abs(pres_tol)
705 ELSEIF (PRESENT(pres_diff) .AND. PRESENT(pres_tol)) THEN
706 conv_p = abs(pres_diff) < abs(pres_tol)
707 END IF
708
709 IF (output_unit > 0) THEN
710
711 tag = "OPT| "
712
713 WRITE (unit=output_unit, fmt="(T2,A)") trim(tag)
714 WRITE (unit=output_unit, fmt="(T2,A,T55,1X,F25.10)") &
715 tag//"Maximum step size", maxdum(1), &
716 tag//"Convergence limit for maximum step size", dxcon
717 IF (conv_dx) THEN
718 WRITE (unit=output_unit, fmt="(T2,A,T77,A4)") &
719 tag//"Maximum step size is converged", " YES"
720 ELSE
721 WRITE (unit=output_unit, fmt="(T2,A,T77,A4)") &
722 tag//"Maximum step size is converged", " NO"
723 END IF
724
725 WRITE (unit=output_unit, fmt="(T2,A)") trim(tag)
726 WRITE (unit=output_unit, fmt="(T2,A,T55,1X,F25.10)") &
727 tag//"RMS step size", maxdum(2), &
728 tag//"Convergence limit for RMS step size", rmsxcon
729 IF (conv_rdx) THEN
730 WRITE (unit=output_unit, fmt="(T2,A,T77,A4)") &
731 tag//"RMS step size is converged", " YES"
732 ELSE
733 WRITE (unit=output_unit, fmt="(T2,A,T77,A4)") &
734 tag//"RMS step size is converged", " NO"
735 END IF
736
737 WRITE (unit=output_unit, fmt="(T2,A)") trim(tag)
738 WRITE (unit=output_unit, fmt="(T2,A,T55,1X,F25.10)") &
739 tag//"Maximum gradient", maxdum(3), &
740 tag//"Convergence limit for maximum gradient", gcon
741 IF (conv_g) THEN
742 WRITE (unit=output_unit, fmt="(T2,A,T77,A4)") &
743 tag//"Maximum gradient is converged", " YES"
744 ELSE
745 WRITE (unit=output_unit, fmt="(T2,A,T77,A4)") &
746 tag//"Maximum gradient is converged", " NO"
747 END IF
748
749 WRITE (unit=output_unit, fmt="(T2,A)") trim(tag)
750 WRITE (unit=output_unit, fmt="(T2,A,T55,1X,F25.10)") &
751 tag//"RMS gradient", maxdum(4), &
752 tag//"Convergence limit for RMS gradient", rmsgcon
753 IF (conv_rg) THEN
754 WRITE (unit=output_unit, fmt="(T2,A,T77,A4)") &
755 tag//"RMS gradient is converged", " YES"
756 ELSE
757 WRITE (unit=output_unit, fmt="(T2,A,T77,A4)") &
758 tag//"RMS gradient is converged", " NO"
759 END IF
760
761 IF (PRESENT(pres_diff) .AND. PRESENT(pres_tol)) THEN
762 WRITE (unit=output_unit, fmt="(T2,A)") trim(tag)
763 IF (PRESENT(pres_diff_constr)) THEN
764 WRITE (unit=output_unit, fmt="(T2,A,T55,1X,F25.10)") &
765 tag//"Pressure deviation without constraint ["// &
766 trim(adjustl(stress_unit))//"]", &
767 cp_unit_from_cp2k(pres_diff, trim(stress_unit))
768 WRITE (unit=output_unit, fmt="(T2,A,T55,1X,F25.10)") &
769 tag//"Pressure deviation with constraint ["// &
770 trim(adjustl(stress_unit))//"]", &
771 cp_unit_from_cp2k(pres_diff_constr, trim(stress_unit))
772 ELSE
773 WRITE (unit=output_unit, fmt="(T2,A,T55,1X,F25.10)") &
774 tag//"Pressure deviation ["//trim(adjustl(stress_unit))//"]", &
775 cp_unit_from_cp2k(pres_diff, trim(stress_unit))
776 END IF
777 WRITE (unit=output_unit, fmt="(T2,A,T55,1X,F25.10)") &
778 tag//"Pressure tolerance ["//trim(adjustl(stress_unit))//"]", &
779 cp_unit_from_cp2k(pres_tol, trim(stress_unit))
780 IF (conv_p) THEN
781 WRITE (unit=output_unit, fmt="(T2,A,T77,A4)") &
782 tag//"Pressure is converged", " YES"
783 ELSE
784 WRITE (unit=output_unit, fmt="(T2,A,T77,A4)") &
785 tag//"Pressure is converged", " NO"
786 END IF
787 END IF
788
789 WRITE (unit=output_unit, fmt="(T2,A)") tag//repeat("*", 74)
790
791 IF (max_memory /= 0) THEN
792 WRITE (unit=output_unit, fmt="(T2,A,T60,1X,I20)") &
793 tag//"Estimated peak process memory after this step [MiB]", &
794 (max_memory + (1024*1024) - 1)/(1024*1024)
795 END IF
796
797 END IF
798
799 IF (conv_dx .AND. conv_rdx .AND. conv_g .AND. conv_rg .AND. conv_p) conv = .true.
800
801 IF ((conv) .AND. (output_unit > 0)) THEN
802 WRITE (unit=output_unit, fmt="(/,T2,A)") repeat("*", 79)
803 WRITE (unit=output_unit, fmt="(T2,A,T25,A,T78,A)") &
804 "***", "GEOMETRY OPTIMIZATION COMPLETED", "***"
805 WRITE (unit=output_unit, fmt="(T2,A)") repeat("*", 79)
806 END IF
807
808 END SUBROUTINE check_converg
809
810! **************************************************************************************************
811!> \brief ...
812!> \param dimer_env ...
813!> \param output_unit ...
814!> \param conv ...
815!> \date 01.2008
816!> \author Luca Bellucci and Teodoro Laino - created [tlaino]
817! **************************************************************************************************
818 SUBROUTINE check_rot_conv(dimer_env, output_unit, conv)
819
820 TYPE(dimer_env_type), POINTER :: dimer_env
821 INTEGER, INTENT(IN) :: output_unit
822 LOGICAL, INTENT(OUT) :: conv
823
824 CHARACTER(LEN=5) :: tag
825
826 conv = (abs(dimer_env%rot%angle2) < dimer_env%rot%angle_tol)
827
828 IF (output_unit > 0) THEN
829 tag = "OPT| "
830 WRITE (unit=output_unit, fmt="(T2,A)") trim(tag)
831 WRITE (unit=output_unit, fmt="(T2,A,T55,1X,F25.10)") &
832 tag//"Predicted angle step size", dimer_env%rot%angle1, &
833 tag//"Effective angle step size", dimer_env%rot%angle2, &
834 tag//"Convergence limit for angle step size", dimer_env%rot%angle_tol
835 IF (conv) THEN
836 WRITE (unit=output_unit, fmt="(T2,A,T77,A4)") &
837 tag//"Angle step size is converged", " YES"
838 ELSE
839 WRITE (unit=output_unit, fmt="(T2,A,T77,A4)") &
840 tag//"Angle step size is converged", " NO"
841 END IF
842 WRITE (unit=output_unit, fmt="(T2,A)") tag//repeat("*", 74)
843 END IF
844
845 IF ((conv) .AND. (output_unit > 0)) THEN
846 WRITE (unit=output_unit, fmt="(/,T2,A)") repeat("*", 79)
847 WRITE (unit=output_unit, fmt="(T2,A,T25,A,T78,A)") &
848 "***", "ROTATION OPTIMIZATION COMPLETED", "***"
849 WRITE (unit=output_unit, fmt="(T2,A)") repeat("*", 79)
850 END IF
851
852 END SUBROUTINE check_rot_conv
853
854! **************************************************************************************************
855!> \brief ...
856!> \param output_unit ...
857!> \param conv ...
858!> \param it ...
859!> \param gopt_env ...
860!> \param x0 ...
861!> \param master ...
862!> \param para_env ...
863!> \param force_env ...
864!> \param motion_section ...
865!> \param root_section ...
866!> \date 11.2007
867!> \author Teodoro Laino [tlaino] - University of Zurich
868! **************************************************************************************************
869 RECURSIVE SUBROUTINE write_final_info(output_unit, conv, it, gopt_env, x0, master, para_env, force_env, &
870 motion_section, root_section)
871 INTEGER, INTENT(IN) :: output_unit
872 LOGICAL, INTENT(IN) :: conv
873 INTEGER, INTENT(INOUT) :: it
874 TYPE(gopt_f_type), POINTER :: gopt_env
875 REAL(kind=dp), DIMENSION(:), POINTER :: x0
876 INTEGER, INTENT(IN) :: master
877 TYPE(mp_para_env_type), POINTER :: para_env
878 TYPE(force_env_type), POINTER :: force_env
879 TYPE(section_vals_type), POINTER :: motion_section, root_section
880
881 CHARACTER(LEN=4) :: constraint_label
882 LOGICAL :: keep_angles, keep_symmetry, &
883 keep_volume
884 REAL(kind=dp) :: etot
885 TYPE(cell_type), POINTER :: cell
886 TYPE(cp_subsys_type), POINTER :: subsys
887 TYPE(particle_list_type), POINTER :: particles
888 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
889
890 CALL force_env_get(force_env, cell=cell, subsys=subsys)
891 CALL cp_subsys_get(subsys=subsys, particles=particles)
892 particle_set => particles%els
893
894 ! Passing gopt_f_type pointer gopt_env to particle_methods where
895 ! write_final_structure is defined causes a circular dependency, so it
896 ! is necessary to get some flags by preprocessing...
897 keep_angles = .true.
898 keep_symmetry = .true.
899 keep_volume = .true.
900 constraint_label = "NONE"
901 IF (gopt_env%type_id == default_cell_method_id) THEN
902 keep_angles = gopt_env%cell_env%keep_angles
903 keep_symmetry = gopt_env%cell_env%keep_symmetry
904 keep_volume = gopt_env%cell_env%keep_volume
905 SELECT CASE (gopt_env%cell_env%constraint_id)
906 CASE (fix_x)
907 constraint_label = " X"
908 CASE (fix_y)
909 constraint_label = " Y"
910 CASE (fix_z)
911 constraint_label = " Z"
912 CASE (fix_xy)
913 constraint_label = " XY"
914 CASE (fix_xz)
915 constraint_label = " XZ"
916 CASE (fix_yz)
917 constraint_label = " YZ"
918 CASE (fix_none)
919 constraint_label = "NONE"
920 END SELECT
921 END IF
922 CALL write_final_structure(particle_set, cell, motion_section, conv, &
923 keep_angles, keep_symmetry, keep_volume, &
924 gopt_env%label, constraint_label)
925
926 IF (conv) THEN
927 it = it + 1
928 CALL write_structure_data(particle_set, cell, motion_section)
929 CALL write_restart(force_env=force_env, root_section=root_section)
930
931 IF (output_unit > 0) &
932 WRITE (unit=output_unit, fmt="(/,T20,' Reevaluating energy at the minimum')")
933
934 CALL cp_eval_at(gopt_env, x0, f=etot, master=master, final_evaluation=.true., &
935 para_env=para_env)
936 CALL write_geo_traj(force_env, root_section, it, etot)
937 END IF
938
939 END SUBROUTINE write_final_info
940
941! **************************************************************************************************
942!> \brief Specific driver for dumping trajectory during a GEO_OPT
943!> \param force_env ...
944!> \param root_section ...
945!> \param it ...
946!> \param etot ...
947!> \date 11.2007
948!> \par History
949!> 09.2010: Output of core and shell positions and forces (MK)
950!> \author Teodoro Laino [tlaino] - University of Zurich
951! **************************************************************************************************
952 SUBROUTINE write_geo_traj(force_env, root_section, it, etot)
953
954 TYPE(force_env_type), POINTER :: force_env
955 TYPE(section_vals_type), POINTER :: root_section
956 INTEGER, INTENT(IN) :: it
957 REAL(kind=dp), INTENT(IN) :: etot
958
959 LOGICAL :: shell_adiabatic, shell_present
960 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
961 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
962 TYPE(cp_subsys_type), POINTER :: subsys
963 TYPE(particle_list_type), POINTER :: core_particles, shell_particles
964
965 NULLIFY (atomic_kinds)
966 NULLIFY (atomic_kind_set)
967 NULLIFY (core_particles)
968 NULLIFY (shell_particles)
969 NULLIFY (subsys)
970
971 CALL write_trajectory(force_env, root_section, it, 0.0_dp, 0.0_dp, etot)
972 ! Print Force
973 CALL write_trajectory(force_env, root_section, it, 0.0_dp, 0.0_dp, etot, "FORCES", middle_name="frc")
974 CALL force_env_get(force_env, subsys=subsys)
975 CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds)
976 atomic_kind_set => atomic_kinds%els
977 CALL get_atomic_kind_set(atomic_kind_set, &
978 shell_present=shell_present, &
979 shell_adiabatic=shell_adiabatic)
980 IF (shell_present) THEN
981 CALL cp_subsys_get(subsys, &
982 core_particles=core_particles, &
983 shell_particles=shell_particles)
984 CALL write_trajectory(force_env, root_section, it=it, time=0.0_dp, dtime=0.0_dp, &
985 etot=etot, pk_name="SHELL_TRAJECTORY", middle_name="shpos", &
986 particles=shell_particles)
987 IF (shell_adiabatic) THEN
988 CALL write_trajectory(force_env, root_section, it=it, time=0.0_dp, dtime=0.0_dp, &
989 etot=etot, pk_name="SHELL_FORCES", middle_name="shfrc", &
990 particles=shell_particles)
991 CALL write_trajectory(force_env, root_section, it=it, time=0.0_dp, dtime=0.0_dp, &
992 etot=etot, pk_name="CORE_TRAJECTORY", middle_name="copos", &
993 particles=core_particles)
994 CALL write_trajectory(force_env, root_section, it=it, time=0.0_dp, dtime=0.0_dp, &
995 etot=etot, pk_name="CORE_FORCES", middle_name="cofrc", &
996 particles=core_particles)
997 END IF
998 END IF
999
1000 END SUBROUTINE write_geo_traj
1001
1002! **************************************************************************************************
1003!> \brief ...
1004!> \param gopt_env ...
1005!> \param output_unit ...
1006!> \param label ...
1007!> \date 01.2008
1008!> \author Teodoro Laino [tlaino] - University of Zurich
1009! **************************************************************************************************
1010 SUBROUTINE print_geo_opt_header(gopt_env, output_unit, label)
1011
1012 TYPE(gopt_f_type), POINTER :: gopt_env
1013 INTEGER, INTENT(IN) :: output_unit
1014 CHARACTER(LEN=*), INTENT(IN) :: label
1015
1016 CHARACTER(LEN=default_string_length) :: my_format, my_label
1017 INTEGER :: ix
1018
1019 IF (output_unit > 0) THEN
1020 WRITE (unit=output_unit, fmt="(/,T2,A)") repeat("*", 79)
1021 IF (gopt_env%dimer_rotation) THEN
1022 my_label = "OPTIMIZING DIMER ROTATION"
1023 ELSE
1024 my_label = "STARTING "//gopt_env%tag(1:8)//" OPTIMIZATION"
1025 END IF
1026
1027 ix = (80 - 7 - len_trim(my_label))/2
1028 ix = ix + 5
1029 my_format = "(T2,A,T"//cp_to_string(ix)//",A,T78,A)"
1030 WRITE (unit=output_unit, fmt=trim(my_format)) "***", trim(my_label), "***"
1031
1032 ix = (80 - 7 - len_trim(label))/2
1033 ix = ix + 5
1034 my_format = "(T2,A,T"//cp_to_string(ix)//",A,T78,A)"
1035 WRITE (unit=output_unit, fmt=trim(my_format)) "***", trim(label), "***"
1036
1037 WRITE (unit=output_unit, fmt="(T2,A)") repeat("*", 79)
1038 CALL m_flush(output_unit)
1039 END IF
1040 END SUBROUTINE print_geo_opt_header
1041
1042! **************************************************************************************************
1043!> \brief ...
1044!> \param gopt_env ...
1045!> \param output_unit ...
1046!> \date 01.2008
1047!> \author Teodoro Laino [tlaino] - University of Zurich
1048! **************************************************************************************************
1049 SUBROUTINE print_geo_opt_nc(gopt_env, output_unit)
1050
1051 TYPE(gopt_f_type), POINTER :: gopt_env
1052 INTEGER, INTENT(IN) :: output_unit
1053
1054 IF (output_unit > 0) THEN
1055 WRITE (unit=output_unit, fmt="(/,T2,A)") &
1056 "*** MAXIMUM NUMBER OF OPTIMIZATION STEPS REACHED ***"
1057 IF (.NOT. gopt_env%dimer_rotation) THEN
1058 WRITE (unit=output_unit, fmt="(T2,A)") &
1059 "*** EXITING GEOMETRY OPTIMIZATION ***"
1060 ELSE
1061 WRITE (unit=output_unit, fmt="(T2,A)") &
1062 "*** EXITING ROTATION OPTIMIZATION ***"
1063 END IF
1064 CALL m_flush(output_unit)
1065 END IF
1066
1067 END SUBROUTINE print_geo_opt_nc
1068
1069! **************************************************************************************************
1070!> \brief Prints information during GEO_OPT common to all optimizers
1071!> \param force_env ...
1072!> \param root_section ...
1073!> \param motion_section ...
1074!> \param its ...
1075!> \param opt_energy ...
1076!> \date 02.2008
1077!> \author Teodoro Laino [tlaino] - University of Zurich
1078!> \version 1.0
1079! **************************************************************************************************
1080 SUBROUTINE geo_opt_io(force_env, root_section, motion_section, its, opt_energy)
1081
1082 TYPE(force_env_type), POINTER :: force_env
1083 TYPE(section_vals_type), POINTER :: root_section, motion_section
1084 INTEGER, INTENT(IN) :: its
1085 REAL(kind=dp), INTENT(IN) :: opt_energy
1086
1087 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
1088 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1089 TYPE(cell_type), POINTER :: cell
1090 TYPE(cp_subsys_type), POINTER :: subsys
1091 TYPE(distribution_1d_type), POINTER :: local_particles
1092 TYPE(mp_para_env_type), POINTER :: para_env
1093 TYPE(particle_list_type), POINTER :: particles
1094 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1095 TYPE(virial_type), POINTER :: virial
1096
1097 NULLIFY (para_env, atomic_kind_set, subsys, particle_set, &
1098 local_particles, atomic_kinds, particles)
1099
1100 ! Write Restart File
1101 CALL write_restart(force_env=force_env, root_section=root_section)
1102
1103 ! Write Trajectory
1104 CALL write_geo_traj(force_env, root_section, its, opt_energy)
1105
1106 ! Write the stress Tensor
1107 CALL force_env_get(force_env, cell=cell, para_env=para_env, &
1108 subsys=subsys)
1109 CALL cp_subsys_get(subsys=subsys, atomic_kinds=atomic_kinds, local_particles=local_particles, &
1110 particles=particles, virial=virial)
1111 atomic_kind_set => atomic_kinds%els
1112 particle_set => particles%els
1113 CALL virial_evaluate(atomic_kind_set, particle_set, local_particles, &
1114 virial, para_env)
1115 CALL write_stress_tensor_to_file(virial, cell, motion_section, its, 0.0_dp)
1116
1117 ! Write the cell
1118 CALL write_simulation_cell(cell, motion_section, its, 0.0_dp)
1119
1120 END SUBROUTINE geo_opt_io
1121
1122! **************************************************************************************************
1123!> \brief Apply coordinate transformations after cell (shape) change
1124!> \param gopt_env ...
1125!> \param cell ...
1126!> \param x ...
1127!> \param update_forces ...
1128!> \date 05.11.2012 (revised version of unbiase_coordinates moved here, MK)
1129!> \author Matthias Krack
1130!> \version 1.0
1131! **************************************************************************************************
1132 SUBROUTINE apply_cell_change(gopt_env, cell, x, update_forces)
1133
1134 TYPE(gopt_f_type), POINTER :: gopt_env
1135 TYPE(cell_type), POINTER :: cell
1136 REAL(kind=dp), DIMENSION(:), POINTER :: x
1137 LOGICAL, INTENT(IN) :: update_forces
1138
1139 INTEGER :: i, iatom, idg, j, natom, nparticle, &
1140 shell_index
1141 REAL(kind=dp) :: fc, fs, mass
1142 REAL(kind=dp), DIMENSION(3) :: s
1143 TYPE(cell_type), POINTER :: cell_ref
1144 TYPE(cp_subsys_type), POINTER :: subsys
1145 TYPE(particle_list_type), POINTER :: core_particles, particles, &
1146 shell_particles
1147
1148 NULLIFY (cell_ref)
1149 NULLIFY (core_particles)
1150 NULLIFY (particles)
1151 NULLIFY (shell_particles)
1152 NULLIFY (subsys)
1153
1154 natom = force_env_get_natom(gopt_env%force_env)
1155 nparticle = force_env_get_nparticle(gopt_env%force_env)
1156 CALL force_env_get(gopt_env%force_env, &
1157 subsys=subsys)
1158 CALL cp_subsys_get(subsys=subsys, &
1159 core_particles=core_particles, &
1160 particles=particles, &
1161 shell_particles=shell_particles)
1162
1163 ! Retrieve the reference cell
1164 CALL cell_create(cell_ref)
1165 CALL cell_copy(cell, cell_ref, tag="CELL_OPT_REF")
1166
1167 ! Load the updated cell information
1168 idg = 3*nparticle
1169 CALL init_cell(cell_ref, hmat=gopt_env%h_ref)
1170 cpassert((SIZE(x) == idg + 6))
1171
1172 IF (update_forces) THEN
1173
1174 ! Transform particle forces back to reference cell
1175 idg = 1
1176 DO iatom = 1, natom
1177 CALL real_to_scaled(s, x(idg:idg + 2), cell)
1178 CALL scaled_to_real(x(idg:idg + 2), s, cell_ref)
1179 idg = idg + 3
1180 END DO
1181
1182 ELSE
1183
1184 ! Update cell
1185 DO i = 1, 3
1186 DO j = 1, i
1187 idg = idg + 1
1188 cell%hmat(j, i) = x(idg)
1189 END DO
1190 END DO
1191 CALL init_cell(cell)
1192 CALL cp_subsys_set(subsys, cell=cell)
1193
1194 ! Retrieve particle coordinates for the current cell
1195 idg = 1
1196 DO iatom = 1, natom
1197 CALL real_to_scaled(s, x(idg:idg + 2), cell_ref)
1198 shell_index = particles%els(iatom)%shell_index
1199 IF (shell_index == 0) THEN
1200 CALL scaled_to_real(particles%els(iatom)%r, s, cell)
1201 ELSE
1202 CALL scaled_to_real(core_particles%els(shell_index)%r, s, cell)
1203 i = 3*(natom + shell_index - 1) + 1
1204 CALL real_to_scaled(s, x(i:i + 2), cell_ref)
1205 CALL scaled_to_real(shell_particles%els(shell_index)%r, s, cell)
1206 ! Update atomic position due to core and shell motion
1207 mass = particles%els(iatom)%atomic_kind%mass
1208 fc = core_particles%els(shell_index)%atomic_kind%shell%mass_core/mass
1209 fs = shell_particles%els(shell_index)%atomic_kind%shell%mass_shell/mass
1210 particles%els(iatom)%r(1:3) = fc*core_particles%els(shell_index)%r(1:3) + &
1211 fs*shell_particles%els(shell_index)%r(1:3)
1212 END IF
1213 idg = idg + 3
1214 END DO
1215 END IF
1216
1217 CALL cell_release(cell_ref)
1218
1219 END SUBROUTINE apply_cell_change
1220
1221END MODULE gopt_f_methods
cp_eval_at
subroutine cp_eval_at(gopt_env, x, f, gradient, master, final_evaluation, para_env)
evaluete the potential energy and its gradients using an array with same dimension as the particle_se...
Definition gopt_f77_methods.F:25
cp_log_handling::cp_to_string
Definition cp_log_handling.F:90
atomic_kind_list_types
represent a simple array based list of the given type
Definition atomic_kind_list_types.F:15
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind_set
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Definition atomic_kind_types.F:223
cell_methods
Handles all functions related to the CELL.
Definition cell_methods.F:15
cell_methods::init_cell
subroutine, public init_cell(cell, hmat, periodic)
Initialise/readjust a simulation cell after hmat has been changed.
Definition cell_methods.F:182
cell_methods::cell_create
subroutine, public cell_create(cell, hmat, periodic, tag)
allocates and initializes a cell
Definition cell_methods.F:92
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cell_types::scaled_to_real
subroutine, public scaled_to_real(r, s, cell)
Transform scaled cell coordinates real coordinates. r=h*s.
Definition cell_types.F:565
cell_types::real_to_scaled
subroutine, public real_to_scaled(s, r, cell)
Transform real to scaled cell coordinates. s=h_inv*r.
Definition cell_types.F:535
cell_types::cell_release
subroutine, public cell_release(cell)
releases the given cell (see doc/ReferenceCounting.html)
Definition cell_types.F:608
cell_types::cell_copy
subroutine, public cell_copy(cell_in, cell_out, tag)
Copy cell variable.
Definition cell_types.F:137
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_subsys_types
types that represent a subsys, i.e. a part of the system
Definition cp_subsys_types.F:16
cp_subsys_types::cp_subsys_set
subroutine, public cp_subsys_set(subsys, atomic_kinds, particles, local_particles, molecules, molecule_kinds, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, results, cell, cell_ref, use_ref_cell)
sets various propreties of the subsys
Definition cp_subsys_types.F:210
cp_subsys_types::cp_subsys_get
subroutine, public cp_subsys_get(subsys, ref_count, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell, cell_ref, use_ref_cell)
returns information about various attributes of the given subsys
Definition cp_subsys_types.F:364
cp_subsys_types::pack_subsys_particles
subroutine, public pack_subsys_particles(subsys, f, r, s, v, fscale, cell)
Pack components of a subsystem particle sets into a single vector.
Definition cp_subsys_types.F:467
cp_units
unit conversion facility
Definition cp_units.F:30
cp_units::cp_unit_from_cp2k
real(kind=dp) function, public cp_unit_from_cp2k(value, unit_str, defaults, power)
converts from the internal cp2k units to the given unit
Definition cp_units.F:1239
dimer_types
Contains types used for a Dimer Method calculations.
Definition dimer_types.F:14
dimer_utils
Contains utilities for a Dimer Method calculations.
Definition dimer_utils.F:14
dimer_utils::update_dimer_vec
subroutine, public update_dimer_vec(dimer_env, motion_section)
Updates the orientation of the dimer vector in the input file.
Definition dimer_utils.F:74
distribution_1d_types
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
Definition distribution_1d_types.F:24
force_env_types
Interface for the force calculations.
Definition force_env_types.F:18
force_env_types::force_env_get_natom
integer function, public force_env_get_natom(force_env)
returns the number of atoms
Definition force_env_types.F:537
force_env_types::use_qmmm
integer, parameter, public use_qmmm
Definition force_env_types.F:83
force_env_types::force_env_get
recursive subroutine, public force_env_get(force_env, in_use, fist_env, qs_env, meta_env, fp_env, subsys, para_env, potential_energy, additional_potential, kinetic_energy, harmonic_shell, kinetic_shell, cell, sub_force_env, qmmm_env, qmmmx_env, eip_env, pwdft_env, globenv, input, force_env_section, method_name_id, root_section, mixed_env, nnp_env, embed_env, ipi_env)
returns various attributes about the force environment
Definition force_env_types.F:342
force_env_types::use_qmmmx
integer, parameter, public use_qmmmx
Definition force_env_types.F:83
force_env_types::force_env_get_nparticle
integer function, public force_env_get_nparticle(force_env)
returns the number of particles in a force environment
Definition force_env_types.F:558
gopt_f_methods
contains a functional that calculates the energy and its derivatives for the geometry optimizer
Definition gopt_f_methods.F:14
gopt_f_methods::print_geo_opt_header
subroutine, public print_geo_opt_header(gopt_env, output_unit, label)
...
Definition gopt_f_methods.F:1011
gopt_f_methods::gopt_f_io_init
subroutine, public gopt_f_io_init(gopt_env, output_unit, opt_energy, wildcard, its, used_time)
Handles the Output during an optimization run.
Definition gopt_f_methods.F:209
gopt_f_methods::gopt_f_create_x0
subroutine, public gopt_f_create_x0(gopt_env, x0)
returns the value of the parameters for the actual configuration
Definition gopt_f_methods.F:132
gopt_f_methods::gopt_f_io_finalize
recursive subroutine, public gopt_f_io_finalize(gopt_env, force_env, x0, conv, its, root_section, para_env, master, output_unit)
Handles the Output at the end of an optimization run.
Definition gopt_f_methods.F:464
gopt_f_methods::apply_cell_change
subroutine, public apply_cell_change(gopt_env, cell, x, update_forces)
Apply coordinate transformations after cell (shape) change.
Definition gopt_f_methods.F:1133
gopt_f_methods::gopt_f_io
subroutine, public gopt_f_io(gopt_env, force_env, root_section, its, opt_energy, output_unit, eold, emin, wildcard, gopt_param, ndf, dx, xi, conv, pred, rat, step, rad, used_time)
Handles the Output during an optimization run.
Definition gopt_f_methods.F:310
gopt_f_methods::print_geo_opt_nc
subroutine, public print_geo_opt_nc(gopt_env, output_unit)
...
Definition gopt_f_methods.F:1050
gopt_f_methods::gopt_f_ii
subroutine, public gopt_f_ii(its, output_unit)
Prints iteration step of the optimization procedure on screen.
Definition gopt_f_methods.F:186
gopt_f_types
contains a functional that calculates the energy and its derivatives for the geometry optimizer
Definition gopt_f_types.F:15
gopt_param_types
contains typo and related routines to handle parameters controlling the GEO_OPT module
Definition gopt_param_types.F:14
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::default_shellcore_method_id
integer, parameter, public default_shellcore_method_id
Definition input_constants.F:549
input_constants::fix_xz
integer, parameter, public fix_xz
Definition input_constants.F:1291
input_constants::default_cell_method_id
integer, parameter, public default_cell_method_id
Definition input_constants.F:549
input_constants::default_minimization_method_id
integer, parameter, public default_minimization_method_id
Definition input_constants.F:549
input_constants::default_ts_method_id
integer, parameter, public default_ts_method_id
Definition input_constants.F:549
input_constants::fix_y
integer, parameter, public fix_y
Definition input_constants.F:1291
input_constants::fix_none
integer, parameter, public fix_none
Definition input_constants.F:1291
input_constants::fix_z
integer, parameter, public fix_z
Definition input_constants.F:1291
input_constants::fix_xy
integer, parameter, public fix_xy
Definition input_constants.F:1291
input_constants::fix_yz
integer, parameter, public fix_yz
Definition input_constants.F:1291
input_constants::fix_x
integer, parameter, public fix_x
Definition input_constants.F:1291
input_cp2k_restarts
Set of routines to dump the restart file of CP2K.
Definition input_cp2k_restarts.F:13
input_cp2k_restarts::write_restart
subroutine, public write_restart(md_env, force_env, root_section, coords, vels, pint_env, helium_env)
checks if a restart needs to be written and does so, updating all necessary fields in the input file....
Definition input_cp2k_restarts.F:112
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1047
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::int_8
integer, parameter, public int_8
Definition kinds.F:54
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
machine
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
machine::m_flush
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition machine.F:124
md_energies
prints all energy info per timestep to the screen or to user defined output files
Definition md_energies.F:16
md_energies::sample_memory
integer(kind=int_8) function, public sample_memory(para_env)
Samples memory usage.
Definition md_energies.F:812
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
motion_utils
Output Utilities for MOTION_SECTION.
Definition motion_utils.F:13
motion_utils::write_simulation_cell
subroutine, public write_simulation_cell(cell, motion_section, itimes, time, pos, act)
Prints the Simulation Cell.
Definition motion_utils.F:621
motion_utils::write_trajectory
subroutine, public write_trajectory(force_env, root_section, it, time, dtime, etot, pk_name, pos, act, middle_name, particles, extended_xmol_title)
Prints the information controlled by the TRAJECTORY section.
Definition motion_utils.F:282
motion_utils::write_stress_tensor_to_file
subroutine, public write_stress_tensor_to_file(virial, cell, motion_section, itimes, time, pos, act)
Prints the Stress Tensor.
Definition motion_utils.F:549
particle_list_types
represent a simple array based list of the given type
Definition particle_list_types.F:15
particle_methods
Define methods related to particle_type.
Definition particle_methods.F:14
particle_methods::write_final_structure
subroutine, public write_final_structure(particle_set, cell, input_section, conv, keep_angles, keep_symmetry, keep_volume, gopt_env_label, constraint_label)
Write the final geometry and cell information to files.
Definition particle_methods.F:1083
particle_methods::write_structure_data
subroutine, public write_structure_data(particle_set, cell, input_section)
Write structure data requested by a separate structure data input section to the output unit....
Definition particle_methods.F:829
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
qmmm_util
Definition qmmm_util.F:13
qmmm_util::apply_qmmm_translate
subroutine, public apply_qmmm_translate(qmmm_env)
Apply translation to the full system in order to center the QM system into the QM box.
Definition qmmm_util.F:374
qmmmx_util
Routines used for force-mixing QM/MM calculations.
Definition qmmmx_util.F:14
qmmmx_util::apply_qmmmx_translate
subroutine, public apply_qmmmx_translate(qmmmx_env)
Apply translation to the full system in order to center the QM system into the QM box.
Definition qmmmx_util.F:75
virial_methods
Definition virial_methods.F:12
virial_methods::virial_evaluate
subroutine, public virial_evaluate(atomic_kind_set, particle_set, local_particles, virial, igroup)
Computes the kinetic part of the pressure tensor and updates the full VIRIAL (PV)
Definition virial_methods.F:83
virial_types
Definition virial_types.F:13
atomic_kind_list_types::atomic_kind_list_type
represent a list of objects
Definition atomic_kind_list_types.F:45
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
cp_subsys_types::cp_subsys_type
represents a system: atoms, molecules, their pos,vel,...
Definition cp_subsys_types.F:89
dimer_types::dimer_env_type
Defines the environment for a Dimer Method calculation.
Definition dimer_types.F:103
distribution_1d_types::distribution_1d_type
structure to store local (to a processor) ordered lists of integers.
Definition distribution_1d_types.F:62
force_env_types::force_env_type
wrapper to abstract the force evaluation of the various methods
Definition force_env_types.F:145
gopt_f_types::gopt_f_type
calculates the potential energy of a system, and its derivatives
Definition gopt_f_types.F:60
gopt_param_types::gopt_param_type
Definition gopt_param_types.F:55
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:743
particle_list_types::particle_list_type
represent a list of objects
Definition particle_list_types.F:46
particle_types::particle_type
Definition particle_types.F:35
virial_types::virial_type
Definition virial_types.F:26