108#include "./base/base_uses.f90"
114 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'almo_scf_optimizer'
122 LOGICAL,
PARAMETER :: debug_mode = .false.
123 LOGICAL,
PARAMETER :: safe_mode = .false.
124 LOGICAL,
PARAMETER :: almo_mathematica = .false.
125 INTEGER,
PARAMETER :: hessian_path_reuse = 1, &
126 hessian_path_assemble = 2
145 CHARACTER(len=*),
PARAMETER :: routinen =
'almo_scf_block_diagonal'
147 INTEGER :: handle, iscf, ispin, nspin, unit_nr
148 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: local_nocc_of_domain
149 LOGICAL :: converged, prepare_to_exit, should_stop, &
150 use_diis, use_prev_as_guess
151 REAL(kind=
dp) :: density_rec, energy_diff, energy_new, energy_old, error_norm, &
152 error_norm_ispin, kts_sum, prev_error_norm, t1, t2, true_mixing_fraction
153 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: local_mu
155 DIMENSION(:) :: almo_diis
157 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:) :: matrix_mixing_old_blk
160 CALL timeset(routinen, handle)
164 IF (logger%para_env%is_source())
THEN
172 use_prev_as_guess = .false.
174 nspin = almo_scf_env%nspins
175 ALLOCATE (local_mu(almo_scf_env%ndomains))
176 ALLOCATE (local_nocc_of_domain(almo_scf_env%ndomains))
179 ALLOCATE (matrix_mixing_old_blk(nspin))
180 ALLOCATE (almo_diis(nspin))
183 template=almo_scf_env%matrix_ks_blk(ispin))
185 sample_err=almo_scf_env%matrix_ks_blk(ispin), &
186 sample_var=almo_scf_env%matrix_s_blk(1), &
188 max_length=optimizer%ndiis)
195 prepare_to_exit = .false.
196 true_mixing_fraction = 0.0_dp
197 error_norm = 1.0e+10_dp
199 IF (unit_nr > 0)
THEN
200 WRITE (unit_nr,
'(T2,A,A,A)') repeat(
"-", 20), &
201 " Optimization of block-diagonal ALMOs ", repeat(
"-", 21)
203 WRITE (unit_nr,
'(T2,A13,A6,A23,A14,A14,A9)')
"Method",
"Iter", &
204 "Total Energy",
"Change",
"Convergence",
"Time"
205 WRITE (unit_nr,
'(T2,A)') repeat(
"-", 79)
221 var=almo_scf_env%matrix_ks_blk(ispin), &
222 err=almo_scf_env%matrix_err_blk(ispin))
227 prev_error_norm = error_norm
230 error_norm_ispin =
dbcsr_maxabs(almo_scf_env%matrix_err_blk(ispin))
231 IF (ispin .EQ. 1) error_norm = error_norm_ispin
232 IF (ispin .GT. 1 .AND. error_norm_ispin .GT. error_norm) &
233 error_norm = error_norm_ispin
236 IF (error_norm .LT. almo_scf_env%eps_prev_guess)
THEN
237 use_prev_as_guess = .true.
239 use_prev_as_guess = .false.
244 IF (error_norm .GT. optimizer%eps_error) converged = .false.
248 start_time=qs_env%start_time, &
249 target_time=qs_env%target_time)
250 IF (should_stop .OR. iscf >= optimizer%max_iter .OR. converged)
THEN
251 prepare_to_exit = .true.
252 IF (iscf == 1) energy_new = energy_old
256 IF (optimizer%early_stopping_on .AND. iscf .EQ. 1) &
257 prepare_to_exit = .false.
259 IF (.NOT. prepare_to_exit)
THEN
262 IF (iscf .NE. 1)
THEN
266 extr_var=almo_scf_env%matrix_ks_blk(ispin))
269 true_mixing_fraction = almo_scf_env%mixing_fraction
271 CALL dbcsr_add(almo_scf_env%matrix_ks_blk(ispin), &
272 matrix_mixing_old_blk(ispin), &
273 true_mixing_fraction, &
274 1.0_dp - true_mixing_fraction)
280 CALL dbcsr_copy(matrix_mixing_old_blk(ispin), &
281 almo_scf_env%matrix_ks_blk(ispin))
285 SELECT CASE (almo_scf_env%almo_update_algorithm)
295 local_nocc_of_domain(:) = almo_scf_env%nocc_of_domain(:, ispin)
296 local_mu(:) = almo_scf_env%mu_of_domain(:, ispin)
300 cpabort(
"Density_matrix_sign has not been tested yet")
312 almo_scf_env%mu_of_domain(:, ispin) = local_mu(:)
319 DO ispin = 1, almo_scf_env%nspins
322 overlap=almo_scf_env%matrix_sigma_blk(ispin), &
323 metric=almo_scf_env%matrix_s_blk(1), &
324 retain_locality=.true., &
325 only_normalize=.false., &
326 nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin), &
327 eps_filter=almo_scf_env%eps_filter, &
328 order_lanczos=almo_scf_env%order_lanczos, &
329 eps_lanczos=almo_scf_env%eps_lanczos, &
330 max_iter_lanczos=almo_scf_env%max_iter_lanczos)
337 DO ispin = 1, almo_scf_env%nspins
340 IF (almo_scf_env%smear)
THEN
342 mo_energies=almo_scf_env%mo_energies(:, ispin), &
343 mu_of_domain=almo_scf_env%mu_of_domain(:, ispin), &
344 real_ne_of_domain=almo_scf_env%real_ne_of_domain(:, ispin), &
345 spin_kts=almo_scf_env%kTS(ispin), &
346 smear_e_temp=almo_scf_env%smear_e_temp, &
347 ndomains=almo_scf_env%ndomains, &
348 nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin))
352 p=almo_scf_env%matrix_p(ispin), &
353 eps_filter=almo_scf_env%eps_filter, &
354 orthog_orbs=.false., &
355 nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin), &
356 s=almo_scf_env%matrix_s(1), &
357 sigma=almo_scf_env%matrix_sigma(ispin), &
358 sigma_inv=almo_scf_env%matrix_sigma_inv(ispin), &
359 use_guess=use_prev_as_guess, &
360 smear=almo_scf_env%smear, &
361 algorithm=almo_scf_env%sigma_inv_algorithm, &
362 inverse_accelerator=almo_scf_env%order_lanczos, &
363 inv_eps_factor=almo_scf_env%matrix_iter_eps_error_factor, &
364 eps_lanczos=almo_scf_env%eps_lanczos, &
365 max_iter_lanczos=almo_scf_env%max_iter_lanczos, &
366 para_env=almo_scf_env%para_env, &
367 blacs_env=almo_scf_env%blacs_env)
371 IF (almo_scf_env%nspins == 1)
THEN
374 IF (almo_scf_env%smear)
THEN
375 almo_scf_env%kTS(1) = almo_scf_env%kTS(1)*2.0_dp
379 IF (almo_scf_env%smear)
THEN
380 kts_sum = sum(almo_scf_env%kTS)
387 almo_scf_env%matrix_p, &
388 almo_scf_env%matrix_ks, &
390 almo_scf_env%eps_filter, &
391 almo_scf_env%mat_distr_aos, &
392 smear=almo_scf_env%smear, &
397 energy_diff = energy_new - energy_old
398 energy_old = energy_new
399 almo_scf_env%almo_scf_energy = energy_new
403 IF (unit_nr > 0)
THEN
404 WRITE (unit_nr,
'(T2,A13,I6,F23.10,E14.5,F14.9,F9.2)')
"ALMO SCF DIIS", &
406 energy_new, energy_diff, error_norm, t2 - t1
410 IF (prepare_to_exit)
EXIT
415 IF (almo_scf_env%smear)
THEN
417 CALL dbcsr_dot(almo_scf_env%matrix_p(ispin), almo_scf_env%matrix_s(1), density_rec)
418 IF (unit_nr > 0)
THEN
419 WRITE (unit_nr,
'(T2,A20,F23.10)')
"Electrons recovered:", density_rec
424 IF (.NOT. converged .AND. (.NOT. optimizer%early_stopping_on))
THEN
425 IF (unit_nr > 0)
THEN
426 cpabort(
"SCF for block-diagonal ALMOs not converged!")
434 DEALLOCATE (almo_diis)
435 DEALLOCATE (matrix_mixing_old_blk)
436 DEALLOCATE (local_mu)
437 DEALLOCATE (local_nocc_of_domain)
439 CALL timestop(handle)
459 CHARACTER(len=*),
PARAMETER :: routinen =
'almo_scf_xalmo_eigensolver'
461 INTEGER :: handle, iscf, ispin, nspin, unit_nr
462 LOGICAL :: converged, prepare_to_exit, should_stop
463 REAL(kind=
dp) :: denergy_tot, density_rec, energy_diff, energy_new, energy_old, error_norm, &
464 error_norm_0, kts_sum, spin_factor, t1, t2
465 REAL(kind=
dp),
DIMENSION(2) :: denergy_spin
467 DIMENSION(:) :: almo_diis
469 TYPE(
dbcsr_type) :: matrix_p_almo_scf_converged
471 DIMENSION(:, :) :: submatrix_mixing_old_blk
473 CALL timeset(routinen, handle)
477 IF (logger%para_env%is_source())
THEN
483 nspin = almo_scf_env%nspins
494 matrix_s=almo_scf_env%matrix_s(1), &
495 subm_s_sqrt=almo_scf_env%domain_s_sqrt(:, ispin), &
496 subm_s_sqrt_inv=almo_scf_env%domain_s_sqrt_inv(:, ispin), &
497 dpattern=almo_scf_env%quench_t(ispin), &
498 map=almo_scf_env%domain_map(ispin), &
499 node_of_domain=almo_scf_env%cpu_of_domain)
514 matrix=almo_scf_env%quench_t(ispin), &
515 submatrix=almo_scf_env%domain_t(:, ispin), &
516 distr_pattern=almo_scf_env%quench_t(ispin), &
517 domain_map=almo_scf_env%domain_map(ispin), &
518 node_of_domain=almo_scf_env%cpu_of_domain, &
523 ALLOCATE (submatrix_mixing_old_blk(almo_scf_env%ndomains, nspin))
525 ALLOCATE (almo_diis(nspin))
552 sample_err=almo_scf_env%domain_s_sqrt(:, ispin), &
554 max_length=optimizer%ndiis)
560 prepare_to_exit = .false.
574 d_var=almo_scf_env%domain_ks_xx(:, ispin), &
575 d_err=almo_scf_env%domain_err(:, ispin))
582 error_norm =
dbcsr_maxabs(almo_scf_env%matrix_err_xx(ispin))
585 IF (error_norm .GT. optimizer%eps_error)
THEN
592 start_time=qs_env%start_time, &
593 target_time=qs_env%target_time)
594 IF (should_stop .OR. iscf >= optimizer%max_iter .OR. converged)
THEN
595 prepare_to_exit = .true.
599 IF (optimizer%early_stopping_on .AND. iscf .EQ. 1) &
600 prepare_to_exit = .false.
602 IF (.NOT. prepare_to_exit)
THEN
605 IF (iscf .NE. 1)
THEN
609 almo_scf_env%mixing_fraction, &
610 almo_scf_env%domain_ks_xx(:, ispin), &
611 1.0_dp - almo_scf_env%mixing_fraction, &
612 submatrix_mixing_old_blk(:, ispin), &
618 d_extr_var=almo_scf_env%domain_ks_xx(:, ispin))
625 almo_scf_env%domain_ks_xx(:, ispin), &
626 submatrix_mixing_old_blk(:, ispin), &
637 IF (iscf .EQ. 1)
THEN
639 template=almo_scf_env%matrix_p(ispin))
640 CALL dbcsr_copy(matrix_p_almo_scf_converged, &
641 almo_scf_env%matrix_p(ispin))
645 IF (almo_scf_env%smear)
THEN
647 mo_energies=almo_scf_env%mo_energies(:, ispin), &
648 mu_of_domain=almo_scf_env%mu_of_domain(:, ispin), &
649 real_ne_of_domain=almo_scf_env%real_ne_of_domain(:, ispin), &
650 spin_kts=almo_scf_env%kTS(ispin), &
651 smear_e_temp=almo_scf_env%smear_e_temp, &
652 ndomains=almo_scf_env%ndomains, &
653 nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin))
658 t=almo_scf_env%matrix_t(ispin), &
659 p=almo_scf_env%matrix_p(ispin), &
660 eps_filter=almo_scf_env%eps_filter, &
661 orthog_orbs=.false., &
662 nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin), &
663 s=almo_scf_env%matrix_s(1), &
664 sigma=almo_scf_env%matrix_sigma(ispin), &
665 sigma_inv=almo_scf_env%matrix_sigma_inv(ispin), &
667 smear=almo_scf_env%smear, &
668 algorithm=almo_scf_env%sigma_inv_algorithm, &
669 inverse_accelerator=almo_scf_env%order_lanczos, &
670 inv_eps_factor=almo_scf_env%matrix_iter_eps_error_factor, &
671 eps_lanczos=almo_scf_env%eps_lanczos, &
672 max_iter_lanczos=almo_scf_env%max_iter_lanczos, &
673 para_env=almo_scf_env%para_env, &
674 blacs_env=almo_scf_env%blacs_env)
675 CALL dbcsr_scale(almo_scf_env%matrix_p(ispin), spin_factor)
677 IF (almo_scf_env%smear)
THEN
678 almo_scf_env%kTS(ispin) = almo_scf_env%kTS(ispin)*spin_factor
683 IF (iscf .EQ. 1)
THEN
685 CALL dbcsr_add(matrix_p_almo_scf_converged, &
686 almo_scf_env%matrix_p(ispin), -1.0_dp, 1.0_dp)
687 CALL dbcsr_dot(almo_scf_env%matrix_ks_0deloc(ispin), &
688 matrix_p_almo_scf_converged, &
695 denergy_tot = denergy_tot + denergy_spin(ispin)
760 IF (iscf .EQ. 1)
THEN
761 CALL energy_lowering_report( &
763 ref_energy=almo_scf_env%almo_scf_energy, &
764 energy_lowering=denergy_tot)
766 energy=almo_scf_env%almo_scf_energy, &
767 energy_singles_corr=denergy_tot)
771 IF (.NOT. almo_scf_env%perturbative_delocalization)
THEN
773 IF (almo_scf_env%smear)
THEN
774 kts_sum = sum(almo_scf_env%kTS)
780 almo_scf_env%matrix_p, &
781 almo_scf_env%matrix_ks, &
783 almo_scf_env%eps_filter, &
784 almo_scf_env%mat_distr_aos, &
785 smear=almo_scf_env%smear, &
791 IF (almo_scf_env%perturbative_delocalization)
THEN
794 CALL almo_dm_to_qs_env(qs_env, almo_scf_env%matrix_p, almo_scf_env%mat_distr_aos)
796 prepare_to_exit = .true.
800 energy_diff = energy_new - energy_old
801 energy_old = energy_new
802 almo_scf_env%almo_scf_energy = energy_new
806 IF (unit_nr > 0)
THEN
807 WRITE (unit_nr,
'(T2,A,I6,F20.9,E11.3,E11.3,E11.3,F8.2)')
"ALMO SCF", &
809 energy_new, energy_diff, error_norm, error_norm_0, t2 - t1
815 IF (prepare_to_exit)
EXIT
820 IF (almo_scf_env%smear)
THEN
822 CALL dbcsr_dot(almo_scf_env%matrix_p(ispin), almo_scf_env%matrix_s(1), density_rec)
823 IF (unit_nr > 0)
THEN
824 WRITE (unit_nr,
'(T2,A20,F23.10)')
"Electrons recovered:", density_rec
829 IF (.NOT. converged .AND. .NOT. optimizer%early_stopping_on)
THEN
830 cpabort(
"SCF for ALMOs on overlapping domains not converged!")
837 DEALLOCATE (almo_diis)
838 DEALLOCATE (submatrix_mixing_old_blk)
840 CALL timestop(handle)
866 matrix_t_in, matrix_t_out, assume_t0_q0x, perturbation_only, &
872 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:), &
873 INTENT(INOUT) :: quench_t, matrix_t_in, matrix_t_out
874 LOGICAL,
INTENT(IN) :: assume_t0_q0x, perturbation_only
875 INTEGER,
INTENT(IN),
OPTIONAL :: special_case
877 CHARACTER(len=*),
PARAMETER :: routinen =
'almo_scf_xalmo_pcg'
879 CHARACTER(LEN=20) :: iter_type
880 INTEGER :: cg_iteration, dim_op, fixed_line_search_niter, handle, idim0, ielem, ispin, &
881 iteration, line_search_iteration, max_iter, my_special_case, ndomains, nmo, nspins, &
882 outer_iteration, outer_max_iter, prec_type, reim, unit_nr
883 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: nocc
884 LOGICAL :: blissful_neglect, converged, just_started, line_search, normalize_orbitals, &
885 optimize_theta, outer_prepare_to_exit, penalty_occ_local, penalty_occ_vol, &
886 prepare_to_exit, reset_conjugator, skip_grad, use_guess
887 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: reim_diag, weights, z2
888 REAL(kind=
dp) :: appr_sec_der, beta, denom, denom2, e0, e1, energy_coeff, energy_diff, &
889 energy_new, energy_old, eps_skip_gradients, fval, g0, g1, grad_norm, grad_norm_frob, &
890 line_search_error, localiz_coeff, localization_obj_function, next_step_size_guess, &
891 penalty_amplitude, penalty_func_new, spin_factor, step_size, t1, t2, tempreal
892 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: grad_norm_spin, &
893 penalty_occ_vol_g_prefactor, &
894 penalty_occ_vol_h_prefactor
897 TYPE(
dbcsr_p_type),
DIMENSION(:),
POINTER :: qs_matrix_s
898 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
POINTER :: op_sm_set_almo, op_sm_set_qs
899 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:) :: ftsiginv, grad, m_sig_sqrti_ii, m_t_in_local, &
900 m_theta, prec_vv, prev_grad, prev_minus_prec_grad, prev_step, siginvtftsiginv, st, step, &
901 stsiginv_0, tempnocc, tempnocc_1, tempoccocc
903 DIMENSION(:, :) :: bad_modes_projector_down, domain_r_down
906 CALL timeset(routinen, handle)
909 IF (
PRESENT(special_case)) my_special_case = special_case
913 IF (logger%para_env%is_source())
THEN
919 nspins = almo_scf_env%nspins
923 blissful_neglect = .false.
925 blissful_neglect = .true.
928 IF (unit_nr > 0)
THEN
930 SELECT CASE (my_special_case)
932 WRITE (unit_nr,
'(T2,A,A,A)') repeat(
"-", 20), &
933 " Optimization of block-diagonal ALMOs ", repeat(
"-", 21)
935 WRITE (unit_nr,
'(T2,A,A,A)') repeat(
"-", 20), &
936 " Optimization of fully delocalized MOs ", repeat(
"-", 20)
938 IF (blissful_neglect)
THEN
939 WRITE (unit_nr,
'(T2,A,A,A)') repeat(
"-", 25), &
940 " LCP optimization of XALMOs ", repeat(
"-", 26)
942 WRITE (unit_nr,
'(T2,A,A,A)') repeat(
"-", 27), &
943 " Optimization of XALMOs ", repeat(
"-", 28)
947 WRITE (unit_nr,
'(T2,A13,A6,A23,A14,A14,A9)')
"Method",
"Iter", &
948 "Objective Function",
"Change",
"Convergence",
"Time"
949 WRITE (unit_nr,
'(T2,A)') repeat(
"-", 79)
954 optimize_theta = almo_scf_env%logical05
955 eps_skip_gradients = almo_scf_env%real01
958 energy_coeff = 1.0_dp
959 localiz_coeff = 0.0_dp
960 penalty_amplitude = 0.0_dp
961 penalty_occ_vol = .false.
963 penalty_occ_local = .false.
965 normalize_orbitals = penalty_occ_vol .OR. penalty_occ_local
966 ALLOCATE (penalty_occ_vol_g_prefactor(nspins))
967 ALLOCATE (penalty_occ_vol_h_prefactor(nspins))
968 penalty_occ_vol_g_prefactor(:) = 0.0_dp
969 penalty_occ_vol_h_prefactor(:) = 0.0_dp
970 penalty_func_new = 0.0_dp
973 prec_type = optimizer%preconditioner
976 fixed_line_search_niter = 0
978 IF (nspins == 1)
THEN
984 ALLOCATE (grad_norm_spin(nspins))
985 ALLOCATE (nocc(nspins))
991 ALLOCATE (m_t_in_local(nspins))
994 template=matrix_t_in(ispin), &
995 matrix_type=dbcsr_type_no_symmetry)
996 CALL dbcsr_copy(m_t_in_local(ispin), matrix_t_in(ispin))
1001 ALLOCATE (m_theta(nspins))
1002 DO ispin = 1, nspins
1004 template=matrix_t_out(ispin), &
1005 matrix_type=dbcsr_type_no_symmetry)
1009 IF (penalty_occ_local)
THEN
1012 matrix_s=qs_matrix_s, &
1015 IF (cell%orthorhombic)
THEN
1020 ALLOCATE (weights(6))
1025 ALLOCATE (op_sm_set_qs(2, dim_op))
1026 ALLOCATE (op_sm_set_almo(2, dim_op))
1028 DO idim0 = 1, dim_op
1029 DO reim = 1,
SIZE(op_sm_set_qs, 1)
1030 NULLIFY (op_sm_set_qs(reim, idim0)%matrix)
1031 ALLOCATE (op_sm_set_qs(reim, idim0)%matrix)
1032 CALL dbcsr_copy(op_sm_set_qs(reim, idim0)%matrix, qs_matrix_s(1)%matrix, &
1034 CALL dbcsr_set(op_sm_set_qs(reim, idim0)%matrix, 0.0_dp)
1035 NULLIFY (op_sm_set_almo(reim, idim0)%matrix)
1036 ALLOCATE (op_sm_set_almo(reim, idim0)%matrix)
1037 CALL dbcsr_copy(op_sm_set_almo(reim, idim0)%matrix, almo_scf_env%matrix_s(1), &
1039 CALL dbcsr_set(op_sm_set_almo(reim, idim0)%matrix, 0.0_dp)
1051 m_t_in=m_t_in_local, &
1052 m_t0=almo_scf_env%matrix_t_blk, &
1053 m_quench_t=quench_t, &
1054 m_overlap=almo_scf_env%matrix_s(1), &
1055 m_sigma_tmpl=almo_scf_env%matrix_sigma_inv, &
1057 xalmo_history=almo_scf_env%xalmo_history, &
1058 assume_t0_q0x=assume_t0_q0x, &
1059 optimize_theta=optimize_theta, &
1060 envelope_amplitude=almo_scf_env%envelope_amplitude, &
1061 eps_filter=almo_scf_env%eps_filter, &
1062 order_lanczos=almo_scf_env%order_lanczos, &
1063 eps_lanczos=almo_scf_env%eps_lanczos, &
1064 max_iter_lanczos=almo_scf_env%max_iter_lanczos, &
1065 nocc_of_domain=almo_scf_env%nocc_of_domain)
1067 ndomains = almo_scf_env%ndomains
1068 ALLOCATE (domain_r_down(ndomains, nspins))
1070 ALLOCATE (bad_modes_projector_down(ndomains, nspins))
1073 ALLOCATE (prec_vv(nspins))
1074 ALLOCATE (siginvtftsiginv(nspins))
1075 ALLOCATE (stsiginv_0(nspins))
1076 ALLOCATE (ftsiginv(nspins))
1077 ALLOCATE (st(nspins))
1078 ALLOCATE (prev_grad(nspins))
1079 ALLOCATE (grad(nspins))
1080 ALLOCATE (prev_step(nspins))
1081 ALLOCATE (step(nspins))
1082 ALLOCATE (prev_minus_prec_grad(nspins))
1083 ALLOCATE (m_sig_sqrti_ii(nspins))
1084 ALLOCATE (tempnocc(nspins))
1085 ALLOCATE (tempnocc_1(nspins))
1086 ALLOCATE (tempoccocc(nspins))
1087 DO ispin = 1, nspins
1091 template=almo_scf_env%matrix_ks(ispin), &
1092 matrix_type=dbcsr_type_no_symmetry)
1094 template=almo_scf_env%matrix_sigma(ispin), &
1095 matrix_type=dbcsr_type_no_symmetry)
1097 template=matrix_t_out(ispin), &
1098 matrix_type=dbcsr_type_no_symmetry)
1100 template=matrix_t_out(ispin), &
1101 matrix_type=dbcsr_type_no_symmetry)
1103 template=matrix_t_out(ispin), &
1104 matrix_type=dbcsr_type_no_symmetry)
1106 template=matrix_t_out(ispin), &
1107 matrix_type=dbcsr_type_no_symmetry)
1109 template=matrix_t_out(ispin), &
1110 matrix_type=dbcsr_type_no_symmetry)
1112 template=matrix_t_out(ispin), &
1113 matrix_type=dbcsr_type_no_symmetry)
1115 template=matrix_t_out(ispin), &
1116 matrix_type=dbcsr_type_no_symmetry)
1118 template=matrix_t_out(ispin), &
1119 matrix_type=dbcsr_type_no_symmetry)
1121 template=almo_scf_env%matrix_sigma_inv(ispin), &
1122 matrix_type=dbcsr_type_no_symmetry)
1124 template=matrix_t_out(ispin), &
1125 matrix_type=dbcsr_type_no_symmetry)
1127 template=matrix_t_out(ispin), &
1128 matrix_type=dbcsr_type_no_symmetry)
1130 template=almo_scf_env%matrix_sigma_inv(ispin), &
1131 matrix_type=dbcsr_type_no_symmetry)
1134 CALL dbcsr_set(prev_step(ispin), 0.0_dp)
1137 nfullrows_total=nocc(ispin))
1144 matrix_s=almo_scf_env%matrix_s(1), &
1145 subm_s_inv=almo_scf_env%domain_s_inv(:, ispin), &
1146 dpattern=quench_t(ispin), &
1147 map=almo_scf_env%domain_map(ispin), &
1148 node_of_domain=almo_scf_env%cpu_of_domain)
1151 matrix_s=almo_scf_env%matrix_s(1), &
1152 subm_s_sqrt=almo_scf_env%domain_s_sqrt(:, ispin), &
1153 subm_s_sqrt_inv=almo_scf_env%domain_s_sqrt_inv(:, ispin), &
1154 dpattern=almo_scf_env%quench_t(ispin), &
1155 map=almo_scf_env%domain_map(ispin), &
1156 node_of_domain=almo_scf_env%cpu_of_domain)
1160 IF (assume_t0_q0x)
THEN
1165 almo_scf_env%matrix_s(1), &
1166 almo_scf_env%matrix_t_blk(ispin), &
1167 0.0_dp, st(ispin), &
1168 filter_eps=almo_scf_env%eps_filter)
1171 almo_scf_env%matrix_sigma_inv_0deloc(ispin), &
1172 0.0_dp, stsiginv_0(ispin), &
1173 filter_eps=almo_scf_env%eps_filter)
1179 matrix_t=almo_scf_env%matrix_t_blk(ispin), &
1180 matrix_sigma_inv=almo_scf_env%matrix_sigma_inv(ispin), &
1181 matrix_s=almo_scf_env%matrix_s(1), &
1182 subm_r_down=domain_r_down(:, ispin), &
1183 dpattern=quench_t(ispin), &
1184 map=almo_scf_env%domain_map(ispin), &
1185 node_of_domain=almo_scf_env%cpu_of_domain, &
1186 filter_eps=almo_scf_env%eps_filter)
1192 IF (penalty_occ_local)
THEN
1196 almo_scf_env%matrix_s(1), &
1197 matrix_t_in(ispin), &
1198 0.0_dp, tempnocc(ispin), &
1199 filter_eps=almo_scf_env%eps_filter)
1202 almo_scf_env%matrix_sigma_inv(ispin), &
1203 0.0_dp, tempnocc_1(ispin), &
1204 filter_eps=almo_scf_env%eps_filter)
1206 DO idim0 = 1,
SIZE(op_sm_set_qs, 2)
1207 DO reim = 1,
SIZE(op_sm_set_qs, 1)
1210 op_sm_set_almo(reim, idim0)%matrix, almo_scf_env%mat_distr_aos)
1213 op_sm_set_almo(reim, idim0)%matrix, &
1214 matrix_t_in(ispin), &
1215 0.0_dp, tempnocc(ispin), &
1216 filter_eps=almo_scf_env%eps_filter)
1219 matrix_t_in(ispin), &
1221 0.0_dp, tempoccocc(ispin), &
1222 filter_eps=almo_scf_env%eps_filter)
1225 tempnocc_1(ispin), &
1226 tempoccocc(ispin), &
1227 0.0_dp, tempnocc(ispin), &
1228 filter_eps=almo_scf_env%eps_filter)
1232 tempnocc_1(ispin), &
1233 0.0_dp, op_sm_set_almo(reim, idim0)%matrix, &
1234 filter_eps=almo_scf_env%eps_filter)
1244 outer_max_iter = optimizer%max_iter_outer_loop
1245 outer_prepare_to_exit = .false.
1248 grad_norm_frob = 0.0_dp
1254 max_iter = optimizer%max_iter
1255 prepare_to_exit = .false.
1256 line_search = .false.
1260 line_search_iteration = 0
1263 energy_diff = 0.0_dp
1264 localization_obj_function = 0.0_dp
1265 line_search_error = 0.0_dp
1271 just_started = (iteration .EQ. 0) .AND. (outer_iteration .EQ. 0)
1273 CALL main_var_to_xalmos_and_loss_func( &
1274 almo_scf_env=almo_scf_env, &
1276 m_main_var_in=m_theta, &
1277 m_t_out=matrix_t_out, &
1278 m_sig_sqrti_ii_out=m_sig_sqrti_ii, &
1279 energy_out=energy_new, &
1280 penalty_out=penalty_func_new, &
1281 m_ftsiginv_out=ftsiginv, &
1282 m_siginvtftsiginv_out=siginvtftsiginv, &
1284 m_stsiginv0_in=stsiginv_0, &
1285 m_quench_t_in=quench_t, &
1286 domain_r_down_in=domain_r_down, &
1287 assume_t0_q0x=assume_t0_q0x, &
1288 just_started=just_started, &
1289 optimize_theta=optimize_theta, &
1290 normalize_orbitals=normalize_orbitals, &
1291 perturbation_only=perturbation_only, &
1292 do_penalty=penalty_occ_vol, &
1293 special_case=my_special_case)
1294 IF (penalty_occ_vol)
THEN
1296 energy_new = energy_new + penalty_func_new
1298 DO ispin = 1, nspins
1299 IF (penalty_occ_vol)
THEN
1300 penalty_occ_vol_g_prefactor(ispin) = &
1301 -2.0_dp*penalty_amplitude*spin_factor*nocc(ispin)
1302 penalty_occ_vol_h_prefactor(ispin) = 0.0_dp
1306 localization_obj_function = 0.0_dp
1308 IF (penalty_occ_local)
THEN
1309 DO ispin = 1, nspins
1312 localization_obj_function = 0.0_dp
1313 CALL dbcsr_get_info(almo_scf_env%matrix_sigma_inv(ispin), nfullrows_total=nmo)
1315 ALLOCATE (reim_diag(nmo))
1319 DO idim0 = 1,
SIZE(op_sm_set_qs, 2)
1323 DO reim = 1,
SIZE(op_sm_set_qs, 1)
1329 op_sm_set_almo(reim, idim0)%matrix, &
1330 matrix_t_out(ispin), &
1331 0.0_dp, tempnocc(ispin), &
1332 filter_eps=almo_scf_env%eps_filter)
1335 matrix_t_out(ispin), &
1337 0.0_dp, tempoccocc(ispin), &
1338 filter_eps=almo_scf_env%eps_filter)
1342 CALL group%sum(reim_diag)
1343 z2(:) = z2(:) + reim_diag(:)*reim_diag(:)
1350 fval = -weights(idim0)*log(abs(z2(ielem)))
1352 fval = weights(idim0) - weights(idim0)*abs(z2(ielem))
1354 fval = weights(idim0) - weights(idim0)*sqrt(abs(z2(ielem)))
1356 localization_obj_function = localization_obj_function + fval
1362 DEALLOCATE (reim_diag)
1364 energy_new = energy_new + localiz_coeff*localization_obj_function
1369 DO ispin = 1, nspins
1371 IF (just_started .AND. almo_mathematica)
THEN
1372 cpwarn_if(ispin .GT. 1,
"Mathematica files will be overwritten")
1373 CALL print_mathematica_matrix(almo_scf_env%matrix_s(1),
"matrixS.dat")
1374 CALL print_mathematica_matrix(almo_scf_env%matrix_ks(ispin),
"matrixF.dat")
1375 CALL print_mathematica_matrix(matrix_t_out(ispin),
"matrixT.dat")
1376 CALL print_mathematica_matrix(quench_t(ispin),
"matrixQ.dat")
1382 IF (line_search_iteration .EQ. 0 .AND. iteration .NE. 0) &
1383 CALL dbcsr_copy(prev_grad(ispin), grad(ispin))
1388 skip_grad = (iteration .GT. 0 .AND. &
1389 fixed_line_search_niter .NE. 0 .AND. &
1390 line_search_iteration .NE. fixed_line_search_niter)
1392 IF (.NOT. skip_grad)
THEN
1394 DO ispin = 1, nspins
1396 CALL compute_gradient( &
1397 m_grad_out=grad(ispin), &
1398 m_ks=almo_scf_env%matrix_ks(ispin), &
1399 m_s=almo_scf_env%matrix_s(1), &
1400 m_t=matrix_t_out(ispin), &
1401 m_t0=almo_scf_env%matrix_t_blk(ispin), &
1402 m_siginv=almo_scf_env%matrix_sigma_inv(ispin), &
1403 m_quench_t=quench_t(ispin), &
1404 m_ftsiginv=ftsiginv(ispin), &
1405 m_siginvtftsiginv=siginvtftsiginv(ispin), &
1407 m_stsiginv0=stsiginv_0(ispin), &
1408 m_theta=m_theta(ispin), &
1409 m_sig_sqrti_ii=m_sig_sqrti_ii(ispin), &
1410 domain_s_inv=almo_scf_env%domain_s_inv(:, ispin), &
1411 domain_r_down=domain_r_down(:, ispin), &
1412 cpu_of_domain=almo_scf_env%cpu_of_domain, &
1413 domain_map=almo_scf_env%domain_map(ispin), &
1414 assume_t0_q0x=assume_t0_q0x, &
1415 optimize_theta=optimize_theta, &
1416 normalize_orbitals=normalize_orbitals, &
1417 penalty_occ_vol=penalty_occ_vol, &
1418 penalty_occ_vol_prefactor=penalty_occ_vol_g_prefactor(ispin), &
1419 envelope_amplitude=almo_scf_env%envelope_amplitude, &
1420 eps_filter=almo_scf_env%eps_filter, &
1421 spin_factor=spin_factor, &
1422 special_case=my_special_case, &
1423 penalty_occ_local=penalty_occ_local, &
1424 op_sm_set=op_sm_set_almo, &
1426 energy_coeff=energy_coeff, &
1427 localiz_coeff=localiz_coeff)
1436 IF (blissful_neglect)
THEN
1437 DO ispin = 1, nspins
1440 IF (iteration .EQ. 0)
THEN
1441 CALL compute_preconditioner( &
1442 domain_prec_out=almo_scf_env%domain_preconditioner(:, ispin), &
1443 bad_modes_projector_down_out=bad_modes_projector_down(:, ispin), &
1444 m_prec_out=prec_vv(ispin), &
1445 m_ks=almo_scf_env%matrix_ks(ispin), &
1446 m_s=almo_scf_env%matrix_s(1), &
1447 m_siginv=almo_scf_env%matrix_sigma_inv(ispin), &
1448 m_quench_t=quench_t(ispin), &
1449 m_ftsiginv=ftsiginv(ispin), &
1450 m_siginvtftsiginv=siginvtftsiginv(ispin), &
1452 para_env=almo_scf_env%para_env, &
1453 blacs_env=almo_scf_env%blacs_env, &
1454 nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin), &
1455 domain_s_inv=almo_scf_env%domain_s_inv(:, ispin), &
1456 domain_s_inv_half=almo_scf_env%domain_s_sqrt_inv(:, ispin), &
1457 domain_s_half=almo_scf_env%domain_s_sqrt(:, ispin), &
1458 domain_r_down=domain_r_down(:, ispin), &
1459 cpu_of_domain=almo_scf_env%cpu_of_domain, &
1460 domain_map=almo_scf_env%domain_map(ispin), &
1461 assume_t0_q0x=assume_t0_q0x, &
1462 penalty_occ_vol=penalty_occ_vol, &
1463 penalty_occ_vol_prefactor=penalty_occ_vol_g_prefactor(ispin), &
1464 eps_filter=almo_scf_env%eps_filter, &
1465 neg_thr=optimizer%neglect_threshold, &
1466 spin_factor=spin_factor, &
1467 skip_inversion=.false., &
1468 special_case=my_special_case)
1472 matrix_in=grad(ispin), &
1473 matrix_out=grad(ispin), &
1474 operator1=almo_scf_env%domain_s_inv(:, ispin), &
1475 operator2=bad_modes_projector_down(:, ispin), &
1476 dpattern=quench_t(ispin), &
1477 map=almo_scf_env%domain_map(ispin), &
1478 node_of_domain=almo_scf_env%cpu_of_domain, &
1480 filter_eps=almo_scf_env%eps_filter)
1487 DO ispin = 1, nspins
1490 grad_norm = maxval(grad_norm_spin)
1492 converged = (grad_norm .LE. optimizer%eps_error)
1493 IF (converged .OR. (iteration .GE. max_iter))
THEN
1494 prepare_to_exit = .true.
1497 IF (optimizer%early_stopping_on .AND. just_started) &
1498 prepare_to_exit = .false.
1500 IF (grad_norm .LT. almo_scf_env%eps_prev_guess) &
1504 IF (.NOT. prepare_to_exit)
THEN
1509 IF (iteration .NE. 0)
THEN
1511 IF (fixed_line_search_niter .EQ. 0)
THEN
1515 IF (.NOT. line_search)
THEN
1517 line_search = .true.
1518 line_search_iteration = line_search_iteration + 1
1524 line_search_error = 0.0_dp
1528 DO ispin = 1, nspins
1530 CALL dbcsr_dot(grad(ispin), step(ispin), tempreal)
1531 line_search_error = line_search_error + tempreal
1532 CALL dbcsr_dot(grad(ispin), grad(ispin), tempreal)
1533 denom = denom + tempreal
1534 CALL dbcsr_dot(step(ispin), step(ispin), tempreal)
1535 denom2 = denom2 + tempreal
1541 line_search_error = line_search_error/sqrt(denom)/sqrt(denom2)
1543 IF (abs(line_search_error) .GT. optimizer%lin_search_eps_error)
THEN
1544 line_search = .true.
1545 line_search_iteration = line_search_iteration + 1
1547 line_search = .false.
1548 line_search_iteration = 0
1549 IF (grad_norm .LT. eps_skip_gradients)
THEN
1550 fixed_line_search_niter = abs(almo_scf_env%integer04)
1558 IF (.NOT. line_search)
THEN
1559 line_search = .true.
1560 line_search_iteration = line_search_iteration + 1
1562 IF (line_search_iteration .EQ. fixed_line_search_niter)
THEN
1563 line_search = .false.
1564 line_search_iteration = 0
1565 line_search_iteration = line_search_iteration + 1
1573 IF (line_search)
THEN
1574 energy_diff = 0.0_dp
1576 energy_diff = energy_new - energy_old
1577 energy_old = energy_new
1581 IF (.NOT. line_search)
THEN
1587 cg_iteration = cg_iteration + 1
1590 DO ispin = 1, nspins
1591 CALL dbcsr_copy(prev_step(ispin), step(ispin))
1595 SELECT CASE (prec_type)
1599 CALL newton_grad_to_step( &
1600 optimizer=almo_scf_env%opt_xalmo_newton_pcg_solver, &
1603 m_s=almo_scf_env%matrix_s(:), &
1604 m_ks=almo_scf_env%matrix_ks(:), &
1605 m_siginv=almo_scf_env%matrix_sigma_inv(:), &
1606 m_quench_t=quench_t(:), &
1607 m_ftsiginv=ftsiginv(:), &
1608 m_siginvtftsiginv=siginvtftsiginv(:), &
1610 m_t=matrix_t_out(:), &
1611 m_sig_sqrti_ii=m_sig_sqrti_ii(:), &
1612 domain_s_inv=almo_scf_env%domain_s_inv(:, :), &
1613 domain_r_down=domain_r_down(:, :), &
1614 domain_map=almo_scf_env%domain_map(:), &
1615 cpu_of_domain=almo_scf_env%cpu_of_domain, &
1616 nocc_of_domain=almo_scf_env%nocc_of_domain(:, :), &
1617 para_env=almo_scf_env%para_env, &
1618 blacs_env=almo_scf_env%blacs_env, &
1619 eps_filter=almo_scf_env%eps_filter, &
1620 optimize_theta=optimize_theta, &
1621 penalty_occ_vol=penalty_occ_vol, &
1622 normalize_orbitals=normalize_orbitals, &
1623 penalty_occ_vol_prefactor=penalty_occ_vol_g_prefactor(:), &
1624 penalty_occ_vol_pf2=penalty_occ_vol_h_prefactor(:), &
1625 special_case=my_special_case &
1631 IF (.NOT. blissful_neglect .AND. &
1632 ((just_started .AND. perturbation_only) .OR. &
1633 (iteration .EQ. 0 .AND. (.NOT. perturbation_only))) &
1637 DO ispin = 1, nspins
1638 CALL compute_preconditioner( &
1639 domain_prec_out=almo_scf_env%domain_preconditioner(:, ispin), &
1640 m_prec_out=prec_vv(ispin), &
1641 m_ks=almo_scf_env%matrix_ks(ispin), &
1642 m_s=almo_scf_env%matrix_s(1), &
1643 m_siginv=almo_scf_env%matrix_sigma_inv(ispin), &
1644 m_quench_t=quench_t(ispin), &
1645 m_ftsiginv=ftsiginv(ispin), &
1646 m_siginvtftsiginv=siginvtftsiginv(ispin), &
1648 para_env=almo_scf_env%para_env, &
1649 blacs_env=almo_scf_env%blacs_env, &
1650 nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin), &
1651 domain_s_inv=almo_scf_env%domain_s_inv(:, ispin), &
1652 domain_r_down=domain_r_down(:, ispin), &
1653 cpu_of_domain=almo_scf_env%cpu_of_domain, &
1654 domain_map=almo_scf_env%domain_map(ispin), &
1655 assume_t0_q0x=assume_t0_q0x, &
1656 penalty_occ_vol=penalty_occ_vol, &
1657 penalty_occ_vol_prefactor=penalty_occ_vol_g_prefactor(ispin), &
1658 eps_filter=almo_scf_env%eps_filter, &
1660 spin_factor=spin_factor, &
1661 skip_inversion=.false., &
1662 special_case=my_special_case)
1673 DO ispin = 1, nspins
1678 0.0_dp, step(ispin), &
1679 filter_eps=almo_scf_env%eps_filter)
1686 IF (optimize_theta)
THEN
1687 cpabort(
"theta is NYI")
1690 DO ispin = 1, nspins
1693 matrix_in=grad(ispin), &
1694 matrix_out=step(ispin), &
1695 operator1=almo_scf_env%domain_preconditioner(:, ispin), &
1696 dpattern=quench_t(ispin), &
1697 map=almo_scf_env%domain_map(ispin), &
1698 node_of_domain=almo_scf_env%cpu_of_domain, &
1700 filter_eps=almo_scf_env%eps_filter)
1720 DO ispin = 1, nspins
1730 IF (iteration .EQ. 0)
THEN
1731 reset_conjugator = .true.
1735 IF (.NOT. reset_conjugator)
THEN
1737 CALL compute_cg_beta( &
1739 reset_conjugator=reset_conjugator, &
1740 conjugator=optimizer%conjugator, &
1742 prev_grad=prev_grad(:), &
1744 prev_step=prev_step(:), &
1745 prev_minus_prec_grad=prev_minus_prec_grad(:) &
1750 IF (reset_conjugator)
THEN
1753 IF (unit_nr > 0 .AND. (.NOT. just_started))
THEN
1754 WRITE (unit_nr,
'(T2,A35)')
"Re-setting conjugator to zero"
1756 reset_conjugator = .false.
1761 DO ispin = 1, nspins
1763 CALL dbcsr_copy(prev_minus_prec_grad(ispin), step(ispin))
1770 CALL dbcsr_add(step(ispin), prev_step(ispin), 1.0_dp, beta)
1777 IF (.NOT. line_search)
THEN
1783 DO ispin = 1, nspins
1784 CALL dbcsr_dot(grad(ispin), step(ispin), tempreal)
1787 IF (iteration .EQ. 0)
THEN
1788 step_size = optimizer%lin_search_step_size_guess
1790 IF (next_step_size_guess .LE. 0.0_dp)
THEN
1791 step_size = optimizer%lin_search_step_size_guess
1794 step_size = next_step_size_guess*1.05_dp
1801 next_step_size_guess = step_size
1803 IF (fixed_line_search_niter .EQ. 0)
THEN
1806 DO ispin = 1, nspins
1807 CALL dbcsr_dot(grad(ispin), step(ispin), tempreal)
1812 appr_sec_der = (g1 - g0)/step_size
1817 step_size = -g1/appr_sec_der
1825 appr_sec_der = 2.0*((e1 - e0)/step_size - g0)/step_size
1826 g1 = appr_sec_der*step_size + g0
1832 step_size = -g1/appr_sec_der
1836 next_step_size_guess = next_step_size_guess + step_size
1840 DO ispin = 1, nspins
1841 CALL dbcsr_add(m_theta(ispin), step(ispin), 1.0_dp, step_size)
1846 IF (line_search)
THEN
1853 IF (unit_nr > 0)
THEN
1854 iter_type = trim(
"ALMO SCF "//iter_type)
1855 WRITE (unit_nr,
'(T2,A13,I6,F23.10,E14.5,F14.9,F9.2)') &
1856 iter_type, iteration, &
1857 energy_new, energy_diff, grad_norm, &
1859 IF (penalty_occ_local .OR. penalty_occ_vol)
THEN
1860 WRITE (unit_nr,
'(T2,A25,F23.10)') &
1861 "Energy component:", (energy_new - penalty_func_new - localization_obj_function)
1863 IF (penalty_occ_local)
THEN
1864 WRITE (unit_nr,
'(T2,A25,F23.10)') &
1865 "Localization component:", localization_obj_function
1867 IF (penalty_occ_vol)
THEN
1868 WRITE (unit_nr,
'(T2,A25,F23.10)') &
1869 "Penalty component:", penalty_func_new
1874 IF (penalty_occ_vol)
THEN
1875 almo_scf_env%almo_scf_energy = energy_new - penalty_func_new - localization_obj_function
1877 almo_scf_env%almo_scf_energy = energy_new - localization_obj_function
1883 iteration = iteration + 1
1884 IF (prepare_to_exit)
EXIT
1888 IF (converged .OR. (outer_iteration .GE. outer_max_iter))
THEN
1889 outer_prepare_to_exit = .true.
1892 outer_iteration = outer_iteration + 1
1893 IF (outer_prepare_to_exit)
EXIT
1897 DO ispin = 1, nspins
1898 IF (converged .AND. almo_mathematica)
THEN
1899 cpwarn_if(ispin .GT. 1,
"Mathematica files will be overwritten")
1900 CALL print_mathematica_matrix(matrix_t_out(ispin),
"matrixTf.dat")
1907 CALL wrap_up_xalmo_scf( &
1909 almo_scf_env=almo_scf_env, &
1910 perturbation_in=perturbation_only, &
1911 m_xalmo_in=matrix_t_out, &
1912 m_quench_in=quench_t, &
1913 energy_inout=energy_new)
1917 DO ispin = 1, nspins
1938 DEALLOCATE (tempnocc)
1939 DEALLOCATE (tempnocc_1)
1940 DEALLOCATE (tempoccocc)
1941 DEALLOCATE (prec_vv)
1942 DEALLOCATE (siginvtftsiginv)
1943 DEALLOCATE (stsiginv_0)
1944 DEALLOCATE (ftsiginv)
1946 DEALLOCATE (prev_grad)
1948 DEALLOCATE (prev_step)
1950 DEALLOCATE (prev_minus_prec_grad)
1951 DEALLOCATE (m_sig_sqrti_ii)
1953 DEALLOCATE (domain_r_down)
1954 DEALLOCATE (bad_modes_projector_down)
1956 DEALLOCATE (penalty_occ_vol_g_prefactor)
1957 DEALLOCATE (penalty_occ_vol_h_prefactor)
1958 DEALLOCATE (grad_norm_spin)
1961 DEALLOCATE (m_theta, m_t_in_local)
1962 IF (penalty_occ_local)
THEN
1963 DO idim0 = 1, dim_op
1964 DO reim = 1,
SIZE(op_sm_set_qs, 1)
1965 DEALLOCATE (op_sm_set_qs(reim, idim0)%matrix)
1966 DEALLOCATE (op_sm_set_almo(reim, idim0)%matrix)
1969 DEALLOCATE (op_sm_set_qs)
1970 DEALLOCATE (op_sm_set_almo)
1971 DEALLOCATE (weights)
1974 IF (.NOT. converged .AND. .NOT. optimizer%early_stopping_on)
THEN
1975 cpabort(
"Optimization not converged! ")
1978 CALL timestop(handle)
1999 matrix_s, matrix_mo_in, matrix_mo_out, &
2000 template_matrix_sigma, overlap_determinant, &
2001 mat_distr_aos, virtuals, eps_filter)
2005 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:), &
2006 INTENT(INOUT) :: matrix_mo_in, matrix_mo_out
2007 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:), &
2008 INTENT(IN) :: template_matrix_sigma
2009 REAL(kind=
dp),
INTENT(INOUT) :: overlap_determinant
2010 INTEGER,
INTENT(IN) :: mat_distr_aos
2011 LOGICAL,
INTENT(IN) :: virtuals
2012 REAL(kind=
dp),
INTENT(IN) :: eps_filter
2014 CHARACTER(len=*),
PARAMETER :: routinen =
'almo_scf_construct_nlmos'
2016 CHARACTER(LEN=30) :: iter_type, print_string
2017 INTEGER :: cg_iteration, dim_op, handle, iatom, idim0, isgf, ispin, iteration, &
2018 line_search_iteration, linear_search_type, max_iter, natom, ncol, nspins, &
2019 outer_iteration, outer_max_iter, prec_type, reim, unit_nr
2020 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: first_sgf, last_sgf, nocc, nsgf
2021 LOGICAL :: converged, d_bfgs, just_started, l_bfgs, &
2022 line_search, outer_prepare_to_exit, &
2023 prepare_to_exit, reset_conjugator
2024 REAL(kind=
dp) :: appr_sec_der, beta, bfgs_rho, bfgs_sum, denom, denom2, e0, e1, g0, g0sign, &
2025 g1, g1sign, grad_norm, line_search_error, localization_obj_function, &
2026 localization_obj_function_ispin, next_step_size_guess, obj_function_ispin, objf_diff, &
2027 objf_new, objf_old, penalty_amplitude, penalty_func_ispin, penalty_func_new, spin_factor, &
2028 step_size, t1, t2, tempreal
2029 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: diagonal, grad_norm_spin, &
2030 penalty_vol_prefactor, &
2031 suggested_vol_penalty, weights
2034 TYPE(
dbcsr_p_type),
DIMENSION(:),
POINTER :: qs_matrix_s
2035 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
POINTER :: op_sm_set_almo, op_sm_set_qs
2036 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:) :: approx_inv_hessian, bfgs_s, bfgs_y, grad, &
2037 m_s0, m_sig_sqrti_ii, m_siginv, m_sigma, m_t_mo_local, m_theta, m_theta_normalized, &
2038 prev_grad, prev_m_theta, prev_minus_prec_grad, prev_step, step, tempnocc1, tempoccocc1, &
2039 tempoccocc2, tempoccocc3
2040 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:, :, :) :: m_b0
2044 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
2046 CALL timeset(routinen, handle)
2050 IF (logger%para_env%is_source())
THEN
2056 nspins =
SIZE(matrix_mo_in)
2058 IF (unit_nr > 0)
THEN
2060 IF (.NOT. virtuals)
THEN
2061 WRITE (unit_nr,
'(T2,A,A,A)') repeat(
"-", 24), &
2062 " Optimization of occupied NLMOs ", repeat(
"-", 23)
2064 WRITE (unit_nr,
'(T2,A,A,A)') repeat(
"-", 24), &
2065 " Optimization of virtual NLMOs ", repeat(
"-", 24)
2068 WRITE (unit_nr,
'(T2,A13,A6,A23,A14,A14,A9)')
"Method",
"Iter", &
2069 "Objective Function",
"Change",
"Convergence",
"Time"
2070 WRITE (unit_nr,
'(T2,A)') repeat(
"-", 79)
2073 NULLIFY (particle_set)
2076 matrix_s=qs_matrix_s, &
2078 particle_set=particle_set, &
2079 qs_kind_set=qs_kind_set)
2081 natom =
SIZE(particle_set, 1)
2082 ALLOCATE (first_sgf(natom))
2083 ALLOCATE (last_sgf(natom))
2084 ALLOCATE (nsgf(natom))
2087 first_sgf=first_sgf, last_sgf=last_sgf, nsgf=nsgf)
2091 ALLOCATE (m_theta(nspins))
2092 DO ispin = 1, nspins
2094 template=template_matrix_sigma(ispin), &
2095 matrix_type=dbcsr_type_no_symmetry)
2101 SELECT CASE (optimizer%opt_penalty%operator_type)
2104 IF (cell%orthorhombic)
THEN
2109 ALLOCATE (weights(6))
2112 ALLOCATE (op_sm_set_qs(2, dim_op))
2113 ALLOCATE (op_sm_set_almo(2, dim_op))
2115 ALLOCATE (m_b0(2, dim_op, nspins))
2116 DO idim0 = 1, dim_op
2117 DO reim = 1,
SIZE(op_sm_set_qs, 1)
2118 NULLIFY (op_sm_set_qs(reim, idim0)%matrix, op_sm_set_almo(reim, idim0)%matrix)
2119 ALLOCATE (op_sm_set_qs(reim, idim0)%matrix)
2120 ALLOCATE (op_sm_set_almo(reim, idim0)%matrix)
2121 CALL dbcsr_copy(op_sm_set_qs(reim, idim0)%matrix, qs_matrix_s(1)%matrix, &
2123 CALL dbcsr_set(op_sm_set_qs(reim, idim0)%matrix, 0.0_dp)
2124 CALL dbcsr_copy(op_sm_set_almo(reim, idim0)%matrix, matrix_s, &
2126 CALL dbcsr_set(op_sm_set_almo(reim, idim0)%matrix, 0.0_dp)
2127 DO ispin = 1, nspins
2129 template=m_theta(ispin), &
2130 matrix_type=dbcsr_type_no_symmetry)
2131 CALL dbcsr_set(m_b0(reim, idim0, ispin), 0.0_dp)
2141 ALLOCATE (weights(dim_op))
2144 ALLOCATE (m_b0(1, dim_op, nspins))
2146 DO idim0 = 1, dim_op
2148 DO ispin = 1, nspins
2150 template=m_theta(ispin), &
2151 matrix_type=dbcsr_type_no_symmetry)
2152 CALL dbcsr_set(m_b0(reim, idim0, ispin), 0.0_dp)
2160 penalty_amplitude = optimizer%opt_penalty%penalty_strength
2165 prec_type = optimizer%preconditioner
2171 IF (l_bfgs .AND. (optimizer%conjugator .NE.
cg_zero))
THEN
2172 cpabort(
"Cannot use conjugators with BFGS")
2175 CALL lbfgs_create(nlmo_lbfgs_history, nspins, nstore=10)
2178 IF (nspins == 1)
THEN
2179 spin_factor = 2.0_dp
2181 spin_factor = 1.0_dp
2184 ALLOCATE (grad_norm_spin(nspins))
2185 ALLOCATE (nocc(nspins))
2186 ALLOCATE (penalty_vol_prefactor(nspins))
2187 ALLOCATE (suggested_vol_penalty(nspins))
2193 ALLOCATE (m_t_mo_local(nspins))
2194 DO ispin = 1, nspins
2196 template=matrix_mo_in(ispin), &
2197 matrix_type=dbcsr_type_no_symmetry)
2198 CALL dbcsr_copy(m_t_mo_local(ispin), matrix_mo_in(ispin))
2201 ALLOCATE (approx_inv_hessian(nspins))
2202 ALLOCATE (m_theta_normalized(nspins))
2203 ALLOCATE (prev_m_theta(nspins))
2204 ALLOCATE (m_s0(nspins))
2205 ALLOCATE (prev_grad(nspins))
2206 ALLOCATE (grad(nspins))
2207 ALLOCATE (prev_step(nspins))
2208 ALLOCATE (step(nspins))
2209 ALLOCATE (prev_minus_prec_grad(nspins))
2210 ALLOCATE (m_sig_sqrti_ii(nspins))
2211 ALLOCATE (m_sigma(nspins))
2212 ALLOCATE (m_siginv(nspins))
2213 ALLOCATE (tempnocc1(nspins))
2214 ALLOCATE (tempoccocc1(nspins))
2215 ALLOCATE (tempoccocc2(nspins))
2216 ALLOCATE (tempoccocc3(nspins))
2217 ALLOCATE (bfgs_y(nspins))
2218 ALLOCATE (bfgs_s(nspins))
2220 DO ispin = 1, nspins
2224 template=matrix_mo_out(ispin), &
2225 matrix_type=dbcsr_type_no_symmetry)
2227 template=m_theta(ispin), &
2228 matrix_type=dbcsr_type_no_symmetry)
2230 template=m_theta(ispin), &
2231 matrix_type=dbcsr_type_no_symmetry)
2233 template=m_theta(ispin), &
2234 matrix_type=dbcsr_type_no_symmetry)
2236 template=m_theta(ispin), &
2237 matrix_type=dbcsr_type_no_symmetry)
2239 template=m_theta(ispin), &
2240 matrix_type=dbcsr_type_no_symmetry)
2242 template=m_theta(ispin), &
2243 matrix_type=dbcsr_type_no_symmetry)
2245 template=m_theta(ispin), &
2246 matrix_type=dbcsr_type_no_symmetry)
2248 template=m_theta(ispin), &
2249 matrix_type=dbcsr_type_no_symmetry)
2251 template=m_theta(ispin), &
2252 matrix_type=dbcsr_type_no_symmetry)
2254 template=m_theta(ispin), &
2255 matrix_type=dbcsr_type_no_symmetry)
2257 template=m_theta(ispin), &
2258 matrix_type=dbcsr_type_no_symmetry)
2260 template=m_theta(ispin), &
2261 matrix_type=dbcsr_type_no_symmetry)
2263 template=m_theta(ispin), &
2264 matrix_type=dbcsr_type_no_symmetry)
2266 template=m_theta(ispin), &
2267 matrix_type=dbcsr_type_no_symmetry)
2269 template=m_theta(ispin), &
2270 matrix_type=dbcsr_type_no_symmetry)
2272 template=m_theta(ispin), &
2273 matrix_type=dbcsr_type_no_symmetry)
2275 template=m_theta(ispin), &
2276 matrix_type=dbcsr_type_no_symmetry)
2279 CALL dbcsr_set(prev_step(ispin), 0.0_dp)
2282 nfullrows_total=nocc(ispin))
2284 penalty_vol_prefactor(ispin) = -penalty_amplitude
2289 m_t_mo_local(ispin), &
2290 0.0_dp, tempnocc1(ispin), &
2291 filter_eps=eps_filter)
2293 m_t_mo_local(ispin), &
2295 0.0_dp, m_s0(ispin), &
2296 filter_eps=eps_filter)
2298 SELECT CASE (optimizer%opt_penalty%operator_type)
2303 DO idim0 = 1,
SIZE(op_sm_set_qs, 2)
2305 DO reim = 1,
SIZE(op_sm_set_qs, 1)
2308 op_sm_set_almo(reim, idim0)%matrix, mat_distr_aos)
2311 op_sm_set_almo(reim, idim0)%matrix, &
2312 m_t_mo_local(ispin), &
2313 0.0_dp, tempnocc1(ispin), &
2314 filter_eps=eps_filter)
2317 m_t_mo_local(ispin), &
2319 0.0_dp, m_b0(reim, idim0, ispin), &
2320 filter_eps=eps_filter)
2322 DEALLOCATE (op_sm_set_qs(reim, idim0)%matrix)
2323 DEALLOCATE (op_sm_set_almo(reim, idim0)%matrix)
2334 isgf = first_sgf(iatom)
2339 m_t_mo_local(ispin), &
2340 0.0_dp, tempnocc1(ispin), &
2341 filter_eps=eps_filter)
2344 m_t_mo_local(ispin), &
2346 0.0_dp, m_b0(1, iatom, ispin), &
2347 first_k=isgf, last_k=isgf + ncol - 1, &
2348 filter_eps=eps_filter)
2352 m_t_mo_local(ispin), &
2353 0.0_dp, tempnocc1(ispin), &
2354 first_k=isgf, last_k=isgf + ncol - 1, &
2355 filter_eps=eps_filter)
2358 m_t_mo_local(ispin), &
2360 1.0_dp, m_b0(1, iatom, ispin), &
2361 filter_eps=eps_filter)
2369 IF (optimizer%opt_penalty%operator_type .EQ.
op_loc_berry)
THEN
2370 DO idim0 = 1,
SIZE(op_sm_set_qs, 2)
2371 DO reim = 1,
SIZE(op_sm_set_qs, 1)
2372 DEALLOCATE (op_sm_set_qs(reim, idim0)%matrix)
2373 DEALLOCATE (op_sm_set_almo(reim, idim0)%matrix)
2376 DEALLOCATE (op_sm_set_qs, op_sm_set_almo)
2380 outer_max_iter = optimizer%max_iter_outer_loop
2381 outer_prepare_to_exit = .false.
2384 penalty_func_new = 0.0_dp
2385 linear_search_type = 1
2386 localization_obj_function = 0.0_dp
2387 penalty_func_new = 0.0_dp
2392 max_iter = optimizer%max_iter
2393 prepare_to_exit = .false.
2394 line_search = .false.
2398 line_search_iteration = 0
2399 obj_function_ispin = 0.0_dp
2403 line_search_error = 0.0_dp
2405 next_step_size_guess = 0.0_dp
2409 just_started = (iteration .EQ. 0) .AND. (outer_iteration .EQ. 0)
2411 DO ispin = 1, nspins
2417 m_s0(ispin), m_theta(ispin), 0.0_dp, &
2418 tempoccocc1(ispin), &
2419 filter_eps=eps_filter)
2420 CALL dbcsr_set(m_sig_sqrti_ii(ispin), 0.0_dp)
2423 m_theta(ispin), tempoccocc1(ispin), 0.0_dp, &
2424 m_sig_sqrti_ii(ispin), &
2425 retain_sparsity=.true.)
2426 ALLOCATE (diagonal(nocc(ispin)))
2428 CALL group%sum(diagonal)
2430 diagonal(:) = 1.0_dp/sqrt(diagonal(:))
2431 CALL dbcsr_set(m_sig_sqrti_ii(ispin), 0.0_dp)
2433 DEALLOCATE (diagonal)
2437 m_sig_sqrti_ii(ispin), &
2438 0.0_dp, m_theta_normalized(ispin), &
2439 filter_eps=eps_filter)
2443 m_t_mo_local(ispin), &
2444 m_theta_normalized(ispin), &
2445 0.0_dp, matrix_mo_out(ispin), &
2446 filter_eps=eps_filter)
2451 localization_obj_function = 0.0_dp
2452 penalty_func_new = 0.0_dp
2453 DO ispin = 1, nspins
2455 CALL compute_obj_nlmos( &
2457 localization_obj_function_ispin=localization_obj_function_ispin, &
2458 penalty_func_ispin=penalty_func_ispin, &
2459 overlap_determinant=overlap_determinant, &
2460 m_sigma=m_sigma(ispin), &
2462 m_b0=m_b0(:, :, ispin), &
2463 m_theta_normalized=m_theta_normalized(ispin), &
2464 template_matrix_mo=matrix_mo_out(ispin), &
2467 just_started=just_started, &
2468 penalty_vol_prefactor=penalty_vol_prefactor(ispin), &
2469 penalty_amplitude=penalty_amplitude, &
2470 eps_filter=eps_filter)
2472 localization_obj_function = localization_obj_function + localization_obj_function_ispin
2473 penalty_func_new = penalty_func_new + penalty_func_ispin
2476 objf_new = penalty_func_new + localization_obj_function
2478 DO ispin = 1, nspins
2482 IF (line_search_iteration .EQ. 0 .AND. iteration .NE. 0)
THEN
2483 CALL dbcsr_copy(prev_grad(ispin), grad(ispin))
2489 DO ispin = 1, nspins
2492 matrix_inverse=m_siginv(ispin), &
2493 matrix=m_sigma(ispin), &
2494 threshold=eps_filter*10.0_dp, &
2495 filter_eps=eps_filter, &
2498 CALL compute_gradient_nlmos( &
2499 m_grad_out=grad(ispin), &
2500 m_b0=m_b0(:, :, ispin), &
2503 m_theta_normalized=m_theta_normalized(ispin), &
2504 m_siginv=m_siginv(ispin), &
2505 m_sig_sqrti_ii=m_sig_sqrti_ii(ispin), &
2506 penalty_vol_prefactor=penalty_vol_prefactor(ispin), &
2507 eps_filter=eps_filter, &
2508 suggested_vol_penalty=suggested_vol_penalty(ispin))
2513 DO ispin = 1, nspins
2516 grad_norm = maxval(grad_norm_spin)
2518 converged = (grad_norm .LE. optimizer%eps_error)
2519 IF (converged .OR. (iteration .GE. max_iter))
THEN
2520 prepare_to_exit = .true.
2524 IF (.NOT. prepare_to_exit)
THEN
2529 IF (iteration .NE. 0)
THEN
2533 IF (.NOT. line_search)
THEN
2535 line_search = .true.
2536 line_search_iteration = line_search_iteration + 1
2542 line_search_error = 0.0_dp
2546 DO ispin = 1, nspins
2548 CALL dbcsr_dot(grad(ispin), step(ispin), tempreal)
2549 line_search_error = line_search_error + tempreal
2550 CALL dbcsr_dot(grad(ispin), grad(ispin), tempreal)
2551 denom = denom + tempreal
2552 CALL dbcsr_dot(step(ispin), step(ispin), tempreal)
2553 denom2 = denom2 + tempreal
2559 line_search_error = line_search_error/sqrt(denom)/sqrt(denom2)
2561 IF (abs(line_search_error) .GT. optimizer%lin_search_eps_error)
THEN
2562 line_search = .true.
2563 line_search_iteration = line_search_iteration + 1
2565 line_search = .false.
2566 line_search_iteration = 0
2573 IF (line_search)
THEN
2576 objf_diff = objf_new - objf_old
2581 IF (.NOT. line_search)
THEN
2583 cg_iteration = cg_iteration + 1
2586 DO ispin = 1, nspins
2587 CALL dbcsr_copy(prev_step(ispin), step(ispin))
2595 DO ispin = 1, nspins
2605 IF (iteration .EQ. 0)
THEN
2611 IF (nspins .GT. 1)
THEN
2612 DO ispin = 2, nspins
2613 CALL dbcsr_copy(approx_inv_hessian(ispin), approx_inv_hessian(1))
2617 ELSE IF (l_bfgs)
THEN
2619 CALL lbfgs_seed(nlmo_lbfgs_history, m_theta, grad)
2620 DO ispin = 1, nspins
2628 DO ispin = 1, nspins
2666 DO ispin = 1, nspins
2670 CALL dbcsr_add(bfgs_y(ispin), prev_grad(ispin), 1.0_dp, -1.0_dp)
2671 CALL dbcsr_copy(bfgs_s(ispin), m_theta(ispin))
2672 CALL dbcsr_add(bfgs_s(ispin), prev_m_theta(ispin), 1.0_dp, -1.0_dp)
2675 CALL dbcsr_dot(grad(ispin), step(ispin), bfgs_rho)
2676 bfgs_rho = 1.0_dp/bfgs_rho
2679 CALL dbcsr_dot(bfgs_y(ispin), bfgs_y(ispin), bfgs_sum)
2682 CALL dbcsr_copy(tempoccocc2(ispin), approx_inv_hessian(ispin))
2686 CALL dbcsr_add(tempoccocc2(ispin), tempoccocc1(ispin), 1.0_dp, bfgs_rho)
2690 approx_inv_hessian(ispin), tempoccocc3(ispin))
2691 CALL dbcsr_add(tempoccocc2(ispin), tempoccocc3(ispin), &
2692 1.0_dp, bfgs_rho*bfgs_rho*bfgs_sum)
2696 approx_inv_hessian(ispin), tempoccocc1(ispin))
2698 CALL dbcsr_add(tempoccocc2(ispin), tempoccocc3(ispin), &
2699 1.0_dp, -2.0_dp*bfgs_rho)
2701 CALL dbcsr_copy(approx_inv_hessian(ispin), tempoccocc2(ispin))
2705 ELSE IF (l_bfgs)
THEN
2713 IF (.NOT. l_bfgs)
THEN
2715 DO ispin = 1, nspins
2718 grad(ispin), step(ispin))
2728 IF (iteration .EQ. 0)
THEN
2729 reset_conjugator = .true.
2733 IF (.NOT. reset_conjugator)
THEN
2734 CALL compute_cg_beta( &
2736 reset_conjugator=reset_conjugator, &
2737 conjugator=optimizer%conjugator, &
2739 prev_grad=prev_grad(:), &
2741 prev_step=prev_step(:), &
2742 prev_minus_prec_grad=prev_minus_prec_grad(:) &
2747 IF (reset_conjugator)
THEN
2750 IF (unit_nr > 0 .AND. (.NOT. just_started))
THEN
2751 WRITE (unit_nr,
'(T2,A35)')
"Re-setting conjugator to zero"
2753 reset_conjugator = .false.
2758 DO ispin = 1, nspins
2760 CALL dbcsr_copy(prev_minus_prec_grad(ispin), step(ispin))
2763 CALL dbcsr_add(step(ispin), prev_step(ispin), 1.0_dp, beta)
2770 IF (.NOT. line_search)
THEN
2776 DO ispin = 1, nspins
2777 CALL dbcsr_dot(grad(ispin), step(ispin), tempreal)
2780 g0sign = sign(1.0_dp, g0)
2781 IF (linear_search_type .EQ. 1)
THEN
2782 IF (iteration .EQ. 0)
THEN
2783 step_size = optimizer%lin_search_step_size_guess
2785 IF (next_step_size_guess .LE. 0.0_dp)
THEN
2786 step_size = optimizer%lin_search_step_size_guess
2789 step_size = optimizer%lin_search_step_size_guess
2793 ELSE IF (linear_search_type .EQ. 2)
THEN
2796 step_size = optimizer%lin_search_step_size_guess
2798 IF (unit_nr > 0)
THEN
2799 WRITE (unit_nr,
'(T21,3A19)')
"Line position",
"Line grad",
"Next line step"
2800 WRITE (unit_nr,
'(T2,A19,3F19.5)')
"Line search", 0.0_dp, g0, step_size
2802 next_step_size_guess = step_size
2806 DO ispin = 1, nspins
2807 CALL dbcsr_dot(grad(ispin), step(ispin), tempreal)
2810 g1sign = sign(1.0_dp, g1)
2811 IF (linear_search_type .EQ. 1)
THEN
2814 appr_sec_der = (g1 - g0)/step_size
2819 step_size = -g1/appr_sec_der
2820 ELSE IF (linear_search_type .EQ. 2)
THEN
2823 IF (g1sign .NE. g0sign)
THEN
2824 step_size = -step_size/2.0;
2826 step_size = step_size*1.5;
2830 IF (unit_nr > 0)
THEN
2831 WRITE (unit_nr,
'(T21,3A19)')
"Line position",
"Line grad",
"Next line step"
2832 WRITE (unit_nr,
'(T2,A19,3F19.5)')
"Line search", next_step_size_guess, g1, step_size
2837 next_step_size_guess = next_step_size_guess + step_size
2841 DO ispin = 1, nspins
2842 IF (.NOT. line_search)
THEN
2844 CALL dbcsr_copy(prev_m_theta(ispin), m_theta(ispin))
2846 CALL dbcsr_add(m_theta(ispin), step(ispin), 1.0_dp, step_size)
2851 IF (line_search)
THEN
2858 IF (unit_nr > 0)
THEN
2859 iter_type = trim(
"NLMO OPT "//iter_type)
2860 WRITE (unit_nr,
'(T2,A13,I6,F23.10,E14.5,F14.9,F9.2)') &
2861 iter_type, iteration, &
2862 objf_new, objf_diff, grad_norm, &
2864 WRITE (unit_nr,
'(T2,A19,F23.10)') &
2865 "Localization:", localization_obj_function
2866 WRITE (unit_nr,
'(T2,A19,F23.10)') &
2867 "Orthogonalization:", penalty_func_new
2871 iteration = iteration + 1
2872 IF (prepare_to_exit)
EXIT
2876 IF (converged .OR. (outer_iteration .GE. outer_max_iter))
THEN
2877 outer_prepare_to_exit = .true.
2880 outer_iteration = outer_iteration + 1
2881 IF (outer_prepare_to_exit)
EXIT
2886 optimizer%opt_penalty%penalty_strength = 0.0_dp
2887 DO ispin = 1, nspins
2888 optimizer%opt_penalty%penalty_strength = optimizer%opt_penalty%penalty_strength + &
2889 (-1.0_dp)*penalty_vol_prefactor(ispin)
2891 optimizer%opt_penalty%penalty_strength = optimizer%opt_penalty%penalty_strength/nspins
2896 iter_type =
"Unconverged"
2899 IF (unit_nr > 0)
THEN
2900 WRITE (unit_nr,
'()')
2901 print_string = trim(iter_type)//
" localization:"
2902 WRITE (unit_nr,
'(T2,A29,F30.10)') &
2903 print_string, localization_obj_function
2904 print_string = trim(iter_type)//
" determinant:"
2905 WRITE (unit_nr,
'(T2,A29,F30.10)') &
2906 print_string, overlap_determinant
2907 print_string = trim(iter_type)//
" penalty strength:"
2908 WRITE (unit_nr,
'(T2,A29,F30.10)') &
2909 print_string, optimizer%opt_penalty%penalty_strength
2916 DO ispin = 1, nspins
2917 DO idim0 = 1,
SIZE(m_b0, 2)
2918 DO reim = 1,
SIZE(m_b0, 1)
2944 DEALLOCATE (grad_norm_spin)
2946 DEALLOCATE (penalty_vol_prefactor)
2947 DEALLOCATE (suggested_vol_penalty)
2949 DEALLOCATE (approx_inv_hessian)
2950 DEALLOCATE (prev_m_theta)
2951 DEALLOCATE (m_theta_normalized)
2953 DEALLOCATE (prev_grad)
2955 DEALLOCATE (prev_step)
2957 DEALLOCATE (prev_minus_prec_grad)
2958 DEALLOCATE (m_sig_sqrti_ii)
2959 DEALLOCATE (m_sigma)
2960 DEALLOCATE (m_siginv)
2961 DEALLOCATE (tempnocc1)
2962 DEALLOCATE (tempoccocc1)
2963 DEALLOCATE (tempoccocc2)
2964 DEALLOCATE (tempoccocc3)
2968 DEALLOCATE (m_theta, m_t_mo_local)
2970 DEALLOCATE (weights)
2971 DEALLOCATE (first_sgf, last_sgf, nsgf)
2973 IF (.NOT. converged)
THEN
2974 cpabort(
"Optimization not converged! ")
2977 CALL timestop(handle)
2999 SUBROUTINE xalmo_analysis(detailed_analysis, eps_filter, m_T_in, m_T0_in, &
3000 m_siginv_in, m_siginv0_in, m_S_in, m_KS0_in, m_quench_t_in, energy_out, &
3001 m_eda_out, m_cta_out)
3003 LOGICAL,
INTENT(IN) :: detailed_analysis
3004 REAL(kind=
dp),
INTENT(IN) :: eps_filter
3005 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(IN) :: m_t_in, m_t0_in, m_siginv_in, &
3006 m_siginv0_in, m_s_in, m_ks0_in, &
3008 REAL(kind=
dp),
INTENT(INOUT) :: energy_out
3009 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(INOUT) :: m_eda_out, m_cta_out
3011 CHARACTER(len=*),
PARAMETER :: routinen =
'xalmo_analysis'
3013 INTEGER :: handle, ispin, nspins
3014 REAL(kind=
dp) :: energy_ispin, spin_factor
3015 TYPE(
dbcsr_type) :: ftsiginv0, fvo0, m_x, siginvtftsiginv0, &
3018 CALL timeset(routinen, handle)
3020 nspins =
SIZE(m_t_in)
3022 IF (nspins == 1)
THEN
3023 spin_factor = 2.0_dp
3025 spin_factor = 1.0_dp
3029 DO ispin = 1, nspins
3033 template=m_t_in(ispin), &
3034 matrix_type=dbcsr_type_no_symmetry)
3036 template=m_t_in(ispin), &
3037 matrix_type=dbcsr_type_no_symmetry)
3039 template=m_t_in(ispin), &
3040 matrix_type=dbcsr_type_no_symmetry)
3042 template=m_t_in(ispin), &
3043 matrix_type=dbcsr_type_no_symmetry)
3045 template=m_siginv0_in(ispin), &
3046 matrix_type=dbcsr_type_no_symmetry)
3049 CALL compute_frequently_used_matrices( &
3050 filter_eps=eps_filter, &
3051 m_t_in=m_t0_in(ispin), &
3052 m_siginv_in=m_siginv0_in(ispin), &
3054 m_f_in=m_ks0_in(ispin), &
3055 m_ftsiginv_out=ftsiginv0, &
3056 m_siginvtftsiginv_out=siginvtftsiginv0, &
3059 CALL dbcsr_copy(fvo0, ftsiginv0, keep_sparsity=.true.)
3064 retain_sparsity=.true.)
3068 CALL dbcsr_add(m_x, m_t_in(ispin), -1.0_dp, 1.0_dp)
3071 energy_out = energy_out + energy_ispin*spin_factor
3073 IF (detailed_analysis)
THEN
3083 m_siginv0_in(ispin), &
3084 0.0_dp, ftsiginv0, &
3085 filter_eps=eps_filter)
3091 filter_eps=eps_filter)
3096 0.0_dp, siginvtftsiginv0, &
3097 filter_eps=eps_filter)
3104 filter_eps=eps_filter)
3108 m_siginv_in(ispin), &
3109 0.0_dp, ftsiginv0, &
3110 filter_eps=eps_filter)
3113 ftsiginv0, m_cta_out(ispin))
3127 CALL timestop(handle)
3129 END SUBROUTINE xalmo_analysis
3146 SUBROUTINE compute_frequently_used_matrices(filter_eps, &
3147 m_T_in, m_siginv_in, m_S_in, m_F_in, m_FTsiginv_out, &
3148 m_siginvTFTsiginv_out, m_ST_out)
3150 REAL(kind=
dp),
INTENT(IN) :: filter_eps
3151 TYPE(
dbcsr_type),
INTENT(IN) :: m_t_in, m_siginv_in, m_s_in, m_f_in
3152 TYPE(
dbcsr_type),
INTENT(INOUT) :: m_ftsiginv_out, m_siginvtftsiginv_out, &
3155 CHARACTER(len=*),
PARAMETER :: routinen =
'compute_frequently_used_matrices'
3160 CALL timeset(routinen, handle)
3164 matrix_type=dbcsr_type_no_symmetry)
3166 template=m_siginv_in, &
3167 matrix_type=dbcsr_type_no_symmetry)
3172 0.0_dp, m_tmp_no_1, &
3173 filter_eps=filter_eps)
3178 0.0_dp, m_ftsiginv_out, &
3179 filter_eps=filter_eps)
3184 0.0_dp, m_tmp_oo_1, &
3185 filter_eps=filter_eps)
3190 0.0_dp, m_siginvtftsiginv_out, &
3191 filter_eps=filter_eps)
3197 filter_eps=filter_eps)
3202 CALL timestop(handle)
3204 END SUBROUTINE compute_frequently_used_matrices
3214 SUBROUTINE split_v_blk(almo_scf_env)
3218 CHARACTER(len=*),
PARAMETER :: routinen =
'split_v_blk'
3220 INTEGER :: discarded_v, handle, iblock_col, &
3221 iblock_col_size, iblock_row, &
3222 iblock_row_size, ispin, retained_v
3223 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: data_p
3226 CALL timeset(routinen, handle)
3228 DO ispin = 1, almo_scf_env%nspins
3231 work_mutable=.true.)
3233 work_mutable=.true.)
3240 row_size=iblock_row_size, col_size=iblock_col_size)
3242 IF (iblock_row .NE. iblock_col)
THEN
3243 cpabort(
"off-diagonal block found")
3246 retained_v = almo_scf_env%nvirt_of_domain(iblock_col, ispin)
3247 discarded_v = almo_scf_env%nvirt_disc_of_domain(iblock_col, ispin)
3248 cpassert(retained_v .GT. 0)
3249 cpassert(discarded_v .GT. 0)
3250 CALL dbcsr_put_block(almo_scf_env%matrix_v_disc_blk(ispin), iblock_row, iblock_col, &
3251 block=data_p(:, (retained_v + 1):iblock_col_size))
3252 CALL dbcsr_put_block(almo_scf_env%matrix_v_blk(ispin), iblock_row, iblock_col, &
3253 block=data_p(:, 1:retained_v))
3263 CALL timestop(handle)
3265 END SUBROUTINE split_v_blk
3274 SUBROUTINE harris_foulkes_correction(almo_scf_env)
3278 CHARACTER(len=*),
PARAMETER :: routinen =
'harris_foulkes_correction'
3279 INTEGER,
PARAMETER :: cayley_transform = 1, dm_ls_step = 2
3281 INTEGER :: algorithm_id, handle, handle1, handle2, handle3, handle4, handle5, handle6, &
3282 handle7, handle8, ispin, iteration, n, nmins, nspin, opt_k_max_iter, &
3283 outer_opt_k_iteration, outer_opt_k_max_iter, unit_nr
3284 INTEGER,
DIMENSION(1) :: fake, nelectron_spin_real
3285 LOGICAL :: converged, line_search, md_in_k_space, outer_opt_k_prepare_to_exit, &
3286 prepare_to_exit, reset_conjugator, reset_step_size, use_cubic_approximation, &
3287 use_quadratic_approximation
3288 REAL(kind=
dp) :: aa, bb, beta, conjugacy_error, conjugacy_error_threshold, &
3289 delta_obj_function, denom, energy_correction_final, frob_matrix, frob_matrix_base, fun0, &
3290 fun1, gfun0, gfun1, grad_norm, grad_norm_frob, kappa, kin_energy, line_search_error, &
3291 line_search_error_threshold, num_threshold, numer, obj_function, quadratic_approx_error, &
3292 quadratic_approx_error_threshold, safety_multiplier, spin_factor, step_size, &
3293 step_size_quadratic_approx, step_size_quadratic_approx2, t1, t1a, t1cholesky, t2, t2a, &
3294 t2cholesky, tau, time_step, x_opt_eps_adaptive, x_opt_eps_adaptive_factor
3295 REAL(kind=
dp),
DIMENSION(1) :: local_mu
3296 REAL(kind=
dp),
DIMENSION(2) :: energy_correction
3297 REAL(kind=
dp),
DIMENSION(3) :: minima
3300 TYPE(
dbcsr_type) :: grad, k_vd_index_down, k_vr_index_down, matrix_k_central, matrix_tmp1, &
3301 matrix_tmp2, prec, prev_grad, prev_minus_prec_grad, prev_step, sigma_oo_curr, &
3302 sigma_oo_curr_inv, sigma_vv_sqrt, sigma_vv_sqrt_guess, sigma_vv_sqrt_inv, &
3303 sigma_vv_sqrt_inv_guess, step, t_curr, tmp1_n_vr, tmp2_n_o, tmp3_vd_vr, tmp4_o_vr, &
3304 tmp_k_blk, vd_fixed, vd_index_sqrt, vd_index_sqrt_inv, velocity, vr_fixed, vr_index_sqrt, &
3306 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:) :: matrix_p_almo_scf_converged
3308 CALL timeset(routinen, handle)
3312 IF (logger%para_env%is_source())
THEN
3318 nspin = almo_scf_env%nspins
3319 energy_correction_final = 0.0_dp
3320 IF (nspin .EQ. 1)
THEN
3321 spin_factor = 2.0_dp
3323 spin_factor = 1.0_dp
3326 IF (almo_scf_env%deloc_use_occ_orbs)
THEN
3327 algorithm_id = cayley_transform
3329 algorithm_id = dm_ls_step
3334 SELECT CASE (algorithm_id)
3335 CASE (cayley_transform)
3339 IF (almo_scf_env%nspins == 1)
THEN
3340 CALL dbcsr_scale(almo_scf_env%matrix_p(1), 1.0_dp/spin_factor)
3346 CALL dbcsr_copy(almo_scf_env%matrix_t(ispin), &
3347 almo_scf_env%matrix_t_blk(ispin))
3354 IF (unit_nr > 0)
THEN
3355 WRITE (unit_nr, *)
"sqrt and inv(sqrt) of MO overlap matrix"
3357 CALL dbcsr_create(almo_scf_env%matrix_sigma_sqrt(ispin), &
3358 template=almo_scf_env%matrix_sigma(ispin), &
3359 matrix_type=dbcsr_type_no_symmetry)
3360 CALL dbcsr_create(almo_scf_env%matrix_sigma_sqrt_inv(ispin), &
3361 template=almo_scf_env%matrix_sigma(ispin), &
3362 matrix_type=dbcsr_type_no_symmetry)
3365 almo_scf_env%matrix_sigma_sqrt_inv(ispin), &
3366 almo_scf_env%matrix_sigma(ispin), &
3367 threshold=almo_scf_env%eps_filter, &
3368 order=almo_scf_env%order_lanczos, &
3369 eps_lanczos=almo_scf_env%eps_lanczos, &
3370 max_iter_lanczos=almo_scf_env%max_iter_lanczos)
3373 CALL dbcsr_create(matrix_tmp1, template=almo_scf_env%matrix_sigma(ispin), &
3374 matrix_type=dbcsr_type_no_symmetry)
3375 CALL dbcsr_create(matrix_tmp2, template=almo_scf_env%matrix_sigma(ispin), &
3376 matrix_type=dbcsr_type_no_symmetry)
3378 CALL dbcsr_multiply(
"N",
"N", 1.0_dp, almo_scf_env%matrix_sigma_sqrt_inv(ispin), &
3379 almo_scf_env%matrix_sigma(ispin), &
3380 0.0_dp, matrix_tmp1, filter_eps=almo_scf_env%eps_filter)
3382 almo_scf_env%matrix_sigma_sqrt_inv(ispin), &
3383 0.0_dp, matrix_tmp2, filter_eps=almo_scf_env%eps_filter)
3388 IF (unit_nr > 0)
THEN
3389 WRITE (unit_nr, *)
"Error for (inv(sqrt(SIG))*SIG*inv(sqrt(SIG))-I)", frob_matrix/frob_matrix_base
3397 IF (almo_scf_env%almo_update_algorithm .EQ.
almo_scf_diag)
THEN
3403 line_search_error_threshold = almo_scf_env%real01
3404 conjugacy_error_threshold = almo_scf_env%real02
3405 quadratic_approx_error_threshold = almo_scf_env%real03
3406 x_opt_eps_adaptive_factor = almo_scf_env%real04
3409 outer_opt_k_max_iter = almo_scf_env%opt_k_outer_max_iter
3410 outer_opt_k_prepare_to_exit = .false.
3411 outer_opt_k_iteration = 0
3413 grad_norm_frob = 0.0_dp
3414 CALL dbcsr_set(almo_scf_env%matrix_x(ispin), 0.0_dp)
3415 IF (almo_scf_env%deloc_truncate_virt .EQ.
virt_full) outer_opt_k_max_iter = 0
3421 psi_out=almo_scf_env%matrix_v(ispin), &
3422 psi_projector=almo_scf_env%matrix_t_blk(ispin), &
3423 metric=almo_scf_env%matrix_s(1), &
3424 project_out=.true., &
3425 psi_projector_orthogonal=.false., &
3426 proj_in_template=almo_scf_env%matrix_ov(ispin), &
3427 eps_filter=almo_scf_env%eps_filter, &
3428 sig_inv_projector=almo_scf_env%matrix_sigma_inv(ispin))
3433 template=almo_scf_env%matrix_v(ispin))
3434 CALL dbcsr_copy(vr_fixed, almo_scf_env%matrix_v(ispin))
3438 template=almo_scf_env%matrix_sigma_vv(ispin), &
3439 matrix_type=dbcsr_type_no_symmetry)
3441 template=almo_scf_env%matrix_sigma_vv(ispin), &
3442 matrix_type=dbcsr_type_no_symmetry)
3444 template=almo_scf_env%matrix_sigma_vv(ispin), &
3445 matrix_type=dbcsr_type_no_symmetry)
3447 template=almo_scf_env%matrix_sigma_vv(ispin), &
3448 matrix_type=dbcsr_type_no_symmetry)
3449 CALL dbcsr_set(sigma_vv_sqrt_guess, 0.0_dp)
3451 CALL dbcsr_filter(sigma_vv_sqrt_guess, almo_scf_env%eps_filter)
3452 CALL dbcsr_set(sigma_vv_sqrt_inv_guess, 0.0_dp)
3454 CALL dbcsr_filter(sigma_vv_sqrt_inv_guess, almo_scf_env%eps_filter)
3457 IF (almo_scf_env%deloc_truncate_virt .NE.
virt_full)
THEN
3476 psi_out=almo_scf_env%matrix_v_disc(ispin), &
3477 psi_projector=almo_scf_env%matrix_t_blk(ispin), &
3478 metric=almo_scf_env%matrix_s(1), &
3479 project_out=.true., &
3480 psi_projector_orthogonal=.false., &
3481 proj_in_template=almo_scf_env%matrix_ov_disc(ispin), &
3482 eps_filter=almo_scf_env%eps_filter, &
3483 sig_inv_projector=almo_scf_env%matrix_sigma_inv(ispin))
3488 template=almo_scf_env%matrix_v_disc(ispin))
3489 CALL dbcsr_copy(vd_fixed, almo_scf_env%matrix_v_disc(ispin))
3493 template=almo_scf_env%matrix_sigma_vv_blk(ispin), &
3494 matrix_type=dbcsr_type_no_symmetry)
3507 template=almo_scf_env%matrix_vv_disc_blk(ispin), &
3508 matrix_type=dbcsr_type_no_symmetry)
3551 template=almo_scf_env%matrix_k_blk(ispin))
3552 CALL dbcsr_copy(grad, almo_scf_env%matrix_k_blk(ispin))
3555 md_in_k_space = almo_scf_env%logical01
3556 IF (md_in_k_space)
THEN
3558 template=almo_scf_env%matrix_k_blk(ispin))
3559 CALL dbcsr_copy(velocity, almo_scf_env%matrix_k_blk(ispin))
3561 time_step = almo_scf_env%opt_k_trial_step_size
3565 template=almo_scf_env%matrix_k_blk(ispin))
3568 template=almo_scf_env%matrix_k_blk(ispin))
3572 template=almo_scf_env%matrix_k_blk(ispin))
3573 CALL dbcsr_copy(prec, almo_scf_env%matrix_k_blk(ispin))
3577 CALL dbcsr_set(almo_scf_env%matrix_k_blk(ispin), 0.0_dp)
3581 template=almo_scf_env%matrix_k_blk(ispin))
3583 almo_scf_env%matrix_k_blk(ispin))
3585 template=almo_scf_env%matrix_k_blk(ispin))
3587 template=almo_scf_env%matrix_k_blk(ispin))
3590 template=almo_scf_env%matrix_t(ispin))
3592 template=almo_scf_env%matrix_sigma(ispin), &
3593 matrix_type=dbcsr_type_no_symmetry)
3595 template=almo_scf_env%matrix_sigma(ispin), &
3596 matrix_type=dbcsr_type_no_symmetry)
3598 template=almo_scf_env%matrix_v(ispin))
3600 template=almo_scf_env%matrix_k_blk(ispin))
3602 template=almo_scf_env%matrix_t(ispin))
3604 template=almo_scf_env%matrix_ov(ispin))
3606 template=almo_scf_env%matrix_k_blk(ispin))
3621 opt_k_max_iter = almo_scf_env%opt_k_max_iter
3624 prepare_to_exit = .false.
3626 line_search = .false.
3627 obj_function = 0.0_dp
3628 conjugacy_error = 0.0_dp
3629 line_search_error = 0.0_dp
3634 step_size_quadratic_approx = 0.0_dp
3635 reset_step_size = .true.
3636 IF (almo_scf_env%deloc_truncate_virt .EQ.
virt_full) opt_k_max_iter = 0
3641 CALL timeset(
'k_opt_vr', handle1)
3643 IF (almo_scf_env%deloc_truncate_virt .NE.
virt_full)
THEN
3647 almo_scf_env%matrix_k_blk(ispin), &
3648 0.0_dp, almo_scf_env%matrix_v(ispin), &
3649 filter_eps=almo_scf_env%eps_filter)
3650 CALL dbcsr_add(almo_scf_env%matrix_v(ispin), vr_fixed, &
3658 CALL get_overlap(bra=almo_scf_env%matrix_v(ispin), &
3659 ket=almo_scf_env%matrix_v(ispin), &
3660 overlap=almo_scf_env%matrix_sigma_vv(ispin), &
3661 metric=almo_scf_env%matrix_s(1), &
3662 retain_overlap_sparsity=.false., &
3663 eps_filter=almo_scf_env%eps_filter)
3666 IF (almo_scf_env%deloc_truncate_virt .EQ.
virt_full)
THEN
3667 CALL timeset(
'cholesky', handle2)
3673 template=almo_scf_env%matrix_sigma_vv(ispin), &
3674 matrix_type=dbcsr_type_no_symmetry)
3678 para_env=almo_scf_env%para_env, &
3679 blacs_env=almo_scf_env%blacs_env)
3680 CALL make_triu(sigma_vv_sqrt)
3681 CALL dbcsr_filter(sigma_vv_sqrt, almo_scf_env%eps_filter)
3684 CALL dbcsr_create(matrix_tmp1, template=almo_scf_env%matrix_sigma_vv(ispin), &
3685 matrix_type=dbcsr_type_no_symmetry)
3689 sigma_vv_sqrt_inv, op=
"SOLVE", pos=
"RIGHT", &
3690 para_env=almo_scf_env%para_env, &
3691 blacs_env=almo_scf_env%blacs_env)
3692 CALL dbcsr_filter(sigma_vv_sqrt_inv, almo_scf_env%eps_filter)
3695 CALL dbcsr_create(matrix_tmp1, template=almo_scf_env%matrix_sigma_vv(ispin), &
3696 matrix_type=dbcsr_type_no_symmetry)
3701 -1.0_dp, matrix_tmp1, filter_eps=almo_scf_env%eps_filter)
3705 IF (unit_nr > 0)
THEN
3706 WRITE (unit_nr, *)
"Error for ( U^T * U - Sig )", &
3707 frob_matrix/frob_matrix_base
3711 0.0_dp, matrix_tmp1, filter_eps=almo_scf_env%eps_filter)
3715 IF (unit_nr > 0)
THEN
3716 WRITE (unit_nr, *)
"Error for ( inv(U) * U - I )", &
3717 frob_matrix/frob_matrix_base
3722 IF (unit_nr > 0)
THEN
3723 WRITE (unit_nr, *)
"Cholesky+inverse wall-time: ", t2cholesky - t1cholesky
3725 CALL timestop(handle2)
3728 sigma_vv_sqrt_inv, &
3729 almo_scf_env%matrix_sigma_vv(ispin), &
3732 threshold=almo_scf_env%eps_filter, &
3733 order=almo_scf_env%order_lanczos, &
3734 eps_lanczos=almo_scf_env%eps_lanczos, &
3735 max_iter_lanczos=almo_scf_env%max_iter_lanczos)
3736 CALL dbcsr_copy(sigma_vv_sqrt_inv_guess, sigma_vv_sqrt_inv)
3737 CALL dbcsr_copy(sigma_vv_sqrt_guess, sigma_vv_sqrt)
3739 CALL dbcsr_create(matrix_tmp1, template=almo_scf_env%matrix_sigma_vv(ispin), &
3740 matrix_type=dbcsr_type_no_symmetry)
3741 CALL dbcsr_create(matrix_tmp2, template=almo_scf_env%matrix_sigma_vv(ispin), &
3742 matrix_type=dbcsr_type_no_symmetry)
3745 almo_scf_env%matrix_sigma_vv(ispin), &
3746 0.0_dp, matrix_tmp1, filter_eps=almo_scf_env%eps_filter)
3748 sigma_vv_sqrt_inv, &
3749 0.0_dp, matrix_tmp2, filter_eps=almo_scf_env%eps_filter)
3754 IF (unit_nr > 0)
THEN
3755 WRITE (unit_nr, *)
"Error for (inv(sqrt(SIGVV))*SIGVV*inv(sqrt(SIGVV))-I)", &
3756 frob_matrix/frob_matrix_base
3763 CALL timestop(handle1)
3767 IF ((iteration .EQ. 0) .AND. (.NOT. line_search) .AND. &
3768 (outer_opt_k_iteration .EQ. 0))
THEN
3769 x_opt_eps_adaptive = &
3770 almo_scf_env%deloc_cayley_eps_convergence
3772 x_opt_eps_adaptive = &
3773 max(abs(almo_scf_env%deloc_cayley_eps_convergence), &
3774 abs(x_opt_eps_adaptive_factor*grad_norm))
3778 para_env=almo_scf_env%para_env, &
3779 blacs_env=almo_scf_env%blacs_env, &
3780 use_occ_orbs=.true., &
3781 use_virt_orbs=.true., &
3782 occ_orbs_orthogonal=.false., &
3783 virt_orbs_orthogonal=.false., &
3784 pp_preconditioner_full=almo_scf_env%deloc_cayley_occ_precond, &
3785 qq_preconditioner_full=almo_scf_env%deloc_cayley_vir_precond, &
3786 tensor_type=almo_scf_env%deloc_cayley_tensor_type, &
3787 neglect_quadratic_term=almo_scf_env%deloc_cayley_linear, &
3788 conjugator=almo_scf_env%deloc_cayley_conjugator, &
3789 max_iter=almo_scf_env%deloc_cayley_max_iter, &
3790 calculate_energy_corr=.true., &
3793 eps_convergence=x_opt_eps_adaptive, &
3794 eps_filter=almo_scf_env%eps_filter, &
3796 q_index_up=sigma_vv_sqrt_inv, &
3797 q_index_down=sigma_vv_sqrt, &
3798 p_index_up=almo_scf_env%matrix_sigma_sqrt_inv(ispin), &
3799 p_index_down=almo_scf_env%matrix_sigma_sqrt(ispin), &
3800 matrix_ks=almo_scf_env%matrix_ks_0deloc(ispin), &
3801 matrix_t=almo_scf_env%matrix_t(ispin), &
3802 matrix_qp_template=almo_scf_env%matrix_vo(ispin), &
3803 matrix_pq_template=almo_scf_env%matrix_ov(ispin), &
3804 matrix_v=almo_scf_env%matrix_v(ispin), &
3805 matrix_x_guess=almo_scf_env%matrix_x(ispin))
3810 energy_correction=energy_correction(ispin), &
3811 copy_matrix_x=almo_scf_env%matrix_x(ispin))
3815 energy_correction(1) = energy_correction(1)*spin_factor
3817 IF (opt_k_max_iter .NE. 0)
THEN
3819 CALL timeset(
'k_opt_t_curr', handle3)
3823 almo_scf_env%matrix_v(ispin), &
3824 almo_scf_env%matrix_x(ispin), &
3826 filter_eps=almo_scf_env%eps_filter)
3827 CALL dbcsr_add(t_curr, almo_scf_env%matrix_t_blk(ispin), &
3836 overlap=sigma_oo_curr, &
3837 metric=almo_scf_env%matrix_s(1), &
3838 retain_overlap_sparsity=.false., &
3839 eps_filter=almo_scf_env%eps_filter)
3840 IF (iteration .EQ. 0)
THEN
3843 threshold=almo_scf_env%eps_filter, &
3844 use_inv_as_guess=.false.)
3848 threshold=almo_scf_env%eps_filter, &
3849 use_inv_as_guess=.true.)
3853 CALL dbcsr_create(matrix_tmp1, template=sigma_oo_curr, &
3854 matrix_type=dbcsr_type_no_symmetry)
3856 sigma_oo_curr_inv, &
3857 0.0_dp, matrix_tmp1, &
3858 filter_eps=almo_scf_env%eps_filter)
3864 IF (unit_nr > 0)
THEN
3865 WRITE (unit_nr, *)
"Error for (SIG*inv(SIG)-I)", &
3866 frob_matrix/frob_matrix_base, frob_matrix_base
3871 CALL dbcsr_create(matrix_tmp1, template=sigma_oo_curr, &
3872 matrix_type=dbcsr_type_no_symmetry)
3875 0.0_dp, matrix_tmp1, &
3876 filter_eps=almo_scf_env%eps_filter)
3882 IF (unit_nr > 0)
THEN
3883 WRITE (unit_nr, *)
"Error for (inv(SIG)*SIG-I)", &
3884 frob_matrix/frob_matrix_base, frob_matrix_base
3889 CALL timestop(handle3)
3890 CALL timeset(
'k_opt_vd', handle4)
3898 sigma_vv_sqrt_inv, &
3899 sigma_vv_sqrt_inv, &
3900 0.0_dp, sigma_vv_sqrt, &
3901 filter_eps=almo_scf_env%eps_filter)
3903 psi_out=almo_scf_env%matrix_v_disc(ispin), &
3904 psi_projector=almo_scf_env%matrix_v(ispin), &
3905 metric=almo_scf_env%matrix_s(1), &
3906 project_out=.false., &
3907 psi_projector_orthogonal=.false., &
3908 proj_in_template=almo_scf_env%matrix_k_tr(ispin), &
3909 eps_filter=almo_scf_env%eps_filter, &
3910 sig_inv_projector=sigma_vv_sqrt)
3912 CALL dbcsr_add(almo_scf_env%matrix_v_disc(ispin), &
3913 vd_fixed, -1.0_dp, +1.0_dp)
3915 CALL timestop(handle4)
3916 CALL timeset(
'k_opt_grad', handle5)
3921 IF (line_search)
THEN
3925 almo_scf_env%matrix_ks_0deloc(ispin), &
3928 filter_eps=almo_scf_env%eps_filter)
3930 sigma_oo_curr_inv, &
3931 almo_scf_env%matrix_x(ispin), &
3932 0.0_dp, tmp4_o_vr, &
3933 filter_eps=almo_scf_env%eps_filter)
3937 0.0_dp, tmp1_n_vr, &
3938 filter_eps=almo_scf_env%eps_filter)
3940 almo_scf_env%matrix_v_disc(ispin), &
3943 retain_sparsity=.true.)
3951 converged = (grad_norm .LT. almo_scf_env%opt_k_eps_convergence)
3952 IF (converged .OR. (iteration .GE. opt_k_max_iter))
THEN
3953 prepare_to_exit = .true.
3955 CALL timestop(handle5)
3957 IF (.NOT. prepare_to_exit)
THEN
3959 CALL timeset(
'k_opt_energy', handle6)
3965 0.0_dp, sigma_oo_curr, &
3966 filter_eps=almo_scf_env%eps_filter)
3967 delta_obj_function = fun0
3968 CALL dbcsr_dot(sigma_oo_curr_inv, sigma_oo_curr, obj_function)
3969 delta_obj_function = obj_function - delta_obj_function
3970 IF (line_search)
THEN
3976 CALL timestop(handle6)
3979 IF (.NOT. line_search)
THEN
3981 CALL timeset(
'k_opt_step', handle7)
3983 IF ((.NOT. md_in_k_space) .AND. &
3984 (iteration .GE. max(0, almo_scf_env%opt_k_prec_iter_start) .AND. &
3985 mod(iteration - almo_scf_env%opt_k_prec_iter_start, &
3986 almo_scf_env%opt_k_prec_iter_freq) .EQ. 0))
THEN
3991 IF (unit_nr > 0)
THEN
3992 WRITE (unit_nr, *)
"Computing preconditioner"
4009 CALL opt_k_create_preconditioner_blk(almo_scf_env, &
4010 almo_scf_env%matrix_v_disc(ispin), &
4022 CALL opt_k_apply_preconditioner_blk(almo_scf_env, &
4028 reset_conjugator = .false.
4030 IF (iteration .LT. max(almo_scf_env%opt_k_conj_iter_start, 1) .OR. &
4031 mod(iteration - almo_scf_env%opt_k_conj_iter_start, &
4032 almo_scf_env%opt_k_conj_iter_freq) .EQ. 0)
THEN
4034 reset_conjugator = .true.
4042 CALL dbcsr_dot(grad, prev_minus_prec_grad, numer)
4043 CALL dbcsr_dot(prev_grad, prev_minus_prec_grad, denom)
4044 conjugacy_error = numer/denom
4046 IF (conjugacy_error .GT. min(0.5_dp, conjugacy_error_threshold))
THEN
4047 reset_conjugator = .true.
4048 IF (unit_nr > 0)
THEN
4049 WRITE (unit_nr, *)
"Lack of progress, conjugacy error is ", conjugacy_error
4054 IF ((iteration .NE. 0) .AND. (.NOT. reset_conjugator))
THEN
4056 CALL dbcsr_dot(prev_grad, prev_step, denom)
4057 line_search_error = numer/denom
4058 IF (line_search_error .GT. line_search_error_threshold)
THEN
4059 reset_conjugator = .true.
4060 IF (unit_nr > 0)
THEN
4061 WRITE (unit_nr, *)
"Bad line search, line search error is ", line_search_error
4069 IF (.NOT. reset_conjugator)
THEN
4071 SELECT CASE (almo_scf_env%opt_k_conjugator)
4074 CALL dbcsr_add(tmp_k_blk, prev_grad, 1.0_dp, -1.0_dp)
4076 CALL dbcsr_dot(tmp_k_blk, prev_step, denom)
4077 beta = -1.0_dp*numer/denom
4085 CALL dbcsr_dot(prev_grad, prev_minus_prec_grad, denom)
4093 CALL dbcsr_dot(prev_grad, prev_minus_prec_grad, denom)
4095 CALL dbcsr_add(tmp_k_blk, prev_grad, 1.0_dp, -1.0_dp)
4104 CALL dbcsr_dot(prev_grad, prev_step, denom)
4107 CALL dbcsr_dot(prev_grad, prev_step, denom)
4112 CALL dbcsr_add(tmp_k_blk, prev_grad, 1.0_dp, -1.0_dp)
4123 CALL dbcsr_add(tmp_k_blk, prev_grad, 1.0_dp, -1.0_dp)
4124 CALL dbcsr_dot(tmp_k_blk, prev_step, denom)
4125 beta = -1.0_dp*numer/denom
4137 CALL dbcsr_add(tmp_k_blk, prev_grad, 1.0_dp, -1.0_dp)
4138 CALL dbcsr_dot(tmp_k_blk, prev_step, denom)
4139 CALL dbcsr_dot(tmp_k_blk, prev_minus_prec_grad, numer)
4140 kappa = -2.0_dp*numer/denom
4142 tau = -1.0_dp*numer/denom
4144 beta = tau - kappa*numer/denom
4148 cpabort(
"illegal conjugator")
4151 IF (beta .LT. 0.0_dp)
THEN
4152 IF (unit_nr > 0)
THEN
4153 WRITE (unit_nr, *)
"Beta is negative, ", beta
4155 reset_conjugator = .true.
4160 IF (md_in_k_space)
THEN
4161 reset_conjugator = .true.
4164 IF (reset_conjugator)
THEN
4169 IF (unit_nr > 0)
THEN
4170 WRITE (unit_nr, *)
"(Re)-setting conjugator to zero"
4179 CALL dbcsr_add(step, prev_step, 1.0_dp, beta)
4181 CALL timestop(handle7)
4185 conjugacy_error = 0.0_dp
4189 IF (line_search)
THEN
4191 line_search_error = gfun1/gfun0
4197 IF (line_search)
THEN
4200 safety_multiplier = 1.0e+1_dp
4201 num_threshold = max(epsilon(1.0_dp), &
4202 safety_multiplier*(almo_scf_env%eps_filter**2)*almo_scf_env%ndomains)
4203 IF (abs(fun1 - fun0 - gfun0*step_size) .LT. num_threshold)
THEN
4204 IF (unit_nr > 0)
THEN
4205 WRITE (unit_nr,
'(T3,A,1X,E17.7)') &
4206 "Numerical accuracy is too low to observe non-linear behavior", &
4207 abs(fun1 - fun0 - gfun0*step_size)
4208 WRITE (unit_nr,
'(T3,A,1X,E17.7,A,1X,E12.3)')
"Error computing ", &
4210 " is smaller than the threshold", num_threshold
4214 IF (abs(gfun0) .LT. num_threshold)
THEN
4215 IF (unit_nr > 0)
THEN
4216 WRITE (unit_nr,
'(T3,A,1X,E17.7,A,1X,E12.3)')
"Linear gradient", &
4218 " is smaller than the threshold", num_threshold
4223 use_quadratic_approximation = .true.
4224 use_cubic_approximation = .false.
4228 step_size_quadratic_approx = -(gfun0*step_size*step_size)/(2.0_dp*(fun1 - fun0 - gfun0*step_size))
4230 step_size_quadratic_approx2 = -(fun1 - fun0 - step_size*gfun1/2.0_dp)/(gfun1 - (fun1 - fun0)/step_size)
4232 IF ((step_size_quadratic_approx .LT. 0.0_dp) .AND. &
4233 (step_size_quadratic_approx2 .LT. 0.0_dp))
THEN
4234 IF (unit_nr > 0)
THEN
4235 WRITE (unit_nr,
'(T3,A,1X,E17.7,1X,E17.7,1X,A)') &
4236 "Quadratic approximation gives negative steps", &
4237 step_size_quadratic_approx, step_size_quadratic_approx2, &
4240 use_cubic_approximation = .true.
4241 use_quadratic_approximation = .false.
4243 IF (step_size_quadratic_approx .LT. 0.0_dp)
THEN
4244 step_size_quadratic_approx = step_size_quadratic_approx2
4246 IF (step_size_quadratic_approx2 .LT. 0.0_dp)
THEN
4247 step_size_quadratic_approx2 = step_size_quadratic_approx
4252 IF (use_quadratic_approximation)
THEN
4253 quadratic_approx_error = abs(step_size_quadratic_approx - &
4254 step_size_quadratic_approx2)/step_size_quadratic_approx
4255 IF (quadratic_approx_error .GT. quadratic_approx_error_threshold)
THEN
4256 IF (unit_nr > 0)
THEN
4257 WRITE (unit_nr,
'(T3,A,1X,E17.7,1X,E17.7,1X,A)')
"Quadratic approximation is poor", &
4258 step_size_quadratic_approx, step_size_quadratic_approx2, &
4259 "Try cubic approximation"
4261 use_cubic_approximation = .true.
4262 use_quadratic_approximation = .false.
4267 IF (use_cubic_approximation)
THEN
4272 bb = (-step_size*gfun1 + 3.0_dp*(fun1 - fun0) - 2.0_dp*step_size*gfun0)/(step_size*step_size)
4273 aa = (gfun1 - 2.0_dp*step_size*bb - gfun0)/(3.0_dp*step_size*step_size)
4275 IF (abs(gfun1 - 2.0_dp*step_size*bb - gfun0) .LT. num_threshold)
THEN
4276 IF (unit_nr > 0)
THEN
4277 WRITE (unit_nr,
'(T3,A,1X,E17.7)') &
4278 "Numerical accuracy is too low to observe cubic behavior", &
4279 abs(gfun1 - 2.0_dp*step_size*bb - gfun0)
4281 use_cubic_approximation = .false.
4282 use_quadratic_approximation = .true.
4284 IF (abs(gfun1) .LT. num_threshold)
THEN
4285 IF (unit_nr > 0)
THEN
4286 WRITE (unit_nr,
'(T3,A,1X,E17.7,A,1X,E12.3)')
"Linear gradient", &
4288 " is smaller than the threshold", num_threshold
4290 use_cubic_approximation = .false.
4291 use_quadratic_approximation = .true.
4296 IF (use_cubic_approximation)
THEN
4300 IF (nmins .LT. 1)
THEN
4301 IF (unit_nr > 0)
THEN
4302 WRITE (unit_nr,
'(T3,A)') &
4303 "Cubic approximation gives zero soultions! Use quadratic approximation"
4305 use_quadratic_approximation = .true.
4306 use_cubic_approximation = .true.
4308 step_size = minima(1)
4309 IF (nmins .GT. 1)
THEN
4310 IF (unit_nr > 0)
THEN
4311 WRITE (unit_nr,
'(T3,A)') &
4312 "More than one solution found! Use quadratic approximation"
4314 use_quadratic_approximation = .true.
4315 use_cubic_approximation = .true.
4320 IF (use_quadratic_approximation)
THEN
4321 IF (unit_nr > 0)
THEN
4322 WRITE (unit_nr,
'(T3,A)')
"Use quadratic approximation"
4324 step_size = (step_size_quadratic_approx + step_size_quadratic_approx2)*0.5_dp
4328 IF (step_size .LT. 0.0_dp)
THEN
4329 cpabort(
"Negative step proposed")
4332 CALL dbcsr_copy(almo_scf_env%matrix_k_blk(ispin), &
4334 CALL dbcsr_add(almo_scf_env%matrix_k_blk(ispin), &
4335 step, 1.0_dp, step_size)
4337 almo_scf_env%matrix_k_blk(ispin))
4338 line_search = .false.
4342 IF (md_in_k_space)
THEN
4345 IF (iteration .NE. 0)
THEN
4347 step, 1.0_dp, 0.5_dp*time_step)
4349 prev_step, 1.0_dp, 0.5_dp*time_step)
4352 kin_energy = 0.5_dp*kin_energy*kin_energy
4355 CALL dbcsr_add(almo_scf_env%matrix_k_blk(ispin), &
4356 velocity, 1.0_dp, time_step)
4357 CALL dbcsr_add(almo_scf_env%matrix_k_blk(ispin), &
4358 step, 1.0_dp, 0.5_dp*time_step*time_step)
4362 IF (reset_step_size)
THEN
4363 step_size = almo_scf_env%opt_k_trial_step_size
4364 reset_step_size = .false.
4366 step_size = step_size*almo_scf_env%opt_k_trial_step_size_multiplier
4368 CALL dbcsr_copy(almo_scf_env%matrix_k_blk(ispin), &
4370 CALL dbcsr_add(almo_scf_env%matrix_k_blk(ispin), &
4371 step, 1.0_dp, step_size)
4372 line_search = .true.
4381 IF (unit_nr > 0)
THEN
4382 IF (md_in_k_space)
THEN
4383 WRITE (unit_nr,
'(T6,A,1X,I5,1X,E12.3,E16.7,F15.9,F15.9,F15.9,E12.3,F15.9,F15.9,F8.3)') &
4384 "K iter CG", iteration, time_step, time_step*iteration, &
4385 energy_correction(ispin), obj_function, delta_obj_function, grad_norm, &
4386 kin_energy, kin_energy + obj_function, beta
4388 IF (line_search .OR. prepare_to_exit)
THEN
4389 WRITE (unit_nr,
'(T6,A,1X,I3,1X,E12.3,F16.10,F16.10,E12.3,E12.3,E12.3,F8.3,F8.3,F10.3)') &
4390 "K iter CG", iteration, step_size, &
4391 energy_correction(ispin), delta_obj_function, grad_norm, &
4392 gfun0, line_search_error, beta, conjugacy_error, t2a - t1a
4395 WRITE (unit_nr,
'(T6,A,1X,I3,1X,E12.3,F16.10,F16.10,E12.3,E12.3,E12.3,F8.3,F8.3,F10.3)') &
4396 "K iter LS", iteration, step_size, &
4397 energy_correction(ispin), delta_obj_function, grad_norm, &
4398 gfun1, line_search_error, beta, conjugacy_error, t2a - t1a
4407 prepare_to_exit = .true.
4410 IF (.NOT. line_search) iteration = iteration + 1
4412 IF (prepare_to_exit)
EXIT
4416 IF (converged .OR. (outer_opt_k_iteration .GE. outer_opt_k_max_iter))
THEN
4417 outer_opt_k_prepare_to_exit = .true.
4420 IF (almo_scf_env%deloc_truncate_virt .NE.
virt_full)
THEN
4422 IF (unit_nr > 0)
THEN
4423 WRITE (unit_nr, *)
"Updating ALMO virtuals"
4426 CALL timeset(
'k_opt_v0_update', handle8)
4430 almo_scf_env%matrix_v_disc_blk(ispin), &
4431 almo_scf_env%matrix_k_blk(ispin), &
4433 filter_eps=almo_scf_env%eps_filter)
4434 CALL dbcsr_add(vr_fixed, almo_scf_env%matrix_v_blk(ispin), &
4439 almo_scf_env%matrix_v_blk(ispin), &
4440 almo_scf_env%matrix_k_blk(ispin), &
4442 filter_eps=almo_scf_env%eps_filter)
4443 CALL dbcsr_add(vd_fixed, almo_scf_env%matrix_v_disc_blk(ispin), &
4449 overlap=k_vr_index_down, &
4450 metric=almo_scf_env%matrix_s_blk(1), &
4451 retain_overlap_sparsity=.false., &
4452 eps_filter=almo_scf_env%eps_filter)
4453 CALL dbcsr_create(vr_index_sqrt_inv, template=k_vr_index_down, &
4454 matrix_type=dbcsr_type_no_symmetry)
4455 CALL dbcsr_create(vr_index_sqrt, template=k_vr_index_down, &
4456 matrix_type=dbcsr_type_no_symmetry)
4458 vr_index_sqrt_inv, &
4460 threshold=almo_scf_env%eps_filter, &
4461 order=almo_scf_env%order_lanczos, &
4462 eps_lanczos=almo_scf_env%eps_lanczos, &
4463 max_iter_lanczos=almo_scf_env%max_iter_lanczos)
4465 CALL dbcsr_create(matrix_tmp1, template=k_vr_index_down, &
4466 matrix_type=dbcsr_type_no_symmetry)
4467 CALL dbcsr_create(matrix_tmp2, template=k_vr_index_down, &
4468 matrix_type=dbcsr_type_no_symmetry)
4472 0.0_dp, matrix_tmp1, filter_eps=almo_scf_env%eps_filter)
4474 vr_index_sqrt_inv, &
4475 0.0_dp, matrix_tmp2, filter_eps=almo_scf_env%eps_filter)
4480 IF (unit_nr > 0)
THEN
4481 WRITE (unit_nr, *)
"Error for (inv(sqrt(SIGVV))*SIGVV*inv(sqrt(SIGVV))-I)", &
4482 frob_matrix/frob_matrix_base
4490 vr_index_sqrt_inv, &
4491 0.0_dp, almo_scf_env%matrix_v_blk(ispin), &
4492 filter_eps=almo_scf_env%eps_filter)
4496 overlap=k_vd_index_down, &
4497 metric=almo_scf_env%matrix_s_blk(1), &
4498 retain_overlap_sparsity=.false., &
4499 eps_filter=almo_scf_env%eps_filter)
4500 CALL dbcsr_create(vd_index_sqrt_inv, template=k_vd_index_down, &
4501 matrix_type=dbcsr_type_no_symmetry)
4502 CALL dbcsr_create(vd_index_sqrt, template=k_vd_index_down, &
4503 matrix_type=dbcsr_type_no_symmetry)
4505 vd_index_sqrt_inv, &
4507 threshold=almo_scf_env%eps_filter, &
4508 order=almo_scf_env%order_lanczos, &
4509 eps_lanczos=almo_scf_env%eps_lanczos, &
4510 max_iter_lanczos=almo_scf_env%max_iter_lanczos)
4512 CALL dbcsr_create(matrix_tmp1, template=k_vd_index_down, &
4513 matrix_type=dbcsr_type_no_symmetry)
4514 CALL dbcsr_create(matrix_tmp2, template=k_vd_index_down, &
4515 matrix_type=dbcsr_type_no_symmetry)
4519 0.0_dp, matrix_tmp1, filter_eps=almo_scf_env%eps_filter)
4521 vd_index_sqrt_inv, &
4522 0.0_dp, matrix_tmp2, filter_eps=almo_scf_env%eps_filter)
4527 IF (unit_nr > 0)
THEN
4528 WRITE (unit_nr, *)
"Error for (inv(sqrt(SIGVV))*SIGVV*inv(sqrt(SIGVV))-I)", &
4529 frob_matrix/frob_matrix_base
4537 vd_index_sqrt_inv, &
4538 0.0_dp, almo_scf_env%matrix_v_disc_blk(ispin), &
4539 filter_eps=almo_scf_env%eps_filter)
4546 CALL timestop(handle8)
4553 IF (almo_scf_env%deloc_truncate_virt .NE.
virt_full)
THEN
4574 IF (md_in_k_space)
THEN
4580 outer_opt_k_iteration = outer_opt_k_iteration + 1
4581 IF (outer_opt_k_prepare_to_exit)
EXIT
4595 IF (.NOT. almo_scf_env%s_sqrt_done)
THEN
4597 IF (unit_nr > 0)
THEN
4598 WRITE (unit_nr, *)
"sqrt and inv(sqrt) of AO overlap matrix"
4601 template=almo_scf_env%matrix_s(1), &
4602 matrix_type=dbcsr_type_no_symmetry)
4604 template=almo_scf_env%matrix_s(1), &
4605 matrix_type=dbcsr_type_no_symmetry)
4608 almo_scf_env%matrix_s_sqrt_inv(1), &
4609 almo_scf_env%matrix_s(1), &
4610 threshold=almo_scf_env%eps_filter, &
4611 order=almo_scf_env%order_lanczos, &
4612 eps_lanczos=almo_scf_env%eps_lanczos, &
4613 max_iter_lanczos=almo_scf_env%max_iter_lanczos)
4616 CALL dbcsr_create(matrix_tmp1, template=almo_scf_env%matrix_s(1), &
4617 matrix_type=dbcsr_type_no_symmetry)
4618 CALL dbcsr_create(matrix_tmp2, template=almo_scf_env%matrix_s(1), &
4619 matrix_type=dbcsr_type_no_symmetry)
4621 CALL dbcsr_multiply(
"N",
"N", 1.0_dp, almo_scf_env%matrix_s_sqrt_inv(1), &
4622 almo_scf_env%matrix_s(1), &
4623 0.0_dp, matrix_tmp1, filter_eps=almo_scf_env%eps_filter)
4624 CALL dbcsr_multiply(
"N",
"N", 1.0_dp, matrix_tmp1, almo_scf_env%matrix_s_sqrt_inv(1), &
4625 0.0_dp, matrix_tmp2, filter_eps=almo_scf_env%eps_filter)
4630 IF (unit_nr > 0)
THEN
4631 WRITE (unit_nr, *)
"Error for (inv(sqrt(S))*S*inv(sqrt(S))-I)", frob_matrix/frob_matrix_base
4638 almo_scf_env%s_sqrt_done = .true.
4646 para_env=almo_scf_env%para_env, &
4647 blacs_env=almo_scf_env%blacs_env, &
4648 use_occ_orbs=.true., &
4649 use_virt_orbs=almo_scf_env%deloc_cayley_use_virt_orbs, &
4650 occ_orbs_orthogonal=.false., &
4651 virt_orbs_orthogonal=almo_scf_env%orthogonal_basis, &
4652 tensor_type=almo_scf_env%deloc_cayley_tensor_type, &
4653 neglect_quadratic_term=almo_scf_env%deloc_cayley_linear, &
4654 calculate_energy_corr=.true., &
4657 pp_preconditioner_full=almo_scf_env%deloc_cayley_occ_precond, &
4658 qq_preconditioner_full=almo_scf_env%deloc_cayley_vir_precond, &
4659 eps_convergence=almo_scf_env%deloc_cayley_eps_convergence, &
4660 eps_filter=almo_scf_env%eps_filter, &
4662 q_index_up=almo_scf_env%matrix_s_sqrt_inv(1), &
4663 q_index_down=almo_scf_env%matrix_s_sqrt(1), &
4664 p_index_up=almo_scf_env%matrix_sigma_sqrt_inv(ispin), &
4665 p_index_down=almo_scf_env%matrix_sigma_sqrt(ispin), &
4666 matrix_ks=almo_scf_env%matrix_ks_0deloc(ispin), &
4667 matrix_p=almo_scf_env%matrix_p(ispin), &
4668 matrix_qp_template=almo_scf_env%matrix_t(ispin), &
4669 matrix_pq_template=almo_scf_env%matrix_t_tr(ispin), &
4670 matrix_t=almo_scf_env%matrix_t(ispin), &
4671 conjugator=almo_scf_env%deloc_cayley_conjugator, &
4672 max_iter=almo_scf_env%deloc_cayley_max_iter)
4680 energy_correction=energy_correction(ispin))
4688 energy_correction(1) = energy_correction(1)*spin_factor
4695 IF (unit_nr > 0)
THEN
4697 WRITE (unit_nr,
'(T2,A,I6,F20.9)')
"ECORR", ispin, &
4698 energy_correction(ispin)
4701 energy_correction_final = energy_correction_final + energy_correction(ispin)
4718 p=almo_scf_env%matrix_p(ispin), &
4719 eps_filter=almo_scf_env%eps_filter, &
4720 orthog_orbs=.false., &
4721 nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin), &
4722 s=almo_scf_env%matrix_s(1), &
4723 sigma=almo_scf_env%matrix_sigma(ispin), &
4724 sigma_inv=almo_scf_env%matrix_sigma_inv(ispin), &
4726 algorithm=almo_scf_env%sigma_inv_algorithm, &
4727 inverse_accelerator=almo_scf_env%order_lanczos, &
4728 inv_eps_factor=almo_scf_env%matrix_iter_eps_error_factor, &
4729 eps_lanczos=almo_scf_env%eps_lanczos, &
4730 max_iter_lanczos=almo_scf_env%max_iter_lanczos, &
4731 para_env=almo_scf_env%para_env, &
4732 blacs_env=almo_scf_env%blacs_env)
4734 IF (almo_scf_env%nspins == 1) &
4743 IF (.NOT. almo_scf_env%s_inv_done)
THEN
4744 IF (unit_nr > 0)
THEN
4745 WRITE (unit_nr, *)
"Inverting AO overlap matrix"
4748 template=almo_scf_env%matrix_s(1), &
4749 matrix_type=dbcsr_type_no_symmetry)
4750 IF (.NOT. almo_scf_env%s_sqrt_done)
THEN
4752 almo_scf_env%matrix_s(1), &
4753 threshold=almo_scf_env%eps_filter)
4755 CALL dbcsr_multiply(
"N",
"N", 1.0_dp, almo_scf_env%matrix_s_sqrt_inv(1), &
4756 almo_scf_env%matrix_s_sqrt_inv(1), &
4757 0.0_dp, almo_scf_env%matrix_s_inv(1), &
4758 filter_eps=almo_scf_env%eps_filter)
4762 CALL dbcsr_create(matrix_tmp1, template=almo_scf_env%matrix_s(1), &
4763 matrix_type=dbcsr_type_no_symmetry)
4764 CALL dbcsr_multiply(
"N",
"N", 1.0_dp, almo_scf_env%matrix_s_inv(1), &
4765 almo_scf_env%matrix_s(1), &
4766 0.0_dp, matrix_tmp1, &
4767 filter_eps=almo_scf_env%eps_filter)
4771 IF (unit_nr > 0)
THEN
4772 WRITE (unit_nr, *)
"Error for (inv(S)*S-I)", &
4773 frob_matrix/frob_matrix_base
4778 almo_scf_env%s_inv_done = .true.
4793 ALLOCATE (matrix_p_almo_scf_converged(nspin))
4795 CALL dbcsr_create(matrix_p_almo_scf_converged(ispin), &
4796 template=almo_scf_env%matrix_p(ispin))
4797 CALL dbcsr_copy(matrix_p_almo_scf_converged(ispin), &
4798 almo_scf_env%matrix_p(ispin))
4804 nelectron_spin_real(1) = almo_scf_env%nelectrons_spin(ispin)
4805 IF (almo_scf_env%nspins == 1) &
4806 nelectron_spin_real(1) = nelectron_spin_real(1)/2
4808 local_mu(1) = sum(almo_scf_env%mu_of_domain(:, ispin))/almo_scf_env%ndomains
4814 cpabort(
"CVS only: density_matrix_sign has not been updated in SVN")
4825 almo_scf_env%mu = local_mu(1)
4836 IF (almo_scf_env%nspins == 1) &
4849 CALL dbcsr_add(matrix_p_almo_scf_converged(ispin), &
4850 almo_scf_env%matrix_p(ispin), -1.0_dp, 1.0_dp)
4851 CALL dbcsr_dot(almo_scf_env%matrix_ks_0deloc(ispin), &
4852 matrix_p_almo_scf_converged(ispin), &
4853 energy_correction(ispin))
4855 energy_correction_final = energy_correction_final + energy_correction(ispin)
4857 IF (unit_nr > 0)
THEN
4859 WRITE (unit_nr,
'(T2,A,I6,F20.9)')
"ECORR", ispin, &
4860 energy_correction(ispin)
4869 DEALLOCATE (matrix_p_almo_scf_converged)
4875 IF (unit_nr > 0)
THEN
4877 WRITE (unit_nr,
'(T2,A,F18.9,F18.9,F18.9,F12.6)')
"ETOT", &
4878 almo_scf_env%almo_scf_energy, &
4879 energy_correction_final, &
4880 almo_scf_env%almo_scf_energy + energy_correction_final, &
4885 CALL timestop(handle)
4887 END SUBROUTINE harris_foulkes_correction
4893 SUBROUTINE make_triu(matrix)
4896 CHARACTER(len=*),
PARAMETER :: routinen =
'make_triu'
4898 INTEGER :: col, handle, i, j, row
4899 REAL(
dp),
DIMENSION(:, :),
POINTER :: block
4902 CALL timeset(routinen, handle)
4907 IF (row > col) block(:, :) = 0.0_dp
4908 IF (row == col)
THEN
4909 DO j = 1,
SIZE(block, 2)
4910 DO i = j + 1,
SIZE(block, 1)
4911 block(i, j) = 0.0_dp
4919 CALL timestop(handle)
4920 END SUBROUTINE make_triu
4942 SUBROUTINE opt_k_create_preconditioner(prec, vd_prop, f, x, oo_inv_x_tr, s, grad, &
4943 vd_blk, t, template_vd_vd_blk, template_vr_vr_blk, template_n_vr, &
4944 spin_factor, eps_filter)
4947 TYPE(
dbcsr_type),
INTENT(IN) :: vd_prop, f, x, oo_inv_x_tr, s
4949 TYPE(
dbcsr_type),
INTENT(IN) :: vd_blk, t, template_vd_vd_blk, &
4950 template_vr_vr_blk, template_n_vr
4951 REAL(kind=
dp),
INTENT(IN) :: spin_factor, eps_filter
4953 CHARACTER(len=*),
PARAMETER :: routinen =
'opt_k_create_preconditioner'
4955 INTEGER :: handle, p_nrows, q_nrows
4956 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: p_diagonal, q_diagonal
4957 TYPE(
dbcsr_type) :: pp_diag, qq_diag, t1, t2, tmp, &
4958 tmp1_n_vr, tmp2_n_vr, tmp_n_vd, &
4959 tmp_vd_vd_blk, tmp_vr_vr_blk
4968 CALL timeset(routinen, handle)
4979 0.0_dp, tmp_n_vd, filter_eps=eps_filter)
4981 template=template_vd_vd_blk)
4982 CALL dbcsr_copy(tmp_vd_vd_blk, template_vd_vd_blk)
4984 0.0_dp, tmp_vd_vd_blk, &
4985 retain_sparsity=.true., &
4986 filter_eps=eps_filter)
4989 ALLOCATE (q_diagonal(q_nrows))
4992 template=template_vd_vd_blk)
4997 0.0_dp, t1, filter_eps=eps_filter)
5000 CALL dbcsr_create(tmp_vr_vr_blk, template=template_vr_vr_blk)
5001 CALL dbcsr_copy(tmp_vr_vr_blk, template_vr_vr_blk)
5003 0.0_dp, tmp_vr_vr_blk, &
5004 retain_sparsity=.true., &
5005 filter_eps=eps_filter)
5008 ALLOCATE (p_diagonal(p_nrows))
5010 CALL dbcsr_create(pp_diag, template=template_vr_vr_blk)
5016 0.0_dp, t2, filter_eps=eps_filter)
5022 0.0_dp, tmp_n_vd, filter_eps=eps_filter)
5024 0.0_dp, tmp_vd_vd_blk, &
5025 retain_sparsity=.true., &
5026 filter_eps=eps_filter)
5033 0.0_dp, t1, filter_eps=eps_filter)
5039 0.0_dp, tmp1_n_vr, filter_eps=eps_filter)
5041 0.0_dp, tmp2_n_vr, filter_eps=eps_filter)
5043 0.0_dp, tmp_vr_vr_blk, &
5044 retain_sparsity=.true., &
5045 filter_eps=eps_filter)
5052 0.0_dp, t2, filter_eps=eps_filter)
5055 CALL dbcsr_add(prec, tmp, 1.0_dp, -1.0_dp)
5060 0.0_dp, tmp_n_vd, filter_eps=eps_filter)
5062 0.0_dp, tmp, retain_sparsity=.true., &
5063 filter_eps=eps_filter)
5068 CALL dbcsr_add(prec, t1, 1.0_dp, 1.0_dp)
5070 CALL inverse_of_elements(prec)
5073 DEALLOCATE (q_diagonal)
5074 DEALLOCATE (p_diagonal)
5086 CALL timestop(handle)
5088 END SUBROUTINE opt_k_create_preconditioner
5103 SUBROUTINE opt_k_create_preconditioner_blk(almo_scf_env, vd_prop, oo_inv_x_tr, &
5104 t_curr, ispin, spin_factor)
5107 TYPE(
dbcsr_type),
INTENT(IN) :: vd_prop, oo_inv_x_tr, t_curr
5108 INTEGER,
INTENT(IN) :: ispin
5109 REAL(kind=
dp),
INTENT(IN) :: spin_factor
5111 CHARACTER(len=*),
PARAMETER :: routinen =
'opt_k_create_preconditioner_blk'
5114 REAL(kind=
dp) :: eps_filter
5115 TYPE(
dbcsr_type) :: opt_k_e_dd, opt_k_e_rr, s_dd_sqrt, &
5116 s_rr_sqrt, t1, tmp, tmp1_n_vr, &
5117 tmp2_n_vr, tmp_n_vd, tmp_vd_vd_blk, &
5122 CALL timeset(routinen, handle)
5124 eps_filter = almo_scf_env%eps_filter
5127 CALL dbcsr_create(tmp_n_vd, template=almo_scf_env%matrix_v_disc(ispin))
5129 template=almo_scf_env%matrix_vv_disc_blk(ispin), &
5130 matrix_type=dbcsr_type_no_symmetry)
5132 almo_scf_env%matrix_s(1), &
5134 0.0_dp, tmp_n_vd, filter_eps=eps_filter)
5136 almo_scf_env%matrix_vv_disc_blk(ispin))
5138 0.0_dp, tmp_vd_vd_blk, &
5139 retain_sparsity=.true.)
5142 template=almo_scf_env%matrix_vv_disc_blk(ispin), &
5143 matrix_type=dbcsr_type_no_symmetry)
5145 almo_scf_env%opt_k_t_dd(ispin), &
5147 threshold=eps_filter, &
5148 order=almo_scf_env%order_lanczos, &
5149 eps_lanczos=almo_scf_env%eps_lanczos, &
5150 max_iter_lanczos=almo_scf_env%max_iter_lanczos)
5154 almo_scf_env%matrix_ks_0deloc(ispin), &
5156 0.0_dp, tmp_n_vd, filter_eps=eps_filter)
5158 almo_scf_env%matrix_vv_disc_blk(ispin))
5160 0.0_dp, tmp_vd_vd_blk, &
5161 retain_sparsity=.true.)
5167 almo_scf_env%opt_k_t_dd(ispin), &
5168 0.0_dp, s_dd_sqrt, filter_eps=eps_filter)
5170 almo_scf_env%opt_k_t_dd(ispin), &
5172 0.0_dp, tmp_vd_vd_blk, filter_eps=eps_filter)
5176 template=almo_scf_env%matrix_vv_disc_blk(ispin))
5179 template=almo_scf_env%matrix_vv_disc_blk(ispin), &
5180 matrix_type=dbcsr_type_no_symmetry)
5188 almo_scf_env%opt_k_t_dd(ispin))
5192 0.0_dp, almo_scf_env%opt_k_t_dd(ispin), &
5193 filter_eps=eps_filter)
5199 template=almo_scf_env%matrix_k_blk_ones(ispin))
5201 almo_scf_env%matrix_k_blk_ones(ispin))
5203 template=almo_scf_env%matrix_k_blk_ones(ispin))
5206 0.0_dp, t1, filter_eps=eps_filter)
5211 template=almo_scf_env%matrix_sigma_vv_blk(ispin), &
5212 matrix_type=dbcsr_type_no_symmetry)
5214 almo_scf_env%matrix_sigma_vv_blk(ispin))
5216 almo_scf_env%matrix_x(ispin), &
5218 0.0_dp, tmp_vr_vr_blk, &
5219 retain_sparsity=.true.)
5223 template=almo_scf_env%matrix_sigma_vv_blk(ispin), &
5224 matrix_type=dbcsr_type_no_symmetry)
5226 almo_scf_env%opt_k_t_rr(ispin), &
5228 threshold=eps_filter, &
5229 order=almo_scf_env%order_lanczos, &
5230 eps_lanczos=almo_scf_env%eps_lanczos, &
5231 max_iter_lanczos=almo_scf_env%max_iter_lanczos)
5235 template=almo_scf_env%matrix_v(ispin))
5237 template=almo_scf_env%matrix_v(ispin))
5239 0.0_dp, tmp1_n_vr, filter_eps=eps_filter)
5241 almo_scf_env%matrix_ks_0deloc(ispin), &
5243 0.0_dp, tmp2_n_vr, filter_eps=eps_filter)
5245 0.0_dp, tmp_vr_vr_blk, &
5246 retain_sparsity=.true.)
5253 almo_scf_env%opt_k_t_rr(ispin), &
5254 0.0_dp, s_rr_sqrt, filter_eps=eps_filter)
5256 almo_scf_env%opt_k_t_rr(ispin), &
5258 0.0_dp, tmp_vr_vr_blk, filter_eps=eps_filter)
5262 template=almo_scf_env%matrix_sigma_vv_blk(ispin))
5265 template=almo_scf_env%matrix_sigma_vv_blk(ispin), &
5266 matrix_type=dbcsr_type_no_symmetry)
5274 almo_scf_env%opt_k_t_rr(ispin))
5278 0.0_dp, almo_scf_env%opt_k_t_rr(ispin), &
5279 filter_eps=eps_filter)
5286 0.0_dp, almo_scf_env%opt_k_denom(ispin), &
5287 filter_eps=eps_filter)
5292 CALL dbcsr_add(almo_scf_env%opt_k_denom(ispin), t1, &
5295 CALL dbcsr_scale(almo_scf_env%opt_k_denom(ispin), &
5298 CALL inverse_of_elements(almo_scf_env%opt_k_denom(ispin))
5302 CALL timestop(handle)
5304 END SUBROUTINE opt_k_create_preconditioner_blk
5318 SUBROUTINE opt_k_apply_preconditioner_blk(almo_scf_env, step, grad, ispin)
5323 INTEGER,
INTENT(IN) :: ispin
5325 CHARACTER(len=*),
PARAMETER :: routinen =
'opt_k_apply_preconditioner_blk'
5328 REAL(kind=
dp) :: eps_filter
5331 CALL timeset(routinen, handle)
5333 eps_filter = almo_scf_env%eps_filter
5335 CALL dbcsr_create(tmp_k, template=almo_scf_env%matrix_k_blk(ispin))
5339 grad, almo_scf_env%opt_k_t_rr(ispin), &
5340 0.0_dp, tmp_k, filter_eps=eps_filter)
5342 almo_scf_env%opt_k_t_dd(ispin), tmp_k, &
5343 0.0_dp, step, filter_eps=eps_filter)
5347 almo_scf_env%opt_k_denom(ispin), tmp_k)
5351 almo_scf_env%opt_k_t_dd(ispin), tmp_k, &
5352 0.0_dp, step, filter_eps=eps_filter)
5354 step, almo_scf_env%opt_k_t_rr(ispin), &
5355 0.0_dp, tmp_k, filter_eps=eps_filter)
5361 CALL timestop(handle)
5363 END SUBROUTINE opt_k_apply_preconditioner_blk
5958 SUBROUTINE compute_gradient(m_grad_out, m_ks, m_s, m_t, m_t0, &
5959 m_siginv, m_quench_t, m_FTsiginv, m_siginvTFTsiginv, m_ST, m_STsiginv0, &
5960 m_theta, domain_s_inv, domain_r_down, &
5961 cpu_of_domain, domain_map, assume_t0_q0x, optimize_theta, &
5962 normalize_orbitals, penalty_occ_vol, penalty_occ_local, &
5963 penalty_occ_vol_prefactor, envelope_amplitude, eps_filter, spin_factor, &
5964 special_case, m_sig_sqrti_ii, op_sm_set, weights, energy_coeff, &
5967 TYPE(
dbcsr_type),
INTENT(INOUT) :: m_grad_out, m_ks, m_s, m_t, m_t0, &
5968 m_siginv, m_quench_t, m_ftsiginv, &
5969 m_siginvtftsiginv, m_st, m_stsiginv0, &
5972 INTENT(IN) :: domain_s_inv, domain_r_down
5973 INTEGER,
DIMENSION(:),
INTENT(IN) :: cpu_of_domain
5975 LOGICAL,
INTENT(IN) :: assume_t0_q0x, optimize_theta, &
5976 normalize_orbitals, penalty_occ_vol
5977 LOGICAL,
INTENT(IN),
OPTIONAL :: penalty_occ_local
5978 REAL(kind=
dp),
INTENT(IN) :: penalty_occ_vol_prefactor, &
5979 envelope_amplitude, eps_filter, &
5981 INTEGER,
INTENT(IN) :: special_case
5982 TYPE(
dbcsr_type),
INTENT(IN),
OPTIONAL :: m_sig_sqrti_ii
5984 POINTER :: op_sm_set
5985 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN),
OPTIONAL :: weights
5986 REAL(kind=
dp),
INTENT(IN),
OPTIONAL :: energy_coeff, localiz_coeff
5988 CHARACTER(len=*),
PARAMETER :: routinen =
'compute_gradient'
5990 INTEGER :: dim0, handle, idim0, nao, reim
5991 LOGICAL :: my_penalty_local
5992 REAL(kind=
dp) :: coeff, energy_g_norm, my_energy_coeff, &
5994 penalty_occ_vol_g_norm
5995 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: tg_diagonal
5996 TYPE(
dbcsr_type) :: m_tmp_no_1, m_tmp_no_2, m_tmp_no_3, &
5997 m_tmp_oo_1, m_tmp_oo_2, temp1, temp2, &
5998 tempnocc1, tempoccocc1
6000 CALL timeset(routinen, handle)
6002 IF (normalize_orbitals .AND. (.NOT.
PRESENT(m_sig_sqrti_ii)))
THEN
6003 cpabort(
"Normalization matrix is required")
6006 my_penalty_local = .false.
6007 my_localiz_coeff = 1.0_dp
6008 my_energy_coeff = 0.0_dp
6009 IF (
PRESENT(localiz_coeff))
THEN
6010 my_localiz_coeff = localiz_coeff
6012 IF (
PRESENT(energy_coeff))
THEN
6013 my_energy_coeff = energy_coeff
6015 IF (
PRESENT(penalty_occ_local))
THEN
6016 my_penalty_local = penalty_occ_local
6025 template=m_quench_t, &
6026 matrix_type=dbcsr_type_no_symmetry)
6028 template=m_quench_t, &
6029 matrix_type=dbcsr_type_no_symmetry)
6031 template=m_quench_t, &
6032 matrix_type=dbcsr_type_no_symmetry)
6034 template=m_siginv, &
6035 matrix_type=dbcsr_type_no_symmetry)
6037 template=m_siginv, &
6038 matrix_type=dbcsr_type_no_symmetry)
6041 matrix_type=dbcsr_type_no_symmetry)
6043 template=m_siginv, &
6044 matrix_type=dbcsr_type_no_symmetry)
6047 matrix_type=dbcsr_type_no_symmetry)
6050 matrix_type=dbcsr_type_no_symmetry)
6067 CALL dbcsr_copy(m_tmp_no_2, m_ftsiginv, keep_sparsity=.true.)
6092 m_siginvtftsiginv, &
6093 1.0_dp, m_tmp_no_2, &
6094 retain_sparsity=.true.)
6098 IF (my_penalty_local)
THEN
6102 DO idim0 = 1,
SIZE(op_sm_set, 2)
6104 DO reim = 1,
SIZE(op_sm_set, 1)
6107 op_sm_set(reim, idim0)%matrix, &
6109 0.0_dp, tempnocc1, &
6110 filter_eps=eps_filter)
6116 0.0_dp, tempoccocc1, &
6117 filter_eps=eps_filter)
6120 ALLOCATE (tg_diagonal(dim0))
6124 DEALLOCATE (tg_diagonal)
6130 filter_eps=eps_filter)
6136 cpabort(
"Localization function is not implemented")
6139 coeff = -weights(idim0)
6141 cpabort(
"Localization function is not implemented")
6144 CALL dbcsr_add(temp2, temp1, 1.0_dp, coeff)
6148 CALL dbcsr_add(m_tmp_no_2, temp2, my_energy_coeff, my_localiz_coeff*4.0_dp)
6152 IF (penalty_occ_vol)
THEN
6162 penalty_occ_vol_prefactor, &
6165 0.0_dp, m_tmp_no_1, &
6166 retain_sparsity=.true.)
6171 CALL dbcsr_add(m_tmp_no_2, m_tmp_no_1, 1.0_dp, 1.0_dp)
6175 IF (normalize_orbitals)
THEN
6200 0.0_dp, m_tmp_no_1, &
6201 retain_sparsity=.true.)
6208 0.0_dp, m_tmp_oo_1, &
6209 retain_sparsity=.true.)
6212 ALLOCATE (tg_diagonal(dim0))
6216 DEALLOCATE (tg_diagonal)
6221 0.0_dp, m_tmp_oo_2, &
6222 filter_eps=eps_filter)
6226 1.0_dp, m_tmp_no_1, &
6227 retain_sparsity=.true.)
6236 IF (assume_t0_q0x)
THEN
6242 0.0_dp, m_tmp_oo_1, &
6243 filter_eps=eps_filter)
6247 1.0_dp, m_grad_out, &
6248 filter_eps=eps_filter)
6250 cpabort(
"Cannot project the zero-order space from itself")
6254 matrix_in=m_tmp_no_1, &
6255 matrix_out=m_grad_out, &
6256 operator2=domain_r_down(:), &
6257 operator1=domain_s_inv(:), &
6258 dpattern=m_quench_t, &
6260 node_of_domain=cpu_of_domain, &
6262 filter_eps=eps_filter, &
6264 use_trimmer=.false.)
6298 IF (optimize_theta)
THEN
6300 CALL dtanh_of_elements(m_tmp_no_2, alpha=1.0_dp/envelope_amplitude)
6327 CALL timestop(handle)
6329 END SUBROUTINE compute_gradient
6339 SUBROUTINE print_mathematica_matrix(matrix, filename)
6342 CHARACTER(len=*),
INTENT(IN) :: filename
6344 CHARACTER(len=*),
PARAMETER :: routinen =
'print_mathematica_matrix'
6346 CHARACTER(LEN=20) :: formatstr, scols
6347 INTEGER :: col, fiunit, handle, hori_offset, jj, &
6348 nblkcols_tot, nblkrows_tot, ncols, &
6349 ncores, nrows, row, unit_nr, &
6351 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: ao_block_sizes, mo_block_sizes
6352 INTEGER,
DIMENSION(:),
POINTER :: ao_blk_sizes, mo_blk_sizes
6354 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: h
6355 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: block_p
6360 CALL timeset(routinen, handle)
6364 IF (logger%para_env%is_source())
THEN
6373 IF (ncores .GT. 1)
THEN
6374 cpabort(
"mathematica files: serial code only")
6377 CALL dbcsr_get_info(matrix, row_blk_size=ao_blk_sizes, col_blk_size=mo_blk_sizes, &
6378 nblkrows_total=nblkrows_tot, nblkcols_total=nblkcols_tot)
6379 cpassert(nblkrows_tot == nblkcols_tot)
6380 ALLOCATE (mo_block_sizes(nblkcols_tot), ao_block_sizes(nblkcols_tot))
6381 mo_block_sizes(:) = mo_blk_sizes(:)
6382 ao_block_sizes(:) = ao_blk_sizes(:)
6386 matrix_type=dbcsr_type_no_symmetry)
6389 ncols = sum(mo_block_sizes)
6390 nrows = sum(ao_block_sizes)
6391 ALLOCATE (h(nrows, ncols))
6395 DO col = 1, nblkcols_tot
6398 DO row = 1, nblkrows_tot
6403 h(vert_offset + 1:vert_offset + ao_block_sizes(row), &
6404 hori_offset + 1:hori_offset + mo_block_sizes(col)) &
6409 vert_offset = vert_offset + ao_block_sizes(row)
6413 hori_offset = hori_offset + mo_block_sizes(col)
6419 IF (unit_nr > 0)
THEN
6420 CALL open_file(filename, unit_number=fiunit, file_status=
'REPLACE')
6421 WRITE (scols,
"(I10)") ncols
6422 formatstr =
"("//trim(scols)//
"E27.17)"
6424 WRITE (fiunit, formatstr) h(jj, :)
6429 DEALLOCATE (mo_block_sizes)
6430 DEALLOCATE (ao_block_sizes)
6433 CALL timestop(handle)
6435 END SUBROUTINE print_mathematica_matrix
6457 SUBROUTINE compute_obj_nlmos(localization_obj_function_ispin, penalty_func_ispin, &
6458 penalty_vol_prefactor, overlap_determinant, m_sigma, nocc, m_B0, &
6459 m_theta_normalized, template_matrix_mo, weights, m_S0, just_started, &
6460 penalty_amplitude, eps_filter)
6462 REAL(kind=
dp),
INTENT(INOUT) :: localization_obj_function_ispin, penalty_func_ispin, &
6463 penalty_vol_prefactor, overlap_determinant
6465 INTEGER,
INTENT(IN) :: nocc
6466 TYPE(
dbcsr_type),
DIMENSION(:, :),
INTENT(IN) :: m_b0
6467 TYPE(
dbcsr_type),
INTENT(IN) :: m_theta_normalized, template_matrix_mo
6468 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: weights
6470 LOGICAL,
INTENT(IN) :: just_started
6471 REAL(kind=
dp),
INTENT(IN) :: penalty_amplitude, eps_filter
6473 CHARACTER(len=*),
PARAMETER :: routinen =
'compute_obj_nlmos'
6475 INTEGER :: handle, idim0, ielem, reim
6476 REAL(kind=
dp) :: det1, fval
6477 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: reim_diag, z2
6478 TYPE(
dbcsr_type) :: tempnocc1, tempoccocc1, tempoccocc2
6481 CALL timeset(routinen, handle)
6484 template=template_matrix_mo, &
6485 matrix_type=dbcsr_type_no_symmetry)
6487 template=m_theta_normalized, &
6488 matrix_type=dbcsr_type_no_symmetry)
6490 template=m_theta_normalized, &
6491 matrix_type=dbcsr_type_no_symmetry)
6493 localization_obj_function_ispin = 0.0_dp
6494 penalty_func_ispin = 0.0_dp
6496 ALLOCATE (reim_diag(nocc))
6500 DO idim0 = 1,
SIZE(m_b0, 2)
6504 DO reim = 1,
SIZE(m_b0, 1)
6507 m_b0(reim, idim0), &
6508 m_theta_normalized, &
6509 0.0_dp, tempoccocc1, &
6510 filter_eps=eps_filter)
6514 m_theta_normalized, &
6516 0.0_dp, tempoccocc2, &
6517 retain_sparsity=.true.)
6521 CALL group%sum(reim_diag)
6522 z2(:) = z2(:) + reim_diag(:)*reim_diag(:)
6529 fval = -weights(idim0)*log(abs(z2(ielem)))
6531 fval = weights(idim0) - weights(idim0)*abs(z2(ielem))
6533 fval = weights(idim0) - weights(idim0)*sqrt(abs(z2(ielem)))
6535 localization_obj_function_ispin = localization_obj_function_ispin + fval
6541 DEALLOCATE (reim_diag)
6545 m_theta_normalized, &
6546 0.0_dp, tempoccocc1, &
6547 filter_eps=eps_filter)
6550 m_theta_normalized, &
6553 filter_eps=eps_filter)
6558 overlap_determinant = det1
6560 IF (just_started .AND. penalty_amplitude .LT. 0.0_dp)
THEN
6561 penalty_vol_prefactor = -(-penalty_amplitude)*localization_obj_function_ispin
6563 penalty_func_ispin = penalty_func_ispin + penalty_vol_prefactor*log(det1)
6569 CALL timestop(handle)
6571 END SUBROUTINE compute_obj_nlmos
6589 SUBROUTINE compute_gradient_nlmos(m_grad_out, m_B0, weights, &
6590 m_S0, m_theta_normalized, m_siginv, m_sig_sqrti_ii, &
6591 penalty_vol_prefactor, eps_filter, suggested_vol_penalty)
6593 TYPE(
dbcsr_type),
INTENT(INOUT) :: m_grad_out
6594 TYPE(
dbcsr_type),
DIMENSION(:, :),
INTENT(IN) :: m_b0
6595 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: weights
6596 TYPE(
dbcsr_type),
INTENT(IN) :: m_s0, m_theta_normalized, m_siginv, &
6598 REAL(kind=
dp),
INTENT(IN) :: penalty_vol_prefactor, eps_filter
6599 REAL(kind=
dp),
INTENT(INOUT) :: suggested_vol_penalty
6601 CHARACTER(len=*),
PARAMETER :: routinen =
'compute_gradient_nlmos'
6603 INTEGER :: dim0, handle, idim0, reim
6604 REAL(kind=
dp) :: norm_loc, norm_vol
6605 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: tg_diagonal, z2
6606 TYPE(
dbcsr_type) :: m_temp_oo_1, m_temp_oo_2, m_temp_oo_3, &
6609 CALL timeset(routinen, handle)
6612 template=m_theta_normalized, &
6613 matrix_type=dbcsr_type_no_symmetry)
6615 template=m_theta_normalized, &
6616 matrix_type=dbcsr_type_no_symmetry)
6618 template=m_theta_normalized, &
6619 matrix_type=dbcsr_type_no_symmetry)
6621 template=m_theta_normalized, &
6622 matrix_type=dbcsr_type_no_symmetry)
6625 ALLOCATE (tg_diagonal(dim0))
6630 DO idim0 = 1,
SIZE(m_b0, 2)
6634 DO reim = 1,
SIZE(m_b0, 1)
6637 m_b0(reim, idim0), &
6638 m_theta_normalized, &
6639 0.0_dp, m_temp_oo_3, &
6640 filter_eps=eps_filter)
6645 m_theta_normalized, &
6647 0.0_dp, m_temp_oo_4, &
6648 filter_eps=eps_filter)
6650 tg_diagonal(:) = 0.0_dp
6655 z2(:) = z2(:) + tg_diagonal(:)*tg_diagonal(:)
6660 1.0_dp, m_temp_oo_2, &
6661 filter_eps=eps_filter)
6669 z2(:) = -weights(idim0)/z2(:)
6671 z2(:) = -weights(idim0)
6673 z2(:) = -weights(idim0)/(2*sqrt(z2(:)))
6683 1.0_dp, m_temp_oo_1, &
6684 filter_eps=eps_filter)
6693 m_theta_normalized, &
6694 0.0_dp, m_temp_oo_2, &
6695 filter_eps=eps_filter)
6703 0.0_dp, m_temp_oo_3, &
6704 filter_eps=eps_filter)
6707 suggested_vol_penalty = norm_loc/norm_vol
6708 CALL dbcsr_add(m_temp_oo_1, m_temp_oo_3, &
6709 1.0_dp, 2.0_dp*penalty_vol_prefactor)
6717 0.0_dp, m_grad_out, &
6718 filter_eps=eps_filter)
6723 m_theta_normalized, &
6725 0.0_dp, m_temp_oo_3, &
6726 filter_eps=eps_filter)
6736 0.0_dp, m_temp_oo_1, &
6737 filter_eps=eps_filter)
6742 1.0_dp, m_grad_out, &
6743 filter_eps=eps_filter)
6745 DEALLOCATE (tg_diagonal)
6751 CALL timestop(handle)
6753 END SUBROUTINE compute_gradient_nlmos
6784 SUBROUTINE compute_xalmos_from_main_var(m_var_in, m_t_out, m_quench_t, &
6785 m_t0, m_oo_template, m_STsiginv0, m_s, m_sig_sqrti_ii_out, domain_r_down, &
6786 domain_s_inv, domain_map, cpu_of_domain, assume_t0_q0x, just_started, &
6787 optimize_theta, normalize_orbitals, envelope_amplitude, eps_filter, &
6788 special_case, nocc_of_domain, order_lanczos, eps_lanczos, max_iter_lanczos)
6791 TYPE(
dbcsr_type),
INTENT(INOUT) :: m_t_out, m_quench_t, m_t0, &
6792 m_oo_template, m_stsiginv0, m_s, &
6795 INTENT(IN) :: domain_r_down, domain_s_inv
6797 INTEGER,
DIMENSION(:),
INTENT(IN) :: cpu_of_domain
6798 LOGICAL,
INTENT(IN) :: assume_t0_q0x, just_started, &
6799 optimize_theta, normalize_orbitals
6800 REAL(kind=
dp),
INTENT(IN) :: envelope_amplitude, eps_filter
6801 INTEGER,
INTENT(IN) :: special_case
6802 INTEGER,
DIMENSION(:),
INTENT(IN) :: nocc_of_domain
6803 INTEGER,
INTENT(IN) :: order_lanczos
6804 REAL(kind=
dp),
INTENT(IN) :: eps_lanczos
6805 INTEGER,
INTENT(IN) :: max_iter_lanczos
6807 CHARACTER(len=*),
PARAMETER :: routinen =
'compute_xalmos_from_main_var'
6809 INTEGER :: handle, unit_nr
6810 REAL(kind=
dp) :: t_norm
6814 CALL timeset(routinen, handle)
6818 IF (logger%para_env%is_source())
THEN
6825 template=m_quench_t, &
6826 matrix_type=dbcsr_type_no_symmetry)
6828 template=m_oo_template, &
6829 matrix_type=dbcsr_type_no_symmetry)
6832 IF (optimize_theta)
THEN
6836 IF (unit_nr > 0)
THEN
6837 WRITE (unit_nr, *)
"Maximum norm of the initial guess: ", t_norm
6838 WRITE (unit_nr, *)
"Maximum allowed amplitude: ", &
6841 IF (t_norm .GT. envelope_amplitude .AND. just_started)
THEN
6842 cpabort(
"Max norm of the initial guess is too large")
6845 CALL tanh_of_elements(m_tmp_no_1, alpha=1.0_dp/envelope_amplitude)
6852 IF (assume_t0_q0x)
THEN
6857 0.0_dp, m_tmp_oo_1, &
6858 filter_eps=eps_filter)
6863 filter_eps=eps_filter)
6865 cpabort(
"cannot use projector with block-daigonal ALMOs")
6869 matrix_in=m_t_out, &
6870 matrix_out=m_tmp_no_1, &
6871 operator1=domain_r_down, &
6872 operator2=domain_s_inv, &
6873 dpattern=m_quench_t, &
6875 node_of_domain=cpu_of_domain, &
6877 filter_eps=eps_filter, &
6878 use_trimmer=.false.)
6883 m_t0, 1.0_dp, 1.0_dp)
6886 IF (normalize_orbitals)
THEN
6889 overlap=m_tmp_oo_1, &
6891 retain_locality=.true., &
6892 only_normalize=.true., &
6893 nocc_of_domain=nocc_of_domain(:), &
6894 eps_filter=eps_filter, &
6895 order_lanczos=order_lanczos, &
6896 eps_lanczos=eps_lanczos, &
6897 max_iter_lanczos=max_iter_lanczos, &
6898 overlap_sqrti=m_sig_sqrti_ii_out)
6906 CALL timestop(handle)
6908 END SUBROUTINE compute_xalmos_from_main_var
6946 SUBROUTINE compute_preconditioner(domain_prec_out, m_prec_out, m_ks, m_s, &
6947 m_siginv, m_quench_t, m_FTsiginv, m_siginvTFTsiginv, m_ST, &
6948 m_STsiginv_out, m_s_vv_out, m_f_vv_out, para_env, &
6949 blacs_env, nocc_of_domain, domain_s_inv, domain_s_inv_half, domain_s_half, &
6950 domain_r_down, cpu_of_domain, &
6951 domain_map, assume_t0_q0x, penalty_occ_vol, penalty_occ_vol_prefactor, &
6952 eps_filter, neg_thr, spin_factor, special_case, bad_modes_projector_down_out, &
6956 INTENT(INOUT) :: domain_prec_out
6957 TYPE(
dbcsr_type),
INTENT(INOUT) :: m_prec_out, m_ks, m_s
6958 TYPE(
dbcsr_type),
INTENT(IN) :: m_siginv, m_quench_t, m_ftsiginv, &
6959 m_siginvtftsiginv, m_st
6960 TYPE(
dbcsr_type),
INTENT(INOUT),
OPTIONAL :: m_stsiginv_out, m_s_vv_out, m_f_vv_out
6963 INTEGER,
DIMENSION(:),
INTENT(IN) :: nocc_of_domain
6965 INTENT(IN) :: domain_s_inv
6967 INTENT(IN),
OPTIONAL :: domain_s_inv_half, domain_s_half
6969 INTENT(IN) :: domain_r_down
6970 INTEGER,
DIMENSION(:),
INTENT(IN) :: cpu_of_domain
6972 LOGICAL,
INTENT(IN) :: assume_t0_q0x, penalty_occ_vol
6973 REAL(kind=
dp),
INTENT(IN) :: penalty_occ_vol_prefactor, eps_filter, &
6974 neg_thr, spin_factor
6975 INTEGER,
INTENT(IN) :: special_case
6977 INTENT(INOUT),
OPTIONAL :: bad_modes_projector_down_out
6978 LOGICAL,
INTENT(IN) :: skip_inversion
6980 CHARACTER(len=*),
PARAMETER :: routinen =
'compute_preconditioner'
6982 INTEGER :: handle, ndim, precond_domain_projector
6983 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: nn_diagonal
6986 CALL timeset(routinen, handle)
6990 matrix_type=dbcsr_type_no_symmetry)
6992 template=m_quench_t, &
6993 matrix_type=dbcsr_type_no_symmetry)
7005 0.0_dp, m_tmp_no_3, &
7006 filter_eps=eps_filter)
7009 IF (
PRESENT(m_stsiginv_out))
THEN
7018 1.0_dp, m_tmp_nn_1, &
7019 filter_eps=eps_filter)
7022 IF (
PRESENT(m_s_vv_out))
THEN
7033 1.0_dp, m_prec_out, &
7034 filter_eps=eps_filter)
7038 1.0_dp, m_prec_out, &
7039 filter_eps=eps_filter)
7042 m_siginvtftsiginv, &
7043 0.0_dp, m_tmp_no_3, &
7044 filter_eps=eps_filter)
7048 1.0_dp, m_prec_out, &
7049 filter_eps=eps_filter)
7051 IF (
PRESENT(m_f_vv_out))
THEN
7057 WRITE (unit_nr, *)
"prefactor0:", penalty_occ_vol_prefactor
7065 CALL dbcsr_add(m_prec_out, m_tmp_nn_1, &
7071 IF (penalty_occ_vol)
THEN
7072 CALL dbcsr_add(m_prec_out, m_tmp_nn_1, &
7073 1.0_dp, penalty_occ_vol_prefactor)
7082 IF (skip_inversion)
THEN
7086 ALLOCATE (nn_diagonal(ndim))
7091 DEALLOCATE (nn_diagonal)
7093 CALL dbcsr_copy(m_prec_out, m_tmp_nn_1, keep_sparsity=.true.)
7098 matrix_in=m_tmp_nn_1, &
7099 matrix_out=m_prec_out, &
7100 nocc=nocc_of_domain(:) &
7107 IF (skip_inversion)
THEN
7113 para_env=para_env, &
7114 blacs_env=blacs_env)
7116 para_env=para_env, &
7117 blacs_env=blacs_env, &
7118 uplo_to_full=.true.)
7126 IF (assume_t0_q0x)
THEN
7127 precond_domain_projector = -1
7129 precond_domain_projector = 0
7134 IF (
PRESENT(bad_modes_projector_down_out))
THEN
7136 matrix_main=m_tmp_nn_1, &
7137 subm_s_inv=domain_s_inv(:), &
7138 subm_s_inv_half=domain_s_inv_half(:), &
7139 subm_s_half=domain_s_half(:), &
7140 subm_r_down=domain_r_down(:), &
7141 matrix_trimmer=m_quench_t, &
7142 dpattern=m_quench_t, &
7144 node_of_domain=cpu_of_domain, &
7146 use_trimmer=.false., &
7147 bad_modes_projector_down=bad_modes_projector_down_out(:), &
7148 eps_zero_eigenvalues=neg_thr, &
7149 my_action=precond_domain_projector, &
7150 skip_inversion=skip_inversion &
7154 matrix_main=m_tmp_nn_1, &
7155 subm_s_inv=domain_s_inv(:), &
7156 subm_r_down=domain_r_down(:), &
7157 matrix_trimmer=m_quench_t, &
7158 dpattern=m_quench_t, &
7160 node_of_domain=cpu_of_domain, &
7162 use_trimmer=.false., &
7164 my_action=precond_domain_projector, &
7165 skip_inversion=skip_inversion &
7235 CALL timestop(handle)
7237 END SUBROUTINE compute_preconditioner
7255 SUBROUTINE compute_cg_beta(beta, numer, denom, reset_conjugator, conjugator, &
7256 grad, prev_grad, step, prev_step, prev_minus_prec_grad)
7258 REAL(kind=
dp),
INTENT(INOUT) :: beta
7259 REAL(kind=
dp),
INTENT(INOUT),
OPTIONAL :: numer, denom
7260 LOGICAL,
INTENT(INOUT) :: reset_conjugator
7261 INTEGER,
INTENT(IN) :: conjugator
7262 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(INOUT) :: grad, prev_grad, step, prev_step
7263 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(INOUT), &
7264 OPTIONAL :: prev_minus_prec_grad
7266 CHARACTER(len=*),
PARAMETER :: routinen =
'compute_cg_beta'
7268 INTEGER :: handle, i, nsize, unit_nr
7269 REAL(kind=
dp) :: den, kappa, my_denom, my_numer, &
7270 my_numer2, my_numer3, num, num2, num3, &
7275 CALL timeset(routinen, handle)
7279 IF (logger%para_env%is_source())
THEN
7285 IF (.NOT.
PRESENT(prev_minus_prec_grad))
THEN
7289 cpabort(
"conjugator needs more input")
7294 IF (
PRESENT(numer) .OR.
PRESENT(denom))
THEN
7298 cpabort(
"cannot return numer/denom")
7313 matrix_type=dbcsr_type_no_symmetry)
7315 SELECT CASE (conjugator)
7318 CALL dbcsr_add(m_tmp_no_1, prev_grad(i), &
7320 CALL dbcsr_dot(m_tmp_no_1, step(i), num)
7321 CALL dbcsr_dot(m_tmp_no_1, prev_step(i), den)
7324 CALL dbcsr_dot(prev_grad(i), prev_minus_prec_grad(i), den)
7326 CALL dbcsr_dot(prev_grad(i), prev_minus_prec_grad(i), den)
7328 CALL dbcsr_add(m_tmp_no_1, prev_grad(i), 1.0_dp, -1.0_dp)
7329 CALL dbcsr_dot(m_tmp_no_1, step(i), num)
7332 CALL dbcsr_dot(prev_grad(i), prev_step(i), den)
7334 CALL dbcsr_dot(prev_grad(i), prev_step(i), den)
7336 CALL dbcsr_add(m_tmp_no_1, prev_grad(i), 1.0_dp, -1.0_dp)
7337 CALL dbcsr_dot(m_tmp_no_1, step(i), num)
7341 CALL dbcsr_add(m_tmp_no_1, prev_grad(i), 1.0_dp, -1.0_dp)
7342 CALL dbcsr_dot(m_tmp_no_1, prev_step(i), den)
7345 CALL dbcsr_add(m_tmp_no_1, prev_grad(i), 1.0_dp, -1.0_dp)
7346 CALL dbcsr_dot(m_tmp_no_1, prev_step(i), den)
7347 CALL dbcsr_dot(m_tmp_no_1, prev_minus_prec_grad(i), num)
7348 CALL dbcsr_dot(m_tmp_no_1, step(i), num2)
7349 CALL dbcsr_dot(prev_step(i), grad(i), num3)
7350 my_numer2 = my_numer2 + num2
7351 my_numer3 = my_numer3 + num3
7356 cpabort(
"illegal conjugator")
7358 my_numer = my_numer + num
7359 my_denom = my_denom + den
7367 SELECT CASE (conjugator)
7369 beta = -1.0_dp*my_numer/my_denom
7371 beta = my_numer/my_denom
7373 kappa = -2.0_dp*my_numer/my_denom
7374 tau = -1.0_dp*my_numer2/my_denom
7375 beta = tau - kappa*my_numer3/my_denom
7379 cpabort(
"illegal conjugator")
7384 IF (beta .LT. 0.0_dp)
THEN
7385 IF (unit_nr > 0)
THEN
7386 WRITE (unit_nr, *)
" Resetting conjugator because beta is negative: ", beta
7388 reset_conjugator = .true.
7391 IF (
PRESENT(numer))
THEN
7394 IF (
PRESENT(denom))
THEN
7398 CALL timestop(handle)
7400 END SUBROUTINE compute_cg_beta
7434 SUBROUTINE newton_grad_to_step(optimizer, m_grad, m_delta, m_s, m_ks, &
7435 m_siginv, m_quench_t, m_FTsiginv, m_siginvTFTsiginv, m_ST, m_t, &
7436 m_sig_sqrti_ii, domain_s_inv, domain_r_down, domain_map, cpu_of_domain, &
7437 nocc_of_domain, para_env, blacs_env, eps_filter, optimize_theta, &
7438 penalty_occ_vol, normalize_orbitals, penalty_occ_vol_prefactor, &
7439 penalty_occ_vol_pf2, special_case)
7442 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(IN) :: m_grad
7443 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(INOUT) :: m_delta, m_s, m_ks, m_siginv, m_quench_t
7444 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(IN) :: m_ftsiginv, m_siginvtftsiginv, m_st, &
7447 INTENT(IN) :: domain_s_inv, domain_r_down
7449 INTEGER,
DIMENSION(:),
INTENT(IN) :: cpu_of_domain
7450 INTEGER,
DIMENSION(:, :),
INTENT(IN) :: nocc_of_domain
7453 REAL(kind=
dp),
INTENT(IN) :: eps_filter
7454 LOGICAL,
INTENT(IN) :: optimize_theta, penalty_occ_vol, &
7456 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: penalty_occ_vol_prefactor, &
7458 INTEGER,
INTENT(IN) :: special_case
7460 CHARACTER(len=*),
PARAMETER :: routinen =
'newton_grad_to_step'
7462 CHARACTER(LEN=20) :: iter_type
7463 INTEGER :: handle, ispin, iteration, max_iter, &
7464 ndomains, nspins, outer_iteration, &
7465 outer_max_iter, unit_nr
7466 LOGICAL :: converged, do_exact_inversion, outer_prepare_to_exit, prepare_to_exit, &
7467 reset_conjugator, use_preconditioner
7468 REAL(kind=
dp) :: alpha, beta, denom, denom_ispin, &
7469 eps_error_target, numer, numer_ispin, &
7470 residue_norm, spin_factor, t1, t2
7471 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: residue_max_norm
7474 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:) :: m_f_vo, m_f_vv, m_hstep, m_prec, &
7475 m_residue, m_residue_prev, m_s_vv, &
7476 m_step, m_stsiginv, m_zet, m_zet_prev
7478 DIMENSION(:, :) :: domain_prec
7480 CALL timeset(routinen, handle)
7484 IF (logger%para_env%is_source())
THEN
7491 IF (optimize_theta)
THEN
7492 cpabort(
"theta is NYI")
7497 outer_max_iter = optimizer%max_iter_outer_loop
7498 max_iter = optimizer%max_iter
7499 eps_error_target = optimizer%eps_error
7503 ndomains =
SIZE(domain_s_inv, 1)
7505 IF (nspins == 1)
THEN
7506 spin_factor = 2.0_dp
7508 spin_factor = 1.0_dp
7511 ALLOCATE (domain_prec(ndomains, nspins))
7515 ALLOCATE (m_residue(nspins))
7516 ALLOCATE (m_residue_prev(nspins))
7517 ALLOCATE (m_step(nspins))
7518 ALLOCATE (m_zet(nspins))
7519 ALLOCATE (m_zet_prev(nspins))
7520 ALLOCATE (m_hstep(nspins))
7521 ALLOCATE (m_prec(nspins))
7522 ALLOCATE (m_s_vv(nspins))
7523 ALLOCATE (m_f_vv(nspins))
7524 ALLOCATE (m_f_vo(nspins))
7525 ALLOCATE (m_stsiginv(nspins))
7527 ALLOCATE (residue_max_norm(nspins))
7530 DO ispin = 1, nspins
7534 template=m_quench_t(ispin), &
7535 matrix_type=dbcsr_type_no_symmetry)
7537 template=m_quench_t(ispin), &
7538 matrix_type=dbcsr_type_no_symmetry)
7540 template=m_quench_t(ispin), &
7541 matrix_type=dbcsr_type_no_symmetry)
7543 template=m_quench_t(ispin), &
7544 matrix_type=dbcsr_type_no_symmetry)
7546 template=m_quench_t(ispin), &
7547 matrix_type=dbcsr_type_no_symmetry)
7549 template=m_quench_t(ispin), &
7550 matrix_type=dbcsr_type_no_symmetry)
7552 template=m_quench_t(ispin), &
7553 matrix_type=dbcsr_type_no_symmetry)
7555 template=m_quench_t(ispin), &
7556 matrix_type=dbcsr_type_no_symmetry)
7558 template=m_ks(ispin), &
7559 matrix_type=dbcsr_type_no_symmetry)
7562 matrix_type=dbcsr_type_no_symmetry)
7564 template=m_ks(ispin), &
7565 matrix_type=dbcsr_type_no_symmetry)
7569 CALL dbcsr_copy(m_f_vo(ispin), m_ftsiginv(ispin))
7572 m_siginvtftsiginv(ispin), &
7573 1.0_dp, m_f_vo(ispin), &
7574 filter_eps=eps_filter)
7583 CALL compute_preconditioner( &
7584 domain_prec_out=domain_prec(:, ispin), &
7585 m_prec_out=m_prec(ispin), &
7588 m_siginv=m_siginv(ispin), &
7589 m_quench_t=m_quench_t(ispin), &
7590 m_ftsiginv=m_ftsiginv(ispin), &
7591 m_siginvtftsiginv=m_siginvtftsiginv(ispin), &
7593 m_stsiginv_out=m_stsiginv(ispin), &
7594 m_s_vv_out=m_s_vv(ispin), &
7595 m_f_vv_out=m_f_vv(ispin), &
7596 para_env=para_env, &
7597 blacs_env=blacs_env, &
7598 nocc_of_domain=nocc_of_domain(:, ispin), &
7599 domain_s_inv=domain_s_inv(:, ispin), &
7600 domain_r_down=domain_r_down(:, ispin), &
7601 cpu_of_domain=cpu_of_domain(:), &
7602 domain_map=domain_map(ispin), &
7603 assume_t0_q0x=.false., &
7604 penalty_occ_vol=penalty_occ_vol, &
7605 penalty_occ_vol_prefactor=penalty_occ_vol_prefactor(ispin), &
7606 eps_filter=eps_filter, &
7608 spin_factor=spin_factor, &
7609 special_case=special_case, &
7610 skip_inversion=.false. &
7616 CALL dbcsr_copy(m_delta(ispin), m_quench_t(ispin))
7619 CALL dbcsr_copy(m_residue(ispin), m_grad(ispin))
7622 do_exact_inversion = .false.
7623 IF (do_exact_inversion)
THEN
7627 CALL dbcsr_copy(m_step(ispin), m_grad(ispin))
7631 CALL hessian_diag_apply( &
7632 matrix_grad=m_step(ispin), &
7633 matrix_step=m_zet(ispin), &
7634 matrix_s_ao=m_s_vv(ispin), &
7635 matrix_f_ao=m_f_vv(ispin), &
7638 matrix_s_mo=m_siginv(ispin), &
7639 matrix_f_mo=m_siginvtftsiginv(ispin), &
7640 matrix_s_vo=m_stsiginv(ispin), &
7641 matrix_f_vo=m_f_vo(ispin), &
7642 quench_t=m_quench_t(ispin), &
7643 spin_factor=spin_factor, &
7644 eps_zero=eps_filter*10.0_dp, &
7645 penalty_occ_vol=penalty_occ_vol, &
7646 penalty_occ_vol_prefactor=penalty_occ_vol_prefactor(ispin), &
7647 penalty_occ_vol_pf2=penalty_occ_vol_pf2(ispin), &
7649 para_env=para_env, &
7650 blacs_env=blacs_env &
7657 IF (use_preconditioner)
THEN
7665 0.0_dp, m_zet(ispin), &
7666 filter_eps=eps_filter)
7671 matrix_in=m_residue(ispin), &
7672 matrix_out=m_zet(ispin), &
7673 operator1=domain_prec(:, ispin), &
7674 dpattern=m_quench_t(ispin), &
7675 map=domain_map(ispin), &
7676 node_of_domain=cpu_of_domain(:), &
7678 filter_eps=eps_filter &
7687 CALL dbcsr_copy(m_zet(ispin), m_residue(ispin))
7698 outer_prepare_to_exit = .false.
7700 residue_norm = 0.0_dp
7705 prepare_to_exit = .false.
7713 CALL apply_hessian( &
7718 m_siginv=m_siginv, &
7719 m_quench_t=m_quench_t, &
7720 m_ftsiginv=m_ftsiginv, &
7721 m_siginvtftsiginv=m_siginvtftsiginv, &
7723 m_stsiginv=m_stsiginv, &
7730 m_sig_sqrti_ii=m_sig_sqrti_ii, &
7731 penalty_occ_vol=penalty_occ_vol, &
7732 normalize_orbitals=normalize_orbitals, &
7733 penalty_occ_vol_prefactor=penalty_occ_vol_prefactor, &
7734 eps_filter=eps_filter, &
7735 path_num=hessian_path_reuse &
7741 DO ispin = 1, nspins
7743 CALL dbcsr_dot(m_residue(ispin), m_zet(ispin), numer_ispin)
7744 CALL dbcsr_dot(m_step(ispin), m_hstep(ispin), denom_ispin)
7746 numer = numer + numer_ispin
7747 denom = denom + denom_ispin
7753 DO ispin = 1, nspins
7756 CALL dbcsr_add(m_delta(ispin), m_step(ispin), 1.0_dp, alpha)
7757 CALL dbcsr_copy(m_residue_prev(ispin), m_residue(ispin))
7758 CALL dbcsr_add(m_residue(ispin), m_hstep(ispin), &
7759 1.0_dp, -1.0_dp*alpha)
7760 residue_max_norm(ispin) =
dbcsr_maxabs(m_residue(ispin))
7765 residue_norm = maxval(residue_max_norm)
7766 converged = (residue_norm .LT. eps_error_target)
7767 IF (converged .OR. (iteration .GE. max_iter))
THEN
7768 prepare_to_exit = .true.
7771 IF (.NOT. prepare_to_exit)
THEN
7773 DO ispin = 1, nspins
7776 CALL dbcsr_copy(m_zet_prev(ispin), m_zet(ispin))
7779 IF (use_preconditioner)
THEN
7791 0.0_dp, m_zet(ispin), &
7792 filter_eps=eps_filter)
7797 matrix_in=m_residue(ispin), &
7798 matrix_out=m_zet(ispin), &
7799 operator1=domain_prec(:, ispin), &
7800 dpattern=m_quench_t(ispin), &
7801 map=domain_map(ispin), &
7802 node_of_domain=cpu_of_domain(:), &
7804 filter_eps=eps_filter &
7813 CALL dbcsr_copy(m_zet(ispin), m_residue(ispin))
7820 CALL compute_cg_beta( &
7822 reset_conjugator=reset_conjugator, &
7825 prev_grad=m_residue_prev, &
7827 prev_step=m_zet_prev)
7829 DO ispin = 1, nspins
7832 CALL dbcsr_add(m_step(ispin), m_zet(ispin), beta, 1.0_dp)
7839 IF (unit_nr > 0)
THEN
7841 iter_type = trim(
"NR STEP")
7842 WRITE (unit_nr,
'(T6,A9,I6,F14.5,F14.5,F15.10,F9.2)') &
7843 iter_type, iteration, &
7844 alpha, beta, residue_norm, &
7849 iteration = iteration + 1
7850 IF (prepare_to_exit)
EXIT
7854 IF (converged .OR. (outer_iteration .GE. outer_max_iter))
THEN
7855 outer_prepare_to_exit = .true.
7858 outer_iteration = outer_iteration + 1
7859 IF (outer_prepare_to_exit)
EXIT
7866 IF (penalty_occ_vol)
THEN
7868 DO ispin = 1, nspins
7871 CALL dbcsr_dot(m_delta(ispin), m_zet(ispin), alpha)
7872 WRITE (unit_nr, *)
"trace(grad.delta): ", alpha
7873 alpha = -1.0_dp/(penalty_occ_vol_pf2(ispin)*alpha - 1.0_dp)
7874 WRITE (unit_nr, *)
"correction alpha: ", alpha
7883 DO ispin = 1, nspins
7887 template=m_siginv(ispin), &
7888 matrix_type=dbcsr_type_no_symmetry)
7890 template=m_siginv(ispin), &
7891 matrix_type=dbcsr_type_no_symmetry)
7895 0.0_dp, m_tmp_oo_1, &
7896 filter_eps=eps_filter)
7900 0.0_dp, m_tmp_oo_2, &
7901 filter_eps=eps_filter)
7902 CALL dbcsr_copy(m_zet(ispin), m_quench_t(ispin))
7906 0.0_dp, m_zet(ispin), &
7907 retain_sparsity=.true.)
7909 WRITE (unit_nr,
"(A50,2F20.10)")
"Occupied-space projection of the step", alpha
7910 CALL dbcsr_add(m_zet(ispin), m_delta(ispin), -1.0_dp, 1.0_dp)
7912 WRITE (unit_nr,
"(A50,2F20.10)")
"Virtual-space projection of the step", alpha
7914 WRITE (unit_nr,
"(A50,2F20.10)")
"Full step", alpha
7921 DO ispin = 1, nspins
7935 DEALLOCATE (domain_prec)
7936 DEALLOCATE (m_residue)
7937 DEALLOCATE (m_residue_prev)
7940 DEALLOCATE (m_zet_prev)
7942 DEALLOCATE (m_hstep)
7946 DEALLOCATE (m_stsiginv)
7947 DEALLOCATE (residue_max_norm)
7949 IF (.NOT. converged)
THEN
7950 cpabort(
"Optimization not converged!")
7955 CALL timestop(handle)
7957 END SUBROUTINE newton_grad_to_step
7985 SUBROUTINE apply_hessian(m_x_in, m_x_out, m_ks, m_s, m_siginv, &
7986 m_quench_t, m_FTsiginv, m_siginvTFTsiginv, m_ST, m_STsiginv, m_s_vv, &
7987 m_ks_vv, m_g_full, m_t, m_sig_sqrti_ii, penalty_occ_vol, &
7988 normalize_orbitals, penalty_occ_vol_prefactor, eps_filter, path_num)
7990 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(INOUT) :: m_x_in, m_x_out, m_ks, m_s
7991 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(IN) :: m_siginv, m_quench_t, m_ftsiginv, &
7992 m_siginvtftsiginv, m_st, m_stsiginv
7993 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(INOUT) :: m_s_vv, m_ks_vv, m_g_full
7994 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(IN) :: m_t, m_sig_sqrti_ii
7995 LOGICAL,
INTENT(IN) :: penalty_occ_vol, normalize_orbitals
7996 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: penalty_occ_vol_prefactor
7997 REAL(kind=
dp),
INTENT(IN) :: eps_filter
7998 INTEGER,
INTENT(IN) :: path_num
8000 CHARACTER(len=*),
PARAMETER :: routinen =
'apply_hessian'
8002 INTEGER :: dim0, handle, ispin, nspins
8003 REAL(kind=
dp) :: penalty_prefactor_local, spin_factor
8004 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: tg_diagonal
8005 TYPE(
dbcsr_type) :: m_tmp_no_1, m_tmp_no_2, m_tmp_oo_1, &
8008 CALL timeset(routinen, handle)
8011 IF (penalty_occ_vol) penalty_prefactor_local = 1._dp
8012 cpassert(
SIZE(m_stsiginv) >= 0)
8013 cpassert(
SIZE(m_siginvtftsiginv) >= 0)
8014 cpassert(
SIZE(m_s) >= 0)
8015 cpassert(
SIZE(m_g_full) >= 0)
8016 cpassert(
SIZE(m_ftsiginv) >= 0)
8017 mark_used(m_siginvtftsiginv)
8018 mark_used(m_stsiginv)
8019 mark_used(m_ftsiginv)
8025 IF (nspins .EQ. 1)
THEN
8026 spin_factor = 2.0_dp
8028 spin_factor = 1.0_dp
8031 DO ispin = 1, nspins
8033 penalty_prefactor_local = penalty_occ_vol_prefactor(ispin)/(2.0_dp*spin_factor)
8036 template=m_siginv(ispin), &
8037 matrix_type=dbcsr_type_no_symmetry)
8039 template=m_quench_t(ispin), &
8040 matrix_type=dbcsr_type_no_symmetry)
8042 template=m_quench_t(ispin), &
8043 matrix_type=dbcsr_type_no_symmetry)
8045 template=m_quench_t(ispin), &
8046 matrix_type=dbcsr_type_no_symmetry)
8049 IF (normalize_orbitals)
THEN
8054 CALL dbcsr_copy(m_tmp_oo_1, m_sig_sqrti_ii(ispin))
8058 0.0_dp, m_tmp_oo_1, &
8059 retain_sparsity=.true.)
8061 ALLOCATE (tg_diagonal(dim0))
8065 DEALLOCATE (tg_diagonal)
8071 1.0_dp, m_tmp_no_1, &
8072 filter_eps=eps_filter)
8075 m_sig_sqrti_ii(ispin), &
8076 0.0_dp, m_tmp_x_in, &
8077 filter_eps=eps_filter)
8085 IF (path_num .EQ. hessian_path_reuse)
THEN
8109 0.0_dp, m_tmp_no_1, &
8110 filter_eps=eps_filter)
8111 CALL dbcsr_copy(m_x_out(ispin), m_quench_t(ispin))
8115 0.0_dp, m_x_out(ispin), &
8116 retain_sparsity=.true.)
8128 CALL dbcsr_copy(m_x_out(ispin), m_quench_t(ispin))
8132 0.0_dp, m_x_out(ispin), &
8133 retain_sparsity=.true.)
8135 CALL dbcsr_copy(m_tmp_no_2, m_quench_t(ispin))
8139 0.0_dp, m_tmp_no_2, &
8140 retain_sparsity=.true.)
8141 CALL dbcsr_add(m_x_out(ispin), m_tmp_no_2, &
8142 1.0_dp, -4.0_dp*penalty_prefactor_local + 1.0_dp)
8213 ELSE IF (path_num .EQ. hessian_path_assemble)
THEN
8218 cpabort(
"path is NYI")
8221 cpabort(
"illegal path")
8225 IF (normalize_orbitals)
THEN
8230 CALL dbcsr_copy(m_tmp_oo_1, m_sig_sqrti_ii(ispin))
8234 0.0_dp, m_tmp_oo_1, &
8235 retain_sparsity=.true.)
8237 ALLOCATE (tg_diagonal(dim0))
8241 DEALLOCATE (tg_diagonal)
8246 1.0_dp, m_x_out(ispin), &
8247 retain_sparsity=.true.)
8251 m_sig_sqrti_ii(ispin), &
8252 0.0_dp, m_x_out(ispin), &
8253 retain_sparsity=.true.)
8271 CALL timestop(handle)
8273 END SUBROUTINE apply_hessian
8298 SUBROUTINE hessian_diag_apply(matrix_grad, matrix_step, matrix_S_ao, &
8299 matrix_F_ao, matrix_S_mo, matrix_F_mo, matrix_S_vo, matrix_F_vo, quench_t, &
8300 penalty_occ_vol, penalty_occ_vol_prefactor, penalty_occ_vol_pf2, &
8301 spin_factor, eps_zero, m_s, para_env, blacs_env)
8303 TYPE(
dbcsr_type),
INTENT(INOUT) :: matrix_grad, matrix_step, matrix_s_ao, &
8304 matrix_f_ao, matrix_s_mo
8306 TYPE(
dbcsr_type),
INTENT(INOUT) :: matrix_s_vo, matrix_f_vo, quench_t
8307 LOGICAL,
INTENT(IN) :: penalty_occ_vol
8308 REAL(kind=
dp),
INTENT(IN) :: penalty_occ_vol_prefactor, &
8309 penalty_occ_vol_pf2, spin_factor, &
8315 CHARACTER(len=*),
PARAMETER :: routinen =
'hessian_diag_apply'
8317 INTEGER :: ao_hori_offset, ao_vert_offset, block_col, block_row, col, h_size, handle, ii, &
8318 info, jj, lev1_hori_offset, lev1_vert_offset, lev2_hori_offset, lev2_vert_offset, lwork, &
8319 nblkcols_tot, nblkrows_tot, ncores, orb_i, orb_j, row, unit_nr, zero_neg_eiv
8320 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: ao_block_sizes, ao_domain_sizes, &
8322 INTEGER,
DIMENSION(:),
POINTER :: ao_blk_sizes, mo_blk_sizes
8323 LOGICAL :: found, found_col, found_row
8324 REAL(kind=
dp) :: penalty_prefactor_local, test_error
8325 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: eigenvalues, grad_vec, step_vec, tmp, &
8327 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: f_ao_block, f_mo_block, h, hinv, &
8328 new_block, s_ao_block, s_mo_block, &
8330 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: block_p
8333 TYPE(
dbcsr_type) :: matrix_f_ao_sym, matrix_f_mo_sym, &
8334 matrix_s_ao_sym, matrix_s_mo_sym
8336 CALL timeset(routinen, handle)
8340 IF (logger%para_env%is_source())
THEN
8347 cpassert(
ASSOCIATED(blacs_env))
8348 cpassert(
ASSOCIATED(para_env))
8349 mark_used(blacs_env)
8359 IF (ncores .GT. 1)
THEN
8360 cpabort(
"serial code only")
8363 CALL dbcsr_get_info(quench_t, row_blk_size=ao_blk_sizes, col_blk_size=mo_blk_sizes, &
8364 nblkrows_total=nblkrows_tot, nblkcols_total=nblkcols_tot)
8365 cpassert(nblkrows_tot == nblkcols_tot)
8366 ALLOCATE (mo_block_sizes(nblkcols_tot), ao_block_sizes(nblkcols_tot))
8367 ALLOCATE (ao_domain_sizes(nblkcols_tot))
8368 mo_block_sizes(:) = mo_blk_sizes(:)
8369 ao_block_sizes(:) = ao_blk_sizes(:)
8370 ao_domain_sizes(:) = 0
8373 template=matrix_s_ao, &
8374 matrix_type=dbcsr_type_no_symmetry)
8376 CALL dbcsr_scale(matrix_s_ao_sym, 2.0_dp*spin_factor)
8379 template=matrix_f_ao, &
8380 matrix_type=dbcsr_type_no_symmetry)
8382 CALL dbcsr_scale(matrix_f_ao_sym, 2.0_dp*spin_factor)
8385 template=matrix_s_mo, &
8386 matrix_type=dbcsr_type_no_symmetry)
8390 template=matrix_f_mo, &
8391 matrix_type=dbcsr_type_no_symmetry)
8394 IF (penalty_occ_vol)
THEN
8395 penalty_prefactor_local = penalty_occ_vol_prefactor/(2.0_dp*spin_factor)
8397 penalty_prefactor_local = 0.0_dp
8400 WRITE (unit_nr, *)
"penalty_prefactor_local: ", penalty_prefactor_local
8401 WRITE (unit_nr, *)
"penalty_prefactor_2: ", penalty_occ_vol_pf2
8411 DO col = 1, nblkcols_tot
8414 DO row = 1, nblkrows_tot
8417 row, col, block_p, found)
8419 ao_domain_sizes(col) = ao_domain_sizes(col) + ao_blk_sizes(row)
8424 h_size = h_size + ao_domain_sizes(col)*mo_block_sizes(col)
8428 ALLOCATE (h(h_size, h_size))
8432 lev1_vert_offset = 0
8434 DO row = 1, nblkcols_tot
8436 lev1_hori_offset = 0
8437 DO col = 1, nblkcols_tot
8440 ALLOCATE (f_ao_block(ao_domain_sizes(row), ao_domain_sizes(col)))
8441 ALLOCATE (s_ao_block(ao_domain_sizes(row), ao_domain_sizes(col)))
8442 ALLOCATE (f_mo_block(mo_block_sizes(row), mo_block_sizes(col)))
8443 ALLOCATE (s_mo_block(mo_block_sizes(row), mo_block_sizes(col)))
8445 f_ao_block(:, :) = 0.0_dp
8446 s_ao_block(:, :) = 0.0_dp
8447 f_mo_block(:, :) = 0.0_dp
8448 s_mo_block(:, :) = 0.0_dp
8453 DO block_row = 1, nblkcols_tot
8456 block_row, row, block_p, found_row)
8460 DO block_col = 1, nblkcols_tot
8463 block_col, col, block_p, found_col)
8467 block_row, block_col, block_p, found)
8470 f_ao_block(ao_vert_offset + 1:ao_vert_offset + ao_block_sizes(block_row), &
8471 ao_hori_offset + 1:ao_hori_offset + ao_block_sizes(block_col)) &
8476 block_row, block_col, block_p, found)
8479 s_ao_block(ao_vert_offset + 1:ao_vert_offset + ao_block_sizes(block_row), &
8480 ao_hori_offset + 1:ao_hori_offset + ao_block_sizes(block_col)) &
8484 ao_hori_offset = ao_hori_offset + ao_block_sizes(block_col)
8490 ao_vert_offset = ao_vert_offset + ao_block_sizes(block_row)
8500 f_mo_block(1:mo_block_sizes(row), 1:mo_block_sizes(col)) = block_p(:, :)
8505 s_mo_block(1:mo_block_sizes(row), 1:mo_block_sizes(col)) = block_p(:, :)
8526 lev2_vert_offset = 0
8527 DO orb_j = 1, mo_block_sizes(row)
8529 lev2_hori_offset = 0
8530 DO orb_i = 1, mo_block_sizes(col)
8531 IF (orb_i .EQ. orb_j .AND. row .EQ. col)
THEN
8532 h(lev1_vert_offset + lev2_vert_offset + 1:lev1_vert_offset + lev2_vert_offset + ao_domain_sizes(row), &
8533 lev1_hori_offset + lev2_hori_offset + 1:lev1_hori_offset + lev2_hori_offset + ao_domain_sizes(col)) &
8535 = f_ao_block(:, :) + s_ao_block(:, :)
8545 lev2_hori_offset = lev2_hori_offset + ao_domain_sizes(col)
8549 lev2_vert_offset = lev2_vert_offset + ao_domain_sizes(row)
8553 lev1_hori_offset = lev1_hori_offset + ao_domain_sizes(col)*mo_block_sizes(col)
8555 DEALLOCATE (f_ao_block)
8556 DEALLOCATE (s_ao_block)
8557 DEALLOCATE (f_mo_block)
8558 DEALLOCATE (s_mo_block)
8562 lev1_vert_offset = lev1_vert_offset + ao_domain_sizes(row)*mo_block_sizes(row)
8746 ALLOCATE (grad_vec(h_size))
8747 grad_vec(:) = 0.0_dp
8748 lev1_vert_offset = 0
8750 DO col = 1, nblkcols_tot
8753 lev2_vert_offset = 0
8754 DO row = 1, nblkrows_tot
8757 row, col, block_p, found_row)
8761 row, col, block_p, found)
8764 DO orb_i = 1, mo_block_sizes(col)
8765 grad_vec(lev1_vert_offset + ao_domain_sizes(col)*(orb_i - 1) + lev2_vert_offset + 1: &
8766 lev1_vert_offset + ao_domain_sizes(col)*(orb_i - 1) + lev2_vert_offset + ao_block_sizes(row)) &
8784 lev2_vert_offset = lev2_vert_offset + ao_block_sizes(row)
8790 lev1_vert_offset = lev1_vert_offset + ao_domain_sizes(col)*mo_block_sizes(col)
8802 ALLOCATE (hinv(h_size, h_size))
8803 hinv(:, :) = h(:, :)
8806 ALLOCATE (eigenvalues(h_size))
8809 ALLOCATE (work(max(1, lwork)))
8810 CALL dsyev(
'V',
'L', h_size, hinv, h_size, eigenvalues, work, lwork, info)
8811 lwork = int(work(1))
8814 ALLOCATE (work(max(1, lwork)))
8815 CALL dsyev(
'V',
'L', h_size, hinv, h_size, eigenvalues, work, lwork, info)
8816 IF (info .NE. 0)
THEN
8817 WRITE (unit_nr, *)
'DSYEV ERROR MESSAGE: ', info
8818 cpabort(
"DSYEV failed")
8823 ALLOCATE (step_vec(h_size))
8825 step_vec(:) = matmul(transpose(hinv), grad_vec)
8844 ALLOCATE (test(h_size, h_size))
8847 WRITE (unit_nr,
"(I10,F20.10,F20.10)") jj, eigenvalues(jj), step_vec(jj)
8848 IF (eigenvalues(jj) .GT. eps_zero)
THEN
8849 test(jj, :) = hinv(:, jj)/eigenvalues(jj)
8851 test(jj, :) = hinv(:, jj)*0.0_dp
8852 zero_neg_eiv = zero_neg_eiv + 1
8855 WRITE (unit_nr, *)
'ZERO OR NEGATIVE EIGENVALUES: ', zero_neg_eiv
8856 DEALLOCATE (step_vec)
8858 ALLOCATE (test2(h_size, h_size))
8859 test2(:, :) = matmul(hinv, test)
8860 hinv(:, :) = test2(:, :)
8861 DEALLOCATE (test, test2)
8882 DEALLOCATE (eigenvalues)
8904 ALLOCATE (test(h_size, h_size))
8905 test(:, :) = matmul(hinv, h)
8907 test(ii, ii) = test(ii, ii) - 1.0_dp
8912 test_error = test_error + test(jj, ii)*test(jj, ii)
8915 WRITE (unit_nr, *)
"Hessian inversion error: ", sqrt(test_error)
8919 ALLOCATE (step_vec(h_size))
8920 ALLOCATE (tmp(h_size))
8921 tmp(:) = matmul(hinv, grad_vec)
8923 step_vec(:) = -1.0_dp*tmp(:)
8925 ALLOCATE (tmpr(h_size))
8926 tmpr(:) = matmul(h, step_vec)
8927 tmp(:) = tmpr(:) + grad_vec(:)
8929 WRITE (unit_nr, *)
"NEWTOV step error: ", maxval(abs(tmp))
8935 DEALLOCATE (grad_vec)
8943 template=matrix_grad, &
8944 matrix_type=dbcsr_type_no_symmetry)
8947 lev1_vert_offset = 0
8949 DO col = 1, nblkcols_tot
8952 lev2_vert_offset = 0
8953 DO row = 1, nblkrows_tot
8956 row, col, block_p, found_row)
8959 ALLOCATE (new_block(ao_block_sizes(row), mo_block_sizes(col)))
8960 DO orb_i = 1, mo_block_sizes(col)
8961 new_block(:, orb_i) = &
8962 step_vec(lev1_vert_offset + ao_domain_sizes(col)*(orb_i - 1) + lev2_vert_offset + 1: &
8963 lev1_vert_offset + ao_domain_sizes(col)*(orb_i - 1) + lev2_vert_offset + ao_block_sizes(row))
8966 DEALLOCATE (new_block)
8967 lev2_vert_offset = lev2_vert_offset + ao_block_sizes(row)
8972 lev1_vert_offset = lev1_vert_offset + ao_domain_sizes(col)*mo_block_sizes(col)
8976 DEALLOCATE (step_vec)
8993 DEALLOCATE (mo_block_sizes, ao_block_sizes)
8994 DEALLOCATE (ao_domain_sizes)
8997 template=quench_t, &
8998 matrix_type=dbcsr_type_no_symmetry)
9003 0.0_dp, matrix_s_ao_sym, &
9004 retain_sparsity=.true.)
9006 template=quench_t, &
9007 matrix_type=dbcsr_type_no_symmetry)
9012 0.0_dp, matrix_f_ao_sym, &
9013 retain_sparsity=.true.)
9014 CALL dbcsr_add(matrix_s_ao_sym, matrix_f_ao_sym, &
9016 CALL dbcsr_scale(matrix_s_ao_sym, 2.0_dp*spin_factor)
9017 CALL dbcsr_add(matrix_s_ao_sym, matrix_grad, &
9020 WRITE (unit_nr, *)
"NEWTOL step error: ", test_error
9024 CALL timestop(handle)
9026 END SUBROUTINE hessian_diag_apply
9046 matrix_t_in, matrix_t_out, perturbation_only, &
9052 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:) :: quench_t, matrix_t_in, matrix_t_out
9053 LOGICAL,
INTENT(IN) :: perturbation_only
9054 INTEGER,
INTENT(IN),
OPTIONAL :: special_case
9056 CHARACTER(len=*),
PARAMETER :: routinen =
'almo_scf_xalmo_trustr'
9058 INTEGER :: handle, ispin, iteration, iteration_type_to_report, my_special_case, ndomains, &
9059 nspins, outer_iteration, prec_type, unit_nr
9060 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: nocc
9061 LOGICAL :: assume_t0_q0x, border_reached, inner_loop_success, normalize_orbitals, &
9062 optimize_theta, penalty_occ_vol, reset_conjugator, same_position, scf_converged
9063 REAL(kind=
dp) :: beta, energy_start, energy_trial, eta, expected_reduction, &
9064 fake_step_size_to_report, grad_norm_ratio, grad_norm_ref, loss_change_to_report, &
9065 loss_start, loss_trial, model_grad_norm, penalty_amplitude, penalty_start, penalty_trial, &
9066 radius_current, radius_max, real_temp, rho, spin_factor, step_norm, step_size, t1, &
9067 t1outer, t2, t2outer, y_scalar
9068 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: grad_norm_spin, &
9069 penalty_occ_vol_g_prefactor, &
9070 penalty_occ_vol_h_prefactor
9073 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:) :: ftsiginv, grad, m_model_bd, m_model_d, &
9074 m_model_hessian, m_model_hessian_inv, m_model_r, m_model_r_prev, m_model_rt, &
9075 m_model_rt_prev, m_sig_sqrti_ii, m_theta, m_theta_trial, prev_step, siginvtftsiginv, st, &
9078 DIMENSION(:, :) :: domain_model_hessian_inv, domain_r_down
9081 CALL timeset(routinen, handle)
9086 IF (
PRESENT(special_case)) my_special_case = special_case
9090 IF (logger%para_env%is_source())
THEN
9097 assume_t0_q0x = .false.
9099 optimize_theta = .false.
9101 nspins = almo_scf_env%nspins
9102 IF (nspins == 1)
THEN
9103 spin_factor = 2.0_dp
9105 spin_factor = 1.0_dp
9108 IF (unit_nr > 0)
THEN
9110 SELECT CASE (my_special_case)
9112 WRITE (unit_nr,
'(T2,A,A,A)') repeat(
"-", 20), &
9113 " Optimization of block-diagonal ALMOs ", repeat(
"-", 21)
9115 WRITE (unit_nr,
'(T2,A,A,A)') repeat(
"-", 20), &
9116 " Optimization of fully delocalized MOs ", repeat(
"-", 20)
9118 WRITE (unit_nr,
'(T2,A,A,A)') repeat(
"-", 27), &
9119 " Optimization of XALMOs ", repeat(
"-", 28)
9122 CALL trust_r_report(unit_nr, &
9127 delta_loss=0.0_dp, &
9129 predicted_reduction=0.0_dp, &
9133 WRITE (unit_nr,
'(T2,A)') repeat(
"-", 79)
9137 penalty_occ_vol = .false.
9140 normalize_orbitals = penalty_occ_vol
9141 penalty_amplitude = 0.0_dp
9142 ALLOCATE (penalty_occ_vol_g_prefactor(nspins))
9143 ALLOCATE (penalty_occ_vol_h_prefactor(nspins))
9144 penalty_occ_vol_g_prefactor(:) = 0.0_dp
9145 penalty_occ_vol_h_prefactor(:) = 0.0_dp
9148 prec_type = optimizer%preconditioner
9150 ALLOCATE (grad_norm_spin(nspins))
9151 ALLOCATE (nocc(nspins))
9155 ALLOCATE (m_theta(nspins))
9156 DO ispin = 1, nspins
9158 template=matrix_t_out(ispin), &
9159 matrix_type=dbcsr_type_no_symmetry)
9164 m_t_in=matrix_t_in, &
9165 m_t0=almo_scf_env%matrix_t_blk, &
9166 m_quench_t=quench_t, &
9167 m_overlap=almo_scf_env%matrix_s(1), &
9168 m_sigma_tmpl=almo_scf_env%matrix_sigma_inv, &
9170 xalmo_history=almo_scf_env%xalmo_history, &
9171 assume_t0_q0x=assume_t0_q0x, &
9172 optimize_theta=optimize_theta, &
9173 envelope_amplitude=almo_scf_env%envelope_amplitude, &
9174 eps_filter=almo_scf_env%eps_filter, &
9175 order_lanczos=almo_scf_env%order_lanczos, &
9176 eps_lanczos=almo_scf_env%eps_lanczos, &
9177 max_iter_lanczos=almo_scf_env%max_iter_lanczos, &
9178 nocc_of_domain=almo_scf_env%nocc_of_domain)
9180 ndomains = almo_scf_env%ndomains
9181 ALLOCATE (domain_r_down(ndomains, nspins))
9183 ALLOCATE (domain_model_hessian_inv(ndomains, nspins))
9186 ALLOCATE (m_model_hessian(nspins))
9187 ALLOCATE (m_model_hessian_inv(nspins))
9188 ALLOCATE (siginvtftsiginv(nspins))
9189 ALLOCATE (stsiginv_0(nspins))
9190 ALLOCATE (ftsiginv(nspins))
9191 ALLOCATE (st(nspins))
9192 ALLOCATE (grad(nspins))
9193 ALLOCATE (prev_step(nspins))
9194 ALLOCATE (step(nspins))
9195 ALLOCATE (m_sig_sqrti_ii(nspins))
9196 ALLOCATE (m_model_r(nspins))
9197 ALLOCATE (m_model_rt(nspins))
9198 ALLOCATE (m_model_d(nspins))
9199 ALLOCATE (m_model_bd(nspins))
9200 ALLOCATE (m_model_r_prev(nspins))
9201 ALLOCATE (m_model_rt_prev(nspins))
9202 ALLOCATE (m_theta_trial(nspins))
9204 DO ispin = 1, nspins
9208 template=almo_scf_env%matrix_ks(ispin), &
9209 matrix_type=dbcsr_type_no_symmetry)
9211 template=almo_scf_env%matrix_ks(ispin), &
9212 matrix_type=dbcsr_type_no_symmetry)
9214 template=almo_scf_env%matrix_sigma(ispin), &
9215 matrix_type=dbcsr_type_no_symmetry)
9217 template=matrix_t_out(ispin), &
9218 matrix_type=dbcsr_type_no_symmetry)
9220 template=matrix_t_out(ispin), &
9221 matrix_type=dbcsr_type_no_symmetry)
9223 template=matrix_t_out(ispin), &
9224 matrix_type=dbcsr_type_no_symmetry)
9226 template=matrix_t_out(ispin), &
9227 matrix_type=dbcsr_type_no_symmetry)
9229 template=matrix_t_out(ispin), &
9230 matrix_type=dbcsr_type_no_symmetry)
9232 template=matrix_t_out(ispin), &
9233 matrix_type=dbcsr_type_no_symmetry)
9235 template=almo_scf_env%matrix_sigma_inv(ispin), &
9236 matrix_type=dbcsr_type_no_symmetry)
9238 template=matrix_t_out(ispin), &
9239 matrix_type=dbcsr_type_no_symmetry)
9241 template=matrix_t_out(ispin), &
9242 matrix_type=dbcsr_type_no_symmetry)
9244 template=matrix_t_out(ispin), &
9245 matrix_type=dbcsr_type_no_symmetry)
9247 template=matrix_t_out(ispin), &
9248 matrix_type=dbcsr_type_no_symmetry)
9250 template=matrix_t_out(ispin), &
9251 matrix_type=dbcsr_type_no_symmetry)
9253 template=matrix_t_out(ispin), &
9254 matrix_type=dbcsr_type_no_symmetry)
9256 template=matrix_t_out(ispin), &
9257 matrix_type=dbcsr_type_no_symmetry)
9260 CALL dbcsr_set(prev_step(ispin), 0.0_dp)
9263 nfullrows_total=nocc(ispin))
9271 matrix_s=almo_scf_env%matrix_s(1), &
9272 subm_s_inv=almo_scf_env%domain_s_inv(:, ispin), &
9273 dpattern=quench_t(ispin), &
9274 map=almo_scf_env%domain_map(ispin), &
9275 node_of_domain=almo_scf_env%cpu_of_domain)
9285 template=almo_scf_env%matrix_s(1), &
9286 matrix_type=dbcsr_type_no_symmetry)
9288 almo_scf_env%matrix_s_blk(1), &
9289 threshold=almo_scf_env%eps_filter, &
9290 filter_eps=almo_scf_env%eps_filter)
9296 template=almo_scf_env%matrix_s(1), &
9297 matrix_type=dbcsr_type_no_symmetry)
9300 para_env=almo_scf_env%para_env, &
9301 blacs_env=almo_scf_env%blacs_env)
9303 para_env=almo_scf_env%para_env, &
9304 blacs_env=almo_scf_env%blacs_env, &
9305 uplo_to_full=.true.)
9310 radius_max = optimizer%max_trust_radius
9311 radius_current = min(optimizer%initial_trust_radius, radius_max)
9313 eta = min(max(optimizer%rho_do_not_update, 0.0_dp), 0.25_dp)
9314 energy_start = 0.0_dp
9315 energy_trial = 0.0_dp
9316 penalty_start = 0.0_dp
9317 penalty_trial = 0.0_dp
9321 same_position = .false.
9324 CALL main_var_to_xalmos_and_loss_func( &
9325 almo_scf_env=almo_scf_env, &
9327 m_main_var_in=m_theta, &
9328 m_t_out=matrix_t_out, &
9329 m_sig_sqrti_ii_out=m_sig_sqrti_ii, &
9330 energy_out=energy_start, &
9331 penalty_out=penalty_start, &
9332 m_ftsiginv_out=ftsiginv, &
9333 m_siginvtftsiginv_out=siginvtftsiginv, &
9335 m_stsiginv0_in=stsiginv_0, &
9336 m_quench_t_in=quench_t, &
9337 domain_r_down_in=domain_r_down, &
9338 assume_t0_q0x=assume_t0_q0x, &
9339 just_started=.true., &
9340 optimize_theta=optimize_theta, &
9341 normalize_orbitals=normalize_orbitals, &
9342 perturbation_only=perturbation_only, &
9343 do_penalty=penalty_occ_vol, &
9344 special_case=my_special_case)
9345 loss_start = energy_start + penalty_start
9347 almo_scf_env%almo_scf_energy = energy_start
9349 DO ispin = 1, nspins
9350 IF (penalty_occ_vol)
THEN
9351 penalty_occ_vol_g_prefactor(ispin) = &
9352 -2.0_dp*penalty_amplitude*spin_factor*nocc(ispin)
9353 penalty_occ_vol_h_prefactor(ispin) = 0.0_dp
9358 scf_converged = .false.
9359 adjust_r_loop:
DO outer_iteration = 1, optimizer%max_iter_outer_loop
9362 border_reached = .false.
9364 DO ispin = 1, nspins
9366 CALL dbcsr_filter(step(ispin), almo_scf_env%eps_filter)
9369 IF (.NOT. same_position)
THEN
9371 DO ispin = 1, nspins
9373 IF (unit_nr > 0 .AND. debug_mode)
WRITE (unit_nr, *)
"...Compute model gradient"
9374 CALL compute_gradient( &
9375 m_grad_out=grad(ispin), &
9376 m_ks=almo_scf_env%matrix_ks(ispin), &
9377 m_s=almo_scf_env%matrix_s(1), &
9378 m_t=matrix_t_out(ispin), &
9379 m_t0=almo_scf_env%matrix_t_blk(ispin), &
9380 m_siginv=almo_scf_env%matrix_sigma_inv(ispin), &
9381 m_quench_t=quench_t(ispin), &
9382 m_ftsiginv=ftsiginv(ispin), &
9383 m_siginvtftsiginv=siginvtftsiginv(ispin), &
9385 m_stsiginv0=stsiginv_0(ispin), &
9386 m_theta=m_theta(ispin), &
9387 m_sig_sqrti_ii=m_sig_sqrti_ii(ispin), &
9388 domain_s_inv=almo_scf_env%domain_s_inv(:, ispin), &
9389 domain_r_down=domain_r_down(:, ispin), &
9390 cpu_of_domain=almo_scf_env%cpu_of_domain, &
9391 domain_map=almo_scf_env%domain_map(ispin), &
9392 assume_t0_q0x=assume_t0_q0x, &
9393 optimize_theta=optimize_theta, &
9394 normalize_orbitals=normalize_orbitals, &
9395 penalty_occ_vol=penalty_occ_vol, &
9396 penalty_occ_vol_prefactor=penalty_occ_vol_g_prefactor(ispin), &
9397 envelope_amplitude=almo_scf_env%envelope_amplitude, &
9398 eps_filter=almo_scf_env%eps_filter, &
9399 spin_factor=spin_factor, &
9400 special_case=my_special_case)
9407 DO ispin = 1, nspins
9412 grad_norm_ref = maxval(grad_norm_spin)
9415 CALL trust_r_report(unit_nr, &
9417 iteration=outer_iteration, &
9419 delta_loss=0.0_dp, &
9420 grad_norm=grad_norm_ref, &
9421 predicted_reduction=0.0_dp, &
9423 radius=radius_current, &
9424 new=.NOT. same_position, &
9425 time=t2outer - t1outer)
9428 IF (grad_norm_ref .LE. optimizer%eps_error)
THEN
9429 scf_converged = .true.
9430 border_reached = .false.
9431 expected_reduction = 0.0_dp
9432 IF (.NOT. (optimizer%early_stopping_on .AND. outer_iteration .EQ. 1)) &
9435 scf_converged = .false.
9438 DO ispin = 1, nspins
9440 CALL dbcsr_copy(m_model_r(ispin), grad(ispin))
9446 IF (unit_nr > 0 .AND. debug_mode)
WRITE (unit_nr, *)
"...Multiply Sinv.r"
9450 0.0_dp, m_model_rt(ispin), &
9451 filter_eps=almo_scf_env%eps_filter)
9455 IF (unit_nr > 0 .AND. debug_mode)
WRITE (unit_nr, *)
"...Multiply Sinv_xx.r"
9457 matrix_in=m_model_r(ispin), &
9458 matrix_out=m_model_rt(ispin), &
9459 operator1=almo_scf_env%domain_s_inv(:, ispin), &
9460 dpattern=quench_t(ispin), &
9461 map=almo_scf_env%domain_map(ispin), &
9462 node_of_domain=almo_scf_env%cpu_of_domain, &
9464 filter_eps=almo_scf_env%eps_filter)
9467 cpabort(
"Unknown XALMO special case")
9470 CALL dbcsr_copy(m_model_d(ispin), m_model_rt(ispin))
9475 IF (.NOT. same_position)
THEN
9477 SELECT CASE (prec_type)
9480 IF (unit_nr > 0 .AND. debug_mode)
WRITE (unit_nr, *)
"...Compute model Hessian"
9481 DO ispin = 1, nspins
9482 CALL compute_preconditioner( &
9483 domain_prec_out=almo_scf_env%domain_preconditioner(:, ispin), &
9484 m_prec_out=m_model_hessian(ispin), &
9485 m_ks=almo_scf_env%matrix_ks(ispin), &
9486 m_s=almo_scf_env%matrix_s(1), &
9487 m_siginv=almo_scf_env%matrix_sigma_inv(ispin), &
9488 m_quench_t=quench_t(ispin), &
9489 m_ftsiginv=ftsiginv(ispin), &
9490 m_siginvtftsiginv=siginvtftsiginv(ispin), &
9492 para_env=almo_scf_env%para_env, &
9493 blacs_env=almo_scf_env%blacs_env, &
9494 nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin), &
9495 domain_s_inv=almo_scf_env%domain_s_inv(:, ispin), &
9496 domain_r_down=domain_r_down(:, ispin), &
9497 cpu_of_domain=almo_scf_env%cpu_of_domain, &
9498 domain_map=almo_scf_env%domain_map(ispin), &
9499 assume_t0_q0x=.false., &
9500 penalty_occ_vol=penalty_occ_vol, &
9501 penalty_occ_vol_prefactor=penalty_occ_vol_g_prefactor(ispin), &
9502 eps_filter=almo_scf_env%eps_filter, &
9504 spin_factor=spin_factor, &
9505 skip_inversion=.true., &
9506 special_case=my_special_case)
9511 cpabort(
"Unknown preconditioner")
9518 CALL fixed_r_report(unit_nr, &
9522 border_reached=.false., &
9524 grad_norm_ratio=0.0_dp, &
9527 IF (unit_nr > 0 .AND. debug_mode)
WRITE (unit_nr, *)
"...Start inner loop"
9530 inner_loop_success = .false.
9532 fixed_r_loop:
DO iteration = 1, optimizer%max_iter
9536 DO ispin = 1, nspins
9543 m_model_hessian(ispin), &
9545 0.0_dp, m_model_bd(ispin), &
9546 filter_eps=almo_scf_env%eps_filter)
9551 matrix_in=m_model_d(ispin), &
9552 matrix_out=m_model_bd(ispin), &
9553 operator1=almo_scf_env%domain_preconditioner(:, ispin), &
9554 dpattern=quench_t(ispin), &
9555 map=almo_scf_env%domain_map(ispin), &
9556 node_of_domain=almo_scf_env%cpu_of_domain, &
9558 filter_eps=almo_scf_env%eps_filter)
9563 CALL dbcsr_dot(m_model_d(ispin), m_model_bd(ispin), real_temp)
9564 y_scalar = y_scalar + real_temp
9567 IF (unit_nr > 0 .AND. debug_mode)
WRITE (unit_nr, *)
"...Curvature: ", y_scalar
9570 IF (y_scalar .LT. 0.0_dp)
THEN
9572 CALL step_size_to_border( &
9573 step_size_out=step_size, &
9574 metric_in=almo_scf_env%matrix_s, &
9576 direction_in=m_model_d, &
9577 trust_radius_in=radius_current, &
9578 quench_t_in=quench_t, &
9579 eps_filter_in=almo_scf_env%eps_filter &
9582 DO ispin = 1, nspins
9583 CALL dbcsr_add(step(ispin), m_model_d(ispin), 1.0_dp, step_size)
9586 border_reached = .true.
9587 inner_loop_success = .true.
9589 CALL predicted_reduction( &
9590 reduction_out=expected_reduction, &
9593 hess_in=m_model_hessian, &
9594 hess_submatrix_in=almo_scf_env%domain_preconditioner, &
9595 quench_t_in=quench_t, &
9596 special_case=my_special_case, &
9597 eps_filter=almo_scf_env%eps_filter, &
9598 domain_map=almo_scf_env%domain_map, &
9599 cpu_of_domain=almo_scf_env%cpu_of_domain &
9603 CALL fixed_r_report(unit_nr, &
9605 iteration=iteration, &
9606 step_size=step_size, &
9607 border_reached=border_reached, &
9608 curvature=y_scalar, &
9609 grad_norm_ratio=expected_reduction, &
9618 DO ispin = 1, nspins
9619 CALL dbcsr_dot(m_model_r(ispin), m_model_rt(ispin), real_temp)
9620 step_size = step_size + real_temp
9622 step_size = step_size/y_scalar
9623 IF (unit_nr > 0 .AND. debug_mode)
WRITE (unit_nr, *)
"...Proposed step size: ", step_size
9626 DO ispin = 1, nspins
9627 CALL dbcsr_copy(prev_step(ispin), step(ispin))
9628 CALL dbcsr_add(step(ispin), m_model_d(ispin), 1.0_dp, step_size)
9632 CALL contravariant_matrix_norm( &
9633 norm_out=step_norm, &
9635 metric_in=almo_scf_env%matrix_s, &
9636 quench_t_in=quench_t, &
9637 eps_filter_in=almo_scf_env%eps_filter &
9639 IF (unit_nr > 0 .AND. debug_mode)
WRITE (unit_nr, *)
"...Step norm: ", step_norm
9642 IF (step_norm .GT. radius_current)
THEN
9644 IF (unit_nr > 0 .AND. debug_mode)
WRITE (unit_nr, *)
"...Norm is too large"
9645 CALL step_size_to_border( &
9646 step_size_out=step_size, &
9647 metric_in=almo_scf_env%matrix_s, &
9648 position_in=prev_step, &
9649 direction_in=m_model_d, &
9650 trust_radius_in=radius_current, &
9651 quench_t_in=quench_t, &
9652 eps_filter_in=almo_scf_env%eps_filter &
9654 IF (unit_nr > 0 .AND. debug_mode)
WRITE (unit_nr, *)
"...Step size to border: ", step_size
9656 DO ispin = 1, nspins
9657 CALL dbcsr_copy(step(ispin), prev_step(ispin))
9658 CALL dbcsr_add(step(ispin), m_model_d(ispin), 1.0_dp, step_size)
9661 IF (debug_mode)
THEN
9663 IF (unit_nr > 0)
WRITE (unit_nr, *)
"...Extra norm evaluation"
9664 CALL contravariant_matrix_norm( &
9665 norm_out=step_norm, &
9667 metric_in=almo_scf_env%matrix_s, &
9668 quench_t_in=quench_t, &
9669 eps_filter_in=almo_scf_env%eps_filter &
9671 IF (unit_nr > 0)
WRITE (unit_nr, *)
"...Step norm: ", step_norm
9672 IF (unit_nr > 0)
WRITE (unit_nr, *)
"...Current radius: ", radius_current
9675 border_reached = .true.
9676 inner_loop_success = .true.
9678 CALL predicted_reduction( &
9679 reduction_out=expected_reduction, &
9682 hess_in=m_model_hessian, &
9683 hess_submatrix_in=almo_scf_env%domain_preconditioner, &
9684 quench_t_in=quench_t, &
9685 special_case=my_special_case, &
9686 eps_filter=almo_scf_env%eps_filter, &
9687 domain_map=almo_scf_env%domain_map, &
9688 cpu_of_domain=almo_scf_env%cpu_of_domain &
9692 CALL fixed_r_report(unit_nr, &
9694 iteration=iteration, &
9695 step_size=step_size, &
9696 border_reached=border_reached, &
9697 curvature=y_scalar, &
9698 grad_norm_ratio=expected_reduction, &
9708 border_reached = .false.
9709 inner_loop_success = .true.
9711 CALL predicted_reduction( &
9712 reduction_out=expected_reduction, &
9715 hess_in=m_model_hessian, &
9716 hess_submatrix_in=almo_scf_env%domain_preconditioner, &
9717 quench_t_in=quench_t, &
9718 special_case=my_special_case, &
9719 eps_filter=almo_scf_env%eps_filter, &
9720 domain_map=almo_scf_env%domain_map, &
9721 cpu_of_domain=almo_scf_env%cpu_of_domain &
9725 CALL fixed_r_report(unit_nr, &
9727 iteration=iteration, &
9728 step_size=step_size, &
9729 border_reached=border_reached, &
9730 curvature=y_scalar, &
9731 grad_norm_ratio=expected_reduction, &
9736 ELSE IF (optimizer%trustr_algorithm .EQ.
trustr_dogleg)
THEN
9739 SELECT CASE (prec_type)
9742 IF (unit_nr > 0 .AND. debug_mode)
WRITE (unit_nr, *)
"...Pseudo-invert model Hessian"
9745 DO ispin = 1, nspins
9747 matrix_in=m_model_hessian(ispin), &
9748 matrix_out=m_model_hessian_inv(ispin), &
9749 nocc=almo_scf_env%nocc_of_domain(:, ispin) &
9756 DO ispin = 1, nspins
9757 CALL dbcsr_copy(m_model_hessian_inv(ispin), &
9758 m_model_hessian(ispin))
9760 para_env=almo_scf_env%para_env, &
9761 blacs_env=almo_scf_env%blacs_env)
9763 para_env=almo_scf_env%para_env, &
9764 blacs_env=almo_scf_env%blacs_env, &
9765 uplo_to_full=.true.)
9767 almo_scf_env%eps_filter)
9772 DO ispin = 1, nspins
9774 matrix_main=m_model_hessian(ispin), &
9775 subm_s_inv=almo_scf_env%domain_s_inv(:, ispin), &
9776 subm_r_down=domain_r_down(:, ispin), &
9777 matrix_trimmer=quench_t(ispin), &
9778 dpattern=quench_t(ispin), &
9779 map=almo_scf_env%domain_map(ispin), &
9780 node_of_domain=almo_scf_env%cpu_of_domain, &
9782 use_trimmer=.false., &
9784 skip_inversion=.false. &
9821 cpabort(
"Unknown preconditioner")
9826 DO ispin = 1, nspins
9833 m_model_hessian_inv(ispin), &
9835 0.0_dp, m_model_bd(ispin), &
9836 filter_eps=almo_scf_env%eps_filter)
9841 matrix_in=m_model_r(ispin), &
9842 matrix_out=m_model_bd(ispin), &
9843 operator1=domain_model_hessian_inv(:, ispin), &
9844 dpattern=quench_t(ispin), &
9845 map=almo_scf_env%domain_map(ispin), &
9846 node_of_domain=almo_scf_env%cpu_of_domain, &
9848 filter_eps=almo_scf_env%eps_filter)
9855 CALL contravariant_matrix_norm( &
9856 norm_out=step_norm, &
9857 matrix_in=m_model_bd, &
9858 metric_in=almo_scf_env%matrix_s, &
9859 quench_t_in=quench_t, &
9860 eps_filter_in=almo_scf_env%eps_filter &
9862 IF (unit_nr > 0 .AND. debug_mode)
WRITE (unit_nr, *)
"...pB norm: ", step_norm
9865 IF (step_norm .LE. radius_current)
THEN
9867 IF (unit_nr > 0 .AND. debug_mode)
WRITE (unit_nr, *)
"...Full dogleg"
9869 border_reached = .false.
9871 DO ispin = 1, nspins
9872 CALL dbcsr_copy(step(ispin), m_model_bd(ispin))
9875 fake_step_size_to_report = 2.0_dp
9876 iteration_type_to_report = 6
9880 IF (unit_nr > 0 .AND. debug_mode)
WRITE (unit_nr, *)
"...pB norm is too large"
9882 border_reached = .true.
9886 DO ispin = 1, nspins
9887 CALL dbcsr_add(m_model_bd(ispin), step(ispin), 1.0_dp, -1.0_dp)
9890 CALL step_size_to_border( &
9891 step_size_out=step_size, &
9892 metric_in=almo_scf_env%matrix_s, &
9894 direction_in=m_model_bd, &
9895 trust_radius_in=radius_current, &
9896 quench_t_in=quench_t, &
9897 eps_filter_in=almo_scf_env%eps_filter &
9899 IF (unit_nr > 0 .AND. debug_mode)
WRITE (unit_nr, *)
"...Step size to border: ", step_size
9900 IF (step_size .GT. 1.0_dp .OR. step_size .LT. 0.0_dp)
THEN
9902 WRITE (unit_nr, *)
"Step size (", step_size,
") must lie inside (0,1)"
9903 cpabort(
"Wrong dog leg step. We should never end up here.")
9906 DO ispin = 1, nspins
9907 CALL dbcsr_add(step(ispin), m_model_bd(ispin), 1.0_dp, step_size)
9910 fake_step_size_to_report = 1.0_dp + step_size
9911 iteration_type_to_report = 7
9915 IF (debug_mode)
THEN
9917 IF (unit_nr > 0)
WRITE (unit_nr, *)
"...Extra norm evaluation"
9918 CALL contravariant_matrix_norm( &
9919 norm_out=step_norm, &
9921 metric_in=almo_scf_env%matrix_s, &
9922 quench_t_in=quench_t, &
9923 eps_filter_in=almo_scf_env%eps_filter &
9925 IF (unit_nr > 0)
WRITE (unit_nr, *)
"...Step norm: ", step_norm
9926 IF (unit_nr > 0)
WRITE (unit_nr, *)
"...Current radius: ", radius_current
9929 CALL predicted_reduction( &
9930 reduction_out=expected_reduction, &
9933 hess_in=m_model_hessian, &
9934 hess_submatrix_in=almo_scf_env%domain_preconditioner, &
9935 quench_t_in=quench_t, &
9936 special_case=my_special_case, &
9937 eps_filter=almo_scf_env%eps_filter, &
9938 domain_map=almo_scf_env%domain_map, &
9939 cpu_of_domain=almo_scf_env%cpu_of_domain &
9942 inner_loop_success = .true.
9945 CALL fixed_r_report(unit_nr, &
9946 iter_type=iteration_type_to_report, &
9947 iteration=iteration, &
9948 step_size=fake_step_size_to_report, &
9949 border_reached=border_reached, &
9950 curvature=y_scalar, &
9951 grad_norm_ratio=expected_reduction, &
9959 DO ispin = 1, nspins
9961 CALL dbcsr_copy(m_model_r_prev(ispin), m_model_r(ispin))
9962 CALL dbcsr_add(m_model_r(ispin), m_model_bd(ispin), &
9967 DO ispin = 1, nspins
9972 model_grad_norm = maxval(grad_norm_spin)
9975 grad_norm_ratio = model_grad_norm/grad_norm_ref
9976 IF (grad_norm_ratio .LT. optimizer%model_grad_norm_ratio)
THEN
9978 border_reached = .false.
9979 inner_loop_success = .true.
9981 CALL predicted_reduction( &
9982 reduction_out=expected_reduction, &
9985 hess_in=m_model_hessian, &
9986 hess_submatrix_in=almo_scf_env%domain_preconditioner, &
9987 quench_t_in=quench_t, &
9988 special_case=my_special_case, &
9989 eps_filter=almo_scf_env%eps_filter, &
9990 domain_map=almo_scf_env%domain_map, &
9991 cpu_of_domain=almo_scf_env%cpu_of_domain &
9995 CALL fixed_r_report(unit_nr, &
9997 iteration=iteration, &
9998 step_size=step_size, &
9999 border_reached=border_reached, &
10000 curvature=y_scalar, &
10001 grad_norm_ratio=expected_reduction, &
10009 DO ispin = 1, nspins
10011 CALL dbcsr_copy(m_model_rt_prev(ispin), m_model_rt(ispin))
10014 DO ispin = 1, nspins
10021 m_model_r(ispin), &
10022 0.0_dp, m_model_rt(ispin), &
10023 filter_eps=almo_scf_env%eps_filter)
10028 matrix_in=m_model_r(ispin), &
10029 matrix_out=m_model_rt(ispin), &
10030 operator1=almo_scf_env%domain_s_inv(:, ispin), &
10031 dpattern=quench_t(ispin), &
10032 map=almo_scf_env%domain_map(ispin), &
10033 node_of_domain=almo_scf_env%cpu_of_domain, &
10035 filter_eps=almo_scf_env%eps_filter)
10041 CALL compute_cg_beta( &
10043 reset_conjugator=reset_conjugator, &
10044 conjugator=optimizer%conjugator, &
10045 grad=m_model_r(:), &
10046 prev_grad=m_model_r_prev(:), &
10047 step=m_model_rt(:), &
10048 prev_step=m_model_rt_prev(:) &
10051 DO ispin = 1, nspins
10053 CALL dbcsr_add(m_model_d(ispin), m_model_rt(ispin), beta, 1.0_dp)
10057 CALL fixed_r_report(unit_nr, &
10059 iteration=iteration, &
10060 step_size=step_size, &
10061 border_reached=border_reached, &
10062 curvature=y_scalar, &
10063 grad_norm_ratio=grad_norm_ratio, &
10067 END DO fixed_r_loop
10072 IF (.NOT. inner_loop_success)
THEN
10073 cpabort(
"Inner loop did not produce solution")
10076 DO ispin = 1, nspins
10078 CALL dbcsr_copy(m_theta_trial(ispin), m_theta(ispin))
10079 CALL dbcsr_add(m_theta_trial(ispin), step(ispin), 1.0_dp, 1.0_dp)
10085 CALL main_var_to_xalmos_and_loss_func( &
10086 almo_scf_env=almo_scf_env, &
10088 m_main_var_in=m_theta_trial, &
10089 m_t_out=matrix_t_out, &
10090 m_sig_sqrti_ii_out=m_sig_sqrti_ii, &
10091 energy_out=energy_trial, &
10092 penalty_out=penalty_trial, &
10093 m_ftsiginv_out=ftsiginv, &
10094 m_siginvtftsiginv_out=siginvtftsiginv, &
10096 m_stsiginv0_in=stsiginv_0, &
10097 m_quench_t_in=quench_t, &
10098 domain_r_down_in=domain_r_down, &
10099 assume_t0_q0x=assume_t0_q0x, &
10100 just_started=.false., &
10101 optimize_theta=optimize_theta, &
10102 normalize_orbitals=normalize_orbitals, &
10103 perturbation_only=perturbation_only, &
10104 do_penalty=penalty_occ_vol, &
10105 special_case=my_special_case)
10106 loss_trial = energy_trial + penalty_trial
10109 rho = (loss_trial - loss_start)/expected_reduction
10110 loss_change_to_report = loss_trial - loss_start
10112 IF (rho < 0.25_dp)
THEN
10113 radius_current = 0.25_dp*radius_current
10115 IF (rho > 0.75_dp .AND. border_reached)
THEN
10116 radius_current = min(2.0_dp*radius_current, radius_max)
10120 IF (rho > eta)
THEN
10121 DO ispin = 1, nspins
10122 CALL dbcsr_copy(m_theta(ispin), m_theta_trial(ispin))
10124 loss_start = loss_trial
10125 energy_start = energy_trial
10126 penalty_start = penalty_trial
10127 same_position = .false.
10129 almo_scf_env%almo_scf_energy = energy_trial
10132 same_position = .true.
10134 almo_scf_env%almo_scf_energy = energy_start
10139 CALL trust_r_report(unit_nr, &
10141 iteration=outer_iteration, &
10143 delta_loss=loss_change_to_report, &
10144 grad_norm=0.0_dp, &
10145 predicted_reduction=expected_reduction, &
10147 radius=radius_current, &
10148 new=.NOT. same_position, &
10149 time=t2outer - t1outer)
10152 END DO adjust_r_loop
10155 IF (scf_converged)
THEN
10157 CALL wrap_up_xalmo_scf( &
10159 almo_scf_env=almo_scf_env, &
10160 perturbation_in=perturbation_only, &
10161 m_xalmo_in=matrix_t_out, &
10162 m_quench_in=quench_t, &
10163 energy_inout=energy_start)
10167 DO ispin = 1, nspins
10195 DEALLOCATE (m_model_hessian)
10196 DEALLOCATE (m_model_hessian_inv)
10197 DEALLOCATE (siginvtftsiginv)
10198 DEALLOCATE (stsiginv_0)
10199 DEALLOCATE (ftsiginv)
10202 DEALLOCATE (prev_step)
10204 DEALLOCATE (m_sig_sqrti_ii)
10205 DEALLOCATE (m_model_r)
10206 DEALLOCATE (m_model_rt)
10207 DEALLOCATE (m_model_d)
10208 DEALLOCATE (m_model_bd)
10209 DEALLOCATE (m_model_r_prev)
10210 DEALLOCATE (m_model_rt_prev)
10211 DEALLOCATE (m_theta_trial)
10213 DEALLOCATE (domain_r_down)
10214 DEALLOCATE (domain_model_hessian_inv)
10216 DEALLOCATE (penalty_occ_vol_g_prefactor)
10217 DEALLOCATE (penalty_occ_vol_h_prefactor)
10218 DEALLOCATE (grad_norm_spin)
10221 DEALLOCATE (m_theta)
10223 IF (.NOT. scf_converged .AND. .NOT. optimizer%early_stopping_on)
THEN
10224 cpabort(
"Optimization not converged! ")
10227 CALL timestop(handle)
10260 SUBROUTINE main_var_to_xalmos_and_loss_func(almo_scf_env, qs_env, m_main_var_in, &
10261 m_t_out, energy_out, penalty_out, m_sig_sqrti_ii_out, m_FTsiginv_out, &
10262 m_siginvTFTsiginv_out, m_ST_out, m_STsiginv0_in, m_quench_t_in, domain_r_down_in, &
10263 assume_t0_q0x, just_started, optimize_theta, normalize_orbitals, perturbation_only, &
10264 do_penalty, special_case)
10268 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(IN) :: m_main_var_in
10269 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(INOUT) :: m_t_out
10270 REAL(kind=
dp),
INTENT(OUT) :: energy_out, penalty_out
10271 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(INOUT) :: m_sig_sqrti_ii_out, m_ftsiginv_out, &
10272 m_siginvtftsiginv_out, m_st_out, &
10273 m_stsiginv0_in, m_quench_t_in
10275 INTENT(IN) :: domain_r_down_in
10276 LOGICAL,
INTENT(IN) :: assume_t0_q0x, just_started, &
10277 optimize_theta, normalize_orbitals, &
10278 perturbation_only, do_penalty
10279 INTEGER,
INTENT(IN) :: special_case
10281 CHARACTER(len=*),
PARAMETER :: routinen =
'main_var_to_xalmos_and_loss_func'
10283 INTEGER :: handle, ispin, nspins
10284 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: nocc
10285 REAL(kind=
dp) :: det1, energy_ispin, penalty_amplitude, &
10288 CALL timeset(routinen, handle)
10290 energy_out = 0.0_dp
10291 penalty_out = 0.0_dp
10293 nspins =
SIZE(m_main_var_in)
10294 IF (nspins == 1)
THEN
10295 spin_factor = 2.0_dp
10297 spin_factor = 1.0_dp
10300 penalty_amplitude = 0.0_dp
10302 ALLOCATE (nocc(nspins))
10303 DO ispin = 1, nspins
10305 nfullrows_total=nocc(ispin))
10308 DO ispin = 1, nspins
10311 CALL compute_xalmos_from_main_var( &
10312 m_var_in=m_main_var_in(ispin), &
10313 m_t_out=m_t_out(ispin), &
10314 m_quench_t=m_quench_t_in(ispin), &
10315 m_t0=almo_scf_env%matrix_t_blk(ispin), &
10316 m_oo_template=almo_scf_env%matrix_sigma_inv(ispin), &
10317 m_stsiginv0=m_stsiginv0_in(ispin), &
10318 m_s=almo_scf_env%matrix_s(1), &
10319 m_sig_sqrti_ii_out=m_sig_sqrti_ii_out(ispin), &
10320 domain_r_down=domain_r_down_in(:, ispin), &
10321 domain_s_inv=almo_scf_env%domain_s_inv(:, ispin), &
10322 domain_map=almo_scf_env%domain_map(ispin), &
10323 cpu_of_domain=almo_scf_env%cpu_of_domain, &
10324 assume_t0_q0x=assume_t0_q0x, &
10325 just_started=just_started, &
10326 optimize_theta=optimize_theta, &
10327 normalize_orbitals=normalize_orbitals, &
10328 envelope_amplitude=almo_scf_env%envelope_amplitude, &
10329 eps_filter=almo_scf_env%eps_filter, &
10330 special_case=special_case, &
10331 nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin), &
10332 order_lanczos=almo_scf_env%order_lanczos, &
10333 eps_lanczos=almo_scf_env%eps_lanczos, &
10334 max_iter_lanczos=almo_scf_env%max_iter_lanczos)
10338 t=m_t_out(ispin), &
10339 p=almo_scf_env%matrix_p(ispin), &
10340 eps_filter=almo_scf_env%eps_filter, &
10341 orthog_orbs=.false., &
10342 nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin), &
10343 s=almo_scf_env%matrix_s(1), &
10344 sigma=almo_scf_env%matrix_sigma(ispin), &
10345 sigma_inv=almo_scf_env%matrix_sigma_inv(ispin), &
10346 use_guess=.false., &
10347 algorithm=almo_scf_env%sigma_inv_algorithm, &
10348 inv_eps_factor=almo_scf_env%matrix_iter_eps_error_factor, &
10349 inverse_accelerator=almo_scf_env%order_lanczos, &
10350 eps_lanczos=almo_scf_env%eps_lanczos, &
10351 max_iter_lanczos=almo_scf_env%max_iter_lanczos, &
10352 para_env=almo_scf_env%para_env, &
10353 blacs_env=almo_scf_env%blacs_env)
10362 IF (perturbation_only)
THEN
10364 IF (just_started)
THEN
10365 DO ispin = 1, nspins
10366 CALL dbcsr_copy(almo_scf_env%matrix_ks(ispin), &
10367 almo_scf_env%matrix_ks_0deloc(ispin))
10373 almo_scf_env%matrix_p, &
10374 almo_scf_env%matrix_ks, &
10376 almo_scf_env%eps_filter, &
10377 almo_scf_env%mat_distr_aos)
10380 penalty_out = 0.0_dp
10381 DO ispin = 1, nspins
10383 CALL compute_frequently_used_matrices( &
10384 filter_eps=almo_scf_env%eps_filter, &
10385 m_t_in=m_t_out(ispin), &
10386 m_siginv_in=almo_scf_env%matrix_sigma_inv(ispin), &
10387 m_s_in=almo_scf_env%matrix_s(1), &
10388 m_f_in=almo_scf_env%matrix_ks(ispin), &
10389 m_ftsiginv_out=m_ftsiginv_out(ispin), &
10390 m_siginvtftsiginv_out=m_siginvtftsiginv_out(ispin), &
10391 m_st_out=m_st_out(ispin))
10393 IF (perturbation_only)
THEN
10395 IF (ispin .EQ. 1) energy_out = 0.0_dp
10396 CALL dbcsr_dot(m_t_out(ispin), m_ftsiginv_out(ispin), energy_ispin)
10397 energy_out = energy_out + energy_ispin*spin_factor
10400 IF (do_penalty)
THEN
10402 CALL determinant(almo_scf_env%matrix_sigma(ispin), det1, &
10403 almo_scf_env%eps_filter)
10404 penalty_out = penalty_out - &
10405 penalty_amplitude*spin_factor*nocc(ispin)*log(det1)
10413 CALL timestop(handle)
10415 END SUBROUTINE main_var_to_xalmos_and_loss_func
10432 SUBROUTINE step_size_to_border(step_size_out, metric_in, position_in, &
10433 direction_in, trust_radius_in, quench_t_in, eps_filter_in)
10435 REAL(kind=
dp),
INTENT(INOUT) :: step_size_out
10436 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(IN) :: metric_in, position_in, direction_in
10437 REAL(kind=
dp),
INTENT(IN) :: trust_radius_in
10438 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(IN) :: quench_t_in
10439 REAL(kind=
dp),
INTENT(IN) :: eps_filter_in
10441 INTEGER :: isol, ispin, nsolutions, &
10442 nsolutions_found, nspins
10443 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: nocc
10444 REAL(kind=
dp) :: discrim_sign, discriminant, solution, &
10445 spin_factor, temp_real
10446 REAL(kind=
dp),
DIMENSION(3) :: coef
10447 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:) :: m_temp_no
10449 step_size_out = 0.0_dp
10451 nspins =
SIZE(position_in)
10452 IF (nspins == 1)
THEN
10453 spin_factor = 2.0_dp
10455 spin_factor = 1.0_dp
10458 ALLOCATE (nocc(nspins))
10459 ALLOCATE (m_temp_no(nspins))
10462 DO ispin = 1, nspins
10465 template=direction_in(ispin))
10468 nfullcols_total=nocc(ispin))
10470 CALL dbcsr_copy(m_temp_no(ispin), quench_t_in(ispin))
10473 position_in(ispin), &
10474 0.0_dp, m_temp_no(ispin), &
10475 retain_sparsity=.true.)
10477 CALL dbcsr_dot(position_in(ispin), m_temp_no(ispin), temp_real)
10478 coef(3) = coef(3) + temp_real/nocc(ispin)
10479 CALL dbcsr_dot(direction_in(ispin), m_temp_no(ispin), temp_real)
10480 coef(2) = coef(2) + 2.0_dp*temp_real/nocc(ispin)
10481 CALL dbcsr_copy(m_temp_no(ispin), quench_t_in(ispin))
10484 direction_in(ispin), &
10485 0.0_dp, m_temp_no(ispin), &
10486 retain_sparsity=.true.)
10488 CALL dbcsr_dot(direction_in(ispin), m_temp_no(ispin), temp_real)
10489 coef(1) = coef(1) + temp_real/nocc(ispin)
10496 DEALLOCATE (m_temp_no)
10498 coef(:) = coef(:)*spin_factor
10499 coef(3) = coef(3) - trust_radius_in*trust_radius_in
10502 discriminant = coef(2)*coef(2) - 4.0_dp*coef(1)*coef(3)
10503 IF (discriminant .GT. tiny(discriminant))
THEN
10505 ELSE IF (discriminant .LT. 0.0_dp)
THEN
10507 cpabort(
"Step to border: no solutions")
10512 discrim_sign = 1.0_dp
10513 nsolutions_found = 0
10514 DO isol = 1, nsolutions
10515 solution = (-coef(2) + discrim_sign*sqrt(discriminant))/(2.0_dp*coef(1))
10516 IF (solution .GT. 0.0_dp)
THEN
10517 nsolutions_found = nsolutions_found + 1
10518 step_size_out = solution
10520 discrim_sign = -discrim_sign
10523 IF (nsolutions_found == 0)
THEN
10524 cpabort(
"Step to border: no positive solutions")
10525 ELSE IF (nsolutions_found == 2)
THEN
10526 cpabort(
"Two positive border steps possible!")
10529 END SUBROUTINE step_size_to_border
10542 SUBROUTINE contravariant_matrix_norm(norm_out, matrix_in, metric_in, &
10543 quench_t_in, eps_filter_in)
10545 REAL(kind=
dp),
INTENT(OUT) :: norm_out
10546 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(IN) :: matrix_in, metric_in, quench_t_in
10547 REAL(kind=
dp),
INTENT(IN) :: eps_filter_in
10549 INTEGER :: ispin, nspins
10550 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: nocc
10551 REAL(kind=
dp) :: my_norm, spin_factor, temp_real
10552 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:) :: m_temp_no
10557 nspins =
SIZE(matrix_in)
10558 IF (nspins == 1)
THEN
10559 spin_factor = 2.0_dp
10561 spin_factor = 1.0_dp
10564 ALLOCATE (nocc(nspins))
10565 ALLOCATE (m_temp_no(nspins))
10568 DO ispin = 1, nspins
10570 CALL dbcsr_create(m_temp_no(ispin), template=matrix_in(ispin))
10573 nfullcols_total=nocc(ispin))
10575 CALL dbcsr_copy(m_temp_no(ispin), quench_t_in(ispin))
10578 matrix_in(ispin), &
10579 0.0_dp, m_temp_no(ispin), &
10580 retain_sparsity=.true.)
10582 CALL dbcsr_dot(matrix_in(ispin), m_temp_no(ispin), temp_real)
10584 my_norm = my_norm + temp_real/nocc(ispin)
10591 DEALLOCATE (m_temp_no)
10593 my_norm = my_norm*spin_factor
10594 norm_out = sqrt(my_norm)
10596 END SUBROUTINE contravariant_matrix_norm
10615 SUBROUTINE predicted_reduction(reduction_out, grad_in, step_in, hess_in, &
10616 hess_submatrix_in, quench_t_in, special_case, eps_filter, domain_map, &
10620 REAL(kind=
dp),
INTENT(INOUT) :: reduction_out
10621 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(INOUT) :: grad_in, step_in, hess_in
10623 INTENT(IN) :: hess_submatrix_in
10624 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(IN) :: quench_t_in
10625 INTEGER,
INTENT(IN) :: special_case
10626 REAL(kind=
dp),
INTENT(IN) :: eps_filter
10628 INTEGER,
DIMENSION(:),
INTENT(IN) :: cpu_of_domain
10630 INTEGER :: ispin, nspins
10631 REAL(kind=
dp) :: my_reduction, spin_factor, temp_real
10632 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:) :: m_temp_no
10634 reduction_out = 0.0_dp
10636 nspins =
SIZE(grad_in)
10637 IF (nspins == 1)
THEN
10638 spin_factor = 2.0_dp
10640 spin_factor = 1.0_dp
10643 ALLOCATE (m_temp_no(nspins))
10645 my_reduction = 0.0_dp
10646 DO ispin = 1, nspins
10648 CALL dbcsr_create(m_temp_no(ispin), template=grad_in(ispin))
10650 CALL dbcsr_dot(step_in(ispin), grad_in(ispin), temp_real)
10651 my_reduction = my_reduction + temp_real
10660 0.0_dp, m_temp_no(ispin), &
10661 filter_eps=eps_filter)
10666 matrix_in=step_in(ispin), &
10667 matrix_out=m_temp_no(ispin), &
10668 operator1=hess_submatrix_in(:, ispin), &
10669 dpattern=quench_t_in(ispin), &
10670 map=domain_map(ispin), &
10671 node_of_domain=cpu_of_domain, &
10673 filter_eps=eps_filter)
10678 CALL dbcsr_dot(step_in(ispin), m_temp_no(ispin), temp_real)
10679 my_reduction = my_reduction + 0.5_dp*temp_real
10686 my_reduction = spin_factor*my_reduction
10688 reduction_out = my_reduction
10690 DEALLOCATE (m_temp_no)
10692 END SUBROUTINE predicted_reduction
10709 SUBROUTINE fixed_r_report(unit_nr, iter_type, iteration, step_size, &
10710 border_reached, curvature, grad_norm_ratio, predicted_reduction, time)
10712 INTEGER,
INTENT(IN) :: unit_nr, iter_type, iteration
10713 REAL(kind=
dp),
INTENT(IN) :: step_size
10714 LOGICAL,
INTENT(IN) :: border_reached
10715 REAL(kind=
dp),
INTENT(IN) :: curvature
10716 REAL(kind=
dp),
INTENT(IN),
OPTIONAL :: grad_norm_ratio, predicted_reduction
10717 REAL(kind=
dp),
INTENT(IN) :: time
10719 CHARACTER(LEN=20) :: iter_type_str
10720 REAL(kind=
dp) :: loss_or_grad_change
10722 loss_or_grad_change = 0.0_dp
10723 IF (
PRESENT(grad_norm_ratio))
THEN
10724 loss_or_grad_change = grad_norm_ratio
10725 ELSE IF (
PRESENT(predicted_reduction))
THEN
10726 loss_or_grad_change = predicted_reduction
10728 cpabort(
"one argument is missing")
10731 SELECT CASE (iter_type)
10733 iter_type_str = trim(
"Ignored")
10735 iter_type_str = trim(
"PCG")
10737 iter_type_str = trim(
"Neg. curvatr.")
10739 iter_type_str = trim(
"Step too long")
10741 iter_type_str = trim(
"Grad. reduced")
10743 iter_type_str = trim(
"Cauchy point")
10745 iter_type_str = trim(
"Full dogleg")
10747 iter_type_str = trim(
"Part. dogleg")
10749 cpabort(
"unknown report type")
10752 IF (unit_nr > 0)
THEN
10754 SELECT CASE (iter_type)
10758 WRITE (unit_nr,
'(T4,A15,A6,A10,A10,A7,A20,A8)') &
10764 "Grad/o.f. reduc", &
10769 WRITE (unit_nr,
'(T4,A15,I6,F10.5,F10.5,L7,F20.10,F8.2)') &
10772 curvature, step_size, border_reached, &
10773 loss_or_grad_change, &
10779 SELECT CASE (iter_type)
10780 CASE (2, 3, 4, 5, 6, 7)
10788 END SUBROUTINE fixed_r_report
10807 SUBROUTINE trust_r_report(unit_nr, iter_type, iteration, radius, &
10808 loss, delta_loss, grad_norm, predicted_reduction, rho, new, time)
10810 INTEGER,
INTENT(IN) :: unit_nr, iter_type, iteration
10811 REAL(kind=
dp),
INTENT(IN) :: radius, loss, delta_loss, grad_norm, &
10812 predicted_reduction, rho
10813 LOGICAL,
INTENT(IN) :: new
10814 REAL(kind=
dp),
INTENT(IN) :: time
10816 CHARACTER(LEN=20) :: iter_status, iter_type_str
10818 SELECT CASE (iter_type)
10820 iter_type_str = trim(
"Iter")
10821 iter_status = trim(
"Stat")
10823 iter_type_str = trim(
"TR INI")
10825 iter_status =
" New"
10827 iter_status =
" Redo"
10830 iter_type_str = trim(
"TR FIN")
10832 iter_status =
" Acc"
10834 iter_status =
" Rej"
10837 cpabort(
"unknown report type")
10840 IF (unit_nr > 0)
THEN
10842 SELECT CASE (iter_type)
10845 WRITE (unit_nr,
'(T2,A6,A5,A6,A22,A10,T67,A7,A6)') &
10849 "Objective Function", &
10853 WRITE (unit_nr,
'(T41,A10,A10,A6)') &
10857 "Change",
"Expct.",
"Rho"
10863 WRITE (unit_nr,
'(T2,A6,A5,I6,F22.10,ES10.2,T67,ES7.0,F6.1)') &
10874 WRITE (unit_nr,
'(T2,A6,A5,I6,F22.10,ES10.2,ES10.2,F6.1,ES7.0,F6.1)') &
10879 delta_loss, predicted_reduction, rho, &
10886 END SUBROUTINE trust_r_report
10894 SUBROUTINE energy_lowering_report(unit_nr, ref_energy, energy_lowering)
10896 INTEGER,
INTENT(IN) :: unit_nr
10897 REAL(kind=
dp),
INTENT(IN) :: ref_energy, energy_lowering
10900 IF (unit_nr > 0)
THEN
10902 WRITE (unit_nr,
'(T2,A35,F25.10)')
"ENERGY OF BLOCK-DIAGONAL ALMOs:", &
10904 WRITE (unit_nr,
'(T2,A35,F25.10)')
"ENERGY LOWERING:", &
10906 WRITE (unit_nr,
'(T2,A35,F25.10)')
"CORRECTED ENERGY:", &
10907 ref_energy + energy_lowering
10911 END SUBROUTINE energy_lowering_report
10923 SUBROUTINE wrap_up_xalmo_scf(qs_env, almo_scf_env, perturbation_in, &
10924 m_xalmo_in, m_quench_in, energy_inout)
10928 LOGICAL,
INTENT(IN) :: perturbation_in
10929 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(IN) :: m_xalmo_in, m_quench_in
10930 REAL(kind=
dp),
INTENT(INOUT) :: energy_inout
10932 CHARACTER(len=*),
PARAMETER :: routinen =
'wrap_up_xalmo_scf'
10934 INTEGER :: eda_unit, handle, ispin, nspins, unit_nr
10936 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:) :: m_temp_no1, m_temp_no2
10939 CALL timeset(routinen, handle)
10943 IF (logger%para_env%is_source())
THEN
10949 nspins = almo_scf_env%nspins
10953 IF (perturbation_in)
THEN
10955 ALLOCATE (m_temp_no1(nspins))
10956 ALLOCATE (m_temp_no2(nspins))
10958 DO ispin = 1, nspins
10959 CALL dbcsr_create(m_temp_no1(ispin), template=m_xalmo_in(ispin))
10960 CALL dbcsr_create(m_temp_no2(ispin), template=m_xalmo_in(ispin))
10965 almo_scf_env%mat_distr_aos)
10970 CALL xalmo_analysis( &
10971 detailed_analysis=almo_scf_env%almo_analysis%do_analysis, &
10972 eps_filter=almo_scf_env%eps_filter, &
10973 m_t_in=m_xalmo_in, &
10974 m_t0_in=almo_scf_env%matrix_t_blk, &
10975 m_siginv_in=almo_scf_env%matrix_sigma_inv, &
10976 m_siginv0_in=almo_scf_env%matrix_sigma_inv_0deloc, &
10977 m_s_in=almo_scf_env%matrix_s, &
10978 m_ks0_in=almo_scf_env%matrix_ks_0deloc, &
10979 m_quench_t_in=m_quench_in, &
10980 energy_out=energy_inout, &
10981 m_eda_out=m_temp_no1, &
10982 m_cta_out=m_temp_no2 &
10985 IF (almo_scf_env%almo_analysis%do_analysis)
THEN
10987 DO ispin = 1, nspins
10990 IF (unit_nr > 0)
THEN
10991 WRITE (unit_nr,
'(T2,A)')
"DECOMPOSITION OF THE DELOCALIZATION ENERGY"
10998 "ALMO_EDA_CT", extension=
".dat", local=.true.)
10999 CALL print_block_sum(m_temp_no1(ispin), eda_unit)
11001 "ALMO_EDA_CT", local=.true.)
11004 IF (unit_nr > 0)
THEN
11005 WRITE (unit_nr,
'(T2,A)')
"DECOMPOSITION OF CHARGE TRANSFER TERMS"
11009 "ALMO_CTA", extension=
".dat", local=.true.)
11010 CALL print_block_sum(m_temp_no2(ispin), eda_unit)
11012 "ALMO_CTA", local=.true.)
11018 CALL energy_lowering_report( &
11020 ref_energy=almo_scf_env%almo_scf_energy, &
11021 energy_lowering=energy_inout)
11023 energy=almo_scf_env%almo_scf_energy, &
11024 energy_singles_corr=energy_inout)
11026 DO ispin = 1, nspins
11031 DEALLOCATE (m_temp_no1)
11032 DEALLOCATE (m_temp_no2)
11037 energy=energy_inout)
11041 CALL timestop(handle)
11043 END SUBROUTINE wrap_up_xalmo_scf
11051 SUBROUTINE tanh_of_elements(matrix, alpha)
11053 REAL(kind=
dp),
INTENT(IN) :: alpha
11055 CHARACTER(len=*),
PARAMETER :: routinen =
'tanh_of_elements'
11058 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: block
11061 CALL timeset(routinen, handle)
11065 block = tanh(alpha*block)
11068 CALL timestop(handle)
11070 END SUBROUTINE tanh_of_elements
11078 SUBROUTINE dtanh_of_elements(matrix, alpha)
11080 REAL(kind=
dp),
INTENT(IN) :: alpha
11082 CHARACTER(len=*),
PARAMETER :: routinen =
'dtanh_of_elements'
11085 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: block
11088 CALL timeset(routinen, handle)
11092 block = alpha*(1.0_dp - tanh(block)**2)
11095 CALL timestop(handle)
11097 END SUBROUTINE dtanh_of_elements
11104 SUBROUTINE inverse_of_elements(matrix)
11107 CHARACTER(len=*),
PARAMETER :: routinen =
'inverse_of_elements'
11110 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: block
11113 CALL timeset(routinen, handle)
11117 block = 1.0_dp/block
11120 CALL timestop(handle)
11122 END SUBROUTINE inverse_of_elements
11129 SUBROUTINE print_block_sum(matrix, unit_nr)
11131 INTEGER,
INTENT(IN) :: unit_nr
11133 CHARACTER(len=*),
PARAMETER :: routinen =
'print_block_sum'
11135 INTEGER :: col, handle, row
11136 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: block
11139 CALL timeset(routinen, handle)
11141 IF (unit_nr > 0)
THEN
11145 WRITE (unit_nr,
'(I6,I6,ES18.9)') row, col, sum(block)
11150 CALL timestop(handle)
11151 END SUBROUTINE print_block_sum
A DIIS implementation for the ALMO-based SCF methods.
subroutine, public almo_scf_diis_release(diis_env)
destroys the diis structure
subroutine, public almo_scf_diis_extrapolate(diis_env, extr_var, d_extr_var)
extrapolates the variable using the saved history
subroutine, public almo_scf_diis_push(diis_env, var, err, d_var, d_err)
adds a variable-error pair to the diis structure
subroutine, public lbfgs_create(history, nspins, nstore)
create history storage for limited memory bfgs
subroutine, public lbfgs_seed(history, variable, gradient)
interface subroutine to store the first variable/gradient pair
subroutine, public lbfgs_release(history)
release the bfgs history
subroutine, public lbfgs_get_direction(history, variable, gradient, direction)
interface subroutine to store a variable/gradient pair and predict direction
Subroutines for ALMO SCF.
subroutine, public construct_domain_preconditioner(matrix_main, subm_s_inv, subm_s_inv_half, subm_s_half, subm_r_down, matrix_trimmer, dpattern, map, node_of_domain, preconditioner, bad_modes_projector_down, use_trimmer, eps_zero_eigenvalues, my_action, skip_inversion)
Constructs preconditioners for each domain -1. projected preconditioner 0. simple preconditioner.
subroutine, public almo_scf_ks_xx_to_tv_xx(almo_scf_env)
ALMOs by diagonalizing the KS domain submatrices computes both the occupied and virtual orbitals.
subroutine, public xalmo_initial_guess(m_guess, m_t_in, m_t0, m_quench_t, m_overlap, m_sigma_tmpl, nspins, xalmo_history, assume_t0_q0x, optimize_theta, envelope_amplitude, eps_filter, order_lanczos, eps_lanczos, max_iter_lanczos, nocc_of_domain)
create the initial guess for XALMOs
subroutine, public almo_scf_p_blk_to_t_blk(almo_scf_env, ionic)
computes occupied ALMOs from the superimposed atomic density blocks
subroutine, public pseudo_invert_diagonal_blk(matrix_in, matrix_out, nocc)
inverts block-diagonal blocks of a dbcsr_matrix
subroutine, public almo_scf_ks_blk_to_tv_blk(almo_scf_env)
computes ALMOs by diagonalizing the projected blocked KS matrix uses the diagonalization code for blo...
subroutine, public apply_domain_operators(matrix_in, matrix_out, operator1, operator2, dpattern, map, node_of_domain, my_action, filter_eps, matrix_trimmer, use_trimmer)
Parallel code for domain specific operations (my_action) 0. out = op1 * in.
subroutine, public construct_domain_r_down(matrix_t, matrix_sigma_inv, matrix_s, subm_r_down, dpattern, map, node_of_domain, filter_eps)
Constructs subblocks of the covariant-covariant projectors (i.e. DM without spin factor)
subroutine, public almo_scf_t_to_proj(t, p, eps_filter, orthog_orbs, nocc_of_domain, s, sigma, sigma_inv, use_guess, smear, algorithm, para_env, blacs_env, eps_lanczos, max_iter_lanczos, inverse_accelerator, inv_eps_factor)
computes the idempotent density matrix from MOs MOs can be either orthogonal or non-orthogonal
subroutine, public construct_domain_s_inv(matrix_s, subm_s_inv, dpattern, map, node_of_domain)
Constructs S_inv block for each domain.
subroutine, public almo_scf_ks_to_ks_blk(almo_scf_env)
computes the projected KS from the total KS matrix also computes the DIIS error vector as a by-produc...
subroutine, public get_overlap(bra, ket, overlap, metric, retain_overlap_sparsity, eps_filter, smear)
Computes the overlap matrix of MO orbitals.
subroutine, public fill_matrix_with_ones(matrix)
Fill all matrix blocks with 1.0_dp.
subroutine, public apply_projector(psi_in, psi_out, psi_projector, metric, project_out, psi_projector_orthogonal, proj_in_template, eps_filter, sig_inv_projector, sig_inv_template)
applies projector to the orbitals |psi_out> = P |psi_in> OR |psi_out> = (1-P) |psi_in>,...
subroutine, public construct_domain_s_sqrt(matrix_s, subm_s_sqrt, subm_s_sqrt_inv, dpattern, map, node_of_domain)
Constructs S^(+1/2) and S^(-1/2) submatrices for each domain.
subroutine, public orthogonalize_mos(ket, overlap, metric, retain_locality, only_normalize, nocc_of_domain, eps_filter, order_lanczos, eps_lanczos, max_iter_lanczos, overlap_sqrti, smear)
orthogonalize MOs
subroutine, public almo_scf_ks_to_ks_xx(almo_scf_env)
builds projected KS matrices for the overlapping domains also computes the DIIS error vector as a by-...
subroutine, public almo_scf_t_rescaling(matrix_t, mo_energies, mu_of_domain, real_ne_of_domain, spin_kts, smear_e_temp, ndomains, nocc_of_domain)
Apply an occupation-rescaling trick to ALMOs for smearing. Partially occupied orbitals are considered...
Optimization routines for all ALMO-based SCF methods.
subroutine, public almo_scf_xalmo_trustr(qs_env, almo_scf_env, optimizer, quench_t, matrix_t_in, matrix_t_out, perturbation_only, special_case)
Optimization of ALMOs using trust region minimizers.
subroutine, public almo_scf_xalmo_pcg(qs_env, almo_scf_env, optimizer, quench_t, matrix_t_in, matrix_t_out, assume_t0_q0x, perturbation_only, special_case)
Optimization of ALMOs using PCG-like minimizers.
subroutine, public almo_scf_xalmo_eigensolver(qs_env, almo_scf_env, optimizer)
An eigensolver-based SCF to optimize extended ALMOs (i.e. ALMOs on overlapping domains)
subroutine, public almo_scf_construct_nlmos(qs_env, optimizer, matrix_s, matrix_mo_in, matrix_mo_out, template_matrix_sigma, overlap_determinant, mat_distr_aos, virtuals, eps_filter)
Optimization of NLMOs using PCG minimizers.
subroutine, public almo_scf_block_diagonal(qs_env, almo_scf_env, optimizer)
An SCF procedure that optimizes block-diagonal ALMOs using DIIS.
Interface between ALMO SCF and QS.
subroutine, public almo_scf_update_ks_energy(qs_env, energy, energy_singles_corr)
update qs_env total energy
subroutine, public almo_dm_to_almo_ks(qs_env, matrix_p, matrix_ks, energy_total, eps_filter, mat_distr_aos, smear, kts_sum)
uses the ALMO density matrix to compute ALMO KS matrix and the new energy
subroutine, public almo_dm_to_qs_env(qs_env, matrix_p, mat_distr_aos)
return density matrix to the qs_env
subroutine, public matrix_qs_to_almo(matrix_qs, matrix_almo, mat_distr_aos)
convert between two types of matrices: QS style to ALMO style
Types for all ALMO-based methods.
Handles all functions related to the CELL.
methods related to the blacs parallel environment
subroutine, public dbcsr_scale(matrix, alpha_scalar)
...
subroutine, public dbcsr_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset)
...
logical function, public dbcsr_iterator_blocks_left(iterator)
...
subroutine, public dbcsr_iterator_stop(iterator)
...
subroutine, public dbcsr_desymmetrize(matrix_a, matrix_b)
...
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
subroutine, public dbcsr_get_block_p(matrix, row, col, block, found, row_size, col_size)
...
subroutine, public dbcsr_multiply(transa, transb, alpha, matrix_a, matrix_b, beta, matrix_c, first_row, last_row, first_column, last_column, first_k, last_k, retain_sparsity, filter_eps, flop)
...
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
subroutine, public dbcsr_work_create(matrix, nblks_guess, sizedata_guess, n, work_mutable)
...
subroutine, public dbcsr_filter(matrix, eps)
...
subroutine, public dbcsr_finalize(matrix)
...
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
subroutine, public dbcsr_set(matrix, alpha)
...
subroutine, public dbcsr_release(matrix)
...
subroutine, public dbcsr_iterator_readonly_start(iterator, matrix, shared, dynamic, dynamic_byrows)
Like dbcsr_iterator_start() but with matrix being INTENT(IN). When invoking this routine,...
subroutine, public dbcsr_put_block(matrix, row, col, block, summation)
...
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
subroutine, public dbcsr_distribution_get(dist, row_dist, col_dist, nrows, ncols, has_threads, group, mynode, numnodes, nprows, npcols, myprow, mypcol, pgrid, subgroups_defined, prow_group, pcol_group)
...
Interface to (sca)lapack for the Cholesky based procedures.
subroutine, public cp_dbcsr_cholesky_decompose(matrix, n, para_env, blacs_env)
used to replace a symmetric positive def. matrix M with its cholesky decomposition U: M = U^T * U,...
subroutine, public cp_dbcsr_cholesky_restore(matrix, neig, matrixb, matrixout, op, pos, transa, para_env, blacs_env)
...
subroutine, public cp_dbcsr_cholesky_invert(matrix, n, para_env, blacs_env, uplo_to_full)
used to replace the cholesky decomposition by the inverse
subroutine, public dbcsr_set_diag(matrix, diag)
Copies the diagonal elements from the given array into the given matrix.
subroutine, public dbcsr_get_diag(matrix, diag)
Copies the diagonal elements from the given matrix into the given array.
subroutine, public dbcsr_add_on_diag(matrix, alpha)
Adds the given scalar to the diagonal of the matrix. Reserves any missing diagonal blocks.
real(dp) function, public dbcsr_maxabs(matrix)
Compute the maxabs norm of a dbcsr matrix.
real(dp) function, public dbcsr_frobenius_norm(matrix)
Compute the frobenius norm of a dbcsr matrix.
subroutine, public dbcsr_dot(matrix_a, matrix_b, trace)
Computes the dot product of two matrices, also known as the trace of their matrix product.
subroutine, public dbcsr_hadamard_product(matrix_a, matrix_b, matrix_c)
Hadamard product: C = A . B (C needs to be different from A and B)
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.
various routines to log and control the output. The idea is that decisions about where to log should ...
recursive integer function, public cp_logger_get_default_unit_nr(logger, local, skip_not_ionode)
asks the default unit number of the given logger. try to use cp_logger_get_unit_nr
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,...
Cayley transformation methods.
subroutine, public analytic_line_search(a, b, c, d, minima, nmins)
Finds real roots of a cubic equation a*x**3 + b*x**2 + c*x + d = 0 and returns only those roots for...
subroutine, public diagonalize_diagonal_blocks(matrix, c, e)
Diagonalizes diagonal blocks of a symmetric dbcsr matrix and returs its eigenvectors.
subroutine, public ct_step_execute(cts_env)
Performs Cayley transformation.
Types for all cayley transformation methods.
subroutine, public ct_step_env_clean(env)
...
subroutine, public ct_step_env_set(env, para_env, blacs_env, use_occ_orbs, use_virt_orbs, tensor_type, occ_orbs_orthogonal, virt_orbs_orthogonal, neglect_quadratic_term, update_p, update_q, eps_convergence, eps_filter, max_iter, p_index_up, p_index_down, q_index_up, q_index_down, matrix_ks, matrix_p, matrix_qp_template, matrix_pq_template, matrix_t, matrix_v, matrix_x_guess, calculate_energy_corr, conjugator, qq_preconditioner_full, pp_preconditioner_full)
...
subroutine, public ct_step_env_init(env)
...
subroutine, public ct_step_env_get(env, use_occ_orbs, use_virt_orbs, tensor_type, occ_orbs_orthogonal, virt_orbs_orthogonal, neglect_quadratic_term, update_p, update_q, eps_convergence, eps_filter, max_iter, p_index_up, p_index_down, q_index_up, q_index_down, matrix_ks, matrix_p, matrix_qp_template, matrix_pq_template, matrix_t, matrix_v, copy_matrix_x, energy_correction, calculate_energy_corr, converged, qq_preconditioner_full, pp_preconditioner_full)
...
Subroutines to handle submatrices.
subroutine, public maxnorm_submatrices(submatrices, norm)
Computes the max norm of the collection of submatrices.
subroutine, public construct_submatrices(matrix, submatrix, distr_pattern, domain_map, node_of_domain, job_type)
Constructs submatrices for each ALMO domain by collecting distributed DBCSR blocks to local arrays.
Types to handle submatrices.
integer, parameter, public select_row
Routines useful for iterative matrix calculations.
recursive subroutine, public determinant(matrix, det, threshold)
Computes the determinant of a symmetric positive definite matrix using the trace of the matrix logari...
subroutine, public invert_hotelling(matrix_inverse, matrix, threshold, use_inv_as_guess, norm_convergence, filter_eps, accelerator_order, max_iter_lanczos, eps_lanczos, silent)
invert a symmetric positive definite matrix by Hotelling's method explicit symmetrization makes this ...
subroutine, public matrix_sqrt_newton_schulz(matrix_sqrt, matrix_sqrt_inv, matrix, threshold, order, eps_lanczos, max_iter_lanczos, symmetrize, converged, iounit)
compute the sqrt of a matrix via the sign function and the corresponding Newton-Schulz iterations the...
Defines the basic variable types.
integer, parameter, public dp
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.
Define methods related to particle_type.
subroutine, public get_particle_set(particle_set, qs_kind_set, first_sgf, last_sgf, nsgf, nmao, basis)
Get the components of a particle set.
Define the data structure for the particle information.
computes preconditioners, and implements methods to apply them currently used in qs_ot
Perform a QUICKSTEP wavefunction optimization (single point)
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs)
Get the QUICKSTEP environment.
Define the quickstep kind type and their sub types.
Some utilities for the construction of the localization environment.
subroutine, public compute_berry_operator(qs_env, cell, op_sm_set, dim_op)
Computes the Berry operator for periodic systems used to define the spread of the MOS Here the matrix...
Localization methods such as 2x2 Jacobi rotations Steepest Decents Conjugate Gradient.
subroutine, public initialize_weights(cell, weights)
...
Type defining parameters related to the simulation cell.
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
type of a logger, at the moment it contains just a print level starting at which level it should be l...
stores all the informations relevant to an mpi environment
Provides all information about a quickstep kind.