Line data Source code
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 : ! **************************************************************************************************
14 : MODULE gopt_f_methods
15 :
16 : USE atomic_kind_list_types, ONLY: atomic_kind_list_type
17 : USE atomic_kind_types, ONLY: atomic_kind_type,&
18 : get_atomic_kind_set
19 : USE bibliography, ONLY: Henkelman1999,&
20 : cite_reference
21 : USE cell_methods, ONLY: cell_create,&
22 : init_cell,&
23 : write_cell
24 : USE cell_opt_utils, ONLY: apply_cell_constraints,&
25 : cell_opt_constraint_label,&
26 : get_dg_dh,&
27 : rescale_new_cell_volume
28 : USE cell_types, ONLY: cell_copy,&
29 : cell_release,&
30 : cell_type,&
31 : real_to_scaled,&
32 : scaled_to_real
33 : USE cp_log_handling, ONLY: cp_to_string
34 : USE cp_subsys_types, ONLY: cp_subsys_get,&
35 : cp_subsys_set,&
36 : cp_subsys_type,&
37 : pack_subsys_particles,&
38 : unpack_subsys_particles
39 : USE cp_units, ONLY: cp_unit_from_cp2k
40 : USE dimer_methods, ONLY: cp_eval_at_ts
41 : USE dimer_types, ONLY: dimer_env_type
42 : USE dimer_utils, ONLY: update_dimer_vec
43 : USE distribution_1d_types, ONLY: distribution_1d_type
44 : USE force_env_methods, ONLY: force_env_calc_energy_force
45 : USE force_env_types, ONLY: force_env_get,&
46 : force_env_get_natom,&
47 : force_env_get_nparticle,&
48 : force_env_type,&
49 : use_qmmm,&
50 : use_qmmmx
51 : USE gopt_f_types, ONLY: gopt_f_type
52 : USE gopt_param_types, ONLY: gopt_param_type
53 : USE input_constants, ONLY: default_cell_method_id,&
54 : default_minimization_method_id,&
55 : default_shellcore_method_id,&
56 : default_ts_method_id,&
57 : fix_none
58 : USE input_cp2k_restarts, ONLY: write_restart
59 : USE input_section_types, ONLY: section_vals_type,&
60 : section_vals_val_get
61 : USE kinds, ONLY: default_string_length,&
62 : dp,&
63 : int_8
64 : USE machine, ONLY: m_flush
65 : USE md_energies, ONLY: sample_memory
66 : USE message_passing, ONLY: mp_para_env_type
67 : USE motion_utils, ONLY: write_simulation_cell,&
68 : write_stress_tensor_to_file,&
69 : write_trajectory
70 : USE particle_list_types, ONLY: particle_list_type
71 : USE particle_methods, ONLY: write_final_structure,&
72 : write_structure_data
73 : USE particle_types, ONLY: particle_type
74 : USE qmmm_util, ONLY: apply_qmmm_translate
75 : USE qmmmx_util, ONLY: apply_qmmmx_translate
76 : USE space_groups, ONLY: identify_space_group,&
77 : print_spgr,&
78 : spgr_apply_rotations_coord,&
79 : spgr_apply_rotations_force,&
80 : spgr_apply_rotations_stress,&
81 : spgr_project_cell_metric,&
82 : spgr_write_stress_tensor
83 : USE space_groups_types, ONLY: spgr_type
84 : USE virial_methods, ONLY: virial_evaluate,&
85 : virial_update
86 : USE virial_types, ONLY: virial_type
87 : #include "../base/base_uses.f90"
88 :
89 : IMPLICIT NONE
90 : PRIVATE
91 :
92 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gopt_f_methods'
93 : INTEGER(KIND=int_8), PARAMETER, PRIVATE :: bytes_per_mib = 1024*1024
94 :
95 : PUBLIC :: cp_eval_at, &
96 : gopt_f_create_x0, &
97 : gopt_f_project_coord, gopt_f_project_step, gopt_f_project_vector, &
98 : gopt_f_space_group_init, gopt_f_space_group_show, &
99 : print_geo_opt_header, print_geo_opt_nc, &
100 : gopt_f_io_init, gopt_f_io, gopt_f_io_finalize, gopt_f_ii, &
101 : apply_cell_change
102 :
103 : CONTAINS
104 :
105 : ! **************************************************************************************************
106 : !> \brief returns the value of the parameters for the actual configuration
107 : !> \param gopt_env the geometry optimization environment you want the info about
108 : !> x0: the parameter vector (is allocated by this routine)
109 : !> \param x0 ...
110 : !> \par History
111 : !> - Cell optimization revised (06.11.2012,MK)
112 : ! **************************************************************************************************
113 1119 : SUBROUTINE gopt_f_create_x0(gopt_env, x0)
114 :
115 : TYPE(gopt_f_type), POINTER :: gopt_env
116 : REAL(KIND=dp), DIMENSION(:), POINTER :: x0
117 :
118 : INTEGER :: i, idg, j, nparticle
119 : TYPE(cell_type), POINTER :: cell
120 : TYPE(cp_subsys_type), POINTER :: subsys
121 :
122 1119 : NULLIFY (cell)
123 1119 : NULLIFY (subsys)
124 :
125 2004 : SELECT CASE (gopt_env%type_id)
126 : CASE (default_minimization_method_id, default_ts_method_id)
127 885 : CALL force_env_get(gopt_env%force_env, subsys=subsys)
128 : ! before starting we handle the case of translating coordinates (QM/MM)
129 885 : CALL apply_qmmm_translation(gopt_env)
130 885 : nparticle = force_env_get_nparticle(gopt_env%force_env)
131 2655 : ALLOCATE (x0(3*nparticle))
132 885 : CALL pack_subsys_particles(subsys=subsys, r=x0)
133 : CASE (default_cell_method_id)
134 234 : CALL force_env_get(gopt_env%force_env, subsys=subsys, cell=cell)
135 : ! Store reference cell
136 6084 : gopt_env%h_ref = cell%hmat
137 : ! before starting we handle the case of translating coordinates (QM/MM)
138 234 : CALL apply_qmmm_translation(gopt_env)
139 234 : nparticle = force_env_get_nparticle(gopt_env%force_env)
140 702 : ALLOCATE (x0(3*nparticle + 6))
141 234 : CALL pack_subsys_particles(subsys=subsys, r=x0)
142 234 : idg = 3*nparticle
143 936 : DO i = 1, 3
144 2340 : DO j = 1, i
145 1404 : idg = idg + 1
146 2106 : x0(idg) = gopt_env%cell_env%opt_cell%hmat(j, i)
147 : END DO
148 : END DO
149 : CASE DEFAULT
150 1119 : CPABORT("Invalid or not yet implemented type of optimization")
151 : END SELECT
152 :
153 1119 : END SUBROUTINE gopt_f_create_x0
154 :
155 : ! **************************************************************************************************
156 : !> \brief Apply the QM/MM and QM/MM/X coordinate translations before packing x0
157 : !> \param gopt_env ...
158 : ! **************************************************************************************************
159 1119 : SUBROUTINE apply_qmmm_translation(gopt_env)
160 : TYPE(gopt_f_type), POINTER :: gopt_env
161 :
162 1119 : IF (gopt_env%force_env%in_use == use_qmmm) THEN
163 36 : CALL apply_qmmm_translate(gopt_env%force_env%qmmm_env)
164 : END IF
165 1119 : IF (gopt_env%force_env%in_use == use_qmmmx) THEN
166 0 : CALL apply_qmmmx_translate(gopt_env%force_env%qmmmx_env)
167 : END IF
168 1119 : END SUBROUTINE apply_qmmm_translation
169 :
170 : ! **************************************************************************************************
171 : !> \brief evaluete the potential energy and its gradients using an array
172 : !> with same dimension as the particle_set
173 : !> \param gopt_env the geometry optimization environment
174 : !> \param x the position where the function should be evaluated
175 : !> \param f the function value
176 : !> \param gradient the value of its gradient
177 : !> \param master ...
178 : !> \param para_env ...
179 : !> \par History
180 : !> CELL OPTIMIZATION: Teodoro Laino [tlaino] - University of Zurich - 03.2008
181 : !> 07.2020 Pierre Cazade [pcazade] Space Group Symmetry
182 : !> \author Teodoro Laino [tlaino] - University of Zurich - 01.2008
183 : ! **************************************************************************************************
184 17045 : SUBROUTINE cp_eval_at(gopt_env, x, f, gradient, master, &
185 : para_env)
186 :
187 : TYPE(gopt_f_type), POINTER :: gopt_env
188 : REAL(KIND=dp), DIMENSION(:), POINTER :: x
189 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: f
190 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: gradient
191 : INTEGER, INTENT(IN) :: master
192 : TYPE(mp_para_env_type), POINTER :: para_env
193 :
194 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_eval_at'
195 :
196 : INTEGER :: handle, idg, idir, ip, nparticle, nsize, &
197 : shell_index
198 : REAL(KIND=dp) :: f_ts
199 : REAL(KIND=dp), DIMENSION(3, 3) :: av_ptens, av_ptens_opt
200 17045 : REAL(KIND=dp), DIMENSION(:), POINTER :: cell_gradient, gradient_ts
201 : TYPE(cell_type), POINTER :: cell
202 : TYPE(cp_subsys_type), POINTER :: subsys
203 : TYPE(particle_list_type), POINTER :: core_particles, particles, &
204 : shell_particles
205 : TYPE(spgr_type), POINTER :: spgr
206 : TYPE(virial_type), POINTER :: virial
207 :
208 17045 : NULLIFY (cell)
209 17045 : NULLIFY (core_particles)
210 17045 : NULLIFY (gradient_ts)
211 17045 : NULLIFY (particles)
212 17045 : NULLIFY (shell_particles)
213 17045 : NULLIFY (subsys)
214 17045 : NULLIFY (virial)
215 17045 : NULLIFY (spgr)
216 :
217 17045 : CALL timeset(routineN, handle)
218 :
219 17045 : CALL force_env_get(gopt_env%force_env, subsys=subsys, cell=cell)
220 : CALL cp_subsys_get(subsys, &
221 : core_particles=core_particles, &
222 : particles=particles, &
223 : shell_particles=shell_particles, &
224 17045 : virial=virial)
225 :
226 17045 : spgr => gopt_env%spgr
227 :
228 25450 : SELECT CASE (gopt_env%type_id)
229 : CASE (default_minimization_method_id, default_ts_method_id)
230 8405 : CALL unpack_subsys_particles(subsys=subsys, r=x)
231 8405 : CALL write_structure_data(particles%els, cell, gopt_env%motion_section)
232 15022 : SELECT CASE (gopt_env%type_id)
233 : CASE (default_minimization_method_id)
234 : ! Geometry Minimization
235 : CALL force_env_calc_energy_force(gopt_env%force_env, &
236 : calc_force=PRESENT(gradient), &
237 6617 : require_consistent_energy_force=gopt_env%require_consistent_energy_force)
238 : ! Possibly take the potential energy
239 6617 : IF (PRESENT(f)) THEN
240 6617 : CALL force_env_get(gopt_env%force_env, potential_energy=f)
241 : END IF
242 : ! Possibly take the gradients
243 6617 : IF (PRESENT(gradient)) THEN
244 6079 : IF (master == para_env%mepos) THEN ! we are on the master
245 5927 : CALL pack_subsys_particles(subsys=subsys, f=gradient, fscale=-1.0_dp)
246 5927 : IF (spgr%keep_space_group) THEN
247 8 : CALL spgr_apply_rotations_force(spgr, gradient)
248 8 : CALL unpack_subsys_particles(subsys=subsys, f=gradient, fscale=-1.0_dp)
249 : END IF
250 : END IF
251 : END IF
252 : CASE (default_ts_method_id)
253 : ! Transition State Optimization
254 5364 : ALLOCATE (gradient_ts(particles%n_els*3))
255 : ! Real calculation of energy and forces for transition state optimization:
256 : ! When doing dimer methods forces have to be always computed since the function
257 : ! to minimize is not the energy but the effective force
258 1788 : CALL cp_eval_at_ts(gopt_env, x, f_ts, gradient_ts, calc_force=.TRUE.)
259 1788 : CALL cite_reference(Henkelman1999)
260 : ! Possibly take the potential energy
261 1788 : IF (PRESENT(f)) f = f_ts
262 : ! Possibly take the gradients
263 1788 : IF (PRESENT(gradient)) THEN
264 844 : IF (master == para_env%mepos) THEN ! we are on the master
265 844 : CPASSERT(ASSOCIATED(gradient))
266 29420 : gradient = gradient_ts
267 : END IF
268 : END IF
269 10193 : DEALLOCATE (gradient_ts)
270 : END SELECT
271 : ! This call is necessary for QM/MM if a Translation is applied
272 : ! this makes the geometry optimizer consistent
273 8405 : CALL unpack_subsys_particles(subsys=subsys, r=x)
274 : CASE (default_cell_method_id)
275 : ! Check for VIRIAL
276 8320 : IF (.NOT. virial%pv_availability) THEN
277 : CALL cp_abort(__LOCATION__, &
278 : "For the CELL_OPT task, the FORCE_EVAL/STRESS_TENSOR "// &
279 : "keyword MUST be defined in the input file for the "// &
280 0 : "evaluation of the stress tensor, but none is found!")
281 : END IF
282 8320 : IF (gopt_env%cell_env%keep_volume) THEN
283 2284 : nparticle = force_env_get_nparticle(gopt_env%force_env)
284 2284 : idg = 3*nparticle
285 2284 : CALL rescale_new_cell_volume(cell%deth, x, idg)
286 : END IF
287 :
288 8320 : CALL apply_cell_change(gopt_env, cell, x, update_forces=.FALSE.)
289 : ! Possibly output the new cell used for the next calculation
290 8320 : CALL write_cell(cell, gopt_env%geo_section)
291 : ! Compute the pressure tensor
292 8320 : BLOCK
293 : TYPE(virial_type) :: virial_avg
294 : CALL force_env_calc_energy_force(gopt_env%force_env, &
295 : calc_force=PRESENT(gradient), &
296 8320 : require_consistent_energy_force=gopt_env%require_consistent_energy_force)
297 : ! Possibly take the potential energy
298 8320 : virial_avg = virial
299 8320 : CALL virial_update(virial_avg, subsys, para_env)
300 8320 : IF (PRESENT(f)) THEN
301 8320 : CALL force_env_get(gopt_env%force_env, potential_energy=f)
302 : END IF
303 : ! Possibly take the gradients
304 2063360 : IF (PRESENT(gradient)) THEN
305 6992 : CPASSERT(ANY(virial_avg%pv_total /= 0))
306 : ! Convert the average ptens
307 90896 : av_ptens(:, :) = virial_avg%pv_total(:, :)/cell%deth
308 6992 : IF (master == para_env%mepos) THEN ! we are on the master
309 5333 : CPASSERT(ASSOCIATED(gradient))
310 5333 : nparticle = force_env_get_nparticle(gopt_env%force_env)
311 5333 : nsize = 3*nparticle
312 5333 : CPASSERT((SIZE(gradient) == nsize + 6))
313 5333 : CALL pack_subsys_particles(subsys=subsys, f=gradient(1:nsize), fscale=-1.0_dp)
314 5333 : CALL apply_cell_change(gopt_env, cell, gradient, update_forces=.TRUE.)
315 5333 : IF (spgr%keep_space_group) THEN
316 556 : CALL spgr_apply_rotations_force(spgr, gradient)
317 556 : CALL spgr_apply_rotations_stress(spgr, cell, av_ptens)
318 556 : CALL spgr_write_stress_tensor(av_ptens, spgr)
319 : END IF
320 5333 : cell_gradient => gradient(nsize + 1:nsize + 6)
321 37331 : cell_gradient = 0.0_dp
322 : av_ptens_opt = MATMUL(gopt_env%cell_env%input_to_opt, &
323 421307 : MATMUL(av_ptens, gopt_env%cell_env%opt_to_input))
324 : CALL get_dg_dh(cell_gradient, av_ptens_opt, gopt_env%cell_env%pres_ext, &
325 : gopt_env%cell_env%opt_cell, gopt_env%cell_env%mtrx, &
326 : keep_angles=gopt_env%cell_env%keep_angles, &
327 : keep_symmetry=gopt_env%cell_env%keep_symmetry, &
328 : pres_int=gopt_env%cell_env%pres_int, &
329 : pres_constr=gopt_env%cell_env%pres_constr, &
330 5333 : constraint_id=gopt_env%cell_env%constraint_id)
331 : END IF
332 : ! some callers expect pres_int to be available on all ranks. Also, here master is not necessarily a single rank.
333 : ! Assume at least master==0
334 6992 : CALL para_env%bcast(gopt_env%cell_env%pres_int, 0)
335 6992 : IF (gopt_env%cell_env%constraint_id /= fix_none) THEN
336 296 : CALL para_env%bcast(gopt_env%cell_env%pres_constr, 0)
337 : END IF
338 : END IF
339 : END BLOCK
340 : CASE (default_shellcore_method_id)
341 : idg = 0
342 31040 : DO ip = 1, particles%n_els
343 30720 : shell_index = particles%els(ip)%shell_index
344 31040 : IF (shell_index /= 0) THEN
345 122880 : DO idir = 1, 3
346 92160 : idg = 3*(shell_index - 1) + idir
347 122880 : shell_particles%els(shell_index)%r(idir) = core_particles%els(ip)%r(idir) - x(idg)
348 : END DO
349 : END IF
350 : END DO
351 320 : CALL write_structure_data(particles%els, cell, gopt_env%motion_section)
352 :
353 : ! Shell-core optimization
354 : CALL force_env_calc_energy_force(gopt_env%force_env, &
355 : calc_force=PRESENT(gradient), &
356 320 : require_consistent_energy_force=gopt_env%require_consistent_energy_force)
357 :
358 : ! Possibly take the potential energy
359 320 : IF (PRESENT(f)) THEN
360 320 : CALL force_env_get(gopt_env%force_env, potential_energy=f)
361 : END IF
362 :
363 : ! Possibly take the gradients
364 320 : IF (PRESENT(gradient)) THEN
365 320 : IF (master == para_env%mepos) THEN ! we are on the master
366 320 : CPASSERT(ASSOCIATED(gradient))
367 320 : idg = 0
368 31040 : DO ip = 1, shell_particles%n_els
369 123200 : DO idir = 1, 3
370 92160 : idg = idg + 1
371 122880 : gradient(idg) = -(core_particles%els(ip)%f(idir) - shell_particles%els(ip)%f(idir))
372 : END DO
373 : END DO
374 : END IF
375 : END IF
376 : CASE DEFAULT
377 17045 : CPABORT("Invalid or not yet implemented type of optimization")
378 : END SELECT
379 :
380 17045 : CALL timestop(handle)
381 :
382 17045 : END SUBROUTINE cp_eval_at
383 :
384 : ! **************************************************************************************************
385 : !> \brief Prints iteration step of the optimization procedure on screen
386 : !> \param its ...
387 : !> \param output_unit ...
388 : !> \author Teodoro Laino [tlaino] - University of Zurich - 03.2008
389 : ! **************************************************************************************************
390 9274 : SUBROUTINE gopt_f_ii(its, output_unit)
391 :
392 : INTEGER, INTENT(IN) :: its, output_unit
393 :
394 9274 : IF (output_unit > 0) THEN
395 4796 : WRITE (UNIT=output_unit, FMT="(/,T2,26('-'))")
396 4796 : WRITE (UNIT=output_unit, FMT="(T2,A,I6)") "OPTIMIZATION STEP: ", its
397 4796 : WRITE (UNIT=output_unit, FMT="(T2,26('-'))")
398 4796 : CALL m_flush(output_unit)
399 : END IF
400 :
401 9274 : END SUBROUTINE gopt_f_ii
402 :
403 : ! **************************************************************************************************
404 : !> \brief Retrieves the reporting metadata used during an optimization step.
405 : !> \param gopt_env Optimization environment supplying print configuration and units
406 : !> \param force_env Force environment used for memory sampling
407 : !> \param max_memory Sampled memory, or zero if memory reporting is disabled
408 : !> \param energy_unit Energy reporting unit
409 : !> \param stress_unit Stress reporting unit
410 : ! **************************************************************************************************
411 10435 : SUBROUTINE get_reporting_metadata(gopt_env, force_env, max_memory, energy_unit, stress_unit)
412 :
413 : TYPE(gopt_f_type), POINTER :: gopt_env
414 : TYPE(force_env_type), POINTER :: force_env
415 : INTEGER(KIND=int_8), INTENT(OUT) :: max_memory
416 : CHARACTER(LEN=default_string_length), INTENT(OUT) :: energy_unit, stress_unit
417 :
418 : LOGICAL :: print_memory
419 : TYPE(mp_para_env_type), POINTER :: para_env
420 :
421 10435 : NULLIFY (para_env)
422 10435 : CALL section_vals_val_get(gopt_env%motion_section, "PRINT%MEMORY_INFO", l_val=print_memory)
423 10435 : max_memory = 0
424 10435 : IF (print_memory) THEN
425 10435 : CALL force_env_get(force_env, para_env=para_env)
426 10435 : max_memory = sample_memory(para_env)
427 : END IF
428 :
429 : CALL section_vals_val_get(gopt_env%force_env%force_env_section, &
430 : "PRINT%PROGRAM_RUN_INFO%ENERGY_UNIT", &
431 10435 : c_val=energy_unit)
432 : CALL section_vals_val_get(gopt_env%force_env%force_env_section, &
433 : "PRINT%STRESS_TENSOR%STRESS_UNIT", &
434 10435 : c_val=stress_unit)
435 :
436 10435 : END SUBROUTINE get_reporting_metadata
437 :
438 : ! **************************************************************************************************
439 : !> \brief Handles the Output during an optimization run
440 : !> \param gopt_env ...
441 : !> \param output_unit ...
442 : !> \param opt_energy ...
443 : !> \param wildcard ...
444 : !> \param its ...
445 : !> \param used_time ...
446 : !> \author Teodoro Laino [tlaino] - University of Zurich - 03.2008
447 : ! **************************************************************************************************
448 1163 : SUBROUTINE gopt_f_io_init(gopt_env, output_unit, opt_energy, wildcard, its, used_time)
449 :
450 : TYPE(gopt_f_type), POINTER :: gopt_env
451 : INTEGER, INTENT(IN) :: output_unit
452 : REAL(KIND=dp) :: opt_energy
453 : CHARACTER(LEN=5) :: wildcard
454 : INTEGER, INTENT(IN) :: its
455 : REAL(KIND=dp) :: used_time
456 :
457 : CHARACTER(LEN=default_string_length) :: energy_unit, stress_unit
458 : INTEGER(KIND=int_8) :: max_memory
459 : REAL(KIND=dp) :: pres_int
460 :
461 1163 : CALL get_reporting_metadata(gopt_env, gopt_env%force_env, max_memory, energy_unit, stress_unit)
462 :
463 2130 : SELECT CASE (gopt_env%type_id)
464 : CASE (default_ts_method_id, default_minimization_method_id)
465 : ! Geometry Optimization (Minimization and Transition State Search)
466 967 : IF (.NOT. gopt_env%dimer_rotation) THEN
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 839 : stress_unit=stress_unit)
475 : ELSE
476 : CALL write_rot_cycle_infos(output_unit, &
477 : it=its, &
478 : etot=opt_energy, &
479 : dimer_env=gopt_env%dimer_env, &
480 : wildcard=wildcard, &
481 : used_time=used_time, &
482 128 : max_memory=max_memory)
483 : END IF
484 : CASE (default_cell_method_id)
485 : ! Cell Optimization
486 176 : pres_int = gopt_env%cell_env%pres_int
487 : CALL write_cycle_infos(output_unit, &
488 : it=its, &
489 : etot=opt_energy, &
490 : pres_int=pres_int, &
491 : wildcard=wildcard, &
492 : used_time=used_time, &
493 : max_memory=max_memory, &
494 : energy_unit=energy_unit, &
495 176 : stress_unit=stress_unit)
496 : CASE (default_shellcore_method_id)
497 : CALL write_cycle_infos(output_unit, &
498 : it=its, &
499 : etot=opt_energy, &
500 : wildcard=wildcard, &
501 : used_time=used_time, &
502 : max_memory=max_memory, &
503 : energy_unit=energy_unit, &
504 1163 : stress_unit=stress_unit)
505 : END SELECT
506 :
507 1163 : END SUBROUTINE gopt_f_io_init
508 :
509 : ! **************************************************************************************************
510 : !> \brief Handles the Output during an optimization run
511 : !> \param gopt_env ...
512 : !> \param force_env ...
513 : !> \param root_section ...
514 : !> \param its ...
515 : !> \param opt_energy ...
516 : !> \param output_unit ...
517 : !> \param eold ...
518 : !> \param emin ...
519 : !> \param wildcard ...
520 : !> \param gopt_param ...
521 : !> \param ndf ...
522 : !> \param dx ...
523 : !> \param xi ...
524 : !> \param conv ...
525 : !> \param pred ...
526 : !> \param rat ...
527 : !> \param step ...
528 : !> \param rad ...
529 : !> \param used_time ...
530 : !> \author Teodoro Laino [tlaino] - University of Zurich - 03.2008
531 : ! **************************************************************************************************
532 18544 : SUBROUTINE gopt_f_io(gopt_env, force_env, root_section, its, opt_energy, &
533 9272 : output_unit, eold, emin, wildcard, gopt_param, ndf, dx, xi, conv, pred, rat, &
534 : step, rad, used_time)
535 :
536 : TYPE(gopt_f_type), POINTER :: gopt_env
537 : TYPE(force_env_type), POINTER :: force_env
538 : TYPE(section_vals_type), POINTER :: root_section
539 : INTEGER, INTENT(IN) :: its
540 : REAL(KIND=dp), INTENT(IN) :: opt_energy
541 : INTEGER, INTENT(IN) :: output_unit
542 : REAL(KIND=dp) :: eold, emin
543 : CHARACTER(LEN=5) :: wildcard
544 : TYPE(gopt_param_type), POINTER :: gopt_param
545 : INTEGER, INTENT(IN), OPTIONAL :: ndf
546 : REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: dx
547 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: xi
548 : LOGICAL, OPTIONAL :: conv
549 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: pred, rat, step, rad
550 : REAL(KIND=dp) :: used_time
551 :
552 : CHARACTER(LEN=default_string_length) :: energy_unit, stress_unit
553 : INTEGER(KIND=int_8) :: max_memory
554 : REAL(KIND=dp) :: pres_diff, pres_diff_constr, pres_int, &
555 : pres_tol
556 :
557 9272 : CALL get_reporting_metadata(gopt_env, force_env, max_memory, energy_unit, stress_unit)
558 :
559 14350 : SELECT CASE (gopt_env%type_id)
560 : CASE (default_ts_method_id, default_minimization_method_id)
561 : ! Geometry Optimization (Minimization and Transition State Search)
562 5078 : IF (.NOT. gopt_env%dimer_rotation) THEN
563 : CALL geo_opt_io(force_env=force_env, root_section=root_section, &
564 4356 : motion_section=gopt_env%motion_section, its=its, opt_energy=opt_energy)
565 : CALL write_cycle_infos(output_unit, &
566 : it=its, &
567 : etot=opt_energy, &
568 : ediff=(opt_energy - eold), &
569 : pred=pred, &
570 : rat=rat, &
571 : step=step, &
572 : rad=rad, &
573 : emin=emin, &
574 : wildcard=wildcard, &
575 : used_time=used_time, &
576 : max_memory=max_memory, &
577 : energy_unit=energy_unit, &
578 4356 : stress_unit=stress_unit)
579 : ! Possibly check convergence
580 4356 : IF (PRESENT(conv)) THEN
581 4356 : CPASSERT(PRESENT(ndf))
582 4356 : CPASSERT(PRESENT(dx))
583 4356 : CPASSERT(PRESENT(xi))
584 : CALL check_and_report_convergence(output_unit, gopt_param, ndf, dx, xi, conv, &
585 4356 : max_memory, stress_unit)
586 : END IF
587 : ELSE
588 722 : CALL update_dimer_vec(gopt_env%dimer_env, gopt_env%motion_section)
589 722 : CALL write_restart(force_env=force_env, root_section=root_section)
590 : CALL write_rot_cycle_infos(output_unit, its, opt_energy, opt_energy - eold, emin, gopt_env%dimer_env, &
591 722 : wildcard=wildcard, used_time=used_time, max_memory=max_memory)
592 : ! Possibly check convergence
593 722 : IF (PRESENT(conv)) THEN
594 722 : CPASSERT(ASSOCIATED(gopt_env%dimer_env))
595 722 : CALL check_rot_conv(gopt_env%dimer_env, output_unit, conv)
596 : END IF
597 : END IF
598 : CASE (default_cell_method_id)
599 : ! Cell Optimization
600 4044 : pres_diff = gopt_env%cell_env%pres_int - gopt_env%cell_env%pres_ext
601 4044 : pres_int = gopt_env%cell_env%pres_int
602 4044 : pres_tol = gopt_env%cell_env%pres_tol
603 : CALL geo_opt_io(force_env=force_env, root_section=root_section, &
604 4044 : motion_section=gopt_env%motion_section, its=its, opt_energy=opt_energy)
605 : CALL write_cycle_infos(output_unit, &
606 : it=its, &
607 : etot=opt_energy, &
608 : ediff=(opt_energy - eold), &
609 : pred=pred, &
610 : rat=rat, &
611 : step=step, &
612 : rad=rad, &
613 : emin=emin, &
614 : pres_int=pres_int, &
615 : wildcard=wildcard, &
616 : used_time=used_time, &
617 : max_memory=max_memory, &
618 : energy_unit=energy_unit, &
619 4044 : stress_unit=stress_unit)
620 : ! Possibly check convergence
621 4044 : IF (PRESENT(conv)) THEN
622 4044 : CPASSERT(PRESENT(ndf))
623 4044 : CPASSERT(PRESENT(dx))
624 4044 : CPASSERT(PRESENT(xi))
625 4044 : IF (gopt_env%cell_env%constraint_id == fix_none) THEN
626 : CALL check_and_report_convergence(output_unit, gopt_param, ndf, dx, xi, conv, &
627 3920 : max_memory, stress_unit, pres_diff, pres_tol)
628 : ELSE
629 124 : pres_diff_constr = gopt_env%cell_env%pres_constr - gopt_env%cell_env%pres_ext
630 : CALL check_and_report_convergence(output_unit, gopt_param, ndf, dx, xi, conv, &
631 : max_memory, stress_unit, pres_diff, pres_tol, &
632 124 : pres_diff_constr)
633 : END IF
634 : END IF
635 : CASE (default_shellcore_method_id)
636 : CALL write_cycle_infos(output_unit, &
637 : it=its, &
638 : etot=opt_energy, &
639 : ediff=(opt_energy - eold), &
640 : pred=pred, &
641 : rat=rat, &
642 : step=step, &
643 : rad=rad, &
644 : emin=emin, &
645 : wildcard=wildcard, &
646 : used_time=used_time, &
647 : max_memory=max_memory, &
648 : energy_unit=energy_unit, &
649 150 : stress_unit=stress_unit)
650 : ! Possibly check convergence
651 9422 : IF (PRESENT(conv)) THEN
652 150 : CPASSERT(PRESENT(ndf))
653 150 : CPASSERT(PRESENT(dx))
654 150 : CPASSERT(PRESENT(xi))
655 : CALL check_and_report_convergence(output_unit, gopt_param, ndf, dx, xi, conv, &
656 150 : max_memory, stress_unit)
657 : END IF
658 : END SELECT
659 :
660 9272 : END SUBROUTINE gopt_f_io
661 :
662 : ! **************************************************************************************************
663 : !> \brief Handles the Output at the end of an optimization run
664 : !> \param gopt_env ...
665 : !> \param force_env ...
666 : !> \param x0 ...
667 : !> \param conv ...
668 : !> \param its ...
669 : !> \param root_section ...
670 : !> \param para_env ...
671 : !> \param master ...
672 : !> \param output_unit ...
673 : !> \author Teodoro Laino [tlaino] - University of Zurich - 03.2008
674 : ! **************************************************************************************************
675 1267 : RECURSIVE SUBROUTINE gopt_f_io_finalize(gopt_env, force_env, x0, conv, its, root_section, &
676 : para_env, master, output_unit)
677 : TYPE(gopt_f_type), POINTER :: gopt_env
678 : TYPE(force_env_type), POINTER :: force_env
679 : REAL(KIND=dp), DIMENSION(:), POINTER :: x0
680 : LOGICAL :: conv
681 : INTEGER :: its
682 : TYPE(section_vals_type), POINTER :: root_section
683 : TYPE(mp_para_env_type), POINTER :: para_env
684 : INTEGER, INTENT(IN) :: master, output_unit
685 :
686 1267 : IF (gopt_env%eval_opt_geo) THEN
687 1247 : IF (.NOT. gopt_env%dimer_rotation) THEN
688 : CALL write_final_info(output_unit, conv, its, gopt_env, x0, master, &
689 1119 : para_env, force_env, gopt_env%motion_section, root_section)
690 : ELSE
691 128 : CALL update_dimer_vec(gopt_env%dimer_env, gopt_env%motion_section)
692 128 : CALL write_restart(force_env=force_env, root_section=root_section)
693 : END IF
694 : END IF
695 :
696 1267 : END SUBROUTINE gopt_f_io_finalize
697 :
698 : ! **************************************************************************************************
699 : !> \brief ...
700 : !> \param output_unit ...
701 : !> \param it ...
702 : !> \param etot ...
703 : !> \param ediff ...
704 : !> \param pred ...
705 : !> \param rat ...
706 : !> \param step ...
707 : !> \param rad ...
708 : !> \param emin ...
709 : !> \param pres_int ...
710 : !> \param wildcard ...
711 : !> \param used_time ...
712 : !> \param max_memory ...
713 : !> \param energy_unit ...
714 : !> \param stress_unit ...
715 : ! **************************************************************************************************
716 9585 : SUBROUTINE write_cycle_infos(output_unit, it, etot, ediff, pred, rat, step, rad, emin, &
717 : pres_int, wildcard, used_time, max_memory, energy_unit, stress_unit)
718 :
719 : INTEGER, INTENT(IN) :: output_unit, it
720 : REAL(KIND=dp), INTENT(IN) :: etot
721 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: ediff, pred, rat, step, rad, emin, &
722 : pres_int
723 : CHARACTER(LEN=5), INTENT(IN) :: wildcard
724 : REAL(KIND=dp), INTENT(IN) :: used_time
725 : INTEGER(KIND=int_8), INTENT(IN) :: max_memory
726 : CHARACTER(LEN=default_string_length), INTENT(IN) :: energy_unit, stress_unit
727 :
728 : CHARACTER(LEN=5) :: tag
729 :
730 9585 : IF (output_unit > 0) THEN
731 4975 : tag = "OPT| "
732 4975 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") tag//REPEAT("*", 74)
733 : WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,I25)") &
734 4975 : tag//"Step number", it
735 : WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,A25)") &
736 4975 : tag//"Optimization method", wildcard
737 : WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
738 4975 : tag//"Total energy ["//TRIM(ADJUSTL(energy_unit))//"]", &
739 9950 : cp_unit_from_cp2k(etot, TRIM(energy_unit))
740 4975 : IF (PRESENT(pres_int)) THEN
741 : WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
742 2110 : tag//"Internal pressure ["//TRIM(ADJUSTL(stress_unit))//"]", &
743 4220 : cp_unit_from_cp2k(pres_int, TRIM(stress_unit))
744 : END IF
745 4975 : IF (PRESENT(ediff)) THEN
746 : WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
747 4434 : tag//"Effective energy change ["//TRIM(ADJUSTL(energy_unit))//"]", &
748 8868 : cp_unit_from_cp2k(ediff, TRIM(energy_unit))
749 : END IF
750 4975 : IF (PRESENT(pred)) THEN
751 : WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
752 2094 : tag//"Predicted energy change ["//TRIM(ADJUSTL(energy_unit))//"]", &
753 4188 : cp_unit_from_cp2k(pred, TRIM(energy_unit))
754 : END IF
755 4975 : IF (PRESENT(rat)) THEN
756 : WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
757 2350 : tag//"Scaling factor", rat
758 : END IF
759 4975 : IF (PRESENT(step)) THEN
760 : WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
761 2350 : tag//"Step size", step
762 : END IF
763 4975 : IF (PRESENT(rad)) THEN
764 : WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
765 2350 : tag//"Trust radius", rad
766 : END IF
767 4975 : IF (PRESENT(emin)) THEN
768 4434 : IF (etot < emin) THEN
769 : WRITE (UNIT=output_unit, FMT="(T2,A,T77,A4)") &
770 4026 : tag//"Decrease in energy", " YES"
771 : ELSE
772 : WRITE (UNIT=output_unit, FMT="(T2,A,T77,A4)") &
773 408 : tag//"Decrease in energy", " NO"
774 : END IF
775 : END IF
776 : WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.3)") &
777 4975 : tag//"Used time [s]", used_time
778 4975 : IF (it == 0) THEN
779 535 : WRITE (UNIT=output_unit, FMT="(T2,A)") tag//REPEAT("*", 74)
780 535 : IF (max_memory /= 0) THEN
781 : WRITE (UNIT=output_unit, FMT="(T2,A,T60,1X,I20)") &
782 535 : tag//"Estimated peak process memory [MiB]", &
783 1070 : (max_memory + bytes_per_mib - 1)/bytes_per_mib
784 : END IF
785 : END IF
786 : END IF
787 :
788 9585 : END SUBROUTINE write_cycle_infos
789 :
790 : ! **************************************************************************************************
791 : !> \brief ...
792 : !> \param output_unit ...
793 : !> \param it ...
794 : !> \param etot ...
795 : !> \param ediff ...
796 : !> \param emin ...
797 : !> \param dimer_env ...
798 : !> \param used_time ...
799 : !> \param wildcard ...
800 : !> \param max_memory ...
801 : !> \date 01.2008
802 : !> \author Luca Bellucci and Teodoro Laino - created [tlaino]
803 : ! **************************************************************************************************
804 850 : SUBROUTINE write_rot_cycle_infos(output_unit, it, etot, ediff, emin, dimer_env, used_time, &
805 : wildcard, max_memory)
806 :
807 : INTEGER, INTENT(IN) :: output_unit, it
808 : REAL(KIND=dp), INTENT(IN) :: etot
809 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: ediff, emin
810 : TYPE(dimer_env_type), POINTER :: dimer_env
811 : REAL(KIND=dp), INTENT(IN) :: used_time
812 : CHARACTER(LEN=5), INTENT(IN) :: wildcard
813 : INTEGER(KIND=int_8), INTENT(IN) :: max_memory
814 :
815 : CHARACTER(LEN=5) :: tag
816 :
817 850 : IF (output_unit > 0) THEN
818 425 : tag = "OPT| "
819 425 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") tag//REPEAT("*", 74)
820 : WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,I25)") &
821 425 : tag//"Rotational step number", it
822 : WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,A25)") &
823 425 : tag//"Optimization method", wildcard
824 : WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
825 425 : tag//"Local curvature", dimer_env%rot%curvature, &
826 850 : tag//"Total rotational force", etot
827 425 : IF (PRESENT(ediff)) THEN
828 : WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
829 361 : tag//"Rotational force change", ediff
830 : END IF
831 425 : IF (PRESENT(emin)) THEN
832 361 : IF (etot < emin) THEN
833 : WRITE (UNIT=output_unit, FMT="(T2,A,T77,A4)") &
834 155 : tag//"Decrease in rotational force", " YES"
835 : ELSE
836 : WRITE (UNIT=output_unit, FMT="(T2,A,T77,A4)") &
837 206 : tag//"Decrease in rotational force", " NO"
838 : END IF
839 : END IF
840 : WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.3)") &
841 425 : tag//"Used time [s]", used_time
842 425 : IF (it == 0) THEN
843 64 : WRITE (UNIT=output_unit, FMT="(T2,A)") tag//REPEAT("*", 74)
844 64 : IF (max_memory /= 0) THEN
845 : WRITE (UNIT=output_unit, FMT="(T2,A,T60,1X,I20)") &
846 64 : tag//"Estimated peak process memory [MiB]", &
847 128 : (max_memory + bytes_per_mib - 1)/bytes_per_mib
848 : END IF
849 : END IF
850 : END IF
851 :
852 850 : END SUBROUTINE write_rot_cycle_infos
853 :
854 : ! **************************************************************************************************
855 : !> \brief Evaluate the geometry/cell optimization convergence criteria (pure decision, no I/O)
856 : !> \param ndf ...
857 : !> \param dr ...
858 : !> \param g ...
859 : !> \param gopt_param ...
860 : !> \param conv ...
861 : !> \param conv_dx ...
862 : !> \param conv_rdx ...
863 : !> \param conv_g ...
864 : !> \param conv_rg ...
865 : !> \param conv_p ...
866 : !> \param max_dr ...
867 : !> \param rms_dr ...
868 : !> \param max_grad ...
869 : !> \param rms_grad ...
870 : !> \param pres_diff ...
871 : !> \param pres_tol ...
872 : !> \param pres_diff_constr ...
873 : ! **************************************************************************************************
874 8550 : SUBROUTINE check_converg(ndf, dr, g, gopt_param, conv, conv_dx, conv_rdx, conv_g, conv_rg, conv_p, &
875 : max_dr, rms_dr, max_grad, rms_grad, pres_diff, pres_tol, pres_diff_constr)
876 :
877 : INTEGER, INTENT(IN) :: ndf
878 : REAL(KIND=dp), INTENT(IN) :: dr(ndf), g(ndf)
879 : TYPE(gopt_param_type), POINTER :: gopt_param
880 : LOGICAL, INTENT(OUT) :: conv, conv_dx, conv_rdx, conv_g, &
881 : conv_rg, conv_p
882 : REAL(KIND=dp), INTENT(OUT) :: max_dr, rms_dr, max_grad, rms_grad
883 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: pres_diff, pres_tol, pres_diff_constr
884 :
885 : INTEGER :: indf
886 : REAL(KIND=dp) :: dumm, dxcon, gcon, rmsgcon, rmsxcon
887 :
888 8550 : dxcon = gopt_param%max_dr
889 8550 : gcon = gopt_param%max_force
890 8550 : rmsgcon = gopt_param%rms_force
891 8550 : rmsxcon = gopt_param%rms_dr
892 :
893 8550 : conv = .FALSE.
894 8550 : conv_dx = .TRUE.
895 8550 : conv_rdx = .TRUE.
896 8550 : conv_g = .TRUE.
897 8550 : conv_rg = .TRUE.
898 8550 : conv_p = .TRUE.
899 :
900 8550 : dumm = 0.0_dp
901 2545146 : DO indf = 1, ndf
902 2536596 : IF (indf == 1) max_dr = ABS(dr(indf))
903 2536596 : dumm = dumm + dr(indf)**2
904 2536596 : IF (ABS(dr(indf)) > dxcon) conv_dx = .FALSE.
905 2545146 : IF (ABS(dr(indf)) > max_dr) max_dr = ABS(dr(indf))
906 : END DO
907 : ! SQRT(dumm/ndf) > rmsxcon
908 8550 : IF (dumm > (rmsxcon*rmsxcon*ndf)) conv_rdx = .FALSE.
909 8550 : rms_dr = SQRT(dumm/ndf)
910 :
911 8550 : dumm = 0.0_dp
912 2545146 : DO indf = 1, ndf
913 2536596 : IF (indf == 1) max_grad = ABS(g(indf))
914 2536596 : dumm = dumm + g(indf)**2
915 2536596 : IF (ABS(g(indf)) > gcon) conv_g = .FALSE.
916 2545146 : IF (ABS(g(indf)) > max_grad) max_grad = ABS(g(indf))
917 : END DO
918 : ! SQRT(dumm/ndf) > rmsgcon
919 8550 : IF (dumm > (rmsgcon*rmsgcon*ndf)) conv_rg = .FALSE.
920 8550 : rms_grad = SQRT(dumm/ndf)
921 :
922 8550 : IF (PRESENT(pres_diff_constr) .AND. PRESENT(pres_tol)) THEN
923 124 : conv_p = ABS(pres_diff_constr) < ABS(pres_tol)
924 8426 : ELSE IF (PRESENT(pres_diff) .AND. PRESENT(pres_tol)) THEN
925 3920 : conv_p = ABS(pres_diff) < ABS(pres_tol)
926 : END IF
927 :
928 8550 : IF (conv_dx .AND. conv_rdx .AND. conv_g .AND. conv_rg .AND. conv_p) conv = .TRUE.
929 :
930 8550 : END SUBROUTINE check_converg
931 :
932 : ! **************************************************************************************************
933 : !> \brief Print one YES/NO convergence flag line
934 : !> \param output_unit ...
935 : !> \param tag ...
936 : !> \param label ...
937 : !> \param flag ...
938 : ! **************************************************************************************************
939 20119 : SUBROUTINE write_conv_flag(output_unit, tag, label, flag)
940 : INTEGER, INTENT(IN) :: output_unit
941 : CHARACTER(LEN=*), INTENT(IN) :: tag, label
942 : LOGICAL, INTENT(IN) :: flag
943 :
944 20119 : IF (flag) THEN
945 3463 : WRITE (UNIT=output_unit, FMT="(T2,A,T77,A4)") tag//label, " YES"
946 : ELSE
947 16656 : WRITE (UNIT=output_unit, FMT="(T2,A,T77,A4)") tag//label, " NO"
948 : END IF
949 20119 : END SUBROUTINE write_conv_flag
950 :
951 : ! **************************************************************************************************
952 : !> \brief Print the per-cycle convergence report (body moved verbatim from check_converg)
953 : !> \param output_unit ...
954 : !> \param gopt_param ...
955 : !> \param conv ...
956 : !> \param conv_dx ...
957 : !> \param conv_rdx ...
958 : !> \param conv_g ...
959 : !> \param conv_rg ...
960 : !> \param conv_p ...
961 : !> \param max_dr ...
962 : !> \param rms_dr ...
963 : !> \param max_grad ...
964 : !> \param rms_grad ...
965 : !> \param max_memory ...
966 : !> \param stress_unit ...
967 : !> \param pres_diff ...
968 : !> \param pres_tol ...
969 : !> \param pres_diff_constr ...
970 : ! **************************************************************************************************
971 8550 : SUBROUTINE write_convergence_infos(output_unit, gopt_param, conv, conv_dx, conv_rdx, conv_g, &
972 : conv_rg, conv_p, max_dr, rms_dr, max_grad, rms_grad, &
973 : max_memory, stress_unit, pres_diff, pres_tol, pres_diff_constr)
974 :
975 : INTEGER, INTENT(IN) :: output_unit
976 : TYPE(gopt_param_type), POINTER :: gopt_param
977 : LOGICAL, INTENT(IN) :: conv, conv_dx, conv_rdx, conv_g, &
978 : conv_rg, conv_p
979 : REAL(KIND=dp), INTENT(IN) :: max_dr, rms_dr, max_grad, rms_grad
980 : INTEGER(KIND=int_8), INTENT(IN) :: max_memory
981 : CHARACTER(LEN=default_string_length), INTENT(IN) :: stress_unit
982 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: pres_diff, pres_tol, pres_diff_constr
983 :
984 : CHARACTER(LEN=5) :: tag
985 : REAL(KIND=dp) :: dxcon, gcon, rmsgcon, rmsxcon
986 :
987 8550 : dxcon = gopt_param%max_dr
988 8550 : gcon = gopt_param%max_force
989 8550 : rmsgcon = gopt_param%rms_force
990 8550 : rmsxcon = gopt_param%rms_dr
991 :
992 8550 : IF (output_unit > 0) THEN
993 :
994 4434 : tag = "OPT| "
995 :
996 4434 : WRITE (UNIT=output_unit, FMT="(T2,A)") TRIM(tag)
997 : WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
998 4434 : tag//"Maximum step size", max_dr, &
999 8868 : tag//"Convergence limit for maximum step size", dxcon
1000 : CALL write_conv_flag(output_unit, tag, "Maximum step size is converged", &
1001 4434 : conv_dx)
1002 :
1003 4434 : WRITE (UNIT=output_unit, FMT="(T2,A)") TRIM(tag)
1004 : WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
1005 4434 : tag//"RMS step size", rms_dr, &
1006 8868 : tag//"Convergence limit for RMS step size", rmsxcon
1007 : CALL write_conv_flag(output_unit, tag, "RMS step size is converged", &
1008 4434 : conv_rdx)
1009 :
1010 4434 : WRITE (UNIT=output_unit, FMT="(T2,A)") TRIM(tag)
1011 : WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
1012 4434 : tag//"Maximum gradient", max_grad, &
1013 8868 : tag//"Convergence limit for maximum gradient", gcon
1014 : CALL write_conv_flag(output_unit, tag, "Maximum gradient is converged", &
1015 4434 : conv_g)
1016 :
1017 4434 : WRITE (UNIT=output_unit, FMT="(T2,A)") TRIM(tag)
1018 : WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
1019 4434 : tag//"RMS gradient", rms_grad, &
1020 8868 : tag//"Convergence limit for RMS gradient", rmsgcon
1021 : CALL write_conv_flag(output_unit, tag, "RMS gradient is converged", &
1022 4434 : conv_rg)
1023 :
1024 4434 : IF (PRESENT(pres_diff) .AND. PRESENT(pres_tol)) THEN
1025 2022 : WRITE (UNIT=output_unit, FMT="(T2,A)") TRIM(tag)
1026 2022 : IF (PRESENT(pres_diff_constr)) THEN
1027 : WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
1028 : tag//"Pressure deviation without constraint ["// &
1029 62 : TRIM(ADJUSTL(stress_unit))//"]", &
1030 124 : cp_unit_from_cp2k(pres_diff, TRIM(stress_unit))
1031 : WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
1032 : tag//"Pressure deviation with constraint ["// &
1033 62 : TRIM(ADJUSTL(stress_unit))//"]", &
1034 124 : cp_unit_from_cp2k(pres_diff_constr, TRIM(stress_unit))
1035 : ELSE
1036 : WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
1037 1960 : tag//"Pressure deviation ["//TRIM(ADJUSTL(stress_unit))//"]", &
1038 3920 : cp_unit_from_cp2k(pres_diff, TRIM(stress_unit))
1039 : END IF
1040 : WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
1041 2022 : tag//"Pressure tolerance ["//TRIM(ADJUSTL(stress_unit))//"]", &
1042 4044 : cp_unit_from_cp2k(pres_tol, TRIM(stress_unit))
1043 2022 : CALL write_conv_flag(output_unit, tag, "Pressure is converged", conv_p)
1044 : END IF
1045 :
1046 4434 : WRITE (UNIT=output_unit, FMT="(T2,A)") tag//REPEAT("*", 74)
1047 :
1048 4434 : IF (max_memory /= 0) THEN
1049 : WRITE (UNIT=output_unit, FMT="(T2,A,T60,1X,I20)") &
1050 4434 : tag//"Estimated peak process memory after this step [MiB]", &
1051 8868 : (max_memory + bytes_per_mib - 1)/bytes_per_mib
1052 : END IF
1053 :
1054 : END IF
1055 :
1056 8550 : IF ((conv) .AND. (output_unit > 0)) THEN
1057 223 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") REPEAT("*", 79)
1058 : WRITE (UNIT=output_unit, FMT="(T2,A,T25,A,T78,A)") &
1059 223 : "***", "GEOMETRY OPTIMIZATION COMPLETED", "***"
1060 223 : WRITE (UNIT=output_unit, FMT="(T2,A)") REPEAT("*", 79)
1061 : END IF
1062 :
1063 8550 : END SUBROUTINE write_convergence_infos
1064 :
1065 : ! **************************************************************************************************
1066 : !> \brief Check convergence and print the per-cycle report (gopt_f_io call sites)
1067 : !> \param output_unit ...
1068 : !> \param gopt_param ...
1069 : !> \param ndf ...
1070 : !> \param dr ...
1071 : !> \param g ...
1072 : !> \param conv ...
1073 : !> \param max_memory ...
1074 : !> \param stress_unit ...
1075 : !> \param pres_diff ...
1076 : !> \param pres_tol ...
1077 : !> \param pres_diff_constr ...
1078 : ! **************************************************************************************************
1079 8550 : SUBROUTINE check_and_report_convergence(output_unit, gopt_param, ndf, dr, g, conv, &
1080 : max_memory, stress_unit, pres_diff, pres_tol, pres_diff_constr)
1081 :
1082 : INTEGER, INTENT(IN) :: output_unit
1083 : TYPE(gopt_param_type), POINTER :: gopt_param
1084 : INTEGER, INTENT(IN) :: ndf
1085 : REAL(KIND=dp), INTENT(IN) :: dr(ndf), g(ndf)
1086 : LOGICAL, INTENT(INOUT) :: conv
1087 : INTEGER(KIND=int_8), INTENT(IN) :: max_memory
1088 : CHARACTER(LEN=default_string_length), INTENT(IN) :: stress_unit
1089 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: pres_diff, pres_tol, pres_diff_constr
1090 :
1091 : LOGICAL :: conv_dx, conv_g, conv_p, conv_rdx, &
1092 : conv_rg
1093 : REAL(KIND=dp) :: max_dr, max_grad, rms_dr, rms_grad
1094 :
1095 : CALL check_converg(ndf, dr, g, gopt_param, conv, conv_dx, conv_rdx, conv_g, conv_rg, &
1096 : conv_p, max_dr, rms_dr, max_grad, rms_grad, pres_diff, pres_tol, &
1097 8550 : pres_diff_constr)
1098 : CALL write_convergence_infos(output_unit, gopt_param, conv, conv_dx, conv_rdx, conv_g, &
1099 : conv_rg, conv_p, max_dr, rms_dr, max_grad, rms_grad, &
1100 8550 : max_memory, stress_unit, pres_diff, pres_tol, pres_diff_constr)
1101 8550 : END SUBROUTINE check_and_report_convergence
1102 :
1103 : ! **************************************************************************************************
1104 : !> \brief ...
1105 : !> \param dimer_env ...
1106 : !> \param output_unit ...
1107 : !> \param conv ...
1108 : !> \date 01.2008
1109 : !> \author Luca Bellucci and Teodoro Laino - created [tlaino]
1110 : ! **************************************************************************************************
1111 722 : SUBROUTINE check_rot_conv(dimer_env, output_unit, conv)
1112 :
1113 : TYPE(dimer_env_type), POINTER :: dimer_env
1114 : INTEGER, INTENT(IN) :: output_unit
1115 : LOGICAL, INTENT(OUT) :: conv
1116 :
1117 : CHARACTER(LEN=5) :: tag
1118 :
1119 722 : conv = (ABS(dimer_env%rot%angle2) < dimer_env%rot%angle_tol)
1120 :
1121 722 : IF (output_unit > 0) THEN
1122 361 : tag = "OPT| "
1123 361 : WRITE (UNIT=output_unit, FMT="(T2,A)") TRIM(tag)
1124 : WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
1125 361 : tag//"Predicted angle step size", dimer_env%rot%angle1, &
1126 361 : tag//"Effective angle step size", dimer_env%rot%angle2, &
1127 722 : tag//"Convergence limit for angle step size", dimer_env%rot%angle_tol
1128 361 : CALL write_conv_flag(output_unit, tag, "Angle step size is converged", conv)
1129 361 : WRITE (UNIT=output_unit, FMT="(T2,A)") tag//REPEAT("*", 74)
1130 : END IF
1131 :
1132 722 : IF ((conv) .AND. (output_unit > 0)) THEN
1133 59 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") REPEAT("*", 79)
1134 : WRITE (UNIT=output_unit, FMT="(T2,A,T25,A,T78,A)") &
1135 59 : "***", "ROTATION OPTIMIZATION COMPLETED", "***"
1136 59 : WRITE (UNIT=output_unit, FMT="(T2,A)") REPEAT("*", 79)
1137 : END IF
1138 :
1139 722 : END SUBROUTINE check_rot_conv
1140 :
1141 : ! **************************************************************************************************
1142 : !> \brief ...
1143 : !> \param output_unit ...
1144 : !> \param conv ...
1145 : !> \param it ...
1146 : !> \param gopt_env ...
1147 : !> \param x0 ...
1148 : !> \param master ...
1149 : !> \param para_env ...
1150 : !> \param force_env ...
1151 : !> \param motion_section ...
1152 : !> \param root_section ...
1153 : !> \date 11.2007
1154 : !> \author Teodoro Laino [tlaino] - University of Zurich
1155 : ! **************************************************************************************************
1156 1119 : RECURSIVE SUBROUTINE write_final_info(output_unit, conv, it, gopt_env, x0, master, para_env, force_env, &
1157 : motion_section, root_section)
1158 : INTEGER, INTENT(IN) :: output_unit
1159 : LOGICAL, INTENT(IN) :: conv
1160 : INTEGER, INTENT(INOUT) :: it
1161 : TYPE(gopt_f_type), POINTER :: gopt_env
1162 : REAL(KIND=dp), DIMENSION(:), POINTER :: x0
1163 : INTEGER, INTENT(IN) :: master
1164 : TYPE(mp_para_env_type), POINTER :: para_env
1165 : TYPE(force_env_type), POINTER :: force_env
1166 : TYPE(section_vals_type), POINTER :: motion_section, root_section
1167 :
1168 : CHARACTER(LEN=default_string_length) :: constraint_label
1169 : LOGICAL :: keep_angles, keep_symmetry, keep_volume
1170 : REAL(KIND=dp) :: etot
1171 : TYPE(cell_type), POINTER :: cell
1172 : TYPE(cp_subsys_type), POINTER :: subsys
1173 : TYPE(particle_list_type), POINTER :: particles
1174 1119 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1175 :
1176 1119 : CALL force_env_get(force_env, cell=cell, subsys=subsys)
1177 1119 : CALL cp_subsys_get(subsys=subsys, particles=particles)
1178 1119 : particle_set => particles%els
1179 :
1180 : ! Passing gopt_f_type pointer gopt_env to particle_methods where
1181 : ! write_final_structure is defined causes a circular dependency, so it
1182 : ! is necessary to get some flags by preprocessing...
1183 1119 : keep_angles = .TRUE.
1184 1119 : keep_symmetry = .TRUE.
1185 1119 : keep_volume = .TRUE.
1186 1119 : constraint_label = "NONE"
1187 1119 : IF (gopt_env%type_id == default_cell_method_id) THEN
1188 234 : keep_angles = gopt_env%cell_env%keep_angles
1189 234 : keep_symmetry = gopt_env%cell_env%keep_symmetry
1190 234 : keep_volume = gopt_env%cell_env%keep_volume
1191 234 : constraint_label = cell_opt_constraint_label(gopt_env%cell_env%constraint_id)
1192 : END IF
1193 : CALL write_final_structure(particle_set, cell, motion_section, conv, &
1194 : keep_angles, keep_symmetry, keep_volume, &
1195 1119 : gopt_env%label, constraint_label)
1196 :
1197 1119 : IF (conv) THEN
1198 382 : it = it + 1
1199 382 : CALL write_structure_data(particle_set, cell, motion_section)
1200 382 : CALL write_restart(force_env=force_env, root_section=root_section)
1201 :
1202 382 : IF (output_unit > 0) THEN
1203 214 : WRITE (UNIT=output_unit, FMT="(/,T20,' Reevaluating energy at the minimum')")
1204 : END IF
1205 :
1206 : CALL cp_eval_at(gopt_env, x0, f=etot, master=master, &
1207 382 : para_env=para_env)
1208 382 : CALL write_geo_traj(force_env, root_section, it, etot)
1209 382 : CALL write_simulation_cell(cell, motion_section, it, 0.0_dp)
1210 : END IF
1211 :
1212 1119 : END SUBROUTINE write_final_info
1213 :
1214 : ! **************************************************************************************************
1215 : !> \brief Specific driver for dumping trajectory during a GEO_OPT
1216 : !> \param force_env ...
1217 : !> \param root_section ...
1218 : !> \param it ...
1219 : !> \param etot ...
1220 : !> \date 11.2007
1221 : !> \par History
1222 : !> 09.2010: Output of core and shell positions and forces (MK)
1223 : !> \author Teodoro Laino [tlaino] - University of Zurich
1224 : ! **************************************************************************************************
1225 17564 : SUBROUTINE write_geo_traj(force_env, root_section, it, etot)
1226 :
1227 : TYPE(force_env_type), POINTER :: force_env
1228 : TYPE(section_vals_type), POINTER :: root_section
1229 : INTEGER, INTENT(IN) :: it
1230 : REAL(KIND=dp), INTENT(IN) :: etot
1231 :
1232 : LOGICAL :: shell_adiabatic, shell_present
1233 : TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
1234 8782 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1235 : TYPE(cp_subsys_type), POINTER :: subsys
1236 : TYPE(particle_list_type), POINTER :: core_particles, shell_particles
1237 :
1238 8782 : NULLIFY (atomic_kinds)
1239 8782 : NULLIFY (atomic_kind_set)
1240 8782 : NULLIFY (core_particles)
1241 8782 : NULLIFY (shell_particles)
1242 8782 : NULLIFY (subsys)
1243 :
1244 8782 : CALL write_trajectory(force_env, root_section, it, 0.0_dp, 0.0_dp, etot)
1245 : ! Print Force
1246 8782 : CALL write_trajectory(force_env, root_section, it, 0.0_dp, 0.0_dp, etot, "FORCES", middle_name="frc")
1247 8782 : CALL force_env_get(force_env, subsys=subsys)
1248 8782 : CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds)
1249 8782 : atomic_kind_set => atomic_kinds%els
1250 : CALL get_atomic_kind_set(atomic_kind_set, &
1251 : shell_present=shell_present, &
1252 8782 : shell_adiabatic=shell_adiabatic)
1253 8782 : IF (shell_present) THEN
1254 : CALL cp_subsys_get(subsys, &
1255 : core_particles=core_particles, &
1256 3414 : shell_particles=shell_particles)
1257 : CALL write_trajectory(force_env, root_section, it=it, time=0.0_dp, dtime=0.0_dp, &
1258 : etot=etot, pk_name="SHELL_TRAJECTORY", middle_name="shpos", &
1259 3414 : particles=shell_particles)
1260 3414 : IF (shell_adiabatic) THEN
1261 : CALL write_trajectory(force_env, root_section, it=it, time=0.0_dp, dtime=0.0_dp, &
1262 : etot=etot, pk_name="SHELL_FORCES", middle_name="shfrc", &
1263 3414 : particles=shell_particles)
1264 : CALL write_trajectory(force_env, root_section, it=it, time=0.0_dp, dtime=0.0_dp, &
1265 : etot=etot, pk_name="CORE_TRAJECTORY", middle_name="copos", &
1266 3414 : particles=core_particles)
1267 : CALL write_trajectory(force_env, root_section, it=it, time=0.0_dp, dtime=0.0_dp, &
1268 : etot=etot, pk_name="CORE_FORCES", middle_name="cofrc", &
1269 3414 : particles=core_particles)
1270 : END IF
1271 : END IF
1272 :
1273 8782 : END SUBROUTINE write_geo_traj
1274 :
1275 : ! **************************************************************************************************
1276 : !> \brief ...
1277 : !> \param gopt_env ...
1278 : !> \param output_unit ...
1279 : !> \param label ...
1280 : !> \date 01.2008
1281 : !> \author Teodoro Laino [tlaino] - University of Zurich
1282 : ! **************************************************************************************************
1283 1267 : SUBROUTINE print_geo_opt_header(gopt_env, output_unit, label)
1284 :
1285 : TYPE(gopt_f_type), POINTER :: gopt_env
1286 : INTEGER, INTENT(IN) :: output_unit
1287 : CHARACTER(LEN=*), INTENT(IN) :: label
1288 :
1289 : CHARACTER(LEN=default_string_length) :: my_format, my_label
1290 : INTEGER :: ix
1291 :
1292 1267 : IF (output_unit > 0) THEN
1293 657 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") REPEAT("*", 79)
1294 657 : IF (gopt_env%dimer_rotation) THEN
1295 64 : my_label = "OPTIMIZING DIMER ROTATION"
1296 : ELSE
1297 593 : my_label = "STARTING "//gopt_env%tag(1:8)//" OPTIMIZATION"
1298 : END IF
1299 :
1300 657 : ix = (80 - 7 - LEN_TRIM(my_label))/2
1301 657 : ix = ix + 5
1302 657 : my_format = "(T2,A,T"//cp_to_string(ix)//",A,T78,A)"
1303 657 : WRITE (UNIT=output_unit, FMT=TRIM(my_format)) "***", TRIM(my_label), "***"
1304 :
1305 657 : ix = (80 - 7 - LEN_TRIM(label))/2
1306 657 : ix = ix + 5
1307 657 : my_format = "(T2,A,T"//cp_to_string(ix)//",A,T78,A)"
1308 657 : WRITE (UNIT=output_unit, FMT=TRIM(my_format)) "***", TRIM(label), "***"
1309 :
1310 657 : WRITE (UNIT=output_unit, FMT="(T2,A)") REPEAT("*", 79)
1311 657 : CALL m_flush(output_unit)
1312 : END IF
1313 1267 : END SUBROUTINE print_geo_opt_header
1314 :
1315 : ! **************************************************************************************************
1316 : !> \brief ...
1317 : !> \param gopt_env ...
1318 : !> \param output_unit ...
1319 : !> \date 01.2008
1320 : !> \author Teodoro Laino [tlaino] - University of Zurich
1321 : ! **************************************************************************************************
1322 735 : SUBROUTINE print_geo_opt_nc(gopt_env, output_unit)
1323 :
1324 : TYPE(gopt_f_type), POINTER :: gopt_env
1325 : INTEGER, INTENT(IN) :: output_unit
1326 :
1327 735 : IF (output_unit > 0) THEN
1328 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
1329 368 : "*** MAXIMUM NUMBER OF OPTIMIZATION STEPS REACHED ***"
1330 368 : IF (.NOT. gopt_env%dimer_rotation) THEN
1331 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
1332 363 : "*** EXITING GEOMETRY OPTIMIZATION ***"
1333 : ELSE
1334 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
1335 5 : "*** EXITING ROTATION OPTIMIZATION ***"
1336 : END IF
1337 368 : CALL m_flush(output_unit)
1338 : END IF
1339 :
1340 735 : END SUBROUTINE print_geo_opt_nc
1341 :
1342 : ! **************************************************************************************************
1343 : !> \brief Prints information during GEO_OPT common to all optimizers
1344 : !> \param force_env ...
1345 : !> \param root_section ...
1346 : !> \param motion_section ...
1347 : !> \param its ...
1348 : !> \param opt_energy ...
1349 : !> \date 02.2008
1350 : !> \author Teodoro Laino [tlaino] - University of Zurich
1351 : !> \version 1.0
1352 : ! **************************************************************************************************
1353 8400 : SUBROUTINE geo_opt_io(force_env, root_section, motion_section, its, opt_energy)
1354 :
1355 : TYPE(force_env_type), POINTER :: force_env
1356 : TYPE(section_vals_type), POINTER :: root_section, motion_section
1357 : INTEGER, INTENT(IN) :: its
1358 : REAL(KIND=dp), INTENT(IN) :: opt_energy
1359 :
1360 : TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
1361 8400 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1362 : TYPE(cell_type), POINTER :: cell
1363 : TYPE(cp_subsys_type), POINTER :: subsys
1364 : TYPE(distribution_1d_type), POINTER :: local_particles
1365 : TYPE(mp_para_env_type), POINTER :: para_env
1366 : TYPE(particle_list_type), POINTER :: particles
1367 8400 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1368 : TYPE(virial_type), POINTER :: virial
1369 :
1370 8400 : NULLIFY (para_env, atomic_kind_set, subsys, particle_set, &
1371 8400 : local_particles, atomic_kinds, particles)
1372 :
1373 : ! Write Restart File
1374 8400 : CALL write_restart(force_env=force_env, root_section=root_section)
1375 :
1376 : ! Write Trajectory
1377 8400 : CALL write_geo_traj(force_env, root_section, its, opt_energy)
1378 :
1379 : ! Write the stress Tensor
1380 : CALL force_env_get(force_env, cell=cell, para_env=para_env, &
1381 8400 : subsys=subsys)
1382 : CALL cp_subsys_get(subsys=subsys, atomic_kinds=atomic_kinds, local_particles=local_particles, &
1383 8400 : particles=particles, virial=virial)
1384 8400 : atomic_kind_set => atomic_kinds%els
1385 8400 : particle_set => particles%els
1386 : CALL virial_evaluate(atomic_kind_set, particle_set, local_particles, &
1387 8400 : virial, para_env)
1388 8400 : CALL write_stress_tensor_to_file(virial, cell, motion_section, its, 0.0_dp)
1389 :
1390 : ! Write the cell
1391 8400 : CALL write_simulation_cell(cell, motion_section, its, 0.0_dp)
1392 :
1393 8400 : END SUBROUTINE geo_opt_io
1394 :
1395 : ! **************************************************************************************************
1396 : !> \brief Apply coordinate transformations after cell (shape) change
1397 : !> \param gopt_env ...
1398 : !> \param cell ...
1399 : !> \param x ...
1400 : !> \param update_forces ...
1401 : !> \date 05.11.2012 (revised version of unbiase_coordinates moved here, MK)
1402 : !> \author Matthias Krack
1403 : !> \version 1.0
1404 : ! **************************************************************************************************
1405 13653 : SUBROUTINE apply_cell_change(gopt_env, cell, x, update_forces)
1406 :
1407 : TYPE(gopt_f_type), POINTER :: gopt_env
1408 : TYPE(cell_type), POINTER :: cell
1409 : REAL(KIND=dp), DIMENSION(:), POINTER :: x
1410 : LOGICAL, INTENT(IN) :: update_forces
1411 :
1412 : INTEGER :: i, iatom, idg, j, natom, nparticle, &
1413 : shell_index
1414 : REAL(KIND=dp) :: fc, fs, mass
1415 : REAL(KIND=dp), DIMENSION(3) :: s
1416 : TYPE(cell_type), POINTER :: cell_ref
1417 : TYPE(cp_subsys_type), POINTER :: subsys
1418 : TYPE(particle_list_type), POINTER :: core_particles, particles, &
1419 : shell_particles
1420 :
1421 13653 : NULLIFY (cell_ref)
1422 13653 : NULLIFY (core_particles)
1423 13653 : NULLIFY (particles)
1424 13653 : NULLIFY (shell_particles)
1425 13653 : NULLIFY (subsys)
1426 :
1427 13653 : natom = force_env_get_natom(gopt_env%force_env)
1428 13653 : nparticle = force_env_get_nparticle(gopt_env%force_env)
1429 : CALL force_env_get(gopt_env%force_env, &
1430 13653 : subsys=subsys)
1431 : CALL cp_subsys_get(subsys=subsys, &
1432 : core_particles=core_particles, &
1433 : particles=particles, &
1434 13653 : shell_particles=shell_particles)
1435 :
1436 : ! Retrieve the reference cell
1437 13653 : CALL cell_create(cell_ref)
1438 13653 : CALL cell_copy(cell, cell_ref, tag="CELL_OPT_REF")
1439 :
1440 : ! Load the updated cell information
1441 13653 : idg = 3*nparticle
1442 13653 : CALL init_cell(cell_ref, hmat=gopt_env%h_ref)
1443 13653 : CPASSERT((SIZE(x) == idg + 6))
1444 :
1445 13653 : IF (update_forces) THEN
1446 :
1447 : ! Transform particle forces back to reference cell
1448 : idg = 1
1449 339687 : DO iatom = 1, natom
1450 334354 : CALL real_to_scaled(s, x(idg:idg + 2), cell)
1451 334354 : CALL scaled_to_real(x(idg:idg + 2), s, cell_ref)
1452 339687 : idg = idg + 3
1453 : END DO
1454 :
1455 : ELSE
1456 :
1457 : ! Update the six independent components in the canonical optimization frame.
1458 108160 : gopt_env%cell_env%opt_cell%hmat = 0.0_dp
1459 33280 : DO i = 1, 3
1460 83200 : DO j = 1, i
1461 49920 : idg = idg + 1
1462 74880 : gopt_env%cell_env%opt_cell%hmat(j, i) = x(idg)
1463 : END DO
1464 : END DO
1465 8320 : CALL init_cell(gopt_env%cell_env%opt_cell)
1466 8320 : IF (gopt_env%spgr%keep_space_group) THEN
1467 798 : CALL spgr_project_cell_metric(gopt_env%spgr, gopt_env%cell_env%opt_cell)
1468 : ! Keep the optimizer variables synchronized with the projected cell.
1469 798 : idg = 3*nparticle
1470 3192 : DO i = 1, 3
1471 7980 : DO j = 1, i
1472 4788 : idg = idg + 1
1473 7182 : x(idg) = gopt_env%cell_env%opt_cell%hmat(j, i)
1474 : END DO
1475 : END DO
1476 : END IF
1477 :
1478 : ! Reconstruct the physical cell in the fixed orientation of the input cell.
1479 : cell%hmat = MATMUL(gopt_env%cell_env%opt_to_input, &
1480 648960 : gopt_env%cell_env%opt_cell%hmat)
1481 8320 : CALL init_cell(cell)
1482 8320 : CALL cp_subsys_set(subsys, cell=cell)
1483 :
1484 : ! Retrieve particle coordinates for the current cell
1485 8320 : idg = 1
1486 635672 : DO iatom = 1, natom
1487 627352 : CALL real_to_scaled(s, x(idg:idg + 2), cell_ref)
1488 627352 : shell_index = particles%els(iatom)%shell_index
1489 627352 : IF (shell_index == 0) THEN
1490 307612 : CALL scaled_to_real(particles%els(iatom)%r, s, cell)
1491 : ELSE
1492 319740 : CALL scaled_to_real(core_particles%els(shell_index)%r, s, cell)
1493 319740 : i = 3*(natom + shell_index - 1) + 1
1494 319740 : CALL real_to_scaled(s, x(i:i + 2), cell_ref)
1495 319740 : CALL scaled_to_real(shell_particles%els(shell_index)%r, s, cell)
1496 : ! Update atomic position due to core and shell motion
1497 319740 : mass = particles%els(iatom)%atomic_kind%mass
1498 319740 : fc = core_particles%els(shell_index)%atomic_kind%shell%mass_core/mass
1499 319740 : fs = shell_particles%els(shell_index)%atomic_kind%shell%mass_shell/mass
1500 : particles%els(iatom)%r(1:3) = fc*core_particles%els(shell_index)%r(1:3) + &
1501 2557920 : fs*shell_particles%els(shell_index)%r(1:3)
1502 : END IF
1503 635672 : idg = idg + 3
1504 : END DO
1505 : END IF
1506 :
1507 13653 : CALL cell_release(cell_ref)
1508 :
1509 13653 : END SUBROUTINE apply_cell_change
1510 :
1511 : ! **************************************************************************************************
1512 : !> \brief Identify the space group, symmetrize x0 and print it (KEEP_SPACE_GROUP)
1513 : !> \param gopt_env ...
1514 : !> \param x0 ...
1515 : !> \param geo_section ...
1516 : !> \param output_unit ...
1517 : ! **************************************************************************************************
1518 1267 : SUBROUTINE gopt_f_space_group_init(gopt_env, x0, geo_section, output_unit)
1519 : TYPE(gopt_f_type), POINTER :: gopt_env
1520 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: x0
1521 : TYPE(section_vals_type), POINTER :: geo_section
1522 : INTEGER, INTENT(IN) :: output_unit
1523 :
1524 : TYPE(cp_subsys_type), POINTER :: subsys
1525 :
1526 1267 : NULLIFY (subsys)
1527 :
1528 1267 : IF (.NOT. gopt_env%spgr%keep_space_group) RETURN
1529 36 : SELECT CASE (gopt_env%type_id)
1530 : CASE (default_minimization_method_id, default_ts_method_id, default_cell_method_id)
1531 18 : CALL force_env_get(gopt_env%force_env, subsys=subsys)
1532 18 : CALL identify_space_group(subsys, geo_section, gopt_env, output_unit)
1533 18 : CALL spgr_apply_rotations_coord(gopt_env%spgr, x0)
1534 18 : CALL print_spgr(gopt_env%spgr)
1535 : CASE DEFAULT
1536 18 : gopt_env%spgr%keep_space_group = .FALSE.
1537 : END SELECT
1538 : END SUBROUTINE gopt_f_space_group_init
1539 :
1540 : ! **************************************************************************************************
1541 : !> \brief Identify and print the space group of the current structure (SHOW_SPACE_GROUP)
1542 : !> \param gopt_env ...
1543 : !> \param geo_section ...
1544 : !> \param output_unit ...
1545 : ! **************************************************************************************************
1546 1267 : SUBROUTINE gopt_f_space_group_show(gopt_env, geo_section, output_unit)
1547 : TYPE(gopt_f_type), POINTER :: gopt_env
1548 : TYPE(section_vals_type), POINTER :: geo_section
1549 : INTEGER, INTENT(IN) :: output_unit
1550 :
1551 : TYPE(cp_subsys_type), POINTER :: subsys
1552 :
1553 1267 : NULLIFY (subsys)
1554 :
1555 1267 : IF (.NOT. gopt_env%spgr%show_space_group) RETURN
1556 6 : CALL force_env_get(gopt_env%force_env, subsys=subsys)
1557 6 : CALL identify_space_group(subsys, geo_section, gopt_env, output_unit)
1558 6 : CALL print_spgr(gopt_env%spgr)
1559 : END SUBROUTINE gopt_f_space_group_show
1560 :
1561 : ! **************************************************************************************************
1562 : !> \brief Project optimizer coordinates onto the space-group-invariant subspace
1563 : !> (KEEP_SPACE_GROUP)
1564 : !> \param gopt_env ...
1565 : !> \param x ...
1566 : ! **************************************************************************************************
1567 28497 : SUBROUTINE gopt_f_project_coord(gopt_env, x)
1568 : TYPE(gopt_f_type), POINTER :: gopt_env
1569 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: x
1570 :
1571 28497 : IF (gopt_env%spgr%keep_space_group) CALL spgr_apply_rotations_coord(gopt_env%spgr, x)
1572 28497 : END SUBROUTINE gopt_f_project_coord
1573 :
1574 : ! **************************************************************************************************
1575 : !> \brief Project a gradient, displacement or search direction (linear, no
1576 : !> translations)
1577 : !> \param gopt_env ...
1578 : !> \param v ...
1579 : ! **************************************************************************************************
1580 41433 : SUBROUTINE gopt_f_project_vector(gopt_env, v)
1581 : TYPE(gopt_f_type), POINTER :: gopt_env
1582 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: v
1583 :
1584 41433 : IF (gopt_env%spgr%keep_space_group) CALL spgr_apply_rotations_force(gopt_env%spgr, v)
1585 41433 : END SUBROUTINE gopt_f_project_vector
1586 :
1587 : ! **************************************************************************************************
1588 : !> \brief Restrict a proposed optimizer step: cell constraints (CELL_OPT),
1589 : !> then space group
1590 : !> \param gopt_env ...
1591 : !> \param dr ...
1592 : ! **************************************************************************************************
1593 4382 : SUBROUTINE gopt_f_project_step(gopt_env, dr)
1594 : TYPE(gopt_f_type), POINTER :: gopt_env
1595 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT), TARGET :: dr
1596 :
1597 4382 : REAL(KIND=dp), DIMENSION(:), POINTER :: cell_step
1598 :
1599 4382 : IF (gopt_env%type_id == default_cell_method_id) THEN
1600 604 : cell_step => dr(SIZE(dr) - 5:SIZE(dr))
1601 : CALL apply_cell_constraints(cell_step, gopt_env%cell_env%opt_cell, &
1602 : gopt_env%cell_env%keep_angles, gopt_env%cell_env%keep_symmetry, &
1603 604 : gopt_env%cell_env%constraint_id)
1604 : END IF
1605 4382 : CALL gopt_f_project_vector(gopt_env, dr)
1606 4382 : END SUBROUTINE gopt_f_project_step
1607 :
1608 : END MODULE gopt_f_methods
|
