Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 20342 : 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 : USE cell_types, ONLY: cell_type
36 : USE cell_methods, ONLY: write_cell
37 : USE message_passing, ONLY: mp_para_env_type
38 : USE cp_subsys_types, ONLY: cp_subsys_get, &
39 : cp_subsys_type, &
40 : pack_subsys_particles, &
41 : unpack_subsys_particles
42 : USE dimer_methods, ONLY: cp_eval_at_ts
43 : USE force_env_methods, ONLY: force_env_calc_energy_force
44 : USE force_env_types, ONLY: force_env_get, &
45 : force_env_get_nparticle
46 : USE geo_opt, ONLY: cp_geo_opt
47 : USE gopt_f_types, ONLY: gopt_f_type
48 : USE gopt_f_methods, ONLY: apply_cell_change
49 : USE input_constants, ONLY: default_minimization_method_id, &
50 : default_ts_method_id, &
51 : default_cell_direct_id, &
52 : default_cell_method_id, &
53 : default_cell_geo_opt_id, &
54 : default_cell_md_id, &
55 : default_shellcore_method_id, &
56 : nvt_ensemble, &
57 : mol_dyn_run, &
58 : geo_opt_run, &
59 : cell_opt_run, &
60 : fix_none
61 : USE input_section_types, ONLY: section_vals_get, &
62 : section_vals_get_subs_vals, &
63 : section_vals_type, &
64 : section_vals_val_get
65 : USE md_run, ONLY: qs_mol_dyn
66 : USE kinds, ONLY: dp, &
67 : default_string_length
68 : USE particle_list_types, ONLY: particle_list_type
69 : USE particle_methods, ONLY: write_structure_data
70 : USE virial_methods, ONLY: virial_update
71 : USE virial_types, ONLY: virial_type
72 : USE cp_log_handling, ONLY: cp_add_default_logger, &
73 : cp_rm_default_logger
74 : USE space_groups_types, ONLY: spgr_type
75 : USE space_groups, ONLY: spgr_apply_rotations_stress, &
76 : spgr_apply_rotations_coord, &
77 : spgr_apply_rotations_force, &
78 : spgr_write_stress_tensor
79 :
80 : #include "../base/base_uses.f90"
81 : IMPLICIT NONE
82 : TYPE(gopt_f_type), POINTER :: gopt_env
83 : REAL(KIND=dp), DIMENSION(:), POINTER :: x
84 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: f
85 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, &
86 : POINTER :: gradient
87 : INTEGER, INTENT(IN) :: master
88 : LOGICAL, INTENT(IN), OPTIONAL :: final_evaluation
89 : TYPE(mp_para_env_type), POINTER :: para_env
90 :
91 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_eval_at'
92 :
93 : INTEGER :: ensemble, handle, idg, idir, ip, &
94 : nparticle, nsize, shell_index
95 : LOGICAL :: explicit, my_final_evaluation
96 : REAL(KIND=dp) :: f_ts, potential_energy
97 : REAL(KIND=dp), DIMENSION(3, 3) :: av_ptens
98 20342 : REAL(KIND=dp), DIMENSION(:), POINTER :: cell_gradient, gradient_ts
99 : TYPE(cell_type), POINTER :: cell
100 : TYPE(cp_subsys_type), POINTER :: subsys
101 : TYPE(particle_list_type), POINTER :: core_particles, particles, &
102 : shell_particles
103 : TYPE(virial_type), POINTER :: virial
104 : TYPE(cp_logger_type), POINTER :: new_logger
105 : CHARACTER(LEN=default_string_length) :: project_name
106 : TYPE(average_quantities_type), POINTER :: averages
107 : TYPE(section_vals_type), POINTER :: work, avgs_section
108 : TYPE(spgr_type), POINTER :: spgr
109 :
110 20342 : NULLIFY (averages)
111 20342 : NULLIFY (cell)
112 20342 : NULLIFY (core_particles)
113 20342 : NULLIFY (gradient_ts)
114 20342 : NULLIFY (particles)
115 20342 : NULLIFY (shell_particles)
116 20342 : NULLIFY (subsys)
117 20342 : NULLIFY (virial)
118 20342 : NULLIFY (new_logger)
119 20342 : NULLIFY (spgr)
120 :
121 20342 : CALL timeset(routineN, handle)
122 :
123 20342 : CALL force_env_get(gopt_env%force_env, subsys=subsys, cell=cell)
124 : CALL cp_subsys_get(subsys, &
125 : core_particles=core_particles, &
126 : particles=particles, &
127 : shell_particles=shell_particles, &
128 20342 : virial=virial)
129 :
130 20342 : spgr => gopt_env%spgr
131 :
132 20342 : my_final_evaluation = .FALSE.
133 20342 : IF (PRESENT(final_evaluation)) my_final_evaluation = final_evaluation
134 :
135 34212 : SELECT CASE (gopt_env%type_id)
136 : CASE (default_minimization_method_id, default_ts_method_id)
137 13870 : CALL unpack_subsys_particles(subsys=subsys, r=x)
138 13870 : CALL write_structure_data(particles%els, cell, gopt_env%motion_section)
139 25924 : SELECT CASE (gopt_env%type_id)
140 : CASE (default_minimization_method_id)
141 : ! Geometry Minimization
142 : CALL force_env_calc_energy_force(gopt_env%force_env, &
143 : calc_force=PRESENT(gradient), &
144 12054 : require_consistent_energy_force=gopt_env%require_consistent_energy_force)
145 : ! Possibly take the potential energy
146 12054 : IF (PRESENT(f)) THEN
147 12054 : CALL force_env_get(gopt_env%force_env, potential_energy=f)
148 : END IF
149 : ! Possibly take the gradients
150 12054 : IF (PRESENT(gradient)) THEN
151 11536 : IF (master == para_env%mepos) THEN ! we are on the master
152 10478 : CALL pack_subsys_particles(subsys=subsys, f=gradient, fscale=-1.0_dp)
153 10478 : IF (spgr%keep_space_group) THEN
154 0 : CALL spgr_apply_rotations_force(spgr, gradient)
155 0 : CALL unpack_subsys_particles(subsys=subsys, f=gradient, fscale=-1.0_dp)
156 : END IF
157 : END IF
158 : END IF
159 : CASE (default_ts_method_id)
160 : ! Transition State Optimization
161 5448 : ALLOCATE (gradient_ts(particles%n_els*3))
162 : ! Real calculation of energy and forces for transition state optimization:
163 : ! When doing dimer methods forces have to be always computed since the function
164 : ! to minimize is not the energy but the effective force
165 1816 : CALL cp_eval_at_ts(gopt_env, x, f_ts, gradient_ts, calc_force=.TRUE.)
166 1816 : CALL cite_reference(Henkelman1999)
167 : ! Possibly take the potential energy
168 1816 : IF (PRESENT(f)) f = f_ts
169 : ! Possibly take the gradients
170 1816 : IF (PRESENT(gradient)) THEN
171 854 : IF (master == para_env%mepos) THEN ! we are on the master
172 854 : CPASSERT(ASSOCIATED(gradient))
173 29452 : gradient = gradient_ts
174 : END IF
175 : END IF
176 15686 : DEALLOCATE (gradient_ts)
177 : END SELECT
178 : ! This call is necessary for QM/MM if a Translation is applied
179 : ! this makes the geometry optimizer consistent
180 13870 : CALL unpack_subsys_particles(subsys=subsys, r=x)
181 : CASE (default_cell_method_id)
182 : ! Check for VIRIAL
183 6132 : IF (.NOT. virial%pv_availability) &
184 : CALL cp_abort(__LOCATION__, &
185 : "Cell optimization requested but FORCE_EVAL%STRESS_TENSOR was not defined! "// &
186 0 : "Activate the evaluation of the stress tensor for cell optimization!")
187 11650 : SELECT CASE (gopt_env%cell_method_id)
188 : CASE (default_cell_direct_id)
189 5178 : CALL apply_cell_change(gopt_env, cell, x, update_forces=.FALSE.)
190 : ! Possibly output the new cell used for the next calculation
191 5178 : CALL write_cell(cell, gopt_env%geo_section)
192 : ! Compute the pressure tensor
193 5178 : BLOCK
194 : TYPE(virial_type) :: virial_avg
195 : CALL force_env_calc_energy_force(gopt_env%force_env, &
196 : calc_force=PRESENT(gradient), &
197 5178 : require_consistent_energy_force=gopt_env%require_consistent_energy_force)
198 : ! Possibly take the potential energy
199 5178 : CALL force_env_get(gopt_env%force_env, potential_energy=potential_energy)
200 5178 : virial_avg = virial
201 5178 : CALL virial_update(virial_avg, subsys, para_env)
202 5178 : IF (PRESENT(f)) THEN
203 5178 : CALL force_env_get(gopt_env%force_env, potential_energy=f)
204 : END IF
205 : ! Possibly take the gradients
206 1284144 : IF (PRESENT(gradient)) THEN
207 4686 : CPASSERT(ANY(virial_avg%pv_total /= 0))
208 : ! Convert the average ptens
209 60918 : av_ptens(:, :) = virial_avg%pv_total(:, :)/cell%deth
210 4686 : IF (master == para_env%mepos) THEN ! we are on the master
211 3166 : CPASSERT(ASSOCIATED(gradient))
212 3166 : nparticle = force_env_get_nparticle(gopt_env%force_env)
213 3166 : nsize = 3*nparticle
214 3166 : CPASSERT((SIZE(gradient) == nsize + 6))
215 3166 : CALL pack_subsys_particles(subsys=subsys, f=gradient(1:nsize), fscale=-1.0_dp)
216 3166 : CALL apply_cell_change(gopt_env, cell, gradient, update_forces=.TRUE.)
217 3166 : IF (spgr%keep_space_group) THEN
218 20 : CALL spgr_apply_rotations_force(spgr, gradient)
219 20 : CALL spgr_apply_rotations_stress(spgr, cell, av_ptens)
220 20 : CALL spgr_write_stress_tensor(av_ptens, spgr)
221 : END IF
222 3166 : cell_gradient => gradient(nsize + 1:nsize + 6)
223 22162 : cell_gradient = 0.0_dp
224 : CALL get_dg_dh(cell_gradient, av_ptens, gopt_env%cell_env%pres_ext, cell, gopt_env%cell_env%mtrx, &
225 : keep_angles=gopt_env%cell_env%keep_angles, &
226 : keep_symmetry=gopt_env%cell_env%keep_symmetry, &
227 : pres_int=gopt_env%cell_env%pres_int, &
228 : pres_constr=gopt_env%cell_env%pres_constr, &
229 3166 : constraint_id=gopt_env%cell_env%constraint_id)
230 : END IF
231 : ! some callers expect pres_int to be available on all ranks. Also, here master is not necessarily a single rank.
232 : ! Assume at least master==0
233 4686 : CALL para_env%bcast(gopt_env%cell_env%pres_int, 0)
234 4686 : IF (gopt_env%cell_env%constraint_id /= fix_none) &
235 28 : CALL para_env%bcast(gopt_env%cell_env%pres_constr, 0)
236 : END IF
237 : END BLOCK
238 : CASE (default_cell_geo_opt_id, default_cell_md_id)
239 954 : CALL apply_cell_change(gopt_env, cell, x, update_forces=.FALSE.)
240 : ! Possibly output the new cell used for the next calculation
241 954 : CALL write_cell(cell, gopt_env%geo_section)
242 : ! Compute the pressure tensor
243 954 : BLOCK
244 : TYPE(virial_type) :: virial_avg
245 954 : IF (my_final_evaluation) THEN
246 : CALL force_env_calc_energy_force(gopt_env%force_env, &
247 : calc_force=PRESENT(gradient), &
248 24 : require_consistent_energy_force=gopt_env%require_consistent_energy_force)
249 24 : IF (PRESENT(f)) THEN
250 24 : CALL force_env_get(gopt_env%force_env, potential_energy=f)
251 : END IF
252 : ELSE
253 1840 : SELECT CASE (gopt_env%cell_method_id)
254 : CASE (default_cell_geo_opt_id)
255 910 : work => section_vals_get_subs_vals(gopt_env%motion_section, "GEO_OPT")
256 910 : CALL section_vals_get(work, explicit=explicit)
257 910 : IF (.NOT. explicit) &
258 : CALL cp_abort(__LOCATION__, &
259 0 : "Cell optimization at 0K was requested. GEO_OPT section MUST be provided in the input file!")
260 : ! Perform a geometry optimization
261 : CALL gopt_new_logger_create(new_logger, gopt_env%force_env%root_section, para_env, &
262 910 : project_name, id_run=geo_opt_run)
263 910 : CALL cp_add_default_logger(new_logger)
264 910 : CALL cp_geo_opt(gopt_env%force_env, gopt_env%globenv, eval_opt_geo=.FALSE.)
265 910 : CALL force_env_get(gopt_env%force_env, potential_energy=potential_energy)
266 910 : virial_avg = virial
267 : CASE (default_cell_md_id)
268 20 : work => section_vals_get_subs_vals(gopt_env%motion_section, "MD")
269 20 : avgs_section => section_vals_get_subs_vals(work, "AVERAGES")
270 20 : CALL section_vals_get(work, explicit=explicit)
271 20 : IF (.NOT. explicit) &
272 : CALL cp_abort( &
273 : __LOCATION__, &
274 0 : "Cell optimization at finite temperature was requested. MD section MUST be provided in the input file!")
275 : ! Only NVT ensemble is allowed..
276 20 : CALL section_vals_val_get(gopt_env%motion_section, "MD%ENSEMBLE", i_val=ensemble)
277 20 : IF (ensemble /= nvt_ensemble) &
278 : CALL cp_abort(__LOCATION__, &
279 0 : "Cell optimization at finite temperature requires the NVT MD ensemble!")
280 : ! Perform a molecular dynamics
281 : CALL gopt_new_logger_create(new_logger, gopt_env%force_env%root_section, para_env, &
282 20 : project_name, id_run=mol_dyn_run)
283 20 : CALL cp_add_default_logger(new_logger)
284 20 : CALL create_averages(averages, avgs_section, virial_avg=.TRUE., force_env=gopt_env%force_env)
285 20 : CALL qs_mol_dyn(gopt_env%force_env, gopt_env%globenv, averages, rm_restart_info=.FALSE.)
286 : ! Retrieve the average of the stress tensor and the average of the potential energy
287 20 : potential_energy = averages%avepot
288 20 : virial_avg = averages%virial
289 20 : CALL release_averages(averages)
290 : CASE DEFAULT
291 2790 : CPABORT("")
292 : END SELECT
293 930 : CALL cp_rm_default_logger()
294 : CALL gopt_new_logger_release(new_logger, gopt_env%force_env%root_section, para_env, project_name, &
295 930 : cell_opt_run)
296 : ! Update the virial
297 930 : CALL virial_update(virial_avg, subsys, para_env)
298 : ! Possibly take give back the potential energy
299 930 : IF (PRESENT(f)) THEN
300 930 : f = potential_energy
301 : END IF
302 : END IF
303 : ! Possibly give back the gradients
304 236592 : IF (PRESENT(gradient)) THEN
305 930 : CPASSERT(ANY(virial_avg%pv_total /= 0))
306 : ! Convert the average ptens
307 12090 : av_ptens(:, :) = virial_avg%pv_total(:, :)/cell%deth
308 930 : IF (master == para_env%mepos) THEN ! we are on the master
309 904 : CPASSERT(ASSOCIATED(gradient))
310 904 : IF (spgr%keep_space_group) THEN
311 0 : CALL spgr_apply_rotations_stress(spgr, cell, av_ptens)
312 0 : CALL spgr_write_stress_tensor(av_ptens, spgr)
313 : END IF
314 : ! Compute the gradients on the cell
315 : CALL get_dg_dh(gradient, av_ptens, gopt_env%cell_env%pres_ext, cell, gopt_env%cell_env%mtrx, &
316 : keep_angles=gopt_env%cell_env%keep_angles, &
317 : keep_symmetry=gopt_env%cell_env%keep_symmetry, &
318 : pres_int=gopt_env%cell_env%pres_int, &
319 : pres_constr=gopt_env%cell_env%pres_constr, &
320 904 : constraint_id=gopt_env%cell_env%constraint_id)
321 : END IF
322 : ! some callers expect pres_int to be available on all ranks. Also, here master is not necessarily a single rank.
323 : ! Assume at least master==0
324 930 : CALL para_env%bcast(gopt_env%cell_env%pres_int, 0)
325 930 : IF (gopt_env%cell_env%constraint_id /= fix_none) &
326 0 : CALL para_env%bcast(gopt_env%cell_env%pres_constr, 0)
327 : END IF
328 : END BLOCK
329 : CASE DEFAULT
330 6132 : CPABORT("")
331 : END SELECT
332 : CASE (default_shellcore_method_id)
333 340 : idg = 0
334 32980 : DO ip = 1, particles%n_els
335 32640 : shell_index = particles%els(ip)%shell_index
336 32980 : IF (shell_index /= 0) THEN
337 130560 : DO idir = 1, 3
338 97920 : idg = 3*(shell_index - 1) + idir
339 130560 : shell_particles%els(shell_index)%r(idir) = core_particles%els(ip)%r(idir) - x(idg)
340 : END DO
341 : END IF
342 : END DO
343 340 : CALL write_structure_data(particles%els, cell, gopt_env%motion_section)
344 :
345 : ! Shell-core optimization
346 : CALL force_env_calc_energy_force(gopt_env%force_env, &
347 : calc_force=PRESENT(gradient), &
348 340 : require_consistent_energy_force=gopt_env%require_consistent_energy_force)
349 :
350 : ! Possibly take the potential energy
351 340 : IF (PRESENT(f)) THEN
352 340 : CALL force_env_get(gopt_env%force_env, potential_energy=f)
353 : END IF
354 :
355 : ! Possibly take the gradients
356 340 : IF (PRESENT(gradient)) THEN
357 340 : IF (master == para_env%mepos) THEN ! we are on the master
358 340 : CPASSERT(ASSOCIATED(gradient))
359 340 : idg = 0
360 32980 : DO ip = 1, shell_particles%n_els
361 130900 : DO idir = 1, 3
362 97920 : idg = idg + 1
363 130560 : gradient(idg) = -(core_particles%els(ip)%f(idir) - shell_particles%els(ip)%f(idir))
364 : END DO
365 : END DO
366 : END IF
367 : END IF
368 : CASE DEFAULT
369 20342 : CPABORT("")
370 : END SELECT
371 :
372 20342 : CALL timestop(handle)
373 :
374 40684 : END SUBROUTINE cp_eval_at
|
