1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Module with utility for Nudged Elastic Band Calculation 8!> \note 9!> Numerical accuracy for parallel runs: 10!> Each replica starts the SCF run from the one optimized 11!> in a previous run. It may happen then energies and derivatives 12!> of a serial run and a parallel run could be slightly different 13!> 'cause of a different starting density matrix. 14!> Exact results are obtained using: 15!> EXTRAPOLATION USE_GUESS in QS section (Teo 09.2006) 16!> \author Teodoro Laino 10.2006 17! ************************************************************************************************** 18MODULE neb_utils 19 USE bibliography, ONLY: E2002,& 20 Elber1987,& 21 Jonsson1998,& 22 Jonsson2000_1,& 23 Jonsson2000_2,& 24 Wales2004,& 25 cite_reference 26 USE colvar_utils, ONLY: eval_colvar,& 27 get_clv_force,& 28 set_colvars_target 29 USE cp_log_handling, ONLY: cp_get_default_logger,& 30 cp_logger_type 31 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 32 cp_print_key_unit_nr 33 USE cp_para_types, ONLY: cp_para_env_type 34 USE cp_parser_methods, ONLY: parser_get_next_line,& 35 parser_get_object 36 USE cp_parser_types, ONLY: cp_parser_type,& 37 parser_create,& 38 parser_release 39 USE f77_interface, ONLY: f_env_add_defaults,& 40 f_env_rm_defaults,& 41 f_env_type,& 42 get_energy,& 43 get_force,& 44 get_pos,& 45 set_pos 46 USE force_env_methods, ONLY: force_env_calc_energy_force 47 USE force_env_types, ONLY: force_env_get 48 USE geo_opt, ONLY: cp_geo_opt 49 USE global_types, ONLY: global_environment_type 50 USE input_constants, ONLY: & 51 band_diis_opt, do_b_neb, do_band_cartesian, do_ci_neb, do_d_neb, do_eb, do_it_neb, do_sm, & 52 pot_neb_fe, pot_neb_full, pot_neb_me 53 USE input_cp2k_check, ONLY: remove_restart_info 54 USE input_section_types, ONLY: section_vals_get,& 55 section_vals_get_subs_vals,& 56 section_vals_type,& 57 section_vals_val_get 58 USE kinds, ONLY: default_path_length,& 59 default_string_length,& 60 dp 61 USE mathlib, ONLY: matvec_3x3 62 USE md_run, ONLY: qs_mol_dyn 63 USE neb_io, ONLY: dump_replica_coordinates,& 64 handle_band_file_names 65 USE neb_md_utils, ONLY: neb_initialize_velocity 66 USE neb_types, ONLY: neb_type,& 67 neb_var_type 68 USE particle_types, ONLY: particle_type 69 USE physcon, ONLY: bohr 70 USE replica_methods, ONLY: rep_env_calc_e_f 71 USE replica_types, ONLY: rep_env_sync,& 72 replica_env_type 73 USE rmsd, ONLY: rmsd3 74#include "../base/base_uses.f90" 75 76 IMPLICIT NONE 77 PRIVATE 78 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'neb_utils' 79 LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE. 80 81 PUBLIC :: build_replica_coords, & 82 neb_calc_energy_forces, & 83 reorient_images, & 84 reparametrize_images, & 85 check_convergence 86 87CONTAINS 88 89! ************************************************************************************************** 90!> \brief Computes the distance between two replica 91!> \param particle_set ... 92!> \param coords ... 93!> \param i0 ... 94!> \param i ... 95!> \param distance ... 96!> \param iw ... 97!> \param rotate ... 98!> \author Teodoro Laino 09.2006 99! ************************************************************************************************** 100 SUBROUTINE neb_replica_distance(particle_set, coords, i0, i, distance, iw, rotate) 101 TYPE(particle_type), DIMENSION(:), OPTIONAL, & 102 POINTER :: particle_set 103 TYPE(neb_var_type), POINTER :: coords 104 INTEGER, INTENT(IN) :: i0, i 105 REAL(KIND=dp), INTENT(OUT) :: distance 106 INTEGER, INTENT(IN) :: iw 107 LOGICAL, INTENT(IN), OPTIONAL :: rotate 108 109 CHARACTER(len=*), PARAMETER :: routineN = 'neb_replica_distance', & 110 routineP = moduleN//':'//routineN 111 112 LOGICAL :: my_rotate 113 114 my_rotate = .FALSE. 115 IF (PRESENT(rotate)) my_rotate = rotate 116 ! The rotation of the replica is enabled exclusively when working in 117 ! cartesian coordinates 118 IF (my_rotate .AND. (coords%in_use == do_band_cartesian)) THEN 119 CPASSERT(PRESENT(particle_set)) 120 CALL rmsd3(particle_set, coords%xyz(:, i), coords%xyz(:, i0), & 121 iw, rotate=my_rotate) 122 END IF 123 distance = SQRT(DOT_PRODUCT(coords%wrk(:, i) - coords%wrk(:, i0), & 124 coords%wrk(:, i) - coords%wrk(:, i0))) 125 126 END SUBROUTINE neb_replica_distance 127 128! ************************************************************************************************** 129!> \brief Constructs or Read the coordinates for all replica 130!> \param neb_section ... 131!> \param particle_set ... 132!> \param coords ... 133!> \param vels ... 134!> \param neb_env ... 135!> \param iw ... 136!> \param globenv ... 137!> \param para_env ... 138!> \author Teodoro Laino 09.2006 139! ************************************************************************************************** 140 SUBROUTINE build_replica_coords(neb_section, particle_set, & 141 coords, vels, neb_env, iw, globenv, para_env) 142 TYPE(section_vals_type), POINTER :: neb_section 143 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 144 TYPE(neb_var_type), POINTER :: coords, vels 145 TYPE(neb_type), POINTER :: neb_env 146 INTEGER, INTENT(IN) :: iw 147 TYPE(global_environment_type), POINTER :: globenv 148 TYPE(cp_para_env_type), POINTER :: para_env 149 150 CHARACTER(len=*), PARAMETER :: routineN = 'build_replica_coords', & 151 routineP = moduleN//':'//routineN 152 153 CHARACTER(LEN=default_path_length) :: filename 154 CHARACTER(LEN=default_string_length) :: dummy_char 155 INTEGER :: handle, i_rep, iatom, ic, input_nr_replica, is, ivar, j, jtarg, n_rep, natom, & 156 neb_nr_replica, nr_replica_to_interpolate, nval, nvar, shell_index 157 LOGICAL :: check, explicit, my_end, skip_vel_section 158 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: distance 159 REAL(KIND=dp), DIMENSION(3) :: r 160 REAL(KIND=dp), DIMENSION(:), POINTER :: initial_colvars, rptr 161 TYPE(cp_parser_type), POINTER :: parser 162 TYPE(section_vals_type), POINTER :: coord_section, replica_section, & 163 vel_section 164 165 NULLIFY (parser) 166 CALL timeset(routineN, handle) 167 CPASSERT(ASSOCIATED(coords)) 168 CPASSERT(ASSOCIATED(vels)) 169 neb_nr_replica = neb_env%number_of_replica 170 replica_section => section_vals_get_subs_vals(neb_section, "REPLICA") 171 CALL section_vals_get(replica_section, n_repetition=input_nr_replica) 172 ! Calculation is aborted if input replicas are more then the requested ones for the BAND.. 173 CPASSERT(input_nr_replica <= neb_nr_replica) 174 ! Read in replicas coordinates 175 skip_vel_section = (input_nr_replica /= neb_nr_replica) 176 IF ((iw > 0) .AND. skip_vel_section) THEN 177 WRITE (iw, '(T2,A)') 'NEB| The number of replica in the input is different from the number', & 178 'NEB| of replica requested for NEB. More Replica will be interpolated.', & 179 'NEB| Therefore the possibly provided velocities will not be read.' 180 END IF 181 ! Further check on velocity section... 182 DO i_rep = 1, input_nr_replica 183 vel_section => section_vals_get_subs_vals(replica_section, "VELOCITY", & 184 i_rep_section=i_rep) 185 CALL section_vals_get(vel_section, explicit=explicit) 186 skip_vel_section = skip_vel_section .OR. (.NOT. explicit) 187 END DO 188 ! Setup cartesian coordinates and COLVAR (if requested) 189 coords%xyz(:, :) = 0.0_dp 190 DO i_rep = 1, input_nr_replica 191 coord_section => section_vals_get_subs_vals(replica_section, "COORD", & 192 i_rep_section=i_rep) 193 CALL section_vals_get(coord_section, explicit=explicit) 194 ! Cartesian Coordinates 195 IF (explicit) THEN 196 CALL section_vals_val_get(coord_section, "_DEFAULT_KEYWORD_", & 197 n_rep_val=natom) 198 CPASSERT((natom == SIZE(particle_set))) 199 DO iatom = 1, natom 200 CALL section_vals_val_get(coord_section, "_DEFAULT_KEYWORD_", & 201 i_rep_val=iatom, r_vals=rptr) 202 ic = 3*(iatom - 1) 203 coords%xyz(ic + 1:ic + 3, i_rep) = rptr(1:3)*bohr 204 ! Initially core and shell positions are set to the atomic positions 205 shell_index = particle_set(iatom)%shell_index 206 IF (shell_index /= 0) THEN 207 is = 3*(natom + shell_index - 1) 208 coords%xyz(is + 1:is + 3, i_rep) = coords%xyz(ic + 1:ic + 3, i_rep) 209 END IF 210 END DO 211 ELSE 212 CALL section_vals_val_get(replica_section, "COORD_FILE_NAME", & 213 i_rep_section=i_rep, c_val=filename) 214 CPASSERT(TRIM(filename) /= "") 215 CALL parser_create(parser, filename, para_env=para_env, parse_white_lines=.TRUE.) 216 CALL parser_get_next_line(parser, 1) 217 ! Start parser 218 CALL parser_get_object(parser, natom) 219 CPASSERT((natom == SIZE(particle_set))) 220 CALL parser_get_next_line(parser, 1) 221 DO iatom = 1, natom 222 ! Atom coordinates 223 CALL parser_get_next_line(parser, 1, at_end=my_end) 224 IF (my_end) & 225 CALL cp_abort(__LOCATION__, & 226 "Number of lines in XYZ format not equal to the number of atoms."// & 227 " Error in XYZ format for REPLICA coordinates. Very probably the"// & 228 " line with title is missing or is empty. Please check the XYZ file and rerun your job!") 229 READ (parser%input_line, *) dummy_char, r(1:3) 230 ic = 3*(iatom - 1) 231 coords%xyz(ic + 1:ic + 3, i_rep) = r(1:3)*bohr 232 ! Initially core and shell positions are set to the atomic positions 233 shell_index = particle_set(iatom)%shell_index 234 IF (shell_index /= 0) THEN 235 is = 3*(natom + shell_index - 1) 236 coords%xyz(is + 1:is + 3, i_rep) = coords%xyz(ic + 1:ic + 3, i_rep) 237 END IF 238 END DO 239 CALL parser_release(parser) 240 END IF 241 ! Collective Variables 242 IF (neb_env%use_colvar) THEN 243 CALL section_vals_val_get(replica_section, "COLLECTIVE", & 244 i_rep_section=i_rep, n_rep_val=n_rep) 245 IF (n_rep /= 0) THEN 246 ! Read the values of the collective variables 247 NULLIFY (initial_colvars) 248 CALL section_vals_val_get(replica_section, "COLLECTIVE", & 249 i_rep_section=i_rep, r_vals=initial_colvars) 250 check = (neb_env%nsize_int == SIZE(initial_colvars)) 251 CPASSERT(check) 252 coords%int(:, i_rep) = initial_colvars 253 ELSE 254 ! Compute the values of the collective variables 255 CALL eval_colvar(neb_env%force_env, coords%xyz(:, i_rep), coords%int(:, i_rep)) 256 END IF 257 END IF 258 ! Dump cartesian and colvar info.. 259 CALL dump_replica_coordinates(particle_set, coords, i_rep, i_rep, iw, neb_env%use_colvar) 260 ! Setup Velocities 261 IF (skip_vel_section) THEN 262 CALL neb_initialize_velocity(vels%wrk, neb_section, particle_set, & 263 i_rep, iw, globenv, neb_env) 264 ELSE 265 vel_section => section_vals_get_subs_vals(replica_section, "VELOCITY", & 266 i_rep_section=i_rep) 267 CALL section_vals_val_get(vel_section, "_DEFAULT_KEYWORD_", & 268 n_rep_val=nval) 269 ! Setup Velocities for collective or cartesian coordinates 270 IF (neb_env%use_colvar) THEN 271 nvar = SIZE(vels%wrk, 1) 272 CPASSERT(nval == nvar) 273 DO ivar = 1, nvar 274 CALL section_vals_val_get(vel_section, "_DEFAULT_KEYWORD_", & 275 i_rep_val=ivar, r_vals=rptr) 276 vels%wrk(ivar, i_rep) = rptr(1) 277 END DO 278 ELSE 279 natom = SIZE(particle_set) 280 CPASSERT(nval == natom) 281 DO iatom = 1, natom 282 CALL section_vals_val_get(vel_section, "_DEFAULT_KEYWORD_", & 283 i_rep_val=iatom, r_vals=rptr) 284 ic = 3*(iatom - 1) 285 vels%wrk(ic + 1:ic + 3, i_rep) = rptr(1:3) 286 ! Initially set shell velocities to core velocity 287 shell_index = particle_set(iatom)%shell_index 288 IF (shell_index /= 0) THEN 289 is = 3*(natom + shell_index - 1) 290 vels%wrk(is + 1:is + 3, i_rep) = vels%wrk(ic + 1:ic + 3, i_rep) 291 END IF 292 END DO 293 END IF 294 END IF 295 END DO ! i_rep 296 ALLOCATE (distance(neb_nr_replica - 1)) 297 IF (input_nr_replica < neb_nr_replica) THEN 298 ! Interpolate missing replicas 299 nr_replica_to_interpolate = neb_nr_replica - input_nr_replica 300 distance = 0.0_dp 301 IF (iw > 0) THEN 302 WRITE (iw, '(T2,A,I0,A)') 'NEB| Interpolating ', nr_replica_to_interpolate, ' missing Replica.' 303 END IF 304 DO WHILE (nr_replica_to_interpolate > 0) 305 ! Compute distances between known images to find the interval 306 ! where to add a new image 307 DO j = 1, input_nr_replica - 1 308 CALL neb_replica_distance(particle_set, coords, j, j + 1, distance(j), iw, & 309 rotate=neb_env%align_frames) 310 END DO 311 jtarg = MAXLOC(distance(1:input_nr_replica), 1) 312 IF (iw > 0) THEN 313 WRITE (iw, '(T2,3(A,I0),A)') 'NEB| Interpolating Nr. ', & 314 nr_replica_to_interpolate, ' missing Replica between Replica Nr. ', & 315 jtarg, ' and ', jtarg + 1, '.' 316 END IF 317 input_nr_replica = input_nr_replica + 1 318 nr_replica_to_interpolate = nr_replica_to_interpolate - 1 319 ! Interpolation is a simple bisection in XYZ 320 coords%xyz(:, jtarg + 2:input_nr_replica) = coords%xyz(:, jtarg + 1:input_nr_replica - 1) 321 coords%xyz(:, jtarg + 1) = (coords%xyz(:, jtarg) + coords%xyz(:, jtarg + 2))/2.0_dp 322 IF (neb_env%use_colvar) THEN 323 ! Interpolation is a simple bisection also in internal coordinates 324 ! in this case the XYZ coordinates need only as a starting point for computing 325 ! the potential energy function. The reference are the internal coordinates 326 ! interpolated here after.. 327 coords%int(:, jtarg + 2:input_nr_replica) = coords%int(:, jtarg + 1:input_nr_replica - 1) 328 coords%int(:, jtarg + 1) = (coords%int(:, jtarg) + coords%int(:, jtarg + 2))/2.0_dp 329 END IF 330 vels%wrk(:, jtarg + 2:input_nr_replica) = vels%wrk(:, jtarg + 1:input_nr_replica - 1) 331 vels%wrk(:, jtarg + 1) = 0.0_dp 332 CALL dump_replica_coordinates(particle_set, coords, jtarg + 1, & 333 input_nr_replica, iw, neb_env%use_colvar) 334 CALL neb_initialize_velocity(vels%wrk, neb_section, particle_set, & 335 jtarg + 1, iw, globenv, neb_env) 336 END DO 337 END IF 338 vels%wrk(:, 1) = 0.0_dp 339 vels%wrk(:, neb_nr_replica) = 0.0_dp 340 ! If we perform a DIIS optimization we don't need velocities 341 IF (neb_env%opt_type == band_diis_opt) vels%wrk = 0.0_dp 342 ! Compute distances between replicas and in case of Cartesian Coordinates 343 ! Rotate the frames in order to minimize the RMSD 344 DO j = 1, input_nr_replica - 1 345 CALL neb_replica_distance(particle_set, coords, j, j + 1, distance(j), iw, & 346 rotate=neb_env%align_frames) 347 END DO 348 DEALLOCATE (distance) 349 350 CALL timestop(handle) 351 352 END SUBROUTINE build_replica_coords 353 354! ************************************************************************************************** 355!> \brief Driver to compute energy and forces within a NEB, 356!> Based on the use of the replica_env 357!> \param rep_env ... 358!> \param neb_env ... 359!> \param coords ... 360!> \param energies ... 361!> \param forces ... 362!> \param particle_set ... 363!> \param output_unit ... 364!> \author Teodoro Laino 09.2006 365! ************************************************************************************************** 366 SUBROUTINE neb_calc_energy_forces(rep_env, neb_env, coords, energies, forces, & 367 particle_set, output_unit) 368 TYPE(replica_env_type), POINTER :: rep_env 369 TYPE(neb_type), OPTIONAL, POINTER :: neb_env 370 TYPE(neb_var_type), POINTER :: coords 371 REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: energies 372 TYPE(neb_var_type), POINTER :: forces 373 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 374 INTEGER, INTENT(IN) :: output_unit 375 376 CHARACTER(len=*), PARAMETER :: routineN = 'neb_calc_energy_forces', & 377 routineP = moduleN//':'//routineN 378 CHARACTER(LEN=1), DIMENSION(3), PARAMETER :: lab = (/"X", "Y", "Z"/) 379 380 INTEGER :: handle, i, irep, j, n_int, n_rep, & 381 n_rep_neb, nsize_wrk 382 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: tangent, tmp_a, tmp_b 383 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: cvalues, Mmatrix, Mmatrix_tmp 384 385 CALL timeset(routineN, handle) 386 n_int = neb_env%nsize_int 387 n_rep_neb = neb_env%number_of_replica 388 n_rep = rep_env%nrep 389 nsize_wrk = coords%size_wrk(1) 390 energies = 0.0_dp 391 ALLOCATE (cvalues(n_int, n_rep)) 392 ALLOCATE (Mmatrix_tmp(n_int*n_int, n_rep)) 393 ALLOCATE (Mmatrix(n_int*n_int, n_rep_neb)) 394 IF (output_unit > 0) WRITE (output_unit, '(/,T2,A)') "NEB| Computing Energies and Forces" 395 DO irep = 1, n_rep_neb, n_rep 396 DO j = 0, n_rep - 1 397 IF (irep + j <= n_rep_neb) THEN 398 ! If the number of replica in replica_env and the number of replica 399 ! used in the NEB does not match, the other replica in replica_env 400 ! just compute energies and forces keeping the fixed coordinates and 401 ! forces 402 rep_env%r(:, j + 1) = coords%xyz(:, irep + j) 403 END IF 404 END DO 405 ! Fix file name for BAND replicas.. Each BAND replica has its own file 406 ! independently from the number of replicas in replica_env.. 407 CALL handle_band_file_names(rep_env, irep, n_rep_neb, neb_env%istep) 408 ! Let's select the potential we want to use for the band calculation 409 SELECT CASE (neb_env%pot_type) 410 CASE (pot_neb_full) 411 ! Full potential Energy 412 CALL rep_env_calc_e_f(rep_env, calc_f=.TRUE.) 413 CASE (pot_neb_fe) 414 ! Free Energy Case 415 CALL perform_replica_md(rep_env, coords, irep, n_rep_neb, cvalues, Mmatrix_tmp) 416 CASE (pot_neb_me) 417 ! Minimum Potential Energy Case 418 CALL perform_replica_geo(rep_env, coords, irep, n_rep_neb, cvalues, Mmatrix_tmp) 419 END SELECT 420 421 DO j = 0, n_rep - 1 422 IF (irep + j <= n_rep_neb) THEN 423 ! Copy back Forces and Energies 424 forces%wrk(:, irep + j) = rep_env%f(1:nsize_wrk, j + 1) 425 energies(irep + j) = rep_env%f(rep_env%ndim + 1, j + 1) 426 SELECT CASE (neb_env%pot_type) 427 CASE (pot_neb_full) 428 ! Dump Info 429 IF (output_unit > 0) THEN 430 WRITE (output_unit, '(T2,A,I5,A,I5,A)') & 431 "NEB| REPLICA Nr.", irep + j, "- Energy and Forces" 432 WRITE (output_unit, '(T2,A,T42,A,9X,F15.9)') & 433 "NEB|", " Total Energy: ", rep_env%f(rep_env%ndim + 1, j + 1) 434 WRITE (output_unit, '(T2,"NEB|",T10,"ATOM",T33,3(9X,A,7X))') lab(1), lab(2), lab(3) 435 DO i = 1, SIZE(particle_set) 436 WRITE (output_unit, '(T2,"NEB|",T12,A,T30,3(2X,F15.9))') & 437 particle_set(i)%atomic_kind%name, & 438 rep_env%f((i - 1)*3 + 1:(i - 1)*3 + 3, j + 1) 439 END DO 440 END IF 441 CASE (pot_neb_fe, pot_neb_me) 442 ! Let's update the cartesian coordinates. This will make 443 ! easier the next evaluation of energies and forces... 444 coords%xyz(:, irep + j) = rep_env%r(1:rep_env%ndim, j + 1) 445 Mmatrix(:, irep + j) = Mmatrix_tmp(:, j + 1) 446 IF (output_unit > 0) THEN 447 WRITE (output_unit, '(/,T2,A,I5,A,I5,A)') & 448 "NEB| REPLICA Nr.", irep + j, "- Energy, Collective Variables, Forces" 449 WRITE (output_unit, '(T2,A,T42,A,9X,F15.9)') & 450 "NEB|", " Total Energy: ", rep_env%f(rep_env%ndim + 1, j + 1) 451 WRITE (output_unit, & 452 '(T2,"NEB|",T10,"CV Nr.",12X,"Expected COLVAR",5X,"Present COLVAR",10X,"Forces")') 453 DO i = 1, n_int 454 WRITE (output_unit, '(T2,"NEB|",T12,I2,7X,3(5X,F15.9))') & 455 i, coords%int(i, irep + j), cvalues(i, j + 1), rep_env%f(i, j + 1) 456 END DO 457 END IF 458 END SELECT 459 END IF 460 END DO 461 END DO 462 DEALLOCATE (cvalues) 463 DEALLOCATE (Mmatrix_tmp) 464 IF (PRESENT(neb_env)) THEN 465 ! First identify the image of the chain with the higher potential energy 466 ! First and last point of the band are never considered 467 neb_env%nr_HE_image = MAXLOC(energies(2:n_rep_neb - 1), 1) + 1 468 ALLOCATE (tangent(nsize_wrk)) 469 ! Then modify image forces accordingly to the scheme chosen for the 470 ! calculation. 471 neb_env%spring_energy = 0.0_dp 472 IF (neb_env%optimize_end_points) THEN 473 ALLOCATE (tmp_a(SIZE(forces%wrk, 1))) 474 ALLOCATE (tmp_b(SIZE(forces%wrk, 1))) 475 tmp_a(:) = forces%wrk(:, 1) 476 tmp_b(:) = forces%wrk(:, SIZE(forces%wrk, 2)) 477 END IF 478 DO i = 2, neb_env%number_of_replica 479 CALL get_tangent(neb_env, coords, i, tangent, energies, output_unit) 480 CALL get_neb_force(neb_env, tangent, coords, i, forces, Mmatrix=Mmatrix, & 481 iw=output_unit) 482 END DO 483 IF (neb_env%optimize_end_points) THEN 484 forces%wrk(:, 1) = tmp_a ! Image A 485 forces%wrk(:, SIZE(forces%wrk, 2)) = tmp_b ! Image B 486 DEALLOCATE (tmp_a) 487 DEALLOCATE (tmp_b) 488 ELSE 489 ! Nullify forces on the two end points images 490 forces%wrk(:, 1) = 0.0_dp ! Image A 491 forces%wrk(:, SIZE(forces%wrk, 2)) = 0.0_dp ! Image B 492 END IF 493 DEALLOCATE (tangent) 494 END IF 495 DEALLOCATE (Mmatrix) 496 CALL timestop(handle) 497 END SUBROUTINE neb_calc_energy_forces 498 499! ************************************************************************************************** 500!> \brief Driver to perform an MD run on each single replica to 501!> compute specifically Free Energies in a NEB scheme 502!> \param rep_env ... 503!> \param coords ... 504!> \param irep ... 505!> \param n_rep_neb ... 506!> \param cvalues ... 507!> \param Mmatrix ... 508!> \author Teodoro Laino 01.2007 509! ************************************************************************************************** 510 SUBROUTINE perform_replica_md(rep_env, coords, irep, n_rep_neb, cvalues, Mmatrix) 511 TYPE(replica_env_type), POINTER :: rep_env 512 TYPE(neb_var_type), POINTER :: coords 513 INTEGER, INTENT(IN) :: irep, n_rep_neb 514 REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: cvalues, Mmatrix 515 516 CHARACTER(len=*), PARAMETER :: routineN = 'perform_replica_md', & 517 routineP = moduleN//':'//routineN 518 519 INTEGER :: handle, handle2, ierr, j, n_el 520 LOGICAL :: explicit 521 TYPE(cp_logger_type), POINTER :: logger 522 TYPE(f_env_type), POINTER :: f_env 523 TYPE(global_environment_type), POINTER :: globenv 524 TYPE(section_vals_type), POINTER :: md_section, root_section 525 526 CALL timeset(routineN, handle) 527 CALL f_env_add_defaults(f_env_id=rep_env%f_env_id, f_env=f_env, & 528 handle=handle2) 529 logger => cp_get_default_logger() 530 CALL force_env_get(f_env%force_env, globenv=globenv, & 531 root_section=root_section) 532 j = rep_env%local_rep_indices(1) - 1 533 n_el = 3*rep_env%nparticle 534 Mmatrix = 0.0_dp 535 ! Syncronize position on the replica procs 536 CALL set_pos(rep_env%f_env_id, rep_env%r(:, j + 1), n_el, ierr) 537 CPASSERT(ierr == 0) 538 ! 539 IF (irep + j <= n_rep_neb) THEN 540 logger%iter_info%iteration(2) = irep + j 541 CALL remove_restart_info(root_section) 542 md_section => section_vals_get_subs_vals(root_section, "MOTION%MD") 543 CALL section_vals_get(md_section, explicit=explicit) 544 CPASSERT(explicit) 545 ! Let's syncronize the target of Collective Variables for this run 546 CALL set_colvars_target(coords%int(:, irep + j), f_env%force_env) 547 ! Do a molecular dynamics and get back the derivative 548 ! of the free energy w.r.t. the colvar and the metric tensor 549 CALL qs_mol_dyn(f_env%force_env, globenv=globenv) 550 ! Collect the equilibrated coordinates 551 CALL get_pos(rep_env%f_env_id, rep_env%r(1:n_el, j + 1), n_el, ierr) 552 CPASSERT(ierr == 0) 553 ! Write he gradients in the colvar coordinates into the replica_env array 554 ! and copy back also the metric tensor.. 555 ! work in progress.. 556 CPABORT("") 557 rep_env%f(:, j + 1) = 0.0_dp 558 Mmatrix = 0.0_dp 559 ELSE 560 rep_env%r(:, j + 1) = 0.0_dp 561 rep_env%f(:, j + 1) = 0.0_dp 562 cvalues(:, j + 1) = 0.0_dp 563 Mmatrix(:, j + 1) = 0.0_dp 564 END IF 565 CALL rep_env_sync(rep_env, rep_env%f) 566 CALL rep_env_sync(rep_env, rep_env%r) 567 CALL rep_env_sync(rep_env, cvalues) 568 CALL rep_env_sync(rep_env, Mmatrix) 569 CALL f_env_rm_defaults(f_env=f_env, ierr=ierr, handle=handle2) 570 CPASSERT(ierr == 0) 571 CALL timestop(handle) 572 END SUBROUTINE perform_replica_md 573 574! ************************************************************************************************** 575!> \brief Driver to perform a GEO_OPT run on each single replica to 576!> compute specifically minimum energies in a collective variable 577!> NEB scheme 578!> \param rep_env ... 579!> \param coords ... 580!> \param irep ... 581!> \param n_rep_neb ... 582!> \param cvalues ... 583!> \param Mmatrix ... 584!> \author Teodoro Laino 05.2007 585! ************************************************************************************************** 586 SUBROUTINE perform_replica_geo(rep_env, coords, irep, n_rep_neb, cvalues, Mmatrix) 587 TYPE(replica_env_type), POINTER :: rep_env 588 TYPE(neb_var_type), POINTER :: coords 589 INTEGER, INTENT(IN) :: irep, n_rep_neb 590 REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: cvalues, Mmatrix 591 592 CHARACTER(len=*), PARAMETER :: routineN = 'perform_replica_geo', & 593 routineP = moduleN//':'//routineN 594 595 INTEGER :: handle, handle2, ierr, j, n_el 596 LOGICAL :: explicit 597 TYPE(cp_logger_type), POINTER :: logger 598 TYPE(f_env_type), POINTER :: f_env 599 TYPE(global_environment_type), POINTER :: globenv 600 TYPE(section_vals_type), POINTER :: geoopt_section, root_section 601 602 CALL timeset(routineN, handle) 603 CALL f_env_add_defaults(f_env_id=rep_env%f_env_id, f_env=f_env, & 604 handle=handle2) 605 logger => cp_get_default_logger() 606 CALL force_env_get(f_env%force_env, globenv=globenv, & 607 root_section=root_section) 608 j = rep_env%local_rep_indices(1) - 1 609 n_el = 3*rep_env%nparticle 610 Mmatrix = 0.0_dp 611 ! Syncronize position on the replica procs 612 CALL set_pos(rep_env%f_env_id, rep_env%r(:, j + 1), n_el, ierr) 613 CPASSERT(ierr == 0) 614 IF (irep + j <= n_rep_neb) THEN 615 logger%iter_info%iteration(2) = irep + j 616 CALL remove_restart_info(root_section) 617 geoopt_section => section_vals_get_subs_vals(root_section, "MOTION%GEO_OPT") 618 CALL section_vals_get(geoopt_section, explicit=explicit) 619 CPASSERT(explicit) 620 ! Let's syncronize the target of Collective Variables for this run 621 CALL set_colvars_target(coords%int(:, irep + j), f_env%force_env) 622 ! Do a geometry optimization.. 623 CALL cp_geo_opt(f_env%force_env, globenv=globenv) 624 ! Once the geometry optimization is ended let's do a single run 625 ! without any constraints/restraints 626 CALL force_env_calc_energy_force(f_env%force_env, & 627 calc_force=.TRUE., skip_external_control=.TRUE.) 628 ! Collect the optimized coordinates 629 CALL get_pos(rep_env%f_env_id, rep_env%r(1:n_el, j + 1), n_el, ierr) 630 CPASSERT(ierr == 0) 631 ! Collect the gradients in cartesian coordinates 632 CALL get_force(rep_env%f_env_id, rep_env%f(1:n_el, j + 1), n_el, ierr) 633 CPASSERT(ierr == 0) 634 ! Copy the energy 635 CALL get_energy(rep_env%f_env_id, rep_env%f(n_el + 1, j + 1), ierr) 636 CPASSERT(ierr == 0) 637 ! The gradients in the colvar coordinates 638 CALL get_clv_force(f_env%force_env, rep_env%f(1:n_el, j + 1), rep_env%r(1:n_el, j + 1), & 639 SIZE(coords%xyz, 1), SIZE(coords%wrk, 1), cvalues(:, j + 1), Mmatrix(:, j + 1)) 640 ELSE 641 rep_env%r(:, j + 1) = 0.0_dp 642 rep_env%f(:, j + 1) = 0.0_dp 643 cvalues(:, j + 1) = 0.0_dp 644 Mmatrix(:, j + 1) = 0.0_dp 645 END IF 646 CALL rep_env_sync(rep_env, rep_env%f) 647 CALL rep_env_sync(rep_env, rep_env%r) 648 CALL rep_env_sync(rep_env, cvalues) 649 CALL rep_env_sync(rep_env, Mmatrix) 650 CALL f_env_rm_defaults(f_env=f_env, ierr=ierr, handle=handle2) 651 CPASSERT(ierr == 0) 652 CALL timestop(handle) 653 END SUBROUTINE perform_replica_geo 654 655! ************************************************************************************************** 656!> \brief Computes the tangent for point i of the NEB chain 657!> \param neb_env ... 658!> \param coords ... 659!> \param i ... 660!> \param tangent ... 661!> \param energies ... 662!> \param iw ... 663!> \author Teodoro Laino 09.2006 664! ************************************************************************************************** 665 SUBROUTINE get_tangent(neb_env, coords, i, tangent, energies, iw) 666 TYPE(neb_type), POINTER :: neb_env 667 TYPE(neb_var_type), POINTER :: coords 668 INTEGER, INTENT(IN) :: i 669 REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: tangent 670 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: energies 671 INTEGER, INTENT(IN) :: iw 672 673 CHARACTER(len=*), PARAMETER :: routineN = 'get_tangent', routineP = moduleN//':'//routineN 674 675 REAL(KINd=dp) :: distance0, distance1, distance2, DVmax, & 676 Dvmin 677 678 CPASSERT(ASSOCIATED(coords)) 679 tangent(:) = 0.0_dp 680 ! For the last point we don't need any tangent.. 681 IF (i == neb_env%number_of_replica) RETURN 682 ! Several kind of tangents implemented... 683 SELECT CASE (neb_env%id_type) 684 CASE (do_eb) 685 tangent(:) = 0.0_dp 686 CASE (do_b_neb) 687 CALL neb_replica_distance(coords=coords, i0=i, i=i - 1, distance=distance1, iw=iw, & 688 rotate=.FALSE.) 689 CALL neb_replica_distance(coords=coords, i0=i + 1, i=i, distance=distance2, iw=iw, & 690 rotate=.FALSE.) 691 tangent(:) = (coords%wrk(:, i) - coords%wrk(:, i - 1))/distance1 + & 692 (coords%wrk(:, i + 1) - coords%wrk(:, i))/distance2 693 CASE (do_it_neb, do_ci_neb, do_d_neb) 694 IF ((energies(i + 1) .GT. energies(i)) .AND. (energies(i) .GT. (energies(i - 1)))) THEN 695 tangent(:) = coords%wrk(:, i + 1) - coords%wrk(:, i) 696 ELSE IF ((energies(i + 1) .LT. energies(i)) .AND. (energies(i) .LT. (energies(i - 1)))) THEN 697 tangent(:) = coords%wrk(:, i) - coords%wrk(:, i - 1) 698 ELSE 699 DVmax = MAX(ABS(energies(i + 1) - energies(i)), ABS(energies(i - 1) - energies(i))) 700 DVmin = MIN(ABS(energies(i + 1) - energies(i)), ABS(energies(i - 1) - energies(i))) 701 IF (energies(i + 1) .GE. energies(i - 1)) THEN 702 tangent(:) = (coords%wrk(:, i + 1) - coords%wrk(:, i))*DVmax + (coords%wrk(:, i) - coords%wrk(:, i - 1))*DVmin 703 ELSE 704 tangent(:) = (coords%wrk(:, i + 1) - coords%wrk(:, i))*DVmin + (coords%wrk(:, i) - coords%wrk(:, i - 1))*DVmax 705 END IF 706 END IF 707 CASE (do_sm) 708 ! String method.. 709 tangent(:) = 0.0_dp 710 END SELECT 711 distance0 = SQRT(DOT_PRODUCT(tangent(:), tangent(:))) 712 IF (distance0 /= 0.0_dp) tangent(:) = tangent(:)/distance0 713 END SUBROUTINE get_tangent 714 715! ************************************************************************************************** 716!> \brief Computes the forces for point i of the NEB chain 717!> \param neb_env ... 718!> \param tangent ... 719!> \param coords ... 720!> \param i ... 721!> \param forces ... 722!> \param tag ... 723!> \param Mmatrix ... 724!> \param iw ... 725!> \author Teodoro Laino 09.2006 726! ************************************************************************************************** 727 RECURSIVE SUBROUTINE get_neb_force(neb_env, tangent, coords, i, forces, tag, Mmatrix, & 728 iw) 729 TYPE(neb_type), POINTER :: neb_env 730 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: tangent 731 TYPE(neb_var_type), POINTER :: coords 732 INTEGER, INTENT(IN) :: i 733 TYPE(neb_var_type), POINTER :: forces 734 INTEGER, INTENT(IN), OPTIONAL :: tag 735 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: Mmatrix 736 INTEGER, INTENT(IN) :: iw 737 738 CHARACTER(len=*), PARAMETER :: routineN = 'get_neb_force', routineP = moduleN//':'//routineN 739 740 INTEGER :: j, my_tag, nsize_wrk 741 REAL(KIND=dp) :: distance1, distance2, fac, tmp 742 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: dtmp1, wrk 743 744 my_tag = neb_env%id_type 745 IF (PRESENT(tag)) my_tag = tag 746 CPASSERT(ASSOCIATED(forces)) 747 CPASSERT(ASSOCIATED(coords)) 748 nsize_wrk = coords%size_wrk(1) 749 ! All methods but not the classical elastic band will skip the force 750 ! calculation for the last frame of the band 751 SELECT CASE (my_tag) 752 CASE (do_b_neb, do_it_neb, do_ci_neb, do_d_neb) 753 IF (i == neb_env%number_of_replica) RETURN 754 CASE (do_sm) 755 ! String Method 756 ! The forces do not require any projection. Reparametrization required 757 ! after the update of the replica. 758 CALL cite_reference(E2002) 759 RETURN 760 END SELECT 761 ! otherwise proceeed normally.. 762 ALLOCATE (wrk(nsize_wrk)) 763 ! Spring Energy 764 CALL neb_replica_distance(coords=coords, i0=i - 1, i=i, distance=distance1, iw=iw, & 765 rotate=.FALSE.) 766 tmp = distance1 - neb_env%avg_distance 767 neb_env%spring_energy = neb_env%spring_energy + 0.5_dp*neb_env%k*tmp**2 768 SELECT CASE (my_tag) 769 CASE (do_eb) 770 CALL cite_reference(Elber1987) 771 ! Elastic band - Hamiltonian formulation according the original Karplus/Elber 772 ! formulation 773 ALLOCATE (dtmp1(nsize_wrk)) 774 ! derivatives of the spring 775 tmp = distance1 - neb_env%avg_distance 776 dtmp1(:) = 1.0_dp/distance1*(coords%wrk(:, i) - coords%wrk(:, i - 1)) 777 wrk(:) = neb_env%k*tmp*dtmp1 778 forces%wrk(:, i) = forces%wrk(:, i) - wrk 779 forces%wrk(:, i - 1) = forces%wrk(:, i - 1) + wrk 780 ! derivatives due to the average length of the spring 781 fac = 1.0_dp/(neb_env%avg_distance*REAL(neb_env%number_of_replica - 1, KIND=dp)) 782 wrk(:) = neb_env%k*fac*(coords%wrk(:, i) - coords%wrk(:, i - 1)) 783 tmp = 0.0_dp 784 DO j = 2, neb_env%number_of_replica 785 CALL neb_replica_distance(coords=coords, i0=j - 1, i=j, distance=distance1, iw=iw, & 786 rotate=.FALSE.) 787 tmp = tmp + distance1 - neb_env%avg_distance 788 END DO 789 forces%wrk(:, i) = forces%wrk(:, i) + wrk*tmp 790 forces%wrk(:, i - 1) = forces%wrk(:, i - 1) - wrk*tmp 791 DEALLOCATE (dtmp1) 792 CASE (do_b_neb) 793 ! Bisection NEB 794 CALL cite_reference(Jonsson1998) 795 wrk(:) = (coords%wrk(:, i + 1) - 2.0_dp*coords%wrk(:, i) + coords%wrk(:, i - 1)) 796 tmp = neb_env%k*DOT_PRODUCT(wrk, tangent) 797 wrk(:) = forces%wrk(:, i) - dot_product_band(neb_env, forces%wrk(:, i), tangent, Mmatrix)*tangent 798 forces%wrk(:, i) = wrk + tmp*tangent 799 CASE (do_it_neb) 800 ! Improved tangent NEB 801 CALL cite_reference(Jonsson2000_1) 802 CALL neb_replica_distance(coords=coords, i0=i, i=i + 1, distance=distance1, iw=iw, & 803 rotate=.FALSE.) 804 CALL neb_replica_distance(coords=coords, i0=i - 1, i=i, distance=distance2, iw=iw, & 805 rotate=.FALSE.) 806 tmp = neb_env%k*(distance1 - distance2) 807 wrk(:) = forces%wrk(:, i) - dot_product_band(neb_env, forces%wrk(:, i), tangent, Mmatrix)*tangent 808 forces%wrk(:, i) = wrk + tmp*tangent 809 CASE (do_ci_neb) 810 ! Climbing Image NEB 811 CALL cite_reference(Jonsson2000_2) 812 IF (neb_env%istep <= neb_env%nsteps_it .OR. i /= neb_env%nr_HE_image) THEN 813 CALL get_neb_force(neb_env, tangent, coords, i, forces, do_it_neb, Mmatrix, iw) 814 ELSE 815 wrk(:) = forces%wrk(:, i) 816 tmp = -2.0_dp*dot_product_band(neb_env, wrk, tangent, Mmatrix) 817 forces%wrk(:, i) = wrk + tmp*tangent 818 END IF 819 CASE (do_d_neb) 820 ! Doubly NEB 821 CALL cite_reference(Wales2004) 822 ALLOCATE (dtmp1(nsize_wrk)) 823 dtmp1(:) = forces%wrk(:, i) - dot_product_band(neb_env, forces%wrk(:, i), tangent, Mmatrix)*tangent 824 forces%wrk(:, i) = dtmp1 825 tmp = SQRT(DOT_PRODUCT(dtmp1, dtmp1)) 826 dtmp1(:) = dtmp1(:)/tmp 827 ! Project out only the spring component interfering with the 828 ! orthogonal gradient of the band 829 wrk(:) = (coords%wrk(:, i + 1) - 2.0_dp*coords%wrk(:, i) + coords%wrk(:, i - 1)) 830 tmp = DOT_PRODUCT(wrk, dtmp1) 831 dtmp1(:) = neb_env%k*(wrk(:) - tmp*dtmp1(:)) 832 forces%wrk(:, i) = forces%wrk(:, i) + dtmp1(:) 833 DEALLOCATE (dtmp1) 834 END SELECT 835 DEALLOCATE (wrk) 836 END SUBROUTINE get_neb_force 837 838! ************************************************************************************************** 839!> \brief Handles the dot_product when using colvar.. in this case 840!> the scalar product needs to take into account the metric 841!> tensor 842!> \param neb_env ... 843!> \param array1 ... 844!> \param array2 ... 845!> \param array3 ... 846!> \return ... 847!> \author Teodoro Laino 09.2006 848! ************************************************************************************************** 849 FUNCTION dot_product_band(neb_env, array1, array2, array3) RESULT(value) 850 TYPE(neb_type), POINTER :: neb_env 851 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: array1, array2 852 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: array3 853 REAL(KIND=dp) :: value 854 855 CHARACTER(len=*), PARAMETER :: routineN = 'dot_product_band', & 856 routineP = moduleN//':'//routineN 857 858 INTEGER :: nsize_int 859 LOGICAL :: check 860 861 IF (neb_env%use_colvar) THEN 862 nsize_int = neb_env%nsize_int 863 check = ((SIZE(array1) /= SIZE(array2)) .OR. & 864 (SIZE(array1) /= nsize_int) .OR. & 865 (SIZE(array3) /= nsize_int*nsize_int)) 866 ! This condition should always be satisfied.. 867 CPASSERT(check) 868 value = DOT_PRODUCT(MATMUL(RESHAPE(array3, (/nsize_int, nsize_int/)), array1), array2) 869 ELSE 870 value = DOT_PRODUCT(array1, array2) 871 END IF 872 END FUNCTION dot_product_band 873 874! ************************************************************************************************** 875!> \brief Reorient iteratively all images of the NEB chain in order to 876!> have always the smaller RMSD between two following images 877!> \param rotate_frames ... 878!> \param particle_set ... 879!> \param coords ... 880!> \param vels ... 881!> \param iw ... 882!> \param distances ... 883!> \param number_of_replica ... 884!> \author Teodoro Laino 09.2006 885! ************************************************************************************************** 886 SUBROUTINE reorient_images(rotate_frames, particle_set, coords, vels, iw, & 887 distances, number_of_replica) 888 LOGICAL, INTENT(IN) :: rotate_frames 889 TYPE(particle_type), DIMENSION(:), OPTIONAL, & 890 POINTER :: particle_set 891 TYPE(neb_var_type), POINTER :: coords, vels 892 INTEGER, INTENT(IN) :: iw 893 REAL(KIND=dp), DIMENSION(:), OPTIONAL :: distances 894 INTEGER, INTENT(IN) :: number_of_replica 895 896 CHARACTER(len=*), PARAMETER :: routineN = 'reorient_images', & 897 routineP = moduleN//':'//routineN 898 899 INTEGER :: i, k, kind 900 LOGICAL :: check 901 REAL(KIND=dp) :: xtmp 902 REAL(KIND=dp), DIMENSION(3) :: tmp 903 REAL(KIND=dp), DIMENSION(3, 3) :: rot 904 905 rot = 0.0_dp 906 rot(1, 1) = 1.0_dp 907 rot(2, 2) = 1.0_dp 908 rot(3, 3) = 1.0_dp 909 DO i = 2, number_of_replica 910 ! The rotation of the replica is enabled exclusively when working in 911 ! cartesian coordinates 912 IF (rotate_frames .AND. (coords%in_use == do_band_cartesian)) THEN 913 CALL rmsd3(particle_set, coords%xyz(:, i), coords%xyz(:, i - 1), iw, & 914 rotate=.TRUE., rot=rot) 915 ! Rotate velocities 916 DO k = 1, SIZE(vels%xyz, 1)/3 917 kind = (k - 1)*3 918 tmp = vels%xyz(kind + 1:kind + 3, i) 919 CALL matvec_3x3(vels%xyz(kind + 1:kind + 3, i), TRANSPOSE(rot), tmp) 920 END DO 921 END IF 922 IF (PRESENT(distances)) THEN 923 check = SIZE(distances) == (number_of_replica - 1) 924 CPASSERT(check) 925 xtmp = DOT_PRODUCT(coords%wrk(:, i) - coords%wrk(:, i - 1), & 926 coords%wrk(:, i) - coords%wrk(:, i - 1)) 927 distances(i - 1) = SQRT(xtmp) 928 END IF 929 END DO 930 END SUBROUTINE reorient_images 931 932! ************************************************************************************************** 933!> \brief Reparametrization of the replica for String Method with splines 934!> \param reparametrize_frames ... 935!> \param spline_order ... 936!> \param smoothing ... 937!> \param coords ... 938!> \param sline ... 939!> \param distances ... 940!> \author Teodoro Laino - Rodolphe Vuilleumier 09.2008 941! ************************************************************************************************** 942 SUBROUTINE reparametrize_images(reparametrize_frames, spline_order, smoothing, & 943 coords, sline, distances) 944 945 LOGICAL, INTENT(IN) :: reparametrize_frames 946 INTEGER, INTENT(IN) :: spline_order 947 REAL(KIND=dp), INTENT(IN) :: smoothing 948 REAL(KIND=dp), DIMENSION(:, :), POINTER :: coords, sline 949 REAL(KIND=dp), DIMENSION(:) :: distances 950 951 CHARACTER(len=*), PARAMETER :: routineN = 'reparametrize_images', & 952 routineP = moduleN//':'//routineN 953 954 INTEGER :: i, j 955 REAL(KIND=dp) :: avg_distance, xtmp 956 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: tmp_coords 957 958 IF (reparametrize_frames) THEN 959 ALLOCATE (tmp_coords(SIZE(coords, 1), SIZE(coords, 2))) 960 tmp_coords(:, :) = coords 961 ! Smoothing 962 DO i = 2, SIZE(coords, 2) - 1 963 coords(:, i) = tmp_coords(:, i)*(1.0_dp - 2.0_dp*smoothing) + & 964 tmp_coords(:, i - 1)*smoothing + tmp_coords(:, i + 1)*smoothing 965 END DO 966 sline = coords - tmp_coords + sline 967 tmp_coords(:, :) = coords 968 ! Reparametrization 969 SELECT CASE (spline_order) 970 CASE (1) 971 ! Compute distances 972 DO i = 2, SIZE(coords, 2) 973 xtmp = DOT_PRODUCT(coords(:, i) - coords(:, i - 1), coords(:, i) - coords(:, i - 1)) 974 distances(i - 1) = SQRT(xtmp) 975 END DO 976 avg_distance = SUM(distances)/REAL(SIZE(coords, 2) - 1, KIND=dp) 977 ! Redistribute frames 978 DO i = 2, SIZE(coords, 2) - 1 979 xtmp = 0.0_dp 980 DO j = 1, SIZE(coords, 2) - 1 981 xtmp = xtmp + distances(j) 982 IF (xtmp > avg_distance*REAL(i - 1, KIND=dp)) THEN 983 xtmp = (xtmp - avg_distance*REAL(i - 1, KIND=dp))/distances(j) 984 coords(:, i) = (1.0_dp - xtmp)*tmp_coords(:, j + 1) + xtmp*tmp_coords(:, j) 985 EXIT 986 END IF 987 END DO 988 END DO 989 ! Re-compute distances 990 DO i = 2, SIZE(coords, 2) 991 xtmp = DOT_PRODUCT(coords(:, i) - coords(:, i - 1), coords(:, i) - coords(:, i - 1)) 992 distances(i - 1) = SQRT(xtmp) 993 END DO 994 CASE DEFAULT 995 CPWARN("String Method: Spline order greater than 1 not implemented.") 996 END SELECT 997 sline = coords - tmp_coords + sline 998 DEALLOCATE (tmp_coords) 999 END IF 1000 END SUBROUTINE reparametrize_images 1001 1002! ************************************************************************************************** 1003!> \brief Checks for convergence criteria during a NEB run 1004!> \param neb_env ... 1005!> \param Dcoords ... 1006!> \param forces ... 1007!> \return ... 1008!> \author Teodoro Laino 10.2006 1009! ************************************************************************************************** 1010 FUNCTION check_convergence(neb_env, Dcoords, forces) RESULT(converged) 1011 TYPE(neb_type), POINTER :: neb_env 1012 TYPE(neb_var_type), POINTER :: Dcoords, forces 1013 LOGICAL :: converged 1014 1015 CHARACTER(len=*), PARAMETER :: routineN = 'check_convergence', & 1016 routineP = moduleN//':'//routineN 1017 1018 CHARACTER(LEN=3), DIMENSION(4) :: labels 1019 INTEGER :: iw 1020 REAL(KIND=dp) :: max_dr, max_force, my_max_dr, & 1021 my_max_force, my_rms_dr, my_rms_force, & 1022 rms_dr, rms_force 1023 TYPE(cp_logger_type), POINTER :: logger 1024 TYPE(section_vals_type), POINTER :: cc_section 1025 1026 NULLIFY (logger, cc_section) 1027 logger => cp_get_default_logger() 1028 cc_section => section_vals_get_subs_vals(neb_env%neb_section, "CONVERGENCE_CONTROL") 1029 CALL section_vals_val_get(cc_section, "MAX_DR", r_val=max_dr) 1030 CALL section_vals_val_get(cc_section, "MAX_FORCE", r_val=max_force) 1031 CALL section_vals_val_get(cc_section, "RMS_DR", r_val=rms_dr) 1032 CALL section_vals_val_get(cc_section, "RMS_FORCE", r_val=rms_force) 1033 converged = .FALSE. 1034 labels = " NO" 1035 my_max_dr = MAXVAL(ABS(Dcoords%wrk)) 1036 my_max_force = MAXVAL(ABS(forces%wrk)) 1037 my_rms_dr = SQRT(SUM(Dcoords%wrk*Dcoords%wrk)/REAL(SIZE(Dcoords%wrk, 1)*SIZE(Dcoords%wrk, 2), KIND=dp)) 1038 my_rms_force = SQRT(SUM(forces%wrk*forces%wrk)/REAL(SIZE(forces%wrk, 1)*SIZE(forces%wrk, 2), KIND=dp)) 1039 IF (my_max_dr < max_dr) labels(1) = "YES" 1040 IF (my_max_force < max_force) labels(2) = "YES" 1041 IF (my_rms_dr < rms_dr) labels(3) = "YES" 1042 IF (my_rms_force < rms_force) labels(4) = "YES" 1043 IF (ALL(labels == "YES")) converged = .TRUE. 1044 1045 iw = cp_print_key_unit_nr(logger, neb_env%neb_section, "CONVERGENCE_INFO", & 1046 extension=".nebLog") 1047 IF (iw > 0) THEN 1048 ! Print convergence info 1049 WRITE (iw, FMT='(A,A)') ' **************************************', & 1050 '*****************************************' 1051 WRITE (iw, FMT='(1X,A,2X,F8.5,5X,"[",F8.5,"]",1X,T76,"(",A,")")') & 1052 'RMS DISPLACEMENT =', my_rms_dr, rms_dr, labels(3), & 1053 'MAX DISPLACEMENT =', my_max_dr, max_dr, labels(1), & 1054 'RMS FORCE =', my_rms_force, rms_force, labels(4), & 1055 'MAX FORCE =', my_max_force, max_force, labels(2) 1056 WRITE (iw, FMT='(A,A)') ' **************************************', & 1057 '*****************************************' 1058 END IF 1059 CALL cp_print_key_finished_output(iw, logger, neb_env%neb_section, & 1060 "CONVERGENCE_INFO") 1061 END FUNCTION check_convergence 1062 1063END MODULE neb_utils 1064