107#include "./base/base_uses.f90"
113 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'almo_scf_optimizer'
121 LOGICAL,
PARAMETER :: debug_mode = .false.
122 LOGICAL,
PARAMETER :: safe_mode = .false.
123 LOGICAL,
PARAMETER :: almo_mathematica = .false.
124 INTEGER,
PARAMETER :: hessian_path_reuse = 1, &
125 hessian_path_assemble = 2
144 CHARACTER(len=*),
PARAMETER :: routinen =
'almo_scf_block_diagonal'
146 INTEGER :: handle, iscf, ispin, nspin, unit_nr
147 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: local_nocc_of_domain
148 LOGICAL :: converged, prepare_to_exit, should_stop, &
149 use_diis, use_prev_as_guess
150 REAL(kind=
dp) :: density_rec, energy_diff, energy_new, energy_old, error_norm, &
151 error_norm_ispin, kts_sum, prev_error_norm, t1, t2, true_mixing_fraction
152 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: local_mu
154 DIMENSION(:) :: almo_diis
156 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:) :: matrix_mixing_old_blk
159 CALL timeset(routinen, handle)
163 IF (logger%para_env%is_source())
THEN
171 use_prev_as_guess = .false.
173 nspin = almo_scf_env%nspins
174 ALLOCATE (local_mu(almo_scf_env%ndomains))
175 ALLOCATE (local_nocc_of_domain(almo_scf_env%ndomains))
178 ALLOCATE (matrix_mixing_old_blk(nspin))
179 ALLOCATE (almo_diis(nspin))
182 template=almo_scf_env%matrix_ks_blk(ispin))
184 sample_err=almo_scf_env%matrix_ks_blk(ispin), &
185 sample_var=almo_scf_env%matrix_s_blk(1), &
187 max_length=optimizer%ndiis)
194 prepare_to_exit = .false.
195 true_mixing_fraction = 0.0_dp
196 error_norm = 1.0e+10_dp
198 IF (unit_nr > 0)
THEN
199 WRITE (unit_nr,
'(T2,A,A,A)') repeat(
"-", 20), &
200 " Optimization of block-diagonal ALMOs ", repeat(
"-", 21)
202 WRITE (unit_nr,
'(T2,A13,A6,A23,A14,A14,A9)')
"Method",
"Iter", &
203 "Total Energy",
"Change",
"Convergence",
"Time"
204 WRITE (unit_nr,
'(T2,A)') repeat(
"-", 79)
220 var=almo_scf_env%matrix_ks_blk(ispin), &
221 err=almo_scf_env%matrix_err_blk(ispin))
226 prev_error_norm = error_norm
228 error_norm_ispin =
dbcsr_maxabs(almo_scf_env%matrix_err_blk(ispin))
229 IF (ispin == 1) error_norm = error_norm_ispin
230 IF (ispin > 1 .AND. error_norm_ispin > error_norm)
THEN
231 error_norm = error_norm_ispin
235 IF (error_norm < almo_scf_env%eps_prev_guess)
THEN
236 use_prev_as_guess = .true.
238 use_prev_as_guess = .false.
243 IF (error_norm > optimizer%eps_error) converged = .false.
247 start_time=qs_env%start_time, &
248 target_time=qs_env%target_time)
249 IF (should_stop .OR. iscf >= optimizer%max_iter .OR. converged)
THEN
250 prepare_to_exit = .true.
251 IF (iscf == 1) energy_new = energy_old
255 IF (optimizer%early_stopping_on .AND. iscf == 1)
THEN
256 prepare_to_exit = .false.
259 IF (.NOT. prepare_to_exit)
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)
297 cpabort(
"Density_matrix_sign has not been tested yet")
298 almo_scf_env%mu_of_domain(:, ispin) = local_mu(:)
305 DO ispin = 1, almo_scf_env%nspins
308 overlap=almo_scf_env%matrix_sigma_blk(ispin), &
309 metric=almo_scf_env%matrix_s_blk(1), &
310 retain_locality=.true., &
311 only_normalize=.false., &
312 nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin), &
313 eps_filter=almo_scf_env%eps_filter, &
314 order_lanczos=almo_scf_env%order_lanczos, &
315 eps_lanczos=almo_scf_env%eps_lanczos, &
316 max_iter_lanczos=almo_scf_env%max_iter_lanczos)
323 DO ispin = 1, almo_scf_env%nspins
326 IF (almo_scf_env%smear)
THEN
328 mo_energies=almo_scf_env%mo_energies(:, ispin), &
329 mu_of_domain=almo_scf_env%mu_of_domain(:, ispin), &
330 real_ne_of_domain=almo_scf_env%real_ne_of_domain(:, ispin), &
331 spin_kts=almo_scf_env%kTS(ispin), &
332 smear_e_temp=almo_scf_env%smear_e_temp, &
333 ndomains=almo_scf_env%ndomains, &
334 nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin))
338 p=almo_scf_env%matrix_p(ispin), &
339 eps_filter=almo_scf_env%eps_filter, &
340 orthog_orbs=.false., &
341 nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin), &
342 s=almo_scf_env%matrix_s(1), &
343 sigma=almo_scf_env%matrix_sigma(ispin), &
344 sigma_inv=almo_scf_env%matrix_sigma_inv(ispin), &
345 use_guess=use_prev_as_guess, &
346 smear=almo_scf_env%smear, &
347 algorithm=almo_scf_env%sigma_inv_algorithm, &
348 inverse_accelerator=almo_scf_env%order_lanczos, &
349 inv_eps_factor=almo_scf_env%matrix_iter_eps_error_factor, &
350 eps_lanczos=almo_scf_env%eps_lanczos, &
351 max_iter_lanczos=almo_scf_env%max_iter_lanczos, &
352 para_env=almo_scf_env%para_env, &
353 blacs_env=almo_scf_env%blacs_env)
357 IF (almo_scf_env%nspins == 1)
THEN
360 IF (almo_scf_env%smear)
THEN
361 almo_scf_env%kTS(1) = almo_scf_env%kTS(1)*2.0_dp
365 IF (almo_scf_env%smear)
THEN
366 kts_sum = sum(almo_scf_env%kTS)
373 almo_scf_env%matrix_p, &
374 almo_scf_env%matrix_ks, &
376 almo_scf_env%eps_filter, &
377 almo_scf_env%mat_distr_aos, &
378 smear=almo_scf_env%smear, &
383 energy_diff = energy_new - energy_old
384 energy_old = energy_new
385 almo_scf_env%almo_scf_energy = energy_new
389 IF (unit_nr > 0)
THEN
390 WRITE (unit_nr,
'(T2,A13,I6,F23.10,E14.5,F14.9,F9.2)')
"ALMO SCF DIIS", &
392 energy_new, energy_diff, error_norm, t2 - t1
396 IF (prepare_to_exit)
EXIT
401 IF (almo_scf_env%smear)
THEN
403 CALL dbcsr_dot(almo_scf_env%matrix_p(ispin), almo_scf_env%matrix_s(1), density_rec)
404 IF (unit_nr > 0)
THEN
405 WRITE (unit_nr,
'(T2,A20,F23.10)')
"Electrons recovered:", density_rec
410 IF (.NOT. converged .AND. (.NOT. optimizer%early_stopping_on))
THEN
411 IF (unit_nr > 0)
THEN
412 cpabort(
"SCF for block-diagonal ALMOs not converged!")
420 DEALLOCATE (almo_diis)
421 DEALLOCATE (matrix_mixing_old_blk)
422 DEALLOCATE (local_mu)
423 DEALLOCATE (local_nocc_of_domain)
425 CALL timestop(handle)
445 CHARACTER(len=*),
PARAMETER :: routinen =
'almo_scf_xalmo_eigensolver'
447 INTEGER :: handle, iscf, ispin, nspin, unit_nr
448 LOGICAL :: converged, prepare_to_exit, should_stop
449 REAL(kind=
dp) :: denergy_tot, density_rec, energy_diff, energy_new, energy_old, error_norm, &
450 error_norm_0, kts_sum, spin_factor, t1, t2
451 REAL(kind=
dp),
DIMENSION(2) :: denergy_spin
453 DIMENSION(:) :: almo_diis
455 TYPE(
dbcsr_type) :: matrix_p_almo_scf_converged
457 DIMENSION(:, :) :: submatrix_mixing_old_blk
459 CALL timeset(routinen, handle)
463 IF (logger%para_env%is_source())
THEN
469 nspin = almo_scf_env%nspins
480 matrix_s=almo_scf_env%matrix_s(1), &
481 subm_s_sqrt=almo_scf_env%domain_s_sqrt(:, ispin), &
482 subm_s_sqrt_inv=almo_scf_env%domain_s_sqrt_inv(:, ispin), &
483 dpattern=almo_scf_env%quench_t(ispin), &
484 map=almo_scf_env%domain_map(ispin), &
485 node_of_domain=almo_scf_env%cpu_of_domain)
493 matrix=almo_scf_env%quench_t(ispin), &
494 submatrix=almo_scf_env%domain_t(:, ispin), &
495 distr_pattern=almo_scf_env%quench_t(ispin), &
496 domain_map=almo_scf_env%domain_map(ispin), &
497 node_of_domain=almo_scf_env%cpu_of_domain, &
502 ALLOCATE (submatrix_mixing_old_blk(almo_scf_env%ndomains, nspin))
504 ALLOCATE (almo_diis(nspin))
510 sample_err=almo_scf_env%domain_s_sqrt(:, ispin), &
512 max_length=optimizer%ndiis)
518 prepare_to_exit = .false.
532 d_var=almo_scf_env%domain_ks_xx(:, ispin), &
533 d_err=almo_scf_env%domain_err(:, ispin))
539 error_norm =
dbcsr_maxabs(almo_scf_env%matrix_err_xx(ispin))
542 IF (error_norm > optimizer%eps_error)
THEN
549 start_time=qs_env%start_time, &
550 target_time=qs_env%target_time)
551 IF (should_stop .OR. iscf >= optimizer%max_iter .OR. converged)
THEN
552 prepare_to_exit = .true.
556 IF (optimizer%early_stopping_on .AND. iscf == 1)
THEN
557 prepare_to_exit = .false.
560 IF (.NOT. prepare_to_exit)
THEN
566 d_extr_var=almo_scf_env%domain_ks_xx(:, ispin))
572 almo_scf_env%domain_ks_xx(:, ispin), &
573 submatrix_mixing_old_blk(:, ispin), &
586 template=almo_scf_env%matrix_p(ispin))
587 CALL dbcsr_copy(matrix_p_almo_scf_converged, &
588 almo_scf_env%matrix_p(ispin))
592 IF (almo_scf_env%smear)
THEN
594 mo_energies=almo_scf_env%mo_energies(:, ispin), &
595 mu_of_domain=almo_scf_env%mu_of_domain(:, ispin), &
596 real_ne_of_domain=almo_scf_env%real_ne_of_domain(:, ispin), &
597 spin_kts=almo_scf_env%kTS(ispin), &
598 smear_e_temp=almo_scf_env%smear_e_temp, &
599 ndomains=almo_scf_env%ndomains, &
600 nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin))
605 t=almo_scf_env%matrix_t(ispin), &
606 p=almo_scf_env%matrix_p(ispin), &
607 eps_filter=almo_scf_env%eps_filter, &
608 orthog_orbs=.false., &
609 nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin), &
610 s=almo_scf_env%matrix_s(1), &
611 sigma=almo_scf_env%matrix_sigma(ispin), &
612 sigma_inv=almo_scf_env%matrix_sigma_inv(ispin), &
614 smear=almo_scf_env%smear, &
615 algorithm=almo_scf_env%sigma_inv_algorithm, &
616 inverse_accelerator=almo_scf_env%order_lanczos, &
617 inv_eps_factor=almo_scf_env%matrix_iter_eps_error_factor, &
618 eps_lanczos=almo_scf_env%eps_lanczos, &
619 max_iter_lanczos=almo_scf_env%max_iter_lanczos, &
620 para_env=almo_scf_env%para_env, &
621 blacs_env=almo_scf_env%blacs_env)
622 CALL dbcsr_scale(almo_scf_env%matrix_p(ispin), spin_factor)
624 IF (almo_scf_env%smear)
THEN
625 almo_scf_env%kTS(ispin) = almo_scf_env%kTS(ispin)*spin_factor
632 CALL dbcsr_add(matrix_p_almo_scf_converged, &
633 almo_scf_env%matrix_p(ispin), -1.0_dp, 1.0_dp)
634 CALL dbcsr_dot(almo_scf_env%matrix_ks_0deloc(ispin), &
635 matrix_p_almo_scf_converged, &
642 denergy_tot = denergy_tot + denergy_spin(ispin)
650 CALL energy_lowering_report( &
652 ref_energy=almo_scf_env%almo_scf_energy, &
653 energy_lowering=denergy_tot)
655 energy=almo_scf_env%almo_scf_energy, &
656 energy_singles_corr=denergy_tot)
660 IF (.NOT. almo_scf_env%perturbative_delocalization)
THEN
662 IF (almo_scf_env%smear)
THEN
663 kts_sum = sum(almo_scf_env%kTS)
669 almo_scf_env%matrix_p, &
670 almo_scf_env%matrix_ks, &
672 almo_scf_env%eps_filter, &
673 almo_scf_env%mat_distr_aos, &
674 smear=almo_scf_env%smear, &
680 IF (almo_scf_env%perturbative_delocalization)
THEN
683 CALL almo_dm_to_qs_env(qs_env, almo_scf_env%matrix_p, almo_scf_env%mat_distr_aos)
685 prepare_to_exit = .true.
689 energy_diff = energy_new - energy_old
690 energy_old = energy_new
691 almo_scf_env%almo_scf_energy = energy_new
695 IF (unit_nr > 0)
THEN
696 WRITE (unit_nr,
'(T2,A,I6,F20.9,E11.3,E11.3,E11.3,F8.2)')
"ALMO SCF", &
698 energy_new, energy_diff, error_norm, error_norm_0, t2 - t1
704 IF (prepare_to_exit)
EXIT
709 IF (almo_scf_env%smear)
THEN
711 CALL dbcsr_dot(almo_scf_env%matrix_p(ispin), almo_scf_env%matrix_s(1), density_rec)
712 IF (unit_nr > 0)
THEN
713 WRITE (unit_nr,
'(T2,A20,F23.10)')
"Electrons recovered:", density_rec
718 IF (.NOT. converged .AND. .NOT. optimizer%early_stopping_on)
THEN
719 cpabort(
"SCF for ALMOs on overlapping domains not converged!")
726 DEALLOCATE (almo_diis)
727 DEALLOCATE (submatrix_mixing_old_blk)
729 CALL timestop(handle)
755 matrix_t_in, matrix_t_out, assume_t0_q0x, perturbation_only, &
761 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:), &
762 INTENT(INOUT) :: quench_t, matrix_t_in, matrix_t_out
763 LOGICAL,
INTENT(IN) :: assume_t0_q0x, perturbation_only
764 INTEGER,
INTENT(IN),
OPTIONAL :: special_case
766 CHARACTER(len=*),
PARAMETER :: routinen =
'almo_scf_xalmo_pcg'
768 CHARACTER(LEN=20) :: iter_type
769 INTEGER :: cg_iteration, dim_op, fixed_line_search_niter, handle, idim0, ielem, ispin, &
770 iteration, line_search_iteration, max_iter, my_special_case, ndomains, nmo, nspins, &
771 outer_iteration, outer_max_iter, prec_type, reim, unit_nr
772 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: nocc
773 LOGICAL :: blissful_neglect, converged, just_started, line_search, normalize_orbitals, &
774 optimize_theta, outer_prepare_to_exit, penalty_occ_local, penalty_occ_vol, &
775 prepare_to_exit, reset_conjugator, skip_grad, use_guess
776 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: reim_diag, weights, z2
777 REAL(kind=
dp) :: appr_sec_der, beta, denom, denom2, e0, e1, energy_coeff, energy_diff, &
778 energy_new, energy_old, eps_skip_gradients, fval, g0, g1, grad_norm, grad_norm_frob, &
779 line_search_error, localiz_coeff, localization_obj_function, next_step_size_guess, &
780 penalty_amplitude, penalty_func_new, spin_factor, step_size, t1, t2, tempreal
781 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: grad_norm_spin, &
782 penalty_occ_vol_g_prefactor, &
783 penalty_occ_vol_h_prefactor
786 TYPE(
dbcsr_p_type),
DIMENSION(:),
POINTER :: qs_matrix_s
787 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
POINTER :: op_sm_set_almo, op_sm_set_qs
788 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:) :: ftsiginv, grad, m_sig_sqrti_ii, m_t_in_local, &
789 m_theta, prec_vv, prev_grad, prev_minus_prec_grad, prev_step, siginvtftsiginv, st, step, &
790 stsiginv_0, tempnocc, tempnocc_1, tempoccocc
792 DIMENSION(:, :) :: bad_modes_projector_down, domain_r_down
795 CALL timeset(routinen, handle)
798 IF (
PRESENT(special_case)) my_special_case = special_case
802 IF (logger%para_env%is_source())
THEN
808 nspins = almo_scf_env%nspins
812 blissful_neglect = .false.
814 blissful_neglect = .true.
817 IF (unit_nr > 0)
THEN
819 SELECT CASE (my_special_case)
821 WRITE (unit_nr,
'(T2,A,A,A)') repeat(
"-", 20), &
822 " Optimization of block-diagonal ALMOs ", repeat(
"-", 21)
824 WRITE (unit_nr,
'(T2,A,A,A)') repeat(
"-", 20), &
825 " Optimization of fully delocalized MOs ", repeat(
"-", 20)
827 IF (blissful_neglect)
THEN
828 WRITE (unit_nr,
'(T2,A,A,A)') repeat(
"-", 25), &
829 " LCP optimization of XALMOs ", repeat(
"-", 26)
831 WRITE (unit_nr,
'(T2,A,A,A)') repeat(
"-", 27), &
832 " Optimization of XALMOs ", repeat(
"-", 28)
836 WRITE (unit_nr,
'(T2,A13,A6,A23,A14,A14,A9)')
"Method",
"Iter", &
837 "Objective Function",
"Change",
"Convergence",
"Time"
838 WRITE (unit_nr,
'(T2,A)') repeat(
"-", 79)
843 optimize_theta = almo_scf_env%logical05
844 eps_skip_gradients = almo_scf_env%real01
847 energy_coeff = 1.0_dp
848 localiz_coeff = 0.0_dp
849 penalty_amplitude = 0.0_dp
850 penalty_occ_vol = .false.
851 penalty_occ_local = .false.
852 normalize_orbitals = penalty_occ_vol .OR. penalty_occ_local
853 ALLOCATE (penalty_occ_vol_g_prefactor(nspins))
854 ALLOCATE (penalty_occ_vol_h_prefactor(nspins))
855 penalty_occ_vol_g_prefactor(:) = 0.0_dp
856 penalty_occ_vol_h_prefactor(:) = 0.0_dp
857 penalty_func_new = 0.0_dp
860 prec_type = optimizer%preconditioner
863 fixed_line_search_niter = 0
865 IF (nspins == 1)
THEN
871 ALLOCATE (grad_norm_spin(nspins))
872 ALLOCATE (nocc(nspins))
878 ALLOCATE (m_t_in_local(nspins))
881 template=matrix_t_in(ispin), &
882 matrix_type=dbcsr_type_no_symmetry)
883 CALL dbcsr_copy(m_t_in_local(ispin), matrix_t_in(ispin))
888 ALLOCATE (m_theta(nspins))
891 template=matrix_t_out(ispin), &
892 matrix_type=dbcsr_type_no_symmetry)
896 IF (penalty_occ_local)
THEN
899 matrix_s=qs_matrix_s, &
902 IF (cell%orthorhombic)
THEN
907 ALLOCATE (weights(6))
912 ALLOCATE (op_sm_set_qs(2, dim_op))
913 ALLOCATE (op_sm_set_almo(2, dim_op))
916 DO reim = 1,
SIZE(op_sm_set_qs, 1)
917 NULLIFY (op_sm_set_qs(reim, idim0)%matrix)
918 ALLOCATE (op_sm_set_qs(reim, idim0)%matrix)
919 CALL dbcsr_copy(op_sm_set_qs(reim, idim0)%matrix, qs_matrix_s(1)%matrix, &
921 NULLIFY (op_sm_set_almo(reim, idim0)%matrix)
922 ALLOCATE (op_sm_set_almo(reim, idim0)%matrix)
923 CALL dbcsr_copy(op_sm_set_almo(reim, idim0)%matrix, almo_scf_env%matrix_s(1), &
925 CALL dbcsr_set(op_sm_set_almo(reim, idim0)%matrix, 0.0_dp)
937 m_t_in=m_t_in_local, &
938 m_t0=almo_scf_env%matrix_t_blk, &
939 m_quench_t=quench_t, &
940 m_overlap=almo_scf_env%matrix_s(1), &
941 m_sigma_tmpl=almo_scf_env%matrix_sigma_inv, &
943 xalmo_history=almo_scf_env%xalmo_history, &
944 assume_t0_q0x=assume_t0_q0x, &
945 optimize_theta=optimize_theta, &
946 envelope_amplitude=almo_scf_env%envelope_amplitude, &
947 eps_filter=almo_scf_env%eps_filter, &
948 order_lanczos=almo_scf_env%order_lanczos, &
949 eps_lanczos=almo_scf_env%eps_lanczos, &
950 max_iter_lanczos=almo_scf_env%max_iter_lanczos, &
951 nocc_of_domain=almo_scf_env%nocc_of_domain)
953 ndomains = almo_scf_env%ndomains
954 ALLOCATE (domain_r_down(ndomains, nspins))
956 ALLOCATE (bad_modes_projector_down(ndomains, nspins))
959 ALLOCATE (prec_vv(nspins))
960 ALLOCATE (siginvtftsiginv(nspins))
961 ALLOCATE (stsiginv_0(nspins))
962 ALLOCATE (ftsiginv(nspins))
963 ALLOCATE (st(nspins))
964 ALLOCATE (prev_grad(nspins))
965 ALLOCATE (grad(nspins))
966 ALLOCATE (prev_step(nspins))
967 ALLOCATE (step(nspins))
968 ALLOCATE (prev_minus_prec_grad(nspins))
969 ALLOCATE (m_sig_sqrti_ii(nspins))
970 ALLOCATE (tempnocc(nspins))
971 ALLOCATE (tempnocc_1(nspins))
972 ALLOCATE (tempoccocc(nspins))
977 template=almo_scf_env%matrix_ks(ispin), &
978 matrix_type=dbcsr_type_no_symmetry)
980 template=almo_scf_env%matrix_sigma(ispin), &
981 matrix_type=dbcsr_type_no_symmetry)
983 template=matrix_t_out(ispin), &
984 matrix_type=dbcsr_type_no_symmetry)
986 template=matrix_t_out(ispin), &
987 matrix_type=dbcsr_type_no_symmetry)
989 template=matrix_t_out(ispin), &
990 matrix_type=dbcsr_type_no_symmetry)
992 template=matrix_t_out(ispin), &
993 matrix_type=dbcsr_type_no_symmetry)
995 template=matrix_t_out(ispin), &
996 matrix_type=dbcsr_type_no_symmetry)
998 template=matrix_t_out(ispin), &
999 matrix_type=dbcsr_type_no_symmetry)
1001 template=matrix_t_out(ispin), &
1002 matrix_type=dbcsr_type_no_symmetry)
1004 template=matrix_t_out(ispin), &
1005 matrix_type=dbcsr_type_no_symmetry)
1007 template=almo_scf_env%matrix_sigma_inv(ispin), &
1008 matrix_type=dbcsr_type_no_symmetry)
1010 template=matrix_t_out(ispin), &
1011 matrix_type=dbcsr_type_no_symmetry)
1013 template=matrix_t_out(ispin), &
1014 matrix_type=dbcsr_type_no_symmetry)
1016 template=almo_scf_env%matrix_sigma_inv(ispin), &
1017 matrix_type=dbcsr_type_no_symmetry)
1020 CALL dbcsr_set(prev_step(ispin), 0.0_dp)
1023 nfullrows_total=nocc(ispin))
1030 matrix_s=almo_scf_env%matrix_s(1), &
1031 subm_s_inv=almo_scf_env%domain_s_inv(:, ispin), &
1032 dpattern=quench_t(ispin), &
1033 map=almo_scf_env%domain_map(ispin), &
1034 node_of_domain=almo_scf_env%cpu_of_domain)
1037 matrix_s=almo_scf_env%matrix_s(1), &
1038 subm_s_sqrt=almo_scf_env%domain_s_sqrt(:, ispin), &
1039 subm_s_sqrt_inv=almo_scf_env%domain_s_sqrt_inv(:, ispin), &
1040 dpattern=almo_scf_env%quench_t(ispin), &
1041 map=almo_scf_env%domain_map(ispin), &
1042 node_of_domain=almo_scf_env%cpu_of_domain)
1046 IF (assume_t0_q0x)
THEN
1051 almo_scf_env%matrix_s(1), &
1052 almo_scf_env%matrix_t_blk(ispin), &
1053 0.0_dp, st(ispin), &
1054 filter_eps=almo_scf_env%eps_filter)
1057 almo_scf_env%matrix_sigma_inv_0deloc(ispin), &
1058 0.0_dp, stsiginv_0(ispin), &
1059 filter_eps=almo_scf_env%eps_filter)
1065 matrix_t=almo_scf_env%matrix_t_blk(ispin), &
1066 matrix_sigma_inv=almo_scf_env%matrix_sigma_inv(ispin), &
1067 matrix_s=almo_scf_env%matrix_s(1), &
1068 subm_r_down=domain_r_down(:, ispin), &
1069 dpattern=quench_t(ispin), &
1070 map=almo_scf_env%domain_map(ispin), &
1071 node_of_domain=almo_scf_env%cpu_of_domain, &
1072 filter_eps=almo_scf_env%eps_filter)
1078 IF (penalty_occ_local)
THEN
1082 almo_scf_env%matrix_s(1), &
1083 matrix_t_in(ispin), &
1084 0.0_dp, tempnocc(ispin), &
1085 filter_eps=almo_scf_env%eps_filter)
1088 almo_scf_env%matrix_sigma_inv(ispin), &
1089 0.0_dp, tempnocc_1(ispin), &
1090 filter_eps=almo_scf_env%eps_filter)
1092 DO idim0 = 1,
SIZE(op_sm_set_qs, 2)
1093 DO reim = 1,
SIZE(op_sm_set_qs, 1)
1096 op_sm_set_almo(reim, idim0)%matrix, almo_scf_env%mat_distr_aos)
1099 op_sm_set_almo(reim, idim0)%matrix, &
1100 matrix_t_in(ispin), &
1101 0.0_dp, tempnocc(ispin), &
1102 filter_eps=almo_scf_env%eps_filter)
1105 matrix_t_in(ispin), &
1107 0.0_dp, tempoccocc(ispin), &
1108 filter_eps=almo_scf_env%eps_filter)
1111 tempnocc_1(ispin), &
1112 tempoccocc(ispin), &
1113 0.0_dp, tempnocc(ispin), &
1114 filter_eps=almo_scf_env%eps_filter)
1118 tempnocc_1(ispin), &
1119 0.0_dp, op_sm_set_almo(reim, idim0)%matrix, &
1120 filter_eps=almo_scf_env%eps_filter)
1130 outer_max_iter = optimizer%max_iter_outer_loop
1131 outer_prepare_to_exit = .false.
1134 grad_norm_frob = 0.0_dp
1140 max_iter = optimizer%max_iter
1141 prepare_to_exit = .false.
1142 line_search = .false.
1146 line_search_iteration = 0
1149 energy_diff = 0.0_dp
1150 localization_obj_function = 0.0_dp
1151 line_search_error = 0.0_dp
1157 just_started = (iteration == 0) .AND. (outer_iteration == 0)
1159 CALL main_var_to_xalmos_and_loss_func( &
1160 almo_scf_env=almo_scf_env, &
1162 m_main_var_in=m_theta, &
1163 m_t_out=matrix_t_out, &
1164 m_sig_sqrti_ii_out=m_sig_sqrti_ii, &
1165 energy_out=energy_new, &
1166 penalty_out=penalty_func_new, &
1167 m_ftsiginv_out=ftsiginv, &
1168 m_siginvtftsiginv_out=siginvtftsiginv, &
1170 m_stsiginv0_in=stsiginv_0, &
1171 m_quench_t_in=quench_t, &
1172 domain_r_down_in=domain_r_down, &
1173 assume_t0_q0x=assume_t0_q0x, &
1174 just_started=just_started, &
1175 optimize_theta=optimize_theta, &
1176 normalize_orbitals=normalize_orbitals, &
1177 perturbation_only=perturbation_only, &
1178 do_penalty=penalty_occ_vol, &
1179 special_case=my_special_case)
1180 IF (penalty_occ_vol)
THEN
1182 energy_new = energy_new + penalty_func_new
1184 DO ispin = 1, nspins
1185 IF (penalty_occ_vol)
THEN
1186 penalty_occ_vol_g_prefactor(ispin) = &
1187 -2.0_dp*penalty_amplitude*spin_factor*nocc(ispin)
1188 penalty_occ_vol_h_prefactor(ispin) = 0.0_dp
1192 localization_obj_function = 0.0_dp
1194 IF (penalty_occ_local)
THEN
1195 DO ispin = 1, nspins
1198 localization_obj_function = 0.0_dp
1199 CALL dbcsr_get_info(almo_scf_env%matrix_sigma_inv(ispin), nfullrows_total=nmo)
1201 ALLOCATE (reim_diag(nmo))
1205 DO idim0 = 1,
SIZE(op_sm_set_qs, 2)
1209 DO reim = 1,
SIZE(op_sm_set_qs, 1)
1212 op_sm_set_almo(reim, idim0)%matrix, &
1213 matrix_t_out(ispin), &
1214 0.0_dp, tempnocc(ispin), &
1215 filter_eps=almo_scf_env%eps_filter)
1218 matrix_t_out(ispin), &
1220 0.0_dp, tempoccocc(ispin), &
1221 filter_eps=almo_scf_env%eps_filter)
1225 CALL group%sum(reim_diag)
1226 z2(:) = z2(:) + reim_diag(:)*reim_diag(:)
1233 fval = -weights(idim0)*log(abs(z2(ielem)))
1235 fval = weights(idim0) - weights(idim0)*abs(z2(ielem))
1237 fval = weights(idim0) - weights(idim0)*sqrt(abs(z2(ielem)))
1239 localization_obj_function = localization_obj_function + fval
1245 DEALLOCATE (reim_diag)
1247 energy_new = energy_new + localiz_coeff*localization_obj_function
1252 DO ispin = 1, nspins
1254 IF (just_started .AND. almo_mathematica)
THEN
1255 cpwarn_if(ispin > 1,
"Mathematica files will be overwritten")
1256 CALL print_mathematica_matrix(almo_scf_env%matrix_s(1),
"matrixS.dat")
1257 CALL print_mathematica_matrix(almo_scf_env%matrix_ks(ispin),
"matrixF.dat")
1258 CALL print_mathematica_matrix(matrix_t_out(ispin),
"matrixT.dat")
1259 CALL print_mathematica_matrix(quench_t(ispin),
"matrixQ.dat")
1265 IF (line_search_iteration == 0 .AND. iteration /= 0)
THEN
1266 CALL dbcsr_copy(prev_grad(ispin), grad(ispin))
1272 skip_grad = (iteration > 0 .AND. &
1273 fixed_line_search_niter /= 0 .AND. &
1274 line_search_iteration /= fixed_line_search_niter)
1276 IF (.NOT. skip_grad)
THEN
1278 DO ispin = 1, nspins
1280 CALL compute_gradient( &
1281 m_grad_out=grad(ispin), &
1282 m_ks=almo_scf_env%matrix_ks(ispin), &
1283 m_s=almo_scf_env%matrix_s(1), &
1284 m_t=matrix_t_out(ispin), &
1285 m_t0=almo_scf_env%matrix_t_blk(ispin), &
1286 m_siginv=almo_scf_env%matrix_sigma_inv(ispin), &
1287 m_quench_t=quench_t(ispin), &
1288 m_ftsiginv=ftsiginv(ispin), &
1289 m_siginvtftsiginv=siginvtftsiginv(ispin), &
1291 m_stsiginv0=stsiginv_0(ispin), &
1292 m_theta=m_theta(ispin), &
1293 m_sig_sqrti_ii=m_sig_sqrti_ii(ispin), &
1294 domain_s_inv=almo_scf_env%domain_s_inv(:, ispin), &
1295 domain_r_down=domain_r_down(:, ispin), &
1296 cpu_of_domain=almo_scf_env%cpu_of_domain, &
1297 domain_map=almo_scf_env%domain_map(ispin), &
1298 assume_t0_q0x=assume_t0_q0x, &
1299 optimize_theta=optimize_theta, &
1300 normalize_orbitals=normalize_orbitals, &
1301 penalty_occ_vol=penalty_occ_vol, &
1302 penalty_occ_vol_prefactor=penalty_occ_vol_g_prefactor(ispin), &
1303 envelope_amplitude=almo_scf_env%envelope_amplitude, &
1304 eps_filter=almo_scf_env%eps_filter, &
1305 spin_factor=spin_factor, &
1306 special_case=my_special_case, &
1307 penalty_occ_local=penalty_occ_local, &
1308 op_sm_set=op_sm_set_almo, &
1310 energy_coeff=energy_coeff, &
1311 localiz_coeff=localiz_coeff)
1320 IF (blissful_neglect)
THEN
1321 DO ispin = 1, nspins
1324 IF (iteration == 0)
THEN
1325 CALL compute_preconditioner( &
1326 domain_prec_out=almo_scf_env%domain_preconditioner(:, ispin), &
1327 bad_modes_projector_down_out=bad_modes_projector_down(:, ispin), &
1328 m_prec_out=prec_vv(ispin), &
1329 m_ks=almo_scf_env%matrix_ks(ispin), &
1330 m_s=almo_scf_env%matrix_s(1), &
1331 m_siginv=almo_scf_env%matrix_sigma_inv(ispin), &
1332 m_quench_t=quench_t(ispin), &
1333 m_ftsiginv=ftsiginv(ispin), &
1334 m_siginvtftsiginv=siginvtftsiginv(ispin), &
1336 para_env=almo_scf_env%para_env, &
1337 blacs_env=almo_scf_env%blacs_env, &
1338 nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin), &
1339 domain_s_inv=almo_scf_env%domain_s_inv(:, ispin), &
1340 domain_s_inv_half=almo_scf_env%domain_s_sqrt_inv(:, ispin), &
1341 domain_s_half=almo_scf_env%domain_s_sqrt(:, ispin), &
1342 domain_r_down=domain_r_down(:, ispin), &
1343 cpu_of_domain=almo_scf_env%cpu_of_domain, &
1344 domain_map=almo_scf_env%domain_map(ispin), &
1345 assume_t0_q0x=assume_t0_q0x, &
1346 penalty_occ_vol=penalty_occ_vol, &
1347 penalty_occ_vol_prefactor=penalty_occ_vol_g_prefactor(ispin), &
1348 eps_filter=almo_scf_env%eps_filter, &
1349 neg_thr=optimizer%neglect_threshold, &
1350 spin_factor=spin_factor, &
1351 skip_inversion=.false., &
1352 special_case=my_special_case)
1356 matrix_in=grad(ispin), &
1357 matrix_out=grad(ispin), &
1358 operator1=almo_scf_env%domain_s_inv(:, ispin), &
1359 operator2=bad_modes_projector_down(:, ispin), &
1360 dpattern=quench_t(ispin), &
1361 map=almo_scf_env%domain_map(ispin), &
1362 node_of_domain=almo_scf_env%cpu_of_domain, &
1364 filter_eps=almo_scf_env%eps_filter)
1371 DO ispin = 1, nspins
1374 grad_norm = maxval(grad_norm_spin)
1376 converged = (grad_norm <= optimizer%eps_error)
1377 IF (converged .OR. (iteration >= max_iter))
THEN
1378 prepare_to_exit = .true.
1381 IF (optimizer%early_stopping_on .AND. just_started)
THEN
1382 prepare_to_exit = .false.
1385 IF (grad_norm < almo_scf_env%eps_prev_guess)
THEN
1390 IF (.NOT. prepare_to_exit)
THEN
1395 IF (iteration /= 0)
THEN
1397 IF (fixed_line_search_niter == 0)
THEN
1401 IF (.NOT. line_search)
THEN
1403 line_search = .true.
1404 line_search_iteration = line_search_iteration + 1
1410 line_search_error = 0.0_dp
1414 DO ispin = 1, nspins
1416 CALL dbcsr_dot(grad(ispin), step(ispin), tempreal)
1417 line_search_error = line_search_error + tempreal
1418 CALL dbcsr_dot(grad(ispin), grad(ispin), tempreal)
1419 denom = denom + tempreal
1420 CALL dbcsr_dot(step(ispin), step(ispin), tempreal)
1421 denom2 = denom2 + tempreal
1427 line_search_error = line_search_error/sqrt(denom)/sqrt(denom2)
1429 IF (abs(line_search_error) > optimizer%lin_search_eps_error)
THEN
1430 line_search = .true.
1431 line_search_iteration = line_search_iteration + 1
1433 line_search = .false.
1434 line_search_iteration = 0
1435 IF (grad_norm < eps_skip_gradients)
THEN
1436 fixed_line_search_niter = abs(almo_scf_env%integer04)
1444 IF (.NOT. line_search)
THEN
1445 line_search = .true.
1446 line_search_iteration = line_search_iteration + 1
1448 IF (line_search_iteration == fixed_line_search_niter)
THEN
1449 line_search = .false.
1450 line_search_iteration = 0
1451 line_search_iteration = line_search_iteration + 1
1459 IF (line_search)
THEN
1460 energy_diff = 0.0_dp
1462 energy_diff = energy_new - energy_old
1463 energy_old = energy_new
1467 IF (.NOT. line_search)
THEN
1469 cg_iteration = cg_iteration + 1
1472 DO ispin = 1, nspins
1473 CALL dbcsr_copy(prev_step(ispin), step(ispin))
1477 SELECT CASE (prec_type)
1481 CALL newton_grad_to_step( &
1482 optimizer=almo_scf_env%opt_xalmo_newton_pcg_solver, &
1485 m_s=almo_scf_env%matrix_s(:), &
1486 m_ks=almo_scf_env%matrix_ks(:), &
1487 m_siginv=almo_scf_env%matrix_sigma_inv(:), &
1488 m_quench_t=quench_t(:), &
1489 m_ftsiginv=ftsiginv(:), &
1490 m_siginvtftsiginv=siginvtftsiginv(:), &
1492 m_t=matrix_t_out(:), &
1493 m_sig_sqrti_ii=m_sig_sqrti_ii(:), &
1494 domain_s_inv=almo_scf_env%domain_s_inv(:, :), &
1495 domain_r_down=domain_r_down(:, :), &
1496 domain_map=almo_scf_env%domain_map(:), &
1497 cpu_of_domain=almo_scf_env%cpu_of_domain, &
1498 nocc_of_domain=almo_scf_env%nocc_of_domain(:, :), &
1499 para_env=almo_scf_env%para_env, &
1500 blacs_env=almo_scf_env%blacs_env, &
1501 eps_filter=almo_scf_env%eps_filter, &
1502 optimize_theta=optimize_theta, &
1503 penalty_occ_vol=penalty_occ_vol, &
1504 normalize_orbitals=normalize_orbitals, &
1505 penalty_occ_vol_prefactor=penalty_occ_vol_g_prefactor(:), &
1506 penalty_occ_vol_pf2=penalty_occ_vol_h_prefactor(:), &
1507 special_case=my_special_case &
1513 IF (.NOT. blissful_neglect .AND. &
1514 ((just_started .AND. perturbation_only) .OR. &
1515 (iteration == 0 .AND. (.NOT. perturbation_only))) &
1519 DO ispin = 1, nspins
1520 CALL compute_preconditioner( &
1521 domain_prec_out=almo_scf_env%domain_preconditioner(:, ispin), &
1522 m_prec_out=prec_vv(ispin), &
1523 m_ks=almo_scf_env%matrix_ks(ispin), &
1524 m_s=almo_scf_env%matrix_s(1), &
1525 m_siginv=almo_scf_env%matrix_sigma_inv(ispin), &
1526 m_quench_t=quench_t(ispin), &
1527 m_ftsiginv=ftsiginv(ispin), &
1528 m_siginvtftsiginv=siginvtftsiginv(ispin), &
1530 para_env=almo_scf_env%para_env, &
1531 blacs_env=almo_scf_env%blacs_env, &
1532 nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin), &
1533 domain_s_inv=almo_scf_env%domain_s_inv(:, ispin), &
1534 domain_r_down=domain_r_down(:, ispin), &
1535 cpu_of_domain=almo_scf_env%cpu_of_domain, &
1536 domain_map=almo_scf_env%domain_map(ispin), &
1537 assume_t0_q0x=assume_t0_q0x, &
1538 penalty_occ_vol=penalty_occ_vol, &
1539 penalty_occ_vol_prefactor=penalty_occ_vol_g_prefactor(ispin), &
1540 eps_filter=almo_scf_env%eps_filter, &
1542 spin_factor=spin_factor, &
1543 skip_inversion=.false., &
1544 special_case=my_special_case)
1551 DO ispin = 1, nspins
1556 0.0_dp, step(ispin), &
1557 filter_eps=almo_scf_env%eps_filter)
1564 IF (optimize_theta)
THEN
1565 cpabort(
"theta is NYI")
1568 DO ispin = 1, nspins
1571 matrix_in=grad(ispin), &
1572 matrix_out=step(ispin), &
1573 operator1=almo_scf_env%domain_preconditioner(:, ispin), &
1574 dpattern=quench_t(ispin), &
1575 map=almo_scf_env%domain_map(ispin), &
1576 node_of_domain=almo_scf_env%cpu_of_domain, &
1578 filter_eps=almo_scf_env%eps_filter)
1588 DO ispin = 1, nspins
1598 IF (iteration == 0)
THEN
1599 reset_conjugator = .true.
1603 IF (.NOT. reset_conjugator)
THEN
1605 CALL compute_cg_beta( &
1607 reset_conjugator=reset_conjugator, &
1608 conjugator=optimizer%conjugator, &
1610 prev_grad=prev_grad(:), &
1612 prev_step=prev_step(:), &
1613 prev_minus_prec_grad=prev_minus_prec_grad(:) &
1618 IF (reset_conjugator)
THEN
1621 IF (unit_nr > 0 .AND. (.NOT. just_started))
THEN
1622 WRITE (unit_nr,
'(T2,A35)')
"Re-setting conjugator to zero"
1624 reset_conjugator = .false.
1629 DO ispin = 1, nspins
1631 CALL dbcsr_copy(prev_minus_prec_grad(ispin), step(ispin))
1634 CALL dbcsr_add(step(ispin), prev_step(ispin), 1.0_dp, beta)
1641 IF (.NOT. line_search)
THEN
1647 DO ispin = 1, nspins
1648 CALL dbcsr_dot(grad(ispin), step(ispin), tempreal)
1651 IF (iteration == 0)
THEN
1652 step_size = optimizer%lin_search_step_size_guess
1654 IF (next_step_size_guess <= 0.0_dp)
THEN
1655 step_size = optimizer%lin_search_step_size_guess
1658 step_size = next_step_size_guess*1.05_dp
1661 next_step_size_guess = step_size
1663 IF (fixed_line_search_niter == 0)
THEN
1666 DO ispin = 1, nspins
1667 CALL dbcsr_dot(grad(ispin), step(ispin), tempreal)
1672 appr_sec_der = (g1 - g0)/step_size
1673 step_size = -g1/appr_sec_der
1681 appr_sec_der = 2.0_dp*((e1 - e0)/step_size - g0)/step_size
1682 g1 = appr_sec_der*step_size + g0
1683 step_size = -g1/appr_sec_der
1687 next_step_size_guess = next_step_size_guess + step_size
1691 DO ispin = 1, nspins
1692 CALL dbcsr_add(m_theta(ispin), step(ispin), 1.0_dp, step_size)
1697 IF (line_search)
THEN
1704 IF (unit_nr > 0)
THEN
1705 iter_type = trim(
"ALMO SCF "//iter_type)
1706 WRITE (unit_nr,
'(T2,A13,I6,F23.10,E14.5,F14.9,F9.2)') &
1707 iter_type, iteration, &
1708 energy_new, energy_diff, grad_norm, &
1710 IF (penalty_occ_local .OR. penalty_occ_vol)
THEN
1711 WRITE (unit_nr,
'(T2,A25,F23.10)') &
1712 "Energy component:", (energy_new - penalty_func_new - localization_obj_function)
1714 IF (penalty_occ_local)
THEN
1715 WRITE (unit_nr,
'(T2,A25,F23.10)') &
1716 "Localization component:", localization_obj_function
1718 IF (penalty_occ_vol)
THEN
1719 WRITE (unit_nr,
'(T2,A25,F23.10)') &
1720 "Penalty component:", penalty_func_new
1725 IF (penalty_occ_vol)
THEN
1726 almo_scf_env%almo_scf_energy = energy_new - penalty_func_new - localization_obj_function
1728 almo_scf_env%almo_scf_energy = energy_new - localization_obj_function
1734 iteration = iteration + 1
1735 IF (prepare_to_exit)
EXIT
1739 IF (converged .OR. (outer_iteration >= outer_max_iter))
THEN
1740 outer_prepare_to_exit = .true.
1743 outer_iteration = outer_iteration + 1
1744 IF (outer_prepare_to_exit)
EXIT
1748 DO ispin = 1, nspins
1749 IF (converged .AND. almo_mathematica)
THEN
1750 cpwarn_if(ispin > 1,
"Mathematica files will be overwritten")
1751 CALL print_mathematica_matrix(matrix_t_out(ispin),
"matrixTf.dat")
1758 CALL wrap_up_xalmo_scf( &
1760 almo_scf_env=almo_scf_env, &
1761 perturbation_in=perturbation_only, &
1762 m_xalmo_in=matrix_t_out, &
1763 m_quench_in=quench_t, &
1764 energy_inout=energy_new)
1768 DO ispin = 1, nspins
1789 DEALLOCATE (tempnocc)
1790 DEALLOCATE (tempnocc_1)
1791 DEALLOCATE (tempoccocc)
1792 DEALLOCATE (prec_vv)
1793 DEALLOCATE (siginvtftsiginv)
1794 DEALLOCATE (stsiginv_0)
1795 DEALLOCATE (ftsiginv)
1797 DEALLOCATE (prev_grad)
1799 DEALLOCATE (prev_step)
1801 DEALLOCATE (prev_minus_prec_grad)
1802 DEALLOCATE (m_sig_sqrti_ii)
1804 DEALLOCATE (domain_r_down)
1805 DEALLOCATE (bad_modes_projector_down)
1807 DEALLOCATE (penalty_occ_vol_g_prefactor)
1808 DEALLOCATE (penalty_occ_vol_h_prefactor)
1809 DEALLOCATE (grad_norm_spin)
1812 DEALLOCATE (m_theta, m_t_in_local)
1813 IF (penalty_occ_local)
THEN
1814 DO idim0 = 1, dim_op
1815 DO reim = 1,
SIZE(op_sm_set_qs, 1)
1816 DEALLOCATE (op_sm_set_qs(reim, idim0)%matrix)
1817 DEALLOCATE (op_sm_set_almo(reim, idim0)%matrix)
1820 DEALLOCATE (op_sm_set_qs)
1821 DEALLOCATE (op_sm_set_almo)
1822 DEALLOCATE (weights)
1825 IF (.NOT. converged .AND. .NOT. optimizer%early_stopping_on)
THEN
1826 cpabort(
"Optimization not converged! ")
1829 CALL timestop(handle)
1850 matrix_s, matrix_mo_in, matrix_mo_out, &
1851 template_matrix_sigma, overlap_determinant, &
1852 mat_distr_aos, virtuals, eps_filter)
1856 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:), &
1857 INTENT(INOUT) :: matrix_mo_in, matrix_mo_out
1858 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:), &
1859 INTENT(IN) :: template_matrix_sigma
1860 REAL(kind=
dp),
INTENT(INOUT) :: overlap_determinant
1861 INTEGER,
INTENT(IN) :: mat_distr_aos
1862 LOGICAL,
INTENT(IN) :: virtuals
1863 REAL(kind=
dp),
INTENT(IN) :: eps_filter
1865 CHARACTER(len=*),
PARAMETER :: routinen =
'almo_scf_construct_nlmos'
1867 CHARACTER(LEN=30) :: iter_type, print_string
1868 INTEGER :: cg_iteration, dim_op, handle, iatom, idim0, isgf, ispin, iteration, &
1869 line_search_iteration, linear_search_type, max_iter, natom, ncol, nspins, &
1870 outer_iteration, outer_max_iter, prec_type, reim, unit_nr
1871 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: first_sgf, last_sgf, nocc, nsgf
1872 LOGICAL :: converged, d_bfgs, just_started, l_bfgs, &
1873 line_search, outer_prepare_to_exit, &
1874 prepare_to_exit, reset_conjugator
1875 REAL(kind=
dp) :: appr_sec_der, beta, bfgs_rho, bfgs_sum, denom, denom2, e0, e1, g0, g0sign, &
1876 g1, g1sign, grad_norm, line_search_error, localization_obj_function, &
1877 localization_obj_function_ispin, next_step_size_guess, obj_function_ispin, objf_diff, &
1878 objf_new, objf_old, penalty_amplitude, penalty_func_ispin, penalty_func_new, spin_factor, &
1879 step_size, t1, t2, tempreal
1880 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: diagonal, grad_norm_spin, &
1881 penalty_vol_prefactor, &
1882 suggested_vol_penalty, weights
1885 TYPE(
dbcsr_p_type),
DIMENSION(:),
POINTER :: qs_matrix_s
1886 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
POINTER :: op_sm_set_almo, op_sm_set_qs
1887 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:) :: approx_inv_hessian, bfgs_s, bfgs_y, grad, &
1888 m_s0, m_sig_sqrti_ii, m_siginv, m_sigma, m_t_mo_local, m_theta, m_theta_normalized, &
1889 prev_grad, prev_m_theta, prev_minus_prec_grad, prev_step, step, tempnocc1, tempoccocc1, &
1890 tempoccocc2, tempoccocc3
1891 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:, :, :) :: m_b0
1895 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
1897 CALL timeset(routinen, handle)
1901 IF (logger%para_env%is_source())
THEN
1907 nspins =
SIZE(matrix_mo_in)
1909 IF (unit_nr > 0)
THEN
1911 IF (.NOT. virtuals)
THEN
1912 WRITE (unit_nr,
'(T2,A,A,A)') repeat(
"-", 24), &
1913 " Optimization of occupied NLMOs ", repeat(
"-", 23)
1915 WRITE (unit_nr,
'(T2,A,A,A)') repeat(
"-", 24), &
1916 " Optimization of virtual NLMOs ", repeat(
"-", 24)
1919 WRITE (unit_nr,
'(T2,A13,A6,A23,A14,A14,A9)')
"Method",
"Iter", &
1920 "Objective Function",
"Change",
"Convergence",
"Time"
1921 WRITE (unit_nr,
'(T2,A)') repeat(
"-", 79)
1924 NULLIFY (particle_set)
1927 matrix_s=qs_matrix_s, &
1929 particle_set=particle_set, &
1930 qs_kind_set=qs_kind_set)
1932 natom =
SIZE(particle_set, 1)
1933 ALLOCATE (first_sgf(natom))
1934 ALLOCATE (last_sgf(natom))
1935 ALLOCATE (nsgf(natom))
1938 first_sgf=first_sgf, last_sgf=last_sgf, nsgf=nsgf)
1942 ALLOCATE (m_theta(nspins))
1943 DO ispin = 1, nspins
1945 template=template_matrix_sigma(ispin), &
1946 matrix_type=dbcsr_type_no_symmetry)
1952 SELECT CASE (optimizer%opt_penalty%operator_type)
1955 IF (cell%orthorhombic)
THEN
1960 ALLOCATE (weights(6))
1963 ALLOCATE (op_sm_set_qs(2, dim_op))
1964 ALLOCATE (op_sm_set_almo(2, dim_op))
1966 ALLOCATE (m_b0(2, dim_op, nspins))
1967 DO idim0 = 1, dim_op
1968 DO reim = 1,
SIZE(op_sm_set_qs, 1)
1969 NULLIFY (op_sm_set_qs(reim, idim0)%matrix, op_sm_set_almo(reim, idim0)%matrix)
1970 ALLOCATE (op_sm_set_qs(reim, idim0)%matrix)
1971 ALLOCATE (op_sm_set_almo(reim, idim0)%matrix)
1972 CALL dbcsr_copy(op_sm_set_qs(reim, idim0)%matrix, qs_matrix_s(1)%matrix, &
1974 CALL dbcsr_copy(op_sm_set_almo(reim, idim0)%matrix, matrix_s, &
1976 CALL dbcsr_set(op_sm_set_almo(reim, idim0)%matrix, 0.0_dp)
1977 DO ispin = 1, nspins
1979 template=m_theta(ispin), &
1980 matrix_type=dbcsr_type_no_symmetry)
1981 CALL dbcsr_set(m_b0(reim, idim0, ispin), 0.0_dp)
1991 ALLOCATE (weights(dim_op))
1994 ALLOCATE (m_b0(1, dim_op, nspins))
1996 DO idim0 = 1, dim_op
1998 DO ispin = 1, nspins
2000 template=m_theta(ispin), &
2001 matrix_type=dbcsr_type_no_symmetry)
2002 CALL dbcsr_set(m_b0(reim, idim0, ispin), 0.0_dp)
2010 penalty_amplitude = optimizer%opt_penalty%penalty_strength
2013 prec_type = optimizer%preconditioner
2019 IF (l_bfgs .AND. (optimizer%conjugator /=
cg_zero))
THEN
2020 cpabort(
"Cannot use conjugators with BFGS")
2023 CALL lbfgs_create(nlmo_lbfgs_history, nspins, nstore=10)
2026 IF (nspins == 1)
THEN
2027 spin_factor = 2.0_dp
2029 spin_factor = 1.0_dp
2032 ALLOCATE (grad_norm_spin(nspins))
2033 ALLOCATE (nocc(nspins))
2034 ALLOCATE (penalty_vol_prefactor(nspins))
2035 ALLOCATE (suggested_vol_penalty(nspins))
2041 ALLOCATE (m_t_mo_local(nspins))
2042 DO ispin = 1, nspins
2044 template=matrix_mo_in(ispin), &
2045 matrix_type=dbcsr_type_no_symmetry)
2046 CALL dbcsr_copy(m_t_mo_local(ispin), matrix_mo_in(ispin))
2049 ALLOCATE (approx_inv_hessian(nspins))
2050 ALLOCATE (m_theta_normalized(nspins))
2051 ALLOCATE (prev_m_theta(nspins))
2052 ALLOCATE (m_s0(nspins))
2053 ALLOCATE (prev_grad(nspins))
2054 ALLOCATE (grad(nspins))
2055 ALLOCATE (prev_step(nspins))
2056 ALLOCATE (step(nspins))
2057 ALLOCATE (prev_minus_prec_grad(nspins))
2058 ALLOCATE (m_sig_sqrti_ii(nspins))
2059 ALLOCATE (m_sigma(nspins))
2060 ALLOCATE (m_siginv(nspins))
2061 ALLOCATE (tempnocc1(nspins))
2062 ALLOCATE (tempoccocc1(nspins))
2063 ALLOCATE (tempoccocc2(nspins))
2064 ALLOCATE (tempoccocc3(nspins))
2065 ALLOCATE (bfgs_y(nspins))
2066 ALLOCATE (bfgs_s(nspins))
2068 DO ispin = 1, nspins
2072 template=matrix_mo_out(ispin), &
2073 matrix_type=dbcsr_type_no_symmetry)
2075 template=m_theta(ispin), &
2076 matrix_type=dbcsr_type_no_symmetry)
2078 template=m_theta(ispin), &
2079 matrix_type=dbcsr_type_no_symmetry)
2081 template=m_theta(ispin), &
2082 matrix_type=dbcsr_type_no_symmetry)
2084 template=m_theta(ispin), &
2085 matrix_type=dbcsr_type_no_symmetry)
2087 template=m_theta(ispin), &
2088 matrix_type=dbcsr_type_no_symmetry)
2090 template=m_theta(ispin), &
2091 matrix_type=dbcsr_type_no_symmetry)
2093 template=m_theta(ispin), &
2094 matrix_type=dbcsr_type_no_symmetry)
2096 template=m_theta(ispin), &
2097 matrix_type=dbcsr_type_no_symmetry)
2099 template=m_theta(ispin), &
2100 matrix_type=dbcsr_type_no_symmetry)
2102 template=m_theta(ispin), &
2103 matrix_type=dbcsr_type_no_symmetry)
2105 template=m_theta(ispin), &
2106 matrix_type=dbcsr_type_no_symmetry)
2108 template=m_theta(ispin), &
2109 matrix_type=dbcsr_type_no_symmetry)
2111 template=m_theta(ispin), &
2112 matrix_type=dbcsr_type_no_symmetry)
2114 template=m_theta(ispin), &
2115 matrix_type=dbcsr_type_no_symmetry)
2117 template=m_theta(ispin), &
2118 matrix_type=dbcsr_type_no_symmetry)
2120 template=m_theta(ispin), &
2121 matrix_type=dbcsr_type_no_symmetry)
2123 template=m_theta(ispin), &
2124 matrix_type=dbcsr_type_no_symmetry)
2127 CALL dbcsr_set(prev_step(ispin), 0.0_dp)
2130 nfullrows_total=nocc(ispin))
2132 penalty_vol_prefactor(ispin) = -penalty_amplitude
2137 m_t_mo_local(ispin), &
2138 0.0_dp, tempnocc1(ispin), &
2139 filter_eps=eps_filter)
2141 m_t_mo_local(ispin), &
2143 0.0_dp, m_s0(ispin), &
2144 filter_eps=eps_filter)
2146 SELECT CASE (optimizer%opt_penalty%operator_type)
2151 DO idim0 = 1,
SIZE(op_sm_set_qs, 2)
2153 DO reim = 1,
SIZE(op_sm_set_qs, 1)
2156 op_sm_set_almo(reim, idim0)%matrix, mat_distr_aos)
2159 op_sm_set_almo(reim, idim0)%matrix, &
2160 m_t_mo_local(ispin), &
2161 0.0_dp, tempnocc1(ispin), &
2162 filter_eps=eps_filter)
2165 m_t_mo_local(ispin), &
2167 0.0_dp, m_b0(reim, idim0, ispin), &
2168 filter_eps=eps_filter)
2170 DEALLOCATE (op_sm_set_qs(reim, idim0)%matrix)
2171 DEALLOCATE (op_sm_set_almo(reim, idim0)%matrix)
2182 isgf = first_sgf(iatom)
2187 m_t_mo_local(ispin), &
2188 0.0_dp, tempnocc1(ispin), &
2189 filter_eps=eps_filter)
2192 m_t_mo_local(ispin), &
2194 0.0_dp, m_b0(1, iatom, ispin), &
2195 first_k=isgf, last_k=isgf + ncol - 1, &
2196 filter_eps=eps_filter)
2200 m_t_mo_local(ispin), &
2201 0.0_dp, tempnocc1(ispin), &
2202 first_k=isgf, last_k=isgf + ncol - 1, &
2203 filter_eps=eps_filter)
2206 m_t_mo_local(ispin), &
2208 1.0_dp, m_b0(1, iatom, ispin), &
2209 filter_eps=eps_filter)
2217 IF (optimizer%opt_penalty%operator_type ==
op_loc_berry)
THEN
2218 DO idim0 = 1,
SIZE(op_sm_set_qs, 2)
2219 DO reim = 1,
SIZE(op_sm_set_qs, 1)
2220 DEALLOCATE (op_sm_set_qs(reim, idim0)%matrix)
2221 DEALLOCATE (op_sm_set_almo(reim, idim0)%matrix)
2224 DEALLOCATE (op_sm_set_qs, op_sm_set_almo)
2228 outer_max_iter = optimizer%max_iter_outer_loop
2229 outer_prepare_to_exit = .false.
2232 penalty_func_new = 0.0_dp
2233 linear_search_type = 1
2234 localization_obj_function = 0.0_dp
2235 penalty_func_new = 0.0_dp
2240 max_iter = optimizer%max_iter
2241 prepare_to_exit = .false.
2242 line_search = .false.
2246 line_search_iteration = 0
2247 obj_function_ispin = 0.0_dp
2251 line_search_error = 0.0_dp
2253 next_step_size_guess = 0.0_dp
2257 just_started = (iteration == 0) .AND. (outer_iteration == 0)
2259 DO ispin = 1, nspins
2265 m_s0(ispin), m_theta(ispin), 0.0_dp, &
2266 tempoccocc1(ispin), &
2267 filter_eps=eps_filter)
2268 CALL dbcsr_set(m_sig_sqrti_ii(ispin), 0.0_dp)
2271 m_theta(ispin), tempoccocc1(ispin), 0.0_dp, &
2272 m_sig_sqrti_ii(ispin), &
2273 retain_sparsity=.true.)
2274 ALLOCATE (diagonal(nocc(ispin)))
2276 CALL group%sum(diagonal)
2278 diagonal(:) = 1.0_dp/sqrt(diagonal(:))
2279 CALL dbcsr_set(m_sig_sqrti_ii(ispin), 0.0_dp)
2281 DEALLOCATE (diagonal)
2285 m_sig_sqrti_ii(ispin), &
2286 0.0_dp, m_theta_normalized(ispin), &
2287 filter_eps=eps_filter)
2291 m_t_mo_local(ispin), &
2292 m_theta_normalized(ispin), &
2293 0.0_dp, matrix_mo_out(ispin), &
2294 filter_eps=eps_filter)
2299 localization_obj_function = 0.0_dp
2300 penalty_func_new = 0.0_dp
2301 DO ispin = 1, nspins
2303 CALL compute_obj_nlmos( &
2304 localization_obj_function_ispin=localization_obj_function_ispin, &
2305 penalty_func_ispin=penalty_func_ispin, &
2306 overlap_determinant=overlap_determinant, &
2307 m_sigma=m_sigma(ispin), &
2309 m_b0=m_b0(:, :, ispin), &
2310 m_theta_normalized=m_theta_normalized(ispin), &
2311 template_matrix_mo=matrix_mo_out(ispin), &
2314 just_started=just_started, &
2315 penalty_vol_prefactor=penalty_vol_prefactor(ispin), &
2316 penalty_amplitude=penalty_amplitude, &
2317 eps_filter=eps_filter)
2319 localization_obj_function = localization_obj_function + localization_obj_function_ispin
2320 penalty_func_new = penalty_func_new + penalty_func_ispin
2323 objf_new = penalty_func_new + localization_obj_function
2325 DO ispin = 1, nspins
2329 IF (line_search_iteration == 0 .AND. iteration /= 0)
THEN
2330 CALL dbcsr_copy(prev_grad(ispin), grad(ispin))
2336 DO ispin = 1, nspins
2339 matrix_inverse=m_siginv(ispin), &
2340 matrix=m_sigma(ispin), &
2341 threshold=eps_filter*10.0_dp, &
2342 filter_eps=eps_filter, &
2345 CALL compute_gradient_nlmos( &
2346 m_grad_out=grad(ispin), &
2347 m_b0=m_b0(:, :, ispin), &
2350 m_theta_normalized=m_theta_normalized(ispin), &
2351 m_siginv=m_siginv(ispin), &
2352 m_sig_sqrti_ii=m_sig_sqrti_ii(ispin), &
2353 penalty_vol_prefactor=penalty_vol_prefactor(ispin), &
2354 eps_filter=eps_filter, &
2355 suggested_vol_penalty=suggested_vol_penalty(ispin))
2360 DO ispin = 1, nspins
2363 grad_norm = maxval(grad_norm_spin)
2365 converged = (grad_norm <= optimizer%eps_error)
2366 IF (converged .OR. (iteration >= max_iter))
THEN
2367 prepare_to_exit = .true.
2371 IF (.NOT. prepare_to_exit)
THEN
2376 IF (iteration /= 0)
THEN
2380 IF (.NOT. line_search)
THEN
2382 line_search = .true.
2383 line_search_iteration = line_search_iteration + 1
2389 line_search_error = 0.0_dp
2393 DO ispin = 1, nspins
2395 CALL dbcsr_dot(grad(ispin), step(ispin), tempreal)
2396 line_search_error = line_search_error + tempreal
2397 CALL dbcsr_dot(grad(ispin), grad(ispin), tempreal)
2398 denom = denom + tempreal
2399 CALL dbcsr_dot(step(ispin), step(ispin), tempreal)
2400 denom2 = denom2 + tempreal
2406 line_search_error = line_search_error/sqrt(denom)/sqrt(denom2)
2408 IF (abs(line_search_error) > optimizer%lin_search_eps_error)
THEN
2409 line_search = .true.
2410 line_search_iteration = line_search_iteration + 1
2412 line_search = .false.
2413 line_search_iteration = 0
2420 IF (line_search)
THEN
2423 objf_diff = objf_new - objf_old
2428 IF (.NOT. line_search)
THEN
2430 cg_iteration = cg_iteration + 1
2433 DO ispin = 1, nspins
2434 CALL dbcsr_copy(prev_step(ispin), step(ispin))
2442 DO ispin = 1, nspins
2452 IF (iteration == 0)
THEN
2458 IF (nspins > 1)
THEN
2459 DO ispin = 2, nspins
2460 CALL dbcsr_copy(approx_inv_hessian(ispin), approx_inv_hessian(1))
2464 ELSE IF (l_bfgs)
THEN
2466 CALL lbfgs_seed(nlmo_lbfgs_history, m_theta, grad)
2467 DO ispin = 1, nspins
2475 DO ispin = 1, nspins
2489 DO ispin = 1, nspins
2493 CALL dbcsr_add(bfgs_y(ispin), prev_grad(ispin), 1.0_dp, -1.0_dp)
2494 CALL dbcsr_copy(bfgs_s(ispin), m_theta(ispin))
2495 CALL dbcsr_add(bfgs_s(ispin), prev_m_theta(ispin), 1.0_dp, -1.0_dp)
2498 CALL dbcsr_dot(grad(ispin), step(ispin), bfgs_rho)
2499 bfgs_rho = 1.0_dp/bfgs_rho
2502 CALL dbcsr_dot(bfgs_y(ispin), bfgs_y(ispin), bfgs_sum)
2505 CALL dbcsr_copy(tempoccocc2(ispin), approx_inv_hessian(ispin))
2509 CALL dbcsr_add(tempoccocc2(ispin), tempoccocc1(ispin), 1.0_dp, bfgs_rho)
2513 approx_inv_hessian(ispin), tempoccocc3(ispin))
2514 CALL dbcsr_add(tempoccocc2(ispin), tempoccocc3(ispin), &
2515 1.0_dp, bfgs_rho*bfgs_rho*bfgs_sum)
2519 approx_inv_hessian(ispin), tempoccocc1(ispin))
2521 CALL dbcsr_add(tempoccocc2(ispin), tempoccocc3(ispin), &
2522 1.0_dp, -2.0_dp*bfgs_rho)
2524 CALL dbcsr_copy(approx_inv_hessian(ispin), tempoccocc2(ispin))
2528 ELSE IF (l_bfgs)
THEN
2536 IF (.NOT. l_bfgs)
THEN
2538 DO ispin = 1, nspins
2541 grad(ispin), step(ispin))
2551 IF (iteration == 0)
THEN
2552 reset_conjugator = .true.
2556 IF (.NOT. reset_conjugator)
THEN
2557 CALL compute_cg_beta( &
2559 reset_conjugator=reset_conjugator, &
2560 conjugator=optimizer%conjugator, &
2562 prev_grad=prev_grad(:), &
2564 prev_step=prev_step(:), &
2565 prev_minus_prec_grad=prev_minus_prec_grad(:) &
2570 IF (reset_conjugator)
THEN
2573 IF (unit_nr > 0 .AND. (.NOT. just_started))
THEN
2574 WRITE (unit_nr,
'(T2,A35)')
"Re-setting conjugator to zero"
2576 reset_conjugator = .false.
2581 DO ispin = 1, nspins
2583 CALL dbcsr_copy(prev_minus_prec_grad(ispin), step(ispin))
2586 CALL dbcsr_add(step(ispin), prev_step(ispin), 1.0_dp, beta)
2593 IF (.NOT. line_search)
THEN
2599 DO ispin = 1, nspins
2600 CALL dbcsr_dot(grad(ispin), step(ispin), tempreal)
2603 g0sign = sign(1.0_dp, g0)
2604 IF (linear_search_type == 1)
THEN
2605 IF (iteration == 0)
THEN
2606 step_size = optimizer%lin_search_step_size_guess
2608 IF (next_step_size_guess <= 0.0_dp)
THEN
2609 step_size = optimizer%lin_search_step_size_guess
2612 step_size = optimizer%lin_search_step_size_guess
2616 ELSE IF (linear_search_type == 2)
THEN
2619 step_size = optimizer%lin_search_step_size_guess
2621 IF (unit_nr > 0)
THEN
2622 WRITE (unit_nr,
'(T21,3A19)')
"Line position",
"Line grad",
"Next line step"
2623 WRITE (unit_nr,
'(T2,A19,3F19.5)')
"Line search", 0.0_dp, g0, step_size
2625 next_step_size_guess = step_size
2629 DO ispin = 1, nspins
2630 CALL dbcsr_dot(grad(ispin), step(ispin), tempreal)
2633 g1sign = sign(1.0_dp, g1)
2634 IF (linear_search_type == 1)
THEN
2637 appr_sec_der = (g1 - g0)/step_size
2638 step_size = -g1/appr_sec_der
2639 ELSE IF (linear_search_type == 2)
THEN
2642 IF (g1sign /= g0sign)
THEN
2643 step_size = -step_size/2.0_dp
2645 step_size = step_size*1.5_dp
2649 IF (unit_nr > 0)
THEN
2650 WRITE (unit_nr,
'(T21,3A19)')
"Line position",
"Line grad",
"Next line step"
2651 WRITE (unit_nr,
'(T2,A19,3F19.5)')
"Line search", next_step_size_guess, g1, step_size
2656 next_step_size_guess = next_step_size_guess + step_size
2660 DO ispin = 1, nspins
2661 IF (.NOT. line_search)
THEN
2663 CALL dbcsr_copy(prev_m_theta(ispin), m_theta(ispin))
2665 CALL dbcsr_add(m_theta(ispin), step(ispin), 1.0_dp, step_size)
2670 IF (line_search)
THEN
2677 IF (unit_nr > 0)
THEN
2678 iter_type = trim(
"NLMO OPT "//iter_type)
2679 WRITE (unit_nr,
'(T2,A13,I6,F23.10,E14.5,F14.9,F9.2)') &
2680 iter_type, iteration, &
2681 objf_new, objf_diff, grad_norm, &
2683 WRITE (unit_nr,
'(T2,A19,F23.10)') &
2684 "Localization:", localization_obj_function
2685 WRITE (unit_nr,
'(T2,A19,F23.10)') &
2686 "Orthogonalization:", penalty_func_new
2690 iteration = iteration + 1
2691 IF (prepare_to_exit)
EXIT
2695 IF (converged .OR. (outer_iteration >= outer_max_iter))
THEN
2696 outer_prepare_to_exit = .true.
2699 outer_iteration = outer_iteration + 1
2700 IF (outer_prepare_to_exit)
EXIT
2705 optimizer%opt_penalty%penalty_strength = 0.0_dp
2706 DO ispin = 1, nspins
2707 optimizer%opt_penalty%penalty_strength = optimizer%opt_penalty%penalty_strength + &
2708 (-1.0_dp)*penalty_vol_prefactor(ispin)
2710 optimizer%opt_penalty%penalty_strength = optimizer%opt_penalty%penalty_strength/nspins
2715 iter_type =
"Unconverged"
2718 IF (unit_nr > 0)
THEN
2719 WRITE (unit_nr,
'()')
2720 print_string = trim(iter_type)//
" localization:"
2721 WRITE (unit_nr,
'(T2,A29,F30.10)') &
2722 print_string, localization_obj_function
2723 print_string = trim(iter_type)//
" determinant:"
2724 WRITE (unit_nr,
'(T2,A29,F30.10)') &
2725 print_string, overlap_determinant
2726 print_string = trim(iter_type)//
" penalty strength:"
2727 WRITE (unit_nr,
'(T2,A29,F30.10)') &
2728 print_string, optimizer%opt_penalty%penalty_strength
2735 DO ispin = 1, nspins
2736 DO idim0 = 1,
SIZE(m_b0, 2)
2737 DO reim = 1,
SIZE(m_b0, 1)
2763 DEALLOCATE (grad_norm_spin)
2765 DEALLOCATE (penalty_vol_prefactor)
2766 DEALLOCATE (suggested_vol_penalty)
2768 DEALLOCATE (approx_inv_hessian)
2769 DEALLOCATE (prev_m_theta)
2770 DEALLOCATE (m_theta_normalized)
2772 DEALLOCATE (prev_grad)
2774 DEALLOCATE (prev_step)
2776 DEALLOCATE (prev_minus_prec_grad)
2777 DEALLOCATE (m_sig_sqrti_ii)
2778 DEALLOCATE (m_sigma)
2779 DEALLOCATE (m_siginv)
2780 DEALLOCATE (tempnocc1)
2781 DEALLOCATE (tempoccocc1)
2782 DEALLOCATE (tempoccocc2)
2783 DEALLOCATE (tempoccocc3)
2787 DEALLOCATE (m_theta, m_t_mo_local)
2789 DEALLOCATE (weights)
2790 DEALLOCATE (first_sgf, last_sgf, nsgf)
2792 IF (.NOT. converged)
THEN
2793 cpabort(
"Optimization not converged! ")
2796 CALL timestop(handle)
2818 SUBROUTINE xalmo_analysis(detailed_analysis, eps_filter, m_T_in, m_T0_in, &
2819 m_siginv_in, m_siginv0_in, m_S_in, m_KS0_in, m_quench_t_in, energy_out, &
2820 m_eda_out, m_cta_out)
2822 LOGICAL,
INTENT(IN) :: detailed_analysis
2823 REAL(kind=
dp),
INTENT(IN) :: eps_filter
2824 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(IN) :: m_t_in, m_t0_in, m_siginv_in, &
2825 m_siginv0_in, m_s_in, m_ks0_in, &
2827 REAL(kind=
dp),
INTENT(INOUT) :: energy_out
2828 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(INOUT) :: m_eda_out, m_cta_out
2830 CHARACTER(len=*),
PARAMETER :: routinen =
'xalmo_analysis'
2832 INTEGER :: handle, ispin, nspins
2833 REAL(kind=
dp) :: energy_ispin, spin_factor
2834 TYPE(
dbcsr_type) :: ftsiginv0, fvo0, m_x, siginvtftsiginv0, &
2837 CALL timeset(routinen, handle)
2839 nspins =
SIZE(m_t_in)
2841 IF (nspins == 1)
THEN
2842 spin_factor = 2.0_dp
2844 spin_factor = 1.0_dp
2848 DO ispin = 1, nspins
2852 template=m_t_in(ispin), &
2853 matrix_type=dbcsr_type_no_symmetry)
2855 template=m_t_in(ispin), &
2856 matrix_type=dbcsr_type_no_symmetry)
2858 template=m_t_in(ispin), &
2859 matrix_type=dbcsr_type_no_symmetry)
2861 template=m_t_in(ispin), &
2862 matrix_type=dbcsr_type_no_symmetry)
2864 template=m_siginv0_in(ispin), &
2865 matrix_type=dbcsr_type_no_symmetry)
2868 CALL compute_frequently_used_matrices( &
2869 filter_eps=eps_filter, &
2870 m_t_in=m_t0_in(ispin), &
2871 m_siginv_in=m_siginv0_in(ispin), &
2873 m_f_in=m_ks0_in(ispin), &
2874 m_ftsiginv_out=ftsiginv0, &
2875 m_siginvtftsiginv_out=siginvtftsiginv0, &
2878 CALL dbcsr_copy(fvo0, ftsiginv0, keep_sparsity=.true.)
2883 retain_sparsity=.true.)
2887 CALL dbcsr_add(m_x, m_t_in(ispin), -1.0_dp, 1.0_dp)
2890 energy_out = energy_out + energy_ispin*spin_factor
2892 IF (detailed_analysis)
THEN
2902 m_siginv0_in(ispin), &
2903 0.0_dp, ftsiginv0, &
2904 filter_eps=eps_filter)
2910 filter_eps=eps_filter)
2915 0.0_dp, siginvtftsiginv0, &
2916 filter_eps=eps_filter)
2923 filter_eps=eps_filter)
2927 m_siginv_in(ispin), &
2928 0.0_dp, ftsiginv0, &
2929 filter_eps=eps_filter)
2932 ftsiginv0, m_cta_out(ispin))
2946 CALL timestop(handle)
2948 END SUBROUTINE xalmo_analysis
2965 SUBROUTINE compute_frequently_used_matrices(filter_eps, &
2966 m_T_in, m_siginv_in, m_S_in, m_F_in, m_FTsiginv_out, &
2967 m_siginvTFTsiginv_out, m_ST_out)
2969 REAL(kind=
dp),
INTENT(IN) :: filter_eps
2970 TYPE(
dbcsr_type),
INTENT(IN) :: m_t_in, m_siginv_in, m_s_in, m_f_in
2971 TYPE(
dbcsr_type),
INTENT(INOUT) :: m_ftsiginv_out, m_siginvtftsiginv_out, &
2974 CHARACTER(len=*),
PARAMETER :: routinen =
'compute_frequently_used_matrices'
2979 CALL timeset(routinen, handle)
2983 matrix_type=dbcsr_type_no_symmetry)
2985 template=m_siginv_in, &
2986 matrix_type=dbcsr_type_no_symmetry)
2991 0.0_dp, m_tmp_no_1, &
2992 filter_eps=filter_eps)
2997 0.0_dp, m_ftsiginv_out, &
2998 filter_eps=filter_eps)
3003 0.0_dp, m_tmp_oo_1, &
3004 filter_eps=filter_eps)
3009 0.0_dp, m_siginvtftsiginv_out, &
3010 filter_eps=filter_eps)
3016 filter_eps=filter_eps)
3021 CALL timestop(handle)
3023 END SUBROUTINE compute_frequently_used_matrices
3033 SUBROUTINE split_v_blk(almo_scf_env)
3037 CHARACTER(len=*),
PARAMETER :: routinen =
'split_v_blk'
3039 INTEGER :: discarded_v, handle, iblock_col, &
3040 iblock_col_size, iblock_row, &
3041 iblock_row_size, ispin, retained_v
3042 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: data_p
3045 CALL timeset(routinen, handle)
3047 DO ispin = 1, almo_scf_env%nspins
3050 work_mutable=.true.)
3052 work_mutable=.true.)
3059 row_size=iblock_row_size, col_size=iblock_col_size)
3061 IF (iblock_row /= iblock_col)
THEN
3062 cpabort(
"off-diagonal block found")
3065 retained_v = almo_scf_env%nvirt_of_domain(iblock_col, ispin)
3066 discarded_v = almo_scf_env%nvirt_disc_of_domain(iblock_col, ispin)
3067 cpassert(retained_v > 0)
3068 cpassert(discarded_v > 0)
3069 CALL dbcsr_put_block(almo_scf_env%matrix_v_disc_blk(ispin), iblock_row, iblock_col, &
3070 block=data_p(:, (retained_v + 1):iblock_col_size))
3071 CALL dbcsr_put_block(almo_scf_env%matrix_v_blk(ispin), iblock_row, iblock_col, &
3072 block=data_p(:, 1:retained_v))
3082 CALL timestop(handle)
3084 END SUBROUTINE split_v_blk
3093 SUBROUTINE harris_foulkes_correction(almo_scf_env)
3097 CHARACTER(len=*),
PARAMETER :: routinen =
'harris_foulkes_correction'
3098 INTEGER,
PARAMETER :: cayley_transform = 1, dm_ls_step = 2
3100 INTEGER :: algorithm_id, handle, handle1, handle2, handle3, handle4, handle5, handle6, &
3101 handle7, handle8, ispin, iteration, n, nmins, nspin, opt_k_max_iter, &
3102 outer_opt_k_iteration, outer_opt_k_max_iter, unit_nr
3103 INTEGER,
DIMENSION(1) :: fake, nelectron_spin_real
3104 LOGICAL :: converged, line_search, md_in_k_space, outer_opt_k_prepare_to_exit, &
3105 prepare_to_exit, reset_conjugator, reset_step_size, use_cubic_approximation, &
3106 use_quadratic_approximation
3107 REAL(kind=
dp) :: aa, bb, beta, conjugacy_error, conjugacy_error_threshold, &
3108 delta_obj_function, denom, energy_correction_final, frob_matrix, frob_matrix_base, fun0, &
3109 fun1, gfun0, gfun1, grad_norm, grad_norm_frob, kappa, kin_energy, line_search_error, &
3110 line_search_error_threshold, num_threshold, numer, obj_function, quadratic_approx_error, &
3111 quadratic_approx_error_threshold, safety_multiplier, spin_factor, step_size, &
3112 step_size_quadratic_approx, step_size_quadratic_approx2, t1, t1a, t1cholesky, t2, t2a, &
3113 t2cholesky, tau, time_step, x_opt_eps_adaptive, x_opt_eps_adaptive_factor
3114 REAL(kind=
dp),
DIMENSION(1) :: local_mu
3115 REAL(kind=
dp),
DIMENSION(2) :: energy_correction
3116 REAL(kind=
dp),
DIMENSION(3) :: minima
3119 TYPE(
dbcsr_type) :: grad, k_vd_index_down, k_vr_index_down, matrix_k_central, matrix_tmp1, &
3120 matrix_tmp2, prec, prev_grad, prev_minus_prec_grad, prev_step, sigma_oo_curr, &
3121 sigma_oo_curr_inv, sigma_vv_sqrt, sigma_vv_sqrt_guess, sigma_vv_sqrt_inv, &
3122 sigma_vv_sqrt_inv_guess, step, t_curr, tmp1_n_vr, tmp2_n_o, tmp3_vd_vr, tmp4_o_vr, &
3123 tmp_k_blk, vd_fixed, vd_index_sqrt, vd_index_sqrt_inv, velocity, vr_fixed, vr_index_sqrt, &
3125 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:) :: matrix_p_almo_scf_converged
3127 CALL timeset(routinen, handle)
3131 IF (logger%para_env%is_source())
THEN
3137 nspin = almo_scf_env%nspins
3138 energy_correction_final = 0.0_dp
3139 IF (nspin == 1)
THEN
3140 spin_factor = 2.0_dp
3142 spin_factor = 1.0_dp
3145 IF (almo_scf_env%deloc_use_occ_orbs)
THEN
3146 algorithm_id = cayley_transform
3148 algorithm_id = dm_ls_step
3153 SELECT CASE (algorithm_id)
3154 CASE (cayley_transform)
3158 IF (almo_scf_env%nspins == 1)
THEN
3159 CALL dbcsr_scale(almo_scf_env%matrix_p(1), 1.0_dp/spin_factor)
3165 CALL dbcsr_copy(almo_scf_env%matrix_t(ispin), &
3166 almo_scf_env%matrix_t_blk(ispin))
3173 IF (unit_nr > 0)
THEN
3174 WRITE (unit_nr, *)
"sqrt and inv(sqrt) of MO overlap matrix"
3176 CALL dbcsr_create(almo_scf_env%matrix_sigma_sqrt(ispin), &
3177 template=almo_scf_env%matrix_sigma(ispin), &
3178 matrix_type=dbcsr_type_no_symmetry)
3179 CALL dbcsr_create(almo_scf_env%matrix_sigma_sqrt_inv(ispin), &
3180 template=almo_scf_env%matrix_sigma(ispin), &
3181 matrix_type=dbcsr_type_no_symmetry)
3184 almo_scf_env%matrix_sigma_sqrt_inv(ispin), &
3185 almo_scf_env%matrix_sigma(ispin), &
3186 threshold=almo_scf_env%eps_filter, &
3187 order=almo_scf_env%order_lanczos, &
3188 eps_lanczos=almo_scf_env%eps_lanczos, &
3189 max_iter_lanczos=almo_scf_env%max_iter_lanczos)
3192 CALL dbcsr_create(matrix_tmp1, template=almo_scf_env%matrix_sigma(ispin), &
3193 matrix_type=dbcsr_type_no_symmetry)
3194 CALL dbcsr_create(matrix_tmp2, template=almo_scf_env%matrix_sigma(ispin), &
3195 matrix_type=dbcsr_type_no_symmetry)
3197 CALL dbcsr_multiply(
"N",
"N", 1.0_dp, almo_scf_env%matrix_sigma_sqrt_inv(ispin), &
3198 almo_scf_env%matrix_sigma(ispin), &
3199 0.0_dp, matrix_tmp1, filter_eps=almo_scf_env%eps_filter)
3201 almo_scf_env%matrix_sigma_sqrt_inv(ispin), &
3202 0.0_dp, matrix_tmp2, filter_eps=almo_scf_env%eps_filter)
3207 IF (unit_nr > 0)
THEN
3208 WRITE (unit_nr, *)
"Error for (inv(sqrt(SIG))*SIG*inv(sqrt(SIG))-I)", frob_matrix/frob_matrix_base
3216 IF (almo_scf_env%almo_update_algorithm ==
almo_scf_diag)
THEN
3222 line_search_error_threshold = almo_scf_env%real01
3223 conjugacy_error_threshold = almo_scf_env%real02
3224 quadratic_approx_error_threshold = almo_scf_env%real03
3225 x_opt_eps_adaptive_factor = almo_scf_env%real04
3228 outer_opt_k_max_iter = almo_scf_env%opt_k_outer_max_iter
3229 outer_opt_k_prepare_to_exit = .false.
3230 outer_opt_k_iteration = 0
3232 grad_norm_frob = 0.0_dp
3233 CALL dbcsr_set(almo_scf_env%matrix_x(ispin), 0.0_dp)
3234 IF (almo_scf_env%deloc_truncate_virt ==
virt_full) outer_opt_k_max_iter = 0
3240 psi_out=almo_scf_env%matrix_v(ispin), &
3241 psi_projector=almo_scf_env%matrix_t_blk(ispin), &
3242 metric=almo_scf_env%matrix_s(1), &
3243 project_out=.true., &
3244 psi_projector_orthogonal=.false., &
3245 proj_in_template=almo_scf_env%matrix_ov(ispin), &
3246 eps_filter=almo_scf_env%eps_filter, &
3247 sig_inv_projector=almo_scf_env%matrix_sigma_inv(ispin))
3251 template=almo_scf_env%matrix_v(ispin))
3252 CALL dbcsr_copy(vr_fixed, almo_scf_env%matrix_v(ispin))
3256 template=almo_scf_env%matrix_sigma_vv(ispin), &
3257 matrix_type=dbcsr_type_no_symmetry)
3259 template=almo_scf_env%matrix_sigma_vv(ispin), &
3260 matrix_type=dbcsr_type_no_symmetry)
3262 template=almo_scf_env%matrix_sigma_vv(ispin), &
3263 matrix_type=dbcsr_type_no_symmetry)
3265 template=almo_scf_env%matrix_sigma_vv(ispin), &
3266 matrix_type=dbcsr_type_no_symmetry)
3267 CALL dbcsr_set(sigma_vv_sqrt_guess, 0.0_dp)
3269 CALL dbcsr_filter(sigma_vv_sqrt_guess, almo_scf_env%eps_filter)
3270 CALL dbcsr_set(sigma_vv_sqrt_inv_guess, 0.0_dp)
3272 CALL dbcsr_filter(sigma_vv_sqrt_inv_guess, almo_scf_env%eps_filter)
3275 IF (almo_scf_env%deloc_truncate_virt /=
virt_full)
THEN
3283 psi_out=almo_scf_env%matrix_v_disc(ispin), &
3284 psi_projector=almo_scf_env%matrix_t_blk(ispin), &
3285 metric=almo_scf_env%matrix_s(1), &
3286 project_out=.true., &
3287 psi_projector_orthogonal=.false., &
3288 proj_in_template=almo_scf_env%matrix_ov_disc(ispin), &
3289 eps_filter=almo_scf_env%eps_filter, &
3290 sig_inv_projector=almo_scf_env%matrix_sigma_inv(ispin))
3295 template=almo_scf_env%matrix_v_disc(ispin))
3296 CALL dbcsr_copy(vd_fixed, almo_scf_env%matrix_v_disc(ispin))
3300 template=almo_scf_env%matrix_sigma_vv_blk(ispin), &
3301 matrix_type=dbcsr_type_no_symmetry)
3305 template=almo_scf_env%matrix_vv_disc_blk(ispin), &
3306 matrix_type=dbcsr_type_no_symmetry)
3312 template=almo_scf_env%matrix_k_blk(ispin))
3313 CALL dbcsr_copy(grad, almo_scf_env%matrix_k_blk(ispin))
3316 md_in_k_space = almo_scf_env%logical01
3317 IF (md_in_k_space)
THEN
3319 template=almo_scf_env%matrix_k_blk(ispin))
3320 CALL dbcsr_copy(velocity, almo_scf_env%matrix_k_blk(ispin))
3322 time_step = almo_scf_env%opt_k_trial_step_size
3326 template=almo_scf_env%matrix_k_blk(ispin))
3329 template=almo_scf_env%matrix_k_blk(ispin))
3333 template=almo_scf_env%matrix_k_blk(ispin))
3334 CALL dbcsr_copy(prec, almo_scf_env%matrix_k_blk(ispin))
3338 CALL dbcsr_set(almo_scf_env%matrix_k_blk(ispin), 0.0_dp)
3342 template=almo_scf_env%matrix_k_blk(ispin))
3344 almo_scf_env%matrix_k_blk(ispin))
3346 template=almo_scf_env%matrix_k_blk(ispin))
3348 template=almo_scf_env%matrix_k_blk(ispin))
3351 template=almo_scf_env%matrix_t(ispin))
3353 template=almo_scf_env%matrix_sigma(ispin), &
3354 matrix_type=dbcsr_type_no_symmetry)
3356 template=almo_scf_env%matrix_sigma(ispin), &
3357 matrix_type=dbcsr_type_no_symmetry)
3359 template=almo_scf_env%matrix_v(ispin))
3361 template=almo_scf_env%matrix_k_blk(ispin))
3363 template=almo_scf_env%matrix_t(ispin))
3365 template=almo_scf_env%matrix_ov(ispin))
3367 template=almo_scf_env%matrix_k_blk(ispin))
3373 opt_k_max_iter = almo_scf_env%opt_k_max_iter
3376 prepare_to_exit = .false.
3378 line_search = .false.
3379 obj_function = 0.0_dp
3380 conjugacy_error = 0.0_dp
3381 line_search_error = 0.0_dp
3386 step_size_quadratic_approx = 0.0_dp
3387 reset_step_size = .true.
3388 IF (almo_scf_env%deloc_truncate_virt ==
virt_full) opt_k_max_iter = 0
3393 CALL timeset(
'k_opt_vr', handle1)
3395 IF (almo_scf_env%deloc_truncate_virt /=
virt_full)
THEN
3399 almo_scf_env%matrix_k_blk(ispin), &
3400 0.0_dp, almo_scf_env%matrix_v(ispin), &
3401 filter_eps=almo_scf_env%eps_filter)
3402 CALL dbcsr_add(almo_scf_env%matrix_v(ispin), vr_fixed, &
3407 CALL get_overlap(bra=almo_scf_env%matrix_v(ispin), &
3408 ket=almo_scf_env%matrix_v(ispin), &
3409 overlap=almo_scf_env%matrix_sigma_vv(ispin), &
3410 metric=almo_scf_env%matrix_s(1), &
3411 retain_overlap_sparsity=.false., &
3412 eps_filter=almo_scf_env%eps_filter)
3415 IF (almo_scf_env%deloc_truncate_virt ==
virt_full)
THEN
3416 CALL timeset(
'cholesky', handle2)
3422 template=almo_scf_env%matrix_sigma_vv(ispin), &
3423 matrix_type=dbcsr_type_no_symmetry)
3427 para_env=almo_scf_env%para_env, &
3428 blacs_env=almo_scf_env%blacs_env)
3429 CALL make_triu(sigma_vv_sqrt)
3430 CALL dbcsr_filter(sigma_vv_sqrt, almo_scf_env%eps_filter)
3433 CALL dbcsr_create(matrix_tmp1, template=almo_scf_env%matrix_sigma_vv(ispin), &
3434 matrix_type=dbcsr_type_no_symmetry)
3438 sigma_vv_sqrt_inv, op=
"SOLVE", pos=
"RIGHT", &
3439 para_env=almo_scf_env%para_env, &
3440 blacs_env=almo_scf_env%blacs_env)
3441 CALL dbcsr_filter(sigma_vv_sqrt_inv, almo_scf_env%eps_filter)
3444 CALL dbcsr_create(matrix_tmp1, template=almo_scf_env%matrix_sigma_vv(ispin), &
3445 matrix_type=dbcsr_type_no_symmetry)
3450 -1.0_dp, matrix_tmp1, filter_eps=almo_scf_env%eps_filter)
3454 IF (unit_nr > 0)
THEN
3455 WRITE (unit_nr, *)
"Error for ( U^T * U - Sig )", &
3456 frob_matrix/frob_matrix_base
3460 0.0_dp, matrix_tmp1, filter_eps=almo_scf_env%eps_filter)
3464 IF (unit_nr > 0)
THEN
3465 WRITE (unit_nr, *)
"Error for ( inv(U) * U - I )", &
3466 frob_matrix/frob_matrix_base
3471 IF (unit_nr > 0)
THEN
3472 WRITE (unit_nr, *)
"Cholesky+inverse wall-time: ", t2cholesky - t1cholesky
3474 CALL timestop(handle2)
3477 sigma_vv_sqrt_inv, &
3478 almo_scf_env%matrix_sigma_vv(ispin), &
3479 threshold=almo_scf_env%eps_filter, &
3480 order=almo_scf_env%order_lanczos, &
3481 eps_lanczos=almo_scf_env%eps_lanczos, &
3482 max_iter_lanczos=almo_scf_env%max_iter_lanczos)
3483 CALL dbcsr_copy(sigma_vv_sqrt_inv_guess, sigma_vv_sqrt_inv)
3484 CALL dbcsr_copy(sigma_vv_sqrt_guess, sigma_vv_sqrt)
3486 CALL dbcsr_create(matrix_tmp1, template=almo_scf_env%matrix_sigma_vv(ispin), &
3487 matrix_type=dbcsr_type_no_symmetry)
3488 CALL dbcsr_create(matrix_tmp2, template=almo_scf_env%matrix_sigma_vv(ispin), &
3489 matrix_type=dbcsr_type_no_symmetry)
3492 almo_scf_env%matrix_sigma_vv(ispin), &
3493 0.0_dp, matrix_tmp1, filter_eps=almo_scf_env%eps_filter)
3495 sigma_vv_sqrt_inv, &
3496 0.0_dp, matrix_tmp2, filter_eps=almo_scf_env%eps_filter)
3501 IF (unit_nr > 0)
THEN
3502 WRITE (unit_nr, *)
"Error for (inv(sqrt(SIGVV))*SIGVV*inv(sqrt(SIGVV))-I)", &
3503 frob_matrix/frob_matrix_base
3510 CALL timestop(handle1)
3514 IF ((iteration == 0) .AND. (.NOT. line_search) .AND. &
3515 (outer_opt_k_iteration == 0))
THEN
3516 x_opt_eps_adaptive = &
3517 almo_scf_env%deloc_cayley_eps_convergence
3519 x_opt_eps_adaptive = &
3520 max(abs(almo_scf_env%deloc_cayley_eps_convergence), &
3521 abs(x_opt_eps_adaptive_factor*grad_norm))
3525 para_env=almo_scf_env%para_env, &
3526 blacs_env=almo_scf_env%blacs_env, &
3527 use_occ_orbs=.true., &
3528 use_virt_orbs=.true., &
3529 occ_orbs_orthogonal=.false., &
3530 virt_orbs_orthogonal=.false., &
3531 pp_preconditioner_full=almo_scf_env%deloc_cayley_occ_precond, &
3532 qq_preconditioner_full=almo_scf_env%deloc_cayley_vir_precond, &
3533 tensor_type=almo_scf_env%deloc_cayley_tensor_type, &
3534 neglect_quadratic_term=almo_scf_env%deloc_cayley_linear, &
3535 conjugator=almo_scf_env%deloc_cayley_conjugator, &
3536 max_iter=almo_scf_env%deloc_cayley_max_iter, &
3537 calculate_energy_corr=.true., &
3540 eps_convergence=x_opt_eps_adaptive, &
3541 eps_filter=almo_scf_env%eps_filter, &
3543 q_index_up=sigma_vv_sqrt_inv, &
3544 q_index_down=sigma_vv_sqrt, &
3545 p_index_up=almo_scf_env%matrix_sigma_sqrt_inv(ispin), &
3546 p_index_down=almo_scf_env%matrix_sigma_sqrt(ispin), &
3547 matrix_ks=almo_scf_env%matrix_ks_0deloc(ispin), &
3548 matrix_t=almo_scf_env%matrix_t(ispin), &
3549 matrix_qp_template=almo_scf_env%matrix_vo(ispin), &
3550 matrix_pq_template=almo_scf_env%matrix_ov(ispin), &
3551 matrix_v=almo_scf_env%matrix_v(ispin), &
3552 matrix_x_guess=almo_scf_env%matrix_x(ispin))
3557 energy_correction=energy_correction(ispin), &
3558 copy_matrix_x=almo_scf_env%matrix_x(ispin))
3562 energy_correction(1) = energy_correction(1)*spin_factor
3564 IF (opt_k_max_iter /= 0)
THEN
3566 CALL timeset(
'k_opt_t_curr', handle3)
3570 almo_scf_env%matrix_v(ispin), &
3571 almo_scf_env%matrix_x(ispin), &
3573 filter_eps=almo_scf_env%eps_filter)
3574 CALL dbcsr_add(t_curr, almo_scf_env%matrix_t_blk(ispin), &
3580 overlap=sigma_oo_curr, &
3581 metric=almo_scf_env%matrix_s(1), &
3582 retain_overlap_sparsity=.false., &
3583 eps_filter=almo_scf_env%eps_filter)
3584 IF (iteration == 0)
THEN
3587 threshold=almo_scf_env%eps_filter, &
3588 use_inv_as_guess=.false.)
3592 threshold=almo_scf_env%eps_filter, &
3593 use_inv_as_guess=.true.)
3596 CALL dbcsr_create(matrix_tmp1, template=sigma_oo_curr, &
3597 matrix_type=dbcsr_type_no_symmetry)
3599 sigma_oo_curr_inv, &
3600 0.0_dp, matrix_tmp1, &
3601 filter_eps=almo_scf_env%eps_filter)
3605 IF (unit_nr > 0)
THEN
3606 WRITE (unit_nr, *)
"Error for (SIG*inv(SIG)-I)", &
3607 frob_matrix/frob_matrix_base, frob_matrix_base
3612 CALL dbcsr_create(matrix_tmp1, template=sigma_oo_curr, &
3613 matrix_type=dbcsr_type_no_symmetry)
3616 0.0_dp, matrix_tmp1, &
3617 filter_eps=almo_scf_env%eps_filter)
3621 IF (unit_nr > 0)
THEN
3622 WRITE (unit_nr, *)
"Error for (inv(SIG)*SIG-I)", &
3623 frob_matrix/frob_matrix_base, frob_matrix_base
3628 CALL timestop(handle3)
3629 CALL timeset(
'k_opt_vd', handle4)
3637 sigma_vv_sqrt_inv, &
3638 sigma_vv_sqrt_inv, &
3639 0.0_dp, sigma_vv_sqrt, &
3640 filter_eps=almo_scf_env%eps_filter)
3642 psi_out=almo_scf_env%matrix_v_disc(ispin), &
3643 psi_projector=almo_scf_env%matrix_v(ispin), &
3644 metric=almo_scf_env%matrix_s(1), &
3645 project_out=.false., &
3646 psi_projector_orthogonal=.false., &
3647 proj_in_template=almo_scf_env%matrix_k_tr(ispin), &
3648 eps_filter=almo_scf_env%eps_filter, &
3649 sig_inv_projector=sigma_vv_sqrt)
3651 CALL dbcsr_add(almo_scf_env%matrix_v_disc(ispin), &
3652 vd_fixed, -1.0_dp, +1.0_dp)
3654 CALL timestop(handle4)
3655 CALL timeset(
'k_opt_grad', handle5)
3660 IF (line_search)
THEN
3664 almo_scf_env%matrix_ks_0deloc(ispin), &
3667 filter_eps=almo_scf_env%eps_filter)
3669 sigma_oo_curr_inv, &
3670 almo_scf_env%matrix_x(ispin), &
3671 0.0_dp, tmp4_o_vr, &
3672 filter_eps=almo_scf_env%eps_filter)
3676 0.0_dp, tmp1_n_vr, &
3677 filter_eps=almo_scf_env%eps_filter)
3679 almo_scf_env%matrix_v_disc(ispin), &
3682 retain_sparsity=.true.)
3689 converged = (grad_norm < almo_scf_env%opt_k_eps_convergence)
3690 IF (converged .OR. (iteration >= opt_k_max_iter))
THEN
3691 prepare_to_exit = .true.
3693 CALL timestop(handle5)
3695 IF (.NOT. prepare_to_exit)
THEN
3697 CALL timeset(
'k_opt_energy', handle6)
3703 0.0_dp, sigma_oo_curr, &
3704 filter_eps=almo_scf_env%eps_filter)
3705 delta_obj_function = fun0
3706 CALL dbcsr_dot(sigma_oo_curr_inv, sigma_oo_curr, obj_function)
3707 delta_obj_function = obj_function - delta_obj_function
3708 IF (line_search)
THEN
3714 CALL timestop(handle6)
3717 IF (.NOT. line_search)
THEN
3719 CALL timeset(
'k_opt_step', handle7)
3721 IF ((.NOT. md_in_k_space) .AND. &
3722 (iteration >= max(0, almo_scf_env%opt_k_prec_iter_start) .AND. &
3723 mod(iteration - almo_scf_env%opt_k_prec_iter_start, &
3724 almo_scf_env%opt_k_prec_iter_freq) == 0))
THEN
3729 IF (unit_nr > 0)
THEN
3730 WRITE (unit_nr, *)
"Computing preconditioner"
3732 CALL opt_k_create_preconditioner_blk(almo_scf_env, &
3733 almo_scf_env%matrix_v_disc(ispin), &
3745 CALL opt_k_apply_preconditioner_blk(almo_scf_env, &
3750 reset_conjugator = .false.
3752 IF (iteration < max(almo_scf_env%opt_k_conj_iter_start, 1) .OR. &
3753 mod(iteration - almo_scf_env%opt_k_conj_iter_start, &
3754 almo_scf_env%opt_k_conj_iter_freq) == 0)
THEN
3756 reset_conjugator = .true.
3761 CALL dbcsr_dot(grad, prev_minus_prec_grad, numer)
3762 CALL dbcsr_dot(prev_grad, prev_minus_prec_grad, denom)
3763 conjugacy_error = numer/denom
3765 IF (conjugacy_error > min(0.5_dp, conjugacy_error_threshold))
THEN
3766 reset_conjugator = .true.
3767 IF (unit_nr > 0)
THEN
3768 WRITE (unit_nr, *)
"Lack of progress, conjugacy error is ", conjugacy_error
3773 IF ((iteration /= 0) .AND. (.NOT. reset_conjugator))
THEN
3775 CALL dbcsr_dot(prev_grad, prev_step, denom)
3776 line_search_error = numer/denom
3777 IF (line_search_error > line_search_error_threshold)
THEN
3778 reset_conjugator = .true.
3779 IF (unit_nr > 0)
THEN
3780 WRITE (unit_nr, *)
"Bad line search, line search error is ", line_search_error
3788 IF (.NOT. reset_conjugator)
THEN
3790 SELECT CASE (almo_scf_env%opt_k_conjugator)
3793 CALL dbcsr_add(tmp_k_blk, prev_grad, 1.0_dp, -1.0_dp)
3795 CALL dbcsr_dot(tmp_k_blk, prev_step, denom)
3796 beta = -1.0_dp*numer/denom
3799 CALL dbcsr_dot(prev_grad, prev_minus_prec_grad, denom)
3802 CALL dbcsr_dot(prev_grad, prev_minus_prec_grad, denom)
3804 CALL dbcsr_add(tmp_k_blk, prev_grad, 1.0_dp, -1.0_dp)
3809 CALL dbcsr_dot(prev_grad, prev_step, denom)
3812 CALL dbcsr_dot(prev_grad, prev_step, denom)
3814 CALL dbcsr_add(tmp_k_blk, prev_grad, 1.0_dp, -1.0_dp)
3820 CALL dbcsr_add(tmp_k_blk, prev_grad, 1.0_dp, -1.0_dp)
3821 CALL dbcsr_dot(tmp_k_blk, prev_step, denom)
3822 beta = -1.0_dp*numer/denom
3825 CALL dbcsr_add(tmp_k_blk, prev_grad, 1.0_dp, -1.0_dp)
3826 CALL dbcsr_dot(tmp_k_blk, prev_step, denom)
3827 CALL dbcsr_dot(tmp_k_blk, prev_minus_prec_grad, numer)
3828 kappa = -2.0_dp*numer/denom
3830 tau = -1.0_dp*numer/denom
3832 beta = tau - kappa*numer/denom
3836 cpabort(
"illegal conjugator")
3839 IF (beta < 0.0_dp)
THEN
3840 IF (unit_nr > 0)
THEN
3841 WRITE (unit_nr, *)
"Beta is negative, ", beta
3843 reset_conjugator = .true.
3848 IF (md_in_k_space)
THEN
3849 reset_conjugator = .true.
3852 IF (reset_conjugator)
THEN
3856 IF (unit_nr > 0)
THEN
3857 WRITE (unit_nr, *)
"(Re)-setting conjugator to zero"
3866 CALL dbcsr_add(step, prev_step, 1.0_dp, beta)
3868 CALL timestop(handle7)
3872 conjugacy_error = 0.0_dp
3876 IF (line_search)
THEN
3878 line_search_error = gfun1/gfun0
3884 IF (line_search)
THEN
3887 safety_multiplier = 1.0e+1_dp
3888 num_threshold = max(epsilon(1.0_dp), &
3889 safety_multiplier*(almo_scf_env%eps_filter**2)*almo_scf_env%ndomains)
3890 IF (abs(fun1 - fun0 - gfun0*step_size) < num_threshold)
THEN
3891 IF (unit_nr > 0)
THEN
3892 WRITE (unit_nr,
'(T3,A,1X,E17.7)') &
3893 "Numerical accuracy is too low to observe non-linear behavior", &
3894 abs(fun1 - fun0 - gfun0*step_size)
3895 WRITE (unit_nr,
'(T3,A,1X,E17.7,A,1X,E12.3)')
"Error computing ", &
3897 " is smaller than the threshold", num_threshold
3899 cpabort(
"Unable to continue with low numerical accuracy")
3901 IF (abs(gfun0) < num_threshold)
THEN
3902 IF (unit_nr > 0)
THEN
3903 WRITE (unit_nr,
'(T3,A,1X,E17.7,A,1X,E12.3)')
"Linear gradient", &
3905 " is smaller than the threshold", num_threshold
3907 cpabort(
"Unable to continue with low numerical accuracy")
3910 use_quadratic_approximation = .true.
3911 use_cubic_approximation = .false.
3915 step_size_quadratic_approx = -(gfun0*step_size*step_size)/(2.0_dp*(fun1 - fun0 - gfun0*step_size))
3917 step_size_quadratic_approx2 = -(fun1 - fun0 - step_size*gfun1/2.0_dp)/(gfun1 - (fun1 - fun0)/step_size)
3919 IF ((step_size_quadratic_approx < 0.0_dp) .AND. &
3920 (step_size_quadratic_approx2 < 0.0_dp))
THEN
3921 IF (unit_nr > 0)
THEN
3922 WRITE (unit_nr,
'(T3,A,1X,E17.7,1X,E17.7,1X,A)') &
3923 "Quadratic approximation gives negative steps", &
3924 step_size_quadratic_approx, step_size_quadratic_approx2, &
3927 use_cubic_approximation = .true.
3928 use_quadratic_approximation = .false.
3930 IF (step_size_quadratic_approx < 0.0_dp)
THEN
3931 step_size_quadratic_approx = step_size_quadratic_approx2
3933 IF (step_size_quadratic_approx2 < 0.0_dp)
THEN
3934 step_size_quadratic_approx2 = step_size_quadratic_approx
3939 IF (use_quadratic_approximation)
THEN
3940 quadratic_approx_error = abs(step_size_quadratic_approx - &
3941 step_size_quadratic_approx2)/step_size_quadratic_approx
3942 IF (quadratic_approx_error > quadratic_approx_error_threshold)
THEN
3943 IF (unit_nr > 0)
THEN
3944 WRITE (unit_nr,
'(T3,A,1X,E17.7,1X,E17.7,1X,A)')
"Quadratic approximation is poor", &
3945 step_size_quadratic_approx, step_size_quadratic_approx2, &
3946 "Try cubic approximation"
3948 use_cubic_approximation = .true.
3949 use_quadratic_approximation = .false.
3954 IF (use_cubic_approximation)
THEN
3959 bb = (-step_size*gfun1 + 3.0_dp*(fun1 - fun0) - 2.0_dp*step_size*gfun0)/(step_size*step_size)
3960 aa = (gfun1 - 2.0_dp*step_size*bb - gfun0)/(3.0_dp*step_size*step_size)
3962 IF (abs(gfun1 - 2.0_dp*step_size*bb - gfun0) < num_threshold)
THEN
3963 IF (unit_nr > 0)
THEN
3964 WRITE (unit_nr,
'(T3,A,1X,E17.7)') &
3965 "Numerical accuracy is too low to observe cubic behavior", &
3966 abs(gfun1 - 2.0_dp*step_size*bb - gfun0)
3968 use_cubic_approximation = .false.
3969 use_quadratic_approximation = .true.
3971 IF (abs(gfun1) < num_threshold)
THEN
3972 IF (unit_nr > 0)
THEN
3973 WRITE (unit_nr,
'(T3,A,1X,E17.7,A,1X,E12.3)')
"Linear gradient", &
3975 " is smaller than the threshold", num_threshold
3977 use_cubic_approximation = .false.
3978 use_quadratic_approximation = .true.
3983 IF (use_cubic_approximation)
THEN
3988 IF (unit_nr > 0)
THEN
3989 WRITE (unit_nr,
'(T3,A)') &
3990 "Cubic approximation gives zero soultions! Use quadratic approximation"
3992 use_quadratic_approximation = .true.
3993 use_cubic_approximation = .true.
3995 step_size = minima(1)
3997 IF (unit_nr > 0)
THEN
3998 WRITE (unit_nr,
'(T3,A)') &
3999 "More than one solution found! Use quadratic approximation"
4001 use_quadratic_approximation = .true.
4002 use_cubic_approximation = .true.
4007 IF (use_quadratic_approximation)
THEN
4008 IF (unit_nr > 0)
THEN
4009 WRITE (unit_nr,
'(T3,A)')
"Use quadratic approximation"
4011 step_size = (step_size_quadratic_approx + step_size_quadratic_approx2)*0.5_dp
4015 IF (step_size < 0.0_dp)
THEN
4016 cpabort(
"Negative step proposed")
4019 CALL dbcsr_copy(almo_scf_env%matrix_k_blk(ispin), &
4021 CALL dbcsr_add(almo_scf_env%matrix_k_blk(ispin), &
4022 step, 1.0_dp, step_size)
4024 almo_scf_env%matrix_k_blk(ispin))
4025 line_search = .false.
4029 IF (md_in_k_space)
THEN
4032 IF (iteration /= 0)
THEN
4034 step, 1.0_dp, 0.5_dp*time_step)
4036 prev_step, 1.0_dp, 0.5_dp*time_step)
4039 kin_energy = 0.5_dp*kin_energy*kin_energy
4042 CALL dbcsr_add(almo_scf_env%matrix_k_blk(ispin), &
4043 velocity, 1.0_dp, time_step)
4044 CALL dbcsr_add(almo_scf_env%matrix_k_blk(ispin), &
4045 step, 1.0_dp, 0.5_dp*time_step*time_step)
4049 IF (reset_step_size)
THEN
4050 step_size = almo_scf_env%opt_k_trial_step_size
4051 reset_step_size = .false.
4053 step_size = step_size*almo_scf_env%opt_k_trial_step_size_multiplier
4055 CALL dbcsr_copy(almo_scf_env%matrix_k_blk(ispin), &
4057 CALL dbcsr_add(almo_scf_env%matrix_k_blk(ispin), &
4058 step, 1.0_dp, step_size)
4059 line_search = .true.
4068 IF (unit_nr > 0)
THEN
4069 IF (md_in_k_space)
THEN
4070 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)') &
4071 "K iter CG", iteration, time_step, time_step*iteration, &
4072 energy_correction(ispin), obj_function, delta_obj_function, grad_norm, &
4073 kin_energy, kin_energy + obj_function, beta
4075 IF (line_search .OR. prepare_to_exit)
THEN
4076 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)') &
4077 "K iter CG", iteration, step_size, &
4078 energy_correction(ispin), delta_obj_function, grad_norm, &
4079 gfun0, line_search_error, beta, conjugacy_error, t2a - t1a
4081 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)') &
4082 "K iter LS", iteration, step_size, &
4083 energy_correction(ispin), delta_obj_function, grad_norm, &
4084 gfun1, line_search_error, beta, conjugacy_error, t2a - t1a
4092 prepare_to_exit = .true.
4095 IF (.NOT. line_search) iteration = iteration + 1
4097 IF (prepare_to_exit)
EXIT
4101 IF (converged .OR. (outer_opt_k_iteration >= outer_opt_k_max_iter))
THEN
4102 outer_opt_k_prepare_to_exit = .true.
4105 IF (almo_scf_env%deloc_truncate_virt /=
virt_full)
THEN
4107 IF (unit_nr > 0)
THEN
4108 WRITE (unit_nr, *)
"Updating ALMO virtuals"
4111 CALL timeset(
'k_opt_v0_update', handle8)
4115 almo_scf_env%matrix_v_disc_blk(ispin), &
4116 almo_scf_env%matrix_k_blk(ispin), &
4118 filter_eps=almo_scf_env%eps_filter)
4119 CALL dbcsr_add(vr_fixed, almo_scf_env%matrix_v_blk(ispin), &
4124 almo_scf_env%matrix_v_blk(ispin), &
4125 almo_scf_env%matrix_k_blk(ispin), &
4127 filter_eps=almo_scf_env%eps_filter)
4128 CALL dbcsr_add(vd_fixed, almo_scf_env%matrix_v_disc_blk(ispin), &
4134 overlap=k_vr_index_down, &
4135 metric=almo_scf_env%matrix_s_blk(1), &
4136 retain_overlap_sparsity=.false., &
4137 eps_filter=almo_scf_env%eps_filter)
4138 CALL dbcsr_create(vr_index_sqrt_inv, template=k_vr_index_down, &
4139 matrix_type=dbcsr_type_no_symmetry)
4140 CALL dbcsr_create(vr_index_sqrt, template=k_vr_index_down, &
4141 matrix_type=dbcsr_type_no_symmetry)
4143 vr_index_sqrt_inv, &
4145 threshold=almo_scf_env%eps_filter, &
4146 order=almo_scf_env%order_lanczos, &
4147 eps_lanczos=almo_scf_env%eps_lanczos, &
4148 max_iter_lanczos=almo_scf_env%max_iter_lanczos)
4150 CALL dbcsr_create(matrix_tmp1, template=k_vr_index_down, &
4151 matrix_type=dbcsr_type_no_symmetry)
4152 CALL dbcsr_create(matrix_tmp2, template=k_vr_index_down, &
4153 matrix_type=dbcsr_type_no_symmetry)
4157 0.0_dp, matrix_tmp1, filter_eps=almo_scf_env%eps_filter)
4159 vr_index_sqrt_inv, &
4160 0.0_dp, matrix_tmp2, filter_eps=almo_scf_env%eps_filter)
4165 IF (unit_nr > 0)
THEN
4166 WRITE (unit_nr, *)
"Error for (inv(sqrt(SIGVV))*SIGVV*inv(sqrt(SIGVV))-I)", &
4167 frob_matrix/frob_matrix_base
4175 vr_index_sqrt_inv, &
4176 0.0_dp, almo_scf_env%matrix_v_blk(ispin), &
4177 filter_eps=almo_scf_env%eps_filter)
4181 overlap=k_vd_index_down, &
4182 metric=almo_scf_env%matrix_s_blk(1), &
4183 retain_overlap_sparsity=.false., &
4184 eps_filter=almo_scf_env%eps_filter)
4185 CALL dbcsr_create(vd_index_sqrt_inv, template=k_vd_index_down, &
4186 matrix_type=dbcsr_type_no_symmetry)
4187 CALL dbcsr_create(vd_index_sqrt, template=k_vd_index_down, &
4188 matrix_type=dbcsr_type_no_symmetry)
4190 vd_index_sqrt_inv, &
4192 threshold=almo_scf_env%eps_filter, &
4193 order=almo_scf_env%order_lanczos, &
4194 eps_lanczos=almo_scf_env%eps_lanczos, &
4195 max_iter_lanczos=almo_scf_env%max_iter_lanczos)
4197 CALL dbcsr_create(matrix_tmp1, template=k_vd_index_down, &
4198 matrix_type=dbcsr_type_no_symmetry)
4199 CALL dbcsr_create(matrix_tmp2, template=k_vd_index_down, &
4200 matrix_type=dbcsr_type_no_symmetry)
4204 0.0_dp, matrix_tmp1, filter_eps=almo_scf_env%eps_filter)
4206 vd_index_sqrt_inv, &
4207 0.0_dp, matrix_tmp2, filter_eps=almo_scf_env%eps_filter)
4212 IF (unit_nr > 0)
THEN
4213 WRITE (unit_nr, *)
"Error for (inv(sqrt(SIGVV))*SIGVV*inv(sqrt(SIGVV))-I)", &
4214 frob_matrix/frob_matrix_base
4222 vd_index_sqrt_inv, &
4223 0.0_dp, almo_scf_env%matrix_v_disc_blk(ispin), &
4224 filter_eps=almo_scf_env%eps_filter)
4231 CALL timestop(handle8)
4238 IF (almo_scf_env%deloc_truncate_virt /=
virt_full)
THEN
4258 IF (md_in_k_space)
THEN
4264 outer_opt_k_iteration = outer_opt_k_iteration + 1
4265 IF (outer_opt_k_prepare_to_exit)
EXIT
4279 IF (.NOT. almo_scf_env%s_sqrt_done)
THEN
4281 IF (unit_nr > 0)
THEN
4282 WRITE (unit_nr, *)
"sqrt and inv(sqrt) of AO overlap matrix"
4285 template=almo_scf_env%matrix_s(1), &
4286 matrix_type=dbcsr_type_no_symmetry)
4288 template=almo_scf_env%matrix_s(1), &
4289 matrix_type=dbcsr_type_no_symmetry)
4292 almo_scf_env%matrix_s_sqrt_inv(1), &
4293 almo_scf_env%matrix_s(1), &
4294 threshold=almo_scf_env%eps_filter, &
4295 order=almo_scf_env%order_lanczos, &
4296 eps_lanczos=almo_scf_env%eps_lanczos, &
4297 max_iter_lanczos=almo_scf_env%max_iter_lanczos)
4300 CALL dbcsr_create(matrix_tmp1, template=almo_scf_env%matrix_s(1), &
4301 matrix_type=dbcsr_type_no_symmetry)
4302 CALL dbcsr_create(matrix_tmp2, template=almo_scf_env%matrix_s(1), &
4303 matrix_type=dbcsr_type_no_symmetry)
4305 CALL dbcsr_multiply(
"N",
"N", 1.0_dp, almo_scf_env%matrix_s_sqrt_inv(1), &
4306 almo_scf_env%matrix_s(1), &
4307 0.0_dp, matrix_tmp1, filter_eps=almo_scf_env%eps_filter)
4308 CALL dbcsr_multiply(
"N",
"N", 1.0_dp, matrix_tmp1, almo_scf_env%matrix_s_sqrt_inv(1), &
4309 0.0_dp, matrix_tmp2, filter_eps=almo_scf_env%eps_filter)
4314 IF (unit_nr > 0)
THEN
4315 WRITE (unit_nr, *)
"Error for (inv(sqrt(S))*S*inv(sqrt(S))-I)", frob_matrix/frob_matrix_base
4322 almo_scf_env%s_sqrt_done = .true.
4330 para_env=almo_scf_env%para_env, &
4331 blacs_env=almo_scf_env%blacs_env, &
4332 use_occ_orbs=.true., &
4333 use_virt_orbs=almo_scf_env%deloc_cayley_use_virt_orbs, &
4334 occ_orbs_orthogonal=.false., &
4335 virt_orbs_orthogonal=almo_scf_env%orthogonal_basis, &
4336 tensor_type=almo_scf_env%deloc_cayley_tensor_type, &
4337 neglect_quadratic_term=almo_scf_env%deloc_cayley_linear, &
4338 calculate_energy_corr=.true., &
4341 pp_preconditioner_full=almo_scf_env%deloc_cayley_occ_precond, &
4342 qq_preconditioner_full=almo_scf_env%deloc_cayley_vir_precond, &
4343 eps_convergence=almo_scf_env%deloc_cayley_eps_convergence, &
4344 eps_filter=almo_scf_env%eps_filter, &
4346 q_index_up=almo_scf_env%matrix_s_sqrt_inv(1), &
4347 q_index_down=almo_scf_env%matrix_s_sqrt(1), &
4348 p_index_up=almo_scf_env%matrix_sigma_sqrt_inv(ispin), &
4349 p_index_down=almo_scf_env%matrix_sigma_sqrt(ispin), &
4350 matrix_ks=almo_scf_env%matrix_ks_0deloc(ispin), &
4351 matrix_p=almo_scf_env%matrix_p(ispin), &
4352 matrix_qp_template=almo_scf_env%matrix_t(ispin), &
4353 matrix_pq_template=almo_scf_env%matrix_t_tr(ispin), &
4354 matrix_t=almo_scf_env%matrix_t(ispin), &
4355 conjugator=almo_scf_env%deloc_cayley_conjugator, &
4356 max_iter=almo_scf_env%deloc_cayley_max_iter)
4364 energy_correction=energy_correction(ispin))
4370 energy_correction(1) = energy_correction(1)*spin_factor
4377 IF (unit_nr > 0)
THEN
4379 WRITE (unit_nr,
'(T2,A,I6,F20.9)')
"ECORR", ispin, &
4380 energy_correction(ispin)
4383 energy_correction_final = energy_correction_final + energy_correction(ispin)
4388 p=almo_scf_env%matrix_p(ispin), &
4389 eps_filter=almo_scf_env%eps_filter, &
4390 orthog_orbs=.false., &
4391 nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin), &
4392 s=almo_scf_env%matrix_s(1), &
4393 sigma=almo_scf_env%matrix_sigma(ispin), &
4394 sigma_inv=almo_scf_env%matrix_sigma_inv(ispin), &
4396 algorithm=almo_scf_env%sigma_inv_algorithm, &
4397 inverse_accelerator=almo_scf_env%order_lanczos, &
4398 inv_eps_factor=almo_scf_env%matrix_iter_eps_error_factor, &
4399 eps_lanczos=almo_scf_env%eps_lanczos, &
4400 max_iter_lanczos=almo_scf_env%max_iter_lanczos, &
4401 para_env=almo_scf_env%para_env, &
4402 blacs_env=almo_scf_env%blacs_env)
4404 IF (almo_scf_env%nspins == 1)
THEN
4414 IF (.NOT. almo_scf_env%s_inv_done)
THEN
4415 IF (unit_nr > 0)
THEN
4416 WRITE (unit_nr, *)
"Inverting AO overlap matrix"
4419 template=almo_scf_env%matrix_s(1), &
4420 matrix_type=dbcsr_type_no_symmetry)
4421 IF (.NOT. almo_scf_env%s_sqrt_done)
THEN
4423 almo_scf_env%matrix_s(1), &
4424 threshold=almo_scf_env%eps_filter)
4426 CALL dbcsr_multiply(
"N",
"N", 1.0_dp, almo_scf_env%matrix_s_sqrt_inv(1), &
4427 almo_scf_env%matrix_s_sqrt_inv(1), &
4428 0.0_dp, almo_scf_env%matrix_s_inv(1), &
4429 filter_eps=almo_scf_env%eps_filter)
4433 CALL dbcsr_create(matrix_tmp1, template=almo_scf_env%matrix_s(1), &
4434 matrix_type=dbcsr_type_no_symmetry)
4435 CALL dbcsr_multiply(
"N",
"N", 1.0_dp, almo_scf_env%matrix_s_inv(1), &
4436 almo_scf_env%matrix_s(1), &
4437 0.0_dp, matrix_tmp1, &
4438 filter_eps=almo_scf_env%eps_filter)
4442 IF (unit_nr > 0)
THEN
4443 WRITE (unit_nr, *)
"Error for (inv(S)*S-I)", &
4444 frob_matrix/frob_matrix_base
4449 almo_scf_env%s_inv_done = .true.
4453 ALLOCATE (matrix_p_almo_scf_converged(nspin))
4455 CALL dbcsr_create(matrix_p_almo_scf_converged(ispin), &
4456 template=almo_scf_env%matrix_p(ispin))
4457 CALL dbcsr_copy(matrix_p_almo_scf_converged(ispin), &
4458 almo_scf_env%matrix_p(ispin))
4464 nelectron_spin_real(1) = almo_scf_env%nelectrons_spin(ispin)
4465 IF (almo_scf_env%nspins == 1)
THEN
4466 nelectron_spin_real(1) = nelectron_spin_real(1)/2
4469 local_mu(1) = sum(almo_scf_env%mu_of_domain(:, ispin))/almo_scf_env%ndomains
4472 cpabort(
"CVS only: density_matrix_sign has not been updated in SVN")
4474 IF (almo_scf_env%nspins == 1)
THEN
4478 CALL dbcsr_add(matrix_p_almo_scf_converged(ispin), &
4479 almo_scf_env%matrix_p(ispin), -1.0_dp, 1.0_dp)
4480 CALL dbcsr_dot(almo_scf_env%matrix_ks_0deloc(ispin), &
4481 matrix_p_almo_scf_converged(ispin), &
4482 energy_correction(ispin))
4484 energy_correction_final = energy_correction_final + energy_correction(ispin)
4486 IF (unit_nr > 0)
THEN
4488 WRITE (unit_nr,
'(T2,A,I6,F20.9)')
"ECORR", ispin, &
4489 energy_correction(ispin)
4498 DEALLOCATE (matrix_p_almo_scf_converged)
4504 IF (unit_nr > 0)
THEN
4506 WRITE (unit_nr,
'(T2,A,F18.9,F18.9,F18.9,F12.6)')
"ETOT", &
4507 almo_scf_env%almo_scf_energy, &
4508 energy_correction_final, &
4509 almo_scf_env%almo_scf_energy + energy_correction_final, &
4514 CALL timestop(handle)
4516 END SUBROUTINE harris_foulkes_correction
4522 SUBROUTINE make_triu(matrix)
4525 CHARACTER(len=*),
PARAMETER :: routinen =
'make_triu'
4527 INTEGER :: col, handle, i, j, row
4528 REAL(
dp),
DIMENSION(:, :),
POINTER :: block
4531 CALL timeset(routinen, handle)
4536 IF (row > col) block(:, :) = 0.0_dp
4537 IF (row == col)
THEN
4538 DO j = 1,
SIZE(block, 2)
4539 DO i = j + 1,
SIZE(block, 1)
4540 block(i, j) = 0.0_dp
4548 CALL timestop(handle)
4549 END SUBROUTINE make_triu
4571 SUBROUTINE opt_k_create_preconditioner(prec, vd_prop, f, x, oo_inv_x_tr, s, grad, &
4572 vd_blk, t, template_vd_vd_blk, template_vr_vr_blk, template_n_vr, &
4573 spin_factor, eps_filter)
4576 TYPE(
dbcsr_type),
INTENT(IN) :: vd_prop, f, x, oo_inv_x_tr, s
4578 TYPE(
dbcsr_type),
INTENT(IN) :: vd_blk, t, template_vd_vd_blk, &
4579 template_vr_vr_blk, template_n_vr
4580 REAL(kind=
dp),
INTENT(IN) :: spin_factor, eps_filter
4582 CHARACTER(len=*),
PARAMETER :: routinen =
'opt_k_create_preconditioner'
4584 INTEGER :: handle, p_nrows, q_nrows
4585 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: p_diagonal, q_diagonal
4586 TYPE(
dbcsr_type) :: pp_diag, qq_diag, t1, t2, tmp, &
4587 tmp1_n_vr, tmp2_n_vr, tmp_n_vd, &
4588 tmp_vd_vd_blk, tmp_vr_vr_blk
4590 CALL timeset(routinen, handle)
4601 0.0_dp, tmp_n_vd, filter_eps=eps_filter)
4603 template=template_vd_vd_blk)
4604 CALL dbcsr_copy(tmp_vd_vd_blk, template_vd_vd_blk)
4606 0.0_dp, tmp_vd_vd_blk, &
4607 retain_sparsity=.true., &
4608 filter_eps=eps_filter)
4611 ALLOCATE (q_diagonal(q_nrows))
4614 template=template_vd_vd_blk)
4619 0.0_dp, t1, filter_eps=eps_filter)
4622 CALL dbcsr_create(tmp_vr_vr_blk, template=template_vr_vr_blk)
4623 CALL dbcsr_copy(tmp_vr_vr_blk, template_vr_vr_blk)
4625 0.0_dp, tmp_vr_vr_blk, &
4626 retain_sparsity=.true., &
4627 filter_eps=eps_filter)
4630 ALLOCATE (p_diagonal(p_nrows))
4632 CALL dbcsr_create(pp_diag, template=template_vr_vr_blk)
4638 0.0_dp, t2, filter_eps=eps_filter)
4644 0.0_dp, tmp_n_vd, filter_eps=eps_filter)
4646 0.0_dp, tmp_vd_vd_blk, &
4647 retain_sparsity=.true., &
4648 filter_eps=eps_filter)
4655 0.0_dp, t1, filter_eps=eps_filter)
4661 0.0_dp, tmp1_n_vr, filter_eps=eps_filter)
4663 0.0_dp, tmp2_n_vr, filter_eps=eps_filter)
4665 0.0_dp, tmp_vr_vr_blk, &
4666 retain_sparsity=.true., &
4667 filter_eps=eps_filter)
4674 0.0_dp, t2, filter_eps=eps_filter)
4677 CALL dbcsr_add(prec, tmp, 1.0_dp, -1.0_dp)
4682 0.0_dp, tmp_n_vd, filter_eps=eps_filter)
4684 0.0_dp, tmp, retain_sparsity=.true., &
4685 filter_eps=eps_filter)
4690 CALL dbcsr_add(prec, t1, 1.0_dp, 1.0_dp)
4692 CALL inverse_of_elements(prec)
4695 DEALLOCATE (q_diagonal)
4696 DEALLOCATE (p_diagonal)
4708 CALL timestop(handle)
4710 END SUBROUTINE opt_k_create_preconditioner
4725 SUBROUTINE opt_k_create_preconditioner_blk(almo_scf_env, vd_prop, oo_inv_x_tr, &
4726 t_curr, ispin, spin_factor)
4729 TYPE(
dbcsr_type),
INTENT(IN) :: vd_prop, oo_inv_x_tr, t_curr
4730 INTEGER,
INTENT(IN) :: ispin
4731 REAL(kind=
dp),
INTENT(IN) :: spin_factor
4733 CHARACTER(len=*),
PARAMETER :: routinen =
'opt_k_create_preconditioner_blk'
4736 REAL(kind=
dp) :: eps_filter
4737 TYPE(
dbcsr_type) :: opt_k_e_dd, opt_k_e_rr, s_dd_sqrt, &
4738 s_rr_sqrt, t1, tmp, tmp1_n_vr, &
4739 tmp2_n_vr, tmp_n_vd, tmp_vd_vd_blk, &
4744 CALL timeset(routinen, handle)
4746 eps_filter = almo_scf_env%eps_filter
4749 CALL dbcsr_create(tmp_n_vd, template=almo_scf_env%matrix_v_disc(ispin))
4751 template=almo_scf_env%matrix_vv_disc_blk(ispin), &
4752 matrix_type=dbcsr_type_no_symmetry)
4754 almo_scf_env%matrix_s(1), &
4756 0.0_dp, tmp_n_vd, filter_eps=eps_filter)
4758 almo_scf_env%matrix_vv_disc_blk(ispin))
4760 0.0_dp, tmp_vd_vd_blk, &
4761 retain_sparsity=.true.)
4764 template=almo_scf_env%matrix_vv_disc_blk(ispin), &
4765 matrix_type=dbcsr_type_no_symmetry)
4767 almo_scf_env%opt_k_t_dd(ispin), &
4769 threshold=eps_filter, &
4770 order=almo_scf_env%order_lanczos, &
4771 eps_lanczos=almo_scf_env%eps_lanczos, &
4772 max_iter_lanczos=almo_scf_env%max_iter_lanczos)
4776 almo_scf_env%matrix_ks_0deloc(ispin), &
4778 0.0_dp, tmp_n_vd, filter_eps=eps_filter)
4780 almo_scf_env%matrix_vv_disc_blk(ispin))
4782 0.0_dp, tmp_vd_vd_blk, &
4783 retain_sparsity=.true.)
4789 almo_scf_env%opt_k_t_dd(ispin), &
4790 0.0_dp, s_dd_sqrt, filter_eps=eps_filter)
4792 almo_scf_env%opt_k_t_dd(ispin), &
4794 0.0_dp, tmp_vd_vd_blk, filter_eps=eps_filter)
4798 template=almo_scf_env%matrix_vv_disc_blk(ispin))
4801 template=almo_scf_env%matrix_vv_disc_blk(ispin), &
4802 matrix_type=dbcsr_type_no_symmetry)
4810 almo_scf_env%opt_k_t_dd(ispin))
4814 0.0_dp, almo_scf_env%opt_k_t_dd(ispin), &
4815 filter_eps=eps_filter)
4821 template=almo_scf_env%matrix_k_blk_ones(ispin))
4823 almo_scf_env%matrix_k_blk_ones(ispin))
4825 template=almo_scf_env%matrix_k_blk_ones(ispin))
4828 0.0_dp, t1, filter_eps=eps_filter)
4833 template=almo_scf_env%matrix_sigma_vv_blk(ispin), &
4834 matrix_type=dbcsr_type_no_symmetry)
4836 almo_scf_env%matrix_sigma_vv_blk(ispin))
4838 almo_scf_env%matrix_x(ispin), &
4840 0.0_dp, tmp_vr_vr_blk, &
4841 retain_sparsity=.true.)
4845 template=almo_scf_env%matrix_sigma_vv_blk(ispin), &
4846 matrix_type=dbcsr_type_no_symmetry)
4848 almo_scf_env%opt_k_t_rr(ispin), &
4850 threshold=eps_filter, &
4851 order=almo_scf_env%order_lanczos, &
4852 eps_lanczos=almo_scf_env%eps_lanczos, &
4853 max_iter_lanczos=almo_scf_env%max_iter_lanczos)
4857 template=almo_scf_env%matrix_v(ispin))
4859 template=almo_scf_env%matrix_v(ispin))
4861 0.0_dp, tmp1_n_vr, filter_eps=eps_filter)
4863 almo_scf_env%matrix_ks_0deloc(ispin), &
4865 0.0_dp, tmp2_n_vr, filter_eps=eps_filter)
4867 0.0_dp, tmp_vr_vr_blk, &
4868 retain_sparsity=.true.)
4875 almo_scf_env%opt_k_t_rr(ispin), &
4876 0.0_dp, s_rr_sqrt, filter_eps=eps_filter)
4878 almo_scf_env%opt_k_t_rr(ispin), &
4880 0.0_dp, tmp_vr_vr_blk, filter_eps=eps_filter)
4884 template=almo_scf_env%matrix_sigma_vv_blk(ispin))
4887 template=almo_scf_env%matrix_sigma_vv_blk(ispin), &
4888 matrix_type=dbcsr_type_no_symmetry)
4896 almo_scf_env%opt_k_t_rr(ispin))
4900 0.0_dp, almo_scf_env%opt_k_t_rr(ispin), &
4901 filter_eps=eps_filter)
4908 0.0_dp, almo_scf_env%opt_k_denom(ispin), &
4909 filter_eps=eps_filter)
4914 CALL dbcsr_add(almo_scf_env%opt_k_denom(ispin), t1, &
4917 CALL dbcsr_scale(almo_scf_env%opt_k_denom(ispin), &
4920 CALL inverse_of_elements(almo_scf_env%opt_k_denom(ispin))
4924 CALL timestop(handle)
4926 END SUBROUTINE opt_k_create_preconditioner_blk
4940 SUBROUTINE opt_k_apply_preconditioner_blk(almo_scf_env, step, grad, ispin)
4945 INTEGER,
INTENT(IN) :: ispin
4947 CHARACTER(len=*),
PARAMETER :: routinen =
'opt_k_apply_preconditioner_blk'
4950 REAL(kind=
dp) :: eps_filter
4953 CALL timeset(routinen, handle)
4955 eps_filter = almo_scf_env%eps_filter
4957 CALL dbcsr_create(tmp_k, template=almo_scf_env%matrix_k_blk(ispin))
4961 grad, almo_scf_env%opt_k_t_rr(ispin), &
4962 0.0_dp, tmp_k, filter_eps=eps_filter)
4964 almo_scf_env%opt_k_t_dd(ispin), tmp_k, &
4965 0.0_dp, step, filter_eps=eps_filter)
4969 almo_scf_env%opt_k_denom(ispin), tmp_k)
4973 almo_scf_env%opt_k_t_dd(ispin), tmp_k, &
4974 0.0_dp, step, filter_eps=eps_filter)
4976 step, almo_scf_env%opt_k_t_rr(ispin), &
4977 0.0_dp, tmp_k, filter_eps=eps_filter)
4983 CALL timestop(handle)
4985 END SUBROUTINE opt_k_apply_preconditioner_blk
5024 SUBROUTINE compute_gradient(m_grad_out, m_ks, m_s, m_t, m_t0, &
5025 m_siginv, m_quench_t, m_FTsiginv, m_siginvTFTsiginv, m_ST, m_STsiginv0, &
5026 m_theta, domain_s_inv, domain_r_down, &
5027 cpu_of_domain, domain_map, assume_t0_q0x, optimize_theta, &
5028 normalize_orbitals, penalty_occ_vol, penalty_occ_local, &
5029 penalty_occ_vol_prefactor, envelope_amplitude, eps_filter, spin_factor, &
5030 special_case, m_sig_sqrti_ii, op_sm_set, weights, energy_coeff, &
5033 TYPE(
dbcsr_type),
INTENT(INOUT) :: m_grad_out, m_ks, m_s, m_t, m_t0, &
5034 m_siginv, m_quench_t, m_ftsiginv, &
5035 m_siginvtftsiginv, m_st, m_stsiginv0, &
5038 INTENT(IN) :: domain_s_inv, domain_r_down
5039 INTEGER,
DIMENSION(:),
INTENT(IN) :: cpu_of_domain
5041 LOGICAL,
INTENT(IN) :: assume_t0_q0x, optimize_theta, &
5042 normalize_orbitals, penalty_occ_vol
5043 LOGICAL,
INTENT(IN),
OPTIONAL :: penalty_occ_local
5044 REAL(kind=
dp),
INTENT(IN) :: penalty_occ_vol_prefactor, &
5045 envelope_amplitude, eps_filter, &
5047 INTEGER,
INTENT(IN) :: special_case
5048 TYPE(
dbcsr_type),
INTENT(IN),
OPTIONAL :: m_sig_sqrti_ii
5050 POINTER :: op_sm_set
5051 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN),
OPTIONAL :: weights
5052 REAL(kind=
dp),
INTENT(IN),
OPTIONAL :: energy_coeff, localiz_coeff
5054 CHARACTER(len=*),
PARAMETER :: routinen =
'compute_gradient'
5056 INTEGER :: dim0, handle, idim0, nao, reim
5057 LOGICAL :: my_penalty_local
5058 REAL(kind=
dp) :: coeff, energy_g_norm, my_energy_coeff, &
5060 penalty_occ_vol_g_norm
5061 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: tg_diagonal
5062 TYPE(
dbcsr_type) :: m_tmp_no_1, m_tmp_no_2, m_tmp_no_3, &
5063 m_tmp_oo_1, m_tmp_oo_2, temp1, temp2, &
5064 tempnocc1, tempoccocc1
5066 CALL timeset(routinen, handle)
5068 IF (normalize_orbitals .AND. (.NOT.
PRESENT(m_sig_sqrti_ii)))
THEN
5069 cpabort(
"Normalization matrix is required")
5072 my_penalty_local = .false.
5073 my_localiz_coeff = 1.0_dp
5074 my_energy_coeff = 0.0_dp
5075 IF (
PRESENT(localiz_coeff))
THEN
5076 my_localiz_coeff = localiz_coeff
5078 IF (
PRESENT(energy_coeff))
THEN
5079 my_energy_coeff = energy_coeff
5081 IF (
PRESENT(penalty_occ_local))
THEN
5082 my_penalty_local = penalty_occ_local
5091 template=m_quench_t, &
5092 matrix_type=dbcsr_type_no_symmetry)
5094 template=m_quench_t, &
5095 matrix_type=dbcsr_type_no_symmetry)
5097 template=m_quench_t, &
5098 matrix_type=dbcsr_type_no_symmetry)
5100 template=m_siginv, &
5101 matrix_type=dbcsr_type_no_symmetry)
5103 template=m_siginv, &
5104 matrix_type=dbcsr_type_no_symmetry)
5107 matrix_type=dbcsr_type_no_symmetry)
5109 template=m_siginv, &
5110 matrix_type=dbcsr_type_no_symmetry)
5113 matrix_type=dbcsr_type_no_symmetry)
5116 matrix_type=dbcsr_type_no_symmetry)
5119 CALL dbcsr_copy(m_tmp_no_2, m_ftsiginv, keep_sparsity=.true.)
5123 m_siginvtftsiginv, &
5124 1.0_dp, m_tmp_no_2, &
5125 retain_sparsity=.true.)
5129 IF (my_penalty_local)
THEN
5133 DO idim0 = 1,
SIZE(op_sm_set, 2)
5135 DO reim = 1,
SIZE(op_sm_set, 1)
5138 op_sm_set(reim, idim0)%matrix, &
5140 0.0_dp, tempnocc1, &
5141 filter_eps=eps_filter)
5147 0.0_dp, tempoccocc1, &
5148 filter_eps=eps_filter)
5151 ALLOCATE (tg_diagonal(dim0))
5155 DEALLOCATE (tg_diagonal)
5161 filter_eps=eps_filter)
5167 cpabort(
"Localization function is not implemented")
5169 coeff = -weights(idim0)
5171 cpabort(
"Localization function is not implemented")
5173 CALL dbcsr_add(temp2, temp1, 1.0_dp, coeff)
5176 CALL dbcsr_add(m_tmp_no_2, temp2, my_energy_coeff, my_localiz_coeff*4.0_dp)
5180 IF (penalty_occ_vol)
THEN
5183 penalty_occ_vol_prefactor, &
5186 0.0_dp, m_tmp_no_1, &
5187 retain_sparsity=.true.)
5191 CALL dbcsr_add(m_tmp_no_2, m_tmp_no_1, 1.0_dp, 1.0_dp)
5195 IF (normalize_orbitals)
THEN
5209 0.0_dp, m_tmp_no_1, &
5210 retain_sparsity=.true.)
5217 0.0_dp, m_tmp_oo_1, &
5218 retain_sparsity=.true.)
5221 ALLOCATE (tg_diagonal(dim0))
5225 DEALLOCATE (tg_diagonal)
5230 0.0_dp, m_tmp_oo_2, &
5231 filter_eps=eps_filter)
5235 1.0_dp, m_tmp_no_1, &
5236 retain_sparsity=.true.)
5245 IF (assume_t0_q0x)
THEN
5251 0.0_dp, m_tmp_oo_1, &
5252 filter_eps=eps_filter)
5256 1.0_dp, m_grad_out, &
5257 filter_eps=eps_filter)
5259 cpabort(
"Cannot project the zero-order space from itself")
5263 matrix_in=m_tmp_no_1, &
5264 matrix_out=m_grad_out, &
5265 operator2=domain_r_down(:), &
5266 operator1=domain_s_inv(:), &
5267 dpattern=m_quench_t, &
5269 node_of_domain=cpu_of_domain, &
5271 filter_eps=eps_filter, &
5273 use_trimmer=.false.)
5279 IF (optimize_theta)
THEN
5281 CALL dtanh_of_elements(m_tmp_no_2, alpha=1.0_dp/envelope_amplitude)
5308 CALL timestop(handle)
5310 END SUBROUTINE compute_gradient
5320 SUBROUTINE print_mathematica_matrix(matrix, filename)
5323 CHARACTER(len=*),
INTENT(IN) :: filename
5325 CHARACTER(len=*),
PARAMETER :: routinen =
'print_mathematica_matrix'
5327 CHARACTER(LEN=20) :: formatstr, scols
5328 INTEGER :: col, fiunit, handle, hori_offset, jj, &
5329 nblkcols_tot, nblkrows_tot, ncols, &
5330 ncores, nrows, row, unit_nr, &
5332 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: ao_block_sizes, mo_block_sizes
5333 INTEGER,
DIMENSION(:),
POINTER :: ao_blk_sizes, mo_blk_sizes
5335 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: h
5336 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: block_p
5341 CALL timeset(routinen, handle)
5345 IF (logger%para_env%is_source())
THEN
5354 IF (ncores > 1)
THEN
5355 cpabort(
"mathematica files: serial code only")
5358 CALL dbcsr_get_info(matrix, row_blk_size=ao_blk_sizes, col_blk_size=mo_blk_sizes, &
5359 nblkrows_total=nblkrows_tot, nblkcols_total=nblkcols_tot)
5360 cpassert(nblkrows_tot == nblkcols_tot)
5361 ALLOCATE (mo_block_sizes(nblkcols_tot), ao_block_sizes(nblkcols_tot))
5362 mo_block_sizes(:) = mo_blk_sizes(:)
5363 ao_block_sizes(:) = ao_blk_sizes(:)
5367 matrix_type=dbcsr_type_no_symmetry)
5370 ncols = sum(mo_block_sizes)
5371 nrows = sum(ao_block_sizes)
5372 ALLOCATE (h(nrows, ncols))
5376 DO col = 1, nblkcols_tot
5379 DO row = 1, nblkrows_tot
5384 h(vert_offset + 1:vert_offset + ao_block_sizes(row), &
5385 hori_offset + 1:hori_offset + mo_block_sizes(col)) &
5390 vert_offset = vert_offset + ao_block_sizes(row)
5394 hori_offset = hori_offset + mo_block_sizes(col)
5400 IF (unit_nr > 0)
THEN
5401 CALL open_file(filename, unit_number=fiunit, file_status=
'REPLACE')
5402 WRITE (scols,
"(I10)") ncols
5403 formatstr =
"("//trim(scols)//
"E27.17)"
5405 WRITE (fiunit, formatstr) h(jj, :)
5410 DEALLOCATE (mo_block_sizes)
5411 DEALLOCATE (ao_block_sizes)
5414 CALL timestop(handle)
5416 END SUBROUTINE print_mathematica_matrix
5438 SUBROUTINE compute_obj_nlmos(localization_obj_function_ispin, penalty_func_ispin, &
5439 penalty_vol_prefactor, overlap_determinant, m_sigma, nocc, m_B0, &
5440 m_theta_normalized, template_matrix_mo, weights, m_S0, just_started, &
5441 penalty_amplitude, eps_filter)
5443 REAL(kind=
dp),
INTENT(INOUT) :: localization_obj_function_ispin, penalty_func_ispin, &
5444 penalty_vol_prefactor, overlap_determinant
5446 INTEGER,
INTENT(IN) :: nocc
5447 TYPE(
dbcsr_type),
DIMENSION(:, :),
INTENT(IN) :: m_b0
5448 TYPE(
dbcsr_type),
INTENT(IN) :: m_theta_normalized, template_matrix_mo
5449 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: weights
5451 LOGICAL,
INTENT(IN) :: just_started
5452 REAL(kind=
dp),
INTENT(IN) :: penalty_amplitude, eps_filter
5454 CHARACTER(len=*),
PARAMETER :: routinen =
'compute_obj_nlmos'
5456 INTEGER :: handle, idim0, ielem, reim
5457 REAL(kind=
dp) :: det1, fval
5458 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: reim_diag, z2
5459 TYPE(
dbcsr_type) :: tempnocc1, tempoccocc1, tempoccocc2
5462 CALL timeset(routinen, handle)
5465 template=template_matrix_mo, &
5466 matrix_type=dbcsr_type_no_symmetry)
5468 template=m_theta_normalized, &
5469 matrix_type=dbcsr_type_no_symmetry)
5471 template=m_theta_normalized, &
5472 matrix_type=dbcsr_type_no_symmetry)
5474 localization_obj_function_ispin = 0.0_dp
5475 penalty_func_ispin = 0.0_dp
5477 ALLOCATE (reim_diag(nocc))
5481 DO idim0 = 1,
SIZE(m_b0, 2)
5485 DO reim = 1,
SIZE(m_b0, 1)
5488 m_b0(reim, idim0), &
5489 m_theta_normalized, &
5490 0.0_dp, tempoccocc1, &
5491 filter_eps=eps_filter)
5495 m_theta_normalized, &
5497 0.0_dp, tempoccocc2, &
5498 retain_sparsity=.true.)
5502 CALL group%sum(reim_diag)
5503 z2(:) = z2(:) + reim_diag(:)*reim_diag(:)
5510 fval = -weights(idim0)*log(abs(z2(ielem)))
5512 fval = weights(idim0) - weights(idim0)*abs(z2(ielem))
5514 fval = weights(idim0) - weights(idim0)*sqrt(abs(z2(ielem)))
5516 localization_obj_function_ispin = localization_obj_function_ispin + fval
5522 DEALLOCATE (reim_diag)
5526 m_theta_normalized, &
5527 0.0_dp, tempoccocc1, &
5528 filter_eps=eps_filter)
5531 m_theta_normalized, &
5534 filter_eps=eps_filter)
5539 overlap_determinant = det1
5541 IF (just_started .AND. penalty_amplitude < 0.0_dp)
THEN
5542 penalty_vol_prefactor = -(-penalty_amplitude)*localization_obj_function_ispin
5544 penalty_func_ispin = penalty_func_ispin + penalty_vol_prefactor*log(det1)
5550 CALL timestop(handle)
5552 END SUBROUTINE compute_obj_nlmos
5570 SUBROUTINE compute_gradient_nlmos(m_grad_out, m_B0, weights, &
5571 m_S0, m_theta_normalized, m_siginv, m_sig_sqrti_ii, &
5572 penalty_vol_prefactor, eps_filter, suggested_vol_penalty)
5574 TYPE(
dbcsr_type),
INTENT(INOUT) :: m_grad_out
5575 TYPE(
dbcsr_type),
DIMENSION(:, :),
INTENT(IN) :: m_b0
5576 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: weights
5577 TYPE(
dbcsr_type),
INTENT(IN) :: m_s0, m_theta_normalized, m_siginv, &
5579 REAL(kind=
dp),
INTENT(IN) :: penalty_vol_prefactor, eps_filter
5580 REAL(kind=
dp),
INTENT(INOUT) :: suggested_vol_penalty
5582 CHARACTER(len=*),
PARAMETER :: routinen =
'compute_gradient_nlmos'
5584 INTEGER :: dim0, handle, idim0, reim
5585 REAL(kind=
dp) :: norm_loc, norm_vol
5586 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: tg_diagonal, z2
5587 TYPE(
dbcsr_type) :: m_temp_oo_1, m_temp_oo_2, m_temp_oo_3, &
5590 CALL timeset(routinen, handle)
5593 template=m_theta_normalized, &
5594 matrix_type=dbcsr_type_no_symmetry)
5596 template=m_theta_normalized, &
5597 matrix_type=dbcsr_type_no_symmetry)
5599 template=m_theta_normalized, &
5600 matrix_type=dbcsr_type_no_symmetry)
5602 template=m_theta_normalized, &
5603 matrix_type=dbcsr_type_no_symmetry)
5606 ALLOCATE (tg_diagonal(dim0))
5611 DO idim0 = 1,
SIZE(m_b0, 2)
5615 DO reim = 1,
SIZE(m_b0, 1)
5618 m_b0(reim, idim0), &
5619 m_theta_normalized, &
5620 0.0_dp, m_temp_oo_3, &
5621 filter_eps=eps_filter)
5626 m_theta_normalized, &
5628 0.0_dp, m_temp_oo_4, &
5629 filter_eps=eps_filter)
5631 tg_diagonal(:) = 0.0_dp
5635 z2(:) = z2(:) + tg_diagonal(:)*tg_diagonal(:)
5640 1.0_dp, m_temp_oo_2, &
5641 filter_eps=eps_filter)
5649 z2(:) = -weights(idim0)/z2(:)
5651 z2(:) = -weights(idim0)
5653 z2(:) = -weights(idim0)/(2*sqrt(z2(:)))
5663 1.0_dp, m_temp_oo_1, &
5664 filter_eps=eps_filter)
5673 m_theta_normalized, &
5674 0.0_dp, m_temp_oo_2, &
5675 filter_eps=eps_filter)
5683 0.0_dp, m_temp_oo_3, &
5684 filter_eps=eps_filter)
5687 suggested_vol_penalty = norm_loc/norm_vol
5688 CALL dbcsr_add(m_temp_oo_1, m_temp_oo_3, &
5689 1.0_dp, 2.0_dp*penalty_vol_prefactor)
5697 0.0_dp, m_grad_out, &
5698 filter_eps=eps_filter)
5703 m_theta_normalized, &
5705 0.0_dp, m_temp_oo_3, &
5706 filter_eps=eps_filter)
5716 0.0_dp, m_temp_oo_1, &
5717 filter_eps=eps_filter)
5722 1.0_dp, m_grad_out, &
5723 filter_eps=eps_filter)
5725 DEALLOCATE (tg_diagonal)
5731 CALL timestop(handle)
5733 END SUBROUTINE compute_gradient_nlmos
5764 SUBROUTINE compute_xalmos_from_main_var(m_var_in, m_t_out, m_quench_t, &
5765 m_t0, m_oo_template, m_STsiginv0, m_s, m_sig_sqrti_ii_out, domain_r_down, &
5766 domain_s_inv, domain_map, cpu_of_domain, assume_t0_q0x, just_started, &
5767 optimize_theta, normalize_orbitals, envelope_amplitude, eps_filter, &
5768 special_case, nocc_of_domain, order_lanczos, eps_lanczos, max_iter_lanczos)
5771 TYPE(
dbcsr_type),
INTENT(INOUT) :: m_t_out, m_quench_t, m_t0, &
5772 m_oo_template, m_stsiginv0, m_s, &
5775 INTENT(IN) :: domain_r_down, domain_s_inv
5777 INTEGER,
DIMENSION(:),
INTENT(IN) :: cpu_of_domain
5778 LOGICAL,
INTENT(IN) :: assume_t0_q0x, just_started, &
5779 optimize_theta, normalize_orbitals
5780 REAL(kind=
dp),
INTENT(IN) :: envelope_amplitude, eps_filter
5781 INTEGER,
INTENT(IN) :: special_case
5782 INTEGER,
DIMENSION(:),
INTENT(IN) :: nocc_of_domain
5783 INTEGER,
INTENT(IN) :: order_lanczos
5784 REAL(kind=
dp),
INTENT(IN) :: eps_lanczos
5785 INTEGER,
INTENT(IN) :: max_iter_lanczos
5787 CHARACTER(len=*),
PARAMETER :: routinen =
'compute_xalmos_from_main_var'
5789 INTEGER :: handle, unit_nr
5790 REAL(kind=
dp) :: t_norm
5794 CALL timeset(routinen, handle)
5798 IF (logger%para_env%is_source())
THEN
5805 template=m_quench_t, &
5806 matrix_type=dbcsr_type_no_symmetry)
5808 template=m_oo_template, &
5809 matrix_type=dbcsr_type_no_symmetry)
5812 IF (optimize_theta)
THEN
5816 IF (unit_nr > 0)
THEN
5817 WRITE (unit_nr, *)
"Maximum norm of the initial guess: ", t_norm
5818 WRITE (unit_nr, *)
"Maximum allowed amplitude: ", &
5821 IF (t_norm > envelope_amplitude .AND. just_started)
THEN
5822 cpabort(
"Max norm of the initial guess is too large")
5825 CALL tanh_of_elements(m_tmp_no_1, alpha=1.0_dp/envelope_amplitude)
5832 IF (assume_t0_q0x)
THEN
5837 0.0_dp, m_tmp_oo_1, &
5838 filter_eps=eps_filter)
5843 filter_eps=eps_filter)
5845 cpabort(
"cannot use projector with block-daigonal ALMOs")
5849 matrix_in=m_t_out, &
5850 matrix_out=m_tmp_no_1, &
5851 operator1=domain_r_down, &
5852 operator2=domain_s_inv, &
5853 dpattern=m_quench_t, &
5855 node_of_domain=cpu_of_domain, &
5857 filter_eps=eps_filter, &
5858 use_trimmer=.false.)
5863 m_t0, 1.0_dp, 1.0_dp)
5866 IF (normalize_orbitals)
THEN
5869 overlap=m_tmp_oo_1, &
5871 retain_locality=.true., &
5872 only_normalize=.true., &
5873 nocc_of_domain=nocc_of_domain(:), &
5874 eps_filter=eps_filter, &
5875 order_lanczos=order_lanczos, &
5876 eps_lanczos=eps_lanczos, &
5877 max_iter_lanczos=max_iter_lanczos, &
5878 overlap_sqrti=m_sig_sqrti_ii_out)
5886 CALL timestop(handle)
5888 END SUBROUTINE compute_xalmos_from_main_var
5926 SUBROUTINE compute_preconditioner(domain_prec_out, m_prec_out, m_ks, m_s, &
5927 m_siginv, m_quench_t, m_FTsiginv, m_siginvTFTsiginv, m_ST, &
5928 m_STsiginv_out, m_s_vv_out, m_f_vv_out, para_env, &
5929 blacs_env, nocc_of_domain, domain_s_inv, domain_s_inv_half, domain_s_half, &
5930 domain_r_down, cpu_of_domain, &
5931 domain_map, assume_t0_q0x, penalty_occ_vol, penalty_occ_vol_prefactor, &
5932 eps_filter, neg_thr, spin_factor, special_case, bad_modes_projector_down_out, &
5936 INTENT(INOUT) :: domain_prec_out
5937 TYPE(
dbcsr_type),
INTENT(INOUT) :: m_prec_out, m_ks, m_s
5938 TYPE(
dbcsr_type),
INTENT(IN) :: m_siginv, m_quench_t, m_ftsiginv, &
5939 m_siginvtftsiginv, m_st
5940 TYPE(
dbcsr_type),
INTENT(INOUT),
OPTIONAL :: m_stsiginv_out, m_s_vv_out, m_f_vv_out
5943 INTEGER,
DIMENSION(:),
INTENT(IN) :: nocc_of_domain
5945 INTENT(IN) :: domain_s_inv
5947 INTENT(IN),
OPTIONAL :: domain_s_inv_half, domain_s_half
5949 INTENT(IN) :: domain_r_down
5950 INTEGER,
DIMENSION(:),
INTENT(IN) :: cpu_of_domain
5952 LOGICAL,
INTENT(IN) :: assume_t0_q0x, penalty_occ_vol
5953 REAL(kind=
dp),
INTENT(IN) :: penalty_occ_vol_prefactor, eps_filter, &
5954 neg_thr, spin_factor
5955 INTEGER,
INTENT(IN) :: special_case
5957 INTENT(INOUT),
OPTIONAL :: bad_modes_projector_down_out
5958 LOGICAL,
INTENT(IN) :: skip_inversion
5960 CHARACTER(len=*),
PARAMETER :: routinen =
'compute_preconditioner'
5962 INTEGER :: handle, ndim, precond_domain_projector
5963 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: nn_diagonal
5966 CALL timeset(routinen, handle)
5970 matrix_type=dbcsr_type_no_symmetry)
5972 template=m_quench_t, &
5973 matrix_type=dbcsr_type_no_symmetry)
5985 0.0_dp, m_tmp_no_3, &
5986 filter_eps=eps_filter)
5989 IF (
PRESENT(m_stsiginv_out))
THEN
5998 1.0_dp, m_tmp_nn_1, &
5999 filter_eps=eps_filter)
6002 IF (
PRESENT(m_s_vv_out))
THEN
6013 1.0_dp, m_prec_out, &
6014 filter_eps=eps_filter)
6018 1.0_dp, m_prec_out, &
6019 filter_eps=eps_filter)
6022 m_siginvtftsiginv, &
6023 0.0_dp, m_tmp_no_3, &
6024 filter_eps=eps_filter)
6028 1.0_dp, m_prec_out, &
6029 filter_eps=eps_filter)
6031 IF (
PRESENT(m_f_vv_out))
THEN
6036 CALL dbcsr_add(m_prec_out, m_tmp_nn_1, &
6042 IF (penalty_occ_vol)
THEN
6043 CALL dbcsr_add(m_prec_out, m_tmp_nn_1, &
6044 1.0_dp, penalty_occ_vol_prefactor)
6052 IF (skip_inversion)
THEN
6056 ALLOCATE (nn_diagonal(ndim))
6061 DEALLOCATE (nn_diagonal)
6063 CALL dbcsr_copy(m_prec_out, m_tmp_nn_1, keep_sparsity=.true.)
6068 matrix_in=m_tmp_nn_1, &
6069 matrix_out=m_prec_out, &
6070 nocc=nocc_of_domain(:) &
6077 IF (skip_inversion)
THEN
6083 para_env=para_env, &
6084 blacs_env=blacs_env)
6086 para_env=para_env, &
6087 blacs_env=blacs_env, &
6088 uplo_to_full=.true.)
6096 IF (assume_t0_q0x)
THEN
6097 precond_domain_projector = -1
6099 precond_domain_projector = 0
6104 IF (
PRESENT(bad_modes_projector_down_out))
THEN
6106 matrix_main=m_tmp_nn_1, &
6107 subm_s_inv=domain_s_inv(:), &
6108 subm_s_inv_half=domain_s_inv_half(:), &
6109 subm_s_half=domain_s_half(:), &
6110 subm_r_down=domain_r_down(:), &
6111 matrix_trimmer=m_quench_t, &
6112 dpattern=m_quench_t, &
6114 node_of_domain=cpu_of_domain, &
6116 use_trimmer=.false., &
6117 bad_modes_projector_down=bad_modes_projector_down_out(:), &
6118 eps_zero_eigenvalues=neg_thr, &
6119 my_action=precond_domain_projector, &
6120 skip_inversion=skip_inversion &
6124 matrix_main=m_tmp_nn_1, &
6125 subm_s_inv=domain_s_inv(:), &
6126 subm_r_down=domain_r_down(:), &
6127 matrix_trimmer=m_quench_t, &
6128 dpattern=m_quench_t, &
6130 node_of_domain=cpu_of_domain, &
6132 use_trimmer=.false., &
6134 my_action=precond_domain_projector, &
6135 skip_inversion=skip_inversion &
6144 CALL timestop(handle)
6146 END SUBROUTINE compute_preconditioner
6164 SUBROUTINE compute_cg_beta(beta, numer, denom, reset_conjugator, conjugator, &
6165 grad, prev_grad, step, prev_step, prev_minus_prec_grad)
6167 REAL(kind=
dp),
INTENT(INOUT) :: beta
6168 REAL(kind=
dp),
INTENT(INOUT),
OPTIONAL :: numer, denom
6169 LOGICAL,
INTENT(INOUT) :: reset_conjugator
6170 INTEGER,
INTENT(IN) :: conjugator
6171 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(INOUT) :: grad, prev_grad, step, prev_step
6172 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(INOUT), &
6173 OPTIONAL :: prev_minus_prec_grad
6175 CHARACTER(len=*),
PARAMETER :: routinen =
'compute_cg_beta'
6177 INTEGER :: handle, i, nsize, unit_nr
6178 REAL(kind=
dp) :: den, kappa, my_denom, my_numer, &
6179 my_numer2, my_numer3, num, num2, num3, &
6184 CALL timeset(routinen, handle)
6188 IF (logger%para_env%is_source())
THEN
6194 IF (.NOT.
PRESENT(prev_minus_prec_grad))
THEN
6198 cpabort(
"conjugator needs more input")
6203 IF (
PRESENT(numer) .OR.
PRESENT(denom))
THEN
6207 cpabort(
"cannot return numer/denom")
6222 matrix_type=dbcsr_type_no_symmetry)
6224 SELECT CASE (conjugator)
6227 CALL dbcsr_add(m_tmp_no_1, prev_grad(i), &
6229 CALL dbcsr_dot(m_tmp_no_1, step(i), num)
6230 CALL dbcsr_dot(m_tmp_no_1, prev_step(i), den)
6233 CALL dbcsr_dot(prev_grad(i), prev_minus_prec_grad(i), den)
6235 CALL dbcsr_dot(prev_grad(i), prev_minus_prec_grad(i), den)
6237 CALL dbcsr_add(m_tmp_no_1, prev_grad(i), 1.0_dp, -1.0_dp)
6238 CALL dbcsr_dot(m_tmp_no_1, step(i), num)
6241 CALL dbcsr_dot(prev_grad(i), prev_step(i), den)
6243 CALL dbcsr_dot(prev_grad(i), prev_step(i), den)
6245 CALL dbcsr_add(m_tmp_no_1, prev_grad(i), 1.0_dp, -1.0_dp)
6246 CALL dbcsr_dot(m_tmp_no_1, step(i), num)
6250 CALL dbcsr_add(m_tmp_no_1, prev_grad(i), 1.0_dp, -1.0_dp)
6251 CALL dbcsr_dot(m_tmp_no_1, prev_step(i), den)
6254 CALL dbcsr_add(m_tmp_no_1, prev_grad(i), 1.0_dp, -1.0_dp)
6255 CALL dbcsr_dot(m_tmp_no_1, prev_step(i), den)
6256 CALL dbcsr_dot(m_tmp_no_1, prev_minus_prec_grad(i), num)
6257 CALL dbcsr_dot(m_tmp_no_1, step(i), num2)
6258 CALL dbcsr_dot(prev_step(i), grad(i), num3)
6259 my_numer2 = my_numer2 + num2
6260 my_numer3 = my_numer3 + num3
6265 cpabort(
"illegal conjugator")
6267 my_numer = my_numer + num
6268 my_denom = my_denom + den
6276 SELECT CASE (conjugator)
6278 beta = -1.0_dp*my_numer/my_denom
6280 beta = my_numer/my_denom
6282 kappa = -2.0_dp*my_numer/my_denom
6283 tau = -1.0_dp*my_numer2/my_denom
6284 beta = tau - kappa*my_numer3/my_denom
6288 cpabort(
"illegal conjugator")
6293 IF (beta < 0.0_dp)
THEN
6294 IF (unit_nr > 0)
THEN
6295 WRITE (unit_nr, *)
" Resetting conjugator because beta is negative: ", beta
6297 reset_conjugator = .true.
6300 IF (
PRESENT(numer))
THEN
6303 IF (
PRESENT(denom))
THEN
6307 CALL timestop(handle)
6309 END SUBROUTINE compute_cg_beta
6343 SUBROUTINE newton_grad_to_step(optimizer, m_grad, m_delta, m_s, m_ks, &
6344 m_siginv, m_quench_t, m_FTsiginv, m_siginvTFTsiginv, m_ST, m_t, &
6345 m_sig_sqrti_ii, domain_s_inv, domain_r_down, domain_map, cpu_of_domain, &
6346 nocc_of_domain, para_env, blacs_env, eps_filter, optimize_theta, &
6347 penalty_occ_vol, normalize_orbitals, penalty_occ_vol_prefactor, &
6348 penalty_occ_vol_pf2, special_case)
6351 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(IN) :: m_grad
6352 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(INOUT) :: m_delta, m_s, m_ks, m_siginv, m_quench_t
6353 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(IN) :: m_ftsiginv, m_siginvtftsiginv, m_st, &
6356 INTENT(IN) :: domain_s_inv, domain_r_down
6358 INTEGER,
DIMENSION(:),
INTENT(IN) :: cpu_of_domain
6359 INTEGER,
DIMENSION(:, :),
INTENT(IN) :: nocc_of_domain
6362 REAL(kind=
dp),
INTENT(IN) :: eps_filter
6363 LOGICAL,
INTENT(IN) :: optimize_theta, penalty_occ_vol, &
6365 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: penalty_occ_vol_prefactor, &
6367 INTEGER,
INTENT(IN) :: special_case
6369 CHARACTER(len=*),
PARAMETER :: routinen =
'newton_grad_to_step'
6371 CHARACTER(LEN=20) :: iter_type
6372 INTEGER :: handle, ispin, iteration, max_iter, &
6373 ndomains, nspins, outer_iteration, &
6374 outer_max_iter, unit_nr
6375 LOGICAL :: converged, do_exact_inversion, outer_prepare_to_exit, prepare_to_exit, &
6376 reset_conjugator, use_preconditioner
6377 REAL(kind=
dp) :: alpha, beta, denom, denom_ispin, &
6378 eps_error_target, numer, numer_ispin, &
6379 residue_norm, spin_factor, t1, t2
6380 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: residue_max_norm
6383 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:) :: m_f_vo, m_f_vv, m_hstep, m_prec, &
6384 m_residue, m_residue_prev, m_s_vv, &
6385 m_step, m_stsiginv, m_zet, m_zet_prev
6387 DIMENSION(:, :) :: domain_prec
6389 CALL timeset(routinen, handle)
6393 IF (logger%para_env%is_source())
THEN
6400 IF (optimize_theta)
THEN
6401 cpabort(
"theta is NYI")
6406 outer_max_iter = optimizer%max_iter_outer_loop
6407 max_iter = optimizer%max_iter
6408 eps_error_target = optimizer%eps_error
6412 ndomains =
SIZE(domain_s_inv, 1)
6414 IF (nspins == 1)
THEN
6415 spin_factor = 2.0_dp
6417 spin_factor = 1.0_dp
6420 ALLOCATE (domain_prec(ndomains, nspins))
6424 ALLOCATE (m_residue(nspins))
6425 ALLOCATE (m_residue_prev(nspins))
6426 ALLOCATE (m_step(nspins))
6427 ALLOCATE (m_zet(nspins))
6428 ALLOCATE (m_zet_prev(nspins))
6429 ALLOCATE (m_hstep(nspins))
6430 ALLOCATE (m_prec(nspins))
6431 ALLOCATE (m_s_vv(nspins))
6432 ALLOCATE (m_f_vv(nspins))
6433 ALLOCATE (m_f_vo(nspins))
6434 ALLOCATE (m_stsiginv(nspins))
6436 ALLOCATE (residue_max_norm(nspins))
6439 DO ispin = 1, nspins
6443 template=m_quench_t(ispin), &
6444 matrix_type=dbcsr_type_no_symmetry)
6446 template=m_quench_t(ispin), &
6447 matrix_type=dbcsr_type_no_symmetry)
6449 template=m_quench_t(ispin), &
6450 matrix_type=dbcsr_type_no_symmetry)
6452 template=m_quench_t(ispin), &
6453 matrix_type=dbcsr_type_no_symmetry)
6455 template=m_quench_t(ispin), &
6456 matrix_type=dbcsr_type_no_symmetry)
6458 template=m_quench_t(ispin), &
6459 matrix_type=dbcsr_type_no_symmetry)
6461 template=m_quench_t(ispin), &
6462 matrix_type=dbcsr_type_no_symmetry)
6464 template=m_quench_t(ispin), &
6465 matrix_type=dbcsr_type_no_symmetry)
6467 template=m_ks(ispin), &
6468 matrix_type=dbcsr_type_no_symmetry)
6471 matrix_type=dbcsr_type_no_symmetry)
6473 template=m_ks(ispin), &
6474 matrix_type=dbcsr_type_no_symmetry)
6478 CALL dbcsr_copy(m_f_vo(ispin), m_ftsiginv(ispin))
6481 m_siginvtftsiginv(ispin), &
6482 1.0_dp, m_f_vo(ispin), &
6483 filter_eps=eps_filter)
6488 CALL compute_preconditioner( &
6489 domain_prec_out=domain_prec(:, ispin), &
6490 m_prec_out=m_prec(ispin), &
6493 m_siginv=m_siginv(ispin), &
6494 m_quench_t=m_quench_t(ispin), &
6495 m_ftsiginv=m_ftsiginv(ispin), &
6496 m_siginvtftsiginv=m_siginvtftsiginv(ispin), &
6498 m_stsiginv_out=m_stsiginv(ispin), &
6499 m_s_vv_out=m_s_vv(ispin), &
6500 m_f_vv_out=m_f_vv(ispin), &
6501 para_env=para_env, &
6502 blacs_env=blacs_env, &
6503 nocc_of_domain=nocc_of_domain(:, ispin), &
6504 domain_s_inv=domain_s_inv(:, ispin), &
6505 domain_r_down=domain_r_down(:, ispin), &
6506 cpu_of_domain=cpu_of_domain(:), &
6507 domain_map=domain_map(ispin), &
6508 assume_t0_q0x=.false., &
6509 penalty_occ_vol=penalty_occ_vol, &
6510 penalty_occ_vol_prefactor=penalty_occ_vol_prefactor(ispin), &
6511 eps_filter=eps_filter, &
6513 spin_factor=spin_factor, &
6514 special_case=special_case, &
6515 skip_inversion=.false. &
6519 CALL dbcsr_copy(m_delta(ispin), m_quench_t(ispin))
6522 CALL dbcsr_copy(m_residue(ispin), m_grad(ispin))
6525 do_exact_inversion = .false.
6526 IF (do_exact_inversion)
THEN
6530 CALL dbcsr_copy(m_step(ispin), m_grad(ispin))
6534 CALL hessian_diag_apply( &
6535 matrix_grad=m_step(ispin), &
6536 matrix_step=m_zet(ispin), &
6537 matrix_s_ao=m_s_vv(ispin), &
6538 matrix_f_ao=m_f_vv(ispin), &
6541 matrix_s_mo=m_siginv(ispin), &
6542 matrix_f_mo=m_siginvtftsiginv(ispin), &
6543 matrix_s_vo=m_stsiginv(ispin), &
6544 matrix_f_vo=m_f_vo(ispin), &
6545 quench_t=m_quench_t(ispin), &
6546 spin_factor=spin_factor, &
6547 eps_zero=eps_filter*10.0_dp, &
6548 penalty_occ_vol=penalty_occ_vol, &
6549 penalty_occ_vol_prefactor=penalty_occ_vol_prefactor(ispin), &
6550 penalty_occ_vol_pf2=penalty_occ_vol_pf2(ispin), &
6552 para_env=para_env, &
6553 blacs_env=blacs_env)
6557 IF (use_preconditioner)
THEN
6565 0.0_dp, m_zet(ispin), &
6566 filter_eps=eps_filter)
6571 matrix_in=m_residue(ispin), &
6572 matrix_out=m_zet(ispin), &
6573 operator1=domain_prec(:, ispin), &
6574 dpattern=m_quench_t(ispin), &
6575 map=domain_map(ispin), &
6576 node_of_domain=cpu_of_domain(:), &
6578 filter_eps=eps_filter)
6584 CALL dbcsr_copy(m_zet(ispin), m_residue(ispin))
6595 outer_prepare_to_exit = .false.
6597 residue_norm = 0.0_dp
6602 prepare_to_exit = .false.
6610 CALL apply_hessian( &
6615 m_siginv=m_siginv, &
6616 m_quench_t=m_quench_t, &
6617 m_ftsiginv=m_ftsiginv, &
6618 m_siginvtftsiginv=m_siginvtftsiginv, &
6620 m_stsiginv=m_stsiginv, &
6627 m_sig_sqrti_ii=m_sig_sqrti_ii, &
6628 penalty_occ_vol=penalty_occ_vol, &
6629 normalize_orbitals=normalize_orbitals, &
6630 penalty_occ_vol_prefactor=penalty_occ_vol_prefactor, &
6631 eps_filter=eps_filter, &
6632 path_num=hessian_path_reuse)
6637 DO ispin = 1, nspins
6639 CALL dbcsr_dot(m_residue(ispin), m_zet(ispin), numer_ispin)
6640 CALL dbcsr_dot(m_step(ispin), m_hstep(ispin), denom_ispin)
6642 numer = numer + numer_ispin
6643 denom = denom + denom_ispin
6649 DO ispin = 1, nspins
6652 CALL dbcsr_add(m_delta(ispin), m_step(ispin), 1.0_dp, alpha)
6653 CALL dbcsr_copy(m_residue_prev(ispin), m_residue(ispin))
6654 CALL dbcsr_add(m_residue(ispin), m_hstep(ispin), &
6655 1.0_dp, -1.0_dp*alpha)
6656 residue_max_norm(ispin) =
dbcsr_maxabs(m_residue(ispin))
6661 residue_norm = maxval(residue_max_norm)
6662 converged = (residue_norm < eps_error_target)
6663 IF (converged .OR. (iteration >= max_iter))
THEN
6664 prepare_to_exit = .true.
6667 IF (.NOT. prepare_to_exit)
THEN
6669 DO ispin = 1, nspins
6672 CALL dbcsr_copy(m_zet_prev(ispin), m_zet(ispin))
6675 IF (use_preconditioner)
THEN
6683 0.0_dp, m_zet(ispin), &
6684 filter_eps=eps_filter)
6689 matrix_in=m_residue(ispin), &
6690 matrix_out=m_zet(ispin), &
6691 operator1=domain_prec(:, ispin), &
6692 dpattern=m_quench_t(ispin), &
6693 map=domain_map(ispin), &
6694 node_of_domain=cpu_of_domain(:), &
6696 filter_eps=eps_filter)
6702 CALL dbcsr_copy(m_zet(ispin), m_residue(ispin))
6709 CALL compute_cg_beta( &
6711 reset_conjugator=reset_conjugator, &
6714 prev_grad=m_residue_prev, &
6716 prev_step=m_zet_prev)
6718 DO ispin = 1, nspins
6721 CALL dbcsr_add(m_step(ispin), m_zet(ispin), beta, 1.0_dp)
6728 IF (unit_nr > 0)
THEN
6729 iter_type = trim(
"NR STEP")
6730 WRITE (unit_nr,
'(T6,A9,I6,F14.5,F14.5,F15.10,F9.2)') &
6731 iter_type, iteration, &
6732 alpha, beta, residue_norm, &
6737 iteration = iteration + 1
6738 IF (prepare_to_exit)
EXIT
6742 IF (converged .OR. (outer_iteration >= outer_max_iter))
THEN
6743 outer_prepare_to_exit = .true.
6746 outer_iteration = outer_iteration + 1
6747 IF (outer_prepare_to_exit)
EXIT
6751 DO ispin = 1, nspins
6755 template=m_siginv(ispin), &
6756 matrix_type=dbcsr_type_no_symmetry)
6758 template=m_siginv(ispin), &
6759 matrix_type=dbcsr_type_no_symmetry)
6763 0.0_dp, m_tmp_oo_1, &
6764 filter_eps=eps_filter)
6768 0.0_dp, m_tmp_oo_2, &
6769 filter_eps=eps_filter)
6770 CALL dbcsr_copy(m_zet(ispin), m_quench_t(ispin))
6774 0.0_dp, m_zet(ispin), &
6775 retain_sparsity=.true.)
6777 WRITE (unit_nr,
"(A50,2F20.10)")
"Occupied-space projection of the step", alpha
6778 CALL dbcsr_add(m_zet(ispin), m_delta(ispin), -1.0_dp, 1.0_dp)
6780 WRITE (unit_nr,
"(A50,2F20.10)")
"Virtual-space projection of the step", alpha
6782 WRITE (unit_nr,
"(A50,2F20.10)")
"Full step", alpha
6789 DO ispin = 1, nspins
6803 DEALLOCATE (domain_prec)
6804 DEALLOCATE (m_residue)
6805 DEALLOCATE (m_residue_prev)
6808 DEALLOCATE (m_zet_prev)
6810 DEALLOCATE (m_hstep)
6814 DEALLOCATE (m_stsiginv)
6815 DEALLOCATE (residue_max_norm)
6817 IF (.NOT. converged)
THEN
6818 cpabort(
"Optimization not converged!")
6823 CALL timestop(handle)
6825 END SUBROUTINE newton_grad_to_step
6853 SUBROUTINE apply_hessian(m_x_in, m_x_out, m_ks, m_s, m_siginv, &
6854 m_quench_t, m_FTsiginv, m_siginvTFTsiginv, m_ST, m_STsiginv, m_s_vv, &
6855 m_ks_vv, m_g_full, m_t, m_sig_sqrti_ii, penalty_occ_vol, &
6856 normalize_orbitals, penalty_occ_vol_prefactor, eps_filter, path_num)
6858 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(INOUT) :: m_x_in, m_x_out, m_ks, m_s
6859 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(IN) :: m_siginv, m_quench_t, m_ftsiginv, &
6860 m_siginvtftsiginv, m_st, m_stsiginv
6861 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(INOUT) :: m_s_vv, m_ks_vv, m_g_full
6862 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(IN) :: m_t, m_sig_sqrti_ii
6863 LOGICAL,
INTENT(IN) :: penalty_occ_vol, normalize_orbitals
6864 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: penalty_occ_vol_prefactor
6865 REAL(kind=
dp),
INTENT(IN) :: eps_filter
6866 INTEGER,
INTENT(IN) :: path_num
6868 CHARACTER(len=*),
PARAMETER :: routinen =
'apply_hessian'
6870 INTEGER :: dim0, handle, ispin, nspins
6871 REAL(kind=
dp) :: penalty_prefactor_local, spin_factor
6872 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: tg_diagonal
6873 TYPE(
dbcsr_type) :: m_tmp_no_1, m_tmp_no_2, m_tmp_oo_1, &
6876 CALL timeset(routinen, handle)
6879 IF (penalty_occ_vol) penalty_prefactor_local = 1._dp
6880 cpassert(
SIZE(m_stsiginv) >= 0)
6881 cpassert(
SIZE(m_siginvtftsiginv) >= 0)
6882 cpassert(
SIZE(m_s) >= 0)
6883 cpassert(
SIZE(m_g_full) >= 0)
6884 cpassert(
SIZE(m_ftsiginv) >= 0)
6885 mark_used(m_siginvtftsiginv)
6886 mark_used(m_stsiginv)
6887 mark_used(m_ftsiginv)
6893 IF (nspins == 1)
THEN
6894 spin_factor = 2.0_dp
6896 spin_factor = 1.0_dp
6899 DO ispin = 1, nspins
6901 penalty_prefactor_local = penalty_occ_vol_prefactor(ispin)/(2.0_dp*spin_factor)
6904 template=m_siginv(ispin), &
6905 matrix_type=dbcsr_type_no_symmetry)
6907 template=m_quench_t(ispin), &
6908 matrix_type=dbcsr_type_no_symmetry)
6910 template=m_quench_t(ispin), &
6911 matrix_type=dbcsr_type_no_symmetry)
6913 template=m_quench_t(ispin), &
6914 matrix_type=dbcsr_type_no_symmetry)
6917 IF (normalize_orbitals)
THEN
6922 CALL dbcsr_copy(m_tmp_oo_1, m_sig_sqrti_ii(ispin))
6926 0.0_dp, m_tmp_oo_1, &
6927 retain_sparsity=.true.)
6929 ALLOCATE (tg_diagonal(dim0))
6933 DEALLOCATE (tg_diagonal)
6939 1.0_dp, m_tmp_no_1, &
6940 filter_eps=eps_filter)
6943 m_sig_sqrti_ii(ispin), &
6944 0.0_dp, m_tmp_x_in, &
6945 filter_eps=eps_filter)
6953 IF (path_num == hessian_path_reuse)
THEN
6958 CALL dbcsr_copy(m_x_out(ispin), m_quench_t(ispin))
6962 0.0_dp, m_x_out(ispin), &
6963 retain_sparsity=.true.)
6965 CALL dbcsr_copy(m_tmp_no_2, m_quench_t(ispin))
6969 0.0_dp, m_tmp_no_2, &
6970 retain_sparsity=.true.)
6971 CALL dbcsr_add(m_x_out(ispin), m_tmp_no_2, &
6972 1.0_dp, -4.0_dp*penalty_prefactor_local + 1.0_dp)
6974 ELSE IF (path_num == hessian_path_assemble)
THEN
6979 cpabort(
"path is NYI")
6982 cpabort(
"illegal path")
6986 IF (normalize_orbitals)
THEN
6991 CALL dbcsr_copy(m_tmp_oo_1, m_sig_sqrti_ii(ispin))
6995 0.0_dp, m_tmp_oo_1, &
6996 retain_sparsity=.true.)
6998 ALLOCATE (tg_diagonal(dim0))
7002 DEALLOCATE (tg_diagonal)
7007 1.0_dp, m_x_out(ispin), &
7008 retain_sparsity=.true.)
7012 m_sig_sqrti_ii(ispin), &
7013 0.0_dp, m_x_out(ispin), &
7014 retain_sparsity=.true.)
7032 CALL timestop(handle)
7034 END SUBROUTINE apply_hessian
7059 SUBROUTINE hessian_diag_apply(matrix_grad, matrix_step, matrix_S_ao, &
7060 matrix_F_ao, matrix_S_mo, matrix_F_mo, matrix_S_vo, matrix_F_vo, quench_t, &
7061 penalty_occ_vol, penalty_occ_vol_prefactor, penalty_occ_vol_pf2, &
7062 spin_factor, eps_zero, m_s, para_env, blacs_env)
7064 TYPE(
dbcsr_type),
INTENT(INOUT) :: matrix_grad, matrix_step, matrix_s_ao, &
7065 matrix_f_ao, matrix_s_mo
7067 TYPE(
dbcsr_type),
INTENT(INOUT) :: matrix_s_vo, matrix_f_vo, quench_t
7068 LOGICAL,
INTENT(IN) :: penalty_occ_vol
7069 REAL(kind=
dp),
INTENT(IN) :: penalty_occ_vol_prefactor, &
7070 penalty_occ_vol_pf2, spin_factor, &
7076 CHARACTER(len=*),
PARAMETER :: routinen =
'hessian_diag_apply'
7078 INTEGER :: ao_hori_offset, ao_vert_offset, block_col, block_row, col, h_size, handle, ii, &
7079 info, jj, lev1_hori_offset, lev1_vert_offset, lev2_hori_offset, lev2_vert_offset, lwork, &
7080 nblkcols_tot, nblkrows_tot, ncores, orb_i, orb_j, row, unit_nr, zero_neg_eiv
7081 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: ao_block_sizes, ao_domain_sizes, &
7083 INTEGER,
DIMENSION(:),
POINTER :: ao_blk_sizes, mo_blk_sizes
7084 LOGICAL :: found, found_col, found_row
7085 REAL(kind=
dp) :: penalty_prefactor_local, test_error
7086 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: eigenvalues, grad_vec, step_vec, tmp, &
7088 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: f_ao_block, f_mo_block, h, hinv, &
7089 new_block, s_ao_block, s_mo_block, &
7091 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: block_p
7094 TYPE(
dbcsr_type) :: matrix_f_ao_sym, matrix_f_mo_sym, &
7095 matrix_s_ao_sym, matrix_s_mo_sym
7097 CALL timeset(routinen, handle)
7101 IF (logger%para_env%is_source())
THEN
7108 cpassert(
ASSOCIATED(blacs_env))
7109 cpassert(
ASSOCIATED(para_env))
7110 mark_used(blacs_env)
7120 IF (ncores > 1)
THEN
7121 cpabort(
"serial code only")
7124 CALL dbcsr_get_info(quench_t, row_blk_size=ao_blk_sizes, col_blk_size=mo_blk_sizes, &
7125 nblkrows_total=nblkrows_tot, nblkcols_total=nblkcols_tot)
7126 cpassert(nblkrows_tot == nblkcols_tot)
7127 ALLOCATE (mo_block_sizes(nblkcols_tot), ao_block_sizes(nblkcols_tot))
7128 ALLOCATE (ao_domain_sizes(nblkcols_tot))
7129 mo_block_sizes(:) = mo_blk_sizes(:)
7130 ao_block_sizes(:) = ao_blk_sizes(:)
7131 ao_domain_sizes(:) = 0
7134 template=matrix_s_ao, &
7135 matrix_type=dbcsr_type_no_symmetry)
7137 CALL dbcsr_scale(matrix_s_ao_sym, 2.0_dp*spin_factor)
7140 template=matrix_f_ao, &
7141 matrix_type=dbcsr_type_no_symmetry)
7143 CALL dbcsr_scale(matrix_f_ao_sym, 2.0_dp*spin_factor)
7146 template=matrix_s_mo, &
7147 matrix_type=dbcsr_type_no_symmetry)
7151 template=matrix_f_mo, &
7152 matrix_type=dbcsr_type_no_symmetry)
7155 IF (penalty_occ_vol)
THEN
7156 penalty_prefactor_local = penalty_occ_vol_prefactor/(2.0_dp*spin_factor)
7158 penalty_prefactor_local = 0.0_dp
7161 WRITE (unit_nr, *)
"penalty_prefactor_local: ", penalty_prefactor_local
7162 WRITE (unit_nr, *)
"penalty_prefactor_2: ", penalty_occ_vol_pf2
7166 DO col = 1, nblkcols_tot
7169 DO row = 1, nblkrows_tot
7172 row, col, block_p, found)
7174 ao_domain_sizes(col) = ao_domain_sizes(col) + ao_blk_sizes(row)
7179 h_size = h_size + ao_domain_sizes(col)*mo_block_sizes(col)
7183 ALLOCATE (h(h_size, h_size))
7187 lev1_vert_offset = 0
7189 DO row = 1, nblkcols_tot
7191 lev1_hori_offset = 0
7192 DO col = 1, nblkcols_tot
7195 ALLOCATE (f_ao_block(ao_domain_sizes(row), ao_domain_sizes(col)))
7196 ALLOCATE (s_ao_block(ao_domain_sizes(row), ao_domain_sizes(col)))
7197 ALLOCATE (f_mo_block(mo_block_sizes(row), mo_block_sizes(col)))
7198 ALLOCATE (s_mo_block(mo_block_sizes(row), mo_block_sizes(col)))
7200 f_ao_block(:, :) = 0.0_dp
7201 s_ao_block(:, :) = 0.0_dp
7202 f_mo_block(:, :) = 0.0_dp
7203 s_mo_block(:, :) = 0.0_dp
7208 DO block_row = 1, nblkcols_tot
7211 block_row, row, block_p, found_row)
7215 DO block_col = 1, nblkcols_tot
7218 block_col, col, block_p, found_col)
7222 block_row, block_col, block_p, found)
7225 f_ao_block(ao_vert_offset + 1:ao_vert_offset + ao_block_sizes(block_row), &
7226 ao_hori_offset + 1:ao_hori_offset + ao_block_sizes(block_col)) &
7231 block_row, block_col, block_p, found)
7234 s_ao_block(ao_vert_offset + 1:ao_vert_offset + ao_block_sizes(block_row), &
7235 ao_hori_offset + 1:ao_hori_offset + ao_block_sizes(block_col)) &
7239 ao_hori_offset = ao_hori_offset + ao_block_sizes(block_col)
7245 ao_vert_offset = ao_vert_offset + ao_block_sizes(block_row)
7255 f_mo_block(1:mo_block_sizes(row), 1:mo_block_sizes(col)) = block_p(:, :)
7260 s_mo_block(1:mo_block_sizes(row), 1:mo_block_sizes(col)) = block_p(:, :)
7264 lev2_vert_offset = 0
7265 DO orb_j = 1, mo_block_sizes(row)
7267 lev2_hori_offset = 0
7268 DO orb_i = 1, mo_block_sizes(col)
7269 IF (orb_i == orb_j .AND. row == col)
THEN
7270 h(lev1_vert_offset + lev2_vert_offset + 1:lev1_vert_offset + lev2_vert_offset + ao_domain_sizes(row), &
7271 lev1_hori_offset + lev2_hori_offset + 1:lev1_hori_offset + lev2_hori_offset + ao_domain_sizes(col)) &
7272 = f_ao_block(:, :) + s_ao_block(:, :)
7275 lev2_hori_offset = lev2_hori_offset + ao_domain_sizes(col)
7279 lev2_vert_offset = lev2_vert_offset + ao_domain_sizes(row)
7283 lev1_hori_offset = lev1_hori_offset + ao_domain_sizes(col)*mo_block_sizes(col)
7285 DEALLOCATE (f_ao_block)
7286 DEALLOCATE (s_ao_block)
7287 DEALLOCATE (f_mo_block)
7288 DEALLOCATE (s_mo_block)
7292 lev1_vert_offset = lev1_vert_offset + ao_domain_sizes(row)*mo_block_sizes(row)
7302 ALLOCATE (grad_vec(h_size))
7303 grad_vec(:) = 0.0_dp
7304 lev1_vert_offset = 0
7306 DO col = 1, nblkcols_tot
7309 lev2_vert_offset = 0
7310 DO row = 1, nblkrows_tot
7313 row, col, block_p, found_row)
7317 row, col, block_p, found)
7320 DO orb_i = 1, mo_block_sizes(col)
7321 grad_vec(lev1_vert_offset + ao_domain_sizes(col)*(orb_i - 1) + lev2_vert_offset + 1: &
7322 lev1_vert_offset + ao_domain_sizes(col)*(orb_i - 1) + lev2_vert_offset + ao_block_sizes(row)) &
7328 lev2_vert_offset = lev2_vert_offset + ao_block_sizes(row)
7334 lev1_vert_offset = lev1_vert_offset + ao_domain_sizes(col)*mo_block_sizes(col)
7340 ALLOCATE (hinv(h_size, h_size))
7341 hinv(:, :) = h(:, :)
7344 ALLOCATE (eigenvalues(h_size))
7347 ALLOCATE (work(max(1, lwork)))
7348 CALL dsyev(
'V',
'L', h_size, hinv, h_size, eigenvalues, work, lwork, info)
7349 lwork = int(work(1))
7352 ALLOCATE (work(max(1, lwork)))
7353 CALL dsyev(
'V',
'L', h_size, hinv, h_size, eigenvalues, work, lwork, info)
7355 WRITE (unit_nr, *)
'DSYEV ERROR MESSAGE: ', info
7356 cpabort(
"DSYEV failed")
7361 ALLOCATE (step_vec(h_size))
7363 step_vec(:) = matmul(transpose(hinv), grad_vec)
7367 ALLOCATE (test(h_size, h_size))
7370 WRITE (unit_nr,
"(I10,F20.10,F20.10)") jj, eigenvalues(jj), step_vec(jj)
7371 IF (eigenvalues(jj) > eps_zero)
THEN
7372 test(jj, :) = hinv(:, jj)/eigenvalues(jj)
7374 test(jj, :) = hinv(:, jj)*0.0_dp
7375 zero_neg_eiv = zero_neg_eiv + 1
7378 WRITE (unit_nr, *)
'ZERO OR NEGATIVE EIGENVALUES: ', zero_neg_eiv
7379 DEALLOCATE (step_vec)
7381 ALLOCATE (test2(h_size, h_size))
7382 test2(:, :) = matmul(hinv, test)
7383 hinv(:, :) = test2(:, :)
7384 DEALLOCATE (test, test2)
7386 DEALLOCATE (eigenvalues)
7389 ALLOCATE (test(h_size, h_size))
7390 test(:, :) = matmul(hinv, h)
7392 test(ii, ii) = test(ii, ii) - 1.0_dp
7397 test_error = test_error + test(jj, ii)*test(jj, ii)
7400 WRITE (unit_nr, *)
"Hessian inversion error: ", sqrt(test_error)
7404 ALLOCATE (step_vec(h_size))
7405 ALLOCATE (tmp(h_size))
7406 tmp(:) = matmul(hinv, grad_vec)
7407 step_vec(:) = -1.0_dp*tmp(:)
7409 ALLOCATE (tmpr(h_size))
7410 tmpr(:) = matmul(h, step_vec)
7411 tmp(:) = tmpr(:) + grad_vec(:)
7413 WRITE (unit_nr, *)
"NEWTOV step error: ", maxval(abs(tmp))
7419 DEALLOCATE (grad_vec)
7427 template=matrix_grad, &
7428 matrix_type=dbcsr_type_no_symmetry)
7431 lev1_vert_offset = 0
7433 DO col = 1, nblkcols_tot
7436 lev2_vert_offset = 0
7437 DO row = 1, nblkrows_tot
7440 row, col, block_p, found_row)
7443 ALLOCATE (new_block(ao_block_sizes(row), mo_block_sizes(col)))
7444 DO orb_i = 1, mo_block_sizes(col)
7445 new_block(:, orb_i) = &
7446 step_vec(lev1_vert_offset + ao_domain_sizes(col)*(orb_i - 1) + lev2_vert_offset + 1: &
7447 lev1_vert_offset + ao_domain_sizes(col)*(orb_i - 1) + lev2_vert_offset + ao_block_sizes(row))
7450 DEALLOCATE (new_block)
7451 lev2_vert_offset = lev2_vert_offset + ao_block_sizes(row)
7456 lev1_vert_offset = lev1_vert_offset + ao_domain_sizes(col)*mo_block_sizes(col)
7460 DEALLOCATE (step_vec)
7464 DEALLOCATE (mo_block_sizes, ao_block_sizes)
7465 DEALLOCATE (ao_domain_sizes)
7468 template=quench_t, &
7469 matrix_type=dbcsr_type_no_symmetry)
7474 0.0_dp, matrix_s_ao_sym, &
7475 retain_sparsity=.true.)
7477 template=quench_t, &
7478 matrix_type=dbcsr_type_no_symmetry)
7483 0.0_dp, matrix_f_ao_sym, &
7484 retain_sparsity=.true.)
7485 CALL dbcsr_add(matrix_s_ao_sym, matrix_f_ao_sym, &
7487 CALL dbcsr_scale(matrix_s_ao_sym, 2.0_dp*spin_factor)
7488 CALL dbcsr_add(matrix_s_ao_sym, matrix_grad, &
7491 WRITE (unit_nr, *)
"NEWTOL step error: ", test_error
7495 CALL timestop(handle)
7497 END SUBROUTINE hessian_diag_apply
7517 matrix_t_in, matrix_t_out, perturbation_only, &
7523 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:) :: quench_t, matrix_t_in, matrix_t_out
7524 LOGICAL,
INTENT(IN) :: perturbation_only
7525 INTEGER,
INTENT(IN),
OPTIONAL :: special_case
7527 CHARACTER(len=*),
PARAMETER :: routinen =
'almo_scf_xalmo_trustr'
7529 INTEGER :: handle, ispin, iteration, iteration_type_to_report, my_special_case, ndomains, &
7530 nspins, outer_iteration, prec_type, unit_nr
7531 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: nocc
7532 LOGICAL :: assume_t0_q0x, border_reached, inner_loop_success, normalize_orbitals, &
7533 optimize_theta, penalty_occ_vol, reset_conjugator, same_position, scf_converged
7534 REAL(kind=
dp) :: beta, energy_start, energy_trial, eta, expected_reduction, &
7535 fake_step_size_to_report, grad_norm_ratio, grad_norm_ref, loss_change_to_report, &
7536 loss_start, loss_trial, model_grad_norm, penalty_amplitude, penalty_start, penalty_trial, &
7537 radius_current, radius_max, real_temp, rho, spin_factor, step_norm, step_size, t1, &
7538 t1outer, t2, t2outer, y_scalar
7539 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: grad_norm_spin, &
7540 penalty_occ_vol_g_prefactor, &
7541 penalty_occ_vol_h_prefactor
7544 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:) :: ftsiginv, grad, m_model_bd, m_model_d, &
7545 m_model_hessian, m_model_hessian_inv, m_model_r, m_model_r_prev, m_model_rt, &
7546 m_model_rt_prev, m_sig_sqrti_ii, m_theta, m_theta_trial, prev_step, siginvtftsiginv, st, &
7549 DIMENSION(:, :) :: domain_model_hessian_inv, domain_r_down
7552 CALL timeset(routinen, handle)
7557 IF (
PRESENT(special_case)) my_special_case = special_case
7561 IF (logger%para_env%is_source())
THEN
7568 assume_t0_q0x = .false.
7570 optimize_theta = .false.
7572 nspins = almo_scf_env%nspins
7573 IF (nspins == 1)
THEN
7574 spin_factor = 2.0_dp
7576 spin_factor = 1.0_dp
7579 IF (unit_nr > 0)
THEN
7581 SELECT CASE (my_special_case)
7583 WRITE (unit_nr,
'(T2,A,A,A)') repeat(
"-", 20), &
7584 " Optimization of block-diagonal ALMOs ", repeat(
"-", 21)
7586 WRITE (unit_nr,
'(T2,A,A,A)') repeat(
"-", 20), &
7587 " Optimization of fully delocalized MOs ", repeat(
"-", 20)
7589 WRITE (unit_nr,
'(T2,A,A,A)') repeat(
"-", 27), &
7590 " Optimization of XALMOs ", repeat(
"-", 28)
7593 CALL trust_r_report(unit_nr, &
7598 delta_loss=0.0_dp, &
7600 predicted_reduction=0.0_dp, &
7604 WRITE (unit_nr,
'(T2,A)') repeat(
"-", 79)
7608 penalty_occ_vol = .false.
7609 normalize_orbitals = penalty_occ_vol
7610 penalty_amplitude = 0.0_dp
7611 ALLOCATE (penalty_occ_vol_g_prefactor(nspins))
7612 ALLOCATE (penalty_occ_vol_h_prefactor(nspins))
7613 penalty_occ_vol_g_prefactor(:) = 0.0_dp
7614 penalty_occ_vol_h_prefactor(:) = 0.0_dp
7617 prec_type = optimizer%preconditioner
7619 ALLOCATE (grad_norm_spin(nspins))
7620 ALLOCATE (nocc(nspins))
7624 ALLOCATE (m_theta(nspins))
7625 DO ispin = 1, nspins
7627 template=matrix_t_out(ispin), &
7628 matrix_type=dbcsr_type_no_symmetry)
7633 m_t_in=matrix_t_in, &
7634 m_t0=almo_scf_env%matrix_t_blk, &
7635 m_quench_t=quench_t, &
7636 m_overlap=almo_scf_env%matrix_s(1), &
7637 m_sigma_tmpl=almo_scf_env%matrix_sigma_inv, &
7639 xalmo_history=almo_scf_env%xalmo_history, &
7640 assume_t0_q0x=assume_t0_q0x, &
7641 optimize_theta=optimize_theta, &
7642 envelope_amplitude=almo_scf_env%envelope_amplitude, &
7643 eps_filter=almo_scf_env%eps_filter, &
7644 order_lanczos=almo_scf_env%order_lanczos, &
7645 eps_lanczos=almo_scf_env%eps_lanczos, &
7646 max_iter_lanczos=almo_scf_env%max_iter_lanczos, &
7647 nocc_of_domain=almo_scf_env%nocc_of_domain)
7649 ndomains = almo_scf_env%ndomains
7650 ALLOCATE (domain_r_down(ndomains, nspins))
7652 ALLOCATE (domain_model_hessian_inv(ndomains, nspins))
7655 ALLOCATE (m_model_hessian(nspins))
7656 ALLOCATE (m_model_hessian_inv(nspins))
7657 ALLOCATE (siginvtftsiginv(nspins))
7658 ALLOCATE (stsiginv_0(nspins))
7659 ALLOCATE (ftsiginv(nspins))
7660 ALLOCATE (st(nspins))
7661 ALLOCATE (grad(nspins))
7662 ALLOCATE (prev_step(nspins))
7663 ALLOCATE (step(nspins))
7664 ALLOCATE (m_sig_sqrti_ii(nspins))
7665 ALLOCATE (m_model_r(nspins))
7666 ALLOCATE (m_model_rt(nspins))
7667 ALLOCATE (m_model_d(nspins))
7668 ALLOCATE (m_model_bd(nspins))
7669 ALLOCATE (m_model_r_prev(nspins))
7670 ALLOCATE (m_model_rt_prev(nspins))
7671 ALLOCATE (m_theta_trial(nspins))
7673 DO ispin = 1, nspins
7677 template=almo_scf_env%matrix_ks(ispin), &
7678 matrix_type=dbcsr_type_no_symmetry)
7680 template=almo_scf_env%matrix_ks(ispin), &
7681 matrix_type=dbcsr_type_no_symmetry)
7683 template=almo_scf_env%matrix_sigma(ispin), &
7684 matrix_type=dbcsr_type_no_symmetry)
7686 template=matrix_t_out(ispin), &
7687 matrix_type=dbcsr_type_no_symmetry)
7689 template=matrix_t_out(ispin), &
7690 matrix_type=dbcsr_type_no_symmetry)
7692 template=matrix_t_out(ispin), &
7693 matrix_type=dbcsr_type_no_symmetry)
7695 template=matrix_t_out(ispin), &
7696 matrix_type=dbcsr_type_no_symmetry)
7698 template=matrix_t_out(ispin), &
7699 matrix_type=dbcsr_type_no_symmetry)
7701 template=matrix_t_out(ispin), &
7702 matrix_type=dbcsr_type_no_symmetry)
7704 template=almo_scf_env%matrix_sigma_inv(ispin), &
7705 matrix_type=dbcsr_type_no_symmetry)
7707 template=matrix_t_out(ispin), &
7708 matrix_type=dbcsr_type_no_symmetry)
7710 template=matrix_t_out(ispin), &
7711 matrix_type=dbcsr_type_no_symmetry)
7713 template=matrix_t_out(ispin), &
7714 matrix_type=dbcsr_type_no_symmetry)
7716 template=matrix_t_out(ispin), &
7717 matrix_type=dbcsr_type_no_symmetry)
7719 template=matrix_t_out(ispin), &
7720 matrix_type=dbcsr_type_no_symmetry)
7722 template=matrix_t_out(ispin), &
7723 matrix_type=dbcsr_type_no_symmetry)
7725 template=matrix_t_out(ispin), &
7726 matrix_type=dbcsr_type_no_symmetry)
7729 CALL dbcsr_set(prev_step(ispin), 0.0_dp)
7732 nfullrows_total=nocc(ispin))
7740 matrix_s=almo_scf_env%matrix_s(1), &
7741 subm_s_inv=almo_scf_env%domain_s_inv(:, ispin), &
7742 dpattern=quench_t(ispin), &
7743 map=almo_scf_env%domain_map(ispin), &
7744 node_of_domain=almo_scf_env%cpu_of_domain)
7754 template=almo_scf_env%matrix_s(1), &
7755 matrix_type=dbcsr_type_no_symmetry)
7757 almo_scf_env%matrix_s_blk(1), &
7758 threshold=almo_scf_env%eps_filter, &
7759 filter_eps=almo_scf_env%eps_filter)
7765 template=almo_scf_env%matrix_s(1), &
7766 matrix_type=dbcsr_type_no_symmetry)
7769 para_env=almo_scf_env%para_env, &
7770 blacs_env=almo_scf_env%blacs_env)
7772 para_env=almo_scf_env%para_env, &
7773 blacs_env=almo_scf_env%blacs_env, &
7774 uplo_to_full=.true.)
7779 radius_max = optimizer%max_trust_radius
7780 radius_current = min(optimizer%initial_trust_radius, radius_max)
7782 eta = min(max(optimizer%rho_do_not_update, 0.0_dp), 0.25_dp)
7783 energy_start = 0.0_dp
7784 energy_trial = 0.0_dp
7785 penalty_start = 0.0_dp
7786 penalty_trial = 0.0_dp
7790 same_position = .false.
7793 CALL main_var_to_xalmos_and_loss_func( &
7794 almo_scf_env=almo_scf_env, &
7796 m_main_var_in=m_theta, &
7797 m_t_out=matrix_t_out, &
7798 m_sig_sqrti_ii_out=m_sig_sqrti_ii, &
7799 energy_out=energy_start, &
7800 penalty_out=penalty_start, &
7801 m_ftsiginv_out=ftsiginv, &
7802 m_siginvtftsiginv_out=siginvtftsiginv, &
7804 m_stsiginv0_in=stsiginv_0, &
7805 m_quench_t_in=quench_t, &
7806 domain_r_down_in=domain_r_down, &
7807 assume_t0_q0x=assume_t0_q0x, &
7808 just_started=.true., &
7809 optimize_theta=optimize_theta, &
7810 normalize_orbitals=normalize_orbitals, &
7811 perturbation_only=perturbation_only, &
7812 do_penalty=penalty_occ_vol, &
7813 special_case=my_special_case)
7814 loss_start = energy_start + penalty_start
7816 almo_scf_env%almo_scf_energy = energy_start
7818 DO ispin = 1, nspins
7819 IF (penalty_occ_vol)
THEN
7820 penalty_occ_vol_g_prefactor(ispin) = &
7821 -2.0_dp*penalty_amplitude*spin_factor*nocc(ispin)
7822 penalty_occ_vol_h_prefactor(ispin) = 0.0_dp
7827 scf_converged = .false.
7828 adjust_r_loop:
DO outer_iteration = 1, optimizer%max_iter_outer_loop
7831 border_reached = .false.
7833 DO ispin = 1, nspins
7835 CALL dbcsr_filter(step(ispin), almo_scf_env%eps_filter)
7838 IF (.NOT. same_position)
THEN
7840 DO ispin = 1, nspins
7842 IF (unit_nr > 0 .AND. debug_mode)
WRITE (unit_nr, *)
"...Compute model gradient"
7843 CALL compute_gradient( &
7844 m_grad_out=grad(ispin), &
7845 m_ks=almo_scf_env%matrix_ks(ispin), &
7846 m_s=almo_scf_env%matrix_s(1), &
7847 m_t=matrix_t_out(ispin), &
7848 m_t0=almo_scf_env%matrix_t_blk(ispin), &
7849 m_siginv=almo_scf_env%matrix_sigma_inv(ispin), &
7850 m_quench_t=quench_t(ispin), &
7851 m_ftsiginv=ftsiginv(ispin), &
7852 m_siginvtftsiginv=siginvtftsiginv(ispin), &
7854 m_stsiginv0=stsiginv_0(ispin), &
7855 m_theta=m_theta(ispin), &
7856 m_sig_sqrti_ii=m_sig_sqrti_ii(ispin), &
7857 domain_s_inv=almo_scf_env%domain_s_inv(:, ispin), &
7858 domain_r_down=domain_r_down(:, ispin), &
7859 cpu_of_domain=almo_scf_env%cpu_of_domain, &
7860 domain_map=almo_scf_env%domain_map(ispin), &
7861 assume_t0_q0x=assume_t0_q0x, &
7862 optimize_theta=optimize_theta, &
7863 normalize_orbitals=normalize_orbitals, &
7864 penalty_occ_vol=penalty_occ_vol, &
7865 penalty_occ_vol_prefactor=penalty_occ_vol_g_prefactor(ispin), &
7866 envelope_amplitude=almo_scf_env%envelope_amplitude, &
7867 eps_filter=almo_scf_env%eps_filter, &
7868 spin_factor=spin_factor, &
7869 special_case=my_special_case)
7876 DO ispin = 1, nspins
7879 grad_norm_ref = maxval(grad_norm_spin)
7882 CALL trust_r_report(unit_nr, &
7884 iteration=outer_iteration, &
7886 delta_loss=0.0_dp, &
7887 grad_norm=grad_norm_ref, &
7888 predicted_reduction=0.0_dp, &
7890 radius=radius_current, &
7891 new=.NOT. same_position, &
7892 time=t2outer - t1outer)
7895 IF (grad_norm_ref <= optimizer%eps_error)
THEN
7896 scf_converged = .true.
7897 border_reached = .false.
7898 expected_reduction = 0.0_dp
7899 IF (.NOT. (optimizer%early_stopping_on .AND. outer_iteration == 1))
THEN
7903 scf_converged = .false.
7906 DO ispin = 1, nspins
7908 CALL dbcsr_copy(m_model_r(ispin), grad(ispin))
7914 IF (unit_nr > 0 .AND. debug_mode)
WRITE (unit_nr, *)
"...Multiply Sinv.r"
7918 0.0_dp, m_model_rt(ispin), &
7919 filter_eps=almo_scf_env%eps_filter)
7923 IF (unit_nr > 0 .AND. debug_mode)
WRITE (unit_nr, *)
"...Multiply Sinv_xx.r"
7925 matrix_in=m_model_r(ispin), &
7926 matrix_out=m_model_rt(ispin), &
7927 operator1=almo_scf_env%domain_s_inv(:, ispin), &
7928 dpattern=quench_t(ispin), &
7929 map=almo_scf_env%domain_map(ispin), &
7930 node_of_domain=almo_scf_env%cpu_of_domain, &
7932 filter_eps=almo_scf_env%eps_filter)
7935 cpabort(
"Unknown XALMO special case")
7938 CALL dbcsr_copy(m_model_d(ispin), m_model_rt(ispin))
7943 IF (.NOT. same_position)
THEN
7945 SELECT CASE (prec_type)
7948 IF (unit_nr > 0 .AND. debug_mode)
WRITE (unit_nr, *)
"...Compute model Hessian"
7949 DO ispin = 1, nspins
7950 CALL compute_preconditioner( &
7951 domain_prec_out=almo_scf_env%domain_preconditioner(:, ispin), &
7952 m_prec_out=m_model_hessian(ispin), &
7953 m_ks=almo_scf_env%matrix_ks(ispin), &
7954 m_s=almo_scf_env%matrix_s(1), &
7955 m_siginv=almo_scf_env%matrix_sigma_inv(ispin), &
7956 m_quench_t=quench_t(ispin), &
7957 m_ftsiginv=ftsiginv(ispin), &
7958 m_siginvtftsiginv=siginvtftsiginv(ispin), &
7960 para_env=almo_scf_env%para_env, &
7961 blacs_env=almo_scf_env%blacs_env, &
7962 nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin), &
7963 domain_s_inv=almo_scf_env%domain_s_inv(:, ispin), &
7964 domain_r_down=domain_r_down(:, ispin), &
7965 cpu_of_domain=almo_scf_env%cpu_of_domain, &
7966 domain_map=almo_scf_env%domain_map(ispin), &
7967 assume_t0_q0x=.false., &
7968 penalty_occ_vol=penalty_occ_vol, &
7969 penalty_occ_vol_prefactor=penalty_occ_vol_g_prefactor(ispin), &
7970 eps_filter=almo_scf_env%eps_filter, &
7972 spin_factor=spin_factor, &
7973 skip_inversion=.true., &
7974 special_case=my_special_case)
7979 cpabort(
"Unknown preconditioner")
7986 CALL fixed_r_report(unit_nr, &
7990 border_reached=.false., &
7992 grad_norm_ratio=0.0_dp, &
7995 IF (unit_nr > 0 .AND. debug_mode)
WRITE (unit_nr, *)
"...Start inner loop"
7998 inner_loop_success = .false.
8000 fixed_r_loop:
DO iteration = 1, optimizer%max_iter
8004 DO ispin = 1, nspins
8011 m_model_hessian(ispin), &
8013 0.0_dp, m_model_bd(ispin), &
8014 filter_eps=almo_scf_env%eps_filter)
8019 matrix_in=m_model_d(ispin), &
8020 matrix_out=m_model_bd(ispin), &
8021 operator1=almo_scf_env%domain_preconditioner(:, ispin), &
8022 dpattern=quench_t(ispin), &
8023 map=almo_scf_env%domain_map(ispin), &
8024 node_of_domain=almo_scf_env%cpu_of_domain, &
8026 filter_eps=almo_scf_env%eps_filter)
8031 CALL dbcsr_dot(m_model_d(ispin), m_model_bd(ispin), real_temp)
8032 y_scalar = y_scalar + real_temp
8035 IF (unit_nr > 0 .AND. debug_mode)
WRITE (unit_nr, *)
"...Curvature: ", y_scalar
8038 IF (y_scalar < 0.0_dp)
THEN
8040 CALL step_size_to_border( &
8041 step_size_out=step_size, &
8042 metric_in=almo_scf_env%matrix_s, &
8044 direction_in=m_model_d, &
8045 trust_radius_in=radius_current, &
8046 quench_t_in=quench_t, &
8047 eps_filter_in=almo_scf_env%eps_filter &
8050 DO ispin = 1, nspins
8051 CALL dbcsr_add(step(ispin), m_model_d(ispin), 1.0_dp, step_size)
8054 border_reached = .true.
8055 inner_loop_success = .true.
8057 CALL predicted_reduction( &
8058 reduction_out=expected_reduction, &
8061 hess_in=m_model_hessian, &
8062 hess_submatrix_in=almo_scf_env%domain_preconditioner, &
8063 quench_t_in=quench_t, &
8064 special_case=my_special_case, &
8065 eps_filter=almo_scf_env%eps_filter, &
8066 domain_map=almo_scf_env%domain_map, &
8067 cpu_of_domain=almo_scf_env%cpu_of_domain &
8071 CALL fixed_r_report(unit_nr, &
8073 iteration=iteration, &
8074 step_size=step_size, &
8075 border_reached=border_reached, &
8076 curvature=y_scalar, &
8077 grad_norm_ratio=expected_reduction, &
8086 DO ispin = 1, nspins
8087 CALL dbcsr_dot(m_model_r(ispin), m_model_rt(ispin), real_temp)
8088 step_size = step_size + real_temp
8090 step_size = step_size/y_scalar
8091 IF (unit_nr > 0 .AND. debug_mode)
WRITE (unit_nr, *)
"...Proposed step size: ", step_size
8094 DO ispin = 1, nspins
8095 CALL dbcsr_copy(prev_step(ispin), step(ispin))
8096 CALL dbcsr_add(step(ispin), m_model_d(ispin), 1.0_dp, step_size)
8100 CALL contravariant_matrix_norm( &
8101 norm_out=step_norm, &
8103 metric_in=almo_scf_env%matrix_s, &
8104 quench_t_in=quench_t, &
8105 eps_filter_in=almo_scf_env%eps_filter &
8107 IF (unit_nr > 0 .AND. debug_mode)
WRITE (unit_nr, *)
"...Step norm: ", step_norm
8110 IF (step_norm > radius_current)
THEN
8112 IF (unit_nr > 0 .AND. debug_mode)
WRITE (unit_nr, *)
"...Norm is too large"
8113 CALL step_size_to_border( &
8114 step_size_out=step_size, &
8115 metric_in=almo_scf_env%matrix_s, &
8116 position_in=prev_step, &
8117 direction_in=m_model_d, &
8118 trust_radius_in=radius_current, &
8119 quench_t_in=quench_t, &
8120 eps_filter_in=almo_scf_env%eps_filter &
8122 IF (unit_nr > 0 .AND. debug_mode)
WRITE (unit_nr, *)
"...Step size to border: ", step_size
8124 DO ispin = 1, nspins
8125 CALL dbcsr_copy(step(ispin), prev_step(ispin))
8126 CALL dbcsr_add(step(ispin), m_model_d(ispin), 1.0_dp, step_size)
8129 IF (debug_mode)
THEN
8131 IF (unit_nr > 0)
WRITE (unit_nr, *)
"...Extra norm evaluation"
8132 CALL contravariant_matrix_norm( &
8133 norm_out=step_norm, &
8135 metric_in=almo_scf_env%matrix_s, &
8136 quench_t_in=quench_t, &
8137 eps_filter_in=almo_scf_env%eps_filter &
8139 IF (unit_nr > 0)
WRITE (unit_nr, *)
"...Step norm: ", step_norm
8140 IF (unit_nr > 0)
WRITE (unit_nr, *)
"...Current radius: ", radius_current
8143 border_reached = .true.
8144 inner_loop_success = .true.
8146 CALL predicted_reduction( &
8147 reduction_out=expected_reduction, &
8150 hess_in=m_model_hessian, &
8151 hess_submatrix_in=almo_scf_env%domain_preconditioner, &
8152 quench_t_in=quench_t, &
8153 special_case=my_special_case, &
8154 eps_filter=almo_scf_env%eps_filter, &
8155 domain_map=almo_scf_env%domain_map, &
8156 cpu_of_domain=almo_scf_env%cpu_of_domain &
8160 CALL fixed_r_report(unit_nr, &
8162 iteration=iteration, &
8163 step_size=step_size, &
8164 border_reached=border_reached, &
8165 curvature=y_scalar, &
8166 grad_norm_ratio=expected_reduction, &
8176 border_reached = .false.
8177 inner_loop_success = .true.
8179 CALL predicted_reduction( &
8180 reduction_out=expected_reduction, &
8183 hess_in=m_model_hessian, &
8184 hess_submatrix_in=almo_scf_env%domain_preconditioner, &
8185 quench_t_in=quench_t, &
8186 special_case=my_special_case, &
8187 eps_filter=almo_scf_env%eps_filter, &
8188 domain_map=almo_scf_env%domain_map, &
8189 cpu_of_domain=almo_scf_env%cpu_of_domain &
8193 CALL fixed_r_report(unit_nr, &
8195 iteration=iteration, &
8196 step_size=step_size, &
8197 border_reached=border_reached, &
8198 curvature=y_scalar, &
8199 grad_norm_ratio=expected_reduction, &
8207 SELECT CASE (prec_type)
8210 IF (unit_nr > 0 .AND. debug_mode)
WRITE (unit_nr, *)
"...Pseudo-invert model Hessian"
8213 DO ispin = 1, nspins
8215 matrix_in=m_model_hessian(ispin), &
8216 matrix_out=m_model_hessian_inv(ispin), &
8217 nocc=almo_scf_env%nocc_of_domain(:, ispin) &
8224 DO ispin = 1, nspins
8225 CALL dbcsr_copy(m_model_hessian_inv(ispin), &
8226 m_model_hessian(ispin))
8228 para_env=almo_scf_env%para_env, &
8229 blacs_env=almo_scf_env%blacs_env)
8231 para_env=almo_scf_env%para_env, &
8232 blacs_env=almo_scf_env%blacs_env, &
8233 uplo_to_full=.true.)
8235 almo_scf_env%eps_filter)
8240 DO ispin = 1, nspins
8242 matrix_main=m_model_hessian(ispin), &
8243 subm_s_inv=almo_scf_env%domain_s_inv(:, ispin), &
8244 subm_r_down=domain_r_down(:, ispin), &
8245 matrix_trimmer=quench_t(ispin), &
8246 dpattern=quench_t(ispin), &
8247 map=almo_scf_env%domain_map(ispin), &
8248 node_of_domain=almo_scf_env%cpu_of_domain, &
8250 use_trimmer=.false., &
8252 skip_inversion=.false. &
8260 cpabort(
"Unknown preconditioner")
8265 DO ispin = 1, nspins
8272 m_model_hessian_inv(ispin), &
8274 0.0_dp, m_model_bd(ispin), &
8275 filter_eps=almo_scf_env%eps_filter)
8280 matrix_in=m_model_r(ispin), &
8281 matrix_out=m_model_bd(ispin), &
8282 operator1=domain_model_hessian_inv(:, ispin), &
8283 dpattern=quench_t(ispin), &
8284 map=almo_scf_env%domain_map(ispin), &
8285 node_of_domain=almo_scf_env%cpu_of_domain, &
8287 filter_eps=almo_scf_env%eps_filter)
8294 CALL contravariant_matrix_norm( &
8295 norm_out=step_norm, &
8296 matrix_in=m_model_bd, &
8297 metric_in=almo_scf_env%matrix_s, &
8298 quench_t_in=quench_t, &
8299 eps_filter_in=almo_scf_env%eps_filter &
8301 IF (unit_nr > 0 .AND. debug_mode)
WRITE (unit_nr, *)
"...pB norm: ", step_norm
8304 IF (step_norm <= radius_current)
THEN
8306 IF (unit_nr > 0 .AND. debug_mode)
WRITE (unit_nr, *)
"...Full dogleg"
8308 border_reached = .false.
8310 DO ispin = 1, nspins
8311 CALL dbcsr_copy(step(ispin), m_model_bd(ispin))
8314 fake_step_size_to_report = 2.0_dp
8315 iteration_type_to_report = 6
8319 IF (unit_nr > 0 .AND. debug_mode)
WRITE (unit_nr, *)
"...pB norm is too large"
8321 border_reached = .true.
8325 DO ispin = 1, nspins
8326 CALL dbcsr_add(m_model_bd(ispin), step(ispin), 1.0_dp, -1.0_dp)
8329 CALL step_size_to_border( &
8330 step_size_out=step_size, &
8331 metric_in=almo_scf_env%matrix_s, &
8333 direction_in=m_model_bd, &
8334 trust_radius_in=radius_current, &
8335 quench_t_in=quench_t, &
8336 eps_filter_in=almo_scf_env%eps_filter &
8338 IF (unit_nr > 0 .AND. debug_mode)
WRITE (unit_nr, *)
"...Step size to border: ", step_size
8339 IF (step_size > 1.0_dp .OR. step_size < 0.0_dp)
THEN
8340 IF (unit_nr > 0)
THEN
8341 WRITE (unit_nr, *)
"Step size (", step_size,
") must lie inside (0,1)"
8343 cpabort(
"Wrong dog leg step. We should never end up here.")
8346 DO ispin = 1, nspins
8347 CALL dbcsr_add(step(ispin), m_model_bd(ispin), 1.0_dp, step_size)
8350 fake_step_size_to_report = 1.0_dp + step_size
8351 iteration_type_to_report = 7
8355 IF (debug_mode)
THEN
8357 IF (unit_nr > 0)
WRITE (unit_nr, *)
"...Extra norm evaluation"
8358 CALL contravariant_matrix_norm( &
8359 norm_out=step_norm, &
8361 metric_in=almo_scf_env%matrix_s, &
8362 quench_t_in=quench_t, &
8363 eps_filter_in=almo_scf_env%eps_filter &
8365 IF (unit_nr > 0)
WRITE (unit_nr, *)
"...Step norm: ", step_norm
8366 IF (unit_nr > 0)
WRITE (unit_nr, *)
"...Current radius: ", radius_current
8369 CALL predicted_reduction( &
8370 reduction_out=expected_reduction, &
8373 hess_in=m_model_hessian, &
8374 hess_submatrix_in=almo_scf_env%domain_preconditioner, &
8375 quench_t_in=quench_t, &
8376 special_case=my_special_case, &
8377 eps_filter=almo_scf_env%eps_filter, &
8378 domain_map=almo_scf_env%domain_map, &
8379 cpu_of_domain=almo_scf_env%cpu_of_domain &
8382 inner_loop_success = .true.
8385 CALL fixed_r_report(unit_nr, &
8386 iter_type=iteration_type_to_report, &
8387 iteration=iteration, &
8388 step_size=fake_step_size_to_report, &
8389 border_reached=border_reached, &
8390 curvature=y_scalar, &
8391 grad_norm_ratio=expected_reduction, &
8399 DO ispin = 1, nspins
8401 CALL dbcsr_copy(m_model_r_prev(ispin), m_model_r(ispin))
8402 CALL dbcsr_add(m_model_r(ispin), m_model_bd(ispin), &
8407 DO ispin = 1, nspins
8410 model_grad_norm = maxval(grad_norm_spin)
8413 grad_norm_ratio = model_grad_norm/grad_norm_ref
8414 IF (grad_norm_ratio < optimizer%model_grad_norm_ratio)
THEN
8416 border_reached = .false.
8417 inner_loop_success = .true.
8419 CALL predicted_reduction( &
8420 reduction_out=expected_reduction, &
8423 hess_in=m_model_hessian, &
8424 hess_submatrix_in=almo_scf_env%domain_preconditioner, &
8425 quench_t_in=quench_t, &
8426 special_case=my_special_case, &
8427 eps_filter=almo_scf_env%eps_filter, &
8428 domain_map=almo_scf_env%domain_map, &
8429 cpu_of_domain=almo_scf_env%cpu_of_domain &
8433 CALL fixed_r_report(unit_nr, &
8435 iteration=iteration, &
8436 step_size=step_size, &
8437 border_reached=border_reached, &
8438 curvature=y_scalar, &
8439 grad_norm_ratio=expected_reduction, &
8447 DO ispin = 1, nspins
8449 CALL dbcsr_copy(m_model_rt_prev(ispin), m_model_rt(ispin))
8452 DO ispin = 1, nspins
8460 0.0_dp, m_model_rt(ispin), &
8461 filter_eps=almo_scf_env%eps_filter)
8466 matrix_in=m_model_r(ispin), &
8467 matrix_out=m_model_rt(ispin), &
8468 operator1=almo_scf_env%domain_s_inv(:, ispin), &
8469 dpattern=quench_t(ispin), &
8470 map=almo_scf_env%domain_map(ispin), &
8471 node_of_domain=almo_scf_env%cpu_of_domain, &
8473 filter_eps=almo_scf_env%eps_filter)
8479 CALL compute_cg_beta( &
8481 reset_conjugator=reset_conjugator, &
8482 conjugator=optimizer%conjugator, &
8483 grad=m_model_r(:), &
8484 prev_grad=m_model_r_prev(:), &
8485 step=m_model_rt(:), &
8486 prev_step=m_model_rt_prev(:) &
8489 DO ispin = 1, nspins
8491 CALL dbcsr_add(m_model_d(ispin), m_model_rt(ispin), beta, 1.0_dp)
8495 CALL fixed_r_report(unit_nr, &
8497 iteration=iteration, &
8498 step_size=step_size, &
8499 border_reached=border_reached, &
8500 curvature=y_scalar, &
8501 grad_norm_ratio=grad_norm_ratio, &
8510 IF (.NOT. inner_loop_success)
THEN
8511 cpabort(
"Inner loop did not produce solution")
8514 DO ispin = 1, nspins
8516 CALL dbcsr_copy(m_theta_trial(ispin), m_theta(ispin))
8517 CALL dbcsr_add(m_theta_trial(ispin), step(ispin), 1.0_dp, 1.0_dp)
8522 CALL main_var_to_xalmos_and_loss_func( &
8523 almo_scf_env=almo_scf_env, &
8525 m_main_var_in=m_theta_trial, &
8526 m_t_out=matrix_t_out, &
8527 m_sig_sqrti_ii_out=m_sig_sqrti_ii, &
8528 energy_out=energy_trial, &
8529 penalty_out=penalty_trial, &
8530 m_ftsiginv_out=ftsiginv, &
8531 m_siginvtftsiginv_out=siginvtftsiginv, &
8533 m_stsiginv0_in=stsiginv_0, &
8534 m_quench_t_in=quench_t, &
8535 domain_r_down_in=domain_r_down, &
8536 assume_t0_q0x=assume_t0_q0x, &
8537 just_started=.false., &
8538 optimize_theta=optimize_theta, &
8539 normalize_orbitals=normalize_orbitals, &
8540 perturbation_only=perturbation_only, &
8541 do_penalty=penalty_occ_vol, &
8542 special_case=my_special_case)
8543 loss_trial = energy_trial + penalty_trial
8545 rho = (loss_trial - loss_start)/expected_reduction
8546 loss_change_to_report = loss_trial - loss_start
8548 IF (rho < 0.25_dp)
THEN
8549 radius_current = 0.25_dp*radius_current
8551 IF (rho > 0.75_dp .AND. border_reached)
THEN
8552 radius_current = min(2.0_dp*radius_current, radius_max)
8557 DO ispin = 1, nspins
8558 CALL dbcsr_copy(m_theta(ispin), m_theta_trial(ispin))
8560 loss_start = loss_trial
8561 energy_start = energy_trial
8562 penalty_start = penalty_trial
8563 same_position = .false.
8565 almo_scf_env%almo_scf_energy = energy_trial
8568 same_position = .true.
8570 almo_scf_env%almo_scf_energy = energy_start
8575 CALL trust_r_report(unit_nr, &
8577 iteration=outer_iteration, &
8579 delta_loss=loss_change_to_report, &
8581 predicted_reduction=expected_reduction, &
8583 radius=radius_current, &
8584 new=.NOT. same_position, &
8585 time=t2outer - t1outer)
8588 END DO adjust_r_loop
8591 IF (scf_converged)
THEN
8593 CALL wrap_up_xalmo_scf( &
8595 almo_scf_env=almo_scf_env, &
8596 perturbation_in=perturbation_only, &
8597 m_xalmo_in=matrix_t_out, &
8598 m_quench_in=quench_t, &
8599 energy_inout=energy_start)
8603 DO ispin = 1, nspins
8631 DEALLOCATE (m_model_hessian)
8632 DEALLOCATE (m_model_hessian_inv)
8633 DEALLOCATE (siginvtftsiginv)
8634 DEALLOCATE (stsiginv_0)
8635 DEALLOCATE (ftsiginv)
8638 DEALLOCATE (prev_step)
8640 DEALLOCATE (m_sig_sqrti_ii)
8641 DEALLOCATE (m_model_r)
8642 DEALLOCATE (m_model_rt)
8643 DEALLOCATE (m_model_d)
8644 DEALLOCATE (m_model_bd)
8645 DEALLOCATE (m_model_r_prev)
8646 DEALLOCATE (m_model_rt_prev)
8647 DEALLOCATE (m_theta_trial)
8649 DEALLOCATE (domain_r_down)
8650 DEALLOCATE (domain_model_hessian_inv)
8652 DEALLOCATE (penalty_occ_vol_g_prefactor)
8653 DEALLOCATE (penalty_occ_vol_h_prefactor)
8654 DEALLOCATE (grad_norm_spin)
8657 DEALLOCATE (m_theta)
8659 IF (.NOT. scf_converged .AND. .NOT. optimizer%early_stopping_on)
THEN
8660 cpabort(
"Optimization not converged! ")
8663 CALL timestop(handle)
8696 SUBROUTINE main_var_to_xalmos_and_loss_func(almo_scf_env, qs_env, m_main_var_in, &
8697 m_t_out, energy_out, penalty_out, m_sig_sqrti_ii_out, m_FTsiginv_out, &
8698 m_siginvTFTsiginv_out, m_ST_out, m_STsiginv0_in, m_quench_t_in, domain_r_down_in, &
8699 assume_t0_q0x, just_started, optimize_theta, normalize_orbitals, perturbation_only, &
8700 do_penalty, special_case)
8704 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(IN) :: m_main_var_in
8705 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(INOUT) :: m_t_out
8706 REAL(kind=
dp),
INTENT(OUT) :: energy_out, penalty_out
8707 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(INOUT) :: m_sig_sqrti_ii_out, m_ftsiginv_out, &
8708 m_siginvtftsiginv_out, m_st_out, &
8709 m_stsiginv0_in, m_quench_t_in
8711 INTENT(IN) :: domain_r_down_in
8712 LOGICAL,
INTENT(IN) :: assume_t0_q0x, just_started, &
8713 optimize_theta, normalize_orbitals, &
8714 perturbation_only, do_penalty
8715 INTEGER,
INTENT(IN) :: special_case
8717 CHARACTER(len=*),
PARAMETER :: routinen =
'main_var_to_xalmos_and_loss_func'
8719 INTEGER :: handle, ispin, nspins
8720 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: nocc
8721 REAL(kind=
dp) :: det1, energy_ispin, penalty_amplitude, &
8724 CALL timeset(routinen, handle)
8727 penalty_out = 0.0_dp
8729 nspins =
SIZE(m_main_var_in)
8730 IF (nspins == 1)
THEN
8731 spin_factor = 2.0_dp
8733 spin_factor = 1.0_dp
8736 penalty_amplitude = 0.0_dp
8738 ALLOCATE (nocc(nspins))
8739 DO ispin = 1, nspins
8741 nfullrows_total=nocc(ispin))
8744 DO ispin = 1, nspins
8747 CALL compute_xalmos_from_main_var( &
8748 m_var_in=m_main_var_in(ispin), &
8749 m_t_out=m_t_out(ispin), &
8750 m_quench_t=m_quench_t_in(ispin), &
8751 m_t0=almo_scf_env%matrix_t_blk(ispin), &
8752 m_oo_template=almo_scf_env%matrix_sigma_inv(ispin), &
8753 m_stsiginv0=m_stsiginv0_in(ispin), &
8754 m_s=almo_scf_env%matrix_s(1), &
8755 m_sig_sqrti_ii_out=m_sig_sqrti_ii_out(ispin), &
8756 domain_r_down=domain_r_down_in(:, ispin), &
8757 domain_s_inv=almo_scf_env%domain_s_inv(:, ispin), &
8758 domain_map=almo_scf_env%domain_map(ispin), &
8759 cpu_of_domain=almo_scf_env%cpu_of_domain, &
8760 assume_t0_q0x=assume_t0_q0x, &
8761 just_started=just_started, &
8762 optimize_theta=optimize_theta, &
8763 normalize_orbitals=normalize_orbitals, &
8764 envelope_amplitude=almo_scf_env%envelope_amplitude, &
8765 eps_filter=almo_scf_env%eps_filter, &
8766 special_case=special_case, &
8767 nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin), &
8768 order_lanczos=almo_scf_env%order_lanczos, &
8769 eps_lanczos=almo_scf_env%eps_lanczos, &
8770 max_iter_lanczos=almo_scf_env%max_iter_lanczos)
8775 p=almo_scf_env%matrix_p(ispin), &
8776 eps_filter=almo_scf_env%eps_filter, &
8777 orthog_orbs=.false., &
8778 nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin), &
8779 s=almo_scf_env%matrix_s(1), &
8780 sigma=almo_scf_env%matrix_sigma(ispin), &
8781 sigma_inv=almo_scf_env%matrix_sigma_inv(ispin), &
8782 use_guess=.false., &
8783 algorithm=almo_scf_env%sigma_inv_algorithm, &
8784 inv_eps_factor=almo_scf_env%matrix_iter_eps_error_factor, &
8785 inverse_accelerator=almo_scf_env%order_lanczos, &
8786 eps_lanczos=almo_scf_env%eps_lanczos, &
8787 max_iter_lanczos=almo_scf_env%max_iter_lanczos, &
8788 para_env=almo_scf_env%para_env, &
8789 blacs_env=almo_scf_env%blacs_env)
8798 IF (perturbation_only)
THEN
8800 IF (just_started)
THEN
8801 DO ispin = 1, nspins
8802 CALL dbcsr_copy(almo_scf_env%matrix_ks(ispin), &
8803 almo_scf_env%matrix_ks_0deloc(ispin))
8809 almo_scf_env%matrix_p, &
8810 almo_scf_env%matrix_ks, &
8812 almo_scf_env%eps_filter, &
8813 almo_scf_env%mat_distr_aos)
8816 penalty_out = 0.0_dp
8817 DO ispin = 1, nspins
8819 CALL compute_frequently_used_matrices( &
8820 filter_eps=almo_scf_env%eps_filter, &
8821 m_t_in=m_t_out(ispin), &
8822 m_siginv_in=almo_scf_env%matrix_sigma_inv(ispin), &
8823 m_s_in=almo_scf_env%matrix_s(1), &
8824 m_f_in=almo_scf_env%matrix_ks(ispin), &
8825 m_ftsiginv_out=m_ftsiginv_out(ispin), &
8826 m_siginvtftsiginv_out=m_siginvtftsiginv_out(ispin), &
8827 m_st_out=m_st_out(ispin))
8829 IF (perturbation_only)
THEN
8831 IF (ispin == 1) energy_out = 0.0_dp
8832 CALL dbcsr_dot(m_t_out(ispin), m_ftsiginv_out(ispin), energy_ispin)
8833 energy_out = energy_out + energy_ispin*spin_factor
8836 IF (do_penalty)
THEN
8838 CALL determinant(almo_scf_env%matrix_sigma(ispin), det1, &
8839 almo_scf_env%eps_filter)
8840 penalty_out = penalty_out - &
8841 penalty_amplitude*spin_factor*nocc(ispin)*log(det1)
8849 CALL timestop(handle)
8851 END SUBROUTINE main_var_to_xalmos_and_loss_func
8868 SUBROUTINE step_size_to_border(step_size_out, metric_in, position_in, &
8869 direction_in, trust_radius_in, quench_t_in, eps_filter_in)
8871 REAL(kind=
dp),
INTENT(INOUT) :: step_size_out
8872 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(IN) :: metric_in, position_in, direction_in
8873 REAL(kind=
dp),
INTENT(IN) :: trust_radius_in
8874 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(IN) :: quench_t_in
8875 REAL(kind=
dp),
INTENT(IN) :: eps_filter_in
8877 INTEGER :: isol, ispin, nsolutions, &
8878 nsolutions_found, nspins
8879 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: nocc
8880 REAL(kind=
dp) :: discrim_sign, discriminant, solution, &
8881 spin_factor, temp_real
8882 REAL(kind=
dp),
DIMENSION(3) :: coef
8883 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:) :: m_temp_no
8885 step_size_out = 0.0_dp
8887 nspins =
SIZE(position_in)
8888 IF (nspins == 1)
THEN
8889 spin_factor = 2.0_dp
8891 spin_factor = 1.0_dp
8894 ALLOCATE (nocc(nspins))
8895 ALLOCATE (m_temp_no(nspins))
8898 DO ispin = 1, nspins
8901 template=direction_in(ispin))
8904 nfullcols_total=nocc(ispin))
8906 CALL dbcsr_copy(m_temp_no(ispin), quench_t_in(ispin))
8909 position_in(ispin), &
8910 0.0_dp, m_temp_no(ispin), &
8911 retain_sparsity=.true.)
8913 CALL dbcsr_dot(position_in(ispin), m_temp_no(ispin), temp_real)
8914 coef(3) = coef(3) + temp_real/nocc(ispin)
8915 CALL dbcsr_dot(direction_in(ispin), m_temp_no(ispin), temp_real)
8916 coef(2) = coef(2) + 2.0_dp*temp_real/nocc(ispin)
8917 CALL dbcsr_copy(m_temp_no(ispin), quench_t_in(ispin))
8920 direction_in(ispin), &
8921 0.0_dp, m_temp_no(ispin), &
8922 retain_sparsity=.true.)
8924 CALL dbcsr_dot(direction_in(ispin), m_temp_no(ispin), temp_real)
8925 coef(1) = coef(1) + temp_real/nocc(ispin)
8932 DEALLOCATE (m_temp_no)
8934 coef(:) = coef(:)*spin_factor
8935 coef(3) = coef(3) - trust_radius_in*trust_radius_in
8938 discriminant = coef(2)*coef(2) - 4.0_dp*coef(1)*coef(3)
8939 IF (discriminant > tiny(discriminant))
THEN
8941 ELSE IF (discriminant < 0.0_dp)
THEN
8943 cpabort(
"Step to border: no solutions")
8948 discrim_sign = 1.0_dp
8949 nsolutions_found = 0
8950 DO isol = 1, nsolutions
8951 solution = (-coef(2) + discrim_sign*sqrt(discriminant))/(2.0_dp*coef(1))
8952 IF (solution > 0.0_dp)
THEN
8953 nsolutions_found = nsolutions_found + 1
8954 step_size_out = solution
8956 discrim_sign = -discrim_sign
8959 IF (nsolutions_found == 0)
THEN
8960 cpabort(
"Step to border: no positive solutions")
8961 ELSE IF (nsolutions_found == 2)
THEN
8962 cpabort(
"Two positive border steps possible!")
8965 END SUBROUTINE step_size_to_border
8978 SUBROUTINE contravariant_matrix_norm(norm_out, matrix_in, metric_in, &
8979 quench_t_in, eps_filter_in)
8981 REAL(kind=
dp),
INTENT(OUT) :: norm_out
8982 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(IN) :: matrix_in, metric_in, quench_t_in
8983 REAL(kind=
dp),
INTENT(IN) :: eps_filter_in
8985 INTEGER :: ispin, nspins
8986 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: nocc
8987 REAL(kind=
dp) :: my_norm, spin_factor, temp_real
8988 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:) :: m_temp_no
8993 nspins =
SIZE(matrix_in)
8994 IF (nspins == 1)
THEN
8995 spin_factor = 2.0_dp
8997 spin_factor = 1.0_dp
9000 ALLOCATE (nocc(nspins))
9001 ALLOCATE (m_temp_no(nspins))
9004 DO ispin = 1, nspins
9006 CALL dbcsr_create(m_temp_no(ispin), template=matrix_in(ispin))
9009 nfullcols_total=nocc(ispin))
9011 CALL dbcsr_copy(m_temp_no(ispin), quench_t_in(ispin))
9015 0.0_dp, m_temp_no(ispin), &
9016 retain_sparsity=.true.)
9018 CALL dbcsr_dot(matrix_in(ispin), m_temp_no(ispin), temp_real)
9020 my_norm = my_norm + temp_real/nocc(ispin)
9027 DEALLOCATE (m_temp_no)
9029 my_norm = my_norm*spin_factor
9030 norm_out = sqrt(my_norm)
9032 END SUBROUTINE contravariant_matrix_norm
9051 SUBROUTINE predicted_reduction(reduction_out, grad_in, step_in, hess_in, &
9052 hess_submatrix_in, quench_t_in, special_case, eps_filter, domain_map, &
9056 REAL(kind=
dp),
INTENT(INOUT) :: reduction_out
9057 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(INOUT) :: grad_in, step_in, hess_in
9059 INTENT(IN) :: hess_submatrix_in
9060 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(IN) :: quench_t_in
9061 INTEGER,
INTENT(IN) :: special_case
9062 REAL(kind=
dp),
INTENT(IN) :: eps_filter
9064 INTEGER,
DIMENSION(:),
INTENT(IN) :: cpu_of_domain
9066 INTEGER :: ispin, nspins
9067 REAL(kind=
dp) :: my_reduction, spin_factor, temp_real
9068 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:) :: m_temp_no
9070 reduction_out = 0.0_dp
9072 nspins =
SIZE(grad_in)
9073 IF (nspins == 1)
THEN
9074 spin_factor = 2.0_dp
9076 spin_factor = 1.0_dp
9079 ALLOCATE (m_temp_no(nspins))
9081 my_reduction = 0.0_dp
9082 DO ispin = 1, nspins
9084 CALL dbcsr_create(m_temp_no(ispin), template=grad_in(ispin))
9086 CALL dbcsr_dot(step_in(ispin), grad_in(ispin), temp_real)
9087 my_reduction = my_reduction + temp_real
9096 0.0_dp, m_temp_no(ispin), &
9097 filter_eps=eps_filter)
9102 matrix_in=step_in(ispin), &
9103 matrix_out=m_temp_no(ispin), &
9104 operator1=hess_submatrix_in(:, ispin), &
9105 dpattern=quench_t_in(ispin), &
9106 map=domain_map(ispin), &
9107 node_of_domain=cpu_of_domain, &
9109 filter_eps=eps_filter)
9114 CALL dbcsr_dot(step_in(ispin), m_temp_no(ispin), temp_real)
9115 my_reduction = my_reduction + 0.5_dp*temp_real
9122 my_reduction = spin_factor*my_reduction
9124 reduction_out = my_reduction
9126 DEALLOCATE (m_temp_no)
9128 END SUBROUTINE predicted_reduction
9145 SUBROUTINE fixed_r_report(unit_nr, iter_type, iteration, step_size, &
9146 border_reached, curvature, grad_norm_ratio, predicted_reduction, time)
9148 INTEGER,
INTENT(IN) :: unit_nr, iter_type, iteration
9149 REAL(kind=
dp),
INTENT(IN) :: step_size
9150 LOGICAL,
INTENT(IN) :: border_reached
9151 REAL(kind=
dp),
INTENT(IN) :: curvature
9152 REAL(kind=
dp),
INTENT(IN),
OPTIONAL :: grad_norm_ratio, predicted_reduction
9153 REAL(kind=
dp),
INTENT(IN) :: time
9155 CHARACTER(LEN=20) :: iter_type_str
9156 REAL(kind=
dp) :: loss_or_grad_change
9158 loss_or_grad_change = 0.0_dp
9159 IF (
PRESENT(grad_norm_ratio))
THEN
9160 loss_or_grad_change = grad_norm_ratio
9161 ELSE IF (
PRESENT(predicted_reduction))
THEN
9162 loss_or_grad_change = predicted_reduction
9164 cpabort(
"one argument is missing")
9167 SELECT CASE (iter_type)
9169 iter_type_str = trim(
"Ignored")
9171 iter_type_str = trim(
"PCG")
9173 iter_type_str = trim(
"Neg. curvatr.")
9175 iter_type_str = trim(
"Step too long")
9177 iter_type_str = trim(
"Grad. reduced")
9179 iter_type_str = trim(
"Cauchy point")
9181 iter_type_str = trim(
"Full dogleg")
9183 iter_type_str = trim(
"Part. dogleg")
9185 cpabort(
"unknown report type")
9188 IF (unit_nr > 0)
THEN
9190 SELECT CASE (iter_type)
9194 WRITE (unit_nr,
'(T4,A15,A6,A10,A10,A7,A20,A8)') &
9200 "Grad/o.f. reduc", &
9205 WRITE (unit_nr,
'(T4,A15,I6,F10.5,F10.5,L7,F20.10,F8.2)') &
9208 curvature, step_size, border_reached, &
9209 loss_or_grad_change, &
9215 SELECT CASE (iter_type)
9216 CASE (2, 3, 4, 5, 6, 7)
9224 END SUBROUTINE fixed_r_report
9243 SUBROUTINE trust_r_report(unit_nr, iter_type, iteration, radius, &
9244 loss, delta_loss, grad_norm, predicted_reduction, rho, new, time)
9246 INTEGER,
INTENT(IN) :: unit_nr, iter_type, iteration
9247 REAL(kind=
dp),
INTENT(IN) :: radius, loss, delta_loss, grad_norm, &
9248 predicted_reduction, rho
9249 LOGICAL,
INTENT(IN) :: new
9250 REAL(kind=
dp),
INTENT(IN) :: time
9252 CHARACTER(LEN=20) :: iter_status, iter_type_str
9254 SELECT CASE (iter_type)
9256 iter_type_str = trim(
"Iter")
9257 iter_status = trim(
"Stat")
9259 iter_type_str = trim(
"TR INI")
9261 iter_status =
" New"
9263 iter_status =
" Redo"
9266 iter_type_str = trim(
"TR FIN")
9268 iter_status =
" Acc"
9270 iter_status =
" Rej"
9273 cpabort(
"unknown report type")
9276 IF (unit_nr > 0)
THEN
9278 SELECT CASE (iter_type)
9281 WRITE (unit_nr,
'(T2,A6,A5,A6,A22,A10,T67,A7,A6)') &
9285 "Objective Function", &
9289 WRITE (unit_nr,
'(T41,A10,A10,A6)') &
9293 "Change",
"Expct.",
"Rho"
9299 WRITE (unit_nr,
'(T2,A6,A5,I6,F22.10,ES10.2,T67,ES7.0,F6.1)') &
9310 WRITE (unit_nr,
'(T2,A6,A5,I6,F22.10,ES10.2,ES10.2,F6.1,ES7.0,F6.1)') &
9315 delta_loss, predicted_reduction, rho, &
9322 END SUBROUTINE trust_r_report
9330 SUBROUTINE energy_lowering_report(unit_nr, ref_energy, energy_lowering)
9332 INTEGER,
INTENT(IN) :: unit_nr
9333 REAL(kind=
dp),
INTENT(IN) :: ref_energy, energy_lowering
9336 IF (unit_nr > 0)
THEN
9338 WRITE (unit_nr,
'(T2,A35,F25.10)')
"ENERGY OF BLOCK-DIAGONAL ALMOs:", &
9340 WRITE (unit_nr,
'(T2,A35,F25.10)')
"ENERGY LOWERING:", &
9342 WRITE (unit_nr,
'(T2,A35,F25.10)')
"CORRECTED ENERGY:", &
9343 ref_energy + energy_lowering
9347 END SUBROUTINE energy_lowering_report
9359 SUBROUTINE wrap_up_xalmo_scf(qs_env, almo_scf_env, perturbation_in, &
9360 m_xalmo_in, m_quench_in, energy_inout)
9364 LOGICAL,
INTENT(IN) :: perturbation_in
9365 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(IN) :: m_xalmo_in, m_quench_in
9366 REAL(kind=
dp),
INTENT(INOUT) :: energy_inout
9368 CHARACTER(len=*),
PARAMETER :: routinen =
'wrap_up_xalmo_scf'
9370 INTEGER :: eda_unit, handle, ispin, nspins, unit_nr
9372 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:) :: m_temp_no1, m_temp_no2
9375 CALL timeset(routinen, handle)
9379 IF (logger%para_env%is_source())
THEN
9385 nspins = almo_scf_env%nspins
9389 IF (perturbation_in)
THEN
9391 ALLOCATE (m_temp_no1(nspins))
9392 ALLOCATE (m_temp_no2(nspins))
9394 DO ispin = 1, nspins
9395 CALL dbcsr_create(m_temp_no1(ispin), template=m_xalmo_in(ispin))
9396 CALL dbcsr_create(m_temp_no2(ispin), template=m_xalmo_in(ispin))
9401 almo_scf_env%mat_distr_aos)
9406 CALL xalmo_analysis( &
9407 detailed_analysis=almo_scf_env%almo_analysis%do_analysis, &
9408 eps_filter=almo_scf_env%eps_filter, &
9409 m_t_in=m_xalmo_in, &
9410 m_t0_in=almo_scf_env%matrix_t_blk, &
9411 m_siginv_in=almo_scf_env%matrix_sigma_inv, &
9412 m_siginv0_in=almo_scf_env%matrix_sigma_inv_0deloc, &
9413 m_s_in=almo_scf_env%matrix_s, &
9414 m_ks0_in=almo_scf_env%matrix_ks_0deloc, &
9415 m_quench_t_in=m_quench_in, &
9416 energy_out=energy_inout, &
9417 m_eda_out=m_temp_no1, &
9418 m_cta_out=m_temp_no2 &
9421 IF (almo_scf_env%almo_analysis%do_analysis)
THEN
9423 DO ispin = 1, nspins
9426 IF (unit_nr > 0)
THEN
9427 WRITE (unit_nr,
'(T2,A)')
"DECOMPOSITION OF THE DELOCALIZATION ENERGY"
9434 "ALMO_EDA_CT", extension=
".dat", local=.true.)
9435 CALL print_block_sum(m_temp_no1(ispin), eda_unit)
9437 "ALMO_EDA_CT", local=.true.)
9440 IF (unit_nr > 0)
THEN
9441 WRITE (unit_nr,
'(T2,A)')
"DECOMPOSITION OF CHARGE TRANSFER TERMS"
9445 "ALMO_CTA", extension=
".dat", local=.true.)
9446 CALL print_block_sum(m_temp_no2(ispin), eda_unit)
9448 "ALMO_CTA", local=.true.)
9454 CALL energy_lowering_report( &
9456 ref_energy=almo_scf_env%almo_scf_energy, &
9457 energy_lowering=energy_inout)
9459 energy=almo_scf_env%almo_scf_energy, &
9460 energy_singles_corr=energy_inout)
9462 DO ispin = 1, nspins
9467 DEALLOCATE (m_temp_no1)
9468 DEALLOCATE (m_temp_no2)
9473 energy=energy_inout)
9477 CALL timestop(handle)
9479 END SUBROUTINE wrap_up_xalmo_scf
9487 SUBROUTINE tanh_of_elements(matrix, alpha)
9489 REAL(kind=
dp),
INTENT(IN) :: alpha
9491 CHARACTER(len=*),
PARAMETER :: routinen =
'tanh_of_elements'
9494 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: block
9497 CALL timeset(routinen, handle)
9501 block = tanh(alpha*block)
9504 CALL timestop(handle)
9506 END SUBROUTINE tanh_of_elements
9514 SUBROUTINE dtanh_of_elements(matrix, alpha)
9516 REAL(kind=
dp),
INTENT(IN) :: alpha
9518 CHARACTER(len=*),
PARAMETER :: routinen =
'dtanh_of_elements'
9521 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: block
9524 CALL timeset(routinen, handle)
9528 block = alpha*(1.0_dp - tanh(block)**2)
9531 CALL timestop(handle)
9533 END SUBROUTINE dtanh_of_elements
9540 SUBROUTINE inverse_of_elements(matrix)
9543 CHARACTER(len=*),
PARAMETER :: routinen =
'inverse_of_elements'
9546 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: block
9549 CALL timeset(routinen, handle)
9553 block = 1.0_dp/block
9556 CALL timestop(handle)
9558 END SUBROUTINE inverse_of_elements
9565 SUBROUTINE print_block_sum(matrix, unit_nr)
9567 INTEGER,
INTENT(IN) :: unit_nr
9569 CHARACTER(len=*),
PARAMETER :: routinen =
'print_block_sum'
9571 INTEGER :: col, handle, row
9572 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: block
9575 CALL timeset(routinen, handle)
9577 IF (unit_nr > 0)
THEN
9581 WRITE (unit_nr,
'(I6,I6,ES18.9)') row, col, sum(block)
9586 CALL timestop(handle)
9587 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 preconditionersimple 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)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)
...
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_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset, transposed)
...
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.
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, ncgf)
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, mimic, 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, sab_cneo, 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, matrix_vhxc, 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, xcint_weights, 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, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_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, do_rixs, tb_tblite)
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 reciprocal-space operators used by Berry localization. The operators contain the contrac...
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.