75 #include "../base/base_uses.f90"
79 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'neb_utils'
80 LOGICAL,
PARAMETER,
PRIVATE :: debug_this_module = .false.
101 SUBROUTINE neb_replica_distance(particle_set, coords, i0, i, distance, iw, rotate)
102 TYPE(particle_type),
DIMENSION(:),
OPTIONAL, &
103 POINTER :: particle_set
104 TYPE(neb_var_type),
POINTER :: coords
105 INTEGER,
INTENT(IN) :: i0, i
106 REAL(KIND=
dp),
INTENT(OUT) :: distance
107 INTEGER,
INTENT(IN) :: iw
108 LOGICAL,
INTENT(IN),
OPTIONAL :: rotate
113 IF (
PRESENT(rotate)) my_rotate = rotate
117 cpassert(
PRESENT(particle_set))
118 CALL rmsd3(particle_set, coords%xyz(:, i), coords%xyz(:, i0), &
119 iw, rotate=my_rotate)
121 distance = sqrt(dot_product(coords%wrk(:, i) - coords%wrk(:, i0), &
122 coords%wrk(:, i) - coords%wrk(:, i0)))
124 END SUBROUTINE neb_replica_distance
139 coords, vels, neb_env, iw, globenv, para_env)
140 TYPE(section_vals_type),
POINTER :: neb_section
141 TYPE(particle_type),
DIMENSION(:),
POINTER :: particle_set
142 TYPE(neb_var_type),
POINTER :: coords, vels
143 TYPE(neb_type),
POINTER :: neb_env
144 INTEGER,
INTENT(IN) :: iw
145 TYPE(global_environment_type),
POINTER :: globenv
146 TYPE(mp_para_env_type),
POINTER :: para_env
148 CHARACTER(len=*),
PARAMETER :: routinen =
'build_replica_coords'
150 CHARACTER(LEN=default_path_length) :: filename
151 INTEGER :: handle, i_rep, iatom, ic, input_nr_replica, is, ivar, j, jtarg, n_rep, natom, &
152 neb_nr_replica, nr_replica_to_interpolate, nval, nvar, shell_index
153 LOGICAL :: check, explicit, skip_vel_section
154 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: distance
155 REAL(kind=
dp),
DIMENSION(3) :: r
156 REAL(kind=
dp),
DIMENSION(:),
POINTER :: initial_colvars, rptr
157 TYPE(section_vals_type),
POINTER :: coord_section, replica_section, &
160 CALL timeset(routinen, handle)
161 cpassert(
ASSOCIATED(coords))
162 cpassert(
ASSOCIATED(vels))
163 neb_nr_replica = neb_env%number_of_replica
167 cpassert(input_nr_replica <= neb_nr_replica)
169 skip_vel_section = (input_nr_replica /= neb_nr_replica)
170 IF ((iw > 0) .AND. skip_vel_section)
THEN
171 WRITE (iw,
'(T2,A)')
'NEB| The number of replica in the input is different from the number', &
172 'NEB| of replica requested for NEB. More Replica will be interpolated.', &
173 'NEB| Therefore the possibly provided velocities will not be read.'
176 DO i_rep = 1, input_nr_replica
180 skip_vel_section = skip_vel_section .OR. (.NOT. explicit)
183 coords%xyz(:, :) = 0.0_dp
184 DO i_rep = 1, input_nr_replica
192 cpassert((natom ==
SIZE(particle_set)))
195 i_rep_val=iatom, r_vals=rptr)
197 coords%xyz(ic + 1:ic + 3, i_rep) = rptr(1:3)*
bohr
199 shell_index = particle_set(iatom)%shell_index
200 IF (shell_index /= 0)
THEN
201 is = 3*(natom + shell_index - 1)
202 coords%xyz(is + 1:is + 3, i_rep) = coords%xyz(ic + 1:ic + 3, i_rep)
208 CHARACTER(LEN=default_string_length) :: dummy_char
209 TYPE(cp_parser_type) :: parser
211 i_rep_section=i_rep, c_val=filename)
212 cpassert(trim(filename) /=
"")
213 CALL parser_create(parser, filename, para_env=para_env, parse_white_lines=.true.)
216 CALL parser_get_object(parser, natom)
217 cpassert((natom ==
SIZE(particle_set)))
223 CALL cp_abort(__location__, &
224 "Number of lines in XYZ format not equal to the number of atoms."// &
225 " Error in XYZ format for REPLICA coordinates. Very probably the"// &
226 " line with title is missing or is empty. Please check the XYZ file and rerun your job!")
227 READ (parser%input_line, *) dummy_char, r(1:3)
229 coords%xyz(ic + 1:ic + 3, i_rep) = r(1:3)*
bohr
231 shell_index = particle_set(iatom)%shell_index
232 IF (shell_index /= 0)
THEN
233 is = 3*(natom + shell_index - 1)
234 coords%xyz(is + 1:is + 3, i_rep) = coords%xyz(ic + 1:ic + 3, i_rep)
241 IF (neb_env%use_colvar)
THEN
243 i_rep_section=i_rep, n_rep_val=n_rep)
246 NULLIFY (initial_colvars)
248 i_rep_section=i_rep, r_vals=initial_colvars)
249 check = (neb_env%nsize_int ==
SIZE(initial_colvars))
251 coords%int(:, i_rep) = initial_colvars
254 CALL eval_colvar(neb_env%force_env, coords%xyz(:, i_rep), coords%int(:, i_rep))
260 IF (skip_vel_section)
THEN
262 i_rep, iw, globenv, neb_env)
269 IF (neb_env%use_colvar)
THEN
270 nvar =
SIZE(vels%wrk, 1)
271 cpassert(nval == nvar)
274 i_rep_val=ivar, r_vals=rptr)
275 vels%wrk(ivar, i_rep) = rptr(1)
278 natom =
SIZE(particle_set)
279 cpassert(nval == natom)
282 i_rep_val=iatom, r_vals=rptr)
284 vels%wrk(ic + 1:ic + 3, i_rep) = rptr(1:3)
286 shell_index = particle_set(iatom)%shell_index
287 IF (shell_index /= 0)
THEN
288 is = 3*(natom + shell_index - 1)
289 vels%wrk(is + 1:is + 3, i_rep) = vels%wrk(ic + 1:ic + 3, i_rep)
295 ALLOCATE (distance(neb_nr_replica - 1))
296 IF (input_nr_replica < neb_nr_replica)
THEN
298 nr_replica_to_interpolate = neb_nr_replica - input_nr_replica
301 WRITE (iw,
'(T2,A,I0,A)')
'NEB| Interpolating ', nr_replica_to_interpolate,
' missing Replica.'
303 DO WHILE (nr_replica_to_interpolate > 0)
306 DO j = 1, input_nr_replica - 1
307 CALL neb_replica_distance(particle_set, coords, j, j + 1, distance(j), iw, &
308 rotate=neb_env%align_frames)
310 jtarg = maxloc(distance(1:input_nr_replica), 1)
312 WRITE (iw,
'(T2,3(A,I0),A)')
'NEB| Interpolating Nr. ', &
313 nr_replica_to_interpolate,
' missing Replica between Replica Nr. ', &
314 jtarg,
' and ', jtarg + 1,
'.'
316 input_nr_replica = input_nr_replica + 1
317 nr_replica_to_interpolate = nr_replica_to_interpolate - 1
319 coords%xyz(:, jtarg + 2:input_nr_replica) = coords%xyz(:, jtarg + 1:input_nr_replica - 1)
320 coords%xyz(:, jtarg + 1) = (coords%xyz(:, jtarg) + coords%xyz(:, jtarg + 2))/2.0_dp
321 IF (neb_env%use_colvar)
THEN
326 coords%int(:, jtarg + 2:input_nr_replica) = coords%int(:, jtarg + 1:input_nr_replica - 1)
327 coords%int(:, jtarg + 1) = (coords%int(:, jtarg) + coords%int(:, jtarg + 2))/2.0_dp
329 vels%wrk(:, jtarg + 2:input_nr_replica) = vels%wrk(:, jtarg + 1:input_nr_replica - 1)
330 vels%wrk(:, jtarg + 1) = 0.0_dp
332 input_nr_replica, iw, neb_env%use_colvar)
334 jtarg + 1, iw, globenv, neb_env)
337 vels%wrk(:, 1) = 0.0_dp
338 vels%wrk(:, neb_nr_replica) = 0.0_dp
343 DO j = 1, input_nr_replica - 1
344 CALL neb_replica_distance(particle_set, coords, j, j + 1, distance(j), iw, &
345 rotate=neb_env%align_frames)
347 DEALLOCATE (distance)
349 CALL timestop(handle)
366 particle_set, output_unit)
367 TYPE(replica_env_type),
POINTER :: rep_env
368 TYPE(neb_type),
OPTIONAL,
POINTER :: neb_env
369 TYPE(neb_var_type),
POINTER :: coords
370 REAL(kind=dp),
DIMENSION(:),
INTENT(OUT) :: energies
371 TYPE(neb_var_type),
POINTER :: forces
372 TYPE(particle_type),
DIMENSION(:),
POINTER :: particle_set
373 INTEGER,
INTENT(IN) :: output_unit
375 CHARACTER(len=*),
PARAMETER :: routinen =
'neb_calc_energy_forces'
376 CHARACTER(LEN=1),
DIMENSION(3),
PARAMETER :: lab = (/
"X",
"Y",
"Z"/)
378 INTEGER :: handle, i, irep, j, n_int, n_rep, &
380 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:) :: tangent, tmp_a, tmp_b
381 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:, :) :: cvalues, mmatrix, mmatrix_tmp
383 CALL timeset(routinen, handle)
384 n_int = neb_env%nsize_int
385 n_rep_neb = neb_env%number_of_replica
387 nsize_wrk = coords%size_wrk(1)
389 ALLOCATE (cvalues(n_int, n_rep))
390 ALLOCATE (mmatrix_tmp(n_int*n_int, n_rep))
391 ALLOCATE (mmatrix(n_int*n_int, n_rep_neb))
392 IF (output_unit > 0)
WRITE (output_unit,
'(/,T2,A)')
"NEB| Computing Energies and Forces"
393 DO irep = 1, n_rep_neb, n_rep
395 IF (irep + j <= n_rep_neb)
THEN
400 rep_env%r(:, j + 1) = coords%xyz(:, irep + j)
405 CALL handle_band_file_names(rep_env, irep, n_rep_neb, neb_env%istep)
407 SELECT CASE (neb_env%pot_type)
410 CALL rep_env_calc_e_f(rep_env, calc_f=.true.)
413 CALL perform_replica_md(rep_env, coords, irep, n_rep_neb, cvalues, mmatrix_tmp)
416 CALL perform_replica_geo(rep_env, coords, irep, n_rep_neb, cvalues, mmatrix_tmp)
420 IF (irep + j <= n_rep_neb)
THEN
422 forces%wrk(:, irep + j) = rep_env%f(1:nsize_wrk, j + 1)
423 energies(irep + j) = rep_env%f(rep_env%ndim + 1, j + 1)
424 SELECT CASE (neb_env%pot_type)
427 IF (output_unit > 0)
THEN
428 WRITE (output_unit,
'(T2,A,I5,A,I5,A)') &
429 "NEB| REPLICA Nr.", irep + j,
"- Energy and Forces"
430 WRITE (output_unit,
'(T2,A,T43,A,T57,F24.12)') &
431 "NEB|",
"Total energy:", rep_env%f(rep_env%ndim + 1, j + 1)
432 WRITE (output_unit,
'(T2,"NEB|",T10,"ATOM",T33,3(9X,A,7X))') lab(1), lab(2), lab(3)
433 DO i = 1,
SIZE(particle_set)
434 WRITE (output_unit,
'(T2,"NEB|",T12,A,T30,3(2X,F15.9))') &
435 particle_set(i)%atomic_kind%name, &
436 rep_env%f((i - 1)*3 + 1:(i - 1)*3 + 3, j + 1)
439 CASE (pot_neb_fe, pot_neb_me)
442 coords%xyz(:, irep + j) = rep_env%r(1:rep_env%ndim, j + 1)
443 mmatrix(:, irep + j) = mmatrix_tmp(:, j + 1)
444 IF (output_unit > 0)
THEN
445 WRITE (output_unit,
'(/,T2,A,I5,A,I5,A)') &
446 "NEB| REPLICA Nr.", irep + j,
"- Energy, Collective Variables, Forces"
447 WRITE (output_unit,
'(T2,A,T43,A,T57,F24.12)') &
448 "NEB|",
"Total energy:", rep_env%f(rep_env%ndim + 1, j + 1)
449 WRITE (output_unit, &
450 '(T2,"NEB|",T10,"CV Nr.",12X,"Expected COLVAR",5X,"Present COLVAR",10X,"Forces")')
452 WRITE (output_unit,
'(T2,"NEB|",T12,I2,7X,3(5X,F15.9))') &
453 i, coords%int(i, irep + j), cvalues(i, j + 1), rep_env%f(i, j + 1)
461 DEALLOCATE (mmatrix_tmp)
462 IF (
PRESENT(neb_env))
THEN
465 neb_env%nr_HE_image = maxloc(energies(2:n_rep_neb - 1), 1) + 1
466 ALLOCATE (tangent(nsize_wrk))
469 neb_env%spring_energy = 0.0_dp
470 IF (neb_env%optimize_end_points)
THEN
471 ALLOCATE (tmp_a(
SIZE(forces%wrk, 1)))
472 ALLOCATE (tmp_b(
SIZE(forces%wrk, 1)))
473 tmp_a(:) = forces%wrk(:, 1)
474 tmp_b(:) = forces%wrk(:,
SIZE(forces%wrk, 2))
476 DO i = 2, neb_env%number_of_replica
477 CALL get_tangent(neb_env, coords, i, tangent, energies, output_unit)
478 CALL get_neb_force(neb_env, tangent, coords, i, forces, mmatrix=mmatrix, &
481 IF (neb_env%optimize_end_points)
THEN
482 forces%wrk(:, 1) = tmp_a
483 forces%wrk(:,
SIZE(forces%wrk, 2)) = tmp_b
488 forces%wrk(:, 1) = 0.0_dp
489 forces%wrk(:,
SIZE(forces%wrk, 2)) = 0.0_dp
494 CALL timestop(handle)
508 SUBROUTINE perform_replica_md(rep_env, coords, irep, n_rep_neb, cvalues, Mmatrix)
509 TYPE(replica_env_type),
POINTER :: rep_env
510 TYPE(neb_var_type),
POINTER :: coords
511 INTEGER,
INTENT(IN) :: irep, n_rep_neb
512 REAL(kind=dp),
DIMENSION(:, :),
INTENT(OUT) :: cvalues, mmatrix
514 CHARACTER(len=*),
PARAMETER :: routinen =
'perform_replica_md'
516 INTEGER :: handle, handle2, ierr, j, n_el
518 TYPE(cp_logger_type),
POINTER :: logger
519 TYPE(f_env_type),
POINTER :: f_env
520 TYPE(global_environment_type),
POINTER :: globenv
521 TYPE(section_vals_type),
POINTER :: md_section, root_section
523 CALL timeset(routinen, handle)
524 CALL f_env_add_defaults(f_env_id=rep_env%f_env_id, f_env=f_env, &
526 logger => cp_get_default_logger()
527 CALL force_env_get(f_env%force_env, globenv=globenv, &
528 root_section=root_section)
529 j = rep_env%local_rep_indices(1) - 1
530 n_el = 3*rep_env%nparticle
533 CALL set_pos(rep_env%f_env_id, rep_env%r(:, j + 1), n_el, ierr)
536 IF (irep + j <= n_rep_neb)
THEN
537 logger%iter_info%iteration(2) = irep + j
538 CALL remove_restart_info(root_section)
539 md_section => section_vals_get_subs_vals(root_section,
"MOTION%MD")
540 CALL section_vals_get(md_section, explicit=explicit)
543 CALL set_colvars_target(coords%int(:, irep + j), f_env%force_env)
546 CALL qs_mol_dyn(f_env%force_env, globenv=globenv)
548 CALL get_pos(rep_env%f_env_id, rep_env%r(1:n_el, j + 1), n_el, ierr)
554 rep_env%f(:, j + 1) = 0.0_dp
557 rep_env%r(:, j + 1) = 0.0_dp
558 rep_env%f(:, j + 1) = 0.0_dp
559 cvalues(:, j + 1) = 0.0_dp
560 mmatrix(:, j + 1) = 0.0_dp
562 CALL rep_env_sync(rep_env, rep_env%f)
563 CALL rep_env_sync(rep_env, rep_env%r)
564 CALL rep_env_sync(rep_env, cvalues)
565 CALL rep_env_sync(rep_env, mmatrix)
566 CALL f_env_rm_defaults(f_env=f_env, ierr=ierr, handle=handle2)
568 CALL timestop(handle)
569 END SUBROUTINE perform_replica_md
583 SUBROUTINE perform_replica_geo(rep_env, coords, irep, n_rep_neb, cvalues, Mmatrix)
584 TYPE(replica_env_type),
POINTER :: rep_env
585 TYPE(neb_var_type),
POINTER :: coords
586 INTEGER,
INTENT(IN) :: irep, n_rep_neb
587 REAL(kind=dp),
DIMENSION(:, :),
INTENT(OUT) :: cvalues, mmatrix
589 CHARACTER(len=*),
PARAMETER :: routinen =
'perform_replica_geo'
591 INTEGER :: handle, handle2, ierr, j, n_el
593 TYPE(cp_logger_type),
POINTER :: logger
594 TYPE(f_env_type),
POINTER :: f_env
595 TYPE(global_environment_type),
POINTER :: globenv
596 TYPE(section_vals_type),
POINTER :: geoopt_section, root_section
598 CALL timeset(routinen, handle)
599 CALL f_env_add_defaults(f_env_id=rep_env%f_env_id, f_env=f_env, &
601 logger => cp_get_default_logger()
602 CALL force_env_get(f_env%force_env, globenv=globenv, &
603 root_section=root_section)
604 j = rep_env%local_rep_indices(1) - 1
605 n_el = 3*rep_env%nparticle
608 CALL set_pos(rep_env%f_env_id, rep_env%r(:, j + 1), n_el, ierr)
610 IF (irep + j <= n_rep_neb)
THEN
611 logger%iter_info%iteration(2) = irep + j
612 CALL remove_restart_info(root_section)
613 geoopt_section => section_vals_get_subs_vals(root_section,
"MOTION%GEO_OPT")
614 CALL section_vals_get(geoopt_section, explicit=explicit)
617 CALL set_colvars_target(coords%int(:, irep + j), f_env%force_env)
619 CALL cp_geo_opt(f_env%force_env, globenv=globenv)
622 CALL force_env_calc_energy_force(f_env%force_env, &
623 calc_force=.true., skip_external_control=.true.)
625 CALL get_pos(rep_env%f_env_id, rep_env%r(1:n_el, j + 1), n_el, ierr)
628 CALL get_force(rep_env%f_env_id, rep_env%f(1:n_el, j + 1), n_el, ierr)
631 CALL get_energy(rep_env%f_env_id, rep_env%f(n_el + 1, j + 1), ierr)
634 CALL get_clv_force(f_env%force_env, rep_env%f(1:n_el, j + 1), rep_env%r(1:n_el, j + 1), &
635 SIZE(coords%xyz, 1),
SIZE(coords%wrk, 1), cvalues(:, j + 1), mmatrix(:, j + 1))
637 rep_env%r(:, j + 1) = 0.0_dp
638 rep_env%f(:, j + 1) = 0.0_dp
639 cvalues(:, j + 1) = 0.0_dp
640 mmatrix(:, j + 1) = 0.0_dp
642 CALL rep_env_sync(rep_env, rep_env%f)
643 CALL rep_env_sync(rep_env, rep_env%r)
644 CALL rep_env_sync(rep_env, cvalues)
645 CALL rep_env_sync(rep_env, mmatrix)
646 CALL f_env_rm_defaults(f_env=f_env, ierr=ierr, handle=handle2)
648 CALL timestop(handle)
649 END SUBROUTINE perform_replica_geo
661 SUBROUTINE get_tangent(neb_env, coords, i, tangent, energies, iw)
662 TYPE(neb_type),
POINTER :: neb_env
663 TYPE(neb_var_type),
POINTER :: coords
664 INTEGER,
INTENT(IN) :: i
665 REAL(kind=dp),
DIMENSION(:),
INTENT(OUT) :: tangent
666 REAL(kind=dp),
DIMENSION(:),
INTENT(IN) :: energies
667 INTEGER,
INTENT(IN) :: iw
669 REAL(kind=dp) :: distance0, distance1, distance2, dvmax, &
672 cpassert(
ASSOCIATED(coords))
675 IF (i == neb_env%number_of_replica)
RETURN
677 SELECT CASE (neb_env%id_type)
681 CALL neb_replica_distance(coords=coords, i0=i, i=i - 1, distance=distance1, iw=iw, &
683 CALL neb_replica_distance(coords=coords, i0=i + 1, i=i, distance=distance2, iw=iw, &
685 tangent(:) = (coords%wrk(:, i) - coords%wrk(:, i - 1))/distance1 + &
686 (coords%wrk(:, i + 1) - coords%wrk(:, i))/distance2
687 CASE (do_it_neb, do_ci_neb, do_d_neb)
688 IF ((energies(i + 1) .GT. energies(i)) .AND. (energies(i) .GT. (energies(i - 1))))
THEN
689 tangent(:) = coords%wrk(:, i + 1) - coords%wrk(:, i)
690 ELSE IF ((energies(i + 1) .LT. energies(i)) .AND. (energies(i) .LT. (energies(i - 1))))
THEN
691 tangent(:) = coords%wrk(:, i) - coords%wrk(:, i - 1)
693 dvmax = max(abs(energies(i + 1) - energies(i)), abs(energies(i - 1) - energies(i)))
694 dvmin = min(abs(energies(i + 1) - energies(i)), abs(energies(i - 1) - energies(i)))
695 IF (energies(i + 1) .GE. energies(i - 1))
THEN
696 tangent(:) = (coords%wrk(:, i + 1) - coords%wrk(:, i))*dvmax + (coords%wrk(:, i) - coords%wrk(:, i - 1))*dvmin
698 tangent(:) = (coords%wrk(:, i + 1) - coords%wrk(:, i))*dvmin + (coords%wrk(:, i) - coords%wrk(:, i - 1))*dvmax
705 distance0 = sqrt(dot_product(tangent(:), tangent(:)))
706 IF (distance0 /= 0.0_dp) tangent(:) = tangent(:)/distance0
707 END SUBROUTINE get_tangent
721 RECURSIVE SUBROUTINE get_neb_force(neb_env, tangent, coords, i, forces, tag, Mmatrix, &
723 TYPE(neb_type),
POINTER :: neb_env
724 REAL(kind=dp),
DIMENSION(:),
INTENT(IN) :: tangent
725 TYPE(neb_var_type),
POINTER :: coords
726 INTEGER,
INTENT(IN) :: i
727 TYPE(neb_var_type),
POINTER :: forces
728 INTEGER,
INTENT(IN),
OPTIONAL :: tag
729 REAL(kind=dp),
DIMENSION(:, :),
INTENT(IN) :: mmatrix
730 INTEGER,
INTENT(IN) :: iw
732 INTEGER :: j, my_tag, nsize_wrk
733 REAL(kind=dp) :: distance1, distance2,
fac, tmp
734 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:) :: dtmp1, wrk
736 my_tag = neb_env%id_type
737 IF (
PRESENT(tag)) my_tag = tag
738 cpassert(
ASSOCIATED(forces))
739 cpassert(
ASSOCIATED(coords))
740 nsize_wrk = coords%size_wrk(1)
744 CASE (do_b_neb, do_it_neb, do_ci_neb, do_d_neb)
745 IF (i == neb_env%number_of_replica)
RETURN
750 CALL cite_reference(e2002)
754 ALLOCATE (wrk(nsize_wrk))
756 CALL neb_replica_distance(coords=coords, i0=i - 1, i=i, distance=distance1, iw=iw, &
758 tmp = distance1 - neb_env%avg_distance
759 neb_env%spring_energy = neb_env%spring_energy + 0.5_dp*neb_env%k*tmp**2
762 CALL cite_reference(elber1987)
765 ALLOCATE (dtmp1(nsize_wrk))
767 tmp = distance1 - neb_env%avg_distance
768 dtmp1(:) = 1.0_dp/distance1*(coords%wrk(:, i) - coords%wrk(:, i - 1))
769 wrk(:) = neb_env%k*tmp*dtmp1
770 forces%wrk(:, i) = forces%wrk(:, i) - wrk
771 forces%wrk(:, i - 1) = forces%wrk(:, i - 1) + wrk
773 fac = 1.0_dp/(neb_env%avg_distance*real(neb_env%number_of_replica - 1, kind=dp))
774 wrk(:) = neb_env%k*
fac*(coords%wrk(:, i) - coords%wrk(:, i - 1))
776 DO j = 2, neb_env%number_of_replica
777 CALL neb_replica_distance(coords=coords, i0=j - 1, i=j, distance=distance1, iw=iw, &
779 tmp = tmp + distance1 - neb_env%avg_distance
781 forces%wrk(:, i) = forces%wrk(:, i) + wrk*tmp
782 forces%wrk(:, i - 1) = forces%wrk(:, i - 1) - wrk*tmp
786 CALL cite_reference(jonsson1998)
787 wrk(:) = (coords%wrk(:, i + 1) - 2.0_dp*coords%wrk(:, i) + coords%wrk(:, i - 1))
788 tmp = neb_env%k*dot_product(wrk, tangent)
789 wrk(:) = forces%wrk(:, i) - dot_product_band(neb_env, forces%wrk(:, i), tangent, mmatrix)*tangent
790 forces%wrk(:, i) = wrk + tmp*tangent
793 CALL cite_reference(jonsson2000_1)
794 CALL neb_replica_distance(coords=coords, i0=i, i=i + 1, distance=distance1, iw=iw, &
796 CALL neb_replica_distance(coords=coords, i0=i - 1, i=i, distance=distance2, iw=iw, &
798 tmp = neb_env%k*(distance1 - distance2)
799 wrk(:) = forces%wrk(:, i) - dot_product_band(neb_env, forces%wrk(:, i), tangent, mmatrix)*tangent
800 forces%wrk(:, i) = wrk + tmp*tangent
803 CALL cite_reference(jonsson2000_2)
804 IF (neb_env%istep <= neb_env%nsteps_it .OR. i /= neb_env%nr_HE_image)
THEN
805 CALL get_neb_force(neb_env, tangent, coords, i, forces, do_it_neb, mmatrix, iw)
807 wrk(:) = forces%wrk(:, i)
808 tmp = -2.0_dp*dot_product_band(neb_env, wrk, tangent, mmatrix)
809 forces%wrk(:, i) = wrk + tmp*tangent
813 CALL cite_reference(wales2004)
814 ALLOCATE (dtmp1(nsize_wrk))
815 dtmp1(:) = forces%wrk(:, i) - dot_product_band(neb_env, forces%wrk(:, i), tangent, mmatrix)*tangent
816 forces%wrk(:, i) = dtmp1
817 tmp = sqrt(dot_product(dtmp1, dtmp1))
818 dtmp1(:) = dtmp1(:)/tmp
821 wrk(:) = (coords%wrk(:, i + 1) - 2.0_dp*coords%wrk(:, i) + coords%wrk(:, i - 1))
822 tmp = dot_product(wrk, dtmp1)
823 dtmp1(:) = neb_env%k*(wrk(:) - tmp*dtmp1(:))
824 forces%wrk(:, i) = forces%wrk(:, i) + dtmp1(:)
828 END SUBROUTINE get_neb_force
841 FUNCTION dot_product_band(neb_env, array1, array2, array3)
RESULT(value)
842 TYPE(neb_type),
POINTER :: neb_env
843 REAL(kind=dp),
DIMENSION(:),
INTENT(IN) :: array1, array2
844 REAL(kind=dp),
DIMENSION(:, :),
INTENT(IN) :: array3
845 REAL(kind=dp) ::
value
850 IF (neb_env%use_colvar)
THEN
851 nsize_int = neb_env%nsize_int
852 check = ((
SIZE(array1) /=
SIZE(array2)) .OR. &
853 (
SIZE(array1) /= nsize_int) .OR. &
854 (
SIZE(array3) /= nsize_int*nsize_int))
857 value = dot_product(matmul(reshape(array3, (/nsize_int, nsize_int/)), array1), array2)
859 value = dot_product(array1, array2)
861 END FUNCTION dot_product_band
876 distances, number_of_replica)
877 LOGICAL,
INTENT(IN) :: rotate_frames
878 TYPE(particle_type),
DIMENSION(:),
OPTIONAL, &
879 POINTER :: particle_set
880 TYPE(neb_var_type),
POINTER :: coords, vels
881 INTEGER,
INTENT(IN) :: iw
882 REAL(kind=dp),
DIMENSION(:),
OPTIONAL :: distances
883 INTEGER,
INTENT(IN) :: number_of_replica
885 INTEGER :: i, k, kind
887 REAL(kind=dp) :: xtmp
888 REAL(kind=dp),
DIMENSION(3) :: tmp
889 REAL(kind=dp),
DIMENSION(3, 3) :: rot
895 DO i = 2, number_of_replica
898 IF (rotate_frames .AND. (coords%in_use == do_band_cartesian))
THEN
899 CALL rmsd3(particle_set, coords%xyz(:, i), coords%xyz(:, i - 1), iw, &
900 rotate=.true., rot=rot)
902 DO k = 1,
SIZE(vels%xyz, 1)/3
904 tmp = vels%xyz(kind + 1:kind + 3, i)
905 vels%xyz(kind + 1:kind + 3, i) = matmul(transpose(rot), tmp)
908 IF (
PRESENT(distances))
THEN
909 check =
SIZE(distances) == (number_of_replica - 1)
911 xtmp = dot_product(coords%wrk(:, i) - coords%wrk(:, i - 1), &
912 coords%wrk(:, i) - coords%wrk(:, i - 1))
913 distances(i - 1) = sqrt(xtmp)
929 coords, sline, distances)
931 LOGICAL,
INTENT(IN) :: reparametrize_frames
932 INTEGER,
INTENT(IN) :: spline_order
933 REAL(kind=dp),
INTENT(IN) :: smoothing
934 REAL(kind=dp),
DIMENSION(:, :),
POINTER :: coords, sline
935 REAL(kind=dp),
DIMENSION(:) :: distances
938 REAL(kind=dp) :: avg_distance, xtmp
939 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:, :) :: tmp_coords
941 IF (reparametrize_frames)
THEN
942 ALLOCATE (tmp_coords(
SIZE(coords, 1),
SIZE(coords, 2)))
943 tmp_coords(:, :) = coords
945 DO i = 2,
SIZE(coords, 2) - 1
946 coords(:, i) = tmp_coords(:, i)*(1.0_dp - 2.0_dp*smoothing) + &
947 tmp_coords(:, i - 1)*smoothing + tmp_coords(:, i + 1)*smoothing
949 sline = coords - tmp_coords + sline
950 tmp_coords(:, :) = coords
952 SELECT CASE (spline_order)
955 DO i = 2,
SIZE(coords, 2)
956 xtmp = dot_product(coords(:, i) - coords(:, i - 1), coords(:, i) - coords(:, i - 1))
957 distances(i - 1) = sqrt(xtmp)
959 avg_distance = sum(distances)/real(
SIZE(coords, 2) - 1, kind=dp)
961 DO i = 2,
SIZE(coords, 2) - 1
963 DO j = 1,
SIZE(coords, 2) - 1
964 xtmp = xtmp + distances(j)
965 IF (xtmp > avg_distance*real(i - 1, kind=dp))
THEN
966 xtmp = (xtmp - avg_distance*real(i - 1, kind=dp))/distances(j)
967 coords(:, i) = (1.0_dp - xtmp)*tmp_coords(:, j + 1) + xtmp*tmp_coords(:, j)
973 DO i = 2,
SIZE(coords, 2)
974 xtmp = dot_product(coords(:, i) - coords(:, i - 1), coords(:, i) - coords(:, i - 1))
975 distances(i - 1) = sqrt(xtmp)
978 cpwarn(
"String Method: Spline order greater than 1 not implemented.")
980 sline = coords - tmp_coords + sline
981 DEALLOCATE (tmp_coords)
994 TYPE(neb_type),
POINTER :: neb_env
995 TYPE(neb_var_type),
POINTER :: dcoords, forces
998 CHARACTER(LEN=3),
DIMENSION(4) :: labels
1000 REAL(kind=dp) :: max_dr, max_force, my_max_dr, &
1001 my_max_force, my_rms_dr, my_rms_force, &
1003 TYPE(cp_logger_type),
POINTER :: logger
1004 TYPE(section_vals_type),
POINTER :: cc_section
1006 NULLIFY (logger, cc_section)
1007 logger => cp_get_default_logger()
1008 cc_section => section_vals_get_subs_vals(neb_env%neb_section,
"CONVERGENCE_CONTROL")
1009 CALL section_vals_val_get(cc_section,
"MAX_DR", r_val=max_dr)
1010 CALL section_vals_val_get(cc_section,
"MAX_FORCE", r_val=max_force)
1011 CALL section_vals_val_get(cc_section,
"RMS_DR", r_val=rms_dr)
1012 CALL section_vals_val_get(cc_section,
"RMS_FORCE", r_val=rms_force)
1015 my_max_dr = maxval(abs(dcoords%wrk))
1016 my_max_force = maxval(abs(forces%wrk))
1017 my_rms_dr = sqrt(sum(dcoords%wrk*dcoords%wrk)/real(
SIZE(dcoords%wrk, 1)*
SIZE(dcoords%wrk, 2), kind=dp))
1018 my_rms_force = sqrt(sum(forces%wrk*forces%wrk)/real(
SIZE(forces%wrk, 1)*
SIZE(forces%wrk, 2), kind=dp))
1019 IF (my_max_dr < max_dr) labels(1) =
"YES"
1020 IF (my_max_force < max_force) labels(2) =
"YES"
1021 IF (my_rms_dr < rms_dr) labels(3) =
"YES"
1022 IF (my_rms_force < rms_force) labels(4) =
"YES"
1023 IF (all(labels ==
"YES")) converged = .true.
1025 iw = cp_print_key_unit_nr(logger, neb_env%neb_section,
"CONVERGENCE_INFO", &
1026 extension=
".nebLog")
1029 WRITE (iw, fmt=
'(A,A)')
' **************************************', &
1030 '*****************************************'
1031 WRITE (iw, fmt=
'(1X,A,2X,F8.5,5X,"[",F8.5,"]",1X,T76,"(",A,")")') &
1032 'RMS DISPLACEMENT =', my_rms_dr, rms_dr, labels(3), &
1033 'MAX DISPLACEMENT =', my_max_dr, max_dr, labels(1), &
1034 'RMS FORCE =', my_rms_force, rms_force, labels(4), &
1035 'MAX FORCE =', my_max_force, max_force, labels(2)
1036 WRITE (iw, fmt=
'(A,A)')
' **************************************', &
1037 '*****************************************'
1039 CALL cp_print_key_finished_output(iw, logger, neb_env%neb_section, &
static GRID_HOST_DEVICE double fac(const int i)
Factorial function, e.g. fac(5) = 5! = 120.
collects all references to literature in CP2K as new algorithms / method are included from literature...
integer, save, public e2002
integer, save, public elber1987
integer, save, public jonsson2000_1
integer, save, public jonsson1998
integer, save, public jonsson2000_2
integer, save, public wales2004
evaluations of colvar for internal coordinates schemes
subroutine, public eval_colvar(force_env, coords, cvalues, Bmatrix, MassI, Amatrix)
Computes the values of colvars and the Wilson matrix B and its invers A.
subroutine, public get_clv_force(force_env, forces, coords, nsize_xyz, nsize_int, cvalues, Mmatrix)
Computes the forces in the frame of collective variables, and additional also the local metric tensor...
subroutine, public set_colvars_target(targets, force_env)
Set the value of target for constraints/restraints.
various routines to log and control the output. The idea is that decisions about where to log should ...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Utility routines to read data from files. Kept as close as possible to the old parser because.
subroutine, public parser_get_next_line(parser, nline, at_end)
Read the next input line and broadcast the input information. Skip (nline-1) lines and skip also all ...
Utility routines to read data from files. Kept as close as possible to the old parser because.
subroutine, public parser_release(parser)
releases the parser
subroutine, public parser_create(parser, file_name, unit_nr, para_env, end_section_label, separator_chars, comment_char, continuation_char, quote_char, section_char, parse_white_lines, initial_variables, apply_preprocessing)
Start a parser run. Initial variables allow to @SET stuff before opening the file.
interface to use cp2k as library
subroutine, public f_env_add_defaults(f_env_id, f_env, handle)
adds the default environments of the f_env to the stack of the defaults, and returns a new error and ...
subroutine, public get_energy(env_id, e_pot, ierr)
returns the energy of the last configuration calculated
subroutine, public get_pos(env_id, pos, n_el, ierr)
gets the positions of the particles
subroutine, public set_pos(env_id, new_pos, n_el, ierr)
sets the positions of the particles
subroutine, public f_env_rm_defaults(f_env, ierr, handle)
removes the default environments of the f_env to the stack of the defaults, and sets ierr accordingly...
subroutine, public get_force(env_id, frc, n_el, ierr)
gets the forces of the particles
Interface for the force calculations.
recursive subroutine, public force_env_calc_energy_force(force_env, calc_force, consistent_energies, skip_external_control, eval_energy_forces, require_consistent_energy_force, linres, calc_stress_tensor)
Interface routine for force and energy calculations.
Interface for the force calculations.
recursive subroutine, public force_env_get(force_env, in_use, fist_env, qs_env, meta_env, fp_env, subsys, para_env, potential_energy, additional_potential, kinetic_energy, harmonic_shell, kinetic_shell, cell, sub_force_env, qmmm_env, qmmmx_env, eip_env, pwdft_env, globenv, input, force_env_section, method_name_id, root_section, mixed_env, nnp_env, embed_env)
returns various attributes about the force environment
performs geometry optimization
subroutine, public cp_geo_opt(force_env, globenv, eval_opt_geo, rm_restart_info)
Main driver to perform geometry optimization.
Define type storing the global information of a run. Keep the amount of stored data small....
Defines the basic variable types.
integer, parameter, public dp
integer, parameter, public default_string_length
integer, parameter, public default_path_length
Perform a molecular dynamics (MD) run using QUICKSTEP.
subroutine, public qs_mol_dyn(force_env, globenv, averages, rm_restart_info, hmc_e_initial, hmc_e_final, mdctrl)
Main driver module for Molecular Dynamics.
Interface to the message passing library MPI.
I/O Module for Nudged Elastic Band Calculation.
subroutine, public dump_replica_coordinates(particle_set, coords, i_rep, ienum, iw, use_colvar)
dump coordinates of a replica NEB
subroutine, public handle_band_file_names(rep_env, irep, n_rep, istep)
Handles the correct file names during a band calculation.
Module with utility to perform MD Nudged Elastic Band Calculation.
subroutine, public neb_initialize_velocity(vels, neb_section, particle_set, i_rep, iw, globenv, neb_env)
Initialize velocities of replica in an MD optimized algorithm within NEB.
Typo for Nudged Elastic Band Calculation.
Module with utility for Nudged Elastic Band Calculation.
subroutine, public neb_calc_energy_forces(rep_env, neb_env, coords, energies, forces, particle_set, output_unit)
Driver to compute energy and forces within a NEB, Based on the use of the replica_env.
logical function, public check_convergence(neb_env, Dcoords, forces)
Checks for convergence criteria during a NEB run.
subroutine, public build_replica_coords(neb_section, particle_set, coords, vels, neb_env, iw, globenv, para_env)
Constructs or Read the coordinates for all replica.
subroutine, public reorient_images(rotate_frames, particle_set, coords, vels, iw, distances, number_of_replica)
Reorient iteratively all images of the NEB chain in order to have always the smaller RMSD between two...
subroutine, public reparametrize_images(reparametrize_frames, spline_order, smoothing, coords, sline, distances)
Reparametrization of the replica for String Method with splines.
Define the data structure for the particle information.
Definition of physical constants:
real(kind=dp), parameter, public bohr
methods to setup replicas of the same system differing only by atom positions and velocities (as used...
subroutine, public rep_env_calc_e_f(rep_env, calc_f)
evaluates the forces
types used to handle many replica of the same system that differ only in atom positions,...
subroutine, public rep_env_sync(rep_env, vals)
sends the data from each replica to all the other on replica j/=i data from replica i overwrites val(...
Defines functions to perform rmsd in 3D.
subroutine, public rmsd3(particle_set, r, r0, output_unit, weights, my_val, rotate, transl, rot, drmsd3)
Computes the RMSD in 3D. Provides also derivatives.