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