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