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