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 16539 : 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_method_id, &
53 : default_shellcore_method_id, &
54 : fix_none
55 : USE input_section_types, ONLY: section_vals_get, &
56 : section_vals_get_subs_vals, &
57 : section_vals_type, &
58 : section_vals_val_get
59 : USE md_run, ONLY: qs_mol_dyn
60 : USE kinds, ONLY: dp
61 : USE particle_list_types, ONLY: particle_list_type
62 : USE particle_methods, ONLY: write_structure_data
63 : USE virial_methods, ONLY: virial_update
64 : USE virial_types, ONLY: virial_type
65 : USE cp_log_handling, ONLY: cp_add_default_logger, &
66 : cp_rm_default_logger
67 : USE space_groups_types, ONLY: spgr_type
68 : USE space_groups, ONLY: spgr_apply_rotations_stress, &
69 : spgr_apply_rotations_coord, &
70 : spgr_apply_rotations_force, &
71 : spgr_write_stress_tensor
72 :
73 : #include "../base/base_uses.f90"
74 : IMPLICIT NONE
75 : TYPE(gopt_f_type), POINTER :: gopt_env
76 : REAL(KIND=dp), DIMENSION(:), POINTER :: x
77 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: f
78 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, &
79 : POINTER :: gradient
80 : INTEGER, INTENT(IN) :: master
81 : LOGICAL, INTENT(IN), OPTIONAL :: final_evaluation
82 : TYPE(mp_para_env_type), POINTER :: para_env
83 :
84 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_eval_at'
85 :
86 : INTEGER :: handle, idg, idir, ip, &
87 : nparticle, nsize, shell_index
88 : LOGICAL :: my_final_evaluation
89 : REAL(KIND=dp) :: f_ts, potential_energy
90 : REAL(KIND=dp), DIMENSION(3, 3) :: av_ptens
91 16539 : REAL(KIND=dp), DIMENSION(:), POINTER :: cell_gradient, gradient_ts
92 : TYPE(cell_type), POINTER :: cell
93 : TYPE(cp_subsys_type), POINTER :: subsys
94 : TYPE(particle_list_type), POINTER :: core_particles, particles, &
95 : shell_particles
96 : TYPE(virial_type), POINTER :: virial
97 : TYPE(cp_logger_type), POINTER :: new_logger
98 : TYPE(average_quantities_type), POINTER :: averages
99 : TYPE(spgr_type), POINTER :: spgr
100 :
101 16539 : NULLIFY (averages)
102 16539 : NULLIFY (cell)
103 16539 : NULLIFY (core_particles)
104 16539 : NULLIFY (gradient_ts)
105 16539 : NULLIFY (particles)
106 16539 : NULLIFY (shell_particles)
107 16539 : NULLIFY (subsys)
108 16539 : NULLIFY (virial)
109 16539 : NULLIFY (new_logger)
110 16539 : NULLIFY (spgr)
111 :
112 16539 : CALL timeset(routineN, handle)
113 :
114 16539 : CALL force_env_get(gopt_env%force_env, subsys=subsys, cell=cell)
115 : CALL cp_subsys_get(subsys, &
116 : core_particles=core_particles, &
117 : particles=particles, &
118 : shell_particles=shell_particles, &
119 16539 : virial=virial)
120 :
121 16539 : spgr => gopt_env%spgr
122 :
123 16539 : my_final_evaluation = .FALSE.
124 16539 : IF (PRESENT(final_evaluation)) my_final_evaluation = final_evaluation
125 :
126 24810 : SELECT CASE (gopt_env%type_id)
127 : CASE (default_minimization_method_id, default_ts_method_id)
128 8271 : CALL unpack_subsys_particles(subsys=subsys, r=x)
129 8271 : CALL write_structure_data(particles%els, cell, gopt_env%motion_section)
130 14726 : SELECT CASE (gopt_env%type_id)
131 : CASE (default_minimization_method_id)
132 : ! Geometry Minimization
133 : CALL force_env_calc_energy_force(gopt_env%force_env, &
134 : calc_force=PRESENT(gradient), &
135 6455 : require_consistent_energy_force=gopt_env%require_consistent_energy_force)
136 : ! Possibly take the potential energy
137 6455 : IF (PRESENT(f)) THEN
138 6455 : CALL force_env_get(gopt_env%force_env, potential_energy=f)
139 : END IF
140 : ! Possibly take the gradients
141 6455 : IF (PRESENT(gradient)) THEN
142 5905 : IF (master == para_env%mepos) THEN ! we are on the master
143 5754 : CALL pack_subsys_particles(subsys=subsys, f=gradient, fscale=-1.0_dp)
144 5754 : IF (spgr%keep_space_group) THEN
145 8 : CALL spgr_apply_rotations_force(spgr, gradient)
146 8 : CALL unpack_subsys_particles(subsys=subsys, f=gradient, fscale=-1.0_dp)
147 : END IF
148 : END IF
149 : END IF
150 : CASE (default_ts_method_id)
151 : ! Transition State Optimization
152 5448 : ALLOCATE (gradient_ts(particles%n_els*3))
153 : ! Real calculation of energy and forces for transition state optimization:
154 : ! When doing dimer methods forces have to be always computed since the function
155 : ! to minimize is not the energy but the effective force
156 1816 : CALL cp_eval_at_ts(gopt_env, x, f_ts, gradient_ts, calc_force=.TRUE.)
157 1816 : CALL cite_reference(Henkelman1999)
158 : ! Possibly take the potential energy
159 1816 : IF (PRESENT(f)) f = f_ts
160 : ! Possibly take the gradients
161 1816 : IF (PRESENT(gradient)) THEN
162 854 : IF (master == para_env%mepos) THEN ! we are on the master
163 854 : CPASSERT(ASSOCIATED(gradient))
164 29452 : gradient = gradient_ts
165 : END IF
166 : END IF
167 10087 : DEALLOCATE (gradient_ts)
168 : END SELECT
169 : ! This call is necessary for QM/MM if a Translation is applied
170 : ! this makes the geometry optimizer consistent
171 8271 : CALL unpack_subsys_particles(subsys=subsys, r=x)
172 : CASE (default_cell_method_id)
173 : ! Check for VIRIAL
174 7928 : IF (.NOT. virial%pv_availability) &
175 : CALL cp_abort(__LOCATION__, &
176 : "For a cell optimization task with CELL_OPT/TYPE "// &
177 : "DIRECT_CELL_OPT, the FORCE_EVAL/STRESS_TENSOR "// &
178 : "keyword MUST be defined in the input file for the "// &
179 0 : "evaluation of the stress tensor, but none is found!")
180 7928 : IF (gopt_env%cell_env%keep_volume) THEN
181 2278 : nparticle = force_env_get_nparticle(gopt_env%force_env)
182 2278 : idg = 3*nparticle
183 2278 : CALL rescale_new_cell_volume(cell%deth, x, idg)
184 : END IF
185 :
186 7928 : CALL apply_cell_change(gopt_env, cell, x, update_forces=.FALSE.)
187 : ! Possibly output the new cell used for the next calculation
188 7928 : CALL write_cell(cell, gopt_env%geo_section)
189 : ! Compute the pressure tensor
190 7928 : BLOCK
191 : TYPE(virial_type) :: virial_avg
192 : CALL force_env_calc_energy_force(gopt_env%force_env, &
193 : calc_force=PRESENT(gradient), &
194 7928 : require_consistent_energy_force=gopt_env%require_consistent_energy_force)
195 : ! Possibly take the potential energy
196 7928 : CALL force_env_get(gopt_env%force_env, potential_energy=potential_energy)
197 7928 : virial_avg = virial
198 7928 : CALL virial_update(virial_avg, subsys, para_env)
199 7928 : IF (PRESENT(f)) THEN
200 7928 : CALL force_env_get(gopt_env%force_env, potential_energy=f)
201 : END IF
202 : ! Possibly take the gradients
203 1966144 : IF (PRESENT(gradient)) THEN
204 6606 : CPASSERT(ANY(virial_avg%pv_total /= 0))
205 : ! Convert the average ptens
206 85878 : av_ptens(:, :) = virial_avg%pv_total(:, :)/cell%deth
207 6606 : IF (master == para_env%mepos) THEN ! we are on the master
208 5062 : CPASSERT(ASSOCIATED(gradient))
209 5062 : nparticle = force_env_get_nparticle(gopt_env%force_env)
210 5062 : nsize = 3*nparticle
211 5062 : CPASSERT((SIZE(gradient) == nsize + 6))
212 5062 : CALL pack_subsys_particles(subsys=subsys, f=gradient(1:nsize), fscale=-1.0_dp)
213 5062 : CALL apply_cell_change(gopt_env, cell, gradient, update_forces=.TRUE.)
214 5062 : IF (spgr%keep_space_group) THEN
215 544 : CALL spgr_apply_rotations_force(spgr, gradient)
216 544 : CALL spgr_apply_rotations_stress(spgr, cell, av_ptens)
217 544 : CALL spgr_write_stress_tensor(av_ptens, spgr)
218 : END IF
219 5062 : cell_gradient => gradient(nsize + 1:nsize + 6)
220 35434 : cell_gradient = 0.0_dp
221 : CALL get_dg_dh(cell_gradient, av_ptens, gopt_env%cell_env%pres_ext, cell, gopt_env%cell_env%mtrx, &
222 : keep_angles=gopt_env%cell_env%keep_angles, &
223 : keep_symmetry=gopt_env%cell_env%keep_symmetry, &
224 : pres_int=gopt_env%cell_env%pres_int, &
225 : pres_constr=gopt_env%cell_env%pres_constr, &
226 5062 : constraint_id=gopt_env%cell_env%constraint_id)
227 : END IF
228 : ! some callers expect pres_int to be available on all ranks. Also, here master is not necessarily a single rank.
229 : ! Assume at least master==0
230 6606 : CALL para_env%bcast(gopt_env%cell_env%pres_int, 0)
231 6606 : IF (gopt_env%cell_env%constraint_id /= fix_none) &
232 24 : CALL para_env%bcast(gopt_env%cell_env%pres_constr, 0)
233 : END IF
234 : END BLOCK
235 : CASE (default_shellcore_method_id)
236 : idg = 0
237 32980 : DO ip = 1, particles%n_els
238 32640 : shell_index = particles%els(ip)%shell_index
239 32980 : IF (shell_index /= 0) THEN
240 130560 : DO idir = 1, 3
241 97920 : idg = 3*(shell_index - 1) + idir
242 130560 : shell_particles%els(shell_index)%r(idir) = core_particles%els(ip)%r(idir) - x(idg)
243 : END DO
244 : END IF
245 : END DO
246 340 : CALL write_structure_data(particles%els, cell, gopt_env%motion_section)
247 :
248 : ! Shell-core optimization
249 : CALL force_env_calc_energy_force(gopt_env%force_env, &
250 : calc_force=PRESENT(gradient), &
251 340 : require_consistent_energy_force=gopt_env%require_consistent_energy_force)
252 :
253 : ! Possibly take the potential energy
254 340 : IF (PRESENT(f)) THEN
255 340 : CALL force_env_get(gopt_env%force_env, potential_energy=f)
256 : END IF
257 :
258 : ! Possibly take the gradients
259 340 : IF (PRESENT(gradient)) THEN
260 340 : IF (master == para_env%mepos) THEN ! we are on the master
261 340 : CPASSERT(ASSOCIATED(gradient))
262 340 : idg = 0
263 32980 : DO ip = 1, shell_particles%n_els
264 130900 : DO idir = 1, 3
265 97920 : idg = idg + 1
266 130560 : gradient(idg) = -(core_particles%els(ip)%f(idir) - shell_particles%els(ip)%f(idir))
267 : END DO
268 : END DO
269 : END IF
270 : END IF
271 : CASE DEFAULT
272 16539 : CPABORT("Invalid or not yet implemented type of optimization")
273 : END SELECT
274 :
275 16539 : CALL timestop(handle)
276 :
277 33078 : END SUBROUTINE cp_eval_at
|
