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