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