43 USE dbcsr_api,
ONLY: &
44 dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_desymmetrize, &
45 dbcsr_filter, dbcsr_frobenius_norm, dbcsr_get_block_p, dbcsr_init_p, &
46 dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, dbcsr_iterator_start, &
47 dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, dbcsr_p_type, dbcsr_release, &
48 dbcsr_scale, dbcsr_set, dbcsr_transposed, dbcsr_type, dbcsr_type_antisymmetric
80 #include "../base/base_uses.f90"
86 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'rt_propagation_methods'
106 TYPE(qs_environment_type),
POINTER :: qs_env
107 TYPE(rt_prop_type),
POINTER :: rtp
108 TYPE(rtp_control_type),
POINTER :: rtp_control
110 CHARACTER(len=*),
PARAMETER :: routinen =
'propagation_step'
112 INTEGER :: aspc_order, handle, i, im, re, unit_nr
113 TYPE(cell_type),
POINTER :: cell
114 TYPE(cp_fm_type),
DIMENSION(:),
POINTER :: delta_mos, mos_new
115 TYPE(cp_logger_type),
POINTER :: logger
116 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: delta_p, h_last_iter, ks_mix, ks_mix_im, &
117 matrix_ks, matrix_ks_im, matrix_s, &
119 TYPE(dft_control_type),
POINTER :: dft_control
121 CALL timeset(routinen, handle)
124 IF (logger%para_env%is_source())
THEN
130 NULLIFY (cell, delta_p, rho_new, delta_mos, mos_new)
131 NULLIFY (ks_mix, ks_mix_im)
133 CALL get_qs_env(qs_env, cell=cell, matrix_s=matrix_s, dft_control=dft_control)
135 IF (rtp%iter == 1)
THEN
138 IF (rtp_control%fixed_ions)
CALL set_ks_env(qs_env%ks_env, s_mstruct_changed=.false.)
143 IF (dft_control%apply_efield_field)
THEN
144 IF (any(cell%perd(1:3) /= 0)) &
145 cpabort(
"Length gauge (efield) and periodicity are not compatible")
147 ELSE IF (rtp_control%velocity_gauge)
THEN
148 IF (dft_control%apply_vector_potential) &
154 IF (.NOT. rtp_control%fixed_ions)
THEN
157 rtp%delta_iter = 100.0_dp
158 rtp%mixing_factor = 1.0_dp
160 aspc_order = rtp_control%aspc_order
161 CALL aspc_extrapolate(rtp, matrix_s, aspc_order)
162 IF (rtp%linear_scaling)
THEN
169 IF (.NOT. rtp_control%fixed_ions)
THEN
172 rtp%converged = .false.
174 IF (rtp%linear_scaling)
THEN
176 CALL get_rtp(rtp=rtp, rho_new=rho_new)
179 DO i = 1,
SIZE(rho_new)
180 CALL dbcsr_init_p(delta_p(i)%matrix)
181 CALL dbcsr_create(delta_p(i)%matrix, template=rho_new(i)%matrix)
182 CALL dbcsr_copy(delta_p(i)%matrix, rho_new(i)%matrix)
186 CALL get_rtp(rtp=rtp, mos_new=mos_new)
187 ALLOCATE (delta_mos(
SIZE(mos_new)))
188 DO i = 1,
SIZE(mos_new)
190 matrix_struct=mos_new(i)%matrix_struct, &
191 name=
"delta_mos"//trim(adjustl(cp_to_string(i))))
192 CALL cp_fm_to_fm(mos_new(i), delta_mos(i))
197 matrix_ks=matrix_ks, &
198 matrix_ks_im=matrix_ks_im)
200 CALL get_rtp(rtp=rtp, h_last_iter=h_last_iter)
202 IF (unit_nr > 0)
THEN
203 WRITE (unit_nr,
'(t3,a,2f16.8)')
"Mixing the Hamiltonians to improve robustness, mixing factor: ", rtp%mixing_factor
207 DO i = 1,
SIZE(matrix_ks)
208 CALL dbcsr_init_p(ks_mix(i)%matrix)
209 CALL dbcsr_create(ks_mix(i)%matrix, template=matrix_ks(1)%matrix)
210 CALL dbcsr_init_p(ks_mix_im(i)%matrix)
211 CALL dbcsr_create(ks_mix_im(i)%matrix, template=matrix_ks(1)%matrix, matrix_type=dbcsr_type_antisymmetric)
213 DO i = 1,
SIZE(matrix_ks)
216 CALL dbcsr_add(ks_mix(i)%matrix, matrix_ks(i)%matrix, 0.0_dp, rtp%mixing_factor)
217 CALL dbcsr_add(ks_mix(i)%matrix, h_last_iter(re)%matrix, 1.0_dp, 1.0_dp - rtp%mixing_factor)
218 IF (rtp%propagate_complex_ks)
THEN
219 CALL dbcsr_add(ks_mix_im(i)%matrix, matrix_ks_im(i)%matrix, 0.0_dp, rtp%mixing_factor)
220 CALL dbcsr_add(ks_mix_im(i)%matrix, h_last_iter(im)%matrix, 1.0_dp, 1.0_dp - rtp%mixing_factor)
223 CALL calc_sinvh(rtp, ks_mix, ks_mix_im, rtp_control)
224 DO i = 1,
SIZE(matrix_ks)
227 CALL dbcsr_copy(h_last_iter(re)%matrix, ks_mix(i)%matrix)
228 IF (rtp%propagate_complex_ks)
THEN
229 CALL dbcsr_copy(h_last_iter(im)%matrix, ks_mix_im(i)%matrix)
235 CALL calc_sinvh(rtp, matrix_ks, matrix_ks_im, rtp_control)
236 DO i = 1,
SIZE(matrix_ks)
239 CALL dbcsr_copy(h_last_iter(re)%matrix, matrix_ks(i)%matrix)
240 IF (rtp%propagate_complex_ks)
THEN
241 CALL dbcsr_copy(h_last_iter(im)%matrix, matrix_ks_im(i)%matrix)
246 CALL compute_propagator_matrix(rtp, rtp_control%propagator)
248 SELECT CASE (rtp_control%mat_exp)
250 IF (rtp%linear_scaling)
THEN
264 CALL step_finalize(qs_env, rtp_control, delta_mos, delta_p)
265 IF (rtp%linear_scaling)
THEN
268 CALL cp_fm_release(delta_mos)
271 CALL timestop(handle)
287 SUBROUTINE step_finalize(qs_env, rtp_control, delta_mos, delta_P)
288 TYPE(qs_environment_type),
POINTER :: qs_env
289 TYPE(rtp_control_type),
POINTER :: rtp_control
290 TYPE(cp_fm_type),
DIMENSION(:),
POINTER :: delta_mos
291 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: delta_p
293 CHARACTER(len=*),
PARAMETER :: routinen =
'step_finalize'
295 INTEGER :: handle, i, ihist
296 TYPE(cp_fm_type),
DIMENSION(:),
POINTER :: mos_new, mos_old
297 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: exp_h_new, exp_h_old, matrix_ks, &
298 matrix_ks_im, rho_new, rho_old, s_mat
299 TYPE(qs_energy_type),
POINTER :: energy
300 TYPE(rt_prop_type),
POINTER :: rtp
302 CALL timeset(routinen, handle)
304 CALL get_qs_env(qs_env=qs_env, rtp=rtp, matrix_s=s_mat, &
305 matrix_ks=matrix_ks, matrix_ks_im=matrix_ks_im, energy=energy)
306 CALL get_rtp(rtp=rtp, exp_h_old=exp_h_old, exp_h_new=exp_h_new)
308 IF (rtp_control%sc_check_start .LT. rtp%iter)
THEN
309 rtp%delta_iter_old = rtp%delta_iter
310 IF (rtp%linear_scaling)
THEN
313 CALL rt_convergence(rtp, s_mat(1)%matrix, delta_mos, rtp%delta_iter)
315 rtp%converged = (rtp%delta_iter .LT. rtp_control%eps_ener)
320 IF (rtp_control%sc_check_start .LT. rtp%iter + 1)
THEN
321 IF (rtp%delta_iter/rtp%delta_iter_old > 0.9)
THEN
322 rtp%mixing_factor = max(rtp%mixing_factor/2.0_dp, 0.125_dp)
328 IF (rtp%converged)
THEN
329 IF (rtp%linear_scaling)
THEN
330 CALL get_rtp(rtp=rtp, rho_old=rho_old, rho_new=rho_new)
331 CALL purify_mcweeny_complex_nonorth(rho_new, s_mat, rtp%filter_eps, rtp%filter_eps_small, &
332 rtp_control%mcweeny_max_iter, rtp_control%mcweeny_eps)
335 DO i = 1,
SIZE(rho_new)
336 CALL dbcsr_copy(rho_old(i)%matrix, rho_new(i)%matrix)
339 CALL get_rtp(rtp=rtp, mos_old=mos_old, mos_new=mos_new)
340 DO i = 1,
SIZE(mos_new)
341 CALL cp_fm_to_fm(mos_new(i), mos_old(i))
344 IF (rtp_control%propagator ==
do_em)
CALL calc_sinvh(rtp, matrix_ks, matrix_ks_im, rtp_control)
345 DO i = 1,
SIZE(exp_h_new)
346 CALL dbcsr_copy(exp_h_old(i)%matrix, exp_h_new(i)%matrix)
348 ihist = mod(rtp%istep, rtp_control%aspc_order) + 1
349 IF (rtp_control%fixed_ions)
THEN
356 rtp%energy_new = energy%total
358 CALL timestop(handle)
360 END SUBROUTINE step_finalize
369 SUBROUTINE compute_propagator_matrix(rtp, propagator)
370 TYPE(rt_prop_type),
POINTER :: rtp
371 INTEGER :: propagator
373 CHARACTER(len=*),
PARAMETER :: routinen =
'compute_propagator_matrix'
376 REAL(kind=
dp) :: dt, prefac
377 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: exp_h_new, exp_h_old, propagator_matrix
379 CALL timeset(routinen, handle)
380 CALL get_rtp(rtp=rtp, exp_h_new=exp_h_new, exp_h_old=exp_h_old, &
381 propagator_matrix=propagator_matrix, dt=dt)
385 DO i = 1,
SIZE(exp_h_new)
386 CALL dbcsr_add(propagator_matrix(i)%matrix, exp_h_new(i)%matrix, 0.0_dp, prefac)
387 IF (propagator ==
do_em) &
388 CALL dbcsr_add(propagator_matrix(i)%matrix, exp_h_old(i)%matrix, 1.0_dp, prefac)
391 CALL timestop(handle)
393 END SUBROUTINE compute_propagator_matrix
405 SUBROUTINE calc_sinvh(rtp, matrix_ks, matrix_ks_im, rtp_control)
406 TYPE(rt_prop_type),
POINTER :: rtp
407 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_ks, matrix_ks_im
408 TYPE(rtp_control_type),
POINTER :: rtp_control
410 CHARACTER(len=*),
PARAMETER :: routinen =
'calc_SinvH'
411 REAL(kind=
dp),
PARAMETER :: one = 1.0_dp, zero = 0.0_dp
413 INTEGER :: handle, im, ispin, re
414 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: exp_h, sinvb, sinvh, sinvh_imag
415 TYPE(dbcsr_type) :: matrix_ks_nosym
416 TYPE(dbcsr_type),
POINTER :: b_mat, s_inv
418 CALL timeset(routinen, handle)
419 CALL get_rtp(rtp=rtp, s_inv=s_inv, exp_h_new=exp_h)
420 DO ispin = 1,
SIZE(matrix_ks)
423 CALL dbcsr_set(exp_h(re)%matrix, zero)
424 CALL dbcsr_set(exp_h(im)%matrix, zero)
426 CALL dbcsr_create(matrix_ks_nosym, template=matrix_ks(1)%matrix, matrix_type=
"N")
429 DO ispin = 1,
SIZE(matrix_ks)
432 CALL dbcsr_desymmetrize(matrix_ks(ispin)%matrix, matrix_ks_nosym)
433 CALL dbcsr_multiply(
"N",
"N", one, s_inv, matrix_ks_nosym, zero, exp_h(im)%matrix, &
434 filter_eps=rtp%filter_eps)
435 IF (.NOT. rtp_control%fixed_ions)
THEN
436 CALL get_rtp(rtp=rtp, sinvh=sinvh)
437 CALL dbcsr_copy(sinvh(ispin)%matrix, exp_h(im)%matrix)
442 IF (rtp%propagate_complex_ks)
THEN
443 DO ispin = 1,
SIZE(matrix_ks)
446 CALL dbcsr_set(matrix_ks_nosym, 0.0_dp)
447 CALL dbcsr_desymmetrize(matrix_ks_im(ispin)%matrix, matrix_ks_nosym)
449 CALL dbcsr_multiply(
"N",
"N", -one, s_inv, matrix_ks_nosym, zero, exp_h(re)%matrix, &
450 filter_eps=rtp%filter_eps)
451 IF (.NOT. rtp_control%fixed_ions)
THEN
452 CALL get_rtp(rtp=rtp, sinvh_imag=sinvh_imag)
454 CALL dbcsr_copy(sinvh_imag(ispin)%matrix, exp_h(re)%matrix)
459 IF (.NOT. rtp_control%fixed_ions)
THEN
460 CALL get_rtp(rtp=rtp, b_mat=b_mat, sinvb=sinvb)
461 CALL dbcsr_set(matrix_ks_nosym, 0.0_dp)
462 CALL dbcsr_multiply(
"N",
"N", one, s_inv, b_mat, zero, matrix_ks_nosym, filter_eps=rtp%filter_eps)
463 DO ispin = 1,
SIZE(matrix_ks)
467 CALL dbcsr_add(exp_h(re)%matrix, matrix_ks_nosym, 1.0_dp, 1.0_dp)
469 CALL dbcsr_copy(sinvb(ispin)%matrix, matrix_ks_nosym)
474 CALL dbcsr_release(matrix_ks_nosym)
475 CALL timestop(handle)
489 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: s_mat
490 TYPE(rt_prop_type),
POINTER :: rtp
492 CHARACTER(len=*),
PARAMETER :: routinen =
's_matrices_create'
493 REAL(kind=
dp),
PARAMETER :: one = 1.0_dp, zero = 0.0_dp
496 TYPE(dbcsr_type),
POINTER :: s_half, s_inv, s_minus_half
498 CALL timeset(routinen, handle)
500 CALL get_rtp(rtp=rtp, s_inv=s_inv)
502 IF (rtp%linear_scaling)
THEN
503 CALL get_rtp(rtp=rtp, s_half=s_half, s_minus_half=s_minus_half)
505 rtp%newton_schulz_order, rtp%lanzcos_threshold, rtp%lanzcos_max_iter)
506 CALL dbcsr_multiply(
"N",
"N", one, s_minus_half, s_minus_half, zero, s_inv, &
507 filter_eps=rtp%filter_eps)
509 CALL dbcsr_copy(s_inv, s_mat(1)%matrix)
511 blacs_env=rtp%ao_ao_fmstruct%context)
513 blacs_env=rtp%ao_ao_fmstruct%context, upper_to_full=.true.)
516 CALL timestop(handle)
527 SUBROUTINE complex_frobenius_norm(frob_norm, mat_re, mat_im)
529 REAL(kind=
dp),
INTENT(out) :: frob_norm
530 TYPE(dbcsr_type),
POINTER :: mat_re, mat_im
532 CHARACTER(len=*),
PARAMETER :: routinen =
'complex_frobenius_norm'
533 REAL(kind=
dp),
PARAMETER :: one = 1.0_dp, zero = 0.0_dp
535 INTEGER :: col_atom, handle, row_atom
537 REAL(
dp),
DIMENSION(:),
POINTER :: block_values, block_values2
538 TYPE(dbcsr_iterator_type) :: iter
539 TYPE(dbcsr_type),
POINTER :: tmp
541 CALL timeset(routinen, handle)
545 CALL dbcsr_create(tmp, template=mat_re)
547 CALL dbcsr_add(tmp, mat_re, zero, one)
548 CALL dbcsr_add(tmp, mat_im, zero, one)
549 CALL dbcsr_set(tmp, zero)
551 CALL dbcsr_iterator_start(iter, tmp)
552 DO WHILE (dbcsr_iterator_blocks_left(iter))
553 CALL dbcsr_iterator_next_block(iter, row_atom, col_atom, block_values)
554 CALL dbcsr_get_block_p(mat_re, row_atom, col_atom, block_values2, found=found)
556 block_values = block_values2*block_values2
558 CALL dbcsr_get_block_p(mat_im, row_atom, col_atom, block_values2, found=found)
560 block_values = block_values + block_values2*block_values2
562 block_values = sqrt(block_values)
564 CALL dbcsr_iterator_stop(iter)
565 frob_norm = dbcsr_frobenius_norm(tmp)
567 CALL dbcsr_deallocate_matrix(tmp)
569 CALL timestop(handle)
571 END SUBROUTINE complex_frobenius_norm
584 SUBROUTINE purify_mcweeny_complex_nonorth(P, s_mat, eps, eps_small, max_iter, threshold)
586 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: p, s_mat
587 REAL(kind=
dp),
INTENT(in) :: eps, eps_small
588 INTEGER,
INTENT(in) :: max_iter
589 REAL(kind=
dp),
INTENT(in) :: threshold
591 CHARACTER(len=*),
PARAMETER :: routinen =
'purify_mcweeny_complex_nonorth'
592 REAL(kind=
dp),
PARAMETER :: one = 1.0_dp, zero = 0.0_dp
594 INTEGER :: handle, i, im,
imax, ispin, re, unit_nr
595 REAL(kind=
dp) :: frob_norm
596 TYPE(cp_logger_type),
POINTER :: logger
597 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: ps, psp, tmp
599 CALL timeset(routinen, handle)
602 IF (logger%para_env%is_source())
THEN
608 NULLIFY (tmp, ps, psp)
613 CALL dbcsr_init_p(ps(i)%matrix)
614 CALL dbcsr_create(ps(i)%matrix, template=p(1)%matrix)
615 CALL dbcsr_init_p(psp(i)%matrix)
616 CALL dbcsr_create(psp(i)%matrix, template=p(1)%matrix)
617 CALL dbcsr_init_p(tmp(i)%matrix)
618 CALL dbcsr_create(tmp(i)%matrix, template=p(1)%matrix)
620 IF (
SIZE(p) == 2)
THEN
621 CALL dbcsr_scale(p(1)%matrix, one/2)
622 CALL dbcsr_scale(p(2)%matrix, one/2)
624 DO ispin = 1,
SIZE(p)/2
627 imax = max(max_iter, 1)
629 CALL dbcsr_multiply(
"N",
"N", one, p(re)%matrix, s_mat(1)%matrix, &
630 zero, ps(re)%matrix, filter_eps=eps_small)
631 CALL dbcsr_multiply(
"N",
"N", one, p(im)%matrix, s_mat(1)%matrix, &
632 zero, ps(im)%matrix, filter_eps=eps_small)
634 p(re)%matrix, p(im)%matrix, zero, psp(re)%matrix, psp(im)%matrix, &
635 filter_eps=eps_small)
636 CALL dbcsr_copy(tmp(re)%matrix, psp(re)%matrix)
637 CALL dbcsr_copy(tmp(im)%matrix, psp(im)%matrix)
638 CALL dbcsr_add(tmp(re)%matrix, p(re)%matrix, 1.0_dp, -1.0_dp)
639 CALL dbcsr_add(tmp(im)%matrix, p(im)%matrix, 1.0_dp, -1.0_dp)
640 CALL complex_frobenius_norm(frob_norm, tmp(re)%matrix, tmp(im)%matrix)
641 IF (unit_nr .GT. 0)
WRITE (unit_nr,
'(t3,a,2f16.8)')
"Deviation from idempotency: ", frob_norm
642 IF (frob_norm .GT. threshold .AND. max_iter > 0)
THEN
643 CALL dbcsr_copy(p(re)%matrix, psp(re)%matrix)
644 CALL dbcsr_copy(p(im)%matrix, psp(im)%matrix)
646 psp(re)%matrix, psp(im)%matrix, 3.0_dp, p(re)%matrix, p(im)%matrix, &
647 filter_eps=eps_small)
648 CALL dbcsr_filter(p(re)%matrix, eps)
649 CALL dbcsr_filter(p(im)%matrix, eps)
651 CALL dbcsr_transposed(tmp(re)%matrix, p(re)%matrix)
652 CALL dbcsr_add(p(re)%matrix, tmp(re)%matrix, one/2, one/2)
653 CALL dbcsr_transposed(tmp(im)%matrix, p(im)%matrix)
654 CALL dbcsr_add(p(im)%matrix, tmp(im)%matrix, one/2, -one/2)
660 CALL dbcsr_transposed(tmp(re)%matrix, p(re)%matrix)
661 CALL dbcsr_add(p(re)%matrix, tmp(re)%matrix, one/2, one/2)
662 CALL dbcsr_transposed(tmp(im)%matrix, p(im)%matrix)
663 CALL dbcsr_add(p(im)%matrix, tmp(im)%matrix, one/2, -one/2)
665 IF (
SIZE(p) == 2)
THEN
666 CALL dbcsr_scale(p(1)%matrix, one*2)
667 CALL dbcsr_scale(p(2)%matrix, one*2)
673 CALL timestop(handle)
675 END SUBROUTINE purify_mcweeny_complex_nonorth
683 SUBROUTINE aspc_extrapolate(rtp, matrix_s, aspc_order)
684 TYPE(rt_prop_type),
POINTER :: rtp
685 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_s
686 INTEGER,
INTENT(in) :: aspc_order
688 CHARACTER(len=*),
PARAMETER :: routinen =
'aspc_extrapolate'
689 COMPLEX(KIND=dp),
PARAMETER :: cone = (1.0_dp, 0.0_dp), &
690 czero = (0.0_dp, 0.0_dp)
691 REAL(kind=
dp),
PARAMETER :: one = 1.0_dp, zero = 0.0_dp
693 INTEGER :: handle, i, iaspc, icol_local, ihist, &
694 imat, k, kdbl, n, naspc, ncol_local, &
696 REAL(kind=
dp) :: alpha
697 TYPE(cp_cfm_type) :: cfm_tmp, cfm_tmp1, csc
698 TYPE(cp_fm_struct_type),
POINTER :: matrix_struct, matrix_struct_new
699 TYPE(cp_fm_type) :: fm_tmp, fm_tmp1, fm_tmp2
700 TYPE(cp_fm_type),
DIMENSION(:),
POINTER :: mos_new
701 TYPE(cp_fm_type),
DIMENSION(:, :),
POINTER :: mo_hist
702 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: rho_new, s_hist
703 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: rho_hist
706 CALL timeset(routinen, handle)
710 IF (rtp%linear_scaling)
THEN
711 CALL get_rtp(rtp=rtp, rho_new=rho_new)
713 CALL get_rtp(rtp=rtp, mos_new=mos_new)
716 naspc = min(rtp%istep, aspc_order)
717 IF (rtp%linear_scaling)
THEN
719 rho_hist => rtp%history%rho_history
722 alpha = (-1.0_dp)**(iaspc + 1)*real(iaspc, kind=
dp)* &
724 ihist = mod(rtp%istep - iaspc, aspc_order) + 1
726 CALL dbcsr_add(rho_new(imat)%matrix, rho_hist(imat, ihist)%matrix, zero, alpha)
728 CALL dbcsr_add(rho_new(imat)%matrix, rho_hist(imat, ihist)%matrix, one, alpha)
733 mo_hist => rtp%history%mo_history
737 alpha = (-1.0_dp)**(iaspc + 1)*real(iaspc, kind=
dp)* &
739 ihist = mod(rtp%istep - iaspc, aspc_order) + 1
748 mo_hist => rtp%history%mo_history
749 s_hist => rtp%history%s_history
750 DO i = 1,
SIZE(mos_new)/2
751 NULLIFY (matrix_struct, matrix_struct_new)
754 mos_new(2*i)%matrix_struct, &
755 mos_new(2*i)%matrix_struct%context, &
769 ncol_local=ncol_local)
772 template_fmstruct=mos_new(2*i)%matrix_struct, &
783 DO icol_local = 1, ncol_local
784 fm_tmp%local_data(:, icol_local) = mos_new(2*i - 1)%local_data(:, icol_local)
785 fm_tmp%local_data(:, icol_local + ncol_local) = mos_new(2*i)%local_data(:, icol_local)
790 DO icol_local = 1, ncol_local
791 cfm_tmp%local_data(:, icol_local) = cmplx(fm_tmp1%local_data(:, icol_local), &
792 fm_tmp1%local_data(:, icol_local + ncol_local),
dp)
793 cfm_tmp1%local_data(:, icol_local) = cmplx(mos_new(2*i - 1)%local_data(:, icol_local), &
794 mos_new(2*i)%local_data(:, icol_local),
dp)
796 CALL parallel_gemm(
'C',
'N', k, k, n, cone, cfm_tmp1, cfm_tmp, czero, csc)
799 DO icol_local = 1, ncol_local
800 mos_new(2*i - 1)%local_data(:, icol_local) = real(cfm_tmp1%local_data(:, icol_local),
dp)
801 mos_new(2*i)%local_data(:, icol_local) = aimag(cfm_tmp1%local_data(:, icol_local))
806 CALL cp_fm_release(fm_tmp)
807 CALL cp_fm_release(fm_tmp1)
808 CALL cp_fm_release(fm_tmp2)
815 CALL timestop(handle)
817 END SUBROUTINE aspc_extrapolate
828 TYPE(rt_prop_type),
POINTER :: rtp
829 TYPE(cp_fm_type),
DIMENSION(:),
POINTER :: mos
830 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: rho
831 TYPE(dbcsr_p_type),
DIMENSION(:),
OPTIONAL, &
837 IF (rtp%linear_scaling)
THEN
839 CALL dbcsr_copy(rtp%history%rho_history(i, ihist)%matrix, rho(i)%matrix)
843 CALL cp_fm_to_fm(mos(i), rtp%history%mo_history(i, ihist))
845 IF (
PRESENT(s_mat))
THEN
846 IF (
ASSOCIATED(rtp%history%s_history(ihist)%matrix))
THEN
848 CALL dbcsr_deallocate_matrix(rtp%history%s_history(ihist)%matrix)
850 ALLOCATE (rtp%history%s_history(ihist)%matrix)
851 CALL dbcsr_copy(rtp%history%s_history(ihist)%matrix, s_mat(1)%matrix)
static int imax(int x, int y)
Returns the larger of two given integer (missing from the C standard)
collects all references to literature in CP2K as new algorithms / method are included from literature...
integer, save, public kuhne2007
integer, save, public kolafa2004
Handles all functions related to the CELL.
Basic linear algebra operations for complex full matrices.
subroutine, public cp_cfm_triangular_multiply(triangular_matrix, matrix_b, side, transa_tr, invert_tr, uplo_tr, unit_diag_tr, n_rows, n_cols, alpha)
Multiplies in place by a triangular matrix: matrix_b = alpha op(triangular_matrix) matrix_b or (if si...
subroutine, public cp_cfm_cholesky_decompose(matrix, n, info_out)
Used to replace a symmetric positive definite matrix M with its Cholesky decomposition U: M = U^T * U...
Represents a complex full matrix distributed on many processors.
subroutine, public cp_cfm_create(matrix, matrix_struct, name)
Creates a new full matrix with the given structure.
subroutine, public cp_cfm_release(matrix)
Releases a full matrix.
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
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_invert(matrix, n, para_env, blacs_env, upper_to_full)
used to replace the cholesky decomposition by the inverse
DBCSR operations in CP2K.
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
basic linear algebra operations for full matrices
subroutine, public cp_fm_scale_and_add(alpha, matrix_a, beta, matrix_b)
calc A <- alpha*A + beta*B optimized for alpha == 1.0 (just add beta*B) and beta == 0....
represent the structure of a full matrix
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
subroutine, public cp_fm_struct_double(fmstruct, struct, context, col, row)
creates a struct with twice the number of blocks on each core. If matrix A has to be multiplied with ...
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
represent a full matrix distributed on many processors
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
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
all routins needed for a nonperiodic electric field
subroutine, public efield_potential_lengh_gauge(qs_env)
Replace the original implementation of the electric-electronic interaction in the length gauge....
Routines useful for iterative matrix calculations.
subroutine, public matrix_sqrt_newton_schulz(matrix_sqrt, matrix_sqrt_inv, matrix, threshold, order, eps_lanczos, max_iter_lanczos, symmetrize, converged)
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
Routines for calculating a complex matrix exponential with dbcsr matrices. Based on the code in matri...
subroutine, public cp_complex_dbcsr_gemm_3(transa, transb, alpha, A_re, A_im, B_re, B_im, beta, C_re, C_im, filter_eps)
Convenience function. Computes the matrix multiplications needed for the multiplication of complex sp...
Collection of simple mathematical functions and subroutines.
elemental real(kind=dp) function, public binomial(n, k)
The binomial coefficient n over k for 0 <= k <= n is calculated, otherwise zero is returned.
basic linear algebra operations for full matrixes
Utility subroutine for qs energy calculation.
subroutine, public qs_energies_init(qs_env, calc_forces)
Refactoring of qs_energies_scf. Driver routine for the initial setup and calculations for a qs energy...
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, 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, rhs)
Get the QUICKSTEP environment.
routines that build the Kohn-Sham matrix (i.e calculate the coulomb and xc parts
subroutine, public qs_ks_update_qs_env(qs_env, calculate_forces, just_energy, print_active)
updates the Kohn Sham matrix of the given qs_env (facility method)
subroutine, public set_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, complex_ks, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, kinetic, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_ks_im_kp, vppl, rho_core, rho_nlcc, rho_nlcc_g, vee, neighbor_list_id, kpoints, sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, task_list, task_list_soft, subsys, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env)
...
Routines for calculating a complex matrix exponential.
subroutine, public propagate_exp(rtp, rtp_control)
performs propagations if explicit matrix exponentials are used ETRS: exp(i*H(t+dt)*dt/2)*exp(i*H(t)*d...
subroutine, public propagate_bch(rtp, rtp_control)
Propagation using the Baker-Campbell-Hausdorff expansion, currently only works for rtp.
subroutine, public propagate_exp_density(rtp, rtp_control)
Propagation of the density matrix instead of the atomic orbitals via a matrix exponential.
subroutine, public propagate_arnoldi(rtp, rtp_control)
computes U_prop*MOs using arnoldi subspace algorithm
Routines for propagating the orbitals.
subroutine, public calc_sinvh(rtp, matrix_ks, matrix_ks_im, rtp_control)
computes S_inv*H, if needed Sinv*B and S_inv*H_imag and store these quantities to the
subroutine, public propagation_step(qs_env, rtp, rtp_control)
performs a single propagation step a(t+Dt)=U(t+Dt,t)*a(0) and calculates the new exponential
subroutine, public put_data_to_history(rtp, mos, rho, s_mat, ihist)
...
subroutine, public s_matrices_create(s_mat, rtp)
calculates the needed overlap-like matrices depending on the way the exponential is calculated,...
Routine for the real time propagation output.
subroutine, public report_density_occupation(filter_eps, rho)
Reports the sparsity pattern of the complex density matrix.
subroutine, public rt_convergence(rtp, matrix_s, delta_mos, delta_eps)
computes the convergence criterion for RTP and EMD
subroutine, public rt_convergence_density(rtp, delta_P, delta_eps)
computes the convergence criterion for RTP and EMD based on the density matrix
Types and set_get for real time propagation depending on runtype and diagonalization method different...
subroutine, public get_rtp(rtp, exp_H_old, exp_H_new, H_last_iter, rho_old, rho_next, rho_new, mos, mos_new, mos_old, mos_next, S_inv, S_half, S_minus_half, B_mat, C_mat, propagator_matrix, mixing, mixing_factor, S_der, dt, nsteps, SinvH, SinvH_imag, SinvB, admm_mos)
...
subroutine, public calc_update_rho_sparse(qs_env)
Copies the density matrix back into the qs_envrhorho_ao.
subroutine, public calc_s_derivs(qs_env)
Calculates dS/dR respectily the velocity weighted derivatves only needed for ehrenfest MD.
subroutine, public calc_update_rho(qs_env)
calculates the density from the complex MOs and passes the density to qs_env.
Routines to perform the RTP in the velocity gauge.
subroutine, public velocity_gauge_ks_matrix(qs_env, subtract_nl_term)
...
subroutine, public update_vector_potential(qs_env, dft_control)
Update the vector potential in the case where a time-dependant electric field is apply.