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 CALL dbcsr_set(op_sm_set_qs(reim, idim0)%matrix, 0.0_dp)
922 NULLIFY (op_sm_set_almo(reim, idim0)%matrix)
923 ALLOCATE (op_sm_set_almo(reim, idim0)%matrix)
924 CALL dbcsr_copy(op_sm_set_almo(reim, idim0)%matrix, almo_scf_env%matrix_s(1), &
926 CALL dbcsr_set(op_sm_set_almo(reim, idim0)%matrix, 0.0_dp)
938 m_t_in=m_t_in_local, &
939 m_t0=almo_scf_env%matrix_t_blk, &
940 m_quench_t=quench_t, &
941 m_overlap=almo_scf_env%matrix_s(1), &
942 m_sigma_tmpl=almo_scf_env%matrix_sigma_inv, &
944 xalmo_history=almo_scf_env%xalmo_history, &
945 assume_t0_q0x=assume_t0_q0x, &
946 optimize_theta=optimize_theta, &
947 envelope_amplitude=almo_scf_env%envelope_amplitude, &
948 eps_filter=almo_scf_env%eps_filter, &
949 order_lanczos=almo_scf_env%order_lanczos, &
950 eps_lanczos=almo_scf_env%eps_lanczos, &
951 max_iter_lanczos=almo_scf_env%max_iter_lanczos, &
952 nocc_of_domain=almo_scf_env%nocc_of_domain)
954 ndomains = almo_scf_env%ndomains
955 ALLOCATE (domain_r_down(ndomains, nspins))
957 ALLOCATE (bad_modes_projector_down(ndomains, nspins))
960 ALLOCATE (prec_vv(nspins))
961 ALLOCATE (siginvtftsiginv(nspins))
962 ALLOCATE (stsiginv_0(nspins))
963 ALLOCATE (ftsiginv(nspins))
964 ALLOCATE (st(nspins))
965 ALLOCATE (prev_grad(nspins))
966 ALLOCATE (grad(nspins))
967 ALLOCATE (prev_step(nspins))
968 ALLOCATE (step(nspins))
969 ALLOCATE (prev_minus_prec_grad(nspins))
970 ALLOCATE (m_sig_sqrti_ii(nspins))
971 ALLOCATE (tempnocc(nspins))
972 ALLOCATE (tempnocc_1(nspins))
973 ALLOCATE (tempoccocc(nspins))
978 template=almo_scf_env%matrix_ks(ispin), &
979 matrix_type=dbcsr_type_no_symmetry)
981 template=almo_scf_env%matrix_sigma(ispin), &
982 matrix_type=dbcsr_type_no_symmetry)
984 template=matrix_t_out(ispin), &
985 matrix_type=dbcsr_type_no_symmetry)
987 template=matrix_t_out(ispin), &
988 matrix_type=dbcsr_type_no_symmetry)
990 template=matrix_t_out(ispin), &
991 matrix_type=dbcsr_type_no_symmetry)
993 template=matrix_t_out(ispin), &
994 matrix_type=dbcsr_type_no_symmetry)
996 template=matrix_t_out(ispin), &
997 matrix_type=dbcsr_type_no_symmetry)
999 template=matrix_t_out(ispin), &
1000 matrix_type=dbcsr_type_no_symmetry)
1002 template=matrix_t_out(ispin), &
1003 matrix_type=dbcsr_type_no_symmetry)
1005 template=matrix_t_out(ispin), &
1006 matrix_type=dbcsr_type_no_symmetry)
1008 template=almo_scf_env%matrix_sigma_inv(ispin), &
1009 matrix_type=dbcsr_type_no_symmetry)
1011 template=matrix_t_out(ispin), &
1012 matrix_type=dbcsr_type_no_symmetry)
1014 template=matrix_t_out(ispin), &
1015 matrix_type=dbcsr_type_no_symmetry)
1017 template=almo_scf_env%matrix_sigma_inv(ispin), &
1018 matrix_type=dbcsr_type_no_symmetry)
1021 CALL dbcsr_set(prev_step(ispin), 0.0_dp)
1024 nfullrows_total=nocc(ispin))
1031 matrix_s=almo_scf_env%matrix_s(1), &
1032 subm_s_inv=almo_scf_env%domain_s_inv(:, ispin), &
1033 dpattern=quench_t(ispin), &
1034 map=almo_scf_env%domain_map(ispin), &
1035 node_of_domain=almo_scf_env%cpu_of_domain)
1038 matrix_s=almo_scf_env%matrix_s(1), &
1039 subm_s_sqrt=almo_scf_env%domain_s_sqrt(:, ispin), &
1040 subm_s_sqrt_inv=almo_scf_env%domain_s_sqrt_inv(:, ispin), &
1041 dpattern=almo_scf_env%quench_t(ispin), &
1042 map=almo_scf_env%domain_map(ispin), &
1043 node_of_domain=almo_scf_env%cpu_of_domain)
1047 IF (assume_t0_q0x)
THEN
1052 almo_scf_env%matrix_s(1), &
1053 almo_scf_env%matrix_t_blk(ispin), &
1054 0.0_dp, st(ispin), &
1055 filter_eps=almo_scf_env%eps_filter)
1058 almo_scf_env%matrix_sigma_inv_0deloc(ispin), &
1059 0.0_dp, stsiginv_0(ispin), &
1060 filter_eps=almo_scf_env%eps_filter)
1066 matrix_t=almo_scf_env%matrix_t_blk(ispin), &
1067 matrix_sigma_inv=almo_scf_env%matrix_sigma_inv(ispin), &
1068 matrix_s=almo_scf_env%matrix_s(1), &
1069 subm_r_down=domain_r_down(:, ispin), &
1070 dpattern=quench_t(ispin), &
1071 map=almo_scf_env%domain_map(ispin), &
1072 node_of_domain=almo_scf_env%cpu_of_domain, &
1073 filter_eps=almo_scf_env%eps_filter)
1079 IF (penalty_occ_local)
THEN
1083 almo_scf_env%matrix_s(1), &
1084 matrix_t_in(ispin), &
1085 0.0_dp, tempnocc(ispin), &
1086 filter_eps=almo_scf_env%eps_filter)
1089 almo_scf_env%matrix_sigma_inv(ispin), &
1090 0.0_dp, tempnocc_1(ispin), &
1091 filter_eps=almo_scf_env%eps_filter)
1093 DO idim0 = 1,
SIZE(op_sm_set_qs, 2)
1094 DO reim = 1,
SIZE(op_sm_set_qs, 1)
1097 op_sm_set_almo(reim, idim0)%matrix, almo_scf_env%mat_distr_aos)
1100 op_sm_set_almo(reim, idim0)%matrix, &
1101 matrix_t_in(ispin), &
1102 0.0_dp, tempnocc(ispin), &
1103 filter_eps=almo_scf_env%eps_filter)
1106 matrix_t_in(ispin), &
1108 0.0_dp, tempoccocc(ispin), &
1109 filter_eps=almo_scf_env%eps_filter)
1112 tempnocc_1(ispin), &
1113 tempoccocc(ispin), &
1114 0.0_dp, tempnocc(ispin), &
1115 filter_eps=almo_scf_env%eps_filter)
1119 tempnocc_1(ispin), &
1120 0.0_dp, op_sm_set_almo(reim, idim0)%matrix, &
1121 filter_eps=almo_scf_env%eps_filter)
1131 outer_max_iter = optimizer%max_iter_outer_loop
1132 outer_prepare_to_exit = .false.
1135 grad_norm_frob = 0.0_dp
1141 max_iter = optimizer%max_iter
1142 prepare_to_exit = .false.
1143 line_search = .false.
1147 line_search_iteration = 0
1150 energy_diff = 0.0_dp
1151 localization_obj_function = 0.0_dp
1152 line_search_error = 0.0_dp
1158 just_started = (iteration == 0) .AND. (outer_iteration == 0)
1160 CALL main_var_to_xalmos_and_loss_func( &
1161 almo_scf_env=almo_scf_env, &
1163 m_main_var_in=m_theta, &
1164 m_t_out=matrix_t_out, &
1165 m_sig_sqrti_ii_out=m_sig_sqrti_ii, &
1166 energy_out=energy_new, &
1167 penalty_out=penalty_func_new, &
1168 m_ftsiginv_out=ftsiginv, &
1169 m_siginvtftsiginv_out=siginvtftsiginv, &
1171 m_stsiginv0_in=stsiginv_0, &
1172 m_quench_t_in=quench_t, &
1173 domain_r_down_in=domain_r_down, &
1174 assume_t0_q0x=assume_t0_q0x, &
1175 just_started=just_started, &
1176 optimize_theta=optimize_theta, &
1177 normalize_orbitals=normalize_orbitals, &
1178 perturbation_only=perturbation_only, &
1179 do_penalty=penalty_occ_vol, &
1180 special_case=my_special_case)
1181 IF (penalty_occ_vol)
THEN
1183 energy_new = energy_new + penalty_func_new
1185 DO ispin = 1, nspins
1186 IF (penalty_occ_vol)
THEN
1187 penalty_occ_vol_g_prefactor(ispin) = &
1188 -2.0_dp*penalty_amplitude*spin_factor*nocc(ispin)
1189 penalty_occ_vol_h_prefactor(ispin) = 0.0_dp
1193 localization_obj_function = 0.0_dp
1195 IF (penalty_occ_local)
THEN
1196 DO ispin = 1, nspins
1199 localization_obj_function = 0.0_dp
1200 CALL dbcsr_get_info(almo_scf_env%matrix_sigma_inv(ispin), nfullrows_total=nmo)
1202 ALLOCATE (reim_diag(nmo))
1206 DO idim0 = 1,
SIZE(op_sm_set_qs, 2)
1210 DO reim = 1,
SIZE(op_sm_set_qs, 1)
1213 op_sm_set_almo(reim, idim0)%matrix, &
1214 matrix_t_out(ispin), &
1215 0.0_dp, tempnocc(ispin), &
1216 filter_eps=almo_scf_env%eps_filter)
1219 matrix_t_out(ispin), &
1221 0.0_dp, tempoccocc(ispin), &
1222 filter_eps=almo_scf_env%eps_filter)
1226 CALL group%sum(reim_diag)
1227 z2(:) = z2(:) + reim_diag(:)*reim_diag(:)
1234 fval = -weights(idim0)*log(abs(z2(ielem)))
1236 fval = weights(idim0) - weights(idim0)*abs(z2(ielem))
1238 fval = weights(idim0) - weights(idim0)*sqrt(abs(z2(ielem)))
1240 localization_obj_function = localization_obj_function + fval
1246 DEALLOCATE (reim_diag)
1248 energy_new = energy_new + localiz_coeff*localization_obj_function
1253 DO ispin = 1, nspins
1255 IF (just_started .AND. almo_mathematica)
THEN
1256 cpwarn_if(ispin > 1,
"Mathematica files will be overwritten")
1257 CALL print_mathematica_matrix(almo_scf_env%matrix_s(1),
"matrixS.dat")
1258 CALL print_mathematica_matrix(almo_scf_env%matrix_ks(ispin),
"matrixF.dat")
1259 CALL print_mathematica_matrix(matrix_t_out(ispin),
"matrixT.dat")
1260 CALL print_mathematica_matrix(quench_t(ispin),
"matrixQ.dat")
1266 IF (line_search_iteration == 0 .AND. iteration /= 0)
THEN
1267 CALL dbcsr_copy(prev_grad(ispin), grad(ispin))
1273 skip_grad = (iteration > 0 .AND. &
1274 fixed_line_search_niter /= 0 .AND. &
1275 line_search_iteration /= fixed_line_search_niter)
1277 IF (.NOT. skip_grad)
THEN
1279 DO ispin = 1, nspins
1281 CALL compute_gradient( &
1282 m_grad_out=grad(ispin), &
1283 m_ks=almo_scf_env%matrix_ks(ispin), &
1284 m_s=almo_scf_env%matrix_s(1), &
1285 m_t=matrix_t_out(ispin), &
1286 m_t0=almo_scf_env%matrix_t_blk(ispin), &
1287 m_siginv=almo_scf_env%matrix_sigma_inv(ispin), &
1288 m_quench_t=quench_t(ispin), &
1289 m_ftsiginv=ftsiginv(ispin), &
1290 m_siginvtftsiginv=siginvtftsiginv(ispin), &
1292 m_stsiginv0=stsiginv_0(ispin), &
1293 m_theta=m_theta(ispin), &
1294 m_sig_sqrti_ii=m_sig_sqrti_ii(ispin), &
1295 domain_s_inv=almo_scf_env%domain_s_inv(:, ispin), &
1296 domain_r_down=domain_r_down(:, ispin), &
1297 cpu_of_domain=almo_scf_env%cpu_of_domain, &
1298 domain_map=almo_scf_env%domain_map(ispin), &
1299 assume_t0_q0x=assume_t0_q0x, &
1300 optimize_theta=optimize_theta, &
1301 normalize_orbitals=normalize_orbitals, &
1302 penalty_occ_vol=penalty_occ_vol, &
1303 penalty_occ_vol_prefactor=penalty_occ_vol_g_prefactor(ispin), &
1304 envelope_amplitude=almo_scf_env%envelope_amplitude, &
1305 eps_filter=almo_scf_env%eps_filter, &
1306 spin_factor=spin_factor, &
1307 special_case=my_special_case, &
1308 penalty_occ_local=penalty_occ_local, &
1309 op_sm_set=op_sm_set_almo, &
1311 energy_coeff=energy_coeff, &
1312 localiz_coeff=localiz_coeff)
1321 IF (blissful_neglect)
THEN
1322 DO ispin = 1, nspins
1325 IF (iteration == 0)
THEN
1326 CALL compute_preconditioner( &
1327 domain_prec_out=almo_scf_env%domain_preconditioner(:, ispin), &
1328 bad_modes_projector_down_out=bad_modes_projector_down(:, ispin), &
1329 m_prec_out=prec_vv(ispin), &
1330 m_ks=almo_scf_env%matrix_ks(ispin), &
1331 m_s=almo_scf_env%matrix_s(1), &
1332 m_siginv=almo_scf_env%matrix_sigma_inv(ispin), &
1333 m_quench_t=quench_t(ispin), &
1334 m_ftsiginv=ftsiginv(ispin), &
1335 m_siginvtftsiginv=siginvtftsiginv(ispin), &
1337 para_env=almo_scf_env%para_env, &
1338 blacs_env=almo_scf_env%blacs_env, &
1339 nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin), &
1340 domain_s_inv=almo_scf_env%domain_s_inv(:, ispin), &
1341 domain_s_inv_half=almo_scf_env%domain_s_sqrt_inv(:, ispin), &
1342 domain_s_half=almo_scf_env%domain_s_sqrt(:, ispin), &
1343 domain_r_down=domain_r_down(:, ispin), &
1344 cpu_of_domain=almo_scf_env%cpu_of_domain, &
1345 domain_map=almo_scf_env%domain_map(ispin), &
1346 assume_t0_q0x=assume_t0_q0x, &
1347 penalty_occ_vol=penalty_occ_vol, &
1348 penalty_occ_vol_prefactor=penalty_occ_vol_g_prefactor(ispin), &
1349 eps_filter=almo_scf_env%eps_filter, &
1350 neg_thr=optimizer%neglect_threshold, &
1351 spin_factor=spin_factor, &
1352 skip_inversion=.false., &
1353 special_case=my_special_case)
1357 matrix_in=grad(ispin), &
1358 matrix_out=grad(ispin), &
1359 operator1=almo_scf_env%domain_s_inv(:, ispin), &
1360 operator2=bad_modes_projector_down(:, ispin), &
1361 dpattern=quench_t(ispin), &
1362 map=almo_scf_env%domain_map(ispin), &
1363 node_of_domain=almo_scf_env%cpu_of_domain, &
1365 filter_eps=almo_scf_env%eps_filter)
1372 DO ispin = 1, nspins
1375 grad_norm = maxval(grad_norm_spin)
1377 converged = (grad_norm <= optimizer%eps_error)
1378 IF (converged .OR. (iteration >= max_iter))
THEN
1379 prepare_to_exit = .true.
1382 IF (optimizer%early_stopping_on .AND. just_started)
THEN
1383 prepare_to_exit = .false.
1386 IF (grad_norm < almo_scf_env%eps_prev_guess)
THEN
1391 IF (.NOT. prepare_to_exit)
THEN
1396 IF (iteration /= 0)
THEN
1398 IF (fixed_line_search_niter == 0)
THEN
1402 IF (.NOT. line_search)
THEN
1404 line_search = .true.
1405 line_search_iteration = line_search_iteration + 1
1411 line_search_error = 0.0_dp
1415 DO ispin = 1, nspins
1417 CALL dbcsr_dot(grad(ispin), step(ispin), tempreal)
1418 line_search_error = line_search_error + tempreal
1419 CALL dbcsr_dot(grad(ispin), grad(ispin), tempreal)
1420 denom = denom + tempreal
1421 CALL dbcsr_dot(step(ispin), step(ispin), tempreal)
1422 denom2 = denom2 + tempreal
1428 line_search_error = line_search_error/sqrt(denom)/sqrt(denom2)
1430 IF (abs(line_search_error) > optimizer%lin_search_eps_error)
THEN
1431 line_search = .true.
1432 line_search_iteration = line_search_iteration + 1
1434 line_search = .false.
1435 line_search_iteration = 0
1436 IF (grad_norm < eps_skip_gradients)
THEN
1437 fixed_line_search_niter = abs(almo_scf_env%integer04)
1445 IF (.NOT. line_search)
THEN
1446 line_search = .true.
1447 line_search_iteration = line_search_iteration + 1
1449 IF (line_search_iteration == fixed_line_search_niter)
THEN
1450 line_search = .false.
1451 line_search_iteration = 0
1452 line_search_iteration = line_search_iteration + 1
1460 IF (line_search)
THEN
1461 energy_diff = 0.0_dp
1463 energy_diff = energy_new - energy_old
1464 energy_old = energy_new
1468 IF (.NOT. line_search)
THEN
1470 cg_iteration = cg_iteration + 1
1473 DO ispin = 1, nspins
1474 CALL dbcsr_copy(prev_step(ispin), step(ispin))
1478 SELECT CASE (prec_type)
1482 CALL newton_grad_to_step( &
1483 optimizer=almo_scf_env%opt_xalmo_newton_pcg_solver, &
1486 m_s=almo_scf_env%matrix_s(:), &
1487 m_ks=almo_scf_env%matrix_ks(:), &
1488 m_siginv=almo_scf_env%matrix_sigma_inv(:), &
1489 m_quench_t=quench_t(:), &
1490 m_ftsiginv=ftsiginv(:), &
1491 m_siginvtftsiginv=siginvtftsiginv(:), &
1493 m_t=matrix_t_out(:), &
1494 m_sig_sqrti_ii=m_sig_sqrti_ii(:), &
1495 domain_s_inv=almo_scf_env%domain_s_inv(:, :), &
1496 domain_r_down=domain_r_down(:, :), &
1497 domain_map=almo_scf_env%domain_map(:), &
1498 cpu_of_domain=almo_scf_env%cpu_of_domain, &
1499 nocc_of_domain=almo_scf_env%nocc_of_domain(:, :), &
1500 para_env=almo_scf_env%para_env, &
1501 blacs_env=almo_scf_env%blacs_env, &
1502 eps_filter=almo_scf_env%eps_filter, &
1503 optimize_theta=optimize_theta, &
1504 penalty_occ_vol=penalty_occ_vol, &
1505 normalize_orbitals=normalize_orbitals, &
1506 penalty_occ_vol_prefactor=penalty_occ_vol_g_prefactor(:), &
1507 penalty_occ_vol_pf2=penalty_occ_vol_h_prefactor(:), &
1508 special_case=my_special_case &
1514 IF (.NOT. blissful_neglect .AND. &
1515 ((just_started .AND. perturbation_only) .OR. &
1516 (iteration == 0 .AND. (.NOT. perturbation_only))) &
1520 DO ispin = 1, nspins
1521 CALL compute_preconditioner( &
1522 domain_prec_out=almo_scf_env%domain_preconditioner(:, ispin), &
1523 m_prec_out=prec_vv(ispin), &
1524 m_ks=almo_scf_env%matrix_ks(ispin), &
1525 m_s=almo_scf_env%matrix_s(1), &
1526 m_siginv=almo_scf_env%matrix_sigma_inv(ispin), &
1527 m_quench_t=quench_t(ispin), &
1528 m_ftsiginv=ftsiginv(ispin), &
1529 m_siginvtftsiginv=siginvtftsiginv(ispin), &
1531 para_env=almo_scf_env%para_env, &
1532 blacs_env=almo_scf_env%blacs_env, &
1533 nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin), &
1534 domain_s_inv=almo_scf_env%domain_s_inv(:, ispin), &
1535 domain_r_down=domain_r_down(:, ispin), &
1536 cpu_of_domain=almo_scf_env%cpu_of_domain, &
1537 domain_map=almo_scf_env%domain_map(ispin), &
1538 assume_t0_q0x=assume_t0_q0x, &
1539 penalty_occ_vol=penalty_occ_vol, &
1540 penalty_occ_vol_prefactor=penalty_occ_vol_g_prefactor(ispin), &
1541 eps_filter=almo_scf_env%eps_filter, &
1543 spin_factor=spin_factor, &
1544 skip_inversion=.false., &
1545 special_case=my_special_case)
1552 DO ispin = 1, nspins
1557 0.0_dp, step(ispin), &
1558 filter_eps=almo_scf_env%eps_filter)
1565 IF (optimize_theta)
THEN
1566 cpabort(
"theta is NYI")
1569 DO ispin = 1, nspins
1572 matrix_in=grad(ispin), &
1573 matrix_out=step(ispin), &
1574 operator1=almo_scf_env%domain_preconditioner(:, ispin), &
1575 dpattern=quench_t(ispin), &
1576 map=almo_scf_env%domain_map(ispin), &
1577 node_of_domain=almo_scf_env%cpu_of_domain, &
1579 filter_eps=almo_scf_env%eps_filter)
1589 DO ispin = 1, nspins
1599 IF (iteration == 0)
THEN
1600 reset_conjugator = .true.
1604 IF (.NOT. reset_conjugator)
THEN
1606 CALL compute_cg_beta( &
1608 reset_conjugator=reset_conjugator, &
1609 conjugator=optimizer%conjugator, &
1611 prev_grad=prev_grad(:), &
1613 prev_step=prev_step(:), &
1614 prev_minus_prec_grad=prev_minus_prec_grad(:) &
1619 IF (reset_conjugator)
THEN
1622 IF (unit_nr > 0 .AND. (.NOT. just_started))
THEN
1623 WRITE (unit_nr,
'(T2,A35)')
"Re-setting conjugator to zero"
1625 reset_conjugator = .false.
1630 DO ispin = 1, nspins
1632 CALL dbcsr_copy(prev_minus_prec_grad(ispin), step(ispin))
1635 CALL dbcsr_add(step(ispin), prev_step(ispin), 1.0_dp, beta)
1642 IF (.NOT. line_search)
THEN
1648 DO ispin = 1, nspins
1649 CALL dbcsr_dot(grad(ispin), step(ispin), tempreal)
1652 IF (iteration == 0)
THEN
1653 step_size = optimizer%lin_search_step_size_guess
1655 IF (next_step_size_guess <= 0.0_dp)
THEN
1656 step_size = optimizer%lin_search_step_size_guess
1659 step_size = next_step_size_guess*1.05_dp
1662 next_step_size_guess = step_size
1664 IF (fixed_line_search_niter == 0)
THEN
1667 DO ispin = 1, nspins
1668 CALL dbcsr_dot(grad(ispin), step(ispin), tempreal)
1673 appr_sec_der = (g1 - g0)/step_size
1674 step_size = -g1/appr_sec_der
1682 appr_sec_der = 2.0*((e1 - e0)/step_size - g0)/step_size
1683 g1 = appr_sec_der*step_size + g0
1684 step_size = -g1/appr_sec_der
1688 next_step_size_guess = next_step_size_guess + step_size
1692 DO ispin = 1, nspins
1693 CALL dbcsr_add(m_theta(ispin), step(ispin), 1.0_dp, step_size)
1698 IF (line_search)
THEN
1705 IF (unit_nr > 0)
THEN
1706 iter_type = trim(
"ALMO SCF "//iter_type)
1707 WRITE (unit_nr,
'(T2,A13,I6,F23.10,E14.5,F14.9,F9.2)') &
1708 iter_type, iteration, &
1709 energy_new, energy_diff, grad_norm, &
1711 IF (penalty_occ_local .OR. penalty_occ_vol)
THEN
1712 WRITE (unit_nr,
'(T2,A25,F23.10)') &
1713 "Energy component:", (energy_new - penalty_func_new - localization_obj_function)
1715 IF (penalty_occ_local)
THEN
1716 WRITE (unit_nr,
'(T2,A25,F23.10)') &
1717 "Localization component:", localization_obj_function
1719 IF (penalty_occ_vol)
THEN
1720 WRITE (unit_nr,
'(T2,A25,F23.10)') &
1721 "Penalty component:", penalty_func_new
1726 IF (penalty_occ_vol)
THEN
1727 almo_scf_env%almo_scf_energy = energy_new - penalty_func_new - localization_obj_function
1729 almo_scf_env%almo_scf_energy = energy_new - localization_obj_function
1735 iteration = iteration + 1
1736 IF (prepare_to_exit)
EXIT
1740 IF (converged .OR. (outer_iteration >= outer_max_iter))
THEN
1741 outer_prepare_to_exit = .true.
1744 outer_iteration = outer_iteration + 1
1745 IF (outer_prepare_to_exit)
EXIT
1749 DO ispin = 1, nspins
1750 IF (converged .AND. almo_mathematica)
THEN
1751 cpwarn_if(ispin > 1,
"Mathematica files will be overwritten")
1752 CALL print_mathematica_matrix(matrix_t_out(ispin),
"matrixTf.dat")
1759 CALL wrap_up_xalmo_scf( &
1761 almo_scf_env=almo_scf_env, &
1762 perturbation_in=perturbation_only, &
1763 m_xalmo_in=matrix_t_out, &
1764 m_quench_in=quench_t, &
1765 energy_inout=energy_new)
1769 DO ispin = 1, nspins
1790 DEALLOCATE (tempnocc)
1791 DEALLOCATE (tempnocc_1)
1792 DEALLOCATE (tempoccocc)
1793 DEALLOCATE (prec_vv)
1794 DEALLOCATE (siginvtftsiginv)
1795 DEALLOCATE (stsiginv_0)
1796 DEALLOCATE (ftsiginv)
1798 DEALLOCATE (prev_grad)
1800 DEALLOCATE (prev_step)
1802 DEALLOCATE (prev_minus_prec_grad)
1803 DEALLOCATE (m_sig_sqrti_ii)
1805 DEALLOCATE (domain_r_down)
1806 DEALLOCATE (bad_modes_projector_down)
1808 DEALLOCATE (penalty_occ_vol_g_prefactor)
1809 DEALLOCATE (penalty_occ_vol_h_prefactor)
1810 DEALLOCATE (grad_norm_spin)
1813 DEALLOCATE (m_theta, m_t_in_local)
1814 IF (penalty_occ_local)
THEN
1815 DO idim0 = 1, dim_op
1816 DO reim = 1,
SIZE(op_sm_set_qs, 1)
1817 DEALLOCATE (op_sm_set_qs(reim, idim0)%matrix)
1818 DEALLOCATE (op_sm_set_almo(reim, idim0)%matrix)
1821 DEALLOCATE (op_sm_set_qs)
1822 DEALLOCATE (op_sm_set_almo)
1823 DEALLOCATE (weights)
1826 IF (.NOT. converged .AND. .NOT. optimizer%early_stopping_on)
THEN
1827 cpabort(
"Optimization not converged! ")
1830 CALL timestop(handle)
1851 matrix_s, matrix_mo_in, matrix_mo_out, &
1852 template_matrix_sigma, overlap_determinant, &
1853 mat_distr_aos, virtuals, eps_filter)
1857 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:), &
1858 INTENT(INOUT) :: matrix_mo_in, matrix_mo_out
1859 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:), &
1860 INTENT(IN) :: template_matrix_sigma
1861 REAL(kind=
dp),
INTENT(INOUT) :: overlap_determinant
1862 INTEGER,
INTENT(IN) :: mat_distr_aos
1863 LOGICAL,
INTENT(IN) :: virtuals
1864 REAL(kind=
dp),
INTENT(IN) :: eps_filter
1866 CHARACTER(len=*),
PARAMETER :: routinen =
'almo_scf_construct_nlmos'
1868 CHARACTER(LEN=30) :: iter_type, print_string
1869 INTEGER :: cg_iteration, dim_op, handle, iatom, idim0, isgf, ispin, iteration, &
1870 line_search_iteration, linear_search_type, max_iter, natom, ncol, nspins, &
1871 outer_iteration, outer_max_iter, prec_type, reim, unit_nr
1872 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: first_sgf, last_sgf, nocc, nsgf
1873 LOGICAL :: converged, d_bfgs, just_started, l_bfgs, &
1874 line_search, outer_prepare_to_exit, &
1875 prepare_to_exit, reset_conjugator
1876 REAL(kind=
dp) :: appr_sec_der, beta, bfgs_rho, bfgs_sum, denom, denom2, e0, e1, g0, g0sign, &
1877 g1, g1sign, grad_norm, line_search_error, localization_obj_function, &
1878 localization_obj_function_ispin, next_step_size_guess, obj_function_ispin, objf_diff, &
1879 objf_new, objf_old, penalty_amplitude, penalty_func_ispin, penalty_func_new, spin_factor, &
1880 step_size, t1, t2, tempreal
1881 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: diagonal, grad_norm_spin, &
1882 penalty_vol_prefactor, &
1883 suggested_vol_penalty, weights
1886 TYPE(
dbcsr_p_type),
DIMENSION(:),
POINTER :: qs_matrix_s
1887 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
POINTER :: op_sm_set_almo, op_sm_set_qs
1888 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:) :: approx_inv_hessian, bfgs_s, bfgs_y, grad, &
1889 m_s0, m_sig_sqrti_ii, m_siginv, m_sigma, m_t_mo_local, m_theta, m_theta_normalized, &
1890 prev_grad, prev_m_theta, prev_minus_prec_grad, prev_step, step, tempnocc1, tempoccocc1, &
1891 tempoccocc2, tempoccocc3
1892 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:, :, :) :: m_b0
1896 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
1898 CALL timeset(routinen, handle)
1902 IF (logger%para_env%is_source())
THEN
1908 nspins =
SIZE(matrix_mo_in)
1910 IF (unit_nr > 0)
THEN
1912 IF (.NOT. virtuals)
THEN
1913 WRITE (unit_nr,
'(T2,A,A,A)') repeat(
"-", 24), &
1914 " Optimization of occupied NLMOs ", repeat(
"-", 23)
1916 WRITE (unit_nr,
'(T2,A,A,A)') repeat(
"-", 24), &
1917 " Optimization of virtual NLMOs ", repeat(
"-", 24)
1920 WRITE (unit_nr,
'(T2,A13,A6,A23,A14,A14,A9)')
"Method",
"Iter", &
1921 "Objective Function",
"Change",
"Convergence",
"Time"
1922 WRITE (unit_nr,
'(T2,A)') repeat(
"-", 79)
1925 NULLIFY (particle_set)
1928 matrix_s=qs_matrix_s, &
1930 particle_set=particle_set, &
1931 qs_kind_set=qs_kind_set)
1933 natom =
SIZE(particle_set, 1)
1934 ALLOCATE (first_sgf(natom))
1935 ALLOCATE (last_sgf(natom))
1936 ALLOCATE (nsgf(natom))
1939 first_sgf=first_sgf, last_sgf=last_sgf, nsgf=nsgf)
1943 ALLOCATE (m_theta(nspins))
1944 DO ispin = 1, nspins
1946 template=template_matrix_sigma(ispin), &
1947 matrix_type=dbcsr_type_no_symmetry)
1953 SELECT CASE (optimizer%opt_penalty%operator_type)
1956 IF (cell%orthorhombic)
THEN
1961 ALLOCATE (weights(6))
1964 ALLOCATE (op_sm_set_qs(2, dim_op))
1965 ALLOCATE (op_sm_set_almo(2, dim_op))
1967 ALLOCATE (m_b0(2, dim_op, nspins))
1968 DO idim0 = 1, dim_op
1969 DO reim = 1,
SIZE(op_sm_set_qs, 1)
1970 NULLIFY (op_sm_set_qs(reim, idim0)%matrix, op_sm_set_almo(reim, idim0)%matrix)
1971 ALLOCATE (op_sm_set_qs(reim, idim0)%matrix)
1972 ALLOCATE (op_sm_set_almo(reim, idim0)%matrix)
1973 CALL dbcsr_copy(op_sm_set_qs(reim, idim0)%matrix, qs_matrix_s(1)%matrix, &
1975 CALL dbcsr_set(op_sm_set_qs(reim, idim0)%matrix, 0.0_dp)
1976 CALL dbcsr_copy(op_sm_set_almo(reim, idim0)%matrix, matrix_s, &
1978 CALL dbcsr_set(op_sm_set_almo(reim, idim0)%matrix, 0.0_dp)
1979 DO ispin = 1, nspins
1981 template=m_theta(ispin), &
1982 matrix_type=dbcsr_type_no_symmetry)
1983 CALL dbcsr_set(m_b0(reim, idim0, ispin), 0.0_dp)
1993 ALLOCATE (weights(dim_op))
1996 ALLOCATE (m_b0(1, dim_op, nspins))
1998 DO idim0 = 1, dim_op
2000 DO ispin = 1, nspins
2002 template=m_theta(ispin), &
2003 matrix_type=dbcsr_type_no_symmetry)
2004 CALL dbcsr_set(m_b0(reim, idim0, ispin), 0.0_dp)
2012 penalty_amplitude = optimizer%opt_penalty%penalty_strength
2015 prec_type = optimizer%preconditioner
2021 IF (l_bfgs .AND. (optimizer%conjugator /=
cg_zero))
THEN
2022 cpabort(
"Cannot use conjugators with BFGS")
2025 CALL lbfgs_create(nlmo_lbfgs_history, nspins, nstore=10)
2028 IF (nspins == 1)
THEN
2029 spin_factor = 2.0_dp
2031 spin_factor = 1.0_dp
2034 ALLOCATE (grad_norm_spin(nspins))
2035 ALLOCATE (nocc(nspins))
2036 ALLOCATE (penalty_vol_prefactor(nspins))
2037 ALLOCATE (suggested_vol_penalty(nspins))
2043 ALLOCATE (m_t_mo_local(nspins))
2044 DO ispin = 1, nspins
2046 template=matrix_mo_in(ispin), &
2047 matrix_type=dbcsr_type_no_symmetry)
2048 CALL dbcsr_copy(m_t_mo_local(ispin), matrix_mo_in(ispin))
2051 ALLOCATE (approx_inv_hessian(nspins))
2052 ALLOCATE (m_theta_normalized(nspins))
2053 ALLOCATE (prev_m_theta(nspins))
2054 ALLOCATE (m_s0(nspins))
2055 ALLOCATE (prev_grad(nspins))
2056 ALLOCATE (grad(nspins))
2057 ALLOCATE (prev_step(nspins))
2058 ALLOCATE (step(nspins))
2059 ALLOCATE (prev_minus_prec_grad(nspins))
2060 ALLOCATE (m_sig_sqrti_ii(nspins))
2061 ALLOCATE (m_sigma(nspins))
2062 ALLOCATE (m_siginv(nspins))
2063 ALLOCATE (tempnocc1(nspins))
2064 ALLOCATE (tempoccocc1(nspins))
2065 ALLOCATE (tempoccocc2(nspins))
2066 ALLOCATE (tempoccocc3(nspins))
2067 ALLOCATE (bfgs_y(nspins))
2068 ALLOCATE (bfgs_s(nspins))
2070 DO ispin = 1, nspins
2074 template=matrix_mo_out(ispin), &
2075 matrix_type=dbcsr_type_no_symmetry)
2077 template=m_theta(ispin), &
2078 matrix_type=dbcsr_type_no_symmetry)
2080 template=m_theta(ispin), &
2081 matrix_type=dbcsr_type_no_symmetry)
2083 template=m_theta(ispin), &
2084 matrix_type=dbcsr_type_no_symmetry)
2086 template=m_theta(ispin), &
2087 matrix_type=dbcsr_type_no_symmetry)
2089 template=m_theta(ispin), &
2090 matrix_type=dbcsr_type_no_symmetry)
2092 template=m_theta(ispin), &
2093 matrix_type=dbcsr_type_no_symmetry)
2095 template=m_theta(ispin), &
2096 matrix_type=dbcsr_type_no_symmetry)
2098 template=m_theta(ispin), &
2099 matrix_type=dbcsr_type_no_symmetry)
2101 template=m_theta(ispin), &
2102 matrix_type=dbcsr_type_no_symmetry)
2104 template=m_theta(ispin), &
2105 matrix_type=dbcsr_type_no_symmetry)
2107 template=m_theta(ispin), &
2108 matrix_type=dbcsr_type_no_symmetry)
2110 template=m_theta(ispin), &
2111 matrix_type=dbcsr_type_no_symmetry)
2113 template=m_theta(ispin), &
2114 matrix_type=dbcsr_type_no_symmetry)
2116 template=m_theta(ispin), &
2117 matrix_type=dbcsr_type_no_symmetry)
2119 template=m_theta(ispin), &
2120 matrix_type=dbcsr_type_no_symmetry)
2122 template=m_theta(ispin), &
2123 matrix_type=dbcsr_type_no_symmetry)
2125 template=m_theta(ispin), &
2126 matrix_type=dbcsr_type_no_symmetry)
2129 CALL dbcsr_set(prev_step(ispin), 0.0_dp)
2132 nfullrows_total=nocc(ispin))
2134 penalty_vol_prefactor(ispin) = -penalty_amplitude
2139 m_t_mo_local(ispin), &
2140 0.0_dp, tempnocc1(ispin), &
2141 filter_eps=eps_filter)
2143 m_t_mo_local(ispin), &
2145 0.0_dp, m_s0(ispin), &
2146 filter_eps=eps_filter)
2148 SELECT CASE (optimizer%opt_penalty%operator_type)
2153 DO idim0 = 1,
SIZE(op_sm_set_qs, 2)
2155 DO reim = 1,
SIZE(op_sm_set_qs, 1)
2158 op_sm_set_almo(reim, idim0)%matrix, mat_distr_aos)
2161 op_sm_set_almo(reim, idim0)%matrix, &
2162 m_t_mo_local(ispin), &
2163 0.0_dp, tempnocc1(ispin), &
2164 filter_eps=eps_filter)
2167 m_t_mo_local(ispin), &
2169 0.0_dp, m_b0(reim, idim0, ispin), &
2170 filter_eps=eps_filter)
2172 DEALLOCATE (op_sm_set_qs(reim, idim0)%matrix)
2173 DEALLOCATE (op_sm_set_almo(reim, idim0)%matrix)
2184 isgf = first_sgf(iatom)
2189 m_t_mo_local(ispin), &
2190 0.0_dp, tempnocc1(ispin), &
2191 filter_eps=eps_filter)
2194 m_t_mo_local(ispin), &
2196 0.0_dp, m_b0(1, iatom, ispin), &
2197 first_k=isgf, last_k=isgf + ncol - 1, &
2198 filter_eps=eps_filter)
2202 m_t_mo_local(ispin), &
2203 0.0_dp, tempnocc1(ispin), &
2204 first_k=isgf, last_k=isgf + ncol - 1, &
2205 filter_eps=eps_filter)
2208 m_t_mo_local(ispin), &
2210 1.0_dp, m_b0(1, iatom, ispin), &
2211 filter_eps=eps_filter)
2219 IF (optimizer%opt_penalty%operator_type ==
op_loc_berry)
THEN
2220 DO idim0 = 1,
SIZE(op_sm_set_qs, 2)
2221 DO reim = 1,
SIZE(op_sm_set_qs, 1)
2222 DEALLOCATE (op_sm_set_qs(reim, idim0)%matrix)
2223 DEALLOCATE (op_sm_set_almo(reim, idim0)%matrix)
2226 DEALLOCATE (op_sm_set_qs, op_sm_set_almo)
2230 outer_max_iter = optimizer%max_iter_outer_loop
2231 outer_prepare_to_exit = .false.
2234 penalty_func_new = 0.0_dp
2235 linear_search_type = 1
2236 localization_obj_function = 0.0_dp
2237 penalty_func_new = 0.0_dp
2242 max_iter = optimizer%max_iter
2243 prepare_to_exit = .false.
2244 line_search = .false.
2248 line_search_iteration = 0
2249 obj_function_ispin = 0.0_dp
2253 line_search_error = 0.0_dp
2255 next_step_size_guess = 0.0_dp
2259 just_started = (iteration == 0) .AND. (outer_iteration == 0)
2261 DO ispin = 1, nspins
2267 m_s0(ispin), m_theta(ispin), 0.0_dp, &
2268 tempoccocc1(ispin), &
2269 filter_eps=eps_filter)
2270 CALL dbcsr_set(m_sig_sqrti_ii(ispin), 0.0_dp)
2273 m_theta(ispin), tempoccocc1(ispin), 0.0_dp, &
2274 m_sig_sqrti_ii(ispin), &
2275 retain_sparsity=.true.)
2276 ALLOCATE (diagonal(nocc(ispin)))
2278 CALL group%sum(diagonal)
2280 diagonal(:) = 1.0_dp/sqrt(diagonal(:))
2281 CALL dbcsr_set(m_sig_sqrti_ii(ispin), 0.0_dp)
2283 DEALLOCATE (diagonal)
2287 m_sig_sqrti_ii(ispin), &
2288 0.0_dp, m_theta_normalized(ispin), &
2289 filter_eps=eps_filter)
2293 m_t_mo_local(ispin), &
2294 m_theta_normalized(ispin), &
2295 0.0_dp, matrix_mo_out(ispin), &
2296 filter_eps=eps_filter)
2301 localization_obj_function = 0.0_dp
2302 penalty_func_new = 0.0_dp
2303 DO ispin = 1, nspins
2305 CALL compute_obj_nlmos( &
2306 localization_obj_function_ispin=localization_obj_function_ispin, &
2307 penalty_func_ispin=penalty_func_ispin, &
2308 overlap_determinant=overlap_determinant, &
2309 m_sigma=m_sigma(ispin), &
2311 m_b0=m_b0(:, :, ispin), &
2312 m_theta_normalized=m_theta_normalized(ispin), &
2313 template_matrix_mo=matrix_mo_out(ispin), &
2316 just_started=just_started, &
2317 penalty_vol_prefactor=penalty_vol_prefactor(ispin), &
2318 penalty_amplitude=penalty_amplitude, &
2319 eps_filter=eps_filter)
2321 localization_obj_function = localization_obj_function + localization_obj_function_ispin
2322 penalty_func_new = penalty_func_new + penalty_func_ispin
2325 objf_new = penalty_func_new + localization_obj_function
2327 DO ispin = 1, nspins
2331 IF (line_search_iteration == 0 .AND. iteration /= 0)
THEN
2332 CALL dbcsr_copy(prev_grad(ispin), grad(ispin))
2338 DO ispin = 1, nspins
2341 matrix_inverse=m_siginv(ispin), &
2342 matrix=m_sigma(ispin), &
2343 threshold=eps_filter*10.0_dp, &
2344 filter_eps=eps_filter, &
2347 CALL compute_gradient_nlmos( &
2348 m_grad_out=grad(ispin), &
2349 m_b0=m_b0(:, :, ispin), &
2352 m_theta_normalized=m_theta_normalized(ispin), &
2353 m_siginv=m_siginv(ispin), &
2354 m_sig_sqrti_ii=m_sig_sqrti_ii(ispin), &
2355 penalty_vol_prefactor=penalty_vol_prefactor(ispin), &
2356 eps_filter=eps_filter, &
2357 suggested_vol_penalty=suggested_vol_penalty(ispin))
2362 DO ispin = 1, nspins
2365 grad_norm = maxval(grad_norm_spin)
2367 converged = (grad_norm <= optimizer%eps_error)
2368 IF (converged .OR. (iteration >= max_iter))
THEN
2369 prepare_to_exit = .true.
2373 IF (.NOT. prepare_to_exit)
THEN
2378 IF (iteration /= 0)
THEN
2382 IF (.NOT. line_search)
THEN
2384 line_search = .true.
2385 line_search_iteration = line_search_iteration + 1
2391 line_search_error = 0.0_dp
2395 DO ispin = 1, nspins
2397 CALL dbcsr_dot(grad(ispin), step(ispin), tempreal)
2398 line_search_error = line_search_error + tempreal
2399 CALL dbcsr_dot(grad(ispin), grad(ispin), tempreal)
2400 denom = denom + tempreal
2401 CALL dbcsr_dot(step(ispin), step(ispin), tempreal)
2402 denom2 = denom2 + tempreal
2408 line_search_error = line_search_error/sqrt(denom)/sqrt(denom2)
2410 IF (abs(line_search_error) > optimizer%lin_search_eps_error)
THEN
2411 line_search = .true.
2412 line_search_iteration = line_search_iteration + 1
2414 line_search = .false.
2415 line_search_iteration = 0
2422 IF (line_search)
THEN
2425 objf_diff = objf_new - objf_old
2430 IF (.NOT. line_search)
THEN
2432 cg_iteration = cg_iteration + 1
2435 DO ispin = 1, nspins
2436 CALL dbcsr_copy(prev_step(ispin), step(ispin))
2444 DO ispin = 1, nspins
2454 IF (iteration == 0)
THEN
2460 IF (nspins > 1)
THEN
2461 DO ispin = 2, nspins
2462 CALL dbcsr_copy(approx_inv_hessian(ispin), approx_inv_hessian(1))
2466 ELSE IF (l_bfgs)
THEN
2468 CALL lbfgs_seed(nlmo_lbfgs_history, m_theta, grad)
2469 DO ispin = 1, nspins
2477 DO ispin = 1, nspins
2491 DO ispin = 1, nspins
2495 CALL dbcsr_add(bfgs_y(ispin), prev_grad(ispin), 1.0_dp, -1.0_dp)
2496 CALL dbcsr_copy(bfgs_s(ispin), m_theta(ispin))
2497 CALL dbcsr_add(bfgs_s(ispin), prev_m_theta(ispin), 1.0_dp, -1.0_dp)
2500 CALL dbcsr_dot(grad(ispin), step(ispin), bfgs_rho)
2501 bfgs_rho = 1.0_dp/bfgs_rho
2504 CALL dbcsr_dot(bfgs_y(ispin), bfgs_y(ispin), bfgs_sum)
2507 CALL dbcsr_copy(tempoccocc2(ispin), approx_inv_hessian(ispin))
2511 CALL dbcsr_add(tempoccocc2(ispin), tempoccocc1(ispin), 1.0_dp, bfgs_rho)
2515 approx_inv_hessian(ispin), tempoccocc3(ispin))
2516 CALL dbcsr_add(tempoccocc2(ispin), tempoccocc3(ispin), &
2517 1.0_dp, bfgs_rho*bfgs_rho*bfgs_sum)
2521 approx_inv_hessian(ispin), tempoccocc1(ispin))
2523 CALL dbcsr_add(tempoccocc2(ispin), tempoccocc3(ispin), &
2524 1.0_dp, -2.0_dp*bfgs_rho)
2526 CALL dbcsr_copy(approx_inv_hessian(ispin), tempoccocc2(ispin))
2530 ELSE IF (l_bfgs)
THEN
2538 IF (.NOT. l_bfgs)
THEN
2540 DO ispin = 1, nspins
2543 grad(ispin), step(ispin))
2553 IF (iteration == 0)
THEN
2554 reset_conjugator = .true.
2558 IF (.NOT. reset_conjugator)
THEN
2559 CALL compute_cg_beta( &
2561 reset_conjugator=reset_conjugator, &
2562 conjugator=optimizer%conjugator, &
2564 prev_grad=prev_grad(:), &
2566 prev_step=prev_step(:), &
2567 prev_minus_prec_grad=prev_minus_prec_grad(:) &
2572 IF (reset_conjugator)
THEN
2575 IF (unit_nr > 0 .AND. (.NOT. just_started))
THEN
2576 WRITE (unit_nr,
'(T2,A35)')
"Re-setting conjugator to zero"
2578 reset_conjugator = .false.
2583 DO ispin = 1, nspins
2585 CALL dbcsr_copy(prev_minus_prec_grad(ispin), step(ispin))
2588 CALL dbcsr_add(step(ispin), prev_step(ispin), 1.0_dp, beta)
2595 IF (.NOT. line_search)
THEN
2601 DO ispin = 1, nspins
2602 CALL dbcsr_dot(grad(ispin), step(ispin), tempreal)
2605 g0sign = sign(1.0_dp, g0)
2606 IF (linear_search_type == 1)
THEN
2607 IF (iteration == 0)
THEN
2608 step_size = optimizer%lin_search_step_size_guess
2610 IF (next_step_size_guess <= 0.0_dp)
THEN
2611 step_size = optimizer%lin_search_step_size_guess
2614 step_size = optimizer%lin_search_step_size_guess
2618 ELSE IF (linear_search_type == 2)
THEN
2621 step_size = optimizer%lin_search_step_size_guess
2623 IF (unit_nr > 0)
THEN
2624 WRITE (unit_nr,
'(T21,3A19)')
"Line position",
"Line grad",
"Next line step"
2625 WRITE (unit_nr,
'(T2,A19,3F19.5)')
"Line search", 0.0_dp, g0, step_size
2627 next_step_size_guess = step_size
2631 DO ispin = 1, nspins
2632 CALL dbcsr_dot(grad(ispin), step(ispin), tempreal)
2635 g1sign = sign(1.0_dp, g1)
2636 IF (linear_search_type == 1)
THEN
2639 appr_sec_der = (g1 - g0)/step_size
2640 step_size = -g1/appr_sec_der
2641 ELSE IF (linear_search_type == 2)
THEN
2644 IF (g1sign /= g0sign)
THEN
2645 step_size = -step_size/2.0
2647 step_size = step_size*1.5
2651 IF (unit_nr > 0)
THEN
2652 WRITE (unit_nr,
'(T21,3A19)')
"Line position",
"Line grad",
"Next line step"
2653 WRITE (unit_nr,
'(T2,A19,3F19.5)')
"Line search", next_step_size_guess, g1, step_size
2658 next_step_size_guess = next_step_size_guess + step_size
2662 DO ispin = 1, nspins
2663 IF (.NOT. line_search)
THEN
2665 CALL dbcsr_copy(prev_m_theta(ispin), m_theta(ispin))
2667 CALL dbcsr_add(m_theta(ispin), step(ispin), 1.0_dp, step_size)
2672 IF (line_search)
THEN
2679 IF (unit_nr > 0)
THEN
2680 iter_type = trim(
"NLMO OPT "//iter_type)
2681 WRITE (unit_nr,
'(T2,A13,I6,F23.10,E14.5,F14.9,F9.2)') &
2682 iter_type, iteration, &
2683 objf_new, objf_diff, grad_norm, &
2685 WRITE (unit_nr,
'(T2,A19,F23.10)') &
2686 "Localization:", localization_obj_function
2687 WRITE (unit_nr,
'(T2,A19,F23.10)') &
2688 "Orthogonalization:", penalty_func_new
2692 iteration = iteration + 1
2693 IF (prepare_to_exit)
EXIT
2697 IF (converged .OR. (outer_iteration >= outer_max_iter))
THEN
2698 outer_prepare_to_exit = .true.
2701 outer_iteration = outer_iteration + 1
2702 IF (outer_prepare_to_exit)
EXIT
2707 optimizer%opt_penalty%penalty_strength = 0.0_dp
2708 DO ispin = 1, nspins
2709 optimizer%opt_penalty%penalty_strength = optimizer%opt_penalty%penalty_strength + &
2710 (-1.0_dp)*penalty_vol_prefactor(ispin)
2712 optimizer%opt_penalty%penalty_strength = optimizer%opt_penalty%penalty_strength/nspins
2717 iter_type =
"Unconverged"
2720 IF (unit_nr > 0)
THEN
2721 WRITE (unit_nr,
'()')
2722 print_string = trim(iter_type)//
" localization:"
2723 WRITE (unit_nr,
'(T2,A29,F30.10)') &
2724 print_string, localization_obj_function
2725 print_string = trim(iter_type)//
" determinant:"
2726 WRITE (unit_nr,
'(T2,A29,F30.10)') &
2727 print_string, overlap_determinant
2728 print_string = trim(iter_type)//
" penalty strength:"
2729 WRITE (unit_nr,
'(T2,A29,F30.10)') &
2730 print_string, optimizer%opt_penalty%penalty_strength
2737 DO ispin = 1, nspins
2738 DO idim0 = 1,
SIZE(m_b0, 2)
2739 DO reim = 1,
SIZE(m_b0, 1)
2765 DEALLOCATE (grad_norm_spin)
2767 DEALLOCATE (penalty_vol_prefactor)
2768 DEALLOCATE (suggested_vol_penalty)
2770 DEALLOCATE (approx_inv_hessian)
2771 DEALLOCATE (prev_m_theta)
2772 DEALLOCATE (m_theta_normalized)
2774 DEALLOCATE (prev_grad)
2776 DEALLOCATE (prev_step)
2778 DEALLOCATE (prev_minus_prec_grad)
2779 DEALLOCATE (m_sig_sqrti_ii)
2780 DEALLOCATE (m_sigma)
2781 DEALLOCATE (m_siginv)
2782 DEALLOCATE (tempnocc1)
2783 DEALLOCATE (tempoccocc1)
2784 DEALLOCATE (tempoccocc2)
2785 DEALLOCATE (tempoccocc3)
2789 DEALLOCATE (m_theta, m_t_mo_local)
2791 DEALLOCATE (weights)
2792 DEALLOCATE (first_sgf, last_sgf, nsgf)
2794 IF (.NOT. converged)
THEN
2795 cpabort(
"Optimization not converged! ")
2798 CALL timestop(handle)
2820 SUBROUTINE xalmo_analysis(detailed_analysis, eps_filter, m_T_in, m_T0_in, &
2821 m_siginv_in, m_siginv0_in, m_S_in, m_KS0_in, m_quench_t_in, energy_out, &
2822 m_eda_out, m_cta_out)
2824 LOGICAL,
INTENT(IN) :: detailed_analysis
2825 REAL(kind=
dp),
INTENT(IN) :: eps_filter
2826 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(IN) :: m_t_in, m_t0_in, m_siginv_in, &
2827 m_siginv0_in, m_s_in, m_ks0_in, &
2829 REAL(kind=
dp),
INTENT(INOUT) :: energy_out
2830 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(INOUT) :: m_eda_out, m_cta_out
2832 CHARACTER(len=*),
PARAMETER :: routinen =
'xalmo_analysis'
2834 INTEGER :: handle, ispin, nspins
2835 REAL(kind=
dp) :: energy_ispin, spin_factor
2836 TYPE(
dbcsr_type) :: ftsiginv0, fvo0, m_x, siginvtftsiginv0, &
2839 CALL timeset(routinen, handle)
2841 nspins =
SIZE(m_t_in)
2843 IF (nspins == 1)
THEN
2844 spin_factor = 2.0_dp
2846 spin_factor = 1.0_dp
2850 DO ispin = 1, nspins
2854 template=m_t_in(ispin), &
2855 matrix_type=dbcsr_type_no_symmetry)
2857 template=m_t_in(ispin), &
2858 matrix_type=dbcsr_type_no_symmetry)
2860 template=m_t_in(ispin), &
2861 matrix_type=dbcsr_type_no_symmetry)
2863 template=m_t_in(ispin), &
2864 matrix_type=dbcsr_type_no_symmetry)
2866 template=m_siginv0_in(ispin), &
2867 matrix_type=dbcsr_type_no_symmetry)
2870 CALL compute_frequently_used_matrices( &
2871 filter_eps=eps_filter, &
2872 m_t_in=m_t0_in(ispin), &
2873 m_siginv_in=m_siginv0_in(ispin), &
2875 m_f_in=m_ks0_in(ispin), &
2876 m_ftsiginv_out=ftsiginv0, &
2877 m_siginvtftsiginv_out=siginvtftsiginv0, &
2880 CALL dbcsr_copy(fvo0, ftsiginv0, keep_sparsity=.true.)
2885 retain_sparsity=.true.)
2889 CALL dbcsr_add(m_x, m_t_in(ispin), -1.0_dp, 1.0_dp)
2892 energy_out = energy_out + energy_ispin*spin_factor
2894 IF (detailed_analysis)
THEN
2904 m_siginv0_in(ispin), &
2905 0.0_dp, ftsiginv0, &
2906 filter_eps=eps_filter)
2912 filter_eps=eps_filter)
2917 0.0_dp, siginvtftsiginv0, &
2918 filter_eps=eps_filter)
2925 filter_eps=eps_filter)
2929 m_siginv_in(ispin), &
2930 0.0_dp, ftsiginv0, &
2931 filter_eps=eps_filter)
2934 ftsiginv0, m_cta_out(ispin))
2948 CALL timestop(handle)
2950 END SUBROUTINE xalmo_analysis
2967 SUBROUTINE compute_frequently_used_matrices(filter_eps, &
2968 m_T_in, m_siginv_in, m_S_in, m_F_in, m_FTsiginv_out, &
2969 m_siginvTFTsiginv_out, m_ST_out)
2971 REAL(kind=
dp),
INTENT(IN) :: filter_eps
2972 TYPE(
dbcsr_type),
INTENT(IN) :: m_t_in, m_siginv_in, m_s_in, m_f_in
2973 TYPE(
dbcsr_type),
INTENT(INOUT) :: m_ftsiginv_out, m_siginvtftsiginv_out, &
2976 CHARACTER(len=*),
PARAMETER :: routinen =
'compute_frequently_used_matrices'
2981 CALL timeset(routinen, handle)
2985 matrix_type=dbcsr_type_no_symmetry)
2987 template=m_siginv_in, &
2988 matrix_type=dbcsr_type_no_symmetry)
2993 0.0_dp, m_tmp_no_1, &
2994 filter_eps=filter_eps)
2999 0.0_dp, m_ftsiginv_out, &
3000 filter_eps=filter_eps)
3005 0.0_dp, m_tmp_oo_1, &
3006 filter_eps=filter_eps)
3011 0.0_dp, m_siginvtftsiginv_out, &
3012 filter_eps=filter_eps)
3018 filter_eps=filter_eps)
3023 CALL timestop(handle)
3025 END SUBROUTINE compute_frequently_used_matrices
3035 SUBROUTINE split_v_blk(almo_scf_env)
3039 CHARACTER(len=*),
PARAMETER :: routinen =
'split_v_blk'
3041 INTEGER :: discarded_v, handle, iblock_col, &
3042 iblock_col_size, iblock_row, &
3043 iblock_row_size, ispin, retained_v
3044 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: data_p
3047 CALL timeset(routinen, handle)
3049 DO ispin = 1, almo_scf_env%nspins
3052 work_mutable=.true.)
3054 work_mutable=.true.)
3061 row_size=iblock_row_size, col_size=iblock_col_size)
3063 IF (iblock_row /= iblock_col)
THEN
3064 cpabort(
"off-diagonal block found")
3067 retained_v = almo_scf_env%nvirt_of_domain(iblock_col, ispin)
3068 discarded_v = almo_scf_env%nvirt_disc_of_domain(iblock_col, ispin)
3069 cpassert(retained_v > 0)
3070 cpassert(discarded_v > 0)
3071 CALL dbcsr_put_block(almo_scf_env%matrix_v_disc_blk(ispin), iblock_row, iblock_col, &
3072 block=data_p(:, (retained_v + 1):iblock_col_size))
3073 CALL dbcsr_put_block(almo_scf_env%matrix_v_blk(ispin), iblock_row, iblock_col, &
3074 block=data_p(:, 1:retained_v))
3084 CALL timestop(handle)
3086 END SUBROUTINE split_v_blk
3095 SUBROUTINE harris_foulkes_correction(almo_scf_env)
3099 CHARACTER(len=*),
PARAMETER :: routinen =
'harris_foulkes_correction'
3100 INTEGER,
PARAMETER :: cayley_transform = 1, dm_ls_step = 2
3102 INTEGER :: algorithm_id, handle, handle1, handle2, handle3, handle4, handle5, handle6, &
3103 handle7, handle8, ispin, iteration, n, nmins, nspin, opt_k_max_iter, &
3104 outer_opt_k_iteration, outer_opt_k_max_iter, unit_nr
3105 INTEGER,
DIMENSION(1) :: fake, nelectron_spin_real
3106 LOGICAL :: converged, line_search, md_in_k_space, outer_opt_k_prepare_to_exit, &
3107 prepare_to_exit, reset_conjugator, reset_step_size, use_cubic_approximation, &
3108 use_quadratic_approximation
3109 REAL(kind=
dp) :: aa, bb, beta, conjugacy_error, conjugacy_error_threshold, &
3110 delta_obj_function, denom, energy_correction_final, frob_matrix, frob_matrix_base, fun0, &
3111 fun1, gfun0, gfun1, grad_norm, grad_norm_frob, kappa, kin_energy, line_search_error, &
3112 line_search_error_threshold, num_threshold, numer, obj_function, quadratic_approx_error, &
3113 quadratic_approx_error_threshold, safety_multiplier, spin_factor, step_size, &
3114 step_size_quadratic_approx, step_size_quadratic_approx2, t1, t1a, t1cholesky, t2, t2a, &
3115 t2cholesky, tau, time_step, x_opt_eps_adaptive, x_opt_eps_adaptive_factor
3116 REAL(kind=
dp),
DIMENSION(1) :: local_mu
3117 REAL(kind=
dp),
DIMENSION(2) :: energy_correction
3118 REAL(kind=
dp),
DIMENSION(3) :: minima
3121 TYPE(
dbcsr_type) :: grad, k_vd_index_down, k_vr_index_down, matrix_k_central, matrix_tmp1, &
3122 matrix_tmp2, prec, prev_grad, prev_minus_prec_grad, prev_step, sigma_oo_curr, &
3123 sigma_oo_curr_inv, sigma_vv_sqrt, sigma_vv_sqrt_guess, sigma_vv_sqrt_inv, &
3124 sigma_vv_sqrt_inv_guess, step, t_curr, tmp1_n_vr, tmp2_n_o, tmp3_vd_vr, tmp4_o_vr, &
3125 tmp_k_blk, vd_fixed, vd_index_sqrt, vd_index_sqrt_inv, velocity, vr_fixed, vr_index_sqrt, &
3127 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:) :: matrix_p_almo_scf_converged
3129 CALL timeset(routinen, handle)
3133 IF (logger%para_env%is_source())
THEN
3139 nspin = almo_scf_env%nspins
3140 energy_correction_final = 0.0_dp
3141 IF (nspin == 1)
THEN
3142 spin_factor = 2.0_dp
3144 spin_factor = 1.0_dp
3147 IF (almo_scf_env%deloc_use_occ_orbs)
THEN
3148 algorithm_id = cayley_transform
3150 algorithm_id = dm_ls_step
3155 SELECT CASE (algorithm_id)
3156 CASE (cayley_transform)
3160 IF (almo_scf_env%nspins == 1)
THEN
3161 CALL dbcsr_scale(almo_scf_env%matrix_p(1), 1.0_dp/spin_factor)
3167 CALL dbcsr_copy(almo_scf_env%matrix_t(ispin), &
3168 almo_scf_env%matrix_t_blk(ispin))
3175 IF (unit_nr > 0)
THEN
3176 WRITE (unit_nr, *)
"sqrt and inv(sqrt) of MO overlap matrix"
3178 CALL dbcsr_create(almo_scf_env%matrix_sigma_sqrt(ispin), &
3179 template=almo_scf_env%matrix_sigma(ispin), &
3180 matrix_type=dbcsr_type_no_symmetry)
3181 CALL dbcsr_create(almo_scf_env%matrix_sigma_sqrt_inv(ispin), &
3182 template=almo_scf_env%matrix_sigma(ispin), &
3183 matrix_type=dbcsr_type_no_symmetry)
3186 almo_scf_env%matrix_sigma_sqrt_inv(ispin), &
3187 almo_scf_env%matrix_sigma(ispin), &
3188 threshold=almo_scf_env%eps_filter, &
3189 order=almo_scf_env%order_lanczos, &
3190 eps_lanczos=almo_scf_env%eps_lanczos, &
3191 max_iter_lanczos=almo_scf_env%max_iter_lanczos)
3194 CALL dbcsr_create(matrix_tmp1, template=almo_scf_env%matrix_sigma(ispin), &
3195 matrix_type=dbcsr_type_no_symmetry)
3196 CALL dbcsr_create(matrix_tmp2, template=almo_scf_env%matrix_sigma(ispin), &
3197 matrix_type=dbcsr_type_no_symmetry)
3199 CALL dbcsr_multiply(
"N",
"N", 1.0_dp, almo_scf_env%matrix_sigma_sqrt_inv(ispin), &
3200 almo_scf_env%matrix_sigma(ispin), &
3201 0.0_dp, matrix_tmp1, filter_eps=almo_scf_env%eps_filter)
3203 almo_scf_env%matrix_sigma_sqrt_inv(ispin), &
3204 0.0_dp, matrix_tmp2, filter_eps=almo_scf_env%eps_filter)
3209 IF (unit_nr > 0)
THEN
3210 WRITE (unit_nr, *)
"Error for (inv(sqrt(SIG))*SIG*inv(sqrt(SIG))-I)", frob_matrix/frob_matrix_base
3218 IF (almo_scf_env%almo_update_algorithm ==
almo_scf_diag)
THEN
3224 line_search_error_threshold = almo_scf_env%real01
3225 conjugacy_error_threshold = almo_scf_env%real02
3226 quadratic_approx_error_threshold = almo_scf_env%real03
3227 x_opt_eps_adaptive_factor = almo_scf_env%real04
3230 outer_opt_k_max_iter = almo_scf_env%opt_k_outer_max_iter
3231 outer_opt_k_prepare_to_exit = .false.
3232 outer_opt_k_iteration = 0
3234 grad_norm_frob = 0.0_dp
3235 CALL dbcsr_set(almo_scf_env%matrix_x(ispin), 0.0_dp)
3236 IF (almo_scf_env%deloc_truncate_virt ==
virt_full) outer_opt_k_max_iter = 0
3242 psi_out=almo_scf_env%matrix_v(ispin), &
3243 psi_projector=almo_scf_env%matrix_t_blk(ispin), &
3244 metric=almo_scf_env%matrix_s(1), &
3245 project_out=.true., &
3246 psi_projector_orthogonal=.false., &
3247 proj_in_template=almo_scf_env%matrix_ov(ispin), &
3248 eps_filter=almo_scf_env%eps_filter, &
3249 sig_inv_projector=almo_scf_env%matrix_sigma_inv(ispin))
3253 template=almo_scf_env%matrix_v(ispin))
3254 CALL dbcsr_copy(vr_fixed, almo_scf_env%matrix_v(ispin))
3258 template=almo_scf_env%matrix_sigma_vv(ispin), &
3259 matrix_type=dbcsr_type_no_symmetry)
3261 template=almo_scf_env%matrix_sigma_vv(ispin), &
3262 matrix_type=dbcsr_type_no_symmetry)
3264 template=almo_scf_env%matrix_sigma_vv(ispin), &
3265 matrix_type=dbcsr_type_no_symmetry)
3267 template=almo_scf_env%matrix_sigma_vv(ispin), &
3268 matrix_type=dbcsr_type_no_symmetry)
3269 CALL dbcsr_set(sigma_vv_sqrt_guess, 0.0_dp)
3271 CALL dbcsr_filter(sigma_vv_sqrt_guess, almo_scf_env%eps_filter)
3272 CALL dbcsr_set(sigma_vv_sqrt_inv_guess, 0.0_dp)
3274 CALL dbcsr_filter(sigma_vv_sqrt_inv_guess, almo_scf_env%eps_filter)
3277 IF (almo_scf_env%deloc_truncate_virt /=
virt_full)
THEN
3285 psi_out=almo_scf_env%matrix_v_disc(ispin), &
3286 psi_projector=almo_scf_env%matrix_t_blk(ispin), &
3287 metric=almo_scf_env%matrix_s(1), &
3288 project_out=.true., &
3289 psi_projector_orthogonal=.false., &
3290 proj_in_template=almo_scf_env%matrix_ov_disc(ispin), &
3291 eps_filter=almo_scf_env%eps_filter, &
3292 sig_inv_projector=almo_scf_env%matrix_sigma_inv(ispin))
3297 template=almo_scf_env%matrix_v_disc(ispin))
3298 CALL dbcsr_copy(vd_fixed, almo_scf_env%matrix_v_disc(ispin))
3302 template=almo_scf_env%matrix_sigma_vv_blk(ispin), &
3303 matrix_type=dbcsr_type_no_symmetry)
3307 template=almo_scf_env%matrix_vv_disc_blk(ispin), &
3308 matrix_type=dbcsr_type_no_symmetry)
3314 template=almo_scf_env%matrix_k_blk(ispin))
3315 CALL dbcsr_copy(grad, almo_scf_env%matrix_k_blk(ispin))
3318 md_in_k_space = almo_scf_env%logical01
3319 IF (md_in_k_space)
THEN
3321 template=almo_scf_env%matrix_k_blk(ispin))
3322 CALL dbcsr_copy(velocity, almo_scf_env%matrix_k_blk(ispin))
3324 time_step = almo_scf_env%opt_k_trial_step_size
3328 template=almo_scf_env%matrix_k_blk(ispin))
3331 template=almo_scf_env%matrix_k_blk(ispin))
3335 template=almo_scf_env%matrix_k_blk(ispin))
3336 CALL dbcsr_copy(prec, almo_scf_env%matrix_k_blk(ispin))
3340 CALL dbcsr_set(almo_scf_env%matrix_k_blk(ispin), 0.0_dp)
3344 template=almo_scf_env%matrix_k_blk(ispin))
3346 almo_scf_env%matrix_k_blk(ispin))
3348 template=almo_scf_env%matrix_k_blk(ispin))
3350 template=almo_scf_env%matrix_k_blk(ispin))
3353 template=almo_scf_env%matrix_t(ispin))
3355 template=almo_scf_env%matrix_sigma(ispin), &
3356 matrix_type=dbcsr_type_no_symmetry)
3358 template=almo_scf_env%matrix_sigma(ispin), &
3359 matrix_type=dbcsr_type_no_symmetry)
3361 template=almo_scf_env%matrix_v(ispin))
3363 template=almo_scf_env%matrix_k_blk(ispin))
3365 template=almo_scf_env%matrix_t(ispin))
3367 template=almo_scf_env%matrix_ov(ispin))
3369 template=almo_scf_env%matrix_k_blk(ispin))
3375 opt_k_max_iter = almo_scf_env%opt_k_max_iter
3378 prepare_to_exit = .false.
3380 line_search = .false.
3381 obj_function = 0.0_dp
3382 conjugacy_error = 0.0_dp
3383 line_search_error = 0.0_dp
3388 step_size_quadratic_approx = 0.0_dp
3389 reset_step_size = .true.
3390 IF (almo_scf_env%deloc_truncate_virt ==
virt_full) opt_k_max_iter = 0
3395 CALL timeset(
'k_opt_vr', handle1)
3397 IF (almo_scf_env%deloc_truncate_virt /=
virt_full)
THEN
3401 almo_scf_env%matrix_k_blk(ispin), &
3402 0.0_dp, almo_scf_env%matrix_v(ispin), &
3403 filter_eps=almo_scf_env%eps_filter)
3404 CALL dbcsr_add(almo_scf_env%matrix_v(ispin), vr_fixed, &
3409 CALL get_overlap(bra=almo_scf_env%matrix_v(ispin), &
3410 ket=almo_scf_env%matrix_v(ispin), &
3411 overlap=almo_scf_env%matrix_sigma_vv(ispin), &
3412 metric=almo_scf_env%matrix_s(1), &
3413 retain_overlap_sparsity=.false., &
3414 eps_filter=almo_scf_env%eps_filter)
3417 IF (almo_scf_env%deloc_truncate_virt ==
virt_full)
THEN
3418 CALL timeset(
'cholesky', handle2)
3424 template=almo_scf_env%matrix_sigma_vv(ispin), &
3425 matrix_type=dbcsr_type_no_symmetry)
3429 para_env=almo_scf_env%para_env, &
3430 blacs_env=almo_scf_env%blacs_env)
3431 CALL make_triu(sigma_vv_sqrt)
3432 CALL dbcsr_filter(sigma_vv_sqrt, almo_scf_env%eps_filter)
3435 CALL dbcsr_create(matrix_tmp1, template=almo_scf_env%matrix_sigma_vv(ispin), &
3436 matrix_type=dbcsr_type_no_symmetry)
3440 sigma_vv_sqrt_inv, op=
"SOLVE", pos=
"RIGHT", &
3441 para_env=almo_scf_env%para_env, &
3442 blacs_env=almo_scf_env%blacs_env)
3443 CALL dbcsr_filter(sigma_vv_sqrt_inv, almo_scf_env%eps_filter)
3446 CALL dbcsr_create(matrix_tmp1, template=almo_scf_env%matrix_sigma_vv(ispin), &
3447 matrix_type=dbcsr_type_no_symmetry)
3452 -1.0_dp, matrix_tmp1, filter_eps=almo_scf_env%eps_filter)
3456 IF (unit_nr > 0)
THEN
3457 WRITE (unit_nr, *)
"Error for ( U^T * U - Sig )", &
3458 frob_matrix/frob_matrix_base
3462 0.0_dp, matrix_tmp1, filter_eps=almo_scf_env%eps_filter)
3466 IF (unit_nr > 0)
THEN
3467 WRITE (unit_nr, *)
"Error for ( inv(U) * U - I )", &
3468 frob_matrix/frob_matrix_base
3473 IF (unit_nr > 0)
THEN
3474 WRITE (unit_nr, *)
"Cholesky+inverse wall-time: ", t2cholesky - t1cholesky
3476 CALL timestop(handle2)
3479 sigma_vv_sqrt_inv, &
3480 almo_scf_env%matrix_sigma_vv(ispin), &
3481 threshold=almo_scf_env%eps_filter, &
3482 order=almo_scf_env%order_lanczos, &
3483 eps_lanczos=almo_scf_env%eps_lanczos, &
3484 max_iter_lanczos=almo_scf_env%max_iter_lanczos)
3485 CALL dbcsr_copy(sigma_vv_sqrt_inv_guess, sigma_vv_sqrt_inv)
3486 CALL dbcsr_copy(sigma_vv_sqrt_guess, sigma_vv_sqrt)
3488 CALL dbcsr_create(matrix_tmp1, template=almo_scf_env%matrix_sigma_vv(ispin), &
3489 matrix_type=dbcsr_type_no_symmetry)
3490 CALL dbcsr_create(matrix_tmp2, template=almo_scf_env%matrix_sigma_vv(ispin), &
3491 matrix_type=dbcsr_type_no_symmetry)
3494 almo_scf_env%matrix_sigma_vv(ispin), &
3495 0.0_dp, matrix_tmp1, filter_eps=almo_scf_env%eps_filter)
3497 sigma_vv_sqrt_inv, &
3498 0.0_dp, matrix_tmp2, filter_eps=almo_scf_env%eps_filter)
3503 IF (unit_nr > 0)
THEN
3504 WRITE (unit_nr, *)
"Error for (inv(sqrt(SIGVV))*SIGVV*inv(sqrt(SIGVV))-I)", &
3505 frob_matrix/frob_matrix_base
3512 CALL timestop(handle1)
3516 IF ((iteration == 0) .AND. (.NOT. line_search) .AND. &
3517 (outer_opt_k_iteration == 0))
THEN
3518 x_opt_eps_adaptive = &
3519 almo_scf_env%deloc_cayley_eps_convergence
3521 x_opt_eps_adaptive = &
3522 max(abs(almo_scf_env%deloc_cayley_eps_convergence), &
3523 abs(x_opt_eps_adaptive_factor*grad_norm))
3527 para_env=almo_scf_env%para_env, &
3528 blacs_env=almo_scf_env%blacs_env, &
3529 use_occ_orbs=.true., &
3530 use_virt_orbs=.true., &
3531 occ_orbs_orthogonal=.false., &
3532 virt_orbs_orthogonal=.false., &
3533 pp_preconditioner_full=almo_scf_env%deloc_cayley_occ_precond, &
3534 qq_preconditioner_full=almo_scf_env%deloc_cayley_vir_precond, &
3535 tensor_type=almo_scf_env%deloc_cayley_tensor_type, &
3536 neglect_quadratic_term=almo_scf_env%deloc_cayley_linear, &
3537 conjugator=almo_scf_env%deloc_cayley_conjugator, &
3538 max_iter=almo_scf_env%deloc_cayley_max_iter, &
3539 calculate_energy_corr=.true., &
3542 eps_convergence=x_opt_eps_adaptive, &
3543 eps_filter=almo_scf_env%eps_filter, &
3545 q_index_up=sigma_vv_sqrt_inv, &
3546 q_index_down=sigma_vv_sqrt, &
3547 p_index_up=almo_scf_env%matrix_sigma_sqrt_inv(ispin), &
3548 p_index_down=almo_scf_env%matrix_sigma_sqrt(ispin), &
3549 matrix_ks=almo_scf_env%matrix_ks_0deloc(ispin), &
3550 matrix_t=almo_scf_env%matrix_t(ispin), &
3551 matrix_qp_template=almo_scf_env%matrix_vo(ispin), &
3552 matrix_pq_template=almo_scf_env%matrix_ov(ispin), &
3553 matrix_v=almo_scf_env%matrix_v(ispin), &
3554 matrix_x_guess=almo_scf_env%matrix_x(ispin))
3559 energy_correction=energy_correction(ispin), &
3560 copy_matrix_x=almo_scf_env%matrix_x(ispin))
3564 energy_correction(1) = energy_correction(1)*spin_factor
3566 IF (opt_k_max_iter /= 0)
THEN
3568 CALL timeset(
'k_opt_t_curr', handle3)
3572 almo_scf_env%matrix_v(ispin), &
3573 almo_scf_env%matrix_x(ispin), &
3575 filter_eps=almo_scf_env%eps_filter)
3576 CALL dbcsr_add(t_curr, almo_scf_env%matrix_t_blk(ispin), &
3582 overlap=sigma_oo_curr, &
3583 metric=almo_scf_env%matrix_s(1), &
3584 retain_overlap_sparsity=.false., &
3585 eps_filter=almo_scf_env%eps_filter)
3586 IF (iteration == 0)
THEN
3589 threshold=almo_scf_env%eps_filter, &
3590 use_inv_as_guess=.false.)
3594 threshold=almo_scf_env%eps_filter, &
3595 use_inv_as_guess=.true.)
3598 CALL dbcsr_create(matrix_tmp1, template=sigma_oo_curr, &
3599 matrix_type=dbcsr_type_no_symmetry)
3601 sigma_oo_curr_inv, &
3602 0.0_dp, matrix_tmp1, &
3603 filter_eps=almo_scf_env%eps_filter)
3607 IF (unit_nr > 0)
THEN
3608 WRITE (unit_nr, *)
"Error for (SIG*inv(SIG)-I)", &
3609 frob_matrix/frob_matrix_base, frob_matrix_base
3614 CALL dbcsr_create(matrix_tmp1, template=sigma_oo_curr, &
3615 matrix_type=dbcsr_type_no_symmetry)
3618 0.0_dp, matrix_tmp1, &
3619 filter_eps=almo_scf_env%eps_filter)
3623 IF (unit_nr > 0)
THEN
3624 WRITE (unit_nr, *)
"Error for (inv(SIG)*SIG-I)", &
3625 frob_matrix/frob_matrix_base, frob_matrix_base
3630 CALL timestop(handle3)
3631 CALL timeset(
'k_opt_vd', handle4)
3639 sigma_vv_sqrt_inv, &
3640 sigma_vv_sqrt_inv, &
3641 0.0_dp, sigma_vv_sqrt, &
3642 filter_eps=almo_scf_env%eps_filter)
3644 psi_out=almo_scf_env%matrix_v_disc(ispin), &
3645 psi_projector=almo_scf_env%matrix_v(ispin), &
3646 metric=almo_scf_env%matrix_s(1), &
3647 project_out=.false., &
3648 psi_projector_orthogonal=.false., &
3649 proj_in_template=almo_scf_env%matrix_k_tr(ispin), &
3650 eps_filter=almo_scf_env%eps_filter, &
3651 sig_inv_projector=sigma_vv_sqrt)
3653 CALL dbcsr_add(almo_scf_env%matrix_v_disc(ispin), &
3654 vd_fixed, -1.0_dp, +1.0_dp)
3656 CALL timestop(handle4)
3657 CALL timeset(
'k_opt_grad', handle5)
3662 IF (line_search)
THEN
3666 almo_scf_env%matrix_ks_0deloc(ispin), &
3669 filter_eps=almo_scf_env%eps_filter)
3671 sigma_oo_curr_inv, &
3672 almo_scf_env%matrix_x(ispin), &
3673 0.0_dp, tmp4_o_vr, &
3674 filter_eps=almo_scf_env%eps_filter)
3678 0.0_dp, tmp1_n_vr, &
3679 filter_eps=almo_scf_env%eps_filter)
3681 almo_scf_env%matrix_v_disc(ispin), &
3684 retain_sparsity=.true.)
3691 converged = (grad_norm < almo_scf_env%opt_k_eps_convergence)
3692 IF (converged .OR. (iteration >= opt_k_max_iter))
THEN
3693 prepare_to_exit = .true.
3695 CALL timestop(handle5)
3697 IF (.NOT. prepare_to_exit)
THEN
3699 CALL timeset(
'k_opt_energy', handle6)
3705 0.0_dp, sigma_oo_curr, &
3706 filter_eps=almo_scf_env%eps_filter)
3707 delta_obj_function = fun0
3708 CALL dbcsr_dot(sigma_oo_curr_inv, sigma_oo_curr, obj_function)
3709 delta_obj_function = obj_function - delta_obj_function
3710 IF (line_search)
THEN
3716 CALL timestop(handle6)
3719 IF (.NOT. line_search)
THEN
3721 CALL timeset(
'k_opt_step', handle7)
3723 IF ((.NOT. md_in_k_space) .AND. &
3724 (iteration >= max(0, almo_scf_env%opt_k_prec_iter_start) .AND. &
3725 mod(iteration - almo_scf_env%opt_k_prec_iter_start, &
3726 almo_scf_env%opt_k_prec_iter_freq) == 0))
THEN
3731 IF (unit_nr > 0)
THEN
3732 WRITE (unit_nr, *)
"Computing preconditioner"
3734 CALL opt_k_create_preconditioner_blk(almo_scf_env, &
3735 almo_scf_env%matrix_v_disc(ispin), &
3747 CALL opt_k_apply_preconditioner_blk(almo_scf_env, &
3752 reset_conjugator = .false.
3754 IF (iteration < max(almo_scf_env%opt_k_conj_iter_start, 1) .OR. &
3755 mod(iteration - almo_scf_env%opt_k_conj_iter_start, &
3756 almo_scf_env%opt_k_conj_iter_freq) == 0)
THEN
3758 reset_conjugator = .true.
3763 CALL dbcsr_dot(grad, prev_minus_prec_grad, numer)
3764 CALL dbcsr_dot(prev_grad, prev_minus_prec_grad, denom)
3765 conjugacy_error = numer/denom
3767 IF (conjugacy_error > min(0.5_dp, conjugacy_error_threshold))
THEN
3768 reset_conjugator = .true.
3769 IF (unit_nr > 0)
THEN
3770 WRITE (unit_nr, *)
"Lack of progress, conjugacy error is ", conjugacy_error
3775 IF ((iteration /= 0) .AND. (.NOT. reset_conjugator))
THEN
3777 CALL dbcsr_dot(prev_grad, prev_step, denom)
3778 line_search_error = numer/denom
3779 IF (line_search_error > line_search_error_threshold)
THEN
3780 reset_conjugator = .true.
3781 IF (unit_nr > 0)
THEN
3782 WRITE (unit_nr, *)
"Bad line search, line search error is ", line_search_error
3790 IF (.NOT. reset_conjugator)
THEN
3792 SELECT CASE (almo_scf_env%opt_k_conjugator)
3795 CALL dbcsr_add(tmp_k_blk, prev_grad, 1.0_dp, -1.0_dp)
3797 CALL dbcsr_dot(tmp_k_blk, prev_step, denom)
3798 beta = -1.0_dp*numer/denom
3801 CALL dbcsr_dot(prev_grad, prev_minus_prec_grad, denom)
3804 CALL dbcsr_dot(prev_grad, prev_minus_prec_grad, denom)
3806 CALL dbcsr_add(tmp_k_blk, prev_grad, 1.0_dp, -1.0_dp)
3811 CALL dbcsr_dot(prev_grad, prev_step, denom)
3814 CALL dbcsr_dot(prev_grad, prev_step, denom)
3816 CALL dbcsr_add(tmp_k_blk, prev_grad, 1.0_dp, -1.0_dp)
3822 CALL dbcsr_add(tmp_k_blk, prev_grad, 1.0_dp, -1.0_dp)
3823 CALL dbcsr_dot(tmp_k_blk, prev_step, denom)
3824 beta = -1.0_dp*numer/denom
3827 CALL dbcsr_add(tmp_k_blk, prev_grad, 1.0_dp, -1.0_dp)
3828 CALL dbcsr_dot(tmp_k_blk, prev_step, denom)
3829 CALL dbcsr_dot(tmp_k_blk, prev_minus_prec_grad, numer)
3830 kappa = -2.0_dp*numer/denom
3832 tau = -1.0_dp*numer/denom
3834 beta = tau - kappa*numer/denom
3838 cpabort(
"illegal conjugator")
3841 IF (beta < 0.0_dp)
THEN
3842 IF (unit_nr > 0)
THEN
3843 WRITE (unit_nr, *)
"Beta is negative, ", beta
3845 reset_conjugator = .true.
3850 IF (md_in_k_space)
THEN
3851 reset_conjugator = .true.
3854 IF (reset_conjugator)
THEN
3858 IF (unit_nr > 0)
THEN
3859 WRITE (unit_nr, *)
"(Re)-setting conjugator to zero"
3868 CALL dbcsr_add(step, prev_step, 1.0_dp, beta)
3870 CALL timestop(handle7)
3874 conjugacy_error = 0.0_dp
3878 IF (line_search)
THEN
3880 line_search_error = gfun1/gfun0
3886 IF (line_search)
THEN
3889 safety_multiplier = 1.0e+1_dp
3890 num_threshold = max(epsilon(1.0_dp), &
3891 safety_multiplier*(almo_scf_env%eps_filter**2)*almo_scf_env%ndomains)
3892 IF (abs(fun1 - fun0 - gfun0*step_size) < num_threshold)
THEN
3893 IF (unit_nr > 0)
THEN
3894 WRITE (unit_nr,
'(T3,A,1X,E17.7)') &
3895 "Numerical accuracy is too low to observe non-linear behavior", &
3896 abs(fun1 - fun0 - gfun0*step_size)
3897 WRITE (unit_nr,
'(T3,A,1X,E17.7,A,1X,E12.3)')
"Error computing ", &
3899 " is smaller than the threshold", num_threshold
3901 cpabort(
"Unable to continue with low numerical accuracy")
3903 IF (abs(gfun0) < num_threshold)
THEN
3904 IF (unit_nr > 0)
THEN
3905 WRITE (unit_nr,
'(T3,A,1X,E17.7,A,1X,E12.3)')
"Linear gradient", &
3907 " is smaller than the threshold", num_threshold
3909 cpabort(
"Unable to continue with low numerical accuracy")
3912 use_quadratic_approximation = .true.
3913 use_cubic_approximation = .false.
3917 step_size_quadratic_approx = -(gfun0*step_size*step_size)/(2.0_dp*(fun1 - fun0 - gfun0*step_size))
3919 step_size_quadratic_approx2 = -(fun1 - fun0 - step_size*gfun1/2.0_dp)/(gfun1 - (fun1 - fun0)/step_size)
3921 IF ((step_size_quadratic_approx < 0.0_dp) .AND. &
3922 (step_size_quadratic_approx2 < 0.0_dp))
THEN
3923 IF (unit_nr > 0)
THEN
3924 WRITE (unit_nr,
'(T3,A,1X,E17.7,1X,E17.7,1X,A)') &
3925 "Quadratic approximation gives negative steps", &
3926 step_size_quadratic_approx, step_size_quadratic_approx2, &
3929 use_cubic_approximation = .true.
3930 use_quadratic_approximation = .false.
3932 IF (step_size_quadratic_approx < 0.0_dp)
THEN
3933 step_size_quadratic_approx = step_size_quadratic_approx2
3935 IF (step_size_quadratic_approx2 < 0.0_dp)
THEN
3936 step_size_quadratic_approx2 = step_size_quadratic_approx
3941 IF (use_quadratic_approximation)
THEN
3942 quadratic_approx_error = abs(step_size_quadratic_approx - &
3943 step_size_quadratic_approx2)/step_size_quadratic_approx
3944 IF (quadratic_approx_error > quadratic_approx_error_threshold)
THEN
3945 IF (unit_nr > 0)
THEN
3946 WRITE (unit_nr,
'(T3,A,1X,E17.7,1X,E17.7,1X,A)')
"Quadratic approximation is poor", &
3947 step_size_quadratic_approx, step_size_quadratic_approx2, &
3948 "Try cubic approximation"
3950 use_cubic_approximation = .true.
3951 use_quadratic_approximation = .false.
3956 IF (use_cubic_approximation)
THEN
3961 bb = (-step_size*gfun1 + 3.0_dp*(fun1 - fun0) - 2.0_dp*step_size*gfun0)/(step_size*step_size)
3962 aa = (gfun1 - 2.0_dp*step_size*bb - gfun0)/(3.0_dp*step_size*step_size)
3964 IF (abs(gfun1 - 2.0_dp*step_size*bb - gfun0) < num_threshold)
THEN
3965 IF (unit_nr > 0)
THEN
3966 WRITE (unit_nr,
'(T3,A,1X,E17.7)') &
3967 "Numerical accuracy is too low to observe cubic behavior", &
3968 abs(gfun1 - 2.0_dp*step_size*bb - gfun0)
3970 use_cubic_approximation = .false.
3971 use_quadratic_approximation = .true.
3973 IF (abs(gfun1) < num_threshold)
THEN
3974 IF (unit_nr > 0)
THEN
3975 WRITE (unit_nr,
'(T3,A,1X,E17.7,A,1X,E12.3)')
"Linear gradient", &
3977 " is smaller than the threshold", num_threshold
3979 use_cubic_approximation = .false.
3980 use_quadratic_approximation = .true.
3985 IF (use_cubic_approximation)
THEN
3990 IF (unit_nr > 0)
THEN
3991 WRITE (unit_nr,
'(T3,A)') &
3992 "Cubic approximation gives zero soultions! Use quadratic approximation"
3994 use_quadratic_approximation = .true.
3995 use_cubic_approximation = .true.
3997 step_size = minima(1)
3999 IF (unit_nr > 0)
THEN
4000 WRITE (unit_nr,
'(T3,A)') &
4001 "More than one solution found! Use quadratic approximation"
4003 use_quadratic_approximation = .true.
4004 use_cubic_approximation = .true.
4009 IF (use_quadratic_approximation)
THEN
4010 IF (unit_nr > 0)
THEN
4011 WRITE (unit_nr,
'(T3,A)')
"Use quadratic approximation"
4013 step_size = (step_size_quadratic_approx + step_size_quadratic_approx2)*0.5_dp
4017 IF (step_size < 0.0_dp)
THEN
4018 cpabort(
"Negative step proposed")
4021 CALL dbcsr_copy(almo_scf_env%matrix_k_blk(ispin), &
4023 CALL dbcsr_add(almo_scf_env%matrix_k_blk(ispin), &
4024 step, 1.0_dp, step_size)
4026 almo_scf_env%matrix_k_blk(ispin))
4027 line_search = .false.
4031 IF (md_in_k_space)
THEN
4034 IF (iteration /= 0)
THEN
4036 step, 1.0_dp, 0.5_dp*time_step)
4038 prev_step, 1.0_dp, 0.5_dp*time_step)
4041 kin_energy = 0.5_dp*kin_energy*kin_energy
4044 CALL dbcsr_add(almo_scf_env%matrix_k_blk(ispin), &
4045 velocity, 1.0_dp, time_step)
4046 CALL dbcsr_add(almo_scf_env%matrix_k_blk(ispin), &
4047 step, 1.0_dp, 0.5_dp*time_step*time_step)
4051 IF (reset_step_size)
THEN
4052 step_size = almo_scf_env%opt_k_trial_step_size
4053 reset_step_size = .false.
4055 step_size = step_size*almo_scf_env%opt_k_trial_step_size_multiplier
4057 CALL dbcsr_copy(almo_scf_env%matrix_k_blk(ispin), &
4059 CALL dbcsr_add(almo_scf_env%matrix_k_blk(ispin), &
4060 step, 1.0_dp, step_size)
4061 line_search = .true.
4070 IF (unit_nr > 0)
THEN
4071 IF (md_in_k_space)
THEN
4072 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)') &
4073 "K iter CG", iteration, time_step, time_step*iteration, &
4074 energy_correction(ispin), obj_function, delta_obj_function, grad_norm, &
4075 kin_energy, kin_energy + obj_function, beta
4077 IF (line_search .OR. prepare_to_exit)
THEN
4078 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)') &
4079 "K iter CG", iteration, step_size, &
4080 energy_correction(ispin), delta_obj_function, grad_norm, &
4081 gfun0, line_search_error, beta, conjugacy_error, t2a - t1a
4083 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)') &
4084 "K iter LS", iteration, step_size, &
4085 energy_correction(ispin), delta_obj_function, grad_norm, &
4086 gfun1, line_search_error, beta, conjugacy_error, t2a - t1a
4094 prepare_to_exit = .true.
4097 IF (.NOT. line_search) iteration = iteration + 1
4099 IF (prepare_to_exit)
EXIT
4103 IF (converged .OR. (outer_opt_k_iteration >= outer_opt_k_max_iter))
THEN
4104 outer_opt_k_prepare_to_exit = .true.
4107 IF (almo_scf_env%deloc_truncate_virt /=
virt_full)
THEN
4109 IF (unit_nr > 0)
THEN
4110 WRITE (unit_nr, *)
"Updating ALMO virtuals"
4113 CALL timeset(
'k_opt_v0_update', handle8)
4117 almo_scf_env%matrix_v_disc_blk(ispin), &
4118 almo_scf_env%matrix_k_blk(ispin), &
4120 filter_eps=almo_scf_env%eps_filter)
4121 CALL dbcsr_add(vr_fixed, almo_scf_env%matrix_v_blk(ispin), &
4126 almo_scf_env%matrix_v_blk(ispin), &
4127 almo_scf_env%matrix_k_blk(ispin), &
4129 filter_eps=almo_scf_env%eps_filter)
4130 CALL dbcsr_add(vd_fixed, almo_scf_env%matrix_v_disc_blk(ispin), &
4136 overlap=k_vr_index_down, &
4137 metric=almo_scf_env%matrix_s_blk(1), &
4138 retain_overlap_sparsity=.false., &
4139 eps_filter=almo_scf_env%eps_filter)
4140 CALL dbcsr_create(vr_index_sqrt_inv, template=k_vr_index_down, &
4141 matrix_type=dbcsr_type_no_symmetry)
4142 CALL dbcsr_create(vr_index_sqrt, template=k_vr_index_down, &
4143 matrix_type=dbcsr_type_no_symmetry)
4145 vr_index_sqrt_inv, &
4147 threshold=almo_scf_env%eps_filter, &
4148 order=almo_scf_env%order_lanczos, &
4149 eps_lanczos=almo_scf_env%eps_lanczos, &
4150 max_iter_lanczos=almo_scf_env%max_iter_lanczos)
4152 CALL dbcsr_create(matrix_tmp1, template=k_vr_index_down, &
4153 matrix_type=dbcsr_type_no_symmetry)
4154 CALL dbcsr_create(matrix_tmp2, template=k_vr_index_down, &
4155 matrix_type=dbcsr_type_no_symmetry)
4159 0.0_dp, matrix_tmp1, filter_eps=almo_scf_env%eps_filter)
4161 vr_index_sqrt_inv, &
4162 0.0_dp, matrix_tmp2, filter_eps=almo_scf_env%eps_filter)
4167 IF (unit_nr > 0)
THEN
4168 WRITE (unit_nr, *)
"Error for (inv(sqrt(SIGVV))*SIGVV*inv(sqrt(SIGVV))-I)", &
4169 frob_matrix/frob_matrix_base
4177 vr_index_sqrt_inv, &
4178 0.0_dp, almo_scf_env%matrix_v_blk(ispin), &
4179 filter_eps=almo_scf_env%eps_filter)
4183 overlap=k_vd_index_down, &
4184 metric=almo_scf_env%matrix_s_blk(1), &
4185 retain_overlap_sparsity=.false., &
4186 eps_filter=almo_scf_env%eps_filter)
4187 CALL dbcsr_create(vd_index_sqrt_inv, template=k_vd_index_down, &
4188 matrix_type=dbcsr_type_no_symmetry)
4189 CALL dbcsr_create(vd_index_sqrt, template=k_vd_index_down, &
4190 matrix_type=dbcsr_type_no_symmetry)
4192 vd_index_sqrt_inv, &
4194 threshold=almo_scf_env%eps_filter, &
4195 order=almo_scf_env%order_lanczos, &
4196 eps_lanczos=almo_scf_env%eps_lanczos, &
4197 max_iter_lanczos=almo_scf_env%max_iter_lanczos)
4199 CALL dbcsr_create(matrix_tmp1, template=k_vd_index_down, &
4200 matrix_type=dbcsr_type_no_symmetry)
4201 CALL dbcsr_create(matrix_tmp2, template=k_vd_index_down, &
4202 matrix_type=dbcsr_type_no_symmetry)
4206 0.0_dp, matrix_tmp1, filter_eps=almo_scf_env%eps_filter)
4208 vd_index_sqrt_inv, &
4209 0.0_dp, matrix_tmp2, filter_eps=almo_scf_env%eps_filter)
4214 IF (unit_nr > 0)
THEN
4215 WRITE (unit_nr, *)
"Error for (inv(sqrt(SIGVV))*SIGVV*inv(sqrt(SIGVV))-I)", &
4216 frob_matrix/frob_matrix_base
4224 vd_index_sqrt_inv, &
4225 0.0_dp, almo_scf_env%matrix_v_disc_blk(ispin), &
4226 filter_eps=almo_scf_env%eps_filter)
4233 CALL timestop(handle8)
4240 IF (almo_scf_env%deloc_truncate_virt /=
virt_full)
THEN
4260 IF (md_in_k_space)
THEN
4266 outer_opt_k_iteration = outer_opt_k_iteration + 1
4267 IF (outer_opt_k_prepare_to_exit)
EXIT
4281 IF (.NOT. almo_scf_env%s_sqrt_done)
THEN
4283 IF (unit_nr > 0)
THEN
4284 WRITE (unit_nr, *)
"sqrt and inv(sqrt) of AO overlap matrix"
4287 template=almo_scf_env%matrix_s(1), &
4288 matrix_type=dbcsr_type_no_symmetry)
4290 template=almo_scf_env%matrix_s(1), &
4291 matrix_type=dbcsr_type_no_symmetry)
4294 almo_scf_env%matrix_s_sqrt_inv(1), &
4295 almo_scf_env%matrix_s(1), &
4296 threshold=almo_scf_env%eps_filter, &
4297 order=almo_scf_env%order_lanczos, &
4298 eps_lanczos=almo_scf_env%eps_lanczos, &
4299 max_iter_lanczos=almo_scf_env%max_iter_lanczos)
4302 CALL dbcsr_create(matrix_tmp1, template=almo_scf_env%matrix_s(1), &
4303 matrix_type=dbcsr_type_no_symmetry)
4304 CALL dbcsr_create(matrix_tmp2, template=almo_scf_env%matrix_s(1), &
4305 matrix_type=dbcsr_type_no_symmetry)
4307 CALL dbcsr_multiply(
"N",
"N", 1.0_dp, almo_scf_env%matrix_s_sqrt_inv(1), &
4308 almo_scf_env%matrix_s(1), &
4309 0.0_dp, matrix_tmp1, filter_eps=almo_scf_env%eps_filter)
4310 CALL dbcsr_multiply(
"N",
"N", 1.0_dp, matrix_tmp1, almo_scf_env%matrix_s_sqrt_inv(1), &
4311 0.0_dp, matrix_tmp2, filter_eps=almo_scf_env%eps_filter)
4316 IF (unit_nr > 0)
THEN
4317 WRITE (unit_nr, *)
"Error for (inv(sqrt(S))*S*inv(sqrt(S))-I)", frob_matrix/frob_matrix_base
4324 almo_scf_env%s_sqrt_done = .true.
4332 para_env=almo_scf_env%para_env, &
4333 blacs_env=almo_scf_env%blacs_env, &
4334 use_occ_orbs=.true., &
4335 use_virt_orbs=almo_scf_env%deloc_cayley_use_virt_orbs, &
4336 occ_orbs_orthogonal=.false., &
4337 virt_orbs_orthogonal=almo_scf_env%orthogonal_basis, &
4338 tensor_type=almo_scf_env%deloc_cayley_tensor_type, &
4339 neglect_quadratic_term=almo_scf_env%deloc_cayley_linear, &
4340 calculate_energy_corr=.true., &
4343 pp_preconditioner_full=almo_scf_env%deloc_cayley_occ_precond, &
4344 qq_preconditioner_full=almo_scf_env%deloc_cayley_vir_precond, &
4345 eps_convergence=almo_scf_env%deloc_cayley_eps_convergence, &
4346 eps_filter=almo_scf_env%eps_filter, &
4348 q_index_up=almo_scf_env%matrix_s_sqrt_inv(1), &
4349 q_index_down=almo_scf_env%matrix_s_sqrt(1), &
4350 p_index_up=almo_scf_env%matrix_sigma_sqrt_inv(ispin), &
4351 p_index_down=almo_scf_env%matrix_sigma_sqrt(ispin), &
4352 matrix_ks=almo_scf_env%matrix_ks_0deloc(ispin), &
4353 matrix_p=almo_scf_env%matrix_p(ispin), &
4354 matrix_qp_template=almo_scf_env%matrix_t(ispin), &
4355 matrix_pq_template=almo_scf_env%matrix_t_tr(ispin), &
4356 matrix_t=almo_scf_env%matrix_t(ispin), &
4357 conjugator=almo_scf_env%deloc_cayley_conjugator, &
4358 max_iter=almo_scf_env%deloc_cayley_max_iter)
4366 energy_correction=energy_correction(ispin))
4372 energy_correction(1) = energy_correction(1)*spin_factor
4379 IF (unit_nr > 0)
THEN
4381 WRITE (unit_nr,
'(T2,A,I6,F20.9)')
"ECORR", ispin, &
4382 energy_correction(ispin)
4385 energy_correction_final = energy_correction_final + energy_correction(ispin)
4390 p=almo_scf_env%matrix_p(ispin), &
4391 eps_filter=almo_scf_env%eps_filter, &
4392 orthog_orbs=.false., &
4393 nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin), &
4394 s=almo_scf_env%matrix_s(1), &
4395 sigma=almo_scf_env%matrix_sigma(ispin), &
4396 sigma_inv=almo_scf_env%matrix_sigma_inv(ispin), &
4398 algorithm=almo_scf_env%sigma_inv_algorithm, &
4399 inverse_accelerator=almo_scf_env%order_lanczos, &
4400 inv_eps_factor=almo_scf_env%matrix_iter_eps_error_factor, &
4401 eps_lanczos=almo_scf_env%eps_lanczos, &
4402 max_iter_lanczos=almo_scf_env%max_iter_lanczos, &
4403 para_env=almo_scf_env%para_env, &
4404 blacs_env=almo_scf_env%blacs_env)
4406 IF (almo_scf_env%nspins == 1)
THEN
4416 IF (.NOT. almo_scf_env%s_inv_done)
THEN
4417 IF (unit_nr > 0)
THEN
4418 WRITE (unit_nr, *)
"Inverting AO overlap matrix"
4421 template=almo_scf_env%matrix_s(1), &
4422 matrix_type=dbcsr_type_no_symmetry)
4423 IF (.NOT. almo_scf_env%s_sqrt_done)
THEN
4425 almo_scf_env%matrix_s(1), &
4426 threshold=almo_scf_env%eps_filter)
4428 CALL dbcsr_multiply(
"N",
"N", 1.0_dp, almo_scf_env%matrix_s_sqrt_inv(1), &
4429 almo_scf_env%matrix_s_sqrt_inv(1), &
4430 0.0_dp, almo_scf_env%matrix_s_inv(1), &
4431 filter_eps=almo_scf_env%eps_filter)
4435 CALL dbcsr_create(matrix_tmp1, template=almo_scf_env%matrix_s(1), &
4436 matrix_type=dbcsr_type_no_symmetry)
4437 CALL dbcsr_multiply(
"N",
"N", 1.0_dp, almo_scf_env%matrix_s_inv(1), &
4438 almo_scf_env%matrix_s(1), &
4439 0.0_dp, matrix_tmp1, &
4440 filter_eps=almo_scf_env%eps_filter)
4444 IF (unit_nr > 0)
THEN
4445 WRITE (unit_nr, *)
"Error for (inv(S)*S-I)", &
4446 frob_matrix/frob_matrix_base
4451 almo_scf_env%s_inv_done = .true.
4455 ALLOCATE (matrix_p_almo_scf_converged(nspin))
4457 CALL dbcsr_create(matrix_p_almo_scf_converged(ispin), &
4458 template=almo_scf_env%matrix_p(ispin))
4459 CALL dbcsr_copy(matrix_p_almo_scf_converged(ispin), &
4460 almo_scf_env%matrix_p(ispin))
4466 nelectron_spin_real(1) = almo_scf_env%nelectrons_spin(ispin)
4467 IF (almo_scf_env%nspins == 1)
THEN
4468 nelectron_spin_real(1) = nelectron_spin_real(1)/2
4471 local_mu(1) = sum(almo_scf_env%mu_of_domain(:, ispin))/almo_scf_env%ndomains
4474 cpabort(
"CVS only: density_matrix_sign has not been updated in SVN")
4476 IF (almo_scf_env%nspins == 1)
THEN
4480 CALL dbcsr_add(matrix_p_almo_scf_converged(ispin), &
4481 almo_scf_env%matrix_p(ispin), -1.0_dp, 1.0_dp)
4482 CALL dbcsr_dot(almo_scf_env%matrix_ks_0deloc(ispin), &
4483 matrix_p_almo_scf_converged(ispin), &
4484 energy_correction(ispin))
4486 energy_correction_final = energy_correction_final + energy_correction(ispin)
4488 IF (unit_nr > 0)
THEN
4490 WRITE (unit_nr,
'(T2,A,I6,F20.9)')
"ECORR", ispin, &
4491 energy_correction(ispin)
4500 DEALLOCATE (matrix_p_almo_scf_converged)
4506 IF (unit_nr > 0)
THEN
4508 WRITE (unit_nr,
'(T2,A,F18.9,F18.9,F18.9,F12.6)')
"ETOT", &
4509 almo_scf_env%almo_scf_energy, &
4510 energy_correction_final, &
4511 almo_scf_env%almo_scf_energy + energy_correction_final, &
4516 CALL timestop(handle)
4518 END SUBROUTINE harris_foulkes_correction
4524 SUBROUTINE make_triu(matrix)
4527 CHARACTER(len=*),
PARAMETER :: routinen =
'make_triu'
4529 INTEGER :: col, handle, i, j, row
4530 REAL(
dp),
DIMENSION(:, :),
POINTER :: block
4533 CALL timeset(routinen, handle)
4538 IF (row > col) block(:, :) = 0.0_dp
4539 IF (row == col)
THEN
4540 DO j = 1,
SIZE(block, 2)
4541 DO i = j + 1,
SIZE(block, 1)
4542 block(i, j) = 0.0_dp
4550 CALL timestop(handle)
4551 END SUBROUTINE make_triu
4573 SUBROUTINE opt_k_create_preconditioner(prec, vd_prop, f, x, oo_inv_x_tr, s, grad, &
4574 vd_blk, t, template_vd_vd_blk, template_vr_vr_blk, template_n_vr, &
4575 spin_factor, eps_filter)
4578 TYPE(
dbcsr_type),
INTENT(IN) :: vd_prop, f, x, oo_inv_x_tr, s
4580 TYPE(
dbcsr_type),
INTENT(IN) :: vd_blk, t, template_vd_vd_blk, &
4581 template_vr_vr_blk, template_n_vr
4582 REAL(kind=
dp),
INTENT(IN) :: spin_factor, eps_filter
4584 CHARACTER(len=*),
PARAMETER :: routinen =
'opt_k_create_preconditioner'
4586 INTEGER :: handle, p_nrows, q_nrows
4587 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: p_diagonal, q_diagonal
4588 TYPE(
dbcsr_type) :: pp_diag, qq_diag, t1, t2, tmp, &
4589 tmp1_n_vr, tmp2_n_vr, tmp_n_vd, &
4590 tmp_vd_vd_blk, tmp_vr_vr_blk
4592 CALL timeset(routinen, handle)
4603 0.0_dp, tmp_n_vd, filter_eps=eps_filter)
4605 template=template_vd_vd_blk)
4606 CALL dbcsr_copy(tmp_vd_vd_blk, template_vd_vd_blk)
4608 0.0_dp, tmp_vd_vd_blk, &
4609 retain_sparsity=.true., &
4610 filter_eps=eps_filter)
4613 ALLOCATE (q_diagonal(q_nrows))
4616 template=template_vd_vd_blk)
4621 0.0_dp, t1, filter_eps=eps_filter)
4624 CALL dbcsr_create(tmp_vr_vr_blk, template=template_vr_vr_blk)
4625 CALL dbcsr_copy(tmp_vr_vr_blk, template_vr_vr_blk)
4627 0.0_dp, tmp_vr_vr_blk, &
4628 retain_sparsity=.true., &
4629 filter_eps=eps_filter)
4632 ALLOCATE (p_diagonal(p_nrows))
4634 CALL dbcsr_create(pp_diag, template=template_vr_vr_blk)
4640 0.0_dp, t2, filter_eps=eps_filter)
4646 0.0_dp, tmp_n_vd, filter_eps=eps_filter)
4648 0.0_dp, tmp_vd_vd_blk, &
4649 retain_sparsity=.true., &
4650 filter_eps=eps_filter)
4657 0.0_dp, t1, filter_eps=eps_filter)
4663 0.0_dp, tmp1_n_vr, filter_eps=eps_filter)
4665 0.0_dp, tmp2_n_vr, filter_eps=eps_filter)
4667 0.0_dp, tmp_vr_vr_blk, &
4668 retain_sparsity=.true., &
4669 filter_eps=eps_filter)
4676 0.0_dp, t2, filter_eps=eps_filter)
4679 CALL dbcsr_add(prec, tmp, 1.0_dp, -1.0_dp)
4684 0.0_dp, tmp_n_vd, filter_eps=eps_filter)
4686 0.0_dp, tmp, retain_sparsity=.true., &
4687 filter_eps=eps_filter)
4692 CALL dbcsr_add(prec, t1, 1.0_dp, 1.0_dp)
4694 CALL inverse_of_elements(prec)
4697 DEALLOCATE (q_diagonal)
4698 DEALLOCATE (p_diagonal)
4710 CALL timestop(handle)
4712 END SUBROUTINE opt_k_create_preconditioner
4727 SUBROUTINE opt_k_create_preconditioner_blk(almo_scf_env, vd_prop, oo_inv_x_tr, &
4728 t_curr, ispin, spin_factor)
4731 TYPE(
dbcsr_type),
INTENT(IN) :: vd_prop, oo_inv_x_tr, t_curr
4732 INTEGER,
INTENT(IN) :: ispin
4733 REAL(kind=
dp),
INTENT(IN) :: spin_factor
4735 CHARACTER(len=*),
PARAMETER :: routinen =
'opt_k_create_preconditioner_blk'
4738 REAL(kind=
dp) :: eps_filter
4739 TYPE(
dbcsr_type) :: opt_k_e_dd, opt_k_e_rr, s_dd_sqrt, &
4740 s_rr_sqrt, t1, tmp, tmp1_n_vr, &
4741 tmp2_n_vr, tmp_n_vd, tmp_vd_vd_blk, &
4746 CALL timeset(routinen, handle)
4748 eps_filter = almo_scf_env%eps_filter
4751 CALL dbcsr_create(tmp_n_vd, template=almo_scf_env%matrix_v_disc(ispin))
4753 template=almo_scf_env%matrix_vv_disc_blk(ispin), &
4754 matrix_type=dbcsr_type_no_symmetry)
4756 almo_scf_env%matrix_s(1), &
4758 0.0_dp, tmp_n_vd, filter_eps=eps_filter)
4760 almo_scf_env%matrix_vv_disc_blk(ispin))
4762 0.0_dp, tmp_vd_vd_blk, &
4763 retain_sparsity=.true.)
4766 template=almo_scf_env%matrix_vv_disc_blk(ispin), &
4767 matrix_type=dbcsr_type_no_symmetry)
4769 almo_scf_env%opt_k_t_dd(ispin), &
4771 threshold=eps_filter, &
4772 order=almo_scf_env%order_lanczos, &
4773 eps_lanczos=almo_scf_env%eps_lanczos, &
4774 max_iter_lanczos=almo_scf_env%max_iter_lanczos)
4778 almo_scf_env%matrix_ks_0deloc(ispin), &
4780 0.0_dp, tmp_n_vd, filter_eps=eps_filter)
4782 almo_scf_env%matrix_vv_disc_blk(ispin))
4784 0.0_dp, tmp_vd_vd_blk, &
4785 retain_sparsity=.true.)
4791 almo_scf_env%opt_k_t_dd(ispin), &
4792 0.0_dp, s_dd_sqrt, filter_eps=eps_filter)
4794 almo_scf_env%opt_k_t_dd(ispin), &
4796 0.0_dp, tmp_vd_vd_blk, filter_eps=eps_filter)
4800 template=almo_scf_env%matrix_vv_disc_blk(ispin))
4803 template=almo_scf_env%matrix_vv_disc_blk(ispin), &
4804 matrix_type=dbcsr_type_no_symmetry)
4812 almo_scf_env%opt_k_t_dd(ispin))
4816 0.0_dp, almo_scf_env%opt_k_t_dd(ispin), &
4817 filter_eps=eps_filter)
4823 template=almo_scf_env%matrix_k_blk_ones(ispin))
4825 almo_scf_env%matrix_k_blk_ones(ispin))
4827 template=almo_scf_env%matrix_k_blk_ones(ispin))
4830 0.0_dp, t1, filter_eps=eps_filter)
4835 template=almo_scf_env%matrix_sigma_vv_blk(ispin), &
4836 matrix_type=dbcsr_type_no_symmetry)
4838 almo_scf_env%matrix_sigma_vv_blk(ispin))
4840 almo_scf_env%matrix_x(ispin), &
4842 0.0_dp, tmp_vr_vr_blk, &
4843 retain_sparsity=.true.)
4847 template=almo_scf_env%matrix_sigma_vv_blk(ispin), &
4848 matrix_type=dbcsr_type_no_symmetry)
4850 almo_scf_env%opt_k_t_rr(ispin), &
4852 threshold=eps_filter, &
4853 order=almo_scf_env%order_lanczos, &
4854 eps_lanczos=almo_scf_env%eps_lanczos, &
4855 max_iter_lanczos=almo_scf_env%max_iter_lanczos)
4859 template=almo_scf_env%matrix_v(ispin))
4861 template=almo_scf_env%matrix_v(ispin))
4863 0.0_dp, tmp1_n_vr, filter_eps=eps_filter)
4865 almo_scf_env%matrix_ks_0deloc(ispin), &
4867 0.0_dp, tmp2_n_vr, filter_eps=eps_filter)
4869 0.0_dp, tmp_vr_vr_blk, &
4870 retain_sparsity=.true.)
4877 almo_scf_env%opt_k_t_rr(ispin), &
4878 0.0_dp, s_rr_sqrt, filter_eps=eps_filter)
4880 almo_scf_env%opt_k_t_rr(ispin), &
4882 0.0_dp, tmp_vr_vr_blk, filter_eps=eps_filter)
4886 template=almo_scf_env%matrix_sigma_vv_blk(ispin))
4889 template=almo_scf_env%matrix_sigma_vv_blk(ispin), &
4890 matrix_type=dbcsr_type_no_symmetry)
4898 almo_scf_env%opt_k_t_rr(ispin))
4902 0.0_dp, almo_scf_env%opt_k_t_rr(ispin), &
4903 filter_eps=eps_filter)
4910 0.0_dp, almo_scf_env%opt_k_denom(ispin), &
4911 filter_eps=eps_filter)
4916 CALL dbcsr_add(almo_scf_env%opt_k_denom(ispin), t1, &
4919 CALL dbcsr_scale(almo_scf_env%opt_k_denom(ispin), &
4922 CALL inverse_of_elements(almo_scf_env%opt_k_denom(ispin))
4926 CALL timestop(handle)
4928 END SUBROUTINE opt_k_create_preconditioner_blk
4942 SUBROUTINE opt_k_apply_preconditioner_blk(almo_scf_env, step, grad, ispin)
4947 INTEGER,
INTENT(IN) :: ispin
4949 CHARACTER(len=*),
PARAMETER :: routinen =
'opt_k_apply_preconditioner_blk'
4952 REAL(kind=
dp) :: eps_filter
4955 CALL timeset(routinen, handle)
4957 eps_filter = almo_scf_env%eps_filter
4959 CALL dbcsr_create(tmp_k, template=almo_scf_env%matrix_k_blk(ispin))
4963 grad, almo_scf_env%opt_k_t_rr(ispin), &
4964 0.0_dp, tmp_k, filter_eps=eps_filter)
4966 almo_scf_env%opt_k_t_dd(ispin), tmp_k, &
4967 0.0_dp, step, filter_eps=eps_filter)
4971 almo_scf_env%opt_k_denom(ispin), tmp_k)
4975 almo_scf_env%opt_k_t_dd(ispin), tmp_k, &
4976 0.0_dp, step, filter_eps=eps_filter)
4978 step, almo_scf_env%opt_k_t_rr(ispin), &
4979 0.0_dp, tmp_k, filter_eps=eps_filter)
4985 CALL timestop(handle)
4987 END SUBROUTINE opt_k_apply_preconditioner_blk
5026 SUBROUTINE compute_gradient(m_grad_out, m_ks, m_s, m_t, m_t0, &
5027 m_siginv, m_quench_t, m_FTsiginv, m_siginvTFTsiginv, m_ST, m_STsiginv0, &
5028 m_theta, domain_s_inv, domain_r_down, &
5029 cpu_of_domain, domain_map, assume_t0_q0x, optimize_theta, &
5030 normalize_orbitals, penalty_occ_vol, penalty_occ_local, &
5031 penalty_occ_vol_prefactor, envelope_amplitude, eps_filter, spin_factor, &
5032 special_case, m_sig_sqrti_ii, op_sm_set, weights, energy_coeff, &
5035 TYPE(
dbcsr_type),
INTENT(INOUT) :: m_grad_out, m_ks, m_s, m_t, m_t0, &
5036 m_siginv, m_quench_t, m_ftsiginv, &
5037 m_siginvtftsiginv, m_st, m_stsiginv0, &
5040 INTENT(IN) :: domain_s_inv, domain_r_down
5041 INTEGER,
DIMENSION(:),
INTENT(IN) :: cpu_of_domain
5043 LOGICAL,
INTENT(IN) :: assume_t0_q0x, optimize_theta, &
5044 normalize_orbitals, penalty_occ_vol
5045 LOGICAL,
INTENT(IN),
OPTIONAL :: penalty_occ_local
5046 REAL(kind=
dp),
INTENT(IN) :: penalty_occ_vol_prefactor, &
5047 envelope_amplitude, eps_filter, &
5049 INTEGER,
INTENT(IN) :: special_case
5050 TYPE(
dbcsr_type),
INTENT(IN),
OPTIONAL :: m_sig_sqrti_ii
5052 POINTER :: op_sm_set
5053 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN),
OPTIONAL :: weights
5054 REAL(kind=
dp),
INTENT(IN),
OPTIONAL :: energy_coeff, localiz_coeff
5056 CHARACTER(len=*),
PARAMETER :: routinen =
'compute_gradient'
5058 INTEGER :: dim0, handle, idim0, nao, reim
5059 LOGICAL :: my_penalty_local
5060 REAL(kind=
dp) :: coeff, energy_g_norm, my_energy_coeff, &
5062 penalty_occ_vol_g_norm
5063 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: tg_diagonal
5064 TYPE(
dbcsr_type) :: m_tmp_no_1, m_tmp_no_2, m_tmp_no_3, &
5065 m_tmp_oo_1, m_tmp_oo_2, temp1, temp2, &
5066 tempnocc1, tempoccocc1
5068 CALL timeset(routinen, handle)
5070 IF (normalize_orbitals .AND. (.NOT.
PRESENT(m_sig_sqrti_ii)))
THEN
5071 cpabort(
"Normalization matrix is required")
5074 my_penalty_local = .false.
5075 my_localiz_coeff = 1.0_dp
5076 my_energy_coeff = 0.0_dp
5077 IF (
PRESENT(localiz_coeff))
THEN
5078 my_localiz_coeff = localiz_coeff
5080 IF (
PRESENT(energy_coeff))
THEN
5081 my_energy_coeff = energy_coeff
5083 IF (
PRESENT(penalty_occ_local))
THEN
5084 my_penalty_local = penalty_occ_local
5093 template=m_quench_t, &
5094 matrix_type=dbcsr_type_no_symmetry)
5096 template=m_quench_t, &
5097 matrix_type=dbcsr_type_no_symmetry)
5099 template=m_quench_t, &
5100 matrix_type=dbcsr_type_no_symmetry)
5102 template=m_siginv, &
5103 matrix_type=dbcsr_type_no_symmetry)
5105 template=m_siginv, &
5106 matrix_type=dbcsr_type_no_symmetry)
5109 matrix_type=dbcsr_type_no_symmetry)
5111 template=m_siginv, &
5112 matrix_type=dbcsr_type_no_symmetry)
5115 matrix_type=dbcsr_type_no_symmetry)
5118 matrix_type=dbcsr_type_no_symmetry)
5121 CALL dbcsr_copy(m_tmp_no_2, m_ftsiginv, keep_sparsity=.true.)
5125 m_siginvtftsiginv, &
5126 1.0_dp, m_tmp_no_2, &
5127 retain_sparsity=.true.)
5131 IF (my_penalty_local)
THEN
5135 DO idim0 = 1,
SIZE(op_sm_set, 2)
5137 DO reim = 1,
SIZE(op_sm_set, 1)
5140 op_sm_set(reim, idim0)%matrix, &
5142 0.0_dp, tempnocc1, &
5143 filter_eps=eps_filter)
5149 0.0_dp, tempoccocc1, &
5150 filter_eps=eps_filter)
5153 ALLOCATE (tg_diagonal(dim0))
5157 DEALLOCATE (tg_diagonal)
5163 filter_eps=eps_filter)
5169 cpabort(
"Localization function is not implemented")
5171 coeff = -weights(idim0)
5173 cpabort(
"Localization function is not implemented")
5175 CALL dbcsr_add(temp2, temp1, 1.0_dp, coeff)
5178 CALL dbcsr_add(m_tmp_no_2, temp2, my_energy_coeff, my_localiz_coeff*4.0_dp)
5182 IF (penalty_occ_vol)
THEN
5185 penalty_occ_vol_prefactor, &
5188 0.0_dp, m_tmp_no_1, &
5189 retain_sparsity=.true.)
5193 CALL dbcsr_add(m_tmp_no_2, m_tmp_no_1, 1.0_dp, 1.0_dp)
5197 IF (normalize_orbitals)
THEN
5211 0.0_dp, m_tmp_no_1, &
5212 retain_sparsity=.true.)
5219 0.0_dp, m_tmp_oo_1, &
5220 retain_sparsity=.true.)
5223 ALLOCATE (tg_diagonal(dim0))
5227 DEALLOCATE (tg_diagonal)
5232 0.0_dp, m_tmp_oo_2, &
5233 filter_eps=eps_filter)
5237 1.0_dp, m_tmp_no_1, &
5238 retain_sparsity=.true.)
5247 IF (assume_t0_q0x)
THEN
5253 0.0_dp, m_tmp_oo_1, &
5254 filter_eps=eps_filter)
5258 1.0_dp, m_grad_out, &
5259 filter_eps=eps_filter)
5261 cpabort(
"Cannot project the zero-order space from itself")
5265 matrix_in=m_tmp_no_1, &
5266 matrix_out=m_grad_out, &
5267 operator2=domain_r_down(:), &
5268 operator1=domain_s_inv(:), &
5269 dpattern=m_quench_t, &
5271 node_of_domain=cpu_of_domain, &
5273 filter_eps=eps_filter, &
5275 use_trimmer=.false.)
5281 IF (optimize_theta)
THEN
5283 CALL dtanh_of_elements(m_tmp_no_2, alpha=1.0_dp/envelope_amplitude)
5310 CALL timestop(handle)
5312 END SUBROUTINE compute_gradient
5322 SUBROUTINE print_mathematica_matrix(matrix, filename)
5325 CHARACTER(len=*),
INTENT(IN) :: filename
5327 CHARACTER(len=*),
PARAMETER :: routinen =
'print_mathematica_matrix'
5329 CHARACTER(LEN=20) :: formatstr, scols
5330 INTEGER :: col, fiunit, handle, hori_offset, jj, &
5331 nblkcols_tot, nblkrows_tot, ncols, &
5332 ncores, nrows, row, unit_nr, &
5334 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: ao_block_sizes, mo_block_sizes
5335 INTEGER,
DIMENSION(:),
POINTER :: ao_blk_sizes, mo_blk_sizes
5337 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: h
5338 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: block_p
5343 CALL timeset(routinen, handle)
5347 IF (logger%para_env%is_source())
THEN
5356 IF (ncores > 1)
THEN
5357 cpabort(
"mathematica files: serial code only")
5360 CALL dbcsr_get_info(matrix, row_blk_size=ao_blk_sizes, col_blk_size=mo_blk_sizes, &
5361 nblkrows_total=nblkrows_tot, nblkcols_total=nblkcols_tot)
5362 cpassert(nblkrows_tot == nblkcols_tot)
5363 ALLOCATE (mo_block_sizes(nblkcols_tot), ao_block_sizes(nblkcols_tot))
5364 mo_block_sizes(:) = mo_blk_sizes(:)
5365 ao_block_sizes(:) = ao_blk_sizes(:)
5369 matrix_type=dbcsr_type_no_symmetry)
5372 ncols = sum(mo_block_sizes)
5373 nrows = sum(ao_block_sizes)
5374 ALLOCATE (h(nrows, ncols))
5378 DO col = 1, nblkcols_tot
5381 DO row = 1, nblkrows_tot
5386 h(vert_offset + 1:vert_offset + ao_block_sizes(row), &
5387 hori_offset + 1:hori_offset + mo_block_sizes(col)) &
5392 vert_offset = vert_offset + ao_block_sizes(row)
5396 hori_offset = hori_offset + mo_block_sizes(col)
5402 IF (unit_nr > 0)
THEN
5403 CALL open_file(filename, unit_number=fiunit, file_status=
'REPLACE')
5404 WRITE (scols,
"(I10)") ncols
5405 formatstr =
"("//trim(scols)//
"E27.17)"
5407 WRITE (fiunit, formatstr) h(jj, :)
5412 DEALLOCATE (mo_block_sizes)
5413 DEALLOCATE (ao_block_sizes)
5416 CALL timestop(handle)
5418 END SUBROUTINE print_mathematica_matrix
5440 SUBROUTINE compute_obj_nlmos(localization_obj_function_ispin, penalty_func_ispin, &
5441 penalty_vol_prefactor, overlap_determinant, m_sigma, nocc, m_B0, &
5442 m_theta_normalized, template_matrix_mo, weights, m_S0, just_started, &
5443 penalty_amplitude, eps_filter)
5445 REAL(kind=
dp),
INTENT(INOUT) :: localization_obj_function_ispin, penalty_func_ispin, &
5446 penalty_vol_prefactor, overlap_determinant
5448 INTEGER,
INTENT(IN) :: nocc
5449 TYPE(
dbcsr_type),
DIMENSION(:, :),
INTENT(IN) :: m_b0
5450 TYPE(
dbcsr_type),
INTENT(IN) :: m_theta_normalized, template_matrix_mo
5451 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: weights
5453 LOGICAL,
INTENT(IN) :: just_started
5454 REAL(kind=
dp),
INTENT(IN) :: penalty_amplitude, eps_filter
5456 CHARACTER(len=*),
PARAMETER :: routinen =
'compute_obj_nlmos'
5458 INTEGER :: handle, idim0, ielem, reim
5459 REAL(kind=
dp) :: det1, fval
5460 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: reim_diag, z2
5461 TYPE(
dbcsr_type) :: tempnocc1, tempoccocc1, tempoccocc2
5464 CALL timeset(routinen, handle)
5467 template=template_matrix_mo, &
5468 matrix_type=dbcsr_type_no_symmetry)
5470 template=m_theta_normalized, &
5471 matrix_type=dbcsr_type_no_symmetry)
5473 template=m_theta_normalized, &
5474 matrix_type=dbcsr_type_no_symmetry)
5476 localization_obj_function_ispin = 0.0_dp
5477 penalty_func_ispin = 0.0_dp
5479 ALLOCATE (reim_diag(nocc))
5483 DO idim0 = 1,
SIZE(m_b0, 2)
5487 DO reim = 1,
SIZE(m_b0, 1)
5490 m_b0(reim, idim0), &
5491 m_theta_normalized, &
5492 0.0_dp, tempoccocc1, &
5493 filter_eps=eps_filter)
5497 m_theta_normalized, &
5499 0.0_dp, tempoccocc2, &
5500 retain_sparsity=.true.)
5504 CALL group%sum(reim_diag)
5505 z2(:) = z2(:) + reim_diag(:)*reim_diag(:)
5512 fval = -weights(idim0)*log(abs(z2(ielem)))
5514 fval = weights(idim0) - weights(idim0)*abs(z2(ielem))
5516 fval = weights(idim0) - weights(idim0)*sqrt(abs(z2(ielem)))
5518 localization_obj_function_ispin = localization_obj_function_ispin + fval
5524 DEALLOCATE (reim_diag)
5528 m_theta_normalized, &
5529 0.0_dp, tempoccocc1, &
5530 filter_eps=eps_filter)
5533 m_theta_normalized, &
5536 filter_eps=eps_filter)
5541 overlap_determinant = det1
5543 IF (just_started .AND. penalty_amplitude < 0.0_dp)
THEN
5544 penalty_vol_prefactor = -(-penalty_amplitude)*localization_obj_function_ispin
5546 penalty_func_ispin = penalty_func_ispin + penalty_vol_prefactor*log(det1)
5552 CALL timestop(handle)
5554 END SUBROUTINE compute_obj_nlmos
5572 SUBROUTINE compute_gradient_nlmos(m_grad_out, m_B0, weights, &
5573 m_S0, m_theta_normalized, m_siginv, m_sig_sqrti_ii, &
5574 penalty_vol_prefactor, eps_filter, suggested_vol_penalty)
5576 TYPE(
dbcsr_type),
INTENT(INOUT) :: m_grad_out
5577 TYPE(
dbcsr_type),
DIMENSION(:, :),
INTENT(IN) :: m_b0
5578 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: weights
5579 TYPE(
dbcsr_type),
INTENT(IN) :: m_s0, m_theta_normalized, m_siginv, &
5581 REAL(kind=
dp),
INTENT(IN) :: penalty_vol_prefactor, eps_filter
5582 REAL(kind=
dp),
INTENT(INOUT) :: suggested_vol_penalty
5584 CHARACTER(len=*),
PARAMETER :: routinen =
'compute_gradient_nlmos'
5586 INTEGER :: dim0, handle, idim0, reim
5587 REAL(kind=
dp) :: norm_loc, norm_vol
5588 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: tg_diagonal, z2
5589 TYPE(
dbcsr_type) :: m_temp_oo_1, m_temp_oo_2, m_temp_oo_3, &
5592 CALL timeset(routinen, handle)
5595 template=m_theta_normalized, &
5596 matrix_type=dbcsr_type_no_symmetry)
5598 template=m_theta_normalized, &
5599 matrix_type=dbcsr_type_no_symmetry)
5601 template=m_theta_normalized, &
5602 matrix_type=dbcsr_type_no_symmetry)
5604 template=m_theta_normalized, &
5605 matrix_type=dbcsr_type_no_symmetry)
5608 ALLOCATE (tg_diagonal(dim0))
5613 DO idim0 = 1,
SIZE(m_b0, 2)
5617 DO reim = 1,
SIZE(m_b0, 1)
5620 m_b0(reim, idim0), &
5621 m_theta_normalized, &
5622 0.0_dp, m_temp_oo_3, &
5623 filter_eps=eps_filter)
5628 m_theta_normalized, &
5630 0.0_dp, m_temp_oo_4, &
5631 filter_eps=eps_filter)
5633 tg_diagonal(:) = 0.0_dp
5637 z2(:) = z2(:) + tg_diagonal(:)*tg_diagonal(:)
5642 1.0_dp, m_temp_oo_2, &
5643 filter_eps=eps_filter)
5651 z2(:) = -weights(idim0)/z2(:)
5653 z2(:) = -weights(idim0)
5655 z2(:) = -weights(idim0)/(2*sqrt(z2(:)))
5665 1.0_dp, m_temp_oo_1, &
5666 filter_eps=eps_filter)
5675 m_theta_normalized, &
5676 0.0_dp, m_temp_oo_2, &
5677 filter_eps=eps_filter)
5685 0.0_dp, m_temp_oo_3, &
5686 filter_eps=eps_filter)
5689 suggested_vol_penalty = norm_loc/norm_vol
5690 CALL dbcsr_add(m_temp_oo_1, m_temp_oo_3, &
5691 1.0_dp, 2.0_dp*penalty_vol_prefactor)
5699 0.0_dp, m_grad_out, &
5700 filter_eps=eps_filter)
5705 m_theta_normalized, &
5707 0.0_dp, m_temp_oo_3, &
5708 filter_eps=eps_filter)
5718 0.0_dp, m_temp_oo_1, &
5719 filter_eps=eps_filter)
5724 1.0_dp, m_grad_out, &
5725 filter_eps=eps_filter)
5727 DEALLOCATE (tg_diagonal)
5733 CALL timestop(handle)
5735 END SUBROUTINE compute_gradient_nlmos
5766 SUBROUTINE compute_xalmos_from_main_var(m_var_in, m_t_out, m_quench_t, &
5767 m_t0, m_oo_template, m_STsiginv0, m_s, m_sig_sqrti_ii_out, domain_r_down, &
5768 domain_s_inv, domain_map, cpu_of_domain, assume_t0_q0x, just_started, &
5769 optimize_theta, normalize_orbitals, envelope_amplitude, eps_filter, &
5770 special_case, nocc_of_domain, order_lanczos, eps_lanczos, max_iter_lanczos)
5773 TYPE(
dbcsr_type),
INTENT(INOUT) :: m_t_out, m_quench_t, m_t0, &
5774 m_oo_template, m_stsiginv0, m_s, &
5777 INTENT(IN) :: domain_r_down, domain_s_inv
5779 INTEGER,
DIMENSION(:),
INTENT(IN) :: cpu_of_domain
5780 LOGICAL,
INTENT(IN) :: assume_t0_q0x, just_started, &
5781 optimize_theta, normalize_orbitals
5782 REAL(kind=
dp),
INTENT(IN) :: envelope_amplitude, eps_filter
5783 INTEGER,
INTENT(IN) :: special_case
5784 INTEGER,
DIMENSION(:),
INTENT(IN) :: nocc_of_domain
5785 INTEGER,
INTENT(IN) :: order_lanczos
5786 REAL(kind=
dp),
INTENT(IN) :: eps_lanczos
5787 INTEGER,
INTENT(IN) :: max_iter_lanczos
5789 CHARACTER(len=*),
PARAMETER :: routinen =
'compute_xalmos_from_main_var'
5791 INTEGER :: handle, unit_nr
5792 REAL(kind=
dp) :: t_norm
5796 CALL timeset(routinen, handle)
5800 IF (logger%para_env%is_source())
THEN
5807 template=m_quench_t, &
5808 matrix_type=dbcsr_type_no_symmetry)
5810 template=m_oo_template, &
5811 matrix_type=dbcsr_type_no_symmetry)
5814 IF (optimize_theta)
THEN
5818 IF (unit_nr > 0)
THEN
5819 WRITE (unit_nr, *)
"Maximum norm of the initial guess: ", t_norm
5820 WRITE (unit_nr, *)
"Maximum allowed amplitude: ", &
5823 IF (t_norm > envelope_amplitude .AND. just_started)
THEN
5824 cpabort(
"Max norm of the initial guess is too large")
5827 CALL tanh_of_elements(m_tmp_no_1, alpha=1.0_dp/envelope_amplitude)
5834 IF (assume_t0_q0x)
THEN
5839 0.0_dp, m_tmp_oo_1, &
5840 filter_eps=eps_filter)
5845 filter_eps=eps_filter)
5847 cpabort(
"cannot use projector with block-daigonal ALMOs")
5851 matrix_in=m_t_out, &
5852 matrix_out=m_tmp_no_1, &
5853 operator1=domain_r_down, &
5854 operator2=domain_s_inv, &
5855 dpattern=m_quench_t, &
5857 node_of_domain=cpu_of_domain, &
5859 filter_eps=eps_filter, &
5860 use_trimmer=.false.)
5865 m_t0, 1.0_dp, 1.0_dp)
5868 IF (normalize_orbitals)
THEN
5871 overlap=m_tmp_oo_1, &
5873 retain_locality=.true., &
5874 only_normalize=.true., &
5875 nocc_of_domain=nocc_of_domain(:), &
5876 eps_filter=eps_filter, &
5877 order_lanczos=order_lanczos, &
5878 eps_lanczos=eps_lanczos, &
5879 max_iter_lanczos=max_iter_lanczos, &
5880 overlap_sqrti=m_sig_sqrti_ii_out)
5888 CALL timestop(handle)
5890 END SUBROUTINE compute_xalmos_from_main_var
5928 SUBROUTINE compute_preconditioner(domain_prec_out, m_prec_out, m_ks, m_s, &
5929 m_siginv, m_quench_t, m_FTsiginv, m_siginvTFTsiginv, m_ST, &
5930 m_STsiginv_out, m_s_vv_out, m_f_vv_out, para_env, &
5931 blacs_env, nocc_of_domain, domain_s_inv, domain_s_inv_half, domain_s_half, &
5932 domain_r_down, cpu_of_domain, &
5933 domain_map, assume_t0_q0x, penalty_occ_vol, penalty_occ_vol_prefactor, &
5934 eps_filter, neg_thr, spin_factor, special_case, bad_modes_projector_down_out, &
5938 INTENT(INOUT) :: domain_prec_out
5939 TYPE(
dbcsr_type),
INTENT(INOUT) :: m_prec_out, m_ks, m_s
5940 TYPE(
dbcsr_type),
INTENT(IN) :: m_siginv, m_quench_t, m_ftsiginv, &
5941 m_siginvtftsiginv, m_st
5942 TYPE(
dbcsr_type),
INTENT(INOUT),
OPTIONAL :: m_stsiginv_out, m_s_vv_out, m_f_vv_out
5945 INTEGER,
DIMENSION(:),
INTENT(IN) :: nocc_of_domain
5947 INTENT(IN) :: domain_s_inv
5949 INTENT(IN),
OPTIONAL :: domain_s_inv_half, domain_s_half
5951 INTENT(IN) :: domain_r_down
5952 INTEGER,
DIMENSION(:),
INTENT(IN) :: cpu_of_domain
5954 LOGICAL,
INTENT(IN) :: assume_t0_q0x, penalty_occ_vol
5955 REAL(kind=
dp),
INTENT(IN) :: penalty_occ_vol_prefactor, eps_filter, &
5956 neg_thr, spin_factor
5957 INTEGER,
INTENT(IN) :: special_case
5959 INTENT(INOUT),
OPTIONAL :: bad_modes_projector_down_out
5960 LOGICAL,
INTENT(IN) :: skip_inversion
5962 CHARACTER(len=*),
PARAMETER :: routinen =
'compute_preconditioner'
5964 INTEGER :: handle, ndim, precond_domain_projector
5965 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: nn_diagonal
5968 CALL timeset(routinen, handle)
5972 matrix_type=dbcsr_type_no_symmetry)
5974 template=m_quench_t, &
5975 matrix_type=dbcsr_type_no_symmetry)
5987 0.0_dp, m_tmp_no_3, &
5988 filter_eps=eps_filter)
5991 IF (
PRESENT(m_stsiginv_out))
THEN
6000 1.0_dp, m_tmp_nn_1, &
6001 filter_eps=eps_filter)
6004 IF (
PRESENT(m_s_vv_out))
THEN
6015 1.0_dp, m_prec_out, &
6016 filter_eps=eps_filter)
6020 1.0_dp, m_prec_out, &
6021 filter_eps=eps_filter)
6024 m_siginvtftsiginv, &
6025 0.0_dp, m_tmp_no_3, &
6026 filter_eps=eps_filter)
6030 1.0_dp, m_prec_out, &
6031 filter_eps=eps_filter)
6033 IF (
PRESENT(m_f_vv_out))
THEN
6038 CALL dbcsr_add(m_prec_out, m_tmp_nn_1, &
6044 IF (penalty_occ_vol)
THEN
6045 CALL dbcsr_add(m_prec_out, m_tmp_nn_1, &
6046 1.0_dp, penalty_occ_vol_prefactor)
6054 IF (skip_inversion)
THEN
6058 ALLOCATE (nn_diagonal(ndim))
6063 DEALLOCATE (nn_diagonal)
6065 CALL dbcsr_copy(m_prec_out, m_tmp_nn_1, keep_sparsity=.true.)
6070 matrix_in=m_tmp_nn_1, &
6071 matrix_out=m_prec_out, &
6072 nocc=nocc_of_domain(:) &
6079 IF (skip_inversion)
THEN
6085 para_env=para_env, &
6086 blacs_env=blacs_env)
6088 para_env=para_env, &
6089 blacs_env=blacs_env, &
6090 uplo_to_full=.true.)
6098 IF (assume_t0_q0x)
THEN
6099 precond_domain_projector = -1
6101 precond_domain_projector = 0
6106 IF (
PRESENT(bad_modes_projector_down_out))
THEN
6108 matrix_main=m_tmp_nn_1, &
6109 subm_s_inv=domain_s_inv(:), &
6110 subm_s_inv_half=domain_s_inv_half(:), &
6111 subm_s_half=domain_s_half(:), &
6112 subm_r_down=domain_r_down(:), &
6113 matrix_trimmer=m_quench_t, &
6114 dpattern=m_quench_t, &
6116 node_of_domain=cpu_of_domain, &
6118 use_trimmer=.false., &
6119 bad_modes_projector_down=bad_modes_projector_down_out(:), &
6120 eps_zero_eigenvalues=neg_thr, &
6121 my_action=precond_domain_projector, &
6122 skip_inversion=skip_inversion &
6126 matrix_main=m_tmp_nn_1, &
6127 subm_s_inv=domain_s_inv(:), &
6128 subm_r_down=domain_r_down(:), &
6129 matrix_trimmer=m_quench_t, &
6130 dpattern=m_quench_t, &
6132 node_of_domain=cpu_of_domain, &
6134 use_trimmer=.false., &
6136 my_action=precond_domain_projector, &
6137 skip_inversion=skip_inversion &
6146 CALL timestop(handle)
6148 END SUBROUTINE compute_preconditioner
6166 SUBROUTINE compute_cg_beta(beta, numer, denom, reset_conjugator, conjugator, &
6167 grad, prev_grad, step, prev_step, prev_minus_prec_grad)
6169 REAL(kind=
dp),
INTENT(INOUT) :: beta
6170 REAL(kind=
dp),
INTENT(INOUT),
OPTIONAL :: numer, denom
6171 LOGICAL,
INTENT(INOUT) :: reset_conjugator
6172 INTEGER,
INTENT(IN) :: conjugator
6173 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(INOUT) :: grad, prev_grad, step, prev_step
6174 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(INOUT), &
6175 OPTIONAL :: prev_minus_prec_grad
6177 CHARACTER(len=*),
PARAMETER :: routinen =
'compute_cg_beta'
6179 INTEGER :: handle, i, nsize, unit_nr
6180 REAL(kind=
dp) :: den, kappa, my_denom, my_numer, &
6181 my_numer2, my_numer3, num, num2, num3, &
6186 CALL timeset(routinen, handle)
6190 IF (logger%para_env%is_source())
THEN
6196 IF (.NOT.
PRESENT(prev_minus_prec_grad))
THEN
6200 cpabort(
"conjugator needs more input")
6205 IF (
PRESENT(numer) .OR.
PRESENT(denom))
THEN
6209 cpabort(
"cannot return numer/denom")
6224 matrix_type=dbcsr_type_no_symmetry)
6226 SELECT CASE (conjugator)
6229 CALL dbcsr_add(m_tmp_no_1, prev_grad(i), &
6231 CALL dbcsr_dot(m_tmp_no_1, step(i), num)
6232 CALL dbcsr_dot(m_tmp_no_1, prev_step(i), den)
6235 CALL dbcsr_dot(prev_grad(i), prev_minus_prec_grad(i), den)
6237 CALL dbcsr_dot(prev_grad(i), prev_minus_prec_grad(i), den)
6239 CALL dbcsr_add(m_tmp_no_1, prev_grad(i), 1.0_dp, -1.0_dp)
6240 CALL dbcsr_dot(m_tmp_no_1, step(i), num)
6243 CALL dbcsr_dot(prev_grad(i), prev_step(i), den)
6245 CALL dbcsr_dot(prev_grad(i), prev_step(i), den)
6247 CALL dbcsr_add(m_tmp_no_1, prev_grad(i), 1.0_dp, -1.0_dp)
6248 CALL dbcsr_dot(m_tmp_no_1, step(i), num)
6252 CALL dbcsr_add(m_tmp_no_1, prev_grad(i), 1.0_dp, -1.0_dp)
6253 CALL dbcsr_dot(m_tmp_no_1, prev_step(i), den)
6256 CALL dbcsr_add(m_tmp_no_1, prev_grad(i), 1.0_dp, -1.0_dp)
6257 CALL dbcsr_dot(m_tmp_no_1, prev_step(i), den)
6258 CALL dbcsr_dot(m_tmp_no_1, prev_minus_prec_grad(i), num)
6259 CALL dbcsr_dot(m_tmp_no_1, step(i), num2)
6260 CALL dbcsr_dot(prev_step(i), grad(i), num3)
6261 my_numer2 = my_numer2 + num2
6262 my_numer3 = my_numer3 + num3
6267 cpabort(
"illegal conjugator")
6269 my_numer = my_numer + num
6270 my_denom = my_denom + den
6278 SELECT CASE (conjugator)
6280 beta = -1.0_dp*my_numer/my_denom
6282 beta = my_numer/my_denom
6284 kappa = -2.0_dp*my_numer/my_denom
6285 tau = -1.0_dp*my_numer2/my_denom
6286 beta = tau - kappa*my_numer3/my_denom
6290 cpabort(
"illegal conjugator")
6295 IF (beta < 0.0_dp)
THEN
6296 IF (unit_nr > 0)
THEN
6297 WRITE (unit_nr, *)
" Resetting conjugator because beta is negative: ", beta
6299 reset_conjugator = .true.
6302 IF (
PRESENT(numer))
THEN
6305 IF (
PRESENT(denom))
THEN
6309 CALL timestop(handle)
6311 END SUBROUTINE compute_cg_beta
6345 SUBROUTINE newton_grad_to_step(optimizer, m_grad, m_delta, m_s, m_ks, &
6346 m_siginv, m_quench_t, m_FTsiginv, m_siginvTFTsiginv, m_ST, m_t, &
6347 m_sig_sqrti_ii, domain_s_inv, domain_r_down, domain_map, cpu_of_domain, &
6348 nocc_of_domain, para_env, blacs_env, eps_filter, optimize_theta, &
6349 penalty_occ_vol, normalize_orbitals, penalty_occ_vol_prefactor, &
6350 penalty_occ_vol_pf2, special_case)
6353 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(IN) :: m_grad
6354 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(INOUT) :: m_delta, m_s, m_ks, m_siginv, m_quench_t
6355 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(IN) :: m_ftsiginv, m_siginvtftsiginv, m_st, &
6358 INTENT(IN) :: domain_s_inv, domain_r_down
6360 INTEGER,
DIMENSION(:),
INTENT(IN) :: cpu_of_domain
6361 INTEGER,
DIMENSION(:, :),
INTENT(IN) :: nocc_of_domain
6364 REAL(kind=
dp),
INTENT(IN) :: eps_filter
6365 LOGICAL,
INTENT(IN) :: optimize_theta, penalty_occ_vol, &
6367 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: penalty_occ_vol_prefactor, &
6369 INTEGER,
INTENT(IN) :: special_case
6371 CHARACTER(len=*),
PARAMETER :: routinen =
'newton_grad_to_step'
6373 CHARACTER(LEN=20) :: iter_type
6374 INTEGER :: handle, ispin, iteration, max_iter, &
6375 ndomains, nspins, outer_iteration, &
6376 outer_max_iter, unit_nr
6377 LOGICAL :: converged, do_exact_inversion, outer_prepare_to_exit, prepare_to_exit, &
6378 reset_conjugator, use_preconditioner
6379 REAL(kind=
dp) :: alpha, beta, denom, denom_ispin, &
6380 eps_error_target, numer, numer_ispin, &
6381 residue_norm, spin_factor, t1, t2
6382 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: residue_max_norm
6385 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:) :: m_f_vo, m_f_vv, m_hstep, m_prec, &
6386 m_residue, m_residue_prev, m_s_vv, &
6387 m_step, m_stsiginv, m_zet, m_zet_prev
6389 DIMENSION(:, :) :: domain_prec
6391 CALL timeset(routinen, handle)
6395 IF (logger%para_env%is_source())
THEN
6402 IF (optimize_theta)
THEN
6403 cpabort(
"theta is NYI")
6408 outer_max_iter = optimizer%max_iter_outer_loop
6409 max_iter = optimizer%max_iter
6410 eps_error_target = optimizer%eps_error
6414 ndomains =
SIZE(domain_s_inv, 1)
6416 IF (nspins == 1)
THEN
6417 spin_factor = 2.0_dp
6419 spin_factor = 1.0_dp
6422 ALLOCATE (domain_prec(ndomains, nspins))
6426 ALLOCATE (m_residue(nspins))
6427 ALLOCATE (m_residue_prev(nspins))
6428 ALLOCATE (m_step(nspins))
6429 ALLOCATE (m_zet(nspins))
6430 ALLOCATE (m_zet_prev(nspins))
6431 ALLOCATE (m_hstep(nspins))
6432 ALLOCATE (m_prec(nspins))
6433 ALLOCATE (m_s_vv(nspins))
6434 ALLOCATE (m_f_vv(nspins))
6435 ALLOCATE (m_f_vo(nspins))
6436 ALLOCATE (m_stsiginv(nspins))
6438 ALLOCATE (residue_max_norm(nspins))
6441 DO ispin = 1, nspins
6445 template=m_quench_t(ispin), &
6446 matrix_type=dbcsr_type_no_symmetry)
6448 template=m_quench_t(ispin), &
6449 matrix_type=dbcsr_type_no_symmetry)
6451 template=m_quench_t(ispin), &
6452 matrix_type=dbcsr_type_no_symmetry)
6454 template=m_quench_t(ispin), &
6455 matrix_type=dbcsr_type_no_symmetry)
6457 template=m_quench_t(ispin), &
6458 matrix_type=dbcsr_type_no_symmetry)
6460 template=m_quench_t(ispin), &
6461 matrix_type=dbcsr_type_no_symmetry)
6463 template=m_quench_t(ispin), &
6464 matrix_type=dbcsr_type_no_symmetry)
6466 template=m_quench_t(ispin), &
6467 matrix_type=dbcsr_type_no_symmetry)
6469 template=m_ks(ispin), &
6470 matrix_type=dbcsr_type_no_symmetry)
6473 matrix_type=dbcsr_type_no_symmetry)
6475 template=m_ks(ispin), &
6476 matrix_type=dbcsr_type_no_symmetry)
6480 CALL dbcsr_copy(m_f_vo(ispin), m_ftsiginv(ispin))
6483 m_siginvtftsiginv(ispin), &
6484 1.0_dp, m_f_vo(ispin), &
6485 filter_eps=eps_filter)
6490 CALL compute_preconditioner( &
6491 domain_prec_out=domain_prec(:, ispin), &
6492 m_prec_out=m_prec(ispin), &
6495 m_siginv=m_siginv(ispin), &
6496 m_quench_t=m_quench_t(ispin), &
6497 m_ftsiginv=m_ftsiginv(ispin), &
6498 m_siginvtftsiginv=m_siginvtftsiginv(ispin), &
6500 m_stsiginv_out=m_stsiginv(ispin), &
6501 m_s_vv_out=m_s_vv(ispin), &
6502 m_f_vv_out=m_f_vv(ispin), &
6503 para_env=para_env, &
6504 blacs_env=blacs_env, &
6505 nocc_of_domain=nocc_of_domain(:, ispin), &
6506 domain_s_inv=domain_s_inv(:, ispin), &
6507 domain_r_down=domain_r_down(:, ispin), &
6508 cpu_of_domain=cpu_of_domain(:), &
6509 domain_map=domain_map(ispin), &
6510 assume_t0_q0x=.false., &
6511 penalty_occ_vol=penalty_occ_vol, &
6512 penalty_occ_vol_prefactor=penalty_occ_vol_prefactor(ispin), &
6513 eps_filter=eps_filter, &
6515 spin_factor=spin_factor, &
6516 special_case=special_case, &
6517 skip_inversion=.false. &
6521 CALL dbcsr_copy(m_delta(ispin), m_quench_t(ispin))
6524 CALL dbcsr_copy(m_residue(ispin), m_grad(ispin))
6527 do_exact_inversion = .false.
6528 IF (do_exact_inversion)
THEN
6532 CALL dbcsr_copy(m_step(ispin), m_grad(ispin))
6536 CALL hessian_diag_apply( &
6537 matrix_grad=m_step(ispin), &
6538 matrix_step=m_zet(ispin), &
6539 matrix_s_ao=m_s_vv(ispin), &
6540 matrix_f_ao=m_f_vv(ispin), &
6543 matrix_s_mo=m_siginv(ispin), &
6544 matrix_f_mo=m_siginvtftsiginv(ispin), &
6545 matrix_s_vo=m_stsiginv(ispin), &
6546 matrix_f_vo=m_f_vo(ispin), &
6547 quench_t=m_quench_t(ispin), &
6548 spin_factor=spin_factor, &
6549 eps_zero=eps_filter*10.0_dp, &
6550 penalty_occ_vol=penalty_occ_vol, &
6551 penalty_occ_vol_prefactor=penalty_occ_vol_prefactor(ispin), &
6552 penalty_occ_vol_pf2=penalty_occ_vol_pf2(ispin), &
6554 para_env=para_env, &
6555 blacs_env=blacs_env)
6559 IF (use_preconditioner)
THEN
6567 0.0_dp, m_zet(ispin), &
6568 filter_eps=eps_filter)
6573 matrix_in=m_residue(ispin), &
6574 matrix_out=m_zet(ispin), &
6575 operator1=domain_prec(:, ispin), &
6576 dpattern=m_quench_t(ispin), &
6577 map=domain_map(ispin), &
6578 node_of_domain=cpu_of_domain(:), &
6580 filter_eps=eps_filter)
6586 CALL dbcsr_copy(m_zet(ispin), m_residue(ispin))
6597 outer_prepare_to_exit = .false.
6599 residue_norm = 0.0_dp
6604 prepare_to_exit = .false.
6612 CALL apply_hessian( &
6617 m_siginv=m_siginv, &
6618 m_quench_t=m_quench_t, &
6619 m_ftsiginv=m_ftsiginv, &
6620 m_siginvtftsiginv=m_siginvtftsiginv, &
6622 m_stsiginv=m_stsiginv, &
6629 m_sig_sqrti_ii=m_sig_sqrti_ii, &
6630 penalty_occ_vol=penalty_occ_vol, &
6631 normalize_orbitals=normalize_orbitals, &
6632 penalty_occ_vol_prefactor=penalty_occ_vol_prefactor, &
6633 eps_filter=eps_filter, &
6634 path_num=hessian_path_reuse)
6639 DO ispin = 1, nspins
6641 CALL dbcsr_dot(m_residue(ispin), m_zet(ispin), numer_ispin)
6642 CALL dbcsr_dot(m_step(ispin), m_hstep(ispin), denom_ispin)
6644 numer = numer + numer_ispin
6645 denom = denom + denom_ispin
6651 DO ispin = 1, nspins
6654 CALL dbcsr_add(m_delta(ispin), m_step(ispin), 1.0_dp, alpha)
6655 CALL dbcsr_copy(m_residue_prev(ispin), m_residue(ispin))
6656 CALL dbcsr_add(m_residue(ispin), m_hstep(ispin), &
6657 1.0_dp, -1.0_dp*alpha)
6658 residue_max_norm(ispin) =
dbcsr_maxabs(m_residue(ispin))
6663 residue_norm = maxval(residue_max_norm)
6664 converged = (residue_norm < eps_error_target)
6665 IF (converged .OR. (iteration >= max_iter))
THEN
6666 prepare_to_exit = .true.
6669 IF (.NOT. prepare_to_exit)
THEN
6671 DO ispin = 1, nspins
6674 CALL dbcsr_copy(m_zet_prev(ispin), m_zet(ispin))
6677 IF (use_preconditioner)
THEN
6685 0.0_dp, m_zet(ispin), &
6686 filter_eps=eps_filter)
6691 matrix_in=m_residue(ispin), &
6692 matrix_out=m_zet(ispin), &
6693 operator1=domain_prec(:, ispin), &
6694 dpattern=m_quench_t(ispin), &
6695 map=domain_map(ispin), &
6696 node_of_domain=cpu_of_domain(:), &
6698 filter_eps=eps_filter)
6704 CALL dbcsr_copy(m_zet(ispin), m_residue(ispin))
6711 CALL compute_cg_beta( &
6713 reset_conjugator=reset_conjugator, &
6716 prev_grad=m_residue_prev, &
6718 prev_step=m_zet_prev)
6720 DO ispin = 1, nspins
6723 CALL dbcsr_add(m_step(ispin), m_zet(ispin), beta, 1.0_dp)
6730 IF (unit_nr > 0)
THEN
6731 iter_type = trim(
"NR STEP")
6732 WRITE (unit_nr,
'(T6,A9,I6,F14.5,F14.5,F15.10,F9.2)') &
6733 iter_type, iteration, &
6734 alpha, beta, residue_norm, &
6739 iteration = iteration + 1
6740 IF (prepare_to_exit)
EXIT
6744 IF (converged .OR. (outer_iteration >= outer_max_iter))
THEN
6745 outer_prepare_to_exit = .true.
6748 outer_iteration = outer_iteration + 1
6749 IF (outer_prepare_to_exit)
EXIT
6753 DO ispin = 1, nspins
6757 template=m_siginv(ispin), &
6758 matrix_type=dbcsr_type_no_symmetry)
6760 template=m_siginv(ispin), &
6761 matrix_type=dbcsr_type_no_symmetry)
6765 0.0_dp, m_tmp_oo_1, &
6766 filter_eps=eps_filter)
6770 0.0_dp, m_tmp_oo_2, &
6771 filter_eps=eps_filter)
6772 CALL dbcsr_copy(m_zet(ispin), m_quench_t(ispin))
6776 0.0_dp, m_zet(ispin), &
6777 retain_sparsity=.true.)
6779 WRITE (unit_nr,
"(A50,2F20.10)")
"Occupied-space projection of the step", alpha
6780 CALL dbcsr_add(m_zet(ispin), m_delta(ispin), -1.0_dp, 1.0_dp)
6782 WRITE (unit_nr,
"(A50,2F20.10)")
"Virtual-space projection of the step", alpha
6784 WRITE (unit_nr,
"(A50,2F20.10)")
"Full step", alpha
6791 DO ispin = 1, nspins
6805 DEALLOCATE (domain_prec)
6806 DEALLOCATE (m_residue)
6807 DEALLOCATE (m_residue_prev)
6810 DEALLOCATE (m_zet_prev)
6812 DEALLOCATE (m_hstep)
6816 DEALLOCATE (m_stsiginv)
6817 DEALLOCATE (residue_max_norm)
6819 IF (.NOT. converged)
THEN
6820 cpabort(
"Optimization not converged!")
6825 CALL timestop(handle)
6827 END SUBROUTINE newton_grad_to_step
6855 SUBROUTINE apply_hessian(m_x_in, m_x_out, m_ks, m_s, m_siginv, &
6856 m_quench_t, m_FTsiginv, m_siginvTFTsiginv, m_ST, m_STsiginv, m_s_vv, &
6857 m_ks_vv, m_g_full, m_t, m_sig_sqrti_ii, penalty_occ_vol, &
6858 normalize_orbitals, penalty_occ_vol_prefactor, eps_filter, path_num)
6860 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(INOUT) :: m_x_in, m_x_out, m_ks, m_s
6861 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(IN) :: m_siginv, m_quench_t, m_ftsiginv, &
6862 m_siginvtftsiginv, m_st, m_stsiginv
6863 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(INOUT) :: m_s_vv, m_ks_vv, m_g_full
6864 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(IN) :: m_t, m_sig_sqrti_ii
6865 LOGICAL,
INTENT(IN) :: penalty_occ_vol, normalize_orbitals
6866 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: penalty_occ_vol_prefactor
6867 REAL(kind=
dp),
INTENT(IN) :: eps_filter
6868 INTEGER,
INTENT(IN) :: path_num
6870 CHARACTER(len=*),
PARAMETER :: routinen =
'apply_hessian'
6872 INTEGER :: dim0, handle, ispin, nspins
6873 REAL(kind=
dp) :: penalty_prefactor_local, spin_factor
6874 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: tg_diagonal
6875 TYPE(
dbcsr_type) :: m_tmp_no_1, m_tmp_no_2, m_tmp_oo_1, &
6878 CALL timeset(routinen, handle)
6881 IF (penalty_occ_vol) penalty_prefactor_local = 1._dp
6882 cpassert(
SIZE(m_stsiginv) >= 0)
6883 cpassert(
SIZE(m_siginvtftsiginv) >= 0)
6884 cpassert(
SIZE(m_s) >= 0)
6885 cpassert(
SIZE(m_g_full) >= 0)
6886 cpassert(
SIZE(m_ftsiginv) >= 0)
6887 mark_used(m_siginvtftsiginv)
6888 mark_used(m_stsiginv)
6889 mark_used(m_ftsiginv)
6895 IF (nspins == 1)
THEN
6896 spin_factor = 2.0_dp
6898 spin_factor = 1.0_dp
6901 DO ispin = 1, nspins
6903 penalty_prefactor_local = penalty_occ_vol_prefactor(ispin)/(2.0_dp*spin_factor)
6906 template=m_siginv(ispin), &
6907 matrix_type=dbcsr_type_no_symmetry)
6909 template=m_quench_t(ispin), &
6910 matrix_type=dbcsr_type_no_symmetry)
6912 template=m_quench_t(ispin), &
6913 matrix_type=dbcsr_type_no_symmetry)
6915 template=m_quench_t(ispin), &
6916 matrix_type=dbcsr_type_no_symmetry)
6919 IF (normalize_orbitals)
THEN
6924 CALL dbcsr_copy(m_tmp_oo_1, m_sig_sqrti_ii(ispin))
6928 0.0_dp, m_tmp_oo_1, &
6929 retain_sparsity=.true.)
6931 ALLOCATE (tg_diagonal(dim0))
6935 DEALLOCATE (tg_diagonal)
6941 1.0_dp, m_tmp_no_1, &
6942 filter_eps=eps_filter)
6945 m_sig_sqrti_ii(ispin), &
6946 0.0_dp, m_tmp_x_in, &
6947 filter_eps=eps_filter)
6955 IF (path_num == hessian_path_reuse)
THEN
6960 CALL dbcsr_copy(m_x_out(ispin), m_quench_t(ispin))
6964 0.0_dp, m_x_out(ispin), &
6965 retain_sparsity=.true.)
6967 CALL dbcsr_copy(m_tmp_no_2, m_quench_t(ispin))
6971 0.0_dp, m_tmp_no_2, &
6972 retain_sparsity=.true.)
6973 CALL dbcsr_add(m_x_out(ispin), m_tmp_no_2, &
6974 1.0_dp, -4.0_dp*penalty_prefactor_local + 1.0_dp)
6976 ELSE IF (path_num == hessian_path_assemble)
THEN
6981 cpabort(
"path is NYI")
6984 cpabort(
"illegal path")
6988 IF (normalize_orbitals)
THEN
6993 CALL dbcsr_copy(m_tmp_oo_1, m_sig_sqrti_ii(ispin))
6997 0.0_dp, m_tmp_oo_1, &
6998 retain_sparsity=.true.)
7000 ALLOCATE (tg_diagonal(dim0))
7004 DEALLOCATE (tg_diagonal)
7009 1.0_dp, m_x_out(ispin), &
7010 retain_sparsity=.true.)
7014 m_sig_sqrti_ii(ispin), &
7015 0.0_dp, m_x_out(ispin), &
7016 retain_sparsity=.true.)
7034 CALL timestop(handle)
7036 END SUBROUTINE apply_hessian
7061 SUBROUTINE hessian_diag_apply(matrix_grad, matrix_step, matrix_S_ao, &
7062 matrix_F_ao, matrix_S_mo, matrix_F_mo, matrix_S_vo, matrix_F_vo, quench_t, &
7063 penalty_occ_vol, penalty_occ_vol_prefactor, penalty_occ_vol_pf2, &
7064 spin_factor, eps_zero, m_s, para_env, blacs_env)
7066 TYPE(
dbcsr_type),
INTENT(INOUT) :: matrix_grad, matrix_step, matrix_s_ao, &
7067 matrix_f_ao, matrix_s_mo
7069 TYPE(
dbcsr_type),
INTENT(INOUT) :: matrix_s_vo, matrix_f_vo, quench_t
7070 LOGICAL,
INTENT(IN) :: penalty_occ_vol
7071 REAL(kind=
dp),
INTENT(IN) :: penalty_occ_vol_prefactor, &
7072 penalty_occ_vol_pf2, spin_factor, &
7078 CHARACTER(len=*),
PARAMETER :: routinen =
'hessian_diag_apply'
7080 INTEGER :: ao_hori_offset, ao_vert_offset, block_col, block_row, col, h_size, handle, ii, &
7081 info, jj, lev1_hori_offset, lev1_vert_offset, lev2_hori_offset, lev2_vert_offset, lwork, &
7082 nblkcols_tot, nblkrows_tot, ncores, orb_i, orb_j, row, unit_nr, zero_neg_eiv
7083 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: ao_block_sizes, ao_domain_sizes, &
7085 INTEGER,
DIMENSION(:),
POINTER :: ao_blk_sizes, mo_blk_sizes
7086 LOGICAL :: found, found_col, found_row
7087 REAL(kind=
dp) :: penalty_prefactor_local, test_error
7088 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: eigenvalues, grad_vec, step_vec, tmp, &
7090 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: f_ao_block, f_mo_block, h, hinv, &
7091 new_block, s_ao_block, s_mo_block, &
7093 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: block_p
7096 TYPE(
dbcsr_type) :: matrix_f_ao_sym, matrix_f_mo_sym, &
7097 matrix_s_ao_sym, matrix_s_mo_sym
7099 CALL timeset(routinen, handle)
7103 IF (logger%para_env%is_source())
THEN
7110 cpassert(
ASSOCIATED(blacs_env))
7111 cpassert(
ASSOCIATED(para_env))
7112 mark_used(blacs_env)
7122 IF (ncores > 1)
THEN
7123 cpabort(
"serial code only")
7126 CALL dbcsr_get_info(quench_t, row_blk_size=ao_blk_sizes, col_blk_size=mo_blk_sizes, &
7127 nblkrows_total=nblkrows_tot, nblkcols_total=nblkcols_tot)
7128 cpassert(nblkrows_tot == nblkcols_tot)
7129 ALLOCATE (mo_block_sizes(nblkcols_tot), ao_block_sizes(nblkcols_tot))
7130 ALLOCATE (ao_domain_sizes(nblkcols_tot))
7131 mo_block_sizes(:) = mo_blk_sizes(:)
7132 ao_block_sizes(:) = ao_blk_sizes(:)
7133 ao_domain_sizes(:) = 0
7136 template=matrix_s_ao, &
7137 matrix_type=dbcsr_type_no_symmetry)
7139 CALL dbcsr_scale(matrix_s_ao_sym, 2.0_dp*spin_factor)
7142 template=matrix_f_ao, &
7143 matrix_type=dbcsr_type_no_symmetry)
7145 CALL dbcsr_scale(matrix_f_ao_sym, 2.0_dp*spin_factor)
7148 template=matrix_s_mo, &
7149 matrix_type=dbcsr_type_no_symmetry)
7153 template=matrix_f_mo, &
7154 matrix_type=dbcsr_type_no_symmetry)
7157 IF (penalty_occ_vol)
THEN
7158 penalty_prefactor_local = penalty_occ_vol_prefactor/(2.0_dp*spin_factor)
7160 penalty_prefactor_local = 0.0_dp
7163 WRITE (unit_nr, *)
"penalty_prefactor_local: ", penalty_prefactor_local
7164 WRITE (unit_nr, *)
"penalty_prefactor_2: ", penalty_occ_vol_pf2
7168 DO col = 1, nblkcols_tot
7171 DO row = 1, nblkrows_tot
7174 row, col, block_p, found)
7176 ao_domain_sizes(col) = ao_domain_sizes(col) + ao_blk_sizes(row)
7181 h_size = h_size + ao_domain_sizes(col)*mo_block_sizes(col)
7185 ALLOCATE (h(h_size, h_size))
7189 lev1_vert_offset = 0
7191 DO row = 1, nblkcols_tot
7193 lev1_hori_offset = 0
7194 DO col = 1, nblkcols_tot
7197 ALLOCATE (f_ao_block(ao_domain_sizes(row), ao_domain_sizes(col)))
7198 ALLOCATE (s_ao_block(ao_domain_sizes(row), ao_domain_sizes(col)))
7199 ALLOCATE (f_mo_block(mo_block_sizes(row), mo_block_sizes(col)))
7200 ALLOCATE (s_mo_block(mo_block_sizes(row), mo_block_sizes(col)))
7202 f_ao_block(:, :) = 0.0_dp
7203 s_ao_block(:, :) = 0.0_dp
7204 f_mo_block(:, :) = 0.0_dp
7205 s_mo_block(:, :) = 0.0_dp
7210 DO block_row = 1, nblkcols_tot
7213 block_row, row, block_p, found_row)
7217 DO block_col = 1, nblkcols_tot
7220 block_col, col, block_p, found_col)
7224 block_row, block_col, block_p, found)
7227 f_ao_block(ao_vert_offset + 1:ao_vert_offset + ao_block_sizes(block_row), &
7228 ao_hori_offset + 1:ao_hori_offset + ao_block_sizes(block_col)) &
7233 block_row, block_col, block_p, found)
7236 s_ao_block(ao_vert_offset + 1:ao_vert_offset + ao_block_sizes(block_row), &
7237 ao_hori_offset + 1:ao_hori_offset + ao_block_sizes(block_col)) &
7241 ao_hori_offset = ao_hori_offset + ao_block_sizes(block_col)
7247 ao_vert_offset = ao_vert_offset + ao_block_sizes(block_row)
7257 f_mo_block(1:mo_block_sizes(row), 1:mo_block_sizes(col)) = block_p(:, :)
7262 s_mo_block(1:mo_block_sizes(row), 1:mo_block_sizes(col)) = block_p(:, :)
7266 lev2_vert_offset = 0
7267 DO orb_j = 1, mo_block_sizes(row)
7269 lev2_hori_offset = 0
7270 DO orb_i = 1, mo_block_sizes(col)
7271 IF (orb_i == orb_j .AND. row == col)
THEN
7272 h(lev1_vert_offset + lev2_vert_offset + 1:lev1_vert_offset + lev2_vert_offset + ao_domain_sizes(row), &
7273 lev1_hori_offset + lev2_hori_offset + 1:lev1_hori_offset + lev2_hori_offset + ao_domain_sizes(col)) &
7274 = f_ao_block(:, :) + s_ao_block(:, :)
7277 lev2_hori_offset = lev2_hori_offset + ao_domain_sizes(col)
7281 lev2_vert_offset = lev2_vert_offset + ao_domain_sizes(row)
7285 lev1_hori_offset = lev1_hori_offset + ao_domain_sizes(col)*mo_block_sizes(col)
7287 DEALLOCATE (f_ao_block)
7288 DEALLOCATE (s_ao_block)
7289 DEALLOCATE (f_mo_block)
7290 DEALLOCATE (s_mo_block)
7294 lev1_vert_offset = lev1_vert_offset + ao_domain_sizes(row)*mo_block_sizes(row)
7304 ALLOCATE (grad_vec(h_size))
7305 grad_vec(:) = 0.0_dp
7306 lev1_vert_offset = 0
7308 DO col = 1, nblkcols_tot
7311 lev2_vert_offset = 0
7312 DO row = 1, nblkrows_tot
7315 row, col, block_p, found_row)
7319 row, col, block_p, found)
7322 DO orb_i = 1, mo_block_sizes(col)
7323 grad_vec(lev1_vert_offset + ao_domain_sizes(col)*(orb_i - 1) + lev2_vert_offset + 1: &
7324 lev1_vert_offset + ao_domain_sizes(col)*(orb_i - 1) + lev2_vert_offset + ao_block_sizes(row)) &
7330 lev2_vert_offset = lev2_vert_offset + ao_block_sizes(row)
7336 lev1_vert_offset = lev1_vert_offset + ao_domain_sizes(col)*mo_block_sizes(col)
7342 ALLOCATE (hinv(h_size, h_size))
7343 hinv(:, :) = h(:, :)
7346 ALLOCATE (eigenvalues(h_size))
7349 ALLOCATE (work(max(1, lwork)))
7350 CALL dsyev(
'V',
'L', h_size, hinv, h_size, eigenvalues, work, lwork, info)
7351 lwork = int(work(1))
7354 ALLOCATE (work(max(1, lwork)))
7355 CALL dsyev(
'V',
'L', h_size, hinv, h_size, eigenvalues, work, lwork, info)
7357 WRITE (unit_nr, *)
'DSYEV ERROR MESSAGE: ', info
7358 cpabort(
"DSYEV failed")
7363 ALLOCATE (step_vec(h_size))
7365 step_vec(:) = matmul(transpose(hinv), grad_vec)
7369 ALLOCATE (test(h_size, h_size))
7372 WRITE (unit_nr,
"(I10,F20.10,F20.10)") jj, eigenvalues(jj), step_vec(jj)
7373 IF (eigenvalues(jj) > eps_zero)
THEN
7374 test(jj, :) = hinv(:, jj)/eigenvalues(jj)
7376 test(jj, :) = hinv(:, jj)*0.0_dp
7377 zero_neg_eiv = zero_neg_eiv + 1
7380 WRITE (unit_nr, *)
'ZERO OR NEGATIVE EIGENVALUES: ', zero_neg_eiv
7381 DEALLOCATE (step_vec)
7383 ALLOCATE (test2(h_size, h_size))
7384 test2(:, :) = matmul(hinv, test)
7385 hinv(:, :) = test2(:, :)
7386 DEALLOCATE (test, test2)
7388 DEALLOCATE (eigenvalues)
7391 ALLOCATE (test(h_size, h_size))
7392 test(:, :) = matmul(hinv, h)
7394 test(ii, ii) = test(ii, ii) - 1.0_dp
7399 test_error = test_error + test(jj, ii)*test(jj, ii)
7402 WRITE (unit_nr, *)
"Hessian inversion error: ", sqrt(test_error)
7406 ALLOCATE (step_vec(h_size))
7407 ALLOCATE (tmp(h_size))
7408 tmp(:) = matmul(hinv, grad_vec)
7409 step_vec(:) = -1.0_dp*tmp(:)
7411 ALLOCATE (tmpr(h_size))
7412 tmpr(:) = matmul(h, step_vec)
7413 tmp(:) = tmpr(:) + grad_vec(:)
7415 WRITE (unit_nr, *)
"NEWTOV step error: ", maxval(abs(tmp))
7421 DEALLOCATE (grad_vec)
7429 template=matrix_grad, &
7430 matrix_type=dbcsr_type_no_symmetry)
7433 lev1_vert_offset = 0
7435 DO col = 1, nblkcols_tot
7438 lev2_vert_offset = 0
7439 DO row = 1, nblkrows_tot
7442 row, col, block_p, found_row)
7445 ALLOCATE (new_block(ao_block_sizes(row), mo_block_sizes(col)))
7446 DO orb_i = 1, mo_block_sizes(col)
7447 new_block(:, orb_i) = &
7448 step_vec(lev1_vert_offset + ao_domain_sizes(col)*(orb_i - 1) + lev2_vert_offset + 1: &
7449 lev1_vert_offset + ao_domain_sizes(col)*(orb_i - 1) + lev2_vert_offset + ao_block_sizes(row))
7452 DEALLOCATE (new_block)
7453 lev2_vert_offset = lev2_vert_offset + ao_block_sizes(row)
7458 lev1_vert_offset = lev1_vert_offset + ao_domain_sizes(col)*mo_block_sizes(col)
7462 DEALLOCATE (step_vec)
7466 DEALLOCATE (mo_block_sizes, ao_block_sizes)
7467 DEALLOCATE (ao_domain_sizes)
7470 template=quench_t, &
7471 matrix_type=dbcsr_type_no_symmetry)
7476 0.0_dp, matrix_s_ao_sym, &
7477 retain_sparsity=.true.)
7479 template=quench_t, &
7480 matrix_type=dbcsr_type_no_symmetry)
7485 0.0_dp, matrix_f_ao_sym, &
7486 retain_sparsity=.true.)
7487 CALL dbcsr_add(matrix_s_ao_sym, matrix_f_ao_sym, &
7489 CALL dbcsr_scale(matrix_s_ao_sym, 2.0_dp*spin_factor)
7490 CALL dbcsr_add(matrix_s_ao_sym, matrix_grad, &
7493 WRITE (unit_nr, *)
"NEWTOL step error: ", test_error
7497 CALL timestop(handle)
7499 END SUBROUTINE hessian_diag_apply
7519 matrix_t_in, matrix_t_out, perturbation_only, &
7525 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:) :: quench_t, matrix_t_in, matrix_t_out
7526 LOGICAL,
INTENT(IN) :: perturbation_only
7527 INTEGER,
INTENT(IN),
OPTIONAL :: special_case
7529 CHARACTER(len=*),
PARAMETER :: routinen =
'almo_scf_xalmo_trustr'
7531 INTEGER :: handle, ispin, iteration, iteration_type_to_report, my_special_case, ndomains, &
7532 nspins, outer_iteration, prec_type, unit_nr
7533 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: nocc
7534 LOGICAL :: assume_t0_q0x, border_reached, inner_loop_success, normalize_orbitals, &
7535 optimize_theta, penalty_occ_vol, reset_conjugator, same_position, scf_converged
7536 REAL(kind=
dp) :: beta, energy_start, energy_trial, eta, expected_reduction, &
7537 fake_step_size_to_report, grad_norm_ratio, grad_norm_ref, loss_change_to_report, &
7538 loss_start, loss_trial, model_grad_norm, penalty_amplitude, penalty_start, penalty_trial, &
7539 radius_current, radius_max, real_temp, rho, spin_factor, step_norm, step_size, t1, &
7540 t1outer, t2, t2outer, y_scalar
7541 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: grad_norm_spin, &
7542 penalty_occ_vol_g_prefactor, &
7543 penalty_occ_vol_h_prefactor
7546 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:) :: ftsiginv, grad, m_model_bd, m_model_d, &
7547 m_model_hessian, m_model_hessian_inv, m_model_r, m_model_r_prev, m_model_rt, &
7548 m_model_rt_prev, m_sig_sqrti_ii, m_theta, m_theta_trial, prev_step, siginvtftsiginv, st, &
7551 DIMENSION(:, :) :: domain_model_hessian_inv, domain_r_down
7554 CALL timeset(routinen, handle)
7559 IF (
PRESENT(special_case)) my_special_case = special_case
7563 IF (logger%para_env%is_source())
THEN
7570 assume_t0_q0x = .false.
7572 optimize_theta = .false.
7574 nspins = almo_scf_env%nspins
7575 IF (nspins == 1)
THEN
7576 spin_factor = 2.0_dp
7578 spin_factor = 1.0_dp
7581 IF (unit_nr > 0)
THEN
7583 SELECT CASE (my_special_case)
7585 WRITE (unit_nr,
'(T2,A,A,A)') repeat(
"-", 20), &
7586 " Optimization of block-diagonal ALMOs ", repeat(
"-", 21)
7588 WRITE (unit_nr,
'(T2,A,A,A)') repeat(
"-", 20), &
7589 " Optimization of fully delocalized MOs ", repeat(
"-", 20)
7591 WRITE (unit_nr,
'(T2,A,A,A)') repeat(
"-", 27), &
7592 " Optimization of XALMOs ", repeat(
"-", 28)
7595 CALL trust_r_report(unit_nr, &
7600 delta_loss=0.0_dp, &
7602 predicted_reduction=0.0_dp, &
7606 WRITE (unit_nr,
'(T2,A)') repeat(
"-", 79)
7610 penalty_occ_vol = .false.
7611 normalize_orbitals = penalty_occ_vol
7612 penalty_amplitude = 0.0_dp
7613 ALLOCATE (penalty_occ_vol_g_prefactor(nspins))
7614 ALLOCATE (penalty_occ_vol_h_prefactor(nspins))
7615 penalty_occ_vol_g_prefactor(:) = 0.0_dp
7616 penalty_occ_vol_h_prefactor(:) = 0.0_dp
7619 prec_type = optimizer%preconditioner
7621 ALLOCATE (grad_norm_spin(nspins))
7622 ALLOCATE (nocc(nspins))
7626 ALLOCATE (m_theta(nspins))
7627 DO ispin = 1, nspins
7629 template=matrix_t_out(ispin), &
7630 matrix_type=dbcsr_type_no_symmetry)
7635 m_t_in=matrix_t_in, &
7636 m_t0=almo_scf_env%matrix_t_blk, &
7637 m_quench_t=quench_t, &
7638 m_overlap=almo_scf_env%matrix_s(1), &
7639 m_sigma_tmpl=almo_scf_env%matrix_sigma_inv, &
7641 xalmo_history=almo_scf_env%xalmo_history, &
7642 assume_t0_q0x=assume_t0_q0x, &
7643 optimize_theta=optimize_theta, &
7644 envelope_amplitude=almo_scf_env%envelope_amplitude, &
7645 eps_filter=almo_scf_env%eps_filter, &
7646 order_lanczos=almo_scf_env%order_lanczos, &
7647 eps_lanczos=almo_scf_env%eps_lanczos, &
7648 max_iter_lanczos=almo_scf_env%max_iter_lanczos, &
7649 nocc_of_domain=almo_scf_env%nocc_of_domain)
7651 ndomains = almo_scf_env%ndomains
7652 ALLOCATE (domain_r_down(ndomains, nspins))
7654 ALLOCATE (domain_model_hessian_inv(ndomains, nspins))
7657 ALLOCATE (m_model_hessian(nspins))
7658 ALLOCATE (m_model_hessian_inv(nspins))
7659 ALLOCATE (siginvtftsiginv(nspins))
7660 ALLOCATE (stsiginv_0(nspins))
7661 ALLOCATE (ftsiginv(nspins))
7662 ALLOCATE (st(nspins))
7663 ALLOCATE (grad(nspins))
7664 ALLOCATE (prev_step(nspins))
7665 ALLOCATE (step(nspins))
7666 ALLOCATE (m_sig_sqrti_ii(nspins))
7667 ALLOCATE (m_model_r(nspins))
7668 ALLOCATE (m_model_rt(nspins))
7669 ALLOCATE (m_model_d(nspins))
7670 ALLOCATE (m_model_bd(nspins))
7671 ALLOCATE (m_model_r_prev(nspins))
7672 ALLOCATE (m_model_rt_prev(nspins))
7673 ALLOCATE (m_theta_trial(nspins))
7675 DO ispin = 1, nspins
7679 template=almo_scf_env%matrix_ks(ispin), &
7680 matrix_type=dbcsr_type_no_symmetry)
7682 template=almo_scf_env%matrix_ks(ispin), &
7683 matrix_type=dbcsr_type_no_symmetry)
7685 template=almo_scf_env%matrix_sigma(ispin), &
7686 matrix_type=dbcsr_type_no_symmetry)
7688 template=matrix_t_out(ispin), &
7689 matrix_type=dbcsr_type_no_symmetry)
7691 template=matrix_t_out(ispin), &
7692 matrix_type=dbcsr_type_no_symmetry)
7694 template=matrix_t_out(ispin), &
7695 matrix_type=dbcsr_type_no_symmetry)
7697 template=matrix_t_out(ispin), &
7698 matrix_type=dbcsr_type_no_symmetry)
7700 template=matrix_t_out(ispin), &
7701 matrix_type=dbcsr_type_no_symmetry)
7703 template=matrix_t_out(ispin), &
7704 matrix_type=dbcsr_type_no_symmetry)
7706 template=almo_scf_env%matrix_sigma_inv(ispin), &
7707 matrix_type=dbcsr_type_no_symmetry)
7709 template=matrix_t_out(ispin), &
7710 matrix_type=dbcsr_type_no_symmetry)
7712 template=matrix_t_out(ispin), &
7713 matrix_type=dbcsr_type_no_symmetry)
7715 template=matrix_t_out(ispin), &
7716 matrix_type=dbcsr_type_no_symmetry)
7718 template=matrix_t_out(ispin), &
7719 matrix_type=dbcsr_type_no_symmetry)
7721 template=matrix_t_out(ispin), &
7722 matrix_type=dbcsr_type_no_symmetry)
7724 template=matrix_t_out(ispin), &
7725 matrix_type=dbcsr_type_no_symmetry)
7727 template=matrix_t_out(ispin), &
7728 matrix_type=dbcsr_type_no_symmetry)
7731 CALL dbcsr_set(prev_step(ispin), 0.0_dp)
7734 nfullrows_total=nocc(ispin))
7742 matrix_s=almo_scf_env%matrix_s(1), &
7743 subm_s_inv=almo_scf_env%domain_s_inv(:, ispin), &
7744 dpattern=quench_t(ispin), &
7745 map=almo_scf_env%domain_map(ispin), &
7746 node_of_domain=almo_scf_env%cpu_of_domain)
7756 template=almo_scf_env%matrix_s(1), &
7757 matrix_type=dbcsr_type_no_symmetry)
7759 almo_scf_env%matrix_s_blk(1), &
7760 threshold=almo_scf_env%eps_filter, &
7761 filter_eps=almo_scf_env%eps_filter)
7767 template=almo_scf_env%matrix_s(1), &
7768 matrix_type=dbcsr_type_no_symmetry)
7771 para_env=almo_scf_env%para_env, &
7772 blacs_env=almo_scf_env%blacs_env)
7774 para_env=almo_scf_env%para_env, &
7775 blacs_env=almo_scf_env%blacs_env, &
7776 uplo_to_full=.true.)
7781 radius_max = optimizer%max_trust_radius
7782 radius_current = min(optimizer%initial_trust_radius, radius_max)
7784 eta = min(max(optimizer%rho_do_not_update, 0.0_dp), 0.25_dp)
7785 energy_start = 0.0_dp
7786 energy_trial = 0.0_dp
7787 penalty_start = 0.0_dp
7788 penalty_trial = 0.0_dp
7792 same_position = .false.
7795 CALL main_var_to_xalmos_and_loss_func( &
7796 almo_scf_env=almo_scf_env, &
7798 m_main_var_in=m_theta, &
7799 m_t_out=matrix_t_out, &
7800 m_sig_sqrti_ii_out=m_sig_sqrti_ii, &
7801 energy_out=energy_start, &
7802 penalty_out=penalty_start, &
7803 m_ftsiginv_out=ftsiginv, &
7804 m_siginvtftsiginv_out=siginvtftsiginv, &
7806 m_stsiginv0_in=stsiginv_0, &
7807 m_quench_t_in=quench_t, &
7808 domain_r_down_in=domain_r_down, &
7809 assume_t0_q0x=assume_t0_q0x, &
7810 just_started=.true., &
7811 optimize_theta=optimize_theta, &
7812 normalize_orbitals=normalize_orbitals, &
7813 perturbation_only=perturbation_only, &
7814 do_penalty=penalty_occ_vol, &
7815 special_case=my_special_case)
7816 loss_start = energy_start + penalty_start
7818 almo_scf_env%almo_scf_energy = energy_start
7820 DO ispin = 1, nspins
7821 IF (penalty_occ_vol)
THEN
7822 penalty_occ_vol_g_prefactor(ispin) = &
7823 -2.0_dp*penalty_amplitude*spin_factor*nocc(ispin)
7824 penalty_occ_vol_h_prefactor(ispin) = 0.0_dp
7829 scf_converged = .false.
7830 adjust_r_loop:
DO outer_iteration = 1, optimizer%max_iter_outer_loop
7833 border_reached = .false.
7835 DO ispin = 1, nspins
7837 CALL dbcsr_filter(step(ispin), almo_scf_env%eps_filter)
7840 IF (.NOT. same_position)
THEN
7842 DO ispin = 1, nspins
7844 IF (unit_nr > 0 .AND. debug_mode)
WRITE (unit_nr, *)
"...Compute model gradient"
7845 CALL compute_gradient( &
7846 m_grad_out=grad(ispin), &
7847 m_ks=almo_scf_env%matrix_ks(ispin), &
7848 m_s=almo_scf_env%matrix_s(1), &
7849 m_t=matrix_t_out(ispin), &
7850 m_t0=almo_scf_env%matrix_t_blk(ispin), &
7851 m_siginv=almo_scf_env%matrix_sigma_inv(ispin), &
7852 m_quench_t=quench_t(ispin), &
7853 m_ftsiginv=ftsiginv(ispin), &
7854 m_siginvtftsiginv=siginvtftsiginv(ispin), &
7856 m_stsiginv0=stsiginv_0(ispin), &
7857 m_theta=m_theta(ispin), &
7858 m_sig_sqrti_ii=m_sig_sqrti_ii(ispin), &
7859 domain_s_inv=almo_scf_env%domain_s_inv(:, ispin), &
7860 domain_r_down=domain_r_down(:, ispin), &
7861 cpu_of_domain=almo_scf_env%cpu_of_domain, &
7862 domain_map=almo_scf_env%domain_map(ispin), &
7863 assume_t0_q0x=assume_t0_q0x, &
7864 optimize_theta=optimize_theta, &
7865 normalize_orbitals=normalize_orbitals, &
7866 penalty_occ_vol=penalty_occ_vol, &
7867 penalty_occ_vol_prefactor=penalty_occ_vol_g_prefactor(ispin), &
7868 envelope_amplitude=almo_scf_env%envelope_amplitude, &
7869 eps_filter=almo_scf_env%eps_filter, &
7870 spin_factor=spin_factor, &
7871 special_case=my_special_case)
7878 DO ispin = 1, nspins
7881 grad_norm_ref = maxval(grad_norm_spin)
7884 CALL trust_r_report(unit_nr, &
7886 iteration=outer_iteration, &
7888 delta_loss=0.0_dp, &
7889 grad_norm=grad_norm_ref, &
7890 predicted_reduction=0.0_dp, &
7892 radius=radius_current, &
7893 new=.NOT. same_position, &
7894 time=t2outer - t1outer)
7897 IF (grad_norm_ref <= optimizer%eps_error)
THEN
7898 scf_converged = .true.
7899 border_reached = .false.
7900 expected_reduction = 0.0_dp
7901 IF (.NOT. (optimizer%early_stopping_on .AND. outer_iteration == 1))
THEN
7905 scf_converged = .false.
7908 DO ispin = 1, nspins
7910 CALL dbcsr_copy(m_model_r(ispin), grad(ispin))
7916 IF (unit_nr > 0 .AND. debug_mode)
WRITE (unit_nr, *)
"...Multiply Sinv.r"
7920 0.0_dp, m_model_rt(ispin), &
7921 filter_eps=almo_scf_env%eps_filter)
7925 IF (unit_nr > 0 .AND. debug_mode)
WRITE (unit_nr, *)
"...Multiply Sinv_xx.r"
7927 matrix_in=m_model_r(ispin), &
7928 matrix_out=m_model_rt(ispin), &
7929 operator1=almo_scf_env%domain_s_inv(:, ispin), &
7930 dpattern=quench_t(ispin), &
7931 map=almo_scf_env%domain_map(ispin), &
7932 node_of_domain=almo_scf_env%cpu_of_domain, &
7934 filter_eps=almo_scf_env%eps_filter)
7937 cpabort(
"Unknown XALMO special case")
7940 CALL dbcsr_copy(m_model_d(ispin), m_model_rt(ispin))
7945 IF (.NOT. same_position)
THEN
7947 SELECT CASE (prec_type)
7950 IF (unit_nr > 0 .AND. debug_mode)
WRITE (unit_nr, *)
"...Compute model Hessian"
7951 DO ispin = 1, nspins
7952 CALL compute_preconditioner( &
7953 domain_prec_out=almo_scf_env%domain_preconditioner(:, ispin), &
7954 m_prec_out=m_model_hessian(ispin), &
7955 m_ks=almo_scf_env%matrix_ks(ispin), &
7956 m_s=almo_scf_env%matrix_s(1), &
7957 m_siginv=almo_scf_env%matrix_sigma_inv(ispin), &
7958 m_quench_t=quench_t(ispin), &
7959 m_ftsiginv=ftsiginv(ispin), &
7960 m_siginvtftsiginv=siginvtftsiginv(ispin), &
7962 para_env=almo_scf_env%para_env, &
7963 blacs_env=almo_scf_env%blacs_env, &
7964 nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin), &
7965 domain_s_inv=almo_scf_env%domain_s_inv(:, ispin), &
7966 domain_r_down=domain_r_down(:, ispin), &
7967 cpu_of_domain=almo_scf_env%cpu_of_domain, &
7968 domain_map=almo_scf_env%domain_map(ispin), &
7969 assume_t0_q0x=.false., &
7970 penalty_occ_vol=penalty_occ_vol, &
7971 penalty_occ_vol_prefactor=penalty_occ_vol_g_prefactor(ispin), &
7972 eps_filter=almo_scf_env%eps_filter, &
7974 spin_factor=spin_factor, &
7975 skip_inversion=.true., &
7976 special_case=my_special_case)
7981 cpabort(
"Unknown preconditioner")
7988 CALL fixed_r_report(unit_nr, &
7992 border_reached=.false., &
7994 grad_norm_ratio=0.0_dp, &
7997 IF (unit_nr > 0 .AND. debug_mode)
WRITE (unit_nr, *)
"...Start inner loop"
8000 inner_loop_success = .false.
8002 fixed_r_loop:
DO iteration = 1, optimizer%max_iter
8006 DO ispin = 1, nspins
8013 m_model_hessian(ispin), &
8015 0.0_dp, m_model_bd(ispin), &
8016 filter_eps=almo_scf_env%eps_filter)
8021 matrix_in=m_model_d(ispin), &
8022 matrix_out=m_model_bd(ispin), &
8023 operator1=almo_scf_env%domain_preconditioner(:, ispin), &
8024 dpattern=quench_t(ispin), &
8025 map=almo_scf_env%domain_map(ispin), &
8026 node_of_domain=almo_scf_env%cpu_of_domain, &
8028 filter_eps=almo_scf_env%eps_filter)
8033 CALL dbcsr_dot(m_model_d(ispin), m_model_bd(ispin), real_temp)
8034 y_scalar = y_scalar + real_temp
8037 IF (unit_nr > 0 .AND. debug_mode)
WRITE (unit_nr, *)
"...Curvature: ", y_scalar
8040 IF (y_scalar < 0.0_dp)
THEN
8042 CALL step_size_to_border( &
8043 step_size_out=step_size, &
8044 metric_in=almo_scf_env%matrix_s, &
8046 direction_in=m_model_d, &
8047 trust_radius_in=radius_current, &
8048 quench_t_in=quench_t, &
8049 eps_filter_in=almo_scf_env%eps_filter &
8052 DO ispin = 1, nspins
8053 CALL dbcsr_add(step(ispin), m_model_d(ispin), 1.0_dp, step_size)
8056 border_reached = .true.
8057 inner_loop_success = .true.
8059 CALL predicted_reduction( &
8060 reduction_out=expected_reduction, &
8063 hess_in=m_model_hessian, &
8064 hess_submatrix_in=almo_scf_env%domain_preconditioner, &
8065 quench_t_in=quench_t, &
8066 special_case=my_special_case, &
8067 eps_filter=almo_scf_env%eps_filter, &
8068 domain_map=almo_scf_env%domain_map, &
8069 cpu_of_domain=almo_scf_env%cpu_of_domain &
8073 CALL fixed_r_report(unit_nr, &
8075 iteration=iteration, &
8076 step_size=step_size, &
8077 border_reached=border_reached, &
8078 curvature=y_scalar, &
8079 grad_norm_ratio=expected_reduction, &
8088 DO ispin = 1, nspins
8089 CALL dbcsr_dot(m_model_r(ispin), m_model_rt(ispin), real_temp)
8090 step_size = step_size + real_temp
8092 step_size = step_size/y_scalar
8093 IF (unit_nr > 0 .AND. debug_mode)
WRITE (unit_nr, *)
"...Proposed step size: ", step_size
8096 DO ispin = 1, nspins
8097 CALL dbcsr_copy(prev_step(ispin), step(ispin))
8098 CALL dbcsr_add(step(ispin), m_model_d(ispin), 1.0_dp, step_size)
8102 CALL contravariant_matrix_norm( &
8103 norm_out=step_norm, &
8105 metric_in=almo_scf_env%matrix_s, &
8106 quench_t_in=quench_t, &
8107 eps_filter_in=almo_scf_env%eps_filter &
8109 IF (unit_nr > 0 .AND. debug_mode)
WRITE (unit_nr, *)
"...Step norm: ", step_norm
8112 IF (step_norm > radius_current)
THEN
8114 IF (unit_nr > 0 .AND. debug_mode)
WRITE (unit_nr, *)
"...Norm is too large"
8115 CALL step_size_to_border( &
8116 step_size_out=step_size, &
8117 metric_in=almo_scf_env%matrix_s, &
8118 position_in=prev_step, &
8119 direction_in=m_model_d, &
8120 trust_radius_in=radius_current, &
8121 quench_t_in=quench_t, &
8122 eps_filter_in=almo_scf_env%eps_filter &
8124 IF (unit_nr > 0 .AND. debug_mode)
WRITE (unit_nr, *)
"...Step size to border: ", step_size
8126 DO ispin = 1, nspins
8127 CALL dbcsr_copy(step(ispin), prev_step(ispin))
8128 CALL dbcsr_add(step(ispin), m_model_d(ispin), 1.0_dp, step_size)
8131 IF (debug_mode)
THEN
8133 IF (unit_nr > 0)
WRITE (unit_nr, *)
"...Extra norm evaluation"
8134 CALL contravariant_matrix_norm( &
8135 norm_out=step_norm, &
8137 metric_in=almo_scf_env%matrix_s, &
8138 quench_t_in=quench_t, &
8139 eps_filter_in=almo_scf_env%eps_filter &
8141 IF (unit_nr > 0)
WRITE (unit_nr, *)
"...Step norm: ", step_norm
8142 IF (unit_nr > 0)
WRITE (unit_nr, *)
"...Current radius: ", radius_current
8145 border_reached = .true.
8146 inner_loop_success = .true.
8148 CALL predicted_reduction( &
8149 reduction_out=expected_reduction, &
8152 hess_in=m_model_hessian, &
8153 hess_submatrix_in=almo_scf_env%domain_preconditioner, &
8154 quench_t_in=quench_t, &
8155 special_case=my_special_case, &
8156 eps_filter=almo_scf_env%eps_filter, &
8157 domain_map=almo_scf_env%domain_map, &
8158 cpu_of_domain=almo_scf_env%cpu_of_domain &
8162 CALL fixed_r_report(unit_nr, &
8164 iteration=iteration, &
8165 step_size=step_size, &
8166 border_reached=border_reached, &
8167 curvature=y_scalar, &
8168 grad_norm_ratio=expected_reduction, &
8178 border_reached = .false.
8179 inner_loop_success = .true.
8181 CALL predicted_reduction( &
8182 reduction_out=expected_reduction, &
8185 hess_in=m_model_hessian, &
8186 hess_submatrix_in=almo_scf_env%domain_preconditioner, &
8187 quench_t_in=quench_t, &
8188 special_case=my_special_case, &
8189 eps_filter=almo_scf_env%eps_filter, &
8190 domain_map=almo_scf_env%domain_map, &
8191 cpu_of_domain=almo_scf_env%cpu_of_domain &
8195 CALL fixed_r_report(unit_nr, &
8197 iteration=iteration, &
8198 step_size=step_size, &
8199 border_reached=border_reached, &
8200 curvature=y_scalar, &
8201 grad_norm_ratio=expected_reduction, &
8209 SELECT CASE (prec_type)
8212 IF (unit_nr > 0 .AND. debug_mode)
WRITE (unit_nr, *)
"...Pseudo-invert model Hessian"
8215 DO ispin = 1, nspins
8217 matrix_in=m_model_hessian(ispin), &
8218 matrix_out=m_model_hessian_inv(ispin), &
8219 nocc=almo_scf_env%nocc_of_domain(:, ispin) &
8226 DO ispin = 1, nspins
8227 CALL dbcsr_copy(m_model_hessian_inv(ispin), &
8228 m_model_hessian(ispin))
8230 para_env=almo_scf_env%para_env, &
8231 blacs_env=almo_scf_env%blacs_env)
8233 para_env=almo_scf_env%para_env, &
8234 blacs_env=almo_scf_env%blacs_env, &
8235 uplo_to_full=.true.)
8237 almo_scf_env%eps_filter)
8242 DO ispin = 1, nspins
8244 matrix_main=m_model_hessian(ispin), &
8245 subm_s_inv=almo_scf_env%domain_s_inv(:, ispin), &
8246 subm_r_down=domain_r_down(:, ispin), &
8247 matrix_trimmer=quench_t(ispin), &
8248 dpattern=quench_t(ispin), &
8249 map=almo_scf_env%domain_map(ispin), &
8250 node_of_domain=almo_scf_env%cpu_of_domain, &
8252 use_trimmer=.false., &
8254 skip_inversion=.false. &
8262 cpabort(
"Unknown preconditioner")
8267 DO ispin = 1, nspins
8274 m_model_hessian_inv(ispin), &
8276 0.0_dp, m_model_bd(ispin), &
8277 filter_eps=almo_scf_env%eps_filter)
8282 matrix_in=m_model_r(ispin), &
8283 matrix_out=m_model_bd(ispin), &
8284 operator1=domain_model_hessian_inv(:, ispin), &
8285 dpattern=quench_t(ispin), &
8286 map=almo_scf_env%domain_map(ispin), &
8287 node_of_domain=almo_scf_env%cpu_of_domain, &
8289 filter_eps=almo_scf_env%eps_filter)
8296 CALL contravariant_matrix_norm( &
8297 norm_out=step_norm, &
8298 matrix_in=m_model_bd, &
8299 metric_in=almo_scf_env%matrix_s, &
8300 quench_t_in=quench_t, &
8301 eps_filter_in=almo_scf_env%eps_filter &
8303 IF (unit_nr > 0 .AND. debug_mode)
WRITE (unit_nr, *)
"...pB norm: ", step_norm
8306 IF (step_norm <= radius_current)
THEN
8308 IF (unit_nr > 0 .AND. debug_mode)
WRITE (unit_nr, *)
"...Full dogleg"
8310 border_reached = .false.
8312 DO ispin = 1, nspins
8313 CALL dbcsr_copy(step(ispin), m_model_bd(ispin))
8316 fake_step_size_to_report = 2.0_dp
8317 iteration_type_to_report = 6
8321 IF (unit_nr > 0 .AND. debug_mode)
WRITE (unit_nr, *)
"...pB norm is too large"
8323 border_reached = .true.
8327 DO ispin = 1, nspins
8328 CALL dbcsr_add(m_model_bd(ispin), step(ispin), 1.0_dp, -1.0_dp)
8331 CALL step_size_to_border( &
8332 step_size_out=step_size, &
8333 metric_in=almo_scf_env%matrix_s, &
8335 direction_in=m_model_bd, &
8336 trust_radius_in=radius_current, &
8337 quench_t_in=quench_t, &
8338 eps_filter_in=almo_scf_env%eps_filter &
8340 IF (unit_nr > 0 .AND. debug_mode)
WRITE (unit_nr, *)
"...Step size to border: ", step_size
8341 IF (step_size > 1.0_dp .OR. step_size < 0.0_dp)
THEN
8342 IF (unit_nr > 0)
THEN
8343 WRITE (unit_nr, *)
"Step size (", step_size,
") must lie inside (0,1)"
8345 cpabort(
"Wrong dog leg step. We should never end up here.")
8348 DO ispin = 1, nspins
8349 CALL dbcsr_add(step(ispin), m_model_bd(ispin), 1.0_dp, step_size)
8352 fake_step_size_to_report = 1.0_dp + step_size
8353 iteration_type_to_report = 7
8357 IF (debug_mode)
THEN
8359 IF (unit_nr > 0)
WRITE (unit_nr, *)
"...Extra norm evaluation"
8360 CALL contravariant_matrix_norm( &
8361 norm_out=step_norm, &
8363 metric_in=almo_scf_env%matrix_s, &
8364 quench_t_in=quench_t, &
8365 eps_filter_in=almo_scf_env%eps_filter &
8367 IF (unit_nr > 0)
WRITE (unit_nr, *)
"...Step norm: ", step_norm
8368 IF (unit_nr > 0)
WRITE (unit_nr, *)
"...Current radius: ", radius_current
8371 CALL predicted_reduction( &
8372 reduction_out=expected_reduction, &
8375 hess_in=m_model_hessian, &
8376 hess_submatrix_in=almo_scf_env%domain_preconditioner, &
8377 quench_t_in=quench_t, &
8378 special_case=my_special_case, &
8379 eps_filter=almo_scf_env%eps_filter, &
8380 domain_map=almo_scf_env%domain_map, &
8381 cpu_of_domain=almo_scf_env%cpu_of_domain &
8384 inner_loop_success = .true.
8387 CALL fixed_r_report(unit_nr, &
8388 iter_type=iteration_type_to_report, &
8389 iteration=iteration, &
8390 step_size=fake_step_size_to_report, &
8391 border_reached=border_reached, &
8392 curvature=y_scalar, &
8393 grad_norm_ratio=expected_reduction, &
8401 DO ispin = 1, nspins
8403 CALL dbcsr_copy(m_model_r_prev(ispin), m_model_r(ispin))
8404 CALL dbcsr_add(m_model_r(ispin), m_model_bd(ispin), &
8409 DO ispin = 1, nspins
8412 model_grad_norm = maxval(grad_norm_spin)
8415 grad_norm_ratio = model_grad_norm/grad_norm_ref
8416 IF (grad_norm_ratio < optimizer%model_grad_norm_ratio)
THEN
8418 border_reached = .false.
8419 inner_loop_success = .true.
8421 CALL predicted_reduction( &
8422 reduction_out=expected_reduction, &
8425 hess_in=m_model_hessian, &
8426 hess_submatrix_in=almo_scf_env%domain_preconditioner, &
8427 quench_t_in=quench_t, &
8428 special_case=my_special_case, &
8429 eps_filter=almo_scf_env%eps_filter, &
8430 domain_map=almo_scf_env%domain_map, &
8431 cpu_of_domain=almo_scf_env%cpu_of_domain &
8435 CALL fixed_r_report(unit_nr, &
8437 iteration=iteration, &
8438 step_size=step_size, &
8439 border_reached=border_reached, &
8440 curvature=y_scalar, &
8441 grad_norm_ratio=expected_reduction, &
8449 DO ispin = 1, nspins
8451 CALL dbcsr_copy(m_model_rt_prev(ispin), m_model_rt(ispin))
8454 DO ispin = 1, nspins
8462 0.0_dp, m_model_rt(ispin), &
8463 filter_eps=almo_scf_env%eps_filter)
8468 matrix_in=m_model_r(ispin), &
8469 matrix_out=m_model_rt(ispin), &
8470 operator1=almo_scf_env%domain_s_inv(:, ispin), &
8471 dpattern=quench_t(ispin), &
8472 map=almo_scf_env%domain_map(ispin), &
8473 node_of_domain=almo_scf_env%cpu_of_domain, &
8475 filter_eps=almo_scf_env%eps_filter)
8481 CALL compute_cg_beta( &
8483 reset_conjugator=reset_conjugator, &
8484 conjugator=optimizer%conjugator, &
8485 grad=m_model_r(:), &
8486 prev_grad=m_model_r_prev(:), &
8487 step=m_model_rt(:), &
8488 prev_step=m_model_rt_prev(:) &
8491 DO ispin = 1, nspins
8493 CALL dbcsr_add(m_model_d(ispin), m_model_rt(ispin), beta, 1.0_dp)
8497 CALL fixed_r_report(unit_nr, &
8499 iteration=iteration, &
8500 step_size=step_size, &
8501 border_reached=border_reached, &
8502 curvature=y_scalar, &
8503 grad_norm_ratio=grad_norm_ratio, &
8512 IF (.NOT. inner_loop_success)
THEN
8513 cpabort(
"Inner loop did not produce solution")
8516 DO ispin = 1, nspins
8518 CALL dbcsr_copy(m_theta_trial(ispin), m_theta(ispin))
8519 CALL dbcsr_add(m_theta_trial(ispin), step(ispin), 1.0_dp, 1.0_dp)
8524 CALL main_var_to_xalmos_and_loss_func( &
8525 almo_scf_env=almo_scf_env, &
8527 m_main_var_in=m_theta_trial, &
8528 m_t_out=matrix_t_out, &
8529 m_sig_sqrti_ii_out=m_sig_sqrti_ii, &
8530 energy_out=energy_trial, &
8531 penalty_out=penalty_trial, &
8532 m_ftsiginv_out=ftsiginv, &
8533 m_siginvtftsiginv_out=siginvtftsiginv, &
8535 m_stsiginv0_in=stsiginv_0, &
8536 m_quench_t_in=quench_t, &
8537 domain_r_down_in=domain_r_down, &
8538 assume_t0_q0x=assume_t0_q0x, &
8539 just_started=.false., &
8540 optimize_theta=optimize_theta, &
8541 normalize_orbitals=normalize_orbitals, &
8542 perturbation_only=perturbation_only, &
8543 do_penalty=penalty_occ_vol, &
8544 special_case=my_special_case)
8545 loss_trial = energy_trial + penalty_trial
8547 rho = (loss_trial - loss_start)/expected_reduction
8548 loss_change_to_report = loss_trial - loss_start
8550 IF (rho < 0.25_dp)
THEN
8551 radius_current = 0.25_dp*radius_current
8553 IF (rho > 0.75_dp .AND. border_reached)
THEN
8554 radius_current = min(2.0_dp*radius_current, radius_max)
8559 DO ispin = 1, nspins
8560 CALL dbcsr_copy(m_theta(ispin), m_theta_trial(ispin))
8562 loss_start = loss_trial
8563 energy_start = energy_trial
8564 penalty_start = penalty_trial
8565 same_position = .false.
8567 almo_scf_env%almo_scf_energy = energy_trial
8570 same_position = .true.
8572 almo_scf_env%almo_scf_energy = energy_start
8577 CALL trust_r_report(unit_nr, &
8579 iteration=outer_iteration, &
8581 delta_loss=loss_change_to_report, &
8583 predicted_reduction=expected_reduction, &
8585 radius=radius_current, &
8586 new=.NOT. same_position, &
8587 time=t2outer - t1outer)
8590 END DO adjust_r_loop
8593 IF (scf_converged)
THEN
8595 CALL wrap_up_xalmo_scf( &
8597 almo_scf_env=almo_scf_env, &
8598 perturbation_in=perturbation_only, &
8599 m_xalmo_in=matrix_t_out, &
8600 m_quench_in=quench_t, &
8601 energy_inout=energy_start)
8605 DO ispin = 1, nspins
8633 DEALLOCATE (m_model_hessian)
8634 DEALLOCATE (m_model_hessian_inv)
8635 DEALLOCATE (siginvtftsiginv)
8636 DEALLOCATE (stsiginv_0)
8637 DEALLOCATE (ftsiginv)
8640 DEALLOCATE (prev_step)
8642 DEALLOCATE (m_sig_sqrti_ii)
8643 DEALLOCATE (m_model_r)
8644 DEALLOCATE (m_model_rt)
8645 DEALLOCATE (m_model_d)
8646 DEALLOCATE (m_model_bd)
8647 DEALLOCATE (m_model_r_prev)
8648 DEALLOCATE (m_model_rt_prev)
8649 DEALLOCATE (m_theta_trial)
8651 DEALLOCATE (domain_r_down)
8652 DEALLOCATE (domain_model_hessian_inv)
8654 DEALLOCATE (penalty_occ_vol_g_prefactor)
8655 DEALLOCATE (penalty_occ_vol_h_prefactor)
8656 DEALLOCATE (grad_norm_spin)
8659 DEALLOCATE (m_theta)
8661 IF (.NOT. scf_converged .AND. .NOT. optimizer%early_stopping_on)
THEN
8662 cpabort(
"Optimization not converged! ")
8665 CALL timestop(handle)
8698 SUBROUTINE main_var_to_xalmos_and_loss_func(almo_scf_env, qs_env, m_main_var_in, &
8699 m_t_out, energy_out, penalty_out, m_sig_sqrti_ii_out, m_FTsiginv_out, &
8700 m_siginvTFTsiginv_out, m_ST_out, m_STsiginv0_in, m_quench_t_in, domain_r_down_in, &
8701 assume_t0_q0x, just_started, optimize_theta, normalize_orbitals, perturbation_only, &
8702 do_penalty, special_case)
8706 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(IN) :: m_main_var_in
8707 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(INOUT) :: m_t_out
8708 REAL(kind=
dp),
INTENT(OUT) :: energy_out, penalty_out
8709 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(INOUT) :: m_sig_sqrti_ii_out, m_ftsiginv_out, &
8710 m_siginvtftsiginv_out, m_st_out, &
8711 m_stsiginv0_in, m_quench_t_in
8713 INTENT(IN) :: domain_r_down_in
8714 LOGICAL,
INTENT(IN) :: assume_t0_q0x, just_started, &
8715 optimize_theta, normalize_orbitals, &
8716 perturbation_only, do_penalty
8717 INTEGER,
INTENT(IN) :: special_case
8719 CHARACTER(len=*),
PARAMETER :: routinen =
'main_var_to_xalmos_and_loss_func'
8721 INTEGER :: handle, ispin, nspins
8722 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: nocc
8723 REAL(kind=
dp) :: det1, energy_ispin, penalty_amplitude, &
8726 CALL timeset(routinen, handle)
8729 penalty_out = 0.0_dp
8731 nspins =
SIZE(m_main_var_in)
8732 IF (nspins == 1)
THEN
8733 spin_factor = 2.0_dp
8735 spin_factor = 1.0_dp
8738 penalty_amplitude = 0.0_dp
8740 ALLOCATE (nocc(nspins))
8741 DO ispin = 1, nspins
8743 nfullrows_total=nocc(ispin))
8746 DO ispin = 1, nspins
8749 CALL compute_xalmos_from_main_var( &
8750 m_var_in=m_main_var_in(ispin), &
8751 m_t_out=m_t_out(ispin), &
8752 m_quench_t=m_quench_t_in(ispin), &
8753 m_t0=almo_scf_env%matrix_t_blk(ispin), &
8754 m_oo_template=almo_scf_env%matrix_sigma_inv(ispin), &
8755 m_stsiginv0=m_stsiginv0_in(ispin), &
8756 m_s=almo_scf_env%matrix_s(1), &
8757 m_sig_sqrti_ii_out=m_sig_sqrti_ii_out(ispin), &
8758 domain_r_down=domain_r_down_in(:, ispin), &
8759 domain_s_inv=almo_scf_env%domain_s_inv(:, ispin), &
8760 domain_map=almo_scf_env%domain_map(ispin), &
8761 cpu_of_domain=almo_scf_env%cpu_of_domain, &
8762 assume_t0_q0x=assume_t0_q0x, &
8763 just_started=just_started, &
8764 optimize_theta=optimize_theta, &
8765 normalize_orbitals=normalize_orbitals, &
8766 envelope_amplitude=almo_scf_env%envelope_amplitude, &
8767 eps_filter=almo_scf_env%eps_filter, &
8768 special_case=special_case, &
8769 nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin), &
8770 order_lanczos=almo_scf_env%order_lanczos, &
8771 eps_lanczos=almo_scf_env%eps_lanczos, &
8772 max_iter_lanczos=almo_scf_env%max_iter_lanczos)
8777 p=almo_scf_env%matrix_p(ispin), &
8778 eps_filter=almo_scf_env%eps_filter, &
8779 orthog_orbs=.false., &
8780 nocc_of_domain=almo_scf_env%nocc_of_domain(:, ispin), &
8781 s=almo_scf_env%matrix_s(1), &
8782 sigma=almo_scf_env%matrix_sigma(ispin), &
8783 sigma_inv=almo_scf_env%matrix_sigma_inv(ispin), &
8784 use_guess=.false., &
8785 algorithm=almo_scf_env%sigma_inv_algorithm, &
8786 inv_eps_factor=almo_scf_env%matrix_iter_eps_error_factor, &
8787 inverse_accelerator=almo_scf_env%order_lanczos, &
8788 eps_lanczos=almo_scf_env%eps_lanczos, &
8789 max_iter_lanczos=almo_scf_env%max_iter_lanczos, &
8790 para_env=almo_scf_env%para_env, &
8791 blacs_env=almo_scf_env%blacs_env)
8800 IF (perturbation_only)
THEN
8802 IF (just_started)
THEN
8803 DO ispin = 1, nspins
8804 CALL dbcsr_copy(almo_scf_env%matrix_ks(ispin), &
8805 almo_scf_env%matrix_ks_0deloc(ispin))
8811 almo_scf_env%matrix_p, &
8812 almo_scf_env%matrix_ks, &
8814 almo_scf_env%eps_filter, &
8815 almo_scf_env%mat_distr_aos)
8818 penalty_out = 0.0_dp
8819 DO ispin = 1, nspins
8821 CALL compute_frequently_used_matrices( &
8822 filter_eps=almo_scf_env%eps_filter, &
8823 m_t_in=m_t_out(ispin), &
8824 m_siginv_in=almo_scf_env%matrix_sigma_inv(ispin), &
8825 m_s_in=almo_scf_env%matrix_s(1), &
8826 m_f_in=almo_scf_env%matrix_ks(ispin), &
8827 m_ftsiginv_out=m_ftsiginv_out(ispin), &
8828 m_siginvtftsiginv_out=m_siginvtftsiginv_out(ispin), &
8829 m_st_out=m_st_out(ispin))
8831 IF (perturbation_only)
THEN
8833 IF (ispin == 1) energy_out = 0.0_dp
8834 CALL dbcsr_dot(m_t_out(ispin), m_ftsiginv_out(ispin), energy_ispin)
8835 energy_out = energy_out + energy_ispin*spin_factor
8838 IF (do_penalty)
THEN
8840 CALL determinant(almo_scf_env%matrix_sigma(ispin), det1, &
8841 almo_scf_env%eps_filter)
8842 penalty_out = penalty_out - &
8843 penalty_amplitude*spin_factor*nocc(ispin)*log(det1)
8851 CALL timestop(handle)
8853 END SUBROUTINE main_var_to_xalmos_and_loss_func
8870 SUBROUTINE step_size_to_border(step_size_out, metric_in, position_in, &
8871 direction_in, trust_radius_in, quench_t_in, eps_filter_in)
8873 REAL(kind=
dp),
INTENT(INOUT) :: step_size_out
8874 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(IN) :: metric_in, position_in, direction_in
8875 REAL(kind=
dp),
INTENT(IN) :: trust_radius_in
8876 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(IN) :: quench_t_in
8877 REAL(kind=
dp),
INTENT(IN) :: eps_filter_in
8879 INTEGER :: isol, ispin, nsolutions, &
8880 nsolutions_found, nspins
8881 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: nocc
8882 REAL(kind=
dp) :: discrim_sign, discriminant, solution, &
8883 spin_factor, temp_real
8884 REAL(kind=
dp),
DIMENSION(3) :: coef
8885 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:) :: m_temp_no
8887 step_size_out = 0.0_dp
8889 nspins =
SIZE(position_in)
8890 IF (nspins == 1)
THEN
8891 spin_factor = 2.0_dp
8893 spin_factor = 1.0_dp
8896 ALLOCATE (nocc(nspins))
8897 ALLOCATE (m_temp_no(nspins))
8900 DO ispin = 1, nspins
8903 template=direction_in(ispin))
8906 nfullcols_total=nocc(ispin))
8908 CALL dbcsr_copy(m_temp_no(ispin), quench_t_in(ispin))
8911 position_in(ispin), &
8912 0.0_dp, m_temp_no(ispin), &
8913 retain_sparsity=.true.)
8915 CALL dbcsr_dot(position_in(ispin), m_temp_no(ispin), temp_real)
8916 coef(3) = coef(3) + temp_real/nocc(ispin)
8917 CALL dbcsr_dot(direction_in(ispin), m_temp_no(ispin), temp_real)
8918 coef(2) = coef(2) + 2.0_dp*temp_real/nocc(ispin)
8919 CALL dbcsr_copy(m_temp_no(ispin), quench_t_in(ispin))
8922 direction_in(ispin), &
8923 0.0_dp, m_temp_no(ispin), &
8924 retain_sparsity=.true.)
8926 CALL dbcsr_dot(direction_in(ispin), m_temp_no(ispin), temp_real)
8927 coef(1) = coef(1) + temp_real/nocc(ispin)
8934 DEALLOCATE (m_temp_no)
8936 coef(:) = coef(:)*spin_factor
8937 coef(3) = coef(3) - trust_radius_in*trust_radius_in
8940 discriminant = coef(2)*coef(2) - 4.0_dp*coef(1)*coef(3)
8941 IF (discriminant > tiny(discriminant))
THEN
8943 ELSE IF (discriminant < 0.0_dp)
THEN
8945 cpabort(
"Step to border: no solutions")
8950 discrim_sign = 1.0_dp
8951 nsolutions_found = 0
8952 DO isol = 1, nsolutions
8953 solution = (-coef(2) + discrim_sign*sqrt(discriminant))/(2.0_dp*coef(1))
8954 IF (solution > 0.0_dp)
THEN
8955 nsolutions_found = nsolutions_found + 1
8956 step_size_out = solution
8958 discrim_sign = -discrim_sign
8961 IF (nsolutions_found == 0)
THEN
8962 cpabort(
"Step to border: no positive solutions")
8963 ELSE IF (nsolutions_found == 2)
THEN
8964 cpabort(
"Two positive border steps possible!")
8967 END SUBROUTINE step_size_to_border
8980 SUBROUTINE contravariant_matrix_norm(norm_out, matrix_in, metric_in, &
8981 quench_t_in, eps_filter_in)
8983 REAL(kind=
dp),
INTENT(OUT) :: norm_out
8984 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(IN) :: matrix_in, metric_in, quench_t_in
8985 REAL(kind=
dp),
INTENT(IN) :: eps_filter_in
8987 INTEGER :: ispin, nspins
8988 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: nocc
8989 REAL(kind=
dp) :: my_norm, spin_factor, temp_real
8990 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:) :: m_temp_no
8995 nspins =
SIZE(matrix_in)
8996 IF (nspins == 1)
THEN
8997 spin_factor = 2.0_dp
8999 spin_factor = 1.0_dp
9002 ALLOCATE (nocc(nspins))
9003 ALLOCATE (m_temp_no(nspins))
9006 DO ispin = 1, nspins
9008 CALL dbcsr_create(m_temp_no(ispin), template=matrix_in(ispin))
9011 nfullcols_total=nocc(ispin))
9013 CALL dbcsr_copy(m_temp_no(ispin), quench_t_in(ispin))
9017 0.0_dp, m_temp_no(ispin), &
9018 retain_sparsity=.true.)
9020 CALL dbcsr_dot(matrix_in(ispin), m_temp_no(ispin), temp_real)
9022 my_norm = my_norm + temp_real/nocc(ispin)
9029 DEALLOCATE (m_temp_no)
9031 my_norm = my_norm*spin_factor
9032 norm_out = sqrt(my_norm)
9034 END SUBROUTINE contravariant_matrix_norm
9053 SUBROUTINE predicted_reduction(reduction_out, grad_in, step_in, hess_in, &
9054 hess_submatrix_in, quench_t_in, special_case, eps_filter, domain_map, &
9058 REAL(kind=
dp),
INTENT(INOUT) :: reduction_out
9059 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(INOUT) :: grad_in, step_in, hess_in
9061 INTENT(IN) :: hess_submatrix_in
9062 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(IN) :: quench_t_in
9063 INTEGER,
INTENT(IN) :: special_case
9064 REAL(kind=
dp),
INTENT(IN) :: eps_filter
9066 INTEGER,
DIMENSION(:),
INTENT(IN) :: cpu_of_domain
9068 INTEGER :: ispin, nspins
9069 REAL(kind=
dp) :: my_reduction, spin_factor, temp_real
9070 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:) :: m_temp_no
9072 reduction_out = 0.0_dp
9074 nspins =
SIZE(grad_in)
9075 IF (nspins == 1)
THEN
9076 spin_factor = 2.0_dp
9078 spin_factor = 1.0_dp
9081 ALLOCATE (m_temp_no(nspins))
9083 my_reduction = 0.0_dp
9084 DO ispin = 1, nspins
9086 CALL dbcsr_create(m_temp_no(ispin), template=grad_in(ispin))
9088 CALL dbcsr_dot(step_in(ispin), grad_in(ispin), temp_real)
9089 my_reduction = my_reduction + temp_real
9098 0.0_dp, m_temp_no(ispin), &
9099 filter_eps=eps_filter)
9104 matrix_in=step_in(ispin), &
9105 matrix_out=m_temp_no(ispin), &
9106 operator1=hess_submatrix_in(:, ispin), &
9107 dpattern=quench_t_in(ispin), &
9108 map=domain_map(ispin), &
9109 node_of_domain=cpu_of_domain, &
9111 filter_eps=eps_filter)
9116 CALL dbcsr_dot(step_in(ispin), m_temp_no(ispin), temp_real)
9117 my_reduction = my_reduction + 0.5_dp*temp_real
9124 my_reduction = spin_factor*my_reduction
9126 reduction_out = my_reduction
9128 DEALLOCATE (m_temp_no)
9130 END SUBROUTINE predicted_reduction
9147 SUBROUTINE fixed_r_report(unit_nr, iter_type, iteration, step_size, &
9148 border_reached, curvature, grad_norm_ratio, predicted_reduction, time)
9150 INTEGER,
INTENT(IN) :: unit_nr, iter_type, iteration
9151 REAL(kind=
dp),
INTENT(IN) :: step_size
9152 LOGICAL,
INTENT(IN) :: border_reached
9153 REAL(kind=
dp),
INTENT(IN) :: curvature
9154 REAL(kind=
dp),
INTENT(IN),
OPTIONAL :: grad_norm_ratio, predicted_reduction
9155 REAL(kind=
dp),
INTENT(IN) :: time
9157 CHARACTER(LEN=20) :: iter_type_str
9158 REAL(kind=
dp) :: loss_or_grad_change
9160 loss_or_grad_change = 0.0_dp
9161 IF (
PRESENT(grad_norm_ratio))
THEN
9162 loss_or_grad_change = grad_norm_ratio
9163 ELSE IF (
PRESENT(predicted_reduction))
THEN
9164 loss_or_grad_change = predicted_reduction
9166 cpabort(
"one argument is missing")
9169 SELECT CASE (iter_type)
9171 iter_type_str = trim(
"Ignored")
9173 iter_type_str = trim(
"PCG")
9175 iter_type_str = trim(
"Neg. curvatr.")
9177 iter_type_str = trim(
"Step too long")
9179 iter_type_str = trim(
"Grad. reduced")
9181 iter_type_str = trim(
"Cauchy point")
9183 iter_type_str = trim(
"Full dogleg")
9185 iter_type_str = trim(
"Part. dogleg")
9187 cpabort(
"unknown report type")
9190 IF (unit_nr > 0)
THEN
9192 SELECT CASE (iter_type)
9196 WRITE (unit_nr,
'(T4,A15,A6,A10,A10,A7,A20,A8)') &
9202 "Grad/o.f. reduc", &
9207 WRITE (unit_nr,
'(T4,A15,I6,F10.5,F10.5,L7,F20.10,F8.2)') &
9210 curvature, step_size, border_reached, &
9211 loss_or_grad_change, &
9217 SELECT CASE (iter_type)
9218 CASE (2, 3, 4, 5, 6, 7)
9226 END SUBROUTINE fixed_r_report
9245 SUBROUTINE trust_r_report(unit_nr, iter_type, iteration, radius, &
9246 loss, delta_loss, grad_norm, predicted_reduction, rho, new, time)
9248 INTEGER,
INTENT(IN) :: unit_nr, iter_type, iteration
9249 REAL(kind=
dp),
INTENT(IN) :: radius, loss, delta_loss, grad_norm, &
9250 predicted_reduction, rho
9251 LOGICAL,
INTENT(IN) :: new
9252 REAL(kind=
dp),
INTENT(IN) :: time
9254 CHARACTER(LEN=20) :: iter_status, iter_type_str
9256 SELECT CASE (iter_type)
9258 iter_type_str = trim(
"Iter")
9259 iter_status = trim(
"Stat")
9261 iter_type_str = trim(
"TR INI")
9263 iter_status =
" New"
9265 iter_status =
" Redo"
9268 iter_type_str = trim(
"TR FIN")
9270 iter_status =
" Acc"
9272 iter_status =
" Rej"
9275 cpabort(
"unknown report type")
9278 IF (unit_nr > 0)
THEN
9280 SELECT CASE (iter_type)
9283 WRITE (unit_nr,
'(T2,A6,A5,A6,A22,A10,T67,A7,A6)') &
9287 "Objective Function", &
9291 WRITE (unit_nr,
'(T41,A10,A10,A6)') &
9295 "Change",
"Expct.",
"Rho"
9301 WRITE (unit_nr,
'(T2,A6,A5,I6,F22.10,ES10.2,T67,ES7.0,F6.1)') &
9312 WRITE (unit_nr,
'(T2,A6,A5,I6,F22.10,ES10.2,ES10.2,F6.1,ES7.0,F6.1)') &
9317 delta_loss, predicted_reduction, rho, &
9324 END SUBROUTINE trust_r_report
9332 SUBROUTINE energy_lowering_report(unit_nr, ref_energy, energy_lowering)
9334 INTEGER,
INTENT(IN) :: unit_nr
9335 REAL(kind=
dp),
INTENT(IN) :: ref_energy, energy_lowering
9338 IF (unit_nr > 0)
THEN
9340 WRITE (unit_nr,
'(T2,A35,F25.10)')
"ENERGY OF BLOCK-DIAGONAL ALMOs:", &
9342 WRITE (unit_nr,
'(T2,A35,F25.10)')
"ENERGY LOWERING:", &
9344 WRITE (unit_nr,
'(T2,A35,F25.10)')
"CORRECTED ENERGY:", &
9345 ref_energy + energy_lowering
9349 END SUBROUTINE energy_lowering_report
9361 SUBROUTINE wrap_up_xalmo_scf(qs_env, almo_scf_env, perturbation_in, &
9362 m_xalmo_in, m_quench_in, energy_inout)
9366 LOGICAL,
INTENT(IN) :: perturbation_in
9367 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(IN) :: m_xalmo_in, m_quench_in
9368 REAL(kind=
dp),
INTENT(INOUT) :: energy_inout
9370 CHARACTER(len=*),
PARAMETER :: routinen =
'wrap_up_xalmo_scf'
9372 INTEGER :: eda_unit, handle, ispin, nspins, unit_nr
9374 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:) :: m_temp_no1, m_temp_no2
9377 CALL timeset(routinen, handle)
9381 IF (logger%para_env%is_source())
THEN
9387 nspins = almo_scf_env%nspins
9391 IF (perturbation_in)
THEN
9393 ALLOCATE (m_temp_no1(nspins))
9394 ALLOCATE (m_temp_no2(nspins))
9396 DO ispin = 1, nspins
9397 CALL dbcsr_create(m_temp_no1(ispin), template=m_xalmo_in(ispin))
9398 CALL dbcsr_create(m_temp_no2(ispin), template=m_xalmo_in(ispin))
9403 almo_scf_env%mat_distr_aos)
9408 CALL xalmo_analysis( &
9409 detailed_analysis=almo_scf_env%almo_analysis%do_analysis, &
9410 eps_filter=almo_scf_env%eps_filter, &
9411 m_t_in=m_xalmo_in, &
9412 m_t0_in=almo_scf_env%matrix_t_blk, &
9413 m_siginv_in=almo_scf_env%matrix_sigma_inv, &
9414 m_siginv0_in=almo_scf_env%matrix_sigma_inv_0deloc, &
9415 m_s_in=almo_scf_env%matrix_s, &
9416 m_ks0_in=almo_scf_env%matrix_ks_0deloc, &
9417 m_quench_t_in=m_quench_in, &
9418 energy_out=energy_inout, &
9419 m_eda_out=m_temp_no1, &
9420 m_cta_out=m_temp_no2 &
9423 IF (almo_scf_env%almo_analysis%do_analysis)
THEN
9425 DO ispin = 1, nspins
9428 IF (unit_nr > 0)
THEN
9429 WRITE (unit_nr,
'(T2,A)')
"DECOMPOSITION OF THE DELOCALIZATION ENERGY"
9436 "ALMO_EDA_CT", extension=
".dat", local=.true.)
9437 CALL print_block_sum(m_temp_no1(ispin), eda_unit)
9439 "ALMO_EDA_CT", local=.true.)
9442 IF (unit_nr > 0)
THEN
9443 WRITE (unit_nr,
'(T2,A)')
"DECOMPOSITION OF CHARGE TRANSFER TERMS"
9447 "ALMO_CTA", extension=
".dat", local=.true.)
9448 CALL print_block_sum(m_temp_no2(ispin), eda_unit)
9450 "ALMO_CTA", local=.true.)
9456 CALL energy_lowering_report( &
9458 ref_energy=almo_scf_env%almo_scf_energy, &
9459 energy_lowering=energy_inout)
9461 energy=almo_scf_env%almo_scf_energy, &
9462 energy_singles_corr=energy_inout)
9464 DO ispin = 1, nspins
9469 DEALLOCATE (m_temp_no1)
9470 DEALLOCATE (m_temp_no2)
9475 energy=energy_inout)
9479 CALL timestop(handle)
9481 END SUBROUTINE wrap_up_xalmo_scf
9489 SUBROUTINE tanh_of_elements(matrix, alpha)
9491 REAL(kind=
dp),
INTENT(IN) :: alpha
9493 CHARACTER(len=*),
PARAMETER :: routinen =
'tanh_of_elements'
9496 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: block
9499 CALL timeset(routinen, handle)
9503 block = tanh(alpha*block)
9506 CALL timestop(handle)
9508 END SUBROUTINE tanh_of_elements
9516 SUBROUTINE dtanh_of_elements(matrix, alpha)
9518 REAL(kind=
dp),
INTENT(IN) :: alpha
9520 CHARACTER(len=*),
PARAMETER :: routinen =
'dtanh_of_elements'
9523 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: block
9526 CALL timeset(routinen, handle)
9530 block = alpha*(1.0_dp - tanh(block)**2)
9533 CALL timestop(handle)
9535 END SUBROUTINE dtanh_of_elements
9542 SUBROUTINE inverse_of_elements(matrix)
9545 CHARACTER(len=*),
PARAMETER :: routinen =
'inverse_of_elements'
9548 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: block
9551 CALL timeset(routinen, handle)
9555 block = 1.0_dp/block
9558 CALL timestop(handle)
9560 END SUBROUTINE inverse_of_elements
9567 SUBROUTINE print_block_sum(matrix, unit_nr)
9569 INTEGER,
INTENT(IN) :: unit_nr
9571 CHARACTER(len=*),
PARAMETER :: routinen =
'print_block_sum'
9573 INTEGER :: col, handle, row
9574 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: block
9577 CALL timeset(routinen, handle)
9579 IF (unit_nr > 0)
THEN
9583 WRITE (unit_nr,
'(I6,I6,ES18.9)') row, col, sum(block)
9588 CALL timestop(handle)
9589 END SUBROUTINE print_block_sum
A DIIS implementation for the ALMO-based SCF methods.
subroutine, public almo_scf_diis_release(diis_env)
destroys the diis structure
subroutine, public almo_scf_diis_extrapolate(diis_env, extr_var, d_extr_var)
extrapolates the variable using the saved history
subroutine, public almo_scf_diis_push(diis_env, var, err, d_var, d_err)
adds a variable-error pair to the diis structure
subroutine, public lbfgs_create(history, nspins, nstore)
create history storage for limited memory bfgs
subroutine, public lbfgs_seed(history, variable, gradient)
interface subroutine to store the first variable/gradient pair
subroutine, public lbfgs_release(history)
release the bfgs history
subroutine, public lbfgs_get_direction(history, variable, gradient, direction)
interface subroutine to store a variable/gradient pair and predict direction
Subroutines for ALMO SCF.
subroutine, public construct_domain_preconditioner(matrix_main, subm_s_inv, subm_s_inv_half, subm_s_half, subm_r_down, matrix_trimmer, dpattern, map, node_of_domain, preconditioner, bad_modes_projector_down, use_trimmer, eps_zero_eigenvalues, my_action, skip_inversion)
Constructs preconditioners for each domain -1. projected preconditioner 0. simple preconditioner.
subroutine, public almo_scf_ks_xx_to_tv_xx(almo_scf_env)
ALMOs by diagonalizing the KS domain submatrices computes both the occupied and virtual orbitals.
subroutine, public xalmo_initial_guess(m_guess, m_t_in, m_t0, m_quench_t, m_overlap, m_sigma_tmpl, nspins, xalmo_history, assume_t0_q0x, optimize_theta, envelope_amplitude, eps_filter, order_lanczos, eps_lanczos, max_iter_lanczos, nocc_of_domain)
create the initial guess for XALMOs
subroutine, public almo_scf_p_blk_to_t_blk(almo_scf_env, ionic)
computes occupied ALMOs from the superimposed atomic density blocks
subroutine, public pseudo_invert_diagonal_blk(matrix_in, matrix_out, nocc)
inverts block-diagonal blocks of a dbcsr_matrix
subroutine, public almo_scf_ks_blk_to_tv_blk(almo_scf_env)
computes ALMOs by diagonalizing the projected blocked KS matrix uses the diagonalization code for blo...
subroutine, public apply_domain_operators(matrix_in, matrix_out, operator1, operator2, dpattern, map, node_of_domain, my_action, filter_eps, matrix_trimmer, use_trimmer)
Parallel code for domain specific operations (my_action) 0. out = op1 * in.
subroutine, public construct_domain_r_down(matrix_t, matrix_sigma_inv, matrix_s, subm_r_down, dpattern, map, node_of_domain, filter_eps)
Constructs subblocks of the covariant-covariant projectors (i.e. DM without spin factor)
subroutine, public almo_scf_t_to_proj(t, p, eps_filter, orthog_orbs, nocc_of_domain, s, sigma, sigma_inv, use_guess, smear, algorithm, para_env, blacs_env, eps_lanczos, max_iter_lanczos, inverse_accelerator, inv_eps_factor)
computes the idempotent density matrix from MOs MOs can be either orthogonal or non-orthogonal
subroutine, public construct_domain_s_inv(matrix_s, subm_s_inv, dpattern, map, node_of_domain)
Constructs S_inv block for each domain.
subroutine, public almo_scf_ks_to_ks_blk(almo_scf_env)
computes the projected KS from the total KS matrix also computes the DIIS error vector as a by-produc...
subroutine, public get_overlap(bra, ket, overlap, metric, retain_overlap_sparsity, eps_filter, smear)
Computes the overlap matrix of MO orbitals.
subroutine, public fill_matrix_with_ones(matrix)
Fill all matrix blocks with 1.0_dp.
subroutine, public apply_projector(psi_in, psi_out, psi_projector, metric, project_out, psi_projector_orthogonal, proj_in_template, eps_filter, sig_inv_projector, sig_inv_template)
applies projector to the orbitals |psi_out> = P |psi_in> OR |psi_out> = (1-P) |psi_in>,...
subroutine, public construct_domain_s_sqrt(matrix_s, subm_s_sqrt, subm_s_sqrt_inv, dpattern, map, node_of_domain)
Constructs S^(+1/2) and S^(-1/2) submatrices for each domain.
subroutine, public orthogonalize_mos(ket, overlap, metric, retain_locality, only_normalize, nocc_of_domain, eps_filter, order_lanczos, eps_lanczos, max_iter_lanczos, overlap_sqrti, smear)
orthogonalize MOs
subroutine, public almo_scf_ks_to_ks_xx(almo_scf_env)
builds projected KS matrices for the overlapping domains also computes the DIIS error vector as a by-...
subroutine, public almo_scf_t_rescaling(matrix_t, mo_energies, mu_of_domain, real_ne_of_domain, spin_kts, smear_e_temp, ndomains, nocc_of_domain)
Apply an occupation-rescaling trick to ALMOs for smearing. Partially occupied orbitals are considered...
Optimization routines for all ALMO-based SCF methods.
subroutine, public almo_scf_xalmo_trustr(qs_env, almo_scf_env, optimizer, quench_t, matrix_t_in, matrix_t_out, perturbation_only, special_case)
Optimization of ALMOs using trust region minimizers.
subroutine, public almo_scf_xalmo_pcg(qs_env, almo_scf_env, optimizer, quench_t, matrix_t_in, matrix_t_out, assume_t0_q0x, perturbation_only, special_case)
Optimization of ALMOs using PCG-like minimizers.
subroutine, public almo_scf_xalmo_eigensolver(qs_env, almo_scf_env, optimizer)
An eigensolver-based SCF to optimize extended ALMOs (i.e. ALMOs on overlapping domains)
subroutine, public almo_scf_construct_nlmos(qs_env, optimizer, matrix_s, matrix_mo_in, matrix_mo_out, template_matrix_sigma, overlap_determinant, mat_distr_aos, virtuals, eps_filter)
Optimization of NLMOs using PCG minimizers.
subroutine, public almo_scf_block_diagonal(qs_env, almo_scf_env, optimizer)
An SCF procedure that optimizes block-diagonal ALMOs using DIIS.
Interface between ALMO SCF and QS.
subroutine, public almo_scf_update_ks_energy(qs_env, energy, energy_singles_corr)
update qs_env total energy
subroutine, public almo_dm_to_almo_ks(qs_env, matrix_p, matrix_ks, energy_total, eps_filter, mat_distr_aos, smear, kts_sum)
uses the ALMO density matrix to compute ALMO KS matrix and the new energy
subroutine, public almo_dm_to_qs_env(qs_env, matrix_p, mat_distr_aos)
return density matrix to the qs_env
subroutine, public matrix_qs_to_almo(matrix_qs, matrix_almo, mat_distr_aos)
convert between two types of matrices: QS style to ALMO style
Types for all ALMO-based methods.
Handles all functions related to the CELL.
methods related to the blacs parallel environment
subroutine, public dbcsr_scale(matrix, alpha_scalar)
...
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 a*x**3 + b*x**2 + c*x + d = 0 and returns only those roots for...
subroutine, public diagonalize_diagonal_blocks(matrix, c, e)
Diagonalizes diagonal blocks of a symmetric dbcsr matrix and returs its eigenvectors.
subroutine, public ct_step_execute(cts_env)
Performs Cayley transformation.
Types for all cayley transformation methods.
subroutine, public ct_step_env_clean(env)
...
subroutine, public ct_step_env_set(env, para_env, blacs_env, use_occ_orbs, use_virt_orbs, tensor_type, occ_orbs_orthogonal, virt_orbs_orthogonal, neglect_quadratic_term, update_p, update_q, eps_convergence, eps_filter, max_iter, p_index_up, p_index_down, q_index_up, q_index_down, matrix_ks, matrix_p, matrix_qp_template, matrix_pq_template, matrix_t, matrix_v, matrix_x_guess, calculate_energy_corr, conjugator, qq_preconditioner_full, pp_preconditioner_full)
...
subroutine, public ct_step_env_init(env)
...
subroutine, public ct_step_env_get(env, use_occ_orbs, use_virt_orbs, tensor_type, occ_orbs_orthogonal, virt_orbs_orthogonal, neglect_quadratic_term, update_p, update_q, eps_convergence, eps_filter, max_iter, p_index_up, p_index_down, q_index_up, q_index_down, matrix_ks, matrix_p, matrix_qp_template, matrix_pq_template, matrix_t, matrix_v, copy_matrix_x, energy_correction, calculate_energy_corr, converged, qq_preconditioner_full, pp_preconditioner_full)
...
Subroutines to handle submatrices.
subroutine, public maxnorm_submatrices(submatrices, norm)
Computes the max norm of the collection of submatrices.
subroutine, public construct_submatrices(matrix, submatrix, distr_pattern, domain_map, node_of_domain, job_type)
Constructs submatrices for each ALMO domain by collecting distributed DBCSR blocks to local arrays.
Types to handle submatrices.
integer, parameter, public select_row
Routines useful for iterative matrix calculations.
recursive subroutine, public determinant(matrix, det, threshold)
Computes the determinant of a symmetric positive definite matrix using the trace of the matrix logari...
subroutine, public invert_hotelling(matrix_inverse, matrix, threshold, use_inv_as_guess, norm_convergence, filter_eps, accelerator_order, max_iter_lanczos, eps_lanczos, silent)
invert a symmetric positive definite matrix by Hotelling's method explicit symmetrization makes this ...
subroutine, public matrix_sqrt_newton_schulz(matrix_sqrt, matrix_sqrt_inv, matrix, threshold, order, eps_lanczos, max_iter_lanczos, symmetrize, converged, iounit)
compute the sqrt of a matrix via the sign function and the corresponding Newton-Schulz iterations the...
Defines the basic variable types.
integer, parameter, public dp
Machine interface based on Fortran 2003 and POSIX.
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
Interface to the message passing library MPI.
Define methods related to particle_type.
subroutine, public get_particle_set(particle_set, qs_kind_set, first_sgf, last_sgf, nsgf, nmao, basis, 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 Berry operator for periodic systems used to define the spread of the MOS Here the matrix...
Localization methods such as 2x2 Jacobi rotations Steepest Decents Conjugate Gradient.
subroutine, public initialize_weights(cell, weights)
...
Type defining parameters related to the simulation cell.
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
type of a logger, at the moment it contains just a print level starting at which level it should be l...
stores all the informations relevant to an mpi environment
Provides all information about a quickstep kind.