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