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