(git:6d276e9)
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! **************************************************************************************************
15
19 USE bibliography, ONLY: henkelman1999,&
20 cite_reference
21 USE cell_methods, ONLY: cell_create,&
22 init_cell,&
24 USE cell_opt_utils, ONLY: get_dg_dh,&
26 USE cell_types, ONLY: cell_copy,&
28 cell_type,&
47 use_qmmm,&
49 USE gopt_f_types, ONLY: gopt_f_type
51 USE input_constants, ONLY: &
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
79 USE virial_types, ONLY: virial_type
80#include "../base/base_uses.f90"
81
82 IMPLICIT NONE
83 PRIVATE
84
85 LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .true.
86 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gopt_f_methods'
87
88 PUBLIC :: cp_eval_at, &
93
94CONTAINS
95
96! **************************************************************************************************
97!> \brief returns the value of the parameters for the actual configuration
98!> \param gopt_env the geometry optimization environment you want the info about
99!> x0: the parameter vector (is allocated by this routine)
100!> \param x0 ...
101!> \par History
102!> - Cell optimization revised (06.11.2012,MK)
103! **************************************************************************************************
104 SUBROUTINE gopt_f_create_x0(gopt_env, x0)
105
106 TYPE(gopt_f_type), POINTER :: gopt_env
107 REAL(kind=dp), DIMENSION(:), POINTER :: x0
108
109 INTEGER :: i, idg, j, nparticle
110 TYPE(cell_type), POINTER :: cell
111 TYPE(cp_subsys_type), POINTER :: subsys
112
113 NULLIFY (cell)
114 NULLIFY (subsys)
115
116 SELECT CASE (gopt_env%type_id)
118 CALL force_env_get(gopt_env%force_env, subsys=subsys)
119 ! before starting we handle the case of translating coordinates (QM/MM)
120 IF (gopt_env%force_env%in_use == use_qmmm) THEN
121 CALL apply_qmmm_translate(gopt_env%force_env%qmmm_env)
122 END IF
123 IF (gopt_env%force_env%in_use == use_qmmmx) THEN
124 CALL apply_qmmmx_translate(gopt_env%force_env%qmmmx_env)
125 END IF
126 nparticle = force_env_get_nparticle(gopt_env%force_env)
127 ALLOCATE (x0(3*nparticle))
128 CALL pack_subsys_particles(subsys=subsys, r=x0)
130 CALL force_env_get(gopt_env%force_env, subsys=subsys, cell=cell)
131 ! Store reference cell
132 gopt_env%h_ref = cell%hmat
133 ! before starting we handle the case of translating coordinates (QM/MM)
134 IF (gopt_env%force_env%in_use == use_qmmm) THEN
135 CALL apply_qmmm_translate(gopt_env%force_env%qmmm_env)
136 END IF
137 IF (gopt_env%force_env%in_use == use_qmmmx) THEN
138 CALL apply_qmmmx_translate(gopt_env%force_env%qmmmx_env)
139 END IF
140 nparticle = force_env_get_nparticle(gopt_env%force_env)
141 ALLOCATE (x0(3*nparticle + 6))
142 CALL pack_subsys_particles(subsys=subsys, r=x0)
143 idg = 3*nparticle
144 DO i = 1, 3
145 DO j = 1, i
146 idg = idg + 1
147 x0(idg) = gopt_env%cell_env%opt_cell%hmat(j, i)
148 END DO
149 END DO
150 CASE DEFAULT
151 cpabort("Invalid or not yet implemented type of optimization")
152 END SELECT
153
154 END SUBROUTINE gopt_f_create_x0
155
156! **************************************************************************************************
157!> \brief evaluete the potential energy and its gradients using an array
158!> with same dimension as the particle_set
159!> \param gopt_env the geometry optimization environment
160!> \param x the position where the function should be evaluated
161!> \param f the function value
162!> \param gradient the value of its gradient
163!> \param master ...
164!> \param final_evaluation ...
165!> \param para_env ...
166!> \par History
167!> CELL OPTIMIZATION: Teodoro Laino [tlaino] - University of Zurich - 03.2008
168!> 07.2020 Pierre Cazade [pcazade] Space Group Symmetry
169!> \author Teodoro Laino [tlaino] - University of Zurich - 01.2008
170! **************************************************************************************************
171 SUBROUTINE cp_eval_at(gopt_env, x, f, gradient, master, &
172 final_evaluation, para_env)
173
174 TYPE(gopt_f_type), POINTER :: gopt_env
175 REAL(kind=dp), DIMENSION(:), POINTER :: x
176 REAL(kind=dp), INTENT(OUT), OPTIONAL :: f
177 REAL(kind=dp), DIMENSION(:), OPTIONAL, POINTER :: gradient
178 INTEGER, INTENT(IN) :: master
179 LOGICAL, INTENT(IN), OPTIONAL :: final_evaluation
180 TYPE(mp_para_env_type), POINTER :: para_env
181
182 CHARACTER(len=*), PARAMETER :: routinen = 'cp_eval_at'
183
184 INTEGER :: handle, idg, idir, ip, nparticle, nsize, &
185 shell_index
186 REAL(kind=dp) :: f_ts
187 REAL(kind=dp), DIMENSION(3, 3) :: av_ptens, av_ptens_opt
188 REAL(kind=dp), DIMENSION(:), POINTER :: cell_gradient, gradient_ts
189 TYPE(cell_type), POINTER :: cell
190 TYPE(cp_subsys_type), POINTER :: subsys
191 TYPE(particle_list_type), POINTER :: core_particles, particles, &
192 shell_particles
193 TYPE(spgr_type), POINTER :: spgr
194 TYPE(virial_type), POINTER :: virial
195
196 mark_used(final_evaluation)
197
198 NULLIFY (cell)
199 NULLIFY (core_particles)
200 NULLIFY (gradient_ts)
201 NULLIFY (particles)
202 NULLIFY (shell_particles)
203 NULLIFY (subsys)
204 NULLIFY (virial)
205 NULLIFY (spgr)
206
207 CALL timeset(routinen, handle)
208
209 CALL force_env_get(gopt_env%force_env, subsys=subsys, cell=cell)
210 CALL cp_subsys_get(subsys, &
211 core_particles=core_particles, &
212 particles=particles, &
213 shell_particles=shell_particles, &
214 virial=virial)
215
216 spgr => gopt_env%spgr
217
218 SELECT CASE (gopt_env%type_id)
220 CALL unpack_subsys_particles(subsys=subsys, r=x)
221 CALL write_structure_data(particles%els, cell, gopt_env%motion_section)
222 SELECT CASE (gopt_env%type_id)
224 ! Geometry Minimization
225 CALL force_env_calc_energy_force(gopt_env%force_env, &
226 calc_force=PRESENT(gradient), &
227 require_consistent_energy_force=gopt_env%require_consistent_energy_force)
228 ! Possibly take the potential energy
229 IF (PRESENT(f)) THEN
230 CALL force_env_get(gopt_env%force_env, potential_energy=f)
231 END IF
232 ! Possibly take the gradients
233 IF (PRESENT(gradient)) THEN
234 IF (master == para_env%mepos) THEN ! we are on the master
235 CALL pack_subsys_particles(subsys=subsys, f=gradient, fscale=-1.0_dp)
236 IF (spgr%keep_space_group) THEN
237 CALL spgr_apply_rotations_force(spgr, gradient)
238 CALL unpack_subsys_particles(subsys=subsys, f=gradient, fscale=-1.0_dp)
239 END IF
240 END IF
241 END IF
243 ! Transition State Optimization
244 ALLOCATE (gradient_ts(particles%n_els*3))
245 ! Real calculation of energy and forces for transition state optimization:
246 ! When doing dimer methods forces have to be always computed since the function
247 ! to minimize is not the energy but the effective force
248 CALL cp_eval_at_ts(gopt_env, x, f_ts, gradient_ts, calc_force=.true.)
249 CALL cite_reference(henkelman1999)
250 ! Possibly take the potential energy
251 IF (PRESENT(f)) f = f_ts
252 ! Possibly take the gradients
253 IF (PRESENT(gradient)) THEN
254 IF (master == para_env%mepos) THEN ! we are on the master
255 cpassert(ASSOCIATED(gradient))
256 gradient = gradient_ts
257 END IF
258 END IF
259 DEALLOCATE (gradient_ts)
260 END SELECT
261 ! This call is necessary for QM/MM if a Translation is applied
262 ! this makes the geometry optimizer consistent
263 CALL unpack_subsys_particles(subsys=subsys, r=x)
265 ! Check for VIRIAL
266 IF (.NOT. virial%pv_availability) THEN
267 CALL cp_abort(__location__, &
268 "For the CELL_OPT task, the FORCE_EVAL/STRESS_TENSOR "// &
269 "keyword MUST be defined in the input file for the "// &
270 "evaluation of the stress tensor, but none is found!")
271 END IF
272 IF (gopt_env%cell_env%keep_volume) THEN
273 nparticle = force_env_get_nparticle(gopt_env%force_env)
274 idg = 3*nparticle
275 CALL rescale_new_cell_volume(cell%deth, x, idg)
276 END IF
277
278 CALL apply_cell_change(gopt_env, cell, x, update_forces=.false.)
279 ! Possibly output the new cell used for the next calculation
280 CALL write_cell(cell, gopt_env%geo_section)
281 ! Compute the pressure tensor
282 block
283 TYPE(virial_type) :: virial_avg
284 CALL force_env_calc_energy_force(gopt_env%force_env, &
285 calc_force=PRESENT(gradient), &
286 require_consistent_energy_force=gopt_env%require_consistent_energy_force)
287 ! Possibly take the potential energy
288 virial_avg = virial
289 CALL virial_update(virial_avg, subsys, para_env)
290 IF (PRESENT(f)) THEN
291 CALL force_env_get(gopt_env%force_env, potential_energy=f)
292 END IF
293 ! Possibly take the gradients
294 IF (PRESENT(gradient)) THEN
295 cpassert(any(virial_avg%pv_total /= 0))
296 ! Convert the average ptens
297 av_ptens(:, :) = virial_avg%pv_total(:, :)/cell%deth
298 IF (master == para_env%mepos) THEN ! we are on the master
299 cpassert(ASSOCIATED(gradient))
300 nparticle = force_env_get_nparticle(gopt_env%force_env)
301 nsize = 3*nparticle
302 cpassert((SIZE(gradient) == nsize + 6))
303 CALL pack_subsys_particles(subsys=subsys, f=gradient(1:nsize), fscale=-1.0_dp)
304 CALL apply_cell_change(gopt_env, cell, gradient, update_forces=.true.)
305 IF (spgr%keep_space_group) THEN
306 CALL spgr_apply_rotations_force(spgr, gradient)
307 CALL spgr_apply_rotations_stress(spgr, cell, av_ptens)
308 CALL spgr_write_stress_tensor(av_ptens, spgr)
309 END IF
310 cell_gradient => gradient(nsize + 1:nsize + 6)
311 cell_gradient = 0.0_dp
312 av_ptens_opt = matmul(gopt_env%cell_env%input_to_opt, &
313 matmul(av_ptens, gopt_env%cell_env%opt_to_input))
314 CALL get_dg_dh(cell_gradient, av_ptens_opt, gopt_env%cell_env%pres_ext, &
315 gopt_env%cell_env%opt_cell, gopt_env%cell_env%mtrx, &
316 keep_angles=gopt_env%cell_env%keep_angles, &
317 keep_symmetry=gopt_env%cell_env%keep_symmetry, &
318 pres_int=gopt_env%cell_env%pres_int, &
319 pres_constr=gopt_env%cell_env%pres_constr, &
320 constraint_id=gopt_env%cell_env%constraint_id)
321 END IF
322 ! some callers expect pres_int to be available on all ranks. Also, here master is not necessarily a single rank.
323 ! Assume at least master==0
324 CALL para_env%bcast(gopt_env%cell_env%pres_int, 0)
325 IF (gopt_env%cell_env%constraint_id /= fix_none) THEN
326 CALL para_env%bcast(gopt_env%cell_env%pres_constr, 0)
327 END IF
328 END IF
329 END block
331 idg = 0
332 DO ip = 1, particles%n_els
333 shell_index = particles%els(ip)%shell_index
334 IF (shell_index /= 0) THEN
335 DO idir = 1, 3
336 idg = 3*(shell_index - 1) + idir
337 shell_particles%els(shell_index)%r(idir) = core_particles%els(ip)%r(idir) - x(idg)
338 END DO
339 END IF
340 END DO
341 CALL write_structure_data(particles%els, cell, gopt_env%motion_section)
342
343 ! Shell-core optimization
344 CALL force_env_calc_energy_force(gopt_env%force_env, &
345 calc_force=PRESENT(gradient), &
346 require_consistent_energy_force=gopt_env%require_consistent_energy_force)
347
348 ! Possibly take the potential energy
349 IF (PRESENT(f)) THEN
350 CALL force_env_get(gopt_env%force_env, potential_energy=f)
351 END IF
352
353 ! Possibly take the gradients
354 IF (PRESENT(gradient)) THEN
355 IF (master == para_env%mepos) THEN ! we are on the master
356 cpassert(ASSOCIATED(gradient))
357 idg = 0
358 DO ip = 1, shell_particles%n_els
359 DO idir = 1, 3
360 idg = idg + 1
361 gradient(idg) = -(core_particles%els(ip)%f(idir) - shell_particles%els(ip)%f(idir))
362 END DO
363 END DO
364 END IF
365 END IF
366 CASE DEFAULT
367 cpabort("Invalid or not yet implemented type of optimization")
368 END SELECT
369
370 CALL timestop(handle)
371
372 END SUBROUTINE cp_eval_at
373
374! **************************************************************************************************
375!> \brief Prints iteration step of the optimization procedure on screen
376!> \param its ...
377!> \param output_unit ...
378!> \author Teodoro Laino [tlaino] - University of Zurich - 03.2008
379! **************************************************************************************************
380 SUBROUTINE gopt_f_ii(its, output_unit)
381
382 INTEGER, INTENT(IN) :: its, output_unit
383
384 IF (output_unit > 0) THEN
385 WRITE (unit=output_unit, fmt="(/,T2,26('-'))")
386 WRITE (unit=output_unit, fmt="(T2,A,I6)") "OPTIMIZATION STEP: ", its
387 WRITE (unit=output_unit, fmt="(T2,26('-'))")
388 CALL m_flush(output_unit)
389 END IF
390
391 END SUBROUTINE gopt_f_ii
392
393! **************************************************************************************************
394!> \brief Handles the Output during an optimization run
395!> \param gopt_env ...
396!> \param output_unit ...
397!> \param opt_energy ...
398!> \param wildcard ...
399!> \param its ...
400!> \param used_time ...
401!> \author Teodoro Laino [tlaino] - University of Zurich - 03.2008
402! **************************************************************************************************
403 SUBROUTINE gopt_f_io_init(gopt_env, output_unit, opt_energy, wildcard, its, used_time)
404
405 TYPE(gopt_f_type), POINTER :: gopt_env
406 INTEGER, INTENT(IN) :: output_unit
407 REAL(kind=dp) :: opt_energy
408 CHARACTER(LEN=5) :: wildcard
409 INTEGER, INTENT(IN) :: its
410 REAL(kind=dp) :: used_time
411
412 CHARACTER(LEN=default_string_length) :: energy_unit, stress_unit
413 INTEGER(KIND=int_8) :: max_memory
414 LOGICAL :: print_memory
415 REAL(kind=dp) :: pres_int
416 TYPE(mp_para_env_type), POINTER :: para_env
417
418 NULLIFY (para_env)
419 CALL section_vals_val_get(gopt_env%motion_section, "PRINT%MEMORY_INFO", l_val=print_memory)
420 max_memory = 0
421 IF (print_memory) THEN
422 CALL force_env_get(gopt_env%force_env, para_env=para_env)
423 max_memory = sample_memory(para_env)
424 END IF
425
426 CALL section_vals_val_get(gopt_env%force_env%force_env_section, &
427 "PRINT%PROGRAM_RUN_INFO%ENERGY_UNIT", &
428 c_val=energy_unit)
429 CALL section_vals_val_get(gopt_env%force_env%force_env_section, &
430 "PRINT%STRESS_TENSOR%STRESS_UNIT", &
431 c_val=stress_unit)
432
433 SELECT CASE (gopt_env%type_id)
434 CASE (default_ts_method_id, default_minimization_method_id)
435 ! Geometry Optimization (Minimization and Transition State Search)
436 IF (.NOT. gopt_env%dimer_rotation) THEN
437 CALL write_cycle_infos(output_unit, &
438 it=its, &
439 etot=opt_energy, &
440 wildcard=wildcard, &
441 used_time=used_time, &
442 max_memory=max_memory, &
443 energy_unit=energy_unit, &
444 stress_unit=stress_unit)
445 ELSE
446 CALL write_rot_cycle_infos(output_unit, &
447 it=its, &
448 etot=opt_energy, &
449 dimer_env=gopt_env%dimer_env, &
450 wildcard=wildcard, &
451 used_time=used_time, &
452 max_memory=max_memory)
453 END IF
454 CASE (default_cell_method_id)
455 ! Cell Optimization
456 pres_int = gopt_env%cell_env%pres_int
457 CALL write_cycle_infos(output_unit, &
458 it=its, &
459 etot=opt_energy, &
460 pres_int=pres_int, &
461 wildcard=wildcard, &
462 used_time=used_time, &
463 max_memory=max_memory, &
464 energy_unit=energy_unit, &
465 stress_unit=stress_unit)
466 CASE (default_shellcore_method_id)
467 CALL write_cycle_infos(output_unit, &
468 it=its, &
469 etot=opt_energy, &
470 wildcard=wildcard, &
471 used_time=used_time, &
472 max_memory=max_memory, &
473 energy_unit=energy_unit, &
474 stress_unit=stress_unit)
475 END SELECT
476
477 END SUBROUTINE gopt_f_io_init
478
479! **************************************************************************************************
480!> \brief Handles the Output during an optimization run
481!> \param gopt_env ...
482!> \param force_env ...
483!> \param root_section ...
484!> \param its ...
485!> \param opt_energy ...
486!> \param output_unit ...
487!> \param eold ...
488!> \param emin ...
489!> \param wildcard ...
490!> \param gopt_param ...
491!> \param ndf ...
492!> \param dx ...
493!> \param xi ...
494!> \param conv ...
495!> \param pred ...
496!> \param rat ...
497!> \param step ...
498!> \param rad ...
499!> \param used_time ...
500!> \author Teodoro Laino [tlaino] - University of Zurich - 03.2008
501! **************************************************************************************************
502 SUBROUTINE gopt_f_io(gopt_env, force_env, root_section, its, opt_energy, &
503 output_unit, eold, emin, wildcard, gopt_param, ndf, dx, xi, conv, pred, rat, &
504 step, rad, used_time)
505
506 TYPE(gopt_f_type), POINTER :: gopt_env
507 TYPE(force_env_type), POINTER :: force_env
508 TYPE(section_vals_type), POINTER :: root_section
509 INTEGER, INTENT(IN) :: its
510 REAL(kind=dp), INTENT(IN) :: opt_energy
511 INTEGER, INTENT(IN) :: output_unit
512 REAL(kind=dp) :: eold, emin
513 CHARACTER(LEN=5) :: wildcard
514 TYPE(gopt_param_type), POINTER :: gopt_param
515 INTEGER, INTENT(IN), OPTIONAL :: ndf
516 REAL(kind=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: dx
517 REAL(kind=dp), DIMENSION(:), OPTIONAL, POINTER :: xi
518 LOGICAL, OPTIONAL :: conv
519 REAL(kind=dp), INTENT(IN), OPTIONAL :: pred, rat, step, rad
520 REAL(kind=dp) :: used_time
521
522 CHARACTER(LEN=default_string_length) :: energy_unit, stress_unit
523 INTEGER(KIND=int_8) :: max_memory
524 LOGICAL :: print_memory
525 REAL(kind=dp) :: pres_diff, pres_diff_constr, pres_int, &
526 pres_tol
527 TYPE(mp_para_env_type), POINTER :: para_env
528
529 NULLIFY (para_env)
530 CALL section_vals_val_get(gopt_env%motion_section, "PRINT%MEMORY_INFO", l_val=print_memory)
531 max_memory = 0
532 IF (print_memory) THEN
533 CALL force_env_get(force_env, para_env=para_env)
534 max_memory = sample_memory(para_env)
535 END IF
536
537 CALL section_vals_val_get(gopt_env%force_env%force_env_section, &
538 "PRINT%PROGRAM_RUN_INFO%ENERGY_UNIT", &
539 c_val=energy_unit)
540 CALL section_vals_val_get(gopt_env%force_env%force_env_section, &
541 "PRINT%STRESS_TENSOR%STRESS_UNIT", &
542 c_val=stress_unit)
543
544 SELECT CASE (gopt_env%type_id)
545 CASE (default_ts_method_id, default_minimization_method_id)
546 ! Geometry Optimization (Minimization and Transition State Search)
547 IF (.NOT. gopt_env%dimer_rotation) THEN
548 CALL geo_opt_io(force_env=force_env, root_section=root_section, &
549 motion_section=gopt_env%motion_section, its=its, opt_energy=opt_energy)
550 CALL write_cycle_infos(output_unit, &
551 it=its, &
552 etot=opt_energy, &
553 ediff=(opt_energy - eold), &
554 pred=pred, &
555 rat=rat, &
556 step=step, &
557 rad=rad, &
558 emin=emin, &
559 wildcard=wildcard, &
560 used_time=used_time, &
561 max_memory=max_memory, &
562 energy_unit=energy_unit, &
563 stress_unit=stress_unit)
564 ! Possibly check convergence
565 IF (PRESENT(conv)) THEN
566 cpassert(PRESENT(ndf))
567 cpassert(PRESENT(dx))
568 cpassert(PRESENT(xi))
569 CALL check_converg(ndf, dx, xi, output_unit, conv, gopt_param, max_memory, stress_unit)
570 END IF
571 ELSE
572 CALL update_dimer_vec(gopt_env%dimer_env, gopt_env%motion_section)
573 CALL write_restart(force_env=force_env, root_section=root_section)
574 CALL write_rot_cycle_infos(output_unit, its, opt_energy, opt_energy - eold, emin, gopt_env%dimer_env, &
575 wildcard=wildcard, used_time=used_time, max_memory=max_memory)
576 ! Possibly check convergence
577 IF (PRESENT(conv)) THEN
578 cpassert(ASSOCIATED(gopt_env%dimer_env))
579 CALL check_rot_conv(gopt_env%dimer_env, output_unit, conv)
580 END IF
581 END IF
582 CASE (default_cell_method_id)
583 ! Cell Optimization
584 pres_diff = gopt_env%cell_env%pres_int - gopt_env%cell_env%pres_ext
585 pres_int = gopt_env%cell_env%pres_int
586 pres_tol = gopt_env%cell_env%pres_tol
587 CALL geo_opt_io(force_env=force_env, root_section=root_section, &
588 motion_section=gopt_env%motion_section, its=its, opt_energy=opt_energy)
589 CALL write_cycle_infos(output_unit, &
590 it=its, &
591 etot=opt_energy, &
592 ediff=(opt_energy - eold), &
593 pred=pred, &
594 rat=rat, &
595 step=step, &
596 rad=rad, &
597 emin=emin, &
598 pres_int=pres_int, &
599 wildcard=wildcard, &
600 used_time=used_time, &
601 max_memory=max_memory, &
602 energy_unit=energy_unit, &
603 stress_unit=stress_unit)
604 ! Possibly check convergence
605 IF (PRESENT(conv)) THEN
606 cpassert(PRESENT(ndf))
607 cpassert(PRESENT(dx))
608 cpassert(PRESENT(xi))
609 IF (gopt_env%cell_env%constraint_id == fix_none) THEN
610 CALL check_converg(ndf, dx, xi, output_unit, conv, gopt_param, max_memory, stress_unit, &
611 pres_diff, pres_tol)
612 ELSE
613 pres_diff_constr = gopt_env%cell_env%pres_constr - gopt_env%cell_env%pres_ext
614 CALL check_converg(ndf, dx, xi, output_unit, conv, gopt_param, max_memory, stress_unit, &
615 pres_diff, pres_tol, pres_diff_constr)
616 END IF
617 END IF
618 CASE (default_shellcore_method_id)
619 CALL write_cycle_infos(output_unit, &
620 it=its, &
621 etot=opt_energy, &
622 ediff=(opt_energy - eold), &
623 pred=pred, &
624 rat=rat, &
625 step=step, &
626 rad=rad, &
627 emin=emin, &
628 wildcard=wildcard, &
629 used_time=used_time, &
630 max_memory=max_memory, &
631 energy_unit=energy_unit, &
632 stress_unit=stress_unit)
633 ! Possibly check convergence
634 IF (PRESENT(conv)) THEN
635 cpassert(PRESENT(ndf))
636 cpassert(PRESENT(dx))
637 cpassert(PRESENT(xi))
638 CALL check_converg(ndf, dx, xi, output_unit, conv, gopt_param, max_memory, stress_unit)
639 END IF
640 END SELECT
641
642 END SUBROUTINE gopt_f_io
643
644! **************************************************************************************************
645!> \brief Handles the Output at the end of an optimization run
646!> \param gopt_env ...
647!> \param force_env ...
648!> \param x0 ...
649!> \param conv ...
650!> \param its ...
651!> \param root_section ...
652!> \param para_env ...
653!> \param master ...
654!> \param output_unit ...
655!> \author Teodoro Laino [tlaino] - University of Zurich - 03.2008
656! **************************************************************************************************
657 RECURSIVE SUBROUTINE gopt_f_io_finalize(gopt_env, force_env, x0, conv, its, root_section, &
658 para_env, master, output_unit)
659 TYPE(gopt_f_type), POINTER :: gopt_env
660 TYPE(force_env_type), POINTER :: force_env
661 REAL(kind=dp), DIMENSION(:), POINTER :: x0
662 LOGICAL :: conv
663 INTEGER :: its
664 TYPE(section_vals_type), POINTER :: root_section
665 TYPE(mp_para_env_type), POINTER :: para_env
666 INTEGER, INTENT(IN) :: master, output_unit
667
668 IF (gopt_env%eval_opt_geo) THEN
669 IF (.NOT. gopt_env%dimer_rotation) THEN
670 CALL write_final_info(output_unit, conv, its, gopt_env, x0, master, &
671 para_env, force_env, gopt_env%motion_section, root_section)
672 ELSE
673 CALL update_dimer_vec(gopt_env%dimer_env, gopt_env%motion_section)
674 CALL write_restart(force_env=force_env, root_section=root_section)
675 END IF
676 END IF
677
678 END SUBROUTINE gopt_f_io_finalize
679
680! **************************************************************************************************
681!> \brief ...
682!> \param output_unit ...
683!> \param it ...
684!> \param etot ...
685!> \param ediff ...
686!> \param pred ...
687!> \param rat ...
688!> \param step ...
689!> \param rad ...
690!> \param emin ...
691!> \param pres_int ...
692!> \param wildcard ...
693!> \param used_time ...
694!> \param max_memory ...
695!> \param energy_unit ...
696!> \param stress_unit ...
697! **************************************************************************************************
698 SUBROUTINE write_cycle_infos(output_unit, it, etot, ediff, pred, rat, step, rad, emin, &
699 pres_int, wildcard, used_time, max_memory, energy_unit, stress_unit)
700
701 INTEGER, INTENT(IN) :: output_unit, it
702 REAL(kind=dp), INTENT(IN) :: etot
703 REAL(kind=dp), INTENT(IN), OPTIONAL :: ediff, pred, rat, step, rad, emin, &
704 pres_int
705 CHARACTER(LEN=5), INTENT(IN) :: wildcard
706 REAL(kind=dp), INTENT(IN) :: used_time
707 INTEGER(KIND=int_8), INTENT(IN) :: max_memory
708 CHARACTER(LEN=default_string_length), INTENT(IN) :: energy_unit, stress_unit
709
710 CHARACTER(LEN=5) :: tag
711
712 IF (output_unit > 0) THEN
713 tag = "OPT| "
714 WRITE (unit=output_unit, fmt="(/,T2,A)") tag//repeat("*", 74)
715 WRITE (unit=output_unit, fmt="(T2,A,T55,1X,I25)") &
716 tag//"Step number", it
717 WRITE (unit=output_unit, fmt="(T2,A,T55,1X,A25)") &
718 tag//"Optimization method", wildcard
719 WRITE (unit=output_unit, fmt="(T2,A,T55,1X,F25.10)") &
720 tag//"Total energy ["//trim(adjustl(energy_unit))//"]", &
721 cp_unit_from_cp2k(etot, trim(energy_unit))
722 IF (PRESENT(pres_int)) THEN
723 WRITE (unit=output_unit, fmt="(T2,A,T55,1X,F25.10)") &
724 tag//"Internal pressure ["//trim(adjustl(stress_unit))//"]", &
725 cp_unit_from_cp2k(pres_int, trim(stress_unit))
726 END IF
727 IF (PRESENT(ediff)) THEN
728 WRITE (unit=output_unit, fmt="(T2,A,T55,1X,F25.10)") &
729 tag//"Effective energy change ["//trim(adjustl(energy_unit))//"]", &
730 cp_unit_from_cp2k(ediff, trim(energy_unit))
731 END IF
732 IF (PRESENT(pred)) THEN
733 WRITE (unit=output_unit, fmt="(T2,A,T55,1X,F25.10)") &
734 tag//"Predicted energy change ["//trim(adjustl(energy_unit))//"]", &
735 cp_unit_from_cp2k(pred, trim(energy_unit))
736 END IF
737 IF (PRESENT(rat)) THEN
738 WRITE (unit=output_unit, fmt="(T2,A,T55,1X,F25.10)") &
739 tag//"Scaling factor", rat
740 END IF
741 IF (PRESENT(step)) THEN
742 WRITE (unit=output_unit, fmt="(T2,A,T55,1X,F25.10)") &
743 tag//"Step size", step
744 END IF
745 IF (PRESENT(rad)) THEN
746 WRITE (unit=output_unit, fmt="(T2,A,T55,1X,F25.10)") &
747 tag//"Trust radius", rad
748 END IF
749 IF (PRESENT(emin)) THEN
750 IF (etot < emin) THEN
751 WRITE (unit=output_unit, fmt="(T2,A,T77,A4)") &
752 tag//"Decrease in energy", " YES"
753 ELSE
754 WRITE (unit=output_unit, fmt="(T2,A,T77,A4)") &
755 tag//"Decrease in energy", " NO"
756 END IF
757 END IF
758 WRITE (unit=output_unit, fmt="(T2,A,T55,1X,F25.3)") &
759 tag//"Used time [s]", used_time
760 IF (it == 0) THEN
761 WRITE (unit=output_unit, fmt="(T2,A)") tag//repeat("*", 74)
762 IF (max_memory /= 0) THEN
763 WRITE (unit=output_unit, fmt="(T2,A,T60,1X,I20)") &
764 tag//"Estimated peak process memory [MiB]", &
765 (max_memory + (1024*1024) - 1)/(1024*1024)
766 END IF
767 END IF
768 END IF
769
770 END SUBROUTINE write_cycle_infos
771
772! **************************************************************************************************
773!> \brief ...
774!> \param output_unit ...
775!> \param it ...
776!> \param etot ...
777!> \param ediff ...
778!> \param emin ...
779!> \param dimer_env ...
780!> \param used_time ...
781!> \param wildcard ...
782!> \param max_memory ...
783!> \date 01.2008
784!> \author Luca Bellucci and Teodoro Laino - created [tlaino]
785! **************************************************************************************************
786 SUBROUTINE write_rot_cycle_infos(output_unit, it, etot, ediff, emin, dimer_env, used_time, &
787 wildcard, max_memory)
788
789 INTEGER, INTENT(IN) :: output_unit, it
790 REAL(kind=dp), INTENT(IN) :: etot
791 REAL(kind=dp), INTENT(IN), OPTIONAL :: ediff, emin
792 TYPE(dimer_env_type), POINTER :: dimer_env
793 REAL(kind=dp), INTENT(IN) :: used_time
794 CHARACTER(LEN=5), INTENT(IN) :: wildcard
795 INTEGER(KIND=int_8), INTENT(IN) :: max_memory
796
797 CHARACTER(LEN=5) :: tag
798
799 IF (output_unit > 0) THEN
800 tag = "OPT| "
801 WRITE (unit=output_unit, fmt="(/,T2,A)") tag//repeat("*", 74)
802 WRITE (unit=output_unit, fmt="(T2,A,T55,1X,I25)") &
803 tag//"Rotational step number", it
804 WRITE (unit=output_unit, fmt="(T2,A,T55,1X,A25)") &
805 tag//"Optimization method", wildcard
806 WRITE (unit=output_unit, fmt="(T2,A,T55,1X,F25.10)") &
807 tag//"Local curvature", dimer_env%rot%curvature, &
808 tag//"Total rotational force", etot
809 IF (PRESENT(ediff)) THEN
810 WRITE (unit=output_unit, fmt="(T2,A,T55,1X,F25.10)") &
811 tag//"Rotational force change", ediff
812 END IF
813 IF (PRESENT(emin)) THEN
814 IF (etot < emin) THEN
815 WRITE (unit=output_unit, fmt="(T2,A,T77,A4)") &
816 tag//"Decrease in rotational force", " YES"
817 ELSE
818 WRITE (unit=output_unit, fmt="(T2,A,T77,A4)") &
819 tag//"Decrease in rotational force", " NO"
820 END IF
821 END IF
822 WRITE (unit=output_unit, fmt="(T2,A,T55,1X,F25.3)") &
823 tag//"Used time [s]", used_time
824 IF (it == 0) THEN
825 WRITE (unit=output_unit, fmt="(T2,A)") tag//repeat("*", 74)
826 IF (max_memory /= 0) THEN
827 WRITE (unit=output_unit, fmt="(T2,A,T60,1X,I20)") &
828 tag//"Estimated peak process memory [MiB]", &
829 (max_memory + (1024*1024) - 1)/(1024*1024)
830 END IF
831 END IF
832 END IF
833
834 END SUBROUTINE write_rot_cycle_infos
835
836! **************************************************************************************************
837!> \brief ...
838!> \param ndf ...
839!> \param dr ...
840!> \param g ...
841!> \param output_unit ...
842!> \param conv ...
843!> \param gopt_param ...
844!> \param max_memory ...
845!> \param stress_unit ...
846!> \param pres_diff ...
847!> \param pres_tol ...
848!> \param pres_diff_constr ...
849! **************************************************************************************************
850 SUBROUTINE check_converg(ndf, dr, g, output_unit, conv, gopt_param, max_memory, stress_unit, &
851 pres_diff, pres_tol, pres_diff_constr)
852
853 INTEGER, INTENT(IN) :: ndf
854 REAL(kind=dp), INTENT(IN) :: dr(ndf), g(ndf)
855 INTEGER, INTENT(IN) :: output_unit
856 LOGICAL, INTENT(OUT) :: conv
857 TYPE(gopt_param_type), POINTER :: gopt_param
858 INTEGER(KIND=int_8), INTENT(IN) :: max_memory
859 CHARACTER(LEN=default_string_length), INTENT(IN) :: stress_unit
860 REAL(kind=dp), INTENT(IN), OPTIONAL :: pres_diff, pres_tol, pres_diff_constr
861
862 CHARACTER(LEN=5) :: tag
863 INTEGER :: indf
864 LOGICAL :: conv_dx, conv_g, conv_p, conv_rdx, &
865 conv_rg
866 REAL(kind=dp) :: dumm, dxcon, gcon, maxdum(4), rmsgcon, &
867 rmsxcon
868
869 dxcon = gopt_param%max_dr
870 gcon = gopt_param%max_force
871 rmsgcon = gopt_param%rms_force
872 rmsxcon = gopt_param%rms_dr
873
874 conv = .false.
875 conv_dx = .true.
876 conv_rdx = .true.
877 conv_g = .true.
878 conv_rg = .true.
879 conv_p = .true.
880
881 dumm = 0.0_dp
882 DO indf = 1, ndf
883 IF (indf == 1) maxdum(1) = abs(dr(indf))
884 dumm = dumm + dr(indf)**2
885 IF (abs(dr(indf)) > dxcon) conv_dx = .false.
886 IF (abs(dr(indf)) > maxdum(1)) maxdum(1) = abs(dr(indf))
887 END DO
888 ! SQRT(dumm/ndf) > rmsxcon
889 IF (dumm > (rmsxcon*rmsxcon*ndf)) conv_rdx = .false.
890 maxdum(2) = sqrt(dumm/ndf)
891
892 dumm = 0.0_dp
893 DO indf = 1, ndf
894 IF (indf == 1) maxdum(3) = abs(g(indf))
895 dumm = dumm + g(indf)**2
896 IF (abs(g(indf)) > gcon) conv_g = .false.
897 IF (abs(g(indf)) > maxdum(3)) maxdum(3) = abs(g(indf))
898 END DO
899 ! SQRT(dumm/ndf) > rmsgcon
900 IF (dumm > (rmsgcon*rmsgcon*ndf)) conv_rg = .false.
901 maxdum(4) = sqrt(dumm/ndf)
902
903 IF (PRESENT(pres_diff_constr) .AND. PRESENT(pres_tol)) THEN
904 conv_p = abs(pres_diff_constr) < abs(pres_tol)
905 ELSE IF (PRESENT(pres_diff) .AND. PRESENT(pres_tol)) THEN
906 conv_p = abs(pres_diff) < abs(pres_tol)
907 END IF
908
909 IF (output_unit > 0) THEN
910
911 tag = "OPT| "
912
913 WRITE (unit=output_unit, fmt="(T2,A)") trim(tag)
914 WRITE (unit=output_unit, fmt="(T2,A,T55,1X,F25.10)") &
915 tag//"Maximum step size", maxdum(1), &
916 tag//"Convergence limit for maximum step size", dxcon
917 IF (conv_dx) THEN
918 WRITE (unit=output_unit, fmt="(T2,A,T77,A4)") &
919 tag//"Maximum step size is converged", " YES"
920 ELSE
921 WRITE (unit=output_unit, fmt="(T2,A,T77,A4)") &
922 tag//"Maximum step size is converged", " NO"
923 END IF
924
925 WRITE (unit=output_unit, fmt="(T2,A)") trim(tag)
926 WRITE (unit=output_unit, fmt="(T2,A,T55,1X,F25.10)") &
927 tag//"RMS step size", maxdum(2), &
928 tag//"Convergence limit for RMS step size", rmsxcon
929 IF (conv_rdx) THEN
930 WRITE (unit=output_unit, fmt="(T2,A,T77,A4)") &
931 tag//"RMS step size is converged", " YES"
932 ELSE
933 WRITE (unit=output_unit, fmt="(T2,A,T77,A4)") &
934 tag//"RMS step size is converged", " NO"
935 END IF
936
937 WRITE (unit=output_unit, fmt="(T2,A)") trim(tag)
938 WRITE (unit=output_unit, fmt="(T2,A,T55,1X,F25.10)") &
939 tag//"Maximum gradient", maxdum(3), &
940 tag//"Convergence limit for maximum gradient", gcon
941 IF (conv_g) THEN
942 WRITE (unit=output_unit, fmt="(T2,A,T77,A4)") &
943 tag//"Maximum gradient is converged", " YES"
944 ELSE
945 WRITE (unit=output_unit, fmt="(T2,A,T77,A4)") &
946 tag//"Maximum gradient is converged", " NO"
947 END IF
948
949 WRITE (unit=output_unit, fmt="(T2,A)") trim(tag)
950 WRITE (unit=output_unit, fmt="(T2,A,T55,1X,F25.10)") &
951 tag//"RMS gradient", maxdum(4), &
952 tag//"Convergence limit for RMS gradient", rmsgcon
953 IF (conv_rg) THEN
954 WRITE (unit=output_unit, fmt="(T2,A,T77,A4)") &
955 tag//"RMS gradient is converged", " YES"
956 ELSE
957 WRITE (unit=output_unit, fmt="(T2,A,T77,A4)") &
958 tag//"RMS gradient is converged", " NO"
959 END IF
960
961 IF (PRESENT(pres_diff) .AND. PRESENT(pres_tol)) THEN
962 WRITE (unit=output_unit, fmt="(T2,A)") trim(tag)
963 IF (PRESENT(pres_diff_constr)) THEN
964 WRITE (unit=output_unit, fmt="(T2,A,T55,1X,F25.10)") &
965 tag//"Pressure deviation without constraint ["// &
966 trim(adjustl(stress_unit))//"]", &
967 cp_unit_from_cp2k(pres_diff, trim(stress_unit))
968 WRITE (unit=output_unit, fmt="(T2,A,T55,1X,F25.10)") &
969 tag//"Pressure deviation with constraint ["// &
970 trim(adjustl(stress_unit))//"]", &
971 cp_unit_from_cp2k(pres_diff_constr, trim(stress_unit))
972 ELSE
973 WRITE (unit=output_unit, fmt="(T2,A,T55,1X,F25.10)") &
974 tag//"Pressure deviation ["//trim(adjustl(stress_unit))//"]", &
975 cp_unit_from_cp2k(pres_diff, trim(stress_unit))
976 END IF
977 WRITE (unit=output_unit, fmt="(T2,A,T55,1X,F25.10)") &
978 tag//"Pressure tolerance ["//trim(adjustl(stress_unit))//"]", &
979 cp_unit_from_cp2k(pres_tol, trim(stress_unit))
980 IF (conv_p) THEN
981 WRITE (unit=output_unit, fmt="(T2,A,T77,A4)") &
982 tag//"Pressure is converged", " YES"
983 ELSE
984 WRITE (unit=output_unit, fmt="(T2,A,T77,A4)") &
985 tag//"Pressure is converged", " NO"
986 END IF
987 END IF
988
989 WRITE (unit=output_unit, fmt="(T2,A)") tag//repeat("*", 74)
990
991 IF (max_memory /= 0) THEN
992 WRITE (unit=output_unit, fmt="(T2,A,T60,1X,I20)") &
993 tag//"Estimated peak process memory after this step [MiB]", &
994 (max_memory + (1024*1024) - 1)/(1024*1024)
995 END IF
996
997 END IF
998
999 IF (conv_dx .AND. conv_rdx .AND. conv_g .AND. conv_rg .AND. conv_p) conv = .true.
1000
1001 IF ((conv) .AND. (output_unit > 0)) THEN
1002 WRITE (unit=output_unit, fmt="(/,T2,A)") repeat("*", 79)
1003 WRITE (unit=output_unit, fmt="(T2,A,T25,A,T78,A)") &
1004 "***", "GEOMETRY OPTIMIZATION COMPLETED", "***"
1005 WRITE (unit=output_unit, fmt="(T2,A)") repeat("*", 79)
1006 END IF
1007
1008 END SUBROUTINE check_converg
1009
1010! **************************************************************************************************
1011!> \brief ...
1012!> \param dimer_env ...
1013!> \param output_unit ...
1014!> \param conv ...
1015!> \date 01.2008
1016!> \author Luca Bellucci and Teodoro Laino - created [tlaino]
1017! **************************************************************************************************
1018 SUBROUTINE check_rot_conv(dimer_env, output_unit, conv)
1019
1020 TYPE(dimer_env_type), POINTER :: dimer_env
1021 INTEGER, INTENT(IN) :: output_unit
1022 LOGICAL, INTENT(OUT) :: conv
1023
1024 CHARACTER(LEN=5) :: tag
1025
1026 conv = (abs(dimer_env%rot%angle2) < dimer_env%rot%angle_tol)
1027
1028 IF (output_unit > 0) THEN
1029 tag = "OPT| "
1030 WRITE (unit=output_unit, fmt="(T2,A)") trim(tag)
1031 WRITE (unit=output_unit, fmt="(T2,A,T55,1X,F25.10)") &
1032 tag//"Predicted angle step size", dimer_env%rot%angle1, &
1033 tag//"Effective angle step size", dimer_env%rot%angle2, &
1034 tag//"Convergence limit for angle step size", dimer_env%rot%angle_tol
1035 IF (conv) THEN
1036 WRITE (unit=output_unit, fmt="(T2,A,T77,A4)") &
1037 tag//"Angle step size is converged", " YES"
1038 ELSE
1039 WRITE (unit=output_unit, fmt="(T2,A,T77,A4)") &
1040 tag//"Angle step size is converged", " NO"
1041 END IF
1042 WRITE (unit=output_unit, fmt="(T2,A)") tag//repeat("*", 74)
1043 END IF
1044
1045 IF ((conv) .AND. (output_unit > 0)) THEN
1046 WRITE (unit=output_unit, fmt="(/,T2,A)") repeat("*", 79)
1047 WRITE (unit=output_unit, fmt="(T2,A,T25,A,T78,A)") &
1048 "***", "ROTATION OPTIMIZATION COMPLETED", "***"
1049 WRITE (unit=output_unit, fmt="(T2,A)") repeat("*", 79)
1050 END IF
1051
1052 END SUBROUTINE check_rot_conv
1053
1054! **************************************************************************************************
1055!> \brief ...
1056!> \param output_unit ...
1057!> \param conv ...
1058!> \param it ...
1059!> \param gopt_env ...
1060!> \param x0 ...
1061!> \param master ...
1062!> \param para_env ...
1063!> \param force_env ...
1064!> \param motion_section ...
1065!> \param root_section ...
1066!> \date 11.2007
1067!> \author Teodoro Laino [tlaino] - University of Zurich
1068! **************************************************************************************************
1069 RECURSIVE SUBROUTINE write_final_info(output_unit, conv, it, gopt_env, x0, master, para_env, force_env, &
1070 motion_section, root_section)
1071 INTEGER, INTENT(IN) :: output_unit
1072 LOGICAL, INTENT(IN) :: conv
1073 INTEGER, INTENT(INOUT) :: it
1074 TYPE(gopt_f_type), POINTER :: gopt_env
1075 REAL(kind=dp), DIMENSION(:), POINTER :: x0
1076 INTEGER, INTENT(IN) :: master
1077 TYPE(mp_para_env_type), POINTER :: para_env
1078 TYPE(force_env_type), POINTER :: force_env
1079 TYPE(section_vals_type), POINTER :: motion_section, root_section
1080
1081 CHARACTER(LEN=4) :: constraint_label
1082 LOGICAL :: keep_angles, keep_symmetry, keep_volume
1083 REAL(kind=dp) :: etot
1084 TYPE(cell_type), POINTER :: cell
1085 TYPE(cp_subsys_type), POINTER :: subsys
1086 TYPE(particle_list_type), POINTER :: particles
1087 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1088
1089 CALL force_env_get(force_env, cell=cell, subsys=subsys)
1090 CALL cp_subsys_get(subsys=subsys, particles=particles)
1091 particle_set => particles%els
1092
1093 ! Passing gopt_f_type pointer gopt_env to particle_methods where
1094 ! write_final_structure is defined causes a circular dependency, so it
1095 ! is necessary to get some flags by preprocessing...
1096 keep_angles = .true.
1097 keep_symmetry = .true.
1098 keep_volume = .true.
1099 constraint_label = "NONE"
1100 IF (gopt_env%type_id == default_cell_method_id) THEN
1101 keep_angles = gopt_env%cell_env%keep_angles
1102 keep_symmetry = gopt_env%cell_env%keep_symmetry
1103 keep_volume = gopt_env%cell_env%keep_volume
1104 SELECT CASE (gopt_env%cell_env%constraint_id)
1105 CASE (fix_x)
1106 constraint_label = " X"
1107 CASE (fix_y)
1108 constraint_label = " Y"
1109 CASE (fix_z)
1110 constraint_label = " Z"
1111 CASE (fix_xy)
1112 constraint_label = " XY"
1113 CASE (fix_xz)
1114 constraint_label = " XZ"
1115 CASE (fix_yz)
1116 constraint_label = " YZ"
1117 CASE (fix_none)
1118 constraint_label = "NONE"
1119 END SELECT
1120 END IF
1121 CALL write_final_structure(particle_set, cell, motion_section, conv, &
1122 keep_angles, keep_symmetry, keep_volume, &
1123 gopt_env%label, constraint_label)
1124
1125 IF (conv) THEN
1126 it = it + 1
1127 CALL write_structure_data(particle_set, cell, motion_section)
1128 CALL write_restart(force_env=force_env, root_section=root_section)
1129
1130 IF (output_unit > 0) THEN
1131 WRITE (unit=output_unit, fmt="(/,T20,' Reevaluating energy at the minimum')")
1132 END IF
1133
1134 CALL cp_eval_at(gopt_env, x0, f=etot, master=master, final_evaluation=.true., &
1135 para_env=para_env)
1136 CALL write_geo_traj(force_env, root_section, it, etot)
1137 END IF
1138
1139 END SUBROUTINE write_final_info
1140
1141! **************************************************************************************************
1142!> \brief Specific driver for dumping trajectory during a GEO_OPT
1143!> \param force_env ...
1144!> \param root_section ...
1145!> \param it ...
1146!> \param etot ...
1147!> \date 11.2007
1148!> \par History
1149!> 09.2010: Output of core and shell positions and forces (MK)
1150!> \author Teodoro Laino [tlaino] - University of Zurich
1151! **************************************************************************************************
1152 SUBROUTINE write_geo_traj(force_env, root_section, it, etot)
1153
1154 TYPE(force_env_type), POINTER :: force_env
1155 TYPE(section_vals_type), POINTER :: root_section
1156 INTEGER, INTENT(IN) :: it
1157 REAL(kind=dp), INTENT(IN) :: etot
1158
1159 LOGICAL :: shell_adiabatic, shell_present
1160 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
1161 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1162 TYPE(cp_subsys_type), POINTER :: subsys
1163 TYPE(particle_list_type), POINTER :: core_particles, shell_particles
1164
1165 NULLIFY (atomic_kinds)
1166 NULLIFY (atomic_kind_set)
1167 NULLIFY (core_particles)
1168 NULLIFY (shell_particles)
1169 NULLIFY (subsys)
1170
1171 CALL write_trajectory(force_env, root_section, it, 0.0_dp, 0.0_dp, etot)
1172 ! Print Force
1173 CALL write_trajectory(force_env, root_section, it, 0.0_dp, 0.0_dp, etot, "FORCES", middle_name="frc")
1174 CALL force_env_get(force_env, subsys=subsys)
1175 CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds)
1176 atomic_kind_set => atomic_kinds%els
1177 CALL get_atomic_kind_set(atomic_kind_set, &
1178 shell_present=shell_present, &
1179 shell_adiabatic=shell_adiabatic)
1180 IF (shell_present) THEN
1181 CALL cp_subsys_get(subsys, &
1182 core_particles=core_particles, &
1183 shell_particles=shell_particles)
1184 CALL write_trajectory(force_env, root_section, it=it, time=0.0_dp, dtime=0.0_dp, &
1185 etot=etot, pk_name="SHELL_TRAJECTORY", middle_name="shpos", &
1186 particles=shell_particles)
1187 IF (shell_adiabatic) THEN
1188 CALL write_trajectory(force_env, root_section, it=it, time=0.0_dp, dtime=0.0_dp, &
1189 etot=etot, pk_name="SHELL_FORCES", middle_name="shfrc", &
1190 particles=shell_particles)
1191 CALL write_trajectory(force_env, root_section, it=it, time=0.0_dp, dtime=0.0_dp, &
1192 etot=etot, pk_name="CORE_TRAJECTORY", middle_name="copos", &
1193 particles=core_particles)
1194 CALL write_trajectory(force_env, root_section, it=it, time=0.0_dp, dtime=0.0_dp, &
1195 etot=etot, pk_name="CORE_FORCES", middle_name="cofrc", &
1196 particles=core_particles)
1197 END IF
1198 END IF
1199
1200 END SUBROUTINE write_geo_traj
1201
1202! **************************************************************************************************
1203!> \brief ...
1204!> \param gopt_env ...
1205!> \param output_unit ...
1206!> \param label ...
1207!> \date 01.2008
1208!> \author Teodoro Laino [tlaino] - University of Zurich
1209! **************************************************************************************************
1210 SUBROUTINE print_geo_opt_header(gopt_env, output_unit, label)
1211
1212 TYPE(gopt_f_type), POINTER :: gopt_env
1213 INTEGER, INTENT(IN) :: output_unit
1214 CHARACTER(LEN=*), INTENT(IN) :: label
1215
1216 CHARACTER(LEN=default_string_length) :: my_format, my_label
1217 INTEGER :: ix
1218
1219 IF (output_unit > 0) THEN
1220 WRITE (unit=output_unit, fmt="(/,T2,A)") repeat("*", 79)
1221 IF (gopt_env%dimer_rotation) THEN
1222 my_label = "OPTIMIZING DIMER ROTATION"
1223 ELSE
1224 my_label = "STARTING "//gopt_env%tag(1:8)//" OPTIMIZATION"
1225 END IF
1226
1227 ix = (80 - 7 - len_trim(my_label))/2
1228 ix = ix + 5
1229 my_format = "(T2,A,T"//cp_to_string(ix)//",A,T78,A)"
1230 WRITE (unit=output_unit, fmt=trim(my_format)) "***", trim(my_label), "***"
1231
1232 ix = (80 - 7 - len_trim(label))/2
1233 ix = ix + 5
1234 my_format = "(T2,A,T"//cp_to_string(ix)//",A,T78,A)"
1235 WRITE (unit=output_unit, fmt=trim(my_format)) "***", trim(label), "***"
1236
1237 WRITE (unit=output_unit, fmt="(T2,A)") repeat("*", 79)
1238 CALL m_flush(output_unit)
1239 END IF
1240 END SUBROUTINE print_geo_opt_header
1241
1242! **************************************************************************************************
1243!> \brief ...
1244!> \param gopt_env ...
1245!> \param output_unit ...
1246!> \date 01.2008
1247!> \author Teodoro Laino [tlaino] - University of Zurich
1248! **************************************************************************************************
1249 SUBROUTINE print_geo_opt_nc(gopt_env, output_unit)
1250
1251 TYPE(gopt_f_type), POINTER :: gopt_env
1252 INTEGER, INTENT(IN) :: output_unit
1253
1254 IF (output_unit > 0) THEN
1255 WRITE (unit=output_unit, fmt="(/,T2,A)") &
1256 "*** MAXIMUM NUMBER OF OPTIMIZATION STEPS REACHED ***"
1257 IF (.NOT. gopt_env%dimer_rotation) THEN
1258 WRITE (unit=output_unit, fmt="(T2,A)") &
1259 "*** EXITING GEOMETRY OPTIMIZATION ***"
1260 ELSE
1261 WRITE (unit=output_unit, fmt="(T2,A)") &
1262 "*** EXITING ROTATION OPTIMIZATION ***"
1263 END IF
1264 CALL m_flush(output_unit)
1265 END IF
1266
1267 END SUBROUTINE print_geo_opt_nc
1268
1269! **************************************************************************************************
1270!> \brief Prints information during GEO_OPT common to all optimizers
1271!> \param force_env ...
1272!> \param root_section ...
1273!> \param motion_section ...
1274!> \param its ...
1275!> \param opt_energy ...
1276!> \date 02.2008
1277!> \author Teodoro Laino [tlaino] - University of Zurich
1278!> \version 1.0
1279! **************************************************************************************************
1280 SUBROUTINE geo_opt_io(force_env, root_section, motion_section, its, opt_energy)
1281
1282 TYPE(force_env_type), POINTER :: force_env
1283 TYPE(section_vals_type), POINTER :: root_section, motion_section
1284 INTEGER, INTENT(IN) :: its
1285 REAL(kind=dp), INTENT(IN) :: opt_energy
1286
1287 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
1288 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1289 TYPE(cell_type), POINTER :: cell
1290 TYPE(cp_subsys_type), POINTER :: subsys
1291 TYPE(distribution_1d_type), POINTER :: local_particles
1292 TYPE(mp_para_env_type), POINTER :: para_env
1293 TYPE(particle_list_type), POINTER :: particles
1294 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1295 TYPE(virial_type), POINTER :: virial
1296
1297 NULLIFY (para_env, atomic_kind_set, subsys, particle_set, &
1298 local_particles, atomic_kinds, particles)
1299
1300 ! Write Restart File
1301 CALL write_restart(force_env=force_env, root_section=root_section)
1302
1303 ! Write Trajectory
1304 CALL write_geo_traj(force_env, root_section, its, opt_energy)
1305
1306 ! Write the stress Tensor
1307 CALL force_env_get(force_env, cell=cell, para_env=para_env, &
1308 subsys=subsys)
1309 CALL cp_subsys_get(subsys=subsys, atomic_kinds=atomic_kinds, local_particles=local_particles, &
1310 particles=particles, virial=virial)
1311 atomic_kind_set => atomic_kinds%els
1312 particle_set => particles%els
1313 CALL virial_evaluate(atomic_kind_set, particle_set, local_particles, &
1314 virial, para_env)
1315 CALL write_stress_tensor_to_file(virial, cell, motion_section, its, 0.0_dp)
1316
1317 ! Write the cell
1318 CALL write_simulation_cell(cell, motion_section, its, 0.0_dp)
1319
1320 END SUBROUTINE geo_opt_io
1321
1322! **************************************************************************************************
1323!> \brief Apply coordinate transformations after cell (shape) change
1324!> \param gopt_env ...
1325!> \param cell ...
1326!> \param x ...
1327!> \param update_forces ...
1328!> \date 05.11.2012 (revised version of unbiase_coordinates moved here, MK)
1329!> \author Matthias Krack
1330!> \version 1.0
1331! **************************************************************************************************
1332 SUBROUTINE apply_cell_change(gopt_env, cell, x, update_forces)
1333
1334 TYPE(gopt_f_type), POINTER :: gopt_env
1335 TYPE(cell_type), POINTER :: cell
1336 REAL(kind=dp), DIMENSION(:), POINTER :: x
1337 LOGICAL, INTENT(IN) :: update_forces
1338
1339 INTEGER :: i, iatom, idg, j, natom, nparticle, &
1340 shell_index
1341 REAL(kind=dp) :: fc, fs, mass
1342 REAL(kind=dp), DIMENSION(3) :: s
1343 TYPE(cell_type), POINTER :: cell_ref
1344 TYPE(cp_subsys_type), POINTER :: subsys
1345 TYPE(particle_list_type), POINTER :: core_particles, particles, &
1346 shell_particles
1347
1348 NULLIFY (cell_ref)
1349 NULLIFY (core_particles)
1350 NULLIFY (particles)
1351 NULLIFY (shell_particles)
1352 NULLIFY (subsys)
1353
1354 natom = force_env_get_natom(gopt_env%force_env)
1355 nparticle = force_env_get_nparticle(gopt_env%force_env)
1356 CALL force_env_get(gopt_env%force_env, &
1357 subsys=subsys)
1358 CALL cp_subsys_get(subsys=subsys, &
1359 core_particles=core_particles, &
1360 particles=particles, &
1361 shell_particles=shell_particles)
1362
1363 ! Retrieve the reference cell
1364 CALL cell_create(cell_ref)
1365 CALL cell_copy(cell, cell_ref, tag="CELL_OPT_REF")
1366
1367 ! Load the updated cell information
1368 idg = 3*nparticle
1369 CALL init_cell(cell_ref, hmat=gopt_env%h_ref)
1370 cpassert((SIZE(x) == idg + 6))
1371
1372 IF (update_forces) THEN
1373
1374 ! Transform particle forces back to reference cell
1375 idg = 1
1376 DO iatom = 1, natom
1377 CALL real_to_scaled(s, x(idg:idg + 2), cell)
1378 CALL scaled_to_real(x(idg:idg + 2), s, cell_ref)
1379 idg = idg + 3
1380 END DO
1381
1382 ELSE
1383
1384 ! Update the six independent components in the canonical optimization frame.
1385 gopt_env%cell_env%opt_cell%hmat = 0.0_dp
1386 DO i = 1, 3
1387 DO j = 1, i
1388 idg = idg + 1
1389 gopt_env%cell_env%opt_cell%hmat(j, i) = x(idg)
1390 END DO
1391 END DO
1392 CALL init_cell(gopt_env%cell_env%opt_cell)
1393 IF (gopt_env%spgr%keep_space_group) THEN
1394 CALL spgr_project_cell_metric(gopt_env%spgr, gopt_env%cell_env%opt_cell)
1395 ! Keep the optimizer variables synchronized with the projected cell.
1396 idg = 3*nparticle
1397 DO i = 1, 3
1398 DO j = 1, i
1399 idg = idg + 1
1400 x(idg) = gopt_env%cell_env%opt_cell%hmat(j, i)
1401 END DO
1402 END DO
1403 END IF
1404
1405 ! Reconstruct the physical cell in the fixed orientation of the input cell.
1406 cell%hmat = matmul(gopt_env%cell_env%opt_to_input, &
1407 gopt_env%cell_env%opt_cell%hmat)
1408 CALL init_cell(cell)
1409 CALL cp_subsys_set(subsys, cell=cell)
1410
1411 ! Retrieve particle coordinates for the current cell
1412 idg = 1
1413 DO iatom = 1, natom
1414 CALL real_to_scaled(s, x(idg:idg + 2), cell_ref)
1415 shell_index = particles%els(iatom)%shell_index
1416 IF (shell_index == 0) THEN
1417 CALL scaled_to_real(particles%els(iatom)%r, s, cell)
1418 ELSE
1419 CALL scaled_to_real(core_particles%els(shell_index)%r, s, cell)
1420 i = 3*(natom + shell_index - 1) + 1
1421 CALL real_to_scaled(s, x(i:i + 2), cell_ref)
1422 CALL scaled_to_real(shell_particles%els(shell_index)%r, s, cell)
1423 ! Update atomic position due to core and shell motion
1424 mass = particles%els(iatom)%atomic_kind%mass
1425 fc = core_particles%els(shell_index)%atomic_kind%shell%mass_core/mass
1426 fs = shell_particles%els(shell_index)%atomic_kind%shell%mass_shell/mass
1427 particles%els(iatom)%r(1:3) = fc*core_particles%els(shell_index)%r(1:3) + &
1428 fs*shell_particles%els(shell_index)%r(1:3)
1429 END IF
1430 idg = idg + 3
1431 END DO
1432 END IF
1433
1434 CALL cell_release(cell_ref)
1435
1436 END SUBROUTINE apply_cell_change
1437
1438END MODULE gopt_f_methods
represent a simple array based list of the given type
Define the atomic kind types and their sub types.
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.
collects all references to literature in CP2K as new algorithms / method are included from literature...
integer, save, public henkelman1999
Handles all functions related to the CELL.
subroutine, public write_cell(cell, subsys_section, tag)
Write the cell parameters to the output unit.
subroutine, public init_cell(cell, hmat, periodic)
Initialise/readjust a simulation cell after hmat has been changed.
subroutine, public cell_create(cell, hmat, periodic, tag)
allocates and initializes a cell
contains a functional that calculates the energy and its derivatives for the geometry optimizer
subroutine, public rescale_new_cell_volume(vol, x, idg)
Rescale x(idg+1:idg+6) according to vol.
subroutine, public get_dg_dh(gradient, av_ptens, pres_ext, cell, mtrx, keep_angles, keep_symmetry, pres_int, pres_constr, constraint_id)
Computes the derivatives for the cell.
Handles all functions related to the CELL.
Definition cell_types.F:15
subroutine, public scaled_to_real(r, s, cell)
Transform scaled cell coordinates real coordinates. r=h*s.
Definition cell_types.F:625
subroutine, public real_to_scaled(s, r, cell)
Transform real to scaled cell coordinates. s=h_inv*r.
Definition cell_types.F:595
subroutine, public cell_release(cell)
releases the given cell (see doc/ReferenceCounting.html)
Definition cell_types.F:668
subroutine, public cell_copy(cell_in, cell_out, tag)
Copy cell variable.
Definition cell_types.F:160
various routines to log and control the output. The idea is that decisions about where to log should ...
types that represent a subsys, i.e. a part of the system
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
subroutine, public unpack_subsys_particles(subsys, f, r, s, v, fscale, cell)
Unpack components of a subsystem particle sets into a single vector.
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
subroutine, public pack_subsys_particles(subsys, f, r, s, v, fscale, cell)
Pack components of a subsystem particle sets into a single vector.
unit conversion facility
Definition cp_units.F:30
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
Contains types used for a Dimer Method calculations.
recursive subroutine, public cp_eval_at_ts(gopt_env, x, f, gradient, calc_force)
Computes the dimer energy/gradients (including the rotation of the dimer)
Contains types used for a Dimer Method calculations.
Definition dimer_types.F:14
Contains utilities for a Dimer Method calculations.
Definition dimer_utils.F:14
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
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
Interface for the force calculations.
recursive subroutine, public force_env_calc_energy_force(force_env, calc_force, consistent_energies, skip_external_control, eval_energy_forces, require_consistent_energy_force, linres, calc_stress_tensor)
Interface routine for force and energy calculations.
Interface for the force calculations.
integer function, public force_env_get_natom(force_env)
returns the number of atoms
integer, parameter, public use_qmmm
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
integer, parameter, public use_qmmmx
integer function, public force_env_get_nparticle(force_env)
returns the number of particles in a force environment
contains a functional that calculates the energy and its derivatives for the geometry optimizer
subroutine, public print_geo_opt_header(gopt_env, output_unit, label)
...
subroutine, public gopt_f_io_init(gopt_env, output_unit, opt_energy, wildcard, its, used_time)
Handles the Output during an optimization run.
subroutine, public gopt_f_create_x0(gopt_env, x0)
returns the value of the parameters for the actual configuration
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.
subroutine, public apply_cell_change(gopt_env, cell, x, update_forces)
Apply coordinate transformations after cell (shape) change.
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.
subroutine, public print_geo_opt_nc(gopt_env, output_unit)
...
subroutine, public 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...
subroutine, public gopt_f_ii(its, output_unit)
Prints iteration step of the optimization procedure on screen.
contains a functional that calculates the energy and its derivatives for the geometry optimizer
contains typo and related routines to handle parameters controlling the GEO_OPT module
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public default_shellcore_method_id
integer, parameter, public fix_xz
integer, parameter, public default_cell_method_id
integer, parameter, public default_minimization_method_id
integer, parameter, public default_ts_method_id
integer, parameter, public fix_y
integer, parameter, public fix_none
integer, parameter, public fix_z
integer, parameter, public fix_xy
integer, parameter, public fix_yz
integer, parameter, public fix_x
Set of routines to dump the restart file of CP2K.
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....
objects that represent the structure of input sections and the data contained in an input section
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
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public int_8
Definition kinds.F:54
integer, parameter, public dp
Definition kinds.F:34
integer, parameter, public default_string_length
Definition kinds.F:57
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition machine.F:124
prints all energy info per timestep to the screen or to user defined output files
Definition md_energies.F:16
integer(kind=int_8) function, public sample_memory(para_env)
Samples memory usage.
Interface to the message passing library MPI.
Output Utilities for MOTION_SECTION.
subroutine, public write_simulation_cell(cell, motion_section, itimes, time, pos, act)
Prints the Simulation Cell.
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.
subroutine, public write_stress_tensor_to_file(virial, cell, motion_section, itimes, time, pos, act)
Prints the Stress Tensor.
represent a simple array based list of the given type
Define methods related to particle_type.
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.
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....
Define the data structure for the particle information.
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
Routines used for force-mixing QM/MM calculations.
Definition qmmmx_util.F:14
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
Space Group Symmetry Type Module (version 1.0, Ferbruary 12, 2021)
Space Group Symmetry Module (version 1.0, January 16, 2020)
subroutine, public spgr_apply_rotations_stress(spgr, cell, stress)
routine applies the rotation matrices to the stress tensor.
subroutine, public spgr_project_cell_metric(spgr, cell)
Project a cell onto the metric preserved by the selected space-group rotations.
subroutine, public spgr_apply_rotations_force(spgr, force)
routine applies the rotation matrices to the forces.
subroutine, public spgr_write_stress_tensor(stress, spgr)
Variable precision output of the symmetrized stress tensor.
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)
subroutine, public virial_update(virial, subsys, para_env)
Updates the virial given the virial and subsys.
Provides all information about an atomic kind.
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
represents a system: atoms, molecules, their pos,vel,...
Defines the environment for a Dimer Method calculation.
structure to store local (to a processor) ordered lists of integers.
wrapper to abstract the force evaluation of the various methods
calculates the potential energy of a system, and its derivatives
stores all the informations relevant to an mpi environment