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