Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief evaluete the potential energy and its gradients using an array
10 : !> with same dimension as the particle_set
11 : !> \param gopt_env the geometry optimization environment
12 : !> \param x the position where the function should be evaluated
13 : !> \param f the function value
14 : !> \param gradient the value of its gradient
15 : !> \param master ...
16 : !> \param final_evaluation ...
17 : !> \param para_env ...
18 : !> \par History
19 : !> CELL OPTIMIZATION: Teodoro Laino [tlaino] - University of Zurich - 03.2008
20 : !> 07.2020 Pierre Cazade [pcazade] Space Group Symmetry
21 : !> \author Teodoro Laino [tlaino] - University of Zurich - 01.2008
22 : ! **************************************************************************************************
23 22715 : SUBROUTINE cp_eval_at(gopt_env, x, f, gradient, master, &
24 : final_evaluation, para_env)
25 :
26 : USE cp_log_handling, ONLY: cp_logger_type
27 : USE averages_types, ONLY: average_quantities_type, &
28 : create_averages, &
29 : release_averages
30 : USE bibliography, ONLY: Henkelman1999, &
31 : cite_reference
32 : USE cell_opt_utils, ONLY: get_dg_dh, &
33 : gopt_new_logger_create, &
34 : gopt_new_logger_release, &
35 : rescale_new_cell_volume
36 : USE cell_types, ONLY: cell_type
37 : USE cell_methods, ONLY: write_cell
38 : USE message_passing, ONLY: mp_para_env_type
39 : USE cp_subsys_types, ONLY: cp_subsys_get, &
40 : cp_subsys_type, &
41 : pack_subsys_particles, &
42 : unpack_subsys_particles
43 : USE dimer_methods, ONLY: cp_eval_at_ts
44 : USE force_env_methods, ONLY: force_env_calc_energy_force
45 : USE force_env_types, ONLY: force_env_get, &
46 : force_env_get_nparticle
47 : USE geo_opt, ONLY: cp_geo_opt
48 : USE gopt_f_types, ONLY: gopt_f_type
49 : USE gopt_f_methods, ONLY: apply_cell_change
50 : USE input_constants, ONLY: default_minimization_method_id, &
51 : default_ts_method_id, &
52 : default_cell_direct_id, &
53 : default_cell_method_id, &
54 : default_cell_geo_opt_id, &
55 : default_cell_md_id, &
56 : default_shellcore_method_id, &
57 : nvt_ensemble, &
58 : mol_dyn_run, &
59 : geo_opt_run, &
60 : cell_opt_run, &
61 : fix_none
62 : USE input_section_types, ONLY: section_vals_get, &
63 : section_vals_get_subs_vals, &
64 : section_vals_type, &
65 : section_vals_val_get
66 : USE md_run, ONLY: qs_mol_dyn
67 : USE kinds, ONLY: dp, &
68 : default_string_length
69 : USE particle_list_types, ONLY: particle_list_type
70 : USE particle_methods, ONLY: write_structure_data
71 : USE virial_methods, ONLY: virial_update
72 : USE virial_types, ONLY: virial_type
73 : USE cp_log_handling, ONLY: cp_add_default_logger, &
74 : cp_rm_default_logger
75 : USE space_groups_types, ONLY: spgr_type
76 : USE space_groups, ONLY: spgr_apply_rotations_stress, &
77 : spgr_apply_rotations_coord, &
78 : spgr_apply_rotations_force, &
79 : spgr_write_stress_tensor
80 :
81 : #include "../base/base_uses.f90"
82 : IMPLICIT NONE
83 : TYPE(gopt_f_type), POINTER :: gopt_env
84 : REAL(KIND=dp), DIMENSION(:), POINTER :: x
85 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: f
86 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, &
87 : POINTER :: gradient
88 : INTEGER, INTENT(IN) :: master
89 : LOGICAL, INTENT(IN), OPTIONAL :: final_evaluation
90 : TYPE(mp_para_env_type), POINTER :: para_env
91 :
92 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_eval_at'
93 :
94 : INTEGER :: ensemble, handle, idg, idir, ip, &
95 : nparticle, nsize, shell_index
96 : LOGICAL :: explicit, my_final_evaluation
97 : REAL(KIND=dp) :: f_ts, potential_energy
98 : REAL(KIND=dp), DIMENSION(3, 3) :: av_ptens
99 22715 : REAL(KIND=dp), DIMENSION(:), POINTER :: cell_gradient, gradient_ts
100 : TYPE(cell_type), POINTER :: cell
101 : TYPE(cp_subsys_type), POINTER :: subsys
102 : TYPE(particle_list_type), POINTER :: core_particles, particles, &
103 : shell_particles
104 : TYPE(virial_type), POINTER :: virial
105 : TYPE(cp_logger_type), POINTER :: new_logger
106 : CHARACTER(LEN=default_string_length) :: project_name
107 : TYPE(average_quantities_type), POINTER :: averages
108 : TYPE(section_vals_type), POINTER :: work, avgs_section
109 : TYPE(spgr_type), POINTER :: spgr
110 :
111 22715 : NULLIFY (averages)
112 22715 : NULLIFY (cell)
113 22715 : NULLIFY (core_particles)
114 22715 : NULLIFY (gradient_ts)
115 22715 : NULLIFY (particles)
116 22715 : NULLIFY (shell_particles)
117 22715 : NULLIFY (subsys)
118 22715 : NULLIFY (virial)
119 22715 : NULLIFY (new_logger)
120 22715 : NULLIFY (spgr)
121 :
122 22715 : CALL timeset(routineN, handle)
123 :
124 22715 : CALL force_env_get(gopt_env%force_env, subsys=subsys, cell=cell)
125 : CALL cp_subsys_get(subsys, &
126 : core_particles=core_particles, &
127 : particles=particles, &
128 : shell_particles=shell_particles, &
129 22715 : virial=virial)
130 :
131 22715 : spgr => gopt_env%spgr
132 :
133 22715 : my_final_evaluation = .FALSE.
134 22715 : IF (PRESENT(final_evaluation)) my_final_evaluation = final_evaluation
135 :
136 36684 : SELECT CASE (gopt_env%type_id)
137 : CASE (default_minimization_method_id, default_ts_method_id)
138 13969 : CALL unpack_subsys_particles(subsys=subsys, r=x)
139 13969 : CALL write_structure_data(particles%els, cell, gopt_env%motion_section)
140 26122 : SELECT CASE (gopt_env%type_id)
141 : CASE (default_minimization_method_id)
142 : ! Geometry Minimization
143 : CALL force_env_calc_energy_force(gopt_env%force_env, &
144 : calc_force=PRESENT(gradient), &
145 12153 : require_consistent_energy_force=gopt_env%require_consistent_energy_force)
146 : ! Possibly take the potential energy
147 12153 : IF (PRESENT(f)) THEN
148 12153 : CALL force_env_get(gopt_env%force_env, potential_energy=f)
149 : END IF
150 : ! Possibly take the gradients
151 12153 : IF (PRESENT(gradient)) THEN
152 11607 : IF (master == para_env%mepos) THEN ! we are on the master
153 10534 : CALL pack_subsys_particles(subsys=subsys, f=gradient, fscale=-1.0_dp)
154 10534 : IF (spgr%keep_space_group) THEN
155 0 : CALL spgr_apply_rotations_force(spgr, gradient)
156 0 : CALL unpack_subsys_particles(subsys=subsys, f=gradient, fscale=-1.0_dp)
157 : END IF
158 : END IF
159 : END IF
160 : CASE (default_ts_method_id)
161 : ! Transition State Optimization
162 5448 : ALLOCATE (gradient_ts(particles%n_els*3))
163 : ! Real calculation of energy and forces for transition state optimization:
164 : ! When doing dimer methods forces have to be always computed since the function
165 : ! to minimize is not the energy but the effective force
166 1816 : CALL cp_eval_at_ts(gopt_env, x, f_ts, gradient_ts, calc_force=.TRUE.)
167 1816 : CALL cite_reference(Henkelman1999)
168 : ! Possibly take the potential energy
169 1816 : IF (PRESENT(f)) f = f_ts
170 : ! Possibly take the gradients
171 1816 : IF (PRESENT(gradient)) THEN
172 854 : IF (master == para_env%mepos) THEN ! we are on the master
173 854 : CPASSERT(ASSOCIATED(gradient))
174 29452 : gradient = gradient_ts
175 : END IF
176 : END IF
177 15785 : DEALLOCATE (gradient_ts)
178 : END SELECT
179 : ! This call is necessary for QM/MM if a Translation is applied
180 : ! this makes the geometry optimizer consistent
181 13969 : CALL unpack_subsys_particles(subsys=subsys, r=x)
182 : CASE (default_cell_method_id)
183 : ! Check for VIRIAL
184 8406 : IF (.NOT. virial%pv_availability) &
185 : CALL cp_abort(__LOCATION__, &
186 : "For a cell optimization task with CELL_OPT/TYPE "// &
187 : "DIRECT_CELL_OPT, the FORCE_EVAL/STRESS_TENSOR "// &
188 : "keyword MUST be defined in the input file for the "// &
189 0 : "evaluation of the stress tensor, but none is found!")
190 8406 : IF (gopt_env%cell_env%keep_volume) THEN
191 2278 : nparticle = force_env_get_nparticle(gopt_env%force_env)
192 4556 : SELECT CASE (gopt_env%cell_method_id)
193 : CASE (default_cell_direct_id)
194 2278 : idg = 3*nparticle
195 : CASE (default_cell_geo_opt_id, default_cell_md_id)
196 2278 : idg = 0
197 : END SELECT
198 2278 : CALL rescale_new_cell_volume(cell%deth, x, idg)
199 : END IF
200 16198 : SELECT CASE (gopt_env%cell_method_id)
201 : CASE (default_cell_direct_id)
202 7452 : CALL apply_cell_change(gopt_env, cell, x, update_forces=.FALSE.)
203 : ! Possibly output the new cell used for the next calculation
204 7452 : CALL write_cell(cell, gopt_env%geo_section)
205 : ! Compute the pressure tensor
206 7452 : BLOCK
207 : TYPE(virial_type) :: virial_avg
208 : CALL force_env_calc_energy_force(gopt_env%force_env, &
209 : calc_force=PRESENT(gradient), &
210 7452 : require_consistent_energy_force=gopt_env%require_consistent_energy_force)
211 : ! Possibly take the potential energy
212 7452 : CALL force_env_get(gopt_env%force_env, potential_energy=potential_energy)
213 7452 : virial_avg = virial
214 7452 : CALL virial_update(virial_avg, subsys, para_env)
215 7452 : IF (PRESENT(f)) THEN
216 7452 : CALL force_env_get(gopt_env%force_env, potential_energy=f)
217 : END IF
218 : ! Possibly take the gradients
219 1848096 : IF (PRESENT(gradient)) THEN
220 6262 : CPASSERT(ANY(virial_avg%pv_total /= 0))
221 : ! Convert the average ptens
222 81406 : av_ptens(:, :) = virial_avg%pv_total(:, :)/cell%deth
223 6262 : IF (master == para_env%mepos) THEN ! we are on the master
224 4744 : CPASSERT(ASSOCIATED(gradient))
225 4744 : nparticle = force_env_get_nparticle(gopt_env%force_env)
226 4744 : nsize = 3*nparticle
227 4744 : CPASSERT((SIZE(gradient) == nsize + 6))
228 4744 : CALL pack_subsys_particles(subsys=subsys, f=gradient(1:nsize), fscale=-1.0_dp)
229 4744 : CALL apply_cell_change(gopt_env, cell, gradient, update_forces=.TRUE.)
230 4744 : IF (spgr%keep_space_group) THEN
231 544 : CALL spgr_apply_rotations_force(spgr, gradient)
232 544 : CALL spgr_apply_rotations_stress(spgr, cell, av_ptens)
233 544 : CALL spgr_write_stress_tensor(av_ptens, spgr)
234 : END IF
235 4744 : cell_gradient => gradient(nsize + 1:nsize + 6)
236 33208 : cell_gradient = 0.0_dp
237 : CALL get_dg_dh(cell_gradient, av_ptens, gopt_env%cell_env%pres_ext, cell, gopt_env%cell_env%mtrx, &
238 : keep_angles=gopt_env%cell_env%keep_angles, &
239 : keep_symmetry=gopt_env%cell_env%keep_symmetry, &
240 : pres_int=gopt_env%cell_env%pres_int, &
241 : pres_constr=gopt_env%cell_env%pres_constr, &
242 4744 : constraint_id=gopt_env%cell_env%constraint_id)
243 : END IF
244 : ! some callers expect pres_int to be available on all ranks. Also, here master is not necessarily a single rank.
245 : ! Assume at least master==0
246 6262 : CALL para_env%bcast(gopt_env%cell_env%pres_int, 0)
247 6262 : IF (gopt_env%cell_env%constraint_id /= fix_none) &
248 24 : CALL para_env%bcast(gopt_env%cell_env%pres_constr, 0)
249 : END IF
250 : END BLOCK
251 : CASE (default_cell_geo_opt_id, default_cell_md_id)
252 954 : CALL apply_cell_change(gopt_env, cell, x, update_forces=.FALSE.)
253 : ! Possibly output the new cell used for the next calculation
254 954 : CALL write_cell(cell, gopt_env%geo_section)
255 : ! Compute the pressure tensor
256 954 : BLOCK
257 : TYPE(virial_type) :: virial_avg
258 954 : IF (my_final_evaluation) THEN
259 : CALL force_env_calc_energy_force(gopt_env%force_env, &
260 : calc_force=PRESENT(gradient), &
261 24 : require_consistent_energy_force=gopt_env%require_consistent_energy_force)
262 24 : IF (PRESENT(f)) THEN
263 24 : CALL force_env_get(gopt_env%force_env, potential_energy=f)
264 : END IF
265 : ELSE
266 1840 : SELECT CASE (gopt_env%cell_method_id)
267 : CASE (default_cell_geo_opt_id)
268 910 : work => section_vals_get_subs_vals(gopt_env%motion_section, "GEO_OPT")
269 910 : CALL section_vals_get(work, explicit=explicit)
270 910 : IF (.NOT. explicit) &
271 : CALL cp_abort(__LOCATION__, &
272 : "For a cell optimization task with CELL_OPT/TYPE "// &
273 : "GEO_OPT, besides the MOTION/CELL_OPT section, the "// &
274 : "MOTION/GEO_OPT section MUST also be provided in "// &
275 : "the input file for the evaluation of the stress "// &
276 0 : "tensor, but none is found!")
277 : ! Perform a geometry optimization
278 : CALL gopt_new_logger_create(new_logger, gopt_env%force_env%root_section, para_env, &
279 910 : project_name, id_run=geo_opt_run)
280 910 : CALL cp_add_default_logger(new_logger)
281 910 : CALL cp_geo_opt(gopt_env%force_env, gopt_env%globenv, eval_opt_geo=.FALSE.)
282 910 : CALL force_env_get(gopt_env%force_env, potential_energy=potential_energy)
283 910 : virial_avg = virial
284 : CASE (default_cell_md_id)
285 20 : work => section_vals_get_subs_vals(gopt_env%motion_section, "MD")
286 20 : avgs_section => section_vals_get_subs_vals(work, "AVERAGES")
287 20 : CALL section_vals_get(work, explicit=explicit)
288 20 : IF (.NOT. explicit) &
289 : CALL cp_abort(__LOCATION__, &
290 : "For a cell optimization task with CELL_OPT/TYPE MD, "// &
291 : "besides the MOTION/CELL_OPT section, the MOTION/MD "// &
292 : "section MUST also be provided in the input file for "// &
293 0 : "the evaluation of the stress tensor, but none is found!")
294 : ! Only NVT ensemble is allowed..
295 20 : CALL section_vals_val_get(gopt_env%motion_section, "MD%ENSEMBLE", i_val=ensemble)
296 20 : IF (ensemble /= nvt_ensemble) &
297 : CALL cp_abort(__LOCATION__, &
298 : "For a cell optimization task with CELL_OPT/TYPE MD, "// &
299 : "the MOTION/MD/ENSEMBLE keyword MUST be set to NVT "// &
300 0 : "and any other choice of ensemble is not supported!")
301 : ! Perform a molecular dynamics
302 : CALL gopt_new_logger_create(new_logger, gopt_env%force_env%root_section, para_env, &
303 20 : project_name, id_run=mol_dyn_run)
304 20 : CALL cp_add_default_logger(new_logger)
305 20 : CALL create_averages(averages, avgs_section, virial_avg=.TRUE., force_env=gopt_env%force_env)
306 20 : CALL qs_mol_dyn(gopt_env%force_env, gopt_env%globenv, averages, rm_restart_info=.FALSE.)
307 : ! Retrieve the average of the stress tensor and the average of the potential energy
308 20 : potential_energy = averages%avepot
309 20 : virial_avg = averages%virial
310 20 : CALL release_averages(averages)
311 : CASE DEFAULT
312 : ! Should never reach this point
313 2790 : CPABORT("Invalid or not yet implemented type of cell optimization")
314 : END SELECT
315 930 : CALL cp_rm_default_logger()
316 : CALL gopt_new_logger_release(new_logger, gopt_env%force_env%root_section, para_env, project_name, &
317 930 : cell_opt_run)
318 : ! Update the virial
319 930 : CALL virial_update(virial_avg, subsys, para_env)
320 : ! Possibly take give back the potential energy
321 930 : IF (PRESENT(f)) THEN
322 930 : f = potential_energy
323 : END IF
324 : END IF
325 : ! Possibly give back the gradients
326 236592 : IF (PRESENT(gradient)) THEN
327 930 : CPASSERT(ANY(virial_avg%pv_total /= 0))
328 : ! Convert the average ptens
329 12090 : av_ptens(:, :) = virial_avg%pv_total(:, :)/cell%deth
330 930 : IF (master == para_env%mepos) THEN ! we are on the master
331 904 : CPASSERT(ASSOCIATED(gradient))
332 904 : IF (spgr%keep_space_group) THEN
333 0 : CALL spgr_apply_rotations_stress(spgr, cell, av_ptens)
334 0 : CALL spgr_write_stress_tensor(av_ptens, spgr)
335 : END IF
336 : ! Compute the gradients on the cell
337 : CALL get_dg_dh(gradient, av_ptens, gopt_env%cell_env%pres_ext, cell, gopt_env%cell_env%mtrx, &
338 : keep_angles=gopt_env%cell_env%keep_angles, &
339 : keep_symmetry=gopt_env%cell_env%keep_symmetry, &
340 : pres_int=gopt_env%cell_env%pres_int, &
341 : pres_constr=gopt_env%cell_env%pres_constr, &
342 904 : constraint_id=gopt_env%cell_env%constraint_id)
343 : END IF
344 : ! some callers expect pres_int to be available on all ranks. Also, here master is not necessarily a single rank.
345 : ! Assume at least master==0
346 930 : CALL para_env%bcast(gopt_env%cell_env%pres_int, 0)
347 930 : IF (gopt_env%cell_env%constraint_id /= fix_none) &
348 0 : CALL para_env%bcast(gopt_env%cell_env%pres_constr, 0)
349 : END IF
350 : END BLOCK
351 : CASE DEFAULT
352 8406 : CPABORT("Invalid or not yet implemented type of cell optimization")
353 : END SELECT
354 : CASE (default_shellcore_method_id)
355 340 : idg = 0
356 32980 : DO ip = 1, particles%n_els
357 32640 : shell_index = particles%els(ip)%shell_index
358 32980 : IF (shell_index /= 0) THEN
359 130560 : DO idir = 1, 3
360 97920 : idg = 3*(shell_index - 1) + idir
361 130560 : shell_particles%els(shell_index)%r(idir) = core_particles%els(ip)%r(idir) - x(idg)
362 : END DO
363 : END IF
364 : END DO
365 340 : CALL write_structure_data(particles%els, cell, gopt_env%motion_section)
366 :
367 : ! Shell-core optimization
368 : CALL force_env_calc_energy_force(gopt_env%force_env, &
369 : calc_force=PRESENT(gradient), &
370 340 : require_consistent_energy_force=gopt_env%require_consistent_energy_force)
371 :
372 : ! Possibly take the potential energy
373 340 : IF (PRESENT(f)) THEN
374 340 : CALL force_env_get(gopt_env%force_env, potential_energy=f)
375 : END IF
376 :
377 : ! Possibly take the gradients
378 340 : IF (PRESENT(gradient)) THEN
379 340 : IF (master == para_env%mepos) THEN ! we are on the master
380 340 : CPASSERT(ASSOCIATED(gradient))
381 340 : idg = 0
382 32980 : DO ip = 1, shell_particles%n_els
383 130900 : DO idir = 1, 3
384 97920 : idg = idg + 1
385 130560 : gradient(idg) = -(core_particles%els(ip)%f(idir) - shell_particles%els(ip)%f(idir))
386 : END DO
387 : END DO
388 : END IF
389 : END IF
390 : CASE DEFAULT
391 22715 : CPABORT("Invalid or not yet implemented type of optimization")
392 : END SELECT
393 :
394 22715 : CALL timestop(handle)
395 :
396 45430 : END SUBROUTINE cp_eval_at
|
