83#include "../base/base_uses.f90"
102 SUBROUTINE cp_eval_at(gopt_env, x, f, gradient, master, &
103 final_evaluation, para_env)
109 TYPE(gopt_f_type),
POINTER :: gopt_env
110 REAL(KIND=
dp),
DIMENSION(:),
POINTER :: x
111 REAL(KIND=
dp),
INTENT(out),
OPTIONAL :: f
112 REAL(KIND=
dp),
DIMENSION(:),
OPTIONAL, &
114 INTEGER,
INTENT(IN) :: master
115 LOGICAL,
INTENT(IN),
OPTIONAL :: final_evaluation
116 TYPE(mp_para_env_type),
POINTER :: para_env
122 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'bfgs_optimizer'
123 LOGICAL,
PARAMETER :: debug_this_module = .true.
140 RECURSIVE SUBROUTINE geoopt_bfgs(force_env, gopt_param, globenv, geo_section, gopt_env, x0)
146 TYPE(gopt_f_type),
POINTER :: gopt_env
147 REAL(kind=dp),
DIMENSION(:),
POINTER :: x0
149 CHARACTER(len=*),
PARAMETER :: routinen =
'geoopt_bfgs'
150 REAL(kind=dp),
PARAMETER :: one = 1.0_dp, zero = 0.0_dp
152 CHARACTER(LEN=5) :: wildcard
153 CHARACTER(LEN=default_path_length) :: hes_filename
154 INTEGER :: handle, hesunit_read, indf, info, &
155 iter_nr, its, maxiter, ndf, nfree, &
157 LOGICAL :: conv, hesrest, ionode, shell_present, &
158 should_stop, use_mod_hes, use_rfo
159 REAL(kind=dp) :: ediff, emin, eold, etot, pred, rad, rat, &
160 step, t_diff, t_now, t_old
161 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:) :: dg, dr, dx, eigval, gold, work, xold
162 REAL(kind=dp),
DIMENSION(:),
POINTER :: g
166 TYPE(
cp_fm_type) :: eigvec_mat, hess_mat, hess_tmp
168 TYPE(mp_para_env_type),
POINTER :: para_env
173 NULLIFY (logger, g, blacs_env, spgr)
175 para_env => force_env%para_env
176 root_section => force_env%root_section
177 spgr => gopt_env%spgr
180 CALL timeset(routinen, handle)
183 ionode = para_env%is_source()
184 maxiter = gopt_param%max_iter
191 SELECT CASE (gopt_env%type_id)
193 cpabort(
"BFGS method not yet working with DIMER")
201 IF (output_unit > 0)
THEN
203 WRITE (unit=output_unit, fmt=
"(/,T2,A,T78,A3)") &
204 "BFGS| Use rational function optimization for step estimation: ",
"YES"
206 WRITE (unit=output_unit, fmt=
"(/,T2,A,T78,A3)") &
207 "BFGS| Use rational function optimization for step estimation: ",
" NO"
209 IF (use_mod_hes)
THEN
210 WRITE (unit=output_unit, fmt=
"(T2,A,T78,A3)") &
211 "BFGS| Use model Hessian for initial guess: ",
"YES"
213 WRITE (unit=output_unit, fmt=
"(T2,A,T78,A3)") &
214 "BFGS| Use model Hessian for initial guess: ",
" NO"
217 WRITE (unit=output_unit, fmt=
"(T2,A,T78,A3)") &
218 "BFGS| Restart Hessian: ",
"YES"
220 WRITE (unit=output_unit, fmt=
"(T2,A,T78,A3)") &
221 "BFGS| Restart Hessian: ",
" NO"
223 WRITE (unit=output_unit, fmt=
"(T2,A,T61,F20.3)") &
224 "BFGS| Trust radius: ", rad
228 nfree = gopt_env%nfree
230 CALL cp_warn(__location__, &
231 "The dimension of the Hessian matrix ("// &
232 trim(adjustl(
cp_to_string(ndf)))//
") is greater than 3000. "// &
233 "The diagonalisation of the full Hessian matrix needed for BFGS "// &
234 "is computationally expensive. You should consider to use the linear "// &
235 "scaling variant L-BFGS instead.")
239 globenv%blacs_repeatable)
241 nrow_global=ndf, ncol_global=ndf)
242 CALL cp_fm_create(hess_mat, fm_struct_hes, name=
"hess_mat")
243 CALL cp_fm_create(hess_tmp, fm_struct_hes, name=
"hess_tmp")
244 CALL cp_fm_create(eigvec_mat, fm_struct_hes, name=
"eigvec_mat")
245 ALLOCATE (eigval(ndf))
251 IF (use_mod_hes)
THEN
252 IF (shell_present)
THEN
253 CALL cp_warn(__location__, &
254 "No model Hessian is available for core-shell models. "// &
255 "A unit matrix is used as the initial Hessian.")
256 use_mod_hes = .false.
259 CALL cp_warn(__location__, &
260 "No model Hessian is available for cell optimizations. "// &
261 "A unit matrix is used as the initial Hessian.")
262 use_mod_hes = .false.
266 IF (use_mod_hes)
THEN
268 CALL construct_initial_hess(gopt_env%force_env, hess_mat)
274 IF (output_unit > 0)
WRITE (output_unit, *) &
275 "BFGS: Matrix diagonalization failed, using unity as model Hessian."
277 DO its = 1,
SIZE(eigval)
278 IF (eigval(its) < 0.1_dp) eigval(its) = 0.1_dp
282 CALL parallel_gemm(
"N",
"T", ndf, ndf, ndf, one, hess_tmp, eigvec_mat, zero, hess_mat)
311 IF (spgr%keep_space_group)
THEN
318 CALL cp_iterate(logger%iter_info, increment=0, iter_nr_out=iter_nr)
322 CALL cp_eval_at(gopt_env, x0, etot, g, gopt_env%force_env%para_env%mepos, &
323 .false., gopt_env%force_env%para_env)
326 IF (spgr%keep_space_group)
THEN
334 t_diff = t_now - t_old
336 CALL gopt_f_io_init(gopt_env, output_unit, etot, wildcard=wildcard, its=iter_nr, used_time=t_diff)
337 DO its = iter_nr + 1, maxiter
338 CALL cp_iterate(logger%iter_info, last=(its == maxiter))
343 IF (((its - iter_nr) == 1) .AND. hesrest)
THEN
346 IF (len_trim(hes_filename) == 0)
THEN
348 hes_filename = trim(logger%iter_info%project_name)//
"-BFGS.Hessian"
350 IF (output_unit > 0)
THEN
351 WRITE (unit=output_unit, fmt=
"(/,T2,A)") &
352 "BFGS| Checking for Hessian restart file <"//trim(adjustl(hes_filename))//
">"
354 CALL open_file(file_name=trim(hes_filename), file_status=
"OLD", &
355 file_form=
"UNFORMATTED", file_action=
"READ", unit_number=hesunit_read)
356 IF (output_unit > 0)
THEN
357 WRITE (unit=output_unit, fmt=
"(T2,A)") &
358 "BFGS| Hessian restart file read"
362 IF (ionode)
CALL close_file(unit_number=hesunit_read)
364 IF ((its - iter_nr) > 1)
THEN
366 IF (spgr%keep_space_group)
THEN
372 dx(indf) = x0(indf) - xold(indf)
373 dg(indf) = g(indf) - gold(indf)
376 CALL bfgs(ndf, dx, dg, hess_mat, work, para_env)
379 IF (spgr%keep_space_group)
THEN
386 CALL write_bfgs_hessian(geo_section, hess_mat, logger)
392 IF (spgr%keep_space_group)
THEN
408 IF (output_unit > 0)
WRITE (output_unit, *) &
409 "BFGS: Matrix diagonalization failed, resetting Hessian to unity."
416 CALL set_hes_eig(ndf, eigval, work)
418 CALL rat_fun_opt(ndf, dg, eigval, work, eigvec_mat, g, para_env)
420 CALL geoopt_get_step(ndf, eigval, eigvec_mat, hess_tmp, dr, g, para_env, use_rfo)
423 IF (spgr%keep_space_group)
THEN
427 CALL trust_radius(ndf, step, rad, rat, dr, output_unit)
433 IF (spgr%keep_space_group)
THEN
437 CALL energy_predict(ndf, work, hess_mat, dr, g, conv, pred, para_env)
441 CALL cp_eval_at(gopt_env, x0, etot, g, gopt_env%force_env%para_env%mepos, &
442 .false., gopt_env%force_env%para_env)
447 IF (spgr%keep_space_group)
THEN
453 IF (should_stop)
EXIT
457 t_diff = t_now - t_old
459 CALL gopt_f_io(gopt_env, force_env, root_section, its, etot, output_unit, &
460 eold, emin, wildcard, gopt_param, ndf, dr, g, conv, pred, rat, &
461 step, rad, used_time=t_diff)
463 IF (conv .OR. (its == maxiter))
EXIT
464 IF (etot < emin) emin = etot
465 IF (use_rfo)
CALL update_trust_rad(rat, rad, step, ediff)
468 IF (its == maxiter .AND. (.NOT. conv))
THEN
474 IF (spgr%show_space_group)
THEN
480 CALL cp_iterate(logger%iter_info, last=.true., increment=0)
481 CALL write_bfgs_hessian(geo_section, hess_mat, logger)
483 gopt_env%force_env%para_env, gopt_env%force_env%para_env%mepos, output_unit)
501 "PRINT%PROGRAM_RUN_INFO")
502 CALL timestop(handle)
516 SUBROUTINE rat_fun_opt(ndf, dg, eigval, work, eigvec_mat, g, para_env)
518 INTEGER,
INTENT(IN) :: ndf
519 REAL(kind=dp),
INTENT(INOUT) :: dg(ndf), eigval(ndf), work(ndf)
521 REAL(kind=dp),
INTENT(INOUT) :: g(ndf)
522 TYPE(mp_para_env_type),
OPTIONAL,
POINTER :: para_env
524 CHARACTER(LEN=*),
PARAMETER :: routinen =
'rat_fun_opt'
525 REAL(kind=dp),
PARAMETER :: one = 1.0_dp
527 INTEGER :: handle, i, indf, iref, iter, j, k, l, &
528 maxit, ncol_local, nrow_local
529 INTEGER,
DIMENSION(:),
POINTER :: col_indices, row_indices
530 LOGICAL :: bisec, fail, set
531 REAL(kind=dp) :: fun, fun1, fun2, fun3, fung, lam1, lam2, &
532 ln, lp, ssize, step, stol
533 REAL(kind=dp),
CONTIGUOUS,
DIMENSION(:, :),
POINTER :: local_data
535 CALL timeset(routinen, handle)
545 CALL cp_fm_get_info(eigvec_mat, row_indices=row_indices, col_indices=col_indices, &
546 local_data=local_data, nrow_local=nrow_local, ncol_local=ncol_local)
552 dg(l) = dg(l) + local_data(i, k)*g(j)
555 CALL para_env%sum(dg)
566 IF (eigval(iref) < 0.0_dp) ln = eigval(iref) - 0.01_dp
574 fun = fun + dg(indf)**2/(ln - eigval(indf))
575 fung = fung - dg(indf)**2/(ln - eigval(indf)**2)
581 IF (abs(step) < stol)
GOTO 200
582 IF (iter >= maxit)
EXIT
589 IF (eigval(iref) < 0.0_dp) lam1 = eigval(iref) - 0.01_dp
592 fun1 = fun1 + dg(indf)**2/(lam1 - eigval(indf))
595 step = abs(lam1)/1000.0_dp
596 IF (step < ssize) step = ssize
599 IF (iter > maxit)
THEN
606 lam2 = lam1 - iter*step
608 fun2 = fun2 + eigval(indf)**2/(lam2 - eigval(indf))
611 IF (fun2*fun1 < 0.0_dp)
THEN
615 IF (iter > maxit)
THEN
621 step = (lam1 + lam2)/2
624 fun3 = fun3 + dg(indf)**2/(step - eigval(indf))
628 IF (abs(step - lam2) < stol)
THEN
633 IF (fun3*fun1 < stol)
THEN
643 IF ((ln > eigval(iref)) .OR. ((ln > 0.0_dp) .AND. &
644 (eigval(iref) > 0.0_dp)))
THEN
646 IF (.NOT. bisec)
GOTO 100
654 IF (fail .AND. .NOT. set)
THEN
657 eigval(indf) = eigval(indf)*work(indf)
667 eigval(indf) = eigval(indf) - ln
672 CALL timestop(handle)
674 END SUBROUTINE rat_fun_opt
685 SUBROUTINE bfgs(ndf, dx, dg, hess_mat, work, para_env)
686 INTEGER,
INTENT(IN) :: ndf
687 REAL(kind=dp),
INTENT(INOUT) :: dx(ndf), dg(ndf)
689 REAL(kind=dp),
INTENT(INOUT) :: work(ndf)
690 TYPE(mp_para_env_type),
OPTIONAL,
POINTER :: para_env
692 CHARACTER(LEN=*),
PARAMETER :: routinen =
'bfgs'
693 REAL(kind=dp),
PARAMETER :: one = 1.0_dp, zero = 0.0_dp
695 INTEGER :: handle, i, j, k, l, ncol_local, &
697 INTEGER,
DIMENSION(:),
POINTER :: col_indices, row_indices
698 REAL(kind=dp) :: ddot, dxw, gdx
699 REAL(kind=dp),
CONTIGUOUS,
DIMENSION(:, :),
POINTER :: local_hes
701 CALL timeset(routinen, handle)
703 CALL cp_fm_get_info(hess_mat, row_indices=row_indices, col_indices=col_indices, &
704 local_data=local_hes, nrow_local=nrow_local, ncol_local=ncol_local)
711 work(j) = work(j) + local_hes(i, k)*dx(l)
715 CALL para_env%sum(work)
717 gdx = ddot(ndf, dg, 1, dx, 1)
719 dxw = ddot(ndf, dx, 1, work, 1)
726 local_hes(i, k) = local_hes(i, k) + gdx*dg(j)*dg(l) - &
731 CALL timestop(handle)
741 SUBROUTINE set_hes_eig(ndf, eigval, work)
742 INTEGER,
INTENT(IN) :: ndf
743 REAL(kind=dp),
INTENT(INOUT) :: eigval(ndf), work(ndf)
745 CHARACTER(LEN=*),
PARAMETER :: routinen =
'set_hes_eig'
746 REAL(kind=dp),
PARAMETER :: max_neg = -0.5_dp, max_pos = 5.0_dp, &
747 min_eig = 0.005_dp, one = 1.0_dp
749 INTEGER :: handle, indf
752 CALL timeset(routinen, handle)
755 IF (eigval(indf) < 0.0_dp) neg = .true.
756 IF (eigval(indf) > 1000.0_dp) eigval(indf) = 1000.0_dp
759 IF (eigval(indf) < 0.0_dp)
THEN
760 IF (eigval(indf) < max_neg)
THEN
761 eigval(indf) = max_neg
762 ELSE IF (eigval(indf) > -min_eig)
THEN
763 eigval(indf) = -min_eig
765 ELSE IF (eigval(indf) < 1000.0_dp)
THEN
766 IF (eigval(indf) < min_eig)
THEN
767 eigval(indf) = min_eig
768 ELSE IF (eigval(indf) > max_pos)
THEN
769 eigval(indf) = max_pos
775 IF (eigval(indf) < 0.0_dp)
THEN
782 CALL timestop(handle)
784 END SUBROUTINE set_hes_eig
797 SUBROUTINE geoopt_get_step(ndf, eigval, eigvec_mat, hess_tmp, dr, g, para_env, use_rfo)
799 INTEGER,
INTENT(IN) :: ndf
800 REAL(kind=dp),
INTENT(INOUT) :: eigval(ndf)
801 TYPE(
cp_fm_type),
INTENT(IN) :: eigvec_mat, hess_tmp
802 REAL(kind=dp),
INTENT(INOUT) :: dr(ndf), g(ndf)
803 TYPE(mp_para_env_type),
OPTIONAL,
POINTER :: para_env
806 REAL(kind=dp),
PARAMETER :: one = 1.0_dp, zero = 0.0_dp
808 INTEGER :: i, indf, j, k, l, ncol_local, nrow_local
809 INTEGER,
DIMENSION(:),
POINTER :: col_indices, row_indices
810 REAL(kind=dp),
CONTIGUOUS,
DIMENSION(:, :),
POINTER :: local_data
817 eigval(indf) = one/eigval(indf)
821 eigval(indf) = one/max(0.0001_dp, eigval(indf))
830 CALL parallel_gemm(
"N",
"T", ndf, ndf, ndf, one, hess_tmp, eigvec_mat, zero, tmp)
837 CALL cp_fm_get_info(hess_tmp, row_indices=row_indices, col_indices=col_indices, &
838 local_data=local_data, nrow_local=nrow_local, ncol_local=ncol_local)
845 dr(j) = dr(j) - local_data(i, k)*g(l)
849 CALL para_env%sum(dr)
851 END SUBROUTINE geoopt_get_step
862 SUBROUTINE trust_radius(ndf, step, rad, rat, dr, output_unit)
863 INTEGER,
INTENT(IN) :: ndf
864 REAL(kind=dp),
INTENT(INOUT) :: step, rad, rat, dr(ndf)
865 INTEGER,
INTENT(IN) :: output_unit
867 CHARACTER(LEN=*),
PARAMETER :: routinen =
'trust_radius'
868 REAL(kind=dp),
PARAMETER :: one = 1.0_dp
871 REAL(kind=dp) :: scal
873 CALL timeset(routinen, handle)
875 step = maxval(abs(dr))
876 scal = max(one, rad/step)
880 CALL dscal(ndf, rat, dr, 1)
882 IF (output_unit > 0)
THEN
883 WRITE (unit=output_unit, fmt=
"(/,T2,A,F8.5)") &
884 " Step is scaled; Scaling factor = ", rat
888 CALL timestop(handle)
890 END SUBROUTINE trust_radius
903 SUBROUTINE energy_predict(ndf, work, hess_mat, dr, g, conv, pred, para_env)
905 INTEGER,
INTENT(IN) :: ndf
906 REAL(kind=dp),
INTENT(INOUT) :: work(ndf)
908 REAL(kind=dp),
INTENT(INOUT) :: dr(ndf), g(ndf)
909 LOGICAL,
INTENT(INOUT) :: conv
910 REAL(kind=dp),
INTENT(INOUT) :: pred
911 TYPE(mp_para_env_type),
POINTER :: para_env
913 CHARACTER(LEN=*),
PARAMETER :: routinen =
'energy_predict'
914 REAL(kind=dp),
PARAMETER :: zero = 0.0_dp
916 INTEGER :: handle, i, j, k, l, ncol_local, &
918 INTEGER,
DIMENSION(:),
POINTER :: col_indices, row_indices
919 REAL(kind=dp) :: ddot, ener1, ener2
920 REAL(kind=dp),
CONTIGUOUS,
DIMENSION(:, :),
POINTER :: local_data
922 CALL timeset(routinen, handle)
924 ener1 = ddot(ndf, g, 1, dr, 1)
926 CALL cp_fm_get_info(hess_mat, row_indices=row_indices, col_indices=col_indices, &
927 local_data=local_data, nrow_local=nrow_local, ncol_local=ncol_local)
934 work(j) = work(j) + local_data(i, k)*dr(l)
938 CALL para_env%sum(work)
939 ener2 = ddot(ndf, dr, 1, work, 1)
940 pred = ener1 + 0.5_dp*ener2
942 CALL timestop(handle)
944 END SUBROUTINE energy_predict
953 SUBROUTINE update_trust_rad(rat, rad, step, ediff)
955 REAL(kind=dp),
INTENT(INOUT) :: rat, rad, step, ediff
957 CHARACTER(LEN=*),
PARAMETER :: routinen =
'update_trust_rad'
958 REAL(kind=dp),
PARAMETER :: max_trust = 1.0_dp, min_trust = 0.1_dp
962 CALL timeset(routinen, handle)
964 IF (rat > 4.0_dp)
THEN
965 IF (ediff < 0.0_dp)
THEN
970 ELSE IF (rat > 2.0_dp)
THEN
971 IF (ediff < 0.0_dp)
THEN
976 ELSE IF (rat > 4.0_dp/3.0_dp)
THEN
977 IF (ediff < 0.0_dp)
THEN
982 ELSE IF (rat > 10.0_dp/9.0_dp)
THEN
983 IF (ediff < 0.0_dp)
THEN
988 ELSE IF (rat > 0.9_dp)
THEN
989 IF (ediff < 0.0_dp)
THEN
994 ELSE IF (rat > 0.75_dp)
THEN
995 IF (ediff < 0.0_dp)
THEN
1000 ELSE IF (rat > 0.5_dp)
THEN
1001 IF (ediff < 0.0_dp)
THEN
1006 ELSE IF (rat > 0.25_dp)
THEN
1007 IF (ediff < 0.0_dp)
THEN
1012 ELSE IF (ediff < 0.0_dp)
THEN
1018 rad = max(rad, min_trust)
1019 rad = min(rad, max_trust)
1020 CALL timestop(handle)
1022 END SUBROUTINE update_trust_rad
1032 SUBROUTINE write_bfgs_hessian(geo_section, hess_mat, logger)
1038 CHARACTER(LEN=*),
PARAMETER :: routinen =
'write_bfgs_hessian'
1040 INTEGER :: handle, hesunit
1042 CALL timeset(routinen, handle)
1045 extension=
".Hessian", file_form=
"UNFORMATTED", file_action=
"WRITE", &
1046 file_position=
"REWIND")
1052 CALL timestop(handle)
1054 END SUBROUTINE write_bfgs_hessian
1062 SUBROUTINE construct_initial_hess(force_env, hess_mat)
1067 INTEGER :: i, iat_col, iat_row, iglobal, iind, j, &
1068 jat_row, jglobal, jind, k, natom, &
1069 ncol_local, nrow_local, z
1070 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: at_row
1071 INTEGER,
DIMENSION(:),
POINTER :: col_indices, row_indices
1072 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:, :) :: d_ij, rho_ij
1073 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: r_ij
1074 REAL(kind=dp),
DIMENSION(3, 3) :: alpha, r0
1075 REAL(kind=dp),
CONTIGUOUS,
DIMENSION(:, :),
POINTER :: fixed, local_data
1082 CALL force_env_get(force_env=force_env, subsys=subsys, cell=cell)
1084 particles=particles)
1086 alpha(1, :) = [1._dp, 0.3949_dp, 0.3949_dp]
1087 alpha(2, :) = [0.3494_dp, 0.2800_dp, 0.2800_dp]
1088 alpha(3, :) = [0.3494_dp, 0.2800_dp, 0.1800_dp]
1090 r0(1, :) = [1.35_dp, 2.10_dp, 2.53_dp]
1091 r0(2, :) = [2.10_dp, 2.87_dp, 3.40_dp]
1092 r0(3, :) = [2.53_dp, 3.40_dp, 3.40_dp]
1094 CALL cp_fm_get_info(hess_mat, row_indices=row_indices, col_indices=col_indices, &
1095 local_data=local_data, nrow_local=nrow_local, ncol_local=ncol_local)
1096 natom = particles%n_els
1097 ALLOCATE (at_row(natom))
1098 ALLOCATE (rho_ij(natom, natom))
1099 ALLOCATE (d_ij(natom, natom))
1100 ALLOCATE (r_ij(natom, natom, 3))
1101 ALLOCATE (fixed(3, natom))
1105 CALL hess_mat%matrix_struct%para_env%min(fixed(i, :))
1112 IF (z <= 10) at_row(i) = 2
1113 IF (z <= 2) at_row(i) = 1
1120 r_ij(j, i, :) =
pbc(particles%els(i)%r, particles%els(j)%r, cell)
1121 r_ij(i, j, :) = -r_ij(j, i, :)
1122 d_ij(j, i) = sqrt(dot_product(r_ij(j, i, :), r_ij(j, i, :)))
1123 d_ij(i, j) = d_ij(j, i)
1124 rho_ij(j, i) = exp(alpha(jat_row, iat_row)*(r0(jat_row, iat_row)**2 - d_ij(j, i)**2))
1125 rho_ij(i, j) = rho_ij(j, i)
1128 DO i = 1, ncol_local
1129 iglobal = col_indices(i)
1130 iind = mod(iglobal - 1, 3) + 1
1131 iat_col = (iglobal + 2)/3
1132 IF (iat_col > natom) cycle
1133 DO j = 1, nrow_local
1134 jglobal = row_indices(j)
1135 jind = mod(jglobal - 1, 3) + 1
1136 iat_row = (jglobal + 2)/3
1137 IF (iat_row > natom) cycle
1138 IF (iat_row /= iat_col)
THEN
1139 IF (d_ij(iat_row, iat_col) < 6.0_dp) &
1140 local_data(j, i) = local_data(j, i) + &
1141 angle_second_deriv(r_ij, d_ij, rho_ij, iind, jind, iat_col, iat_row, natom)
1143 local_data(j, i) = local_data(j, i) + &
1144 angle_second_deriv(r_ij, d_ij, rho_ij, iind, jind, iat_col, iat_row, natom)
1146 IF (iat_col /= iat_row)
THEN
1147 IF (d_ij(iat_row, iat_col) < 6.0_dp) &
1148 local_data(j, i) = local_data(j, i) - &
1149 dist_second_deriv(r_ij(iat_col, iat_row, :), &
1150 iind, jind, d_ij(iat_row, iat_col), rho_ij(iat_row, iat_col))
1153 IF (k == iat_col) cycle
1154 IF (d_ij(iat_row, k) < 6.0_dp) &
1155 local_data(j, i) = local_data(j, i) + &
1156 dist_second_deriv(r_ij(iat_col, k, :), &
1157 iind, jind, d_ij(iat_row, k), rho_ij(iat_row, k))
1160 IF (fixed(jind, iat_row) < 0.5_dp .OR. fixed(iind, iat_col) < 0.5_dp)
THEN
1161 local_data(j, i) = 0.0_dp
1162 IF (jind == iind .AND. iat_row == iat_col) local_data(j, i) = 1.0_dp
1172 END SUBROUTINE construct_initial_hess
1183 FUNCTION dist_second_deriv(r1, i, j, d, rho)
RESULT(deriv)
1184 REAL(kind=dp),
DIMENSION(3) :: r1
1186 REAL(kind=dp) :: d, rho, deriv
1188 deriv = 0.45_dp*rho*(r1(i)*r1(j))/d**2
1189 END FUNCTION dist_second_deriv
1203 FUNCTION angle_second_deriv(r_ij, d_ij, rho_ij, idir, jdir, iat_der, jat_der, natom)
RESULT(deriv)
1204 REAL(kind=dp),
DIMENSION(:, :, :) :: r_ij
1205 REAL(kind=dp),
DIMENSION(:, :) :: d_ij, rho_ij
1206 INTEGER :: idir, jdir, iat_der, jat_der, natom
1207 REAL(kind=dp) :: deriv
1209 INTEGER :: i, iat, idr, j, jat, jdr
1210 REAL(kind=dp) :: d12, d23, d31, d_mat(3, 2), denom1, &
1211 denom2, denom3, ka1, ka2, ka3, rho12, &
1212 rho23, rho31, rsst1, rsst2, rsst3
1213 REAL(kind=dp),
DIMENSION(3) :: r12, r23, r31
1216 IF (iat_der == jat_der)
THEN
1218 IF (rho_ij(iat_der, i) < 0.00001) cycle
1220 IF (rho_ij(iat_der, j) < 0.00001) cycle
1221 IF (i == iat_der .OR. j == iat_der) cycle
1222 IF (iat_der < i .OR. iat_der > j)
THEN
1223 r12 = r_ij(iat_der, i, :); r23 = r_ij(i, j, :); r31 = r_ij(j, iat_der, :)
1224 d12 = d_ij(iat_der, i); d23 = d_ij(i, j); d31 = d_ij(j, iat_der)
1225 rho12 = rho_ij(iat_der, i); rho23 = rho_ij(i, j); rho31 = rho_ij(j, iat_der)
1227 r12 = r_ij(iat_der, j, :); r23 = r_ij(j, i, :); r31 = r_ij(i, iat_der, :)
1228 d12 = d_ij(iat_der, j); d23 = d_ij(j, i); d31 = d_ij(i, iat_der)
1229 rho12 = rho_ij(iat_der, j); rho23 = rho_ij(j, i); rho31 = rho_ij(i, iat_der)
1231 ka1 = 0.15_dp*rho12*rho23; ka2 = 0.15_dp*rho23*rho31; ka3 = 0.15_dp*rho31*rho12
1232 rsst1 = dot_product(r12, r23); rsst2 = dot_product(r23, r31); rsst3 = dot_product(r31, r12)
1233 denom1 = 1.0_dp - rsst1**2/(d12**2*d23**2); denom2 = 1.0_dp - rsst2**2/(d23**2*d31**2)
1234 denom3 = 1.0_dp - rsst3**2/(d31**2*d12**2)
1235 denom1 = sign(1.0_dp, denom1)*max(abs(denom1), 0.01_dp)
1236 denom2 = sign(1.0_dp, denom2)*max(abs(denom2), 0.01_dp)
1237 denom3 = sign(1.0_dp, denom3)*max(abs(denom3), 0.01_dp)
1238 d_mat(1, 1) = r23(idir)/(d12*d23) - rsst1*r12(idir)/(d12**3*d23)
1239 d_mat(1, 2) = r23(jdir)/(d12*d23) - rsst1*r12(jdir)/(d12**3*d23)
1240 d_mat(2, 1) = -r23(idir)/(d23*d31) + rsst2*r31(idir)/(d23*d31**3)
1241 d_mat(2, 2) = -r23(jdir)/(d23*d31) + rsst2*r31(jdir)/(d23*d31**3)
1242 d_mat(3, 1) = (r31(idir) - r12(idir))/(d31*d12) + rsst3*r31(idir)/(d31**3*d12) - &
1243 rsst3*r12(idir)/(d31*d12**3)
1244 d_mat(3, 2) = (r31(jdir) - r12(jdir))/(d31*d12) + rsst3*r31(jdir)/(d31**3*d12) - &
1245 rsst3*r12(jdir)/(d31*d12**3)
1246 IF (abs(denom1) <= 0.011_dp) d_mat(1, 1) = 0.0_dp
1247 IF (abs(denom2) <= 0.011_dp) d_mat(2, 1) = 0.0_dp
1248 IF (abs(denom3) <= 0.011_dp) d_mat(3, 1) = 0.0_dp
1249 deriv = deriv + ka1*d_mat(1, 1)*d_mat(1, 2)/denom1 + &
1250 ka2*d_mat(2, 1)*d_mat(2, 2)/denom2 + &
1251 ka3*d_mat(3, 1)*d_mat(3, 2)/denom3
1257 IF (i == iat_der .OR. i == jat_der) cycle
1258 IF (jat_der < iat_der)
THEN
1259 iat = jat_der; jat = iat_der; idr = jdir; jdr = idir
1261 iat = iat_der; jat = jat_der; idr = idir; jdr = jdir
1263 IF (jat < i .OR. iat > i)
THEN
1264 r12 = r_ij(iat, jat, :); r23 = r_ij(jat, i, :); r31 = r_ij(i, iat, :)
1265 d12 = d_ij(iat, jat); d23 = d_ij(jat, i); d31 = d_ij(i, iat)
1266 rho12 = rho_ij(iat, jat); rho23 = rho_ij(jat, i); rho31 = rho_ij(i, iat)
1268 r12 = r_ij(iat, i, :); r23 = r_ij(i, jat, :); r31 = r_ij(jat, iat, :)
1269 d12 = d_ij(iat, i); d23 = d_ij(i, jat); d31 = d_ij(jat, iat)
1270 rho12 = rho_ij(iat, i); rho23 = rho_ij(i, jat); rho31 = rho_ij(jat, iat)
1272 ka1 = 0.15_dp*rho12*rho23; ka2 = 0.15_dp*rho23*rho31; ka3 = 0.15_dp*rho31*rho12
1273 rsst1 = dot_product(r12, r23); rsst2 = dot_product(r23, r31); rsst3 = dot_product(r31, r12)
1274 denom1 = 1.0_dp - rsst1**2/(d12**2*d23**2); denom2 = 1.0_dp - rsst2**2/(d23**2*d31**2)
1275 denom3 = 1.0_dp - rsst3**2/(d31**2*d12**2)
1276 denom1 = sign(1.0_dp, denom1)*max(abs(denom1), 0.01_dp)
1277 denom2 = sign(1.0_dp, denom2)*max(abs(denom2), 0.01_dp)
1278 denom3 = sign(1.0_dp, denom3)*max(abs(denom3), 0.01_dp)
1279 d_mat(1, 1) = r23(idr)/(d12*d23) - rsst1*r12(idr)/(d12**3*d23)
1280 d_mat(2, 1) = -r23(idr)/(d23*d31) + rsst2*r31(idr)/(d23*d31**3)
1281 d_mat(3, 1) = (r31(idr) - r12(idr))/(d31*d12) + rsst3*r31(idr)/(d31**3*d12) - &
1282 rsst3*r12(idr)/(d31*d12**3)
1283 IF (jat < i .OR. iat > i)
THEN
1284 d_mat(1, 2) = (r12(jdr) - r23(jdr))/(d12*d23) + rsst1*r12(jdr)/(d12**3*d23) - &
1285 rsst1*r23(jdr)/(d12*d23**3)
1286 d_mat(2, 2) = r31(jdr)/(d23*d31) - rsst2*r23(jdr)/(d23**3*d31)
1287 d_mat(3, 2) = -r31(jdr)/(d31*d12) + rsst3*r12(jdr)/(d31*d12**3)
1289 d_mat(1, 2) = -r12(jdr)/(d12*d23) + rsst1*r23(jdr)/(d12*d23**3)
1290 d_mat(2, 2) = (r23(jdr) - r31(jdr))/(d23*d31) + rsst2*r23(jdr)/(d23**3*d31) - &
1291 rsst2*r31(jdr)/(d23*d31**3)
1292 d_mat(3, 2) = r12(jdr)/(d31*d12) - rsst3*r31(jdr)/(d31**3*d12)
1294 IF (abs(denom1) <= 0.011_dp) d_mat(1, 1) = 0.0_dp
1295 IF (abs(denom2) <= 0.011_dp) d_mat(2, 1) = 0.0_dp
1296 IF (abs(denom3) <= 0.011_dp) d_mat(3, 1) = 0.0_dp
1298 deriv = deriv + ka1*d_mat(1, 1)*d_mat(1, 2)/denom1 + &
1299 ka2*d_mat(2, 1)*d_mat(2, 2)/denom2 + &
1300 ka3*d_mat(3, 1)*d_mat(3, 2)/denom3
1303 deriv = 0.25_dp*deriv
1305 END FUNCTION angle_second_deriv
subroutine cp_eval_at(gopt_env, x, f, gradient, master, final_evaluation, para_env)
evaluete the potential energy and its gradients using an array with same dimension as the particle_se...
represent a simple array based list of the given type
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Routines for Geometry optimization using BFGS algorithm.
recursive subroutine, public geoopt_bfgs(force_env, gopt_param, globenv, geo_section, gopt_env, x0)
Main driver for BFGS geometry optimizations.
collects all references to literature in CP2K as new algorithms / method are included from literature...
integer, save, public lindh1995
Handles all functions related to the CELL.
subroutine, public fix_atom_control(force_env, w)
allows for fix atom constraints
methods related to the blacs parallel environment
subroutine, public cp_blacs_env_release(blacs_env)
releases the given blacs_env
subroutine, public cp_blacs_env_create(blacs_env, para_env, blacs_grid_layout, blacs_repeatable, row_major, grid_2d)
allocates and initializes a type that represent a blacs context
Routines to handle the external control of CP2K.
subroutine, public external_control(should_stop, flag, globenv, target_time, start_time, force_check)
External manipulations during a run : when the <PROJECT_NAME>.EXIT_$runtype command is sent the progr...
Utility routines to open and close files. Tracking of preconnections.
subroutine, public open_file(file_name, file_status, file_form, file_action, file_position, file_pad, unit_number, debug, skip_get_unit_number, file_access)
Opens the requested file using a free unit number.
subroutine, public close_file(unit_number, file_status, keep_preconnection)
Close an open file given by its logical unit number. Optionally, keep the file and unit preconnected.
Basic linear algebra operations for full matrices.
subroutine, public cp_fm_column_scale(matrixa, scaling)
scales column i of matrix a with scaling(i)
subroutine, public cp_fm_transpose(matrix, matrixt)
transposes a matrix matrixt = matrix ^ T
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
subroutine, public choose_eigv_solver(matrix, eigenvectors, eigenvalues, info)
Choose the Eigensolver depending on which library is available ELPA seems to be unstable for small sy...
represent the structure of a full matrix
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
represent a full matrix distributed on many processors
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp, nrow, ncol, set_zero)
creates a new full matrix with the given structure
subroutine, public cp_fm_write_unformatted(fm, unit)
...
subroutine, public cp_fm_read_unformatted(fm, unit)
...
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
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,...
integer, parameter, public cp_p_file
subroutine, public cp_iterate(iteration_info, last, iter_nr, increment, iter_nr_out)
adds one to the actual iteration
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
types that represent a subsys, i.e. a part of the system
subroutine, public cp_subsys_get(subsys, ref_count, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell, cell_ref, use_ref_cell)
returns information about various attributes of the given subsys
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, ipi_env)
returns various attributes about the force environment
Define type storing the global information of a run. Keep the amount of stored data small....
contains a functional that calculates the energy and its derivatives for the geometry optimizer
subroutine, public print_geo_opt_header(gopt_env, output_unit, label)
...
subroutine, public gopt_f_io_init(gopt_env, output_unit, opt_energy, wildcard, its, used_time)
Handles the Output during an optimization run.
recursive subroutine, public gopt_f_io_finalize(gopt_env, force_env, x0, conv, its, root_section, para_env, master, output_unit)
Handles the Output at the end of an optimization run.
subroutine, public gopt_f_io(gopt_env, force_env, root_section, its, opt_energy, output_unit, eold, emin, wildcard, gopt_param, ndf, dx, xi, conv, pred, rat, step, rad, used_time)
Handles the Output during an optimization run.
subroutine, public print_geo_opt_nc(gopt_env, output_unit)
...
subroutine, public gopt_f_ii(its, output_unit)
Prints iteration step of the optimization procedure on screen.
contains a functional that calculates the energy and its derivatives for the geometry optimizer
contains typo and related routines to handle parameters controlling the GEO_OPT module
Defines the basic variable types.
integer, parameter, public dp
integer, parameter, public default_path_length
Machine interface based on Fortran 2003 and POSIX.
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
Interface to the message passing library MPI.
basic linear algebra operations for full matrixes
represent a simple array based list of the given type
Space Group Symmetry Type Module (version 1.0, Ferbruary 12, 2021)
Space Group Symmetry Module (version 1.0, January 16, 2020)
subroutine, public print_spgr(spgr)
routine prints Space Group Information.
subroutine, public spgr_apply_rotations_coord(spgr, coord)
routine applies the rotation matrices to the coordinates.
subroutine, public identify_space_group(subsys, geo_section, gopt_env, iunit)
routine indentifies the space group and finds rotation matrices.
subroutine, public spgr_apply_rotations_force(spgr, force)
routine applies the rotation matrices to the forces.
represent a list of objects
Type defining parameters related to the simulation cell.
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
keeps the information about the structure of a full matrix
type of a logger, at the moment it contains just a print level starting at which level it should be l...
represents a system: atoms, molecules, their pos,vel,...
wrapper to abstract the force evaluation of the various methods
contains the initially parsed file and the initial parallel environment
calculates the potential energy of a system, and its derivatives
stores all the informations relevant to an mpi environment
represent a list of objects