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