62 USE iso_c_binding,
ONLY: c_null_char
154#include "./base/base_uses.f90"
159 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'qs_active_space_methods'
164 INTEGER :: unit_nr = -1, bra_start = -1, ket_start = -1
166 PROCEDURE :: func => eri_fcidump_print_func
170 INTEGER :: bra_start = 0, ket_start = 0
171 REAL(KIND=
dp) :: checksum = 0.0_dp
174 PROCEDURE :: func => eri_fcidump_checksum_func
175 END TYPE eri_fcidump_checksum
186 CLASS(eri_fcidump_checksum) :: this
187 INTEGER,
INTENT(IN) :: bra_start, ket_start
188 this%bra_start = bra_start
189 this%ket_start = ket_start
199 CHARACTER(len=*),
PARAMETER :: routinen =
'active_space_main'
201 CHARACTER(len=10) :: cshell, lnam(5)
202 CHARACTER(len=default_path_length) :: qcschema_filename
203 INTEGER :: as_solver, eri_method, eri_operator, eri_print, group_size, handle, i, iatom, &
204 ishell, isp, ispin, iw, j, jm, m, max_orb_ind, mselect, n1, n2, nao, natom, nel, &
205 nelec_active, nelec_inactive, nelec_total, nmo, nmo_active, nmo_available, nmo_inactive, &
206 nmo_inactive_remaining, nmo_occ, nmo_virtual, nn1, nn2, nrow_global, nspins
207 INTEGER,
DIMENSION(5) :: nshell
208 INTEGER,
DIMENSION(:),
POINTER :: invals
209 LOGICAL :: do_kpoints, ex_operator, ex_perd, &
210 explicit, isolated, stop_after_print, &
212 REAL(kind=
dp) :: eri_eps_filter, eri_eps_grid, eri_eps_int, eri_gpw_cutoff, eri_op_param, &
213 eri_rcut, eri_rel_cutoff, fel, focc, maxocc, nze_percentage
214 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: eigenvalues
215 REAL(kind=
dp),
DIMENSION(:),
POINTER :: evals_virtual
221 TYPE(
cp_fm_type),
POINTER :: fm_target_active, fm_target_inactive, &
222 fmat, mo_coeff, mo_ref, mo_target
224 TYPE(dbcsr_csr_type),
POINTER :: eri_mat
225 TYPE(
dbcsr_p_type),
DIMENSION(:),
POINTER :: ks_matrix, rho_ao, s_matrix
229 TYPE(
mo_set_type),
POINTER :: mo_set, mo_set_active, mo_set_inactive
244 IF (.NOT. explicit)
RETURN
245 CALL timeset(routinen, handle)
251 WRITE (iw,
'(/,T2,A)') &
252 '!-----------------------------------------------------------------------------!'
253 WRITE (iw,
'(T26,A)')
"Generate Active Space Hamiltonian"
254 WRITE (iw,
'(T27,A)')
"Interface to CAS-CI and DMRG-CI"
255 WRITE (iw,
'(T2,A)') &
256 '!-----------------------------------------------------------------------------!'
260 CALL get_qs_env(qs_env, do_kpoints=do_kpoints)
262 CALL cp_abort(__location__,
"K-points not allowd for Active Space Interface")
265 NULLIFY (active_space_env)
267 active_space_env%energy_total = 0.0_dp
268 active_space_env%energy_ref = 0.0_dp
269 active_space_env%energy_inactive = 0.0_dp
270 active_space_env%energy_active = 0.0_dp
276 active_space_env%fcidump = .true.
281 active_space_env%qcschema = .true.
282 active_space_env%qcschema_filename = qcschema_filename
288 SELECT CASE (active_space_env%model)
290 WRITE (iw,
'(T3,A)')
"Hartree-Fock model for interaction Hamiltonian"
292 WRITE (iw,
'(T3,A)')
"Range-separated DFT model for interaction Hamiltonian"
294 WRITE (iw,
'(T3,A)')
"DMFT model for interaction Hamiltonian"
296 cpabort(
"Unknown Model")
302 active_space_env%molecule = isolated
304 IF (active_space_env%molecule)
THEN
305 WRITE (iw,
'(T3,A)')
"System is treated without periodicity"
312 CALL get_qs_env(qs_env, nelectron_total=nelec_total)
314 IF (nelec_active <= 0) cpabort(
"Specify a positive number of active electrons.")
315 IF (nelec_active > nelec_total) cpabort(
"More active electrons than total electrons.")
317 nelec_inactive = nelec_total - nelec_active
318 IF (mod(nelec_inactive, 2) /= 0)
THEN
319 cpabort(
"The remaining number of inactive electrons has to be even.")
322 CALL get_qs_env(qs_env, dft_control=dft_control)
323 nspins = dft_control%nspins
325 WRITE (iw,
'(T3,A,T70,I10)')
"Total number of electrons", nelec_total
326 WRITE (iw,
'(T3,A,T70,I10)')
"Number of inactive electrons", nelec_inactive
327 WRITE (iw,
'(T3,A,T70,I10)')
"Number of active electrons", nelec_active
330 active_space_env%nelec_active = nelec_active
331 active_space_env%nelec_inactive = nelec_inactive
332 active_space_env%nelec_total = nelec_total
333 active_space_env%nspins = nspins
334 active_space_env%multiplicity = dft_control%multiplicity
338 IF (.NOT. explicit)
THEN
339 CALL cp_abort(__location__,
"Number of Active Orbitals has to be specified.")
341 active_space_env%nmo_active = nmo_active
343 nmo_inactive = nelec_inactive/2
344 active_space_env%nmo_inactive = nmo_inactive
348 SELECT CASE (mselect)
350 cpabort(
"Unknown orbital selection method")
352 WRITE (iw,
'(/,T3,A)') &
353 "Active space orbitals selected using energy ordered canonical orbitals"
355 WRITE (iw,
'(/,T3,A)') &
356 "Active space orbitals selected using projected Wannier orbitals"
358 WRITE (iw,
'(/,T3,A)') &
359 "Active space orbitals selected using modified atomic orbitals (MAO)"
361 WRITE (iw,
'(/,T3,A)') &
362 "Active space orbitals selected manually"
365 WRITE (iw,
'(T3,A,T70,I10)')
"Number of inactive orbitals", nmo_inactive
366 WRITE (iw,
'(T3,A,T70,I10)')
"Number of active orbitals", nmo_active
373 IF (iatom <= 0 .OR. iatom > natom)
THEN
375 WRITE (iw,
'(/,T3,A,I3)')
"ERROR: SUBSPACE_ATOM number is not valid", iatom
377 cpabort(
"Select a valid SUBSPACE_ATOM")
384 cshell = adjustl(cshell)
388 IF (cshell(n1:n1) ==
" ")
THEN
392 READ (cshell(n1:),
"(I1,A1)") nshell(i), lnam(i)
398 SELECT CASE (mselect)
400 cpabort(
"Unknown orbital selection method")
405 nmo_occ = nmo_inactive + nmo_active
408 ALLOCATE (active_space_env%inactive_orbitals(nmo_inactive, nspins))
410 DO i = 1, nmo_inactive
411 active_space_env%inactive_orbitals(i, ispin) = i
416 ALLOCATE (active_space_env%active_orbitals(nmo_active, nspins))
419 active_space_env%active_orbitals(i, ispin) = nmo_inactive + i
428 IF (nspins > 1) maxocc = 1.0_dp
429 ALLOCATE (active_space_env%mos_active(nspins))
430 ALLOCATE (active_space_env%mos_inactive(nspins))
432 CALL get_mo_set(mos(ispin), mo_coeff=mo_ref, nao=nao)
433 CALL cp_fm_get_info(mo_ref, context=context, para_env=para_env, nrow_global=nrow_global)
435 CALL allocate_mo_set(active_space_env%mos_active(ispin), nao, nmo_occ, 0, 0.0_dp, maxocc, 0.0_dp)
437 nrow_global=nrow_global, ncol_global=nmo_occ)
438 CALL init_mo_set(active_space_env%mos_active(ispin), fm_struct=fm_struct_tmp, name=
"Active Space MO")
440 IF (nspins == 2)
THEN
441 nel = nelec_inactive/2
445 CALL allocate_mo_set(active_space_env%mos_inactive(ispin), nao, nmo_occ, nel, &
446 REAL(nel, kind=
dp), maxocc, 0.0_dp)
448 nrow_global=nrow_global, ncol_global=nmo_occ)
449 CALL init_mo_set(active_space_env%mos_inactive(ispin), fm_struct=fm_struct_tmp, name=
"Inactive Space MO")
454 IF (dft_control%restricted)
THEN
455 cpabort(
"Unclear how we define MOs in the restricted case ... stopping")
457 IF (dft_control%do_admm)
THEN
458 cpabort(
"ADMM currently not possible for canonical orbital options")
461 ALLOCATE (eigenvalues(nmo_occ, nspins))
463 CALL get_qs_env(qs_env, matrix_ks=ks_matrix, matrix_s=s_matrix, scf_control=scf_control)
467 WRITE (iw,
'(/,T3,A)')
"Calculating virtual MOs..."
473 CALL get_mo_set(mos(ispin), mo_coeff=mo_ref, nmo=nmo_available)
476 nmo_virtual = nmo_occ - nmo_available
477 nmo_virtual = max(nmo_virtual, 0)
479 NULLIFY (evals_virtual)
480 ALLOCATE (evals_virtual(nmo_virtual))
483 nrow_global=nrow_global)
486 nrow_global=nrow_global, ncol_global=nmo_virtual)
487 CALL cp_fm_create(mo_virtual, fm_struct_tmp, name=
"virtual")
491 NULLIFY (local_preconditioner)
494 CALL ot_eigensolver(matrix_h=ks_matrix(ispin)%matrix, matrix_s=s_matrix(1)%matrix, &
495 matrix_c_fm=mo_virtual, matrix_orthogonal_space_fm=mo_ref, &
496 eps_gradient=scf_control%eps_lumos, &
498 iter_max=scf_control%max_iter_lumos, &
499 size_ortho_space=nmo_available)
509 evals_arg=eigenvalues(:, ispin), do_rotation=.true.)
512 mo_set => active_space_env%mos_inactive(ispin)
514 DO i = 1,
SIZE(active_space_env%inactive_orbitals, 1)
515 m = active_space_env%inactive_orbitals(i, ispin)
517 mo_set%eigenvalues(m) = eigenvalues(m, ispin)
519 mo_set%occupation_numbers(m) = 1.0
521 mo_set%occupation_numbers(m) = 2.0
526 mo_set => active_space_env%mos_active(ispin)
529 IF (nspins == 2)
THEN
531 nel = (nelec_active + active_space_env%multiplicity - 1)/2
533 nel = (nelec_active - active_space_env%multiplicity + 1)/2
538 mo_set%nelectron = nel
539 mo_set%n_el_f = real(nel, kind=
dp)
541 m = active_space_env%active_orbitals(i, ispin)
542 IF (m > nmo_available)
THEN
543 CALL cp_fm_to_fm(mo_virtual, mo_target, 1, m - nmo_available, m)
544 eigenvalues(m, ispin) = evals_virtual(m - nmo_available)
545 mo_set%occupation_numbers(m) = 0.0
548 mo_set%occupation_numbers(m) = mos(ispin)%occupation_numbers(m)
550 mo_set%eigenvalues(m) = eigenvalues(m, ispin)
553 DEALLOCATE (evals_virtual)
560 WRITE (iw,
'(/,T3,A,I3,T66,A)')
"Canonical Orbital Selection for spin", ispin, &
562 DO i = 1, nmo_inactive, 4
563 jm = min(3, nmo_inactive - i)
564 WRITE (iw,
'(T3,4(F14.6,A5))') (eigenvalues(i + j, ispin),
" [I]", j=0, jm)
566 DO i = nmo_inactive + 1, nmo_inactive + nmo_active, 4
567 jm = min(3, nmo_inactive + nmo_active - i)
568 WRITE (iw,
'(T3,4(F14.6,A5))') (eigenvalues(i + j, ispin),
" [A]", j=0, jm)
570 WRITE (iw,
'(/,T3,A,I3)')
"Active Orbital Indices for spin", ispin
571 DO i = 1,
SIZE(active_space_env%active_orbitals, 1), 4
572 jm = min(3,
SIZE(active_space_env%active_orbitals, 1) - i)
573 WRITE (iw,
'(T3,4(I4))') (active_space_env%active_orbitals(i + j, ispin), j=0, jm)
577 DEALLOCATE (eigenvalues)
582 IF (dft_control%restricted)
THEN
583 cpabort(
"Unclear how we define MOs in the restricted case ... stopping")
585 IF (dft_control%do_admm)
THEN
586 cpabort(
"ADMM currently not possible for manual orbitals selection")
590 IF (.NOT. explicit)
THEN
591 CALL cp_abort(__location__,
"Manual orbital selection requires to explicitly "// &
592 "set the active orbital indices via ACTIVE_ORBITAL_INDICES")
595 IF (nspins == 1)
THEN
596 cpassert(
SIZE(invals) == nmo_active)
598 cpassert(
SIZE(invals) == 2*nmo_active)
600 ALLOCATE (active_space_env%inactive_orbitals(nmo_inactive, nspins))
601 ALLOCATE (active_space_env%active_orbitals(nmo_active, nspins))
605 active_space_env%active_orbitals(i, ispin) = invals(i + (ispin - 1)*nmo_active)
612 max_orb_ind = maxval(invals)
614 IF (nspins > 1) maxocc = 1.0_dp
615 ALLOCATE (active_space_env%mos_active(nspins))
616 ALLOCATE (active_space_env%mos_inactive(nspins))
619 CALL get_mo_set(mos(ispin), mo_coeff=mo_ref, nao=nao)
620 CALL cp_fm_get_info(mo_ref, context=context, para_env=para_env, nrow_global=nrow_global)
621 CALL allocate_mo_set(active_space_env%mos_active(ispin), nao, max_orb_ind, 0, 0.0_dp, maxocc, 0.0_dp)
623 nrow_global=nrow_global, ncol_global=max_orb_ind)
624 CALL init_mo_set(active_space_env%mos_active(ispin), fm_struct=fm_struct_tmp, name=
"Active Space MO")
628 IF (nspins == 2)
THEN
629 nel = nelec_inactive/2
633 CALL allocate_mo_set(active_space_env%mos_inactive(ispin), nao, max_orb_ind, nel, real(nel, kind=
dp), maxocc, 0.0_dp)
635 nrow_global=nrow_global, ncol_global=max_orb_ind)
636 CALL init_mo_set(active_space_env%mos_inactive(ispin), fm_struct=fm_struct_tmp, name=
"Inactive Space MO")
638 active_space_env%mos_inactive(ispin)%occupation_numbers = 0.0_dp
642 ALLOCATE (eigenvalues(max_orb_ind, nspins))
644 CALL get_qs_env(qs_env, matrix_ks=ks_matrix, matrix_s=s_matrix, scf_control=scf_control)
648 WRITE (iw,
'(/,T3,A)')
"Calculating virtual MOs..."
651 CALL get_mo_set(mos(ispin), mo_coeff=mo_ref, nmo=nmo_available)
652 nmo_virtual = max_orb_ind - nmo_available
653 nmo_virtual = max(nmo_virtual, 0)
655 NULLIFY (evals_virtual)
656 ALLOCATE (evals_virtual(nmo_virtual))
659 nrow_global=nrow_global)
662 nrow_global=nrow_global, ncol_global=nmo_virtual)
663 CALL cp_fm_create(mo_virtual, fm_struct_tmp, name=
"virtual")
667 NULLIFY (local_preconditioner)
669 CALL ot_eigensolver(matrix_h=ks_matrix(ispin)%matrix, matrix_s=s_matrix(1)%matrix, &
670 matrix_c_fm=mo_virtual, matrix_orthogonal_space_fm=mo_ref, &
671 eps_gradient=scf_control%eps_lumos, &
673 iter_max=scf_control%max_iter_lumos, &
674 size_ortho_space=nmo_available)
684 evals_arg=eigenvalues(:, ispin), do_rotation=.true.)
686 mo_set_active => active_space_env%mos_active(ispin)
687 CALL get_mo_set(mo_set_active, mo_coeff=fm_target_active)
688 mo_set_inactive => active_space_env%mos_inactive(ispin)
689 CALL get_mo_set(mo_set_inactive, mo_coeff=fm_target_inactive)
692 nmo_inactive_remaining = nmo_inactive
693 DO i = 1, max_orb_ind
695 IF (any(active_space_env%active_orbitals(:, ispin) == i))
THEN
696 IF (i > nmo_available)
THEN
697 CALL cp_fm_to_fm(mo_virtual, fm_target_active, 1, i - nmo_available, i)
698 eigenvalues(i, ispin) = evals_virtual(i - nmo_available)
699 mo_set_active%occupation_numbers(i) = 0.0
701 CALL cp_fm_to_fm(fm_dummy, fm_target_active, 1, i, i)
702 mo_set_active%occupation_numbers(i) = mos(ispin)%occupation_numbers(i)
704 mo_set_active%eigenvalues(i) = eigenvalues(i, ispin)
706 ELSEIF (nmo_inactive_remaining > 0)
THEN
707 CALL cp_fm_to_fm(fm_dummy, fm_target_inactive, 1, i, i)
709 active_space_env%inactive_orbitals(nmo_inactive - nmo_inactive_remaining + 1, ispin) = i
710 mo_set_inactive%eigenvalues(i) = eigenvalues(i, ispin)
711 mo_set_inactive%occupation_numbers(i) = mos(ispin)%occupation_numbers(i)
713 IF (nmo_inactive_remaining == 1)
THEN
714 mo_set_inactive%homo = i
715 mo_set_inactive%lfomo = i + 1
717 nmo_inactive_remaining = nmo_inactive_remaining - 1
724 DEALLOCATE (evals_virtual)
731 WRITE (iw,
'(/,T3,A,I3,T66,A)')
"Orbital Energies and Selection for spin", ispin,
"[atomic units]"
733 DO i = 1, max_orb_ind, 4
734 jm = min(3, max_orb_ind - i)
735 WRITE (iw,
'(T4)', advance=
"no")
737 IF (any(active_space_env%active_orbitals(:, ispin) == i + j))
THEN
738 WRITE (iw,
'(T3,F12.6,A5)', advance=
"no") eigenvalues(i + j, ispin),
" [A]"
739 ELSEIF (any(active_space_env%inactive_orbitals(:, ispin) == i + j))
THEN
740 WRITE (iw,
'(T3,F12.6,A5)', advance=
"no") eigenvalues(i + j, ispin),
" [I]"
742 WRITE (iw,
'(T3,F12.6,A5)', advance=
"no") eigenvalues(i + j, ispin),
" [V]"
747 WRITE (iw,
'(/,T3,A,I3)')
"Active Orbital Indices for spin", ispin
748 DO i = 1,
SIZE(active_space_env%active_orbitals, 1), 4
749 jm = min(3,
SIZE(active_space_env%active_orbitals, 1) - i)
750 WRITE (iw,
'(T3,4(I4))') (active_space_env%active_orbitals(i + j, ispin), j=0, jm)
754 DEALLOCATE (eigenvalues)
758 NULLIFY (loc_section, loc_print)
760 cpassert(
ASSOCIATED(loc_section))
763 cpabort(
"not yet available")
767 cpabort(
"not yet available")
777 CALL print_orbital_cubes(print_orb, qs_env, active_space_env%mos_active)
779 IF (stop_after_print)
THEN
782 WRITE (iw,
'(/,T2,A)') &
783 '!----------------- Early End of Active Space Interface -----------------------!'
786 CALL timestop(handle)
795 cpassert(
ASSOCIATED(rho_ao))
800 mo_set => active_space_env%mos_inactive(ispin)
802 active_space_env%pmat_inactive(ispin)%matrix => denmat
809 active_space_env%eri%method = eri_method
811 active_space_env%eri%operator = eri_operator
813 active_space_env%eri%operator_parameter = eri_op_param
815 active_space_env%eri%cutoff_radius = eri_rcut
818 active_space_env%eri%eps_integral = eri_eps_int
819 IF (active_space_env%molecule)
THEN
822 IF (sum(cell%perd) /= 0)
THEN
823 cpabort(
"Active space option ISOLATED_SYSTEM requires non-periodic setting")
826 IF (sum(invals) /= 0)
THEN
827 cpabort(
"Active space option ISOLATED_SYSTEM requires non-periodic setting")
830 active_space_env%eri%periodicity(1:3) = 0
831 ELSE IF (ex_perd)
THEN
832 IF (
SIZE(invals) == 1)
THEN
833 active_space_env%eri%periodicity(1:3) = invals(1)
835 active_space_env%eri%periodicity(1:3) = invals(1:3)
839 active_space_env%eri%periodicity(1:3) = cell%perd(1:3)
842 WRITE (iw,
'(/,T3,A)')
"Calculation of Electron Repulsion Integrals"
843 SELECT CASE (eri_method)
845 WRITE (iw,
'(T3,A,T50,A)')
"Integration method",
"GPW Fourier transform over MOs"
847 WRITE (iw,
'(T3,A,T44,A)')
"Integration method",
"Half transformed integrals from GPW"
849 cpabort(
"Unknown ERI method")
851 SELECT CASE (eri_operator)
853 WRITE (iw,
'(T3,A,T66,A)')
"ERI operator",
"Coulomb <1/R>"
855 WRITE (iw,
'(T3,A,T59,A)')
"ERI operator",
"Yukawa <EXP(-a*R)/R>"
856 WRITE (iw,
'(T3,A,T66,F14.3)')
"ERI operator parameter", eri_op_param
858 WRITE (iw,
'(T3,A,T53,A)')
"ERI operator",
"Error function <ERF(a*R)/R>"
859 WRITE (iw,
'(T3,A,T66,F14.3)')
"ERI operator parameter", eri_op_param
861 WRITE (iw,
'(T3,A,T45,A)')
"ERI operator",
"Compl. error function <ERFC(a*R)/R>"
862 WRITE (iw,
'(T3,A,T66,F14.3)')
"ERI operator parameter", eri_op_param
864 WRITE (iw,
'(T3,A,T41,A)')
"ERI operator",
"Gaussian attenuated <EXP(-a*R^2)/R>"
865 WRITE (iw,
'(T3,A,T66,F14.3)')
"ERI operator parameter", eri_op_param
867 WRITE (iw,
'(T3,A,T53,A)')
"ERI operator",
"Truncated Coulomb <H(a-R)/R>"
868 WRITE (iw,
'(T3,A,T66,F14.3)')
"ERI operator parameter", eri_op_param
870 cpabort(
"Unknown ERI operator")
872 WRITE (iw,
'(T3,A,T68,E12.4)')
"Accuracy of ERI", eri_eps_int
873 WRITE (iw,
'(T3,A,T71,3I3)')
"Periodicity", active_space_env%eri%periodicity(1:3)
874 IF (product(active_space_env%eri%periodicity(1:3)) == 0)
THEN
875 IF (eri_rcut > 0.0_dp)
WRITE (iw,
'(T3,A,T65,F14.6)')
"Periodicity (Cutoff)", eri_rcut
878 WRITE (iw,
'(T3,A,T68,I12)')
"Total Number of ERI", (nmo_active**4)/8
880 WRITE (iw,
'(T3,A,T68,I12)')
"Total Number of ERI (aa|aa)", (nmo_active**4)/8
881 WRITE (iw,
'(T3,A,T68,I12)')
"Total Number of ERI (bb|bb)", (nmo_active**4)/8
882 WRITE (iw,
'(T3,A,T68,I12)')
"Total Number of ERI (aa|bb)", (nmo_active**4)/4
888 m = (nspins*(nspins + 1))/2
889 ALLOCATE (active_space_env%eri%eri(m))
891 CALL get_mo_set(active_space_env%mos_active(1), nmo=nmo)
892 ALLOCATE (active_space_env%eri%eri(i)%csr_mat)
893 eri_mat => active_space_env%eri%eri(i)%csr_mat
904 nn1 = (n1*(n1 + 1))/2
905 nn2 = (n2*(n2 + 1))/2
906 CALL dbcsr_csr_create(eri_mat, nn1, nn2, 0_int_8, 0, 0, para_env%get_handle())
907 active_space_env%eri%norb = nmo
910 SELECT CASE (eri_method)
913 active_space_env%eri%eri_gpw%eps_grid = eri_eps_grid
915 active_space_env%eri%eri_gpw%eps_filter = eri_eps_filter
917 active_space_env%eri%eri_gpw%cutoff = eri_gpw_cutoff
919 active_space_env%eri%eri_gpw%rel_cutoff = eri_rel_cutoff
921 active_space_env%eri%eri_gpw%print_level = eri_print
923 active_space_env%eri%eri_gpw%store_wfn = store_wfn
925 active_space_env%eri%eri_gpw%group_size = group_size
927 active_space_env%eri%eri_gpw%redo_poisson = .true.
930 WRITE (iw,
'(/,T2,A,T71,F10.1)')
"ERI_GPW| Energy cutoff [Ry]", eri_gpw_cutoff
931 WRITE (iw,
'(T2,A,T71,F10.1)')
"ERI_GPW| Relative energy cutoff [Ry]", eri_rel_cutoff
934 CALL calculate_eri_gpw(active_space_env%mos_active, active_space_env%active_orbitals, active_space_env%eri, qs_env, iw)
937 cpabort(
"Unknown ERI method")
940 DO isp = 1,
SIZE(active_space_env%eri%eri)
941 eri_mat => active_space_env%eri%eri(isp)%csr_mat
942 nze_percentage = 100.0_dp*(real(eri_mat%nze_total, kind=
dp) &
943 /real(eri_mat%nrows_total, kind=
dp))/real(eri_mat%ncols_total, kind=
dp)
944 WRITE (iw,
'(/,T2,A,I2,T30,A,T68,I12)')
"ERI_GPW| Spinmatrix:", isp, &
945 "Number of CSR non-zero elements:", eri_mat%nze_total
946 WRITE (iw,
'(T2,A,I2,T30,A,T68,F12.4)')
"ERI_GPW| Spinmatrix:", isp, &
947 "Percentage CSR non-zero elements:", nze_percentage
948 WRITE (iw,
'(T2,A,I2,T30,A,T68,I12)')
"ERI_GPW| Spinmatrix:", isp, &
949 "nrows_total", eri_mat%nrows_total
950 WRITE (iw,
'(T2,A,I2,T30,A,T68,I12)')
"ERI_GPW| Spinmatrix:", isp, &
951 "ncols_total", eri_mat%ncols_total
952 WRITE (iw,
'(T2,A,I2,T30,A,T68,I12)')
"ERI_GPW| Spinmatrix:", isp, &
953 "nrows_local", eri_mat%nrows_local
958 nspins = active_space_env%nspins
959 ALLOCATE (active_space_env%p_active(nspins))
961 mo_set => active_space_env%mos_active(isp)
962 CALL get_mo_set(mo_set, mo_coeff=mo_coeff, nmo=nmo)
963 CALL create_subspace_matrix(mo_coeff, active_space_env%p_active(isp), nmo)
965 SELECT CASE (mselect)
967 cpabort(
"Unknown orbital selection method")
970 IF (nspins == 2) focc = 1.0_dp
972 fmat => active_space_env%p_active(isp)
974 IF (nspins == 2)
THEN
976 nel = (active_space_env%nelec_active + active_space_env%multiplicity - 1)/2
978 nel = (active_space_env%nelec_active - active_space_env%multiplicity + 1)/2
981 nel = active_space_env%nelec_active
984 m = active_space_env%active_orbitals(i, isp)
985 fel = min(focc, real(nel, kind=
dp))
987 nel = nel - nint(fel)
992 cpabort(
"NOT IMPLEMENTED")
994 cpabort(
"NOT IMPLEMENTED")
998 CALL calculate_operators(active_space_env%mos_active, qs_env, active_space_env)
1001 active_space_env%energy_ref = energy%total
1003 CALL subspace_fock_matrix(active_space_env)
1006 CALL set_qs_env(qs_env, active_space=active_space_env)
1010 SELECT CASE (as_solver)
1013 WRITE (iw,
'(/,T3,A)')
"No active space solver specified, skipping embedding calculation"
1016 CALL rsdft_embedding(qs_env, active_space_env, as_input)
1018 cpabort(
"Unknown active space solver")
1022 IF (active_space_env%fcidump)
CALL fcidump(active_space_env, as_input)
1025 IF (active_space_env%qcschema)
THEN
1032 WRITE (iw,
'(/,T2,A)') &
1033 '!-------------------- End of Active Space Interface --------------------------!'
1036 CALL timestop(handle)
1048 SUBROUTINE calculate_operators(mos, qs_env, active_space_env)
1050 TYPE(
mo_set_type),
DIMENSION(:),
INTENT(IN) :: mos
1054 CHARACTER(len=*),
PARAMETER :: routinen =
'calculate_operators'
1056 INTEGER :: handle, is, nmo, nspins
1058 TYPE(
dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_vxc
1059 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
POINTER :: h_matrix, ks_matrix
1061 CALL timeset(routinen, handle)
1063 nspins = active_space_env%nspins
1067 CALL get_qs_env(qs_env, matrix_ks_kp=ks_matrix)
1068 IF (
SIZE(ks_matrix, 2) > 1)
THEN
1069 cpabort(
"No k-points allowed at this point")
1071 ALLOCATE (active_space_env%ks_sub(nspins))
1073 CALL get_mo_set(mo_set=mos(is), mo_coeff=mo_coeff, nmo=nmo)
1074 CALL subspace_operator(mo_coeff, nmo, ks_matrix(is, 1)%matrix, active_space_env%ks_sub(is))
1080 NULLIFY (matrix_vxc)
1081 CALL get_qs_env(qs_env, matrix_vxc=matrix_vxc)
1082 IF (
ASSOCIATED(matrix_vxc))
THEN
1083 ALLOCATE (active_space_env%vxc_sub(nspins))
1085 CALL get_mo_set(mo_set=mos(is), mo_coeff=mo_coeff, nmo=nmo)
1086 CALL subspace_operator(mo_coeff, nmo, matrix_vxc(is)%matrix, active_space_env%vxc_sub(is))
1094 CALL get_qs_env(qs_env=qs_env, matrix_h_kp=h_matrix)
1095 IF (
SIZE(h_matrix, 2) > 1)
THEN
1096 cpabort(
"No k-points allowed at this point")
1098 ALLOCATE (active_space_env%h_sub(nspins))
1100 CALL get_mo_set(mo_set=mos(is), mo_coeff=mo_coeff, nmo=nmo)
1101 CALL subspace_operator(mo_coeff, nmo, h_matrix(1, 1)%matrix, active_space_env%h_sub(is))
1104 CALL timestop(handle)
1106 END SUBROUTINE calculate_operators
1117 SUBROUTINE subspace_operator(orbitals, nmo, op_matrix, op_sub)
1120 INTEGER,
INTENT(IN) :: nmo
1124 CHARACTER(len=*),
PARAMETER :: routinen =
'subspace_operator'
1126 INTEGER :: handle, ncol, nrow
1129 CALL timeset(routinen, handle)
1131 CALL cp_fm_get_info(matrix=orbitals, ncol_global=ncol, nrow_global=nrow)
1132 cpassert(nmo <= ncol)
1136 CALL cp_fm_create(vectors, orbitals%matrix_struct,
"vectors")
1138 CALL create_subspace_matrix(orbitals, op_sub, nmo)
1142 CALL parallel_gemm(
'T',
'N', nmo, nmo, nrow, 1.0_dp, orbitals, vectors, 0.0_dp, op_sub)
1148 CALL timestop(handle)
1150 END SUBROUTINE subspace_operator
1160 SUBROUTINE create_subspace_matrix(orbitals, op_sub, n)
1164 INTEGER,
INTENT(IN) :: n
1172 para_env=orbitals%matrix_struct%para_env, &
1173 context=orbitals%matrix_struct%context)
1174 CALL cp_fm_create(op_sub, fm_struct, name=
"Subspace operator")
1179 END SUBROUTINE create_subspace_matrix
1191 SUBROUTINE calculate_eri_gpw(mos, orbitals, eri_env, qs_env, iw)
1193 TYPE(
mo_set_type),
DIMENSION(:),
INTENT(IN) :: mos
1194 INTEGER,
DIMENSION(:, :),
POINTER :: orbitals
1197 INTEGER,
INTENT(IN) :: iw
1199 CHARACTER(len=*),
PARAMETER :: routinen =
'calculate_eri_gpw'
1201 INTEGER :: col_local, color, handle, i1, i2, i3, i4, i_multigrid, icount2, intcount, &
1202 irange(2), irange_sub(2), isp, isp1, isp2, ispin, iwa1, iwa12, iwa2, iwb1, iwb12, iwb2, &
1203 iwbs, iwbt, iwfn, n_multigrid, ncol_global, ncol_local, nmm, nmo, nmo1, nmo2, &
1204 nrow_global, nrow_local, nspins, number_of_subgroups, nx, row_local, stored_integrals
1205 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: eri_index
1206 INTEGER,
DIMENSION(:),
POINTER :: col_indices, row_indices
1207 LOGICAL :: print1, print2, &
1208 skip_load_balance_distributed
1209 REAL(kind=
dp) :: dvol, erint, pair_int, &
1210 progression_factor, rc, rsize
1211 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: eri
1216 TYPE(
cp_fm_type),
ALLOCATABLE,
DIMENSION(:) :: fm_matrix_pq_rnu, fm_matrix_pq_rs, &
1218 TYPE(
cp_fm_type),
POINTER :: mo_coeff, mo_coeff1, mo_coeff2
1220 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:) :: matrix_pq_rnu, mo_coeff_as
1225 POINTER :: sab_orb_sub
1233 DIMENSION(:, :),
TARGET :: wfn_a
1236 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
1240 CALL timeset(routinen, handle)
1243 SELECT CASE (eri_env%eri_gpw%print_level)
1264 CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
1265 IF (eri_env%eri_gpw%group_size <= 1) eri_env%eri_gpw%group_size = para_env%num_pe
1266 IF (mod(para_env%num_pe, eri_env%eri_gpw%group_size) /= 0) &
1267 cpabort(
"Group size must be a divisor of the total number of processes!")
1269 IF (eri_env%eri_gpw%group_size == para_env%num_pe)
THEN
1270 para_env_sub => para_env
1271 CALL para_env_sub%retain()
1273 blacs_env_sub => blacs_env
1274 CALL blacs_env_sub%retain()
1276 number_of_subgroups = 1
1279 number_of_subgroups = para_env%num_pe/eri_env%eri_gpw%group_size
1280 color = para_env%mepos/eri_env%eri_gpw%group_size
1281 ALLOCATE (para_env_sub)
1282 CALL para_env_sub%from_split(para_env, color)
1284 NULLIFY (blacs_env_sub)
1290 CALL get_qs_env(qs_env, dft_control=dft_control)
1291 ALLOCATE (qs_control)
1292 qs_control_old => dft_control%qs_control
1293 qs_control = qs_control_old
1294 dft_control%qs_control => qs_control
1295 progression_factor = qs_control%progression_factor
1296 n_multigrid =
SIZE(qs_control%e_cutoff)
1299 ALLOCATE (qs_control%e_cutoff(n_multigrid))
1301 qs_control%cutoff = eri_env%eri_gpw%cutoff*0.5_dp
1302 qs_control%e_cutoff(1) = qs_control%cutoff
1303 DO i_multigrid = 2, n_multigrid
1304 qs_control%e_cutoff(i_multigrid) = qs_control%e_cutoff(i_multigrid - 1) &
1307 qs_control%relative_cutoff = eri_env%eri_gpw%rel_cutoff*0.5_dp
1312 CALL create_mat_munu(mat_munu, qs_env, eri_env%eri_gpw%eps_grid, blacs_env_sub, sab_orb_sub=sab_orb_sub, &
1313 do_alloc_blocks_from_nbl=.true., dbcsr_sym_type=dbcsr_type_symmetric)
1318 NULLIFY (pw_env_sub)
1320 CALL pw_env_rebuild(pw_env_sub, qs_env, external_para_env=para_env_sub)
1321 CALL pw_env_get(pw_env_sub, auxbas_pw_pool=auxbas_pw_pool, &
1322 poisson_env=poisson_env)
1323 IF (eri_env%eri_gpw%redo_poisson)
THEN
1325 IF (sum(eri_env%periodicity) /= 0)
THEN
1330 poisson_env%parameters%periodic = eri_env%periodicity
1332 IF (eri_env%cutoff_radius > 0.0_dp)
THEN
1333 poisson_env%green_fft%radius = eri_env%cutoff_radius
1336 rc = cell%hmat(1, 1)
1338 rc = min(rc, 0.5_dp*cell%hmat(iwa1, iwa1))
1340 poisson_env%green_fft%radius = rc
1343 CALL pw_eri_green_create(poisson_env%green_fft, eri_env)
1347 IF (sum(cell%perd) /= sum(eri_env%periodicity))
THEN
1348 IF (sum(eri_env%periodicity) /= 0)
THEN
1349 WRITE (unit=iw, fmt=
"(/,T2,A,T51,A30)") &
1350 "ERI_GPW| Switching Poisson solver to",
"PERIODIC"
1352 WRITE (unit=iw, fmt=
"(/,T2,A,T51,A30)") &
1353 "ERI_GPW| Switching Poisson solver to",
"ANALYTIC"
1357 SELECT CASE (poisson_env%green_fft%method)
1359 WRITE (unit=iw, fmt=
"(T2,A,T51,A30)") &
1360 "ERI_GPW| Poisson Greens function",
"PERIODIC"
1362 WRITE (unit=iw, fmt=
"(T2,A,T51,A30)") &
1363 "ERI_GPW| Poisson Greens function",
"ANALYTIC"
1365 cpabort(
"Wrong Greens function setup")
1372 NULLIFY (task_list_sub)
1373 skip_load_balance_distributed = dft_control%qs_control%skip_load_balance_distributed
1376 reorder_rs_grid_ranks=.true., soft_valid=.false., &
1377 skip_load_balance_distributed=skip_load_balance_distributed, &
1378 pw_env_external=pw_env_sub, sab_orb_external=sab_orb_sub)
1383 ALLOCATE (mo_coeff_as(nspins), matrix_pq_rnu(nspins), &
1384 fm_matrix_pq_rnu(nspins), fm_mo_coeff_as(nspins), fm_matrix_pq_rs(nspins))
1385 DO ispin = 1, nspins
1387 REAL(kind=
dp),
DIMENSION(:, :),
ALLOCATABLE :: c
1390 CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff)
1394 c(:, minval(orbitals(:, ispin)):maxval(orbitals(:, ispin))), &
1395 mat_munu%matrix, gd_array, eri_env%eri_gpw%eps_filter)
1400 CALL dbcsr_create(matrix_pq_rnu(ispin), template=mo_coeff_as(ispin))
1401 CALL dbcsr_set(matrix_pq_rnu(ispin), 0.0_dp)
1403 CALL dbcsr_get_info(matrix_pq_rnu(ispin), nfullrows_total=nrow_global, nfullcols_total=ncol_global)
1407 nrow_global=nrow_global, ncol_global=ncol_global)
1416 nrow_global=ncol_global, ncol_global=ncol_global)
1424 CALL auxbas_pw_pool%create_pw(wfn_r)
1425 CALL auxbas_pw_pool%create_pw(rho_g)
1426 CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, cell=cell, &
1427 particle_set=particle_set, atomic_kind_set=atomic_kind_set)
1429 IF (eri_env%eri_gpw%store_wfn)
THEN
1433 DO ispin = 1, nspins
1436 rsize = real(
SIZE(wfn_r%array), kind=
dp)*nx
1438 IF (print1 .AND. iw > 0)
THEN
1439 rsize = rsize*8._dp/1000000._dp
1440 WRITE (iw,
"(T4,'ERI_GPW|',' Store active orbitals on real space grid ',T63,F12.3,' MB')") rsize
1442 ALLOCATE (wfn_a(nmo, nspins))
1443 DO ispin = 1, nspins
1444 CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nmo=nmo)
1445 DO i1 = 1,
SIZE(orbitals, 1)
1446 iwfn = orbitals(i1, ispin)
1447 CALL auxbas_pw_pool%create_pw(wfn_a(iwfn, ispin))
1449 qs_kind_set, cell, dft_control, particle_set, pw_env_sub)
1450 IF (print2 .AND. iw > 0)
THEN
1451 WRITE (iw,
"(T4,'ERI_GPW|',' Orbital stored ',I4,' Spin ',i1)") iwfn, ispin
1457 ALLOCATE (wfn1, wfn2)
1458 CALL auxbas_pw_pool%create_pw(wfn1)
1459 CALL auxbas_pw_pool%create_pw(wfn2)
1461 ALLOCATE (wfn3, wfn4)
1462 CALL auxbas_pw_pool%create_pw(wfn3)
1463 CALL auxbas_pw_pool%create_pw(wfn4)
1468 CALL auxbas_pw_pool%create_pw(rho_r)
1469 CALL auxbas_pw_pool%create_pw(pot_g)
1474 dvol = rho_r%pw_grid%dvol
1475 CALL mp_group%set_handle(eri_env%eri(1)%csr_mat%mp_group%get_handle())
1478 stored_integrals = 0
1480 CALL get_mo_set(mo_set=mos(isp1), nmo=nmo1, mo_coeff=mo_coeff1)
1481 nmm = (nmo1*(nmo1 + 1))/2
1487 irange_sub =
get_limit(nmm, number_of_subgroups, color)
1488 DO i1 = 1,
SIZE(orbitals, 1)
1489 iwa1 = orbitals(i1, isp1)
1490 IF (eri_env%eri_gpw%store_wfn)
THEN
1491 wfn1 => wfn_a(iwa1, isp1)
1494 qs_kind_set, cell, dft_control, particle_set, pw_env_sub)
1496 DO i2 = i1,
SIZE(orbitals, 1)
1497 iwa2 = orbitals(i2, isp1)
1500 IF (iwa12 < irange_sub(1) .OR. iwa12 > irange_sub(2)) cycle
1501 IF (iwa12 >= irange(1) .AND. iwa12 <= irange(2))
THEN
1502 iwa12 = iwa12 - irange(1) + 1
1506 IF (eri_env%eri_gpw%store_wfn)
THEN
1507 wfn2 => wfn_a(iwa2, isp1)
1510 qs_kind_set, cell, dft_control, particle_set, pw_env_sub)
1518 WRITE (iw,
"(T4,'ERI_GPW| Integral rho_ab ',T32,2I4,' [',I1,']',T58,G20.14)") &
1519 iwa1, iwa2, isp1, erint
1525 IF (pair_int < eri_env%eps_integral) cycle
1530 DO isp2 = isp1, nspins
1532 nx = (nmo2*(nmo2 + 1))/2
1533 ALLOCATE (eri(nx), eri_index(nx))
1535 CALL integrate_v_rspace(rho_r, hmat=mat_munu, qs_env=qs_env, &
1536 calculate_forces=.false., compute_tau=.false., gapw=.false., &
1537 pw_env_external=pw_env_sub, task_list_external=task_list_sub)
1539 CALL dbcsr_multiply(
"N",
"N", 1.0_dp, mat_munu%matrix, mo_coeff_as(isp2), &
1540 0.0_dp, matrix_pq_rnu(isp2), filter_eps=eri_env%eri_gpw%eps_filter)
1543 CALL cp_fm_get_info(fm_matrix_pq_rnu(isp2), ncol_global=ncol_global, nrow_global=nrow_global)
1545 CALL parallel_gemm(
"T",
"N", ncol_global, ncol_global, nrow_global, 0.5_dp, &
1546 fm_matrix_pq_rnu(isp2), fm_mo_coeff_as(isp2), &
1547 0.0_dp, fm_matrix_pq_rs(isp2))
1548 CALL parallel_gemm(
"T",
"N", ncol_global, ncol_global, nrow_global, 0.5_dp, &
1549 fm_mo_coeff_as(isp2), fm_matrix_pq_rnu(isp2), &
1550 1.0_dp, fm_matrix_pq_rs(isp2))
1552 CALL cp_fm_get_info(fm_matrix_pq_rs(isp2), ncol_local=ncol_local, nrow_local=nrow_local, &
1553 col_indices=col_indices, row_indices=row_indices)
1556 DO col_local = 1, ncol_local
1557 iwb2 = orbitals(col_indices(col_local), isp2)
1558 DO row_local = 1, nrow_local
1559 iwb1 = orbitals(row_indices(row_local), isp2)
1561 IF (iwb1 <= iwb2)
THEN
1563 erint = fm_matrix_pq_rs(isp2)%local_data(row_local, col_local)
1564 IF (abs(erint) > eri_env%eps_integral .AND. (iwa12 <= iwb12 .OR. isp1 /= isp2))
THEN
1565 icount2 = icount2 + 1
1566 eri(icount2) = erint
1567 eri_index(icount2) = iwb12
1572 stored_integrals = stored_integrals + icount2
1574 isp = (isp1 - 1)*isp2 + (isp2 - isp1 + 1)
1575 CALL update_csr_matrix(eri_env%eri(isp)%csr_mat, icount2, eri, eri_index, iwa12)
1577 DEALLOCATE (eri, eri_index)
1580 DO isp2 = isp1, nspins
1581 CALL get_mo_set(mo_set=mos(isp2), nmo=nmo2, mo_coeff=mo_coeff2)
1582 nx = (nmo2*(nmo2 + 1))/2
1583 ALLOCATE (eri(nx), eri_index(nx))
1586 IF (isp1 == isp2) iwbs = i1
1587 isp = (isp1 - 1)*isp2 + (isp2 - isp1 + 1)
1588 DO i3 = iwbs,
SIZE(orbitals, 1)
1589 iwb1 = orbitals(i3, isp2)
1590 IF (eri_env%eri_gpw%store_wfn)
THEN
1591 wfn3 => wfn_a(iwb1, isp2)
1594 qs_kind_set, cell, dft_control, particle_set, pw_env_sub)
1599 IF (isp1 == isp2 .AND. i1 == i3) iwbt = i2
1600 DO i4 = iwbt,
SIZE(orbitals, 1)
1601 iwb2 = orbitals(i4, isp2)
1602 IF (eri_env%eri_gpw%store_wfn)
THEN
1603 wfn4 => wfn_a(iwb2, isp2)
1606 qs_kind_set, cell, dft_control, particle_set, pw_env_sub)
1610 icount2 = icount2 + 1
1611 eri(icount2) = erint
1616 CALL wfn_r%pw_grid%para%group%sum(eri)
1621 IF (isp1 == isp2) iwbs = i1
1622 isp = (isp1 - 1)*isp2 + (isp2 - isp1 + 1)
1623 DO i3 = iwbs,
SIZE(orbitals, 1)
1624 iwb1 = orbitals(i3, isp2)
1626 IF (isp1 == isp2 .AND. i1 == i3) iwbt = i2
1627 DO i4 = iwbt,
SIZE(orbitals, 1)
1628 iwb2 = orbitals(i4, isp2)
1629 intcount = intcount + 1
1630 erint = eri(intcount)
1631 IF (abs(erint) > eri_env%eps_integral)
THEN
1632 icount2 = icount2 + 1
1633 eri(icount2) = erint
1634 eri_index(icount2) = eri_index(intcount)
1638 stored_integrals = stored_integrals + icount2
1640 CALL update_csr_matrix(eri_env%eri(isp)%csr_mat, icount2, eri, eri_index, iwa12)
1642 DEALLOCATE (eri, eri_index)
1645 cpabort(
"Unknown option")
1651 IF (print1 .AND. iw > 0)
THEN
1652 WRITE (iw,
"(T4,'ERI_GPW|',' Number of Integrals stored ',T68,I10)") stored_integrals
1655 IF (eri_env%eri_gpw%store_wfn)
THEN
1656 DO ispin = 1, nspins
1657 DO i1 = 1,
SIZE(orbitals, 1)
1658 iwfn = orbitals(i1, ispin)
1659 CALL wfn_a(iwfn, ispin)%release()
1666 DEALLOCATE (wfn1, wfn2)
1670 DEALLOCATE (wfn3, wfn4)
1673 CALL auxbas_pw_pool%give_back_pw(wfn_r)
1674 CALL auxbas_pw_pool%give_back_pw(rho_g)
1675 CALL auxbas_pw_pool%give_back_pw(rho_r)
1676 CALL auxbas_pw_pool%give_back_pw(pot_g)
1680 DO ispin = 1, nspins
1687 DEALLOCATE (mo_coeff_as, matrix_pq_rnu, &
1688 fm_mo_coeff_as, fm_matrix_pq_rnu, fm_matrix_pq_rs)
1691 DEALLOCATE (mat_munu%matrix)
1697 dft_control%qs_control => qs_control_old
1698 DEALLOCATE (qs_control%e_cutoff)
1699 DEALLOCATE (qs_control)
1701 CALL timestop(handle)
1703 END SUBROUTINE calculate_eri_gpw
1712 SUBROUTINE pw_eri_green_create(green, eri_env)
1714 TYPE(greens_fn_type),
INTENT(INOUT) :: green
1715 TYPE(eri_type) :: eri_env
1718 REAL(kind=dp) :: a, ea, g2, g3d, ga, gg, rg, rlength
1721 associate(gf => green%influence_fn, grid => green%influence_fn%pw_grid)
1722 SELECT CASE (green%method)
1725 SELECT CASE (eri_env%operator)
1726 CASE (eri_operator_coulomb)
1727 DO ig = grid%first_gne0, grid%ngpts_cut_local
1729 gf%array(ig) = fourpi/g2
1731 IF (grid%have_g0) gf%array(1) = 0.0_dp
1732 CASE (eri_operator_yukawa)
1733 a = eri_env%operator_parameter**2
1734 DO ig = grid%first_gne0, grid%ngpts_cut_local
1736 gf%array(ig) = fourpi/(a + g2)
1738 IF (grid%have_g0) gf%array(1) = fourpi/a
1739 CASE (eri_operator_erf)
1740 a = eri_env%operator_parameter**2
1741 DO ig = grid%first_gne0, grid%ngpts_cut_local
1744 gf%array(ig) = fourpi/g2*exp(ga)
1746 IF (grid%have_g0) gf%array(1) = 0.0_dp
1747 CASE (eri_operator_erfc)
1748 a = eri_env%operator_parameter**2
1749 DO ig = grid%first_gne0, grid%ngpts_cut_local
1752 gf%array(ig) = fourpi/g2*(1._dp - exp(ga))
1754 IF (grid%have_g0) gf%array(1) = 0.25_dp*fourpi/a
1755 CASE (eri_operator_trunc)
1756 a = eri_env%operator_parameter
1757 DO ig = grid%first_gne0, grid%ngpts_cut_local
1760 IF (ga >= 0.005_dp)
THEN
1761 gf%array(ig) = fourpi/g2*(1.0_dp - cos(ga))
1763 gf%array(ig) = fourpi/g2*ga**2/2.0_dp*(1.0_dp - ga**2/12.0_dp)
1766 IF (grid%have_g0) gf%array(1) = 0.0_dp
1767 CASE (eri_operator_gaussian)
1775 SELECT CASE (eri_env%operator)
1776 CASE (eri_operator_coulomb)
1777 rlength = green%radius
1778 DO ig = grid%first_gne0, grid%ngpts_cut_local
1782 gf%array(ig) = g3d*(1.0_dp - cos(rlength*gg))
1784 IF (grid%have_g0) gf%array(1) = 0.5_dp*fourpi*rlength*rlength
1785 CASE (eri_operator_yukawa)
1786 rlength = green%radius
1787 a = eri_env%operator_parameter
1788 ea = exp(-a*rlength)
1789 DO ig = grid%first_gne0, grid%ngpts_cut_local
1792 g3d = fourpi/(a*a + g2)
1794 gf%array(ig) = g3d*(1.0_dp - ea*(cos(rg) + a/gg*sin(rg)))
1796 IF (grid%have_g0) gf%array(1) = fourpi/(a*a)*(1.0_dp - ea*(1._dp + a*rlength))
1797 CASE (eri_operator_erf)
1798 rlength = green%radius
1799 a = eri_env%operator_parameter**2
1800 DO ig = grid%first_gne0, grid%ngpts_cut_local
1804 gf%array(ig) = fourpi/g2*exp(ga)*(1.0_dp - cos(rlength*gg))
1806 IF (grid%have_g0) gf%array(1) = 0.5_dp*fourpi*rlength*rlength
1807 CASE (eri_operator_erfc)
1808 rlength = green%radius
1809 a = eri_env%operator_parameter**2
1810 DO ig = grid%first_gne0, grid%ngpts_cut_local
1814 gf%array(ig) = fourpi/g2*(1._dp - exp(ga))*(1.0_dp - cos(rlength*gg))
1816 IF (grid%have_g0) gf%array(1) = 0._dp
1817 CASE (eri_operator_trunc)
1818 a = eri_env%operator_parameter
1819 DO ig = grid%first_gne0, grid%ngpts_cut_local
1822 IF (ga >= 0.005_dp)
THEN
1823 gf%array(ig) = fourpi/g2*(1.0_dp - cos(ga))
1825 gf%array(ig) = fourpi/g2*ga**2/2.0_dp*(1.0_dp - ga**2/12.0_dp)
1828 IF (grid%have_g0) gf%array(1) = 0.0_dp
1829 CASE (eri_operator_gaussian)
1840 END SUBROUTINE pw_eri_green_create
1852 SUBROUTINE update_csr_matrix(csr_mat, nnz, rdat, rind, irow)
1854 TYPE(dbcsr_csr_type),
INTENT(INOUT) :: csr_mat
1855 INTEGER,
INTENT(IN) :: nnz
1856 REAL(kind=dp),
DIMENSION(:),
INTENT(IN) :: rdat
1857 INTEGER,
DIMENSION(:),
INTENT(IN) :: rind
1858 INTEGER,
INTENT(IN) :: irow
1860 INTEGER :: k, nrow, nze, nze_new
1863 nze = csr_mat%nze_local
1866 CALL reallocate(csr_mat%nzval_local%r_dp, 1, nze_new)
1867 csr_mat%nzval_local%r_dp(nze + 1:nze_new) = rdat(1:nnz)
1869 CALL reallocate(csr_mat%colind_local, 1, nze_new)
1870 csr_mat%colind_local(nze + 1:nze_new) = rind(1:nnz)
1872 nrow = csr_mat%nrows_local
1873 CALL reallocate(csr_mat%rowptr_local, 1, irow + 1)
1874 csr_mat%rowptr_local(nrow + 1:irow) = nze + 1
1875 csr_mat%rowptr_local(irow + 1) = nze_new + 1
1877 CALL reallocate(csr_mat%nzerow_local, 1, irow)
1878 DO k = nrow + 1, irow
1879 csr_mat%nzerow_local(k) = csr_mat%rowptr_local(k + 1) - csr_mat%rowptr_local(k)
1881 csr_mat%nrows_local = irow
1882 csr_mat%nze_local = csr_mat%nze_local + nnz
1884 csr_mat%nze_total = csr_mat%nze_total + nnz
1885 csr_mat%has_indices = .true.
1887 END SUBROUTINE update_csr_matrix
1895 SUBROUTINE print_orbital_cubes(input, qs_env, mos)
1896 TYPE(section_vals_type),
POINTER :: input
1897 TYPE(qs_environment_type),
POINTER :: qs_env
1898 TYPE(mo_set_type),
DIMENSION(:),
INTENT(IN) :: mos
1900 CHARACTER(LEN=default_path_length) :: filebody, filename, title
1901 INTEGER :: i, imo, isp, nmo, str(3), unit_nr
1902 INTEGER,
DIMENSION(:),
POINTER :: alist, blist, istride
1903 LOGICAL :: do_mo, explicit_a, explicit_b
1904 TYPE(atomic_kind_type),
DIMENSION(:),
POINTER :: atomic_kind_set
1905 TYPE(cell_type),
POINTER :: cell
1906 TYPE(cp_fm_type),
POINTER :: mo_coeff
1907 TYPE(dft_control_type),
POINTER :: dft_control
1908 TYPE(mp_para_env_type),
POINTER :: para_env
1909 TYPE(particle_list_type),
POINTER :: particles
1910 TYPE(particle_type),
DIMENSION(:),
POINTER :: particle_set
1911 TYPE(pw_c1d_gs_type) :: wf_g
1912 TYPE(pw_env_type),
POINTER :: pw_env
1913 TYPE(pw_pool_type),
POINTER :: auxbas_pw_pool
1914 TYPE(pw_r3d_rs_type) :: wf_r
1915 TYPE(qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
1916 TYPE(qs_subsys_type),
POINTER :: subsys
1917 TYPE(section_vals_type),
POINTER :: dft_section, scf_input
1919 CALL section_vals_val_get(input,
"FILENAME", c_val=filebody)
1920 CALL section_vals_val_get(input,
"STRIDE", i_vals=istride)
1921 IF (
SIZE(istride) == 1)
THEN
1922 str(1:3) = istride(1)
1923 ELSEIF (
SIZE(istride) == 3)
THEN
1924 str(1:3) = istride(1:3)
1926 cpabort(
"STRIDE arguments inconsistent")
1928 CALL section_vals_val_get(input,
"ALIST", i_vals=alist, explicit=explicit_a)
1929 CALL section_vals_val_get(input,
"BLIST", i_vals=blist, explicit=explicit_b)
1931 CALL get_qs_env(qs_env=qs_env, &
1932 dft_control=dft_control, &
1933 para_env=para_env, &
1935 atomic_kind_set=atomic_kind_set, &
1936 qs_kind_set=qs_kind_set, &
1938 particle_set=particle_set, &
1942 CALL qs_subsys_get(subsys, particles=particles)
1944 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
1945 CALL auxbas_pw_pool%create_pw(wf_r)
1946 CALL auxbas_pw_pool%create_pw(wf_g)
1948 dft_section => section_vals_get_subs_vals(scf_input,
"DFT")
1950 DO isp = 1,
SIZE(mos)
1951 CALL get_mo_set(mo_set=mos(isp), mo_coeff=mo_coeff, nmo=nmo)
1953 IF (
SIZE(mos) > 1)
THEN
1956 CALL write_mo_set_to_output_unit(mos(isp), atomic_kind_set, qs_kind_set, particle_set, &
1957 dft_section, 4, 0, final_mos=.true., spin=
"ALPHA")
1959 CALL write_mo_set_to_output_unit(mos(isp), atomic_kind_set, qs_kind_set, particle_set, &
1960 dft_section, 4, 0, final_mos=.true., spin=
"BETA")
1962 cpabort(
"Invalid spin")
1965 CALL write_mo_set_to_output_unit(mos(isp), atomic_kind_set, qs_kind_set, particle_set, &
1966 dft_section, 4, 0, final_mos=.true.)
1970 IF (isp == 1 .AND. explicit_a)
THEN
1971 IF (alist(1) == -1)
THEN
1975 DO i = 1,
SIZE(alist)
1976 IF (imo == alist(i)) do_mo = .true.
1979 ELSE IF (isp == 2 .AND. explicit_b)
THEN
1980 IF (blist(1) == -1)
THEN
1984 DO i = 1,
SIZE(blist)
1985 IF (imo == blist(i)) do_mo = .true.
1991 IF (.NOT. do_mo) cycle
1992 CALL calculate_wavefunction(mo_coeff, imo, wf_r, wf_g, atomic_kind_set, &
1993 qs_kind_set, cell, dft_control, particle_set, pw_env)
1994 IF (para_env%is_source())
THEN
1995 WRITE (filename,
'(A,A1,I4.4,A1,I1.1,A)') trim(filebody),
"_", imo,
"_", isp,
".cube"
1996 CALL open_file(filename, unit_number=unit_nr, file_status=
"UNKNOWN", file_action=
"WRITE")
1997 WRITE (title, *)
"Active Orbital ", imo,
" spin ", isp
2001 CALL cp_pw_to_cube(wf_r, unit_nr, title, particles=particles, stride=istride)
2002 IF (para_env%is_source())
THEN
2003 CALL close_file(unit_nr)
2008 CALL auxbas_pw_pool%give_back_pw(wf_r)
2009 CALL auxbas_pw_pool%give_back_pw(wf_g)
2011 END SUBROUTINE print_orbital_cubes
2020 SUBROUTINE fcidump(active_space_env, as_input)
2022 TYPE(active_space_type),
POINTER :: active_space_env
2023 TYPE(section_vals_type),
POINTER :: as_input
2025 INTEGER :: i, i1, i2, i3, i4, isym, iw, m1, m2, &
2027 REAL(kind=dp) :: checksum, esub
2028 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:, :) :: fmat
2029 TYPE(cp_logger_type),
POINTER :: logger
2030 TYPE(eri_fcidump_checksum) :: eri_checksum
2034 logger => cp_get_default_logger()
2035 iw = cp_print_key_unit_nr(logger, as_input,
"FCIDUMP", &
2036 extension=
".fcidump", file_status=
"REPLACE", file_action=
"WRITE", file_form=
"FORMATTED")
2038 nspins = active_space_env%nspins
2039 norb =
SIZE(active_space_env%active_orbitals, 1)
2040 IF (nspins == 1)
THEN
2041 associate(ms2 => active_space_env%multiplicity, &
2042 nelec => active_space_env%nelec_active)
2045 WRITE (iw,
"(A,A,I4,A,I4,A,I2,A)")
"&FCI",
" NORB=", norb,
",NELEC=", nelec,
",MS2=", ms2,
","
2047 WRITE (iw,
"(A,1000(I1,','))")
" ORBSYM=", (isym, i=1, norb)
2049 WRITE (iw,
"(A,I1,A)")
" ISYM=", isym,
","
2050 WRITE (iw,
"(A)")
" /"
2054 CALL active_space_env%eri%eri_foreach(1, active_space_env%active_orbitals, &
2055 eri_fcidump_print(iw, 1, 1), 1, 1)
2056 CALL eri_checksum%set(1, 1)
2057 CALL active_space_env%eri%eri_foreach(1, active_space_env%active_orbitals, eri_checksum, 1, 1)
2061 nmo = active_space_env%eri%norb
2062 ALLOCATE (fmat(nmo, nmo))
2063 CALL replicate_and_symmetrize_matrix(nmo, active_space_env%fock_sub(1), fmat)
2066 DO m1 = 1,
SIZE(active_space_env%active_orbitals, 1)
2067 i1 = active_space_env%active_orbitals(m1, 1)
2068 DO m2 = m1,
SIZE(active_space_env%active_orbitals, 1)
2069 i2 = active_space_env%active_orbitals(m2, 1)
2070 checksum = checksum + abs(fmat(i1, i2))
2071 WRITE (iw,
"(ES23.16,4I4)") fmat(i1, i2), m1, m2, i3, i4
2077 esub = active_space_env%energy_inactive
2078 i1 = 0; i2 = 0; i3 = 0; i4 = 0
2079 checksum = checksum + abs(esub)
2080 IF (iw > 0)
WRITE (iw,
"(ES23.16,4I4)") esub, i1, i2, i3, i4
2084 associate(ms2 => active_space_env%multiplicity, &
2085 nelec => active_space_env%nelec_active)
2088 WRITE (iw,
"(A,A,I4,A,I4,A,I2,A)")
"&FCI",
" NORB=", norb,
",NELEC=", nelec,
",MS2=", ms2,
","
2090 WRITE (iw,
"(A,1000(I1,','))")
" ORBSYM=", (isym, i=1, norb)
2092 WRITE (iw,
"(A,I1,A)")
" ISYM=", isym,
","
2093 WRITE (iw,
"(A,I1,A)")
" UHF=", 1,
","
2094 WRITE (iw,
"(A)")
" /"
2099 CALL active_space_env%eri%eri_foreach(1, active_space_env%active_orbitals, &
2100 eri_fcidump_print(iw, 1, 1), 1, 1)
2101 CALL eri_checksum%set(1, 1)
2102 CALL active_space_env%eri%eri_foreach(1, active_space_env%active_orbitals, eri_checksum, 1, 1)
2104 CALL active_space_env%eri%eri_foreach(2, active_space_env%active_orbitals, &
2105 eri_fcidump_print(iw, 1, norb + 1), 1, 2)
2106 CALL eri_checksum%set(1, norb + 1)
2107 CALL active_space_env%eri%eri_foreach(2, active_space_env%active_orbitals, eri_checksum, 1, 2)
2109 CALL active_space_env%eri%eri_foreach(3, active_space_env%active_orbitals, &
2110 eri_fcidump_print(iw, norb + 1, norb + 1), 2, 2)
2111 CALL eri_checksum%set(norb + 1, norb + 1)
2112 CALL active_space_env%eri%eri_foreach(3, active_space_env%active_orbitals, eri_checksum, 2, 2)
2115 nmo = active_space_env%eri%norb
2116 ALLOCATE (fmat(nmo, nmo))
2117 CALL replicate_and_symmetrize_matrix(nmo, active_space_env%fock_sub(1), fmat)
2121 i1 = active_space_env%active_orbitals(m1, 1)
2123 i2 = active_space_env%active_orbitals(m2, 1)
2124 checksum = checksum + abs(fmat(i1, i2))
2125 WRITE (iw,
"(ES23.16,4I4)") fmat(i1, i2), m1, m2, i3, i4
2131 ALLOCATE (fmat(nmo, nmo))
2132 CALL replicate_and_symmetrize_matrix(nmo, active_space_env%fock_sub(2), fmat)
2135 DO m1 = 1,
SIZE(active_space_env%active_orbitals, 1)
2136 i1 = active_space_env%active_orbitals(m1, 2)
2137 DO m2 = m1,
SIZE(active_space_env%active_orbitals, 1)
2138 i2 = active_space_env%active_orbitals(m2, 2)
2139 checksum = checksum + abs(fmat(i1, i2))
2140 WRITE (iw,
"(ES23.16,4I4)") fmat(i1, i2), m1 + norb, m2 + norb, i3, i4
2146 esub = active_space_env%energy_inactive
2147 i1 = 0; i2 = 0; i3 = 0; i4 = 0
2148 checksum = checksum + abs(esub)
2149 IF (iw > 0)
WRITE (iw,
"(ES23.16,4I4)") esub, i1, i2, i3, i4
2153 CALL cp_print_key_finished_output(iw, logger, as_input,
"FCIDUMP")
2156 iw = cp_logger_get_default_io_unit(logger)
2157 IF (iw > 0)
WRITE (iw,
'(T4,A,T66,F12.8)')
"FCIDUMP| Checksum:", eri_checksum%checksum + checksum
2160 END SUBROUTINE fcidump
2168 SUBROUTINE replicate_and_symmetrize_matrix(norb, distributed_matrix, replicated_matrix)
2169 INTEGER,
INTENT(IN) :: norb
2170 TYPE(cp_fm_type),
INTENT(IN) :: distributed_matrix
2171 REAL(dp),
DIMENSION(:, :),
INTENT(INOUT) :: replicated_matrix
2176 replicated_matrix(:, :) = 0.0_dp
2179 CALL cp_fm_get_element(distributed_matrix, i1, i2, mval)
2180 replicated_matrix(i1, i2) = mval
2181 replicated_matrix(i2, i1) = mval
2184 END SUBROUTINE replicate_and_symmetrize_matrix
2192 SUBROUTINE subspace_fock_matrix(active_space_env)
2194 TYPE(active_space_type),
POINTER :: active_space_env
2196 INTEGER :: i1, i2, is, norb, nspins
2197 REAL(kind=dp) :: eeri, eref, esub, mval
2198 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:, :) :: ks_a_mat, ks_a_ref, ks_b_mat, ks_b_ref, &
2199 ks_mat, ks_ref, p_a_mat, p_b_mat, p_mat
2200 TYPE(cp_fm_type),
POINTER :: matrix, mo_coef
2201 TYPE(dbcsr_csr_type),
POINTER :: eri, eri_aa, eri_ab, eri_bb
2203 eref = active_space_env%energy_ref
2204 nspins = active_space_env%nspins
2206 IF (nspins == 1)
THEN
2207 CALL get_mo_set(active_space_env%mos_active(1), nmo=norb, mo_coeff=mo_coef)
2212 ALLOCATE (ks_mat(norb, norb), ks_ref(norb, norb), p_mat(norb, norb))
2216 CALL replicate_and_symmetrize_matrix(norb, active_space_env%ks_sub(1), ks_mat)
2217 CALL replicate_and_symmetrize_matrix(norb, active_space_env%p_active(1), p_mat)
2220 eri => active_space_env%eri%eri(1)%csr_mat
2221 CALL build_subspace_fock_matrix(active_space_env%active_orbitals, eri, p_mat, ks_ref, active_space_env%eri%method)
2225 eeri = 0.5_dp*sum(ks_ref*p_mat)
2231 esub = eref - sum(ks_mat(1:norb, 1:norb)*p_mat(1:norb, 1:norb)) + eeri
2234 ks_mat(1:norb, 1:norb) = ks_mat(1:norb, 1:norb) - ks_ref(1:norb, 1:norb)
2237 active_space_env%energy_inactive = esub
2239 CALL cp_fm_release(active_space_env%fock_sub)
2240 ALLOCATE (active_space_env%fock_sub(nspins))
2242 matrix => active_space_env%ks_sub(is)
2243 CALL cp_fm_create(active_space_env%fock_sub(is), matrix%matrix_struct, &
2244 name=
"Active Fock operator")
2246 matrix => active_space_env%fock_sub(1)
2249 mval = ks_mat(i1, i2)
2250 CALL cp_fm_set_element(matrix, i1, i2, mval)
2255 CALL get_mo_set(active_space_env%mos_active(1), nmo=norb)
2260 ALLOCATE (ks_a_mat(norb, norb), ks_b_mat(norb, norb), &
2261 & ks_a_ref(norb, norb), ks_b_ref(norb, norb), &
2262 & p_a_mat(norb, norb), p_b_mat(norb, norb))
2263 ks_a_ref(:, :) = 0.0_dp; ks_b_ref(:, :) = 0.0_dp
2265 CALL replicate_and_symmetrize_matrix(norb, active_space_env%p_active(1), p_a_mat)
2266 CALL replicate_and_symmetrize_matrix(norb, active_space_env%p_active(2), p_b_mat)
2267 CALL replicate_and_symmetrize_matrix(norb, active_space_env%ks_sub(1), ks_a_mat)
2268 CALL replicate_and_symmetrize_matrix(norb, active_space_env%ks_sub(2), ks_b_mat)
2271 eri_aa => active_space_env%eri%eri(1)%csr_mat
2272 eri_ab => active_space_env%eri%eri(2)%csr_mat
2273 eri_bb => active_space_env%eri%eri(3)%csr_mat
2274 CALL build_subspace_spin_fock_matrix(active_space_env%active_orbitals, eri_aa, eri_ab, p_a_mat, p_b_mat, ks_a_ref, &
2275 tr_mixed_eri=.false., eri_method=active_space_env%eri%method)
2276 CALL build_subspace_spin_fock_matrix(active_space_env%active_orbitals, eri_bb, eri_ab, p_b_mat, p_a_mat, ks_b_ref, &
2277 tr_mixed_eri=.true., eri_method=active_space_env%eri%method)
2281 eeri = 0.5_dp*(sum(ks_a_ref*p_a_mat) + sum(ks_b_ref*p_b_mat))
2282 esub = eref - sum(ks_a_mat*p_a_mat) - sum(ks_b_mat*p_b_mat) + eeri
2283 ks_a_mat(:, :) = ks_a_mat(:, :) - ks_a_ref(:, :)
2284 ks_b_mat(:, :) = ks_b_mat(:, :) - ks_b_ref(:, :)
2286 active_space_env%energy_inactive = esub
2288 CALL cp_fm_release(active_space_env%fock_sub)
2289 ALLOCATE (active_space_env%fock_sub(nspins))
2291 matrix => active_space_env%ks_sub(is)
2292 CALL cp_fm_create(active_space_env%fock_sub(is), matrix%matrix_struct, &
2293 name=
"Active Fock operator")
2296 matrix => active_space_env%fock_sub(1)
2299 mval = ks_a_mat(i1, i2)
2300 CALL cp_fm_set_element(matrix, i1, i2, mval)
2303 matrix => active_space_env%fock_sub(2)
2306 mval = ks_b_mat(i1, i2)
2307 CALL cp_fm_set_element(matrix, i1, i2, mval)
2313 END SUBROUTINE subspace_fock_matrix
2323 SUBROUTINE build_subspace_fock_matrix(active_orbitals, eri, p_mat, ks_ref, eri_method)
2324 INTEGER,
DIMENSION(:, :),
INTENT(IN) :: active_orbitals
2325 TYPE(dbcsr_csr_type),
INTENT(IN) :: eri
2326 REAL(dp),
DIMENSION(:, :),
INTENT(IN) :: p_mat
2327 REAL(dp),
DIMENSION(:, :),
INTENT(INOUT) :: ks_ref
2328 INTEGER,
INTENT(IN) :: eri_method
2330 INTEGER :: i1, i12, i12l, i2, i3, i34, i34l, i4, &
2331 irptr, m1, m2, nindex, nmo_total, norb
2332 INTEGER,
DIMENSION(2) :: irange
2334 TYPE(mp_comm_type) :: mp_group
2337 norb =
SIZE(active_orbitals, 1)
2338 nmo_total =
SIZE(p_mat, 1)
2339 nindex = (nmo_total*(nmo_total + 1))/2
2340 CALL mp_group%set_handle(eri%mp_group%get_handle())
2341 IF (eri_method == eri_method_gpw_ht)
THEN
2342 irange = [1, nindex]
2344 irange = get_irange_csr(nindex, mp_group)
2347 i1 = active_orbitals(m1, 1)
2349 i2 = active_orbitals(m2, 1)
2350 i12 = csr_idx_to_combined(i1, i2, nmo_total)
2351 IF (i12 >= irange(1) .AND. i12 <= irange(2))
THEN
2352 i12l = i12 - irange(1) + 1
2353 irptr = eri%rowptr_local(i12l) - 1
2354 DO i34l = 1, eri%nzerow_local(i12l)
2355 i34 = eri%colind_local(irptr + i34l)
2356 CALL csr_idx_from_combined(i34, nmo_total, i3, i4)
2357 erint = eri%nzval_local%r_dp(irptr + i34l)
2359 ks_ref(i1, i2) = ks_ref(i1, i2) + erint*p_mat(i3, i4)
2361 ks_ref(i1, i2) = ks_ref(i1, i2) + erint*p_mat(i3, i4)
2363 IF (i12 /= i34)
THEN
2364 ks_ref(i3, i4) = ks_ref(i3, i4) + erint*p_mat(i1, i2)
2366 ks_ref(i3, i4) = ks_ref(i3, i4) + erint*p_mat(i1, i2)
2370 erint = -0.5_dp*erint
2371 ks_ref(i1, i3) = ks_ref(i1, i3) + erint*p_mat(i2, i4)
2373 ks_ref(i2, i3) = ks_ref(i2, i3) + erint*p_mat(i1, i4)
2376 ks_ref(i1, i4) = ks_ref(i1, i4) + erint*p_mat(i2, i3)
2378 IF (i1 /= i2 .AND. i3 /= i4)
THEN
2379 ks_ref(i2, i4) = ks_ref(i2, i4) + erint*p_mat(i1, i3)
2387 i1 = active_orbitals(m1, 1)
2389 i2 = active_orbitals(m2, 1)
2390 ks_ref(i2, i1) = ks_ref(i1, i2)
2393 CALL mp_group%sum(ks_ref)
2395 END SUBROUTINE build_subspace_fock_matrix
2408 SUBROUTINE build_subspace_spin_fock_matrix(active_orbitals, eri_aa, eri_ab, p_a_mat, p_b_mat, ks_a_ref, tr_mixed_eri, eri_method)
2409 INTEGER,
DIMENSION(:, :),
INTENT(IN) :: active_orbitals
2410 TYPE(dbcsr_csr_type),
INTENT(IN) :: eri_aa, eri_ab
2411 REAL(dp),
DIMENSION(:, :),
INTENT(IN) :: p_a_mat, p_b_mat
2412 REAL(dp),
DIMENSION(:, :),
INTENT(INOUT) :: ks_a_ref
2413 LOGICAL,
INTENT(IN) :: tr_mixed_eri
2414 INTEGER,
INTENT(IN) :: eri_method
2416 INTEGER :: i1, i12, i12l, i2, i3, i34, i34l, i4, &
2417 irptr, m1, m2, nindex, nmo_total, &
2419 INTEGER,
DIMENSION(2) :: irange
2421 TYPE(mp_comm_type) :: mp_group
2423 norb =
SIZE(active_orbitals, 1)
2424 nmo_total =
SIZE(p_a_mat, 1)
2425 nindex = (nmo_total*(nmo_total + 1))/2
2426 CALL mp_group%set_handle(eri_aa%mp_group%get_handle())
2427 IF (eri_method == eri_method_gpw_ht)
THEN
2428 irange = [1, nindex]
2430 irange = get_irange_csr(nindex, mp_group)
2432 IF (tr_mixed_eri)
THEN
2440 i1 = active_orbitals(m1, spin1)
2442 i2 = active_orbitals(m2, spin1)
2443 i12 = csr_idx_to_combined(i1, i2, nmo_total)
2444 IF (i12 >= irange(1) .AND. i12 <= irange(2))
THEN
2445 i12l = i12 - irange(1) + 1
2446 irptr = eri_aa%rowptr_local(i12l) - 1
2447 DO i34l = 1, eri_aa%nzerow_local(i12l)
2448 i34 = eri_aa%colind_local(irptr + i34l)
2449 CALL csr_idx_from_combined(i34, nmo_total, i3, i4)
2450 erint = eri_aa%nzval_local%r_dp(irptr + i34l)
2453 ks_a_ref(i1, i2) = ks_a_ref(i1, i2) + erint*p_a_mat(i3, i4)
2454 IF (i12 /= i34)
THEN
2456 ks_a_ref(i3, i4) = ks_a_ref(i3, i4) + erint*p_a_mat(i1, i2)
2459 erint = -1.0_dp*erint
2461 ks_a_ref(i1, i3) = ks_a_ref(i1, i3) + erint*p_a_mat(i2, i4)
2464 ks_a_ref(i2, i3) = ks_a_ref(i2, i3) + erint*p_a_mat(i1, i4)
2468 ks_a_ref(i1, i4) = ks_a_ref(i1, i4) + erint*p_a_mat(i2, i3)
2470 IF (i1 /= i2 .AND. i3 /= i4)
THEN
2472 ks_a_ref(i2, i4) = ks_a_ref(i2, i4) + erint*p_a_mat(i1, i3)
2480 CALL mp_group%set_handle(eri_ab%mp_group%get_handle())
2481 IF (eri_method == eri_method_gpw_ht)
THEN
2482 irange = [1, nindex]
2484 irange = get_irange_csr(nindex, mp_group)
2487 i1 = active_orbitals(m1, 1)
2489 i2 = active_orbitals(m2, 1)
2490 i12 = csr_idx_to_combined(i1, i2, nmo_total)
2491 IF (i12 >= irange(1) .AND. i12 <= irange(2))
THEN
2492 i12l = i12 - irange(1) + 1
2493 irptr = eri_ab%rowptr_local(i12l) - 1
2494 DO i34l = 1, eri_ab%nzerow_local(i12l)
2495 i34 = eri_ab%colind_local(irptr + i34l)
2496 CALL csr_idx_from_combined(i34, nmo_total, i3, i4)
2497 erint = eri_ab%nzval_local%r_dp(irptr + i34l)
2499 IF (tr_mixed_eri)
THEN
2501 ks_a_ref(i3, i4) = ks_a_ref(i3, i4) + erint*p_b_mat(i1, i2)
2504 ks_a_ref(i1, i2) = ks_a_ref(i1, i2) + erint*p_b_mat(i3, i4)
2512 i1 = active_orbitals(m1, spin1)
2514 i2 = active_orbitals(m2, spin1)
2515 ks_a_ref(i2, i1) = ks_a_ref(i1, i2)
2518 CALL mp_group%set_handle(eri_aa%mp_group%get_handle())
2519 CALL mp_group%sum(ks_a_ref)
2521 END SUBROUTINE build_subspace_spin_fock_matrix
2533 SUBROUTINE create_pro_basis(pro_basis_set, zval, ishell, nshell, lnam)
2534 TYPE(gto_basis_set_type),
POINTER :: pro_basis_set
2535 INTEGER,
INTENT(IN) :: zval, ishell
2536 INTEGER,
DIMENSION(:),
INTENT(IN) :: nshell
2537 CHARACTER(len=*),
DIMENSION(:),
INTENT(IN) :: lnam
2539 CHARACTER(len=6),
DIMENSION(:),
POINTER :: sym
2541 INTEGER,
DIMENSION(4, 7) :: ne
2542 INTEGER,
DIMENSION(:),
POINTER :: lq, nq
2543 REAL(kind=dp),
DIMENSION(:),
POINTER :: zet
2544 TYPE(sto_basis_set_type),
POINTER :: sto_basis_set
2546 cpassert(.NOT.
ASSOCIATED(pro_basis_set))
2547 NULLIFY (sto_basis_set)
2554 ne(l, i) = ptable(zval)%e_conv(l - 1) - 2*nj*(i - l)
2555 ne(l, i) = max(ne(l, i), 0)
2556 ne(l, i) = min(ne(l, i), 2*nj)
2559 ALLOCATE (nq(ishell), lq(ishell), zet(ishell), sym(ishell))
2562 SELECT CASE (lnam(i))
2572 cpabort(
"Wrong l QN")
2575 zet(i) = srules(zval, ne, nq(1), lq(1))
2577 CALL allocate_sto_basis_set(sto_basis_set)
2578 CALL set_sto_basis_set(sto_basis_set, nshell=1, nq=nq, lq=lq, zet=zet, symbol=sym)
2579 CALL create_gto_from_sto_basis(sto_basis_set, pro_basis_set, 6)
2580 pro_basis_set%norm_type = 2
2581 CALL init_orb_basis_set(pro_basis_set)
2582 CALL deallocate_sto_basis_set(sto_basis_set)
2584 END SUBROUTINE create_pro_basis
2591 SUBROUTINE update_density_ao(active_space_env, rho_ao)
2592 TYPE(active_space_type),
POINTER :: active_space_env
2593 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: rho_ao
2595 INTEGER :: ispin, nao, nmo, nspins
2596 TYPE(cp_fm_type) :: r, u
2597 TYPE(cp_fm_type),
POINTER :: c_active, p_active_mo
2598 TYPE(dbcsr_type),
POINTER :: p_inactive_ao
2599 TYPE(mo_set_type),
DIMENSION(:),
POINTER :: mos_active
2603 nspins = active_space_env%nspins
2604 mos_active => active_space_env%mos_active
2605 DO ispin = 1, nspins
2607 p_inactive_ao => active_space_env%pmat_inactive(ispin)%matrix
2610 CALL dbcsr_copy(rho_ao(ispin)%matrix, p_inactive_ao)
2613 p_active_mo => active_space_env%p_active(ispin)
2616 CALL cp_fm_create(r, p_active_mo%matrix_struct)
2617 CALL cp_fm_to_fm(p_active_mo, r)
2620 CALL get_mo_set(mos_active(ispin), mo_coeff=c_active, nao=nao, nmo=nmo)
2621 CALL cp_fm_create(u, c_active%matrix_struct)
2622 CALL parallel_gemm(
"N",
"N", nao, nmo, nmo, 1.0_dp, c_active, r, 0.0_dp, u)
2624 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_ao(ispin)%matrix, &
2625 matrix_v=u, matrix_g=c_active, ncol=nmo, alpha=1.0_dp)
2627 CALL cp_fm_release(r)
2628 CALL cp_fm_release(u)
2631 END SUBROUTINE update_density_ao
2643 LOGICAL FUNCTION eri_fcidump_print_func(this, i, j, k, l, val)
RESULT(cont)
2644 CLASS(eri_fcidump_print),
INTENT(inout) :: this
2645 INTEGER,
INTENT(in) :: i, j, k, l
2646 REAL(kind=dp),
INTENT(in) :: val
2649 IF (this%unit_nr > 0)
THEN
2650 WRITE (this%unit_nr,
"(ES23.16,4I4)") val, i + this%bra_start - 1, j + this%bra_start - 1, &
2651 & k + this%ket_start - 1, l + this%ket_start - 1
2655 END FUNCTION eri_fcidump_print_func
2667 LOGICAL FUNCTION eri_fcidump_checksum_func(this, i, j, k, l, val)
RESULT(cont)
2668 CLASS(eri_fcidump_checksum),
INTENT(inout) :: this
2669 INTEGER,
INTENT(in) :: i, j, k, l
2670 REAL(kind=dp),
INTENT(in) :: val
2676 this%checksum = this%checksum + abs(val)
2679 END FUNCTION eri_fcidump_checksum_func
2688 SUBROUTINE update_active_density(p_act_mo, active_space_env, ispin)
2689 REAL(kind=dp),
DIMENSION(:) :: p_act_mo
2690 TYPE(active_space_type),
POINTER :: active_space_env
2693 INTEGER :: i1, i2, m1, m2, nmo_active
2694 REAL(kind=dp) :: mval
2695 TYPE(cp_fm_type),
POINTER :: p_active
2697 p_active => active_space_env%p_active(ispin)
2698 nmo_active = active_space_env%nmo_active
2700 DO i1 = 1, nmo_active
2701 m1 = active_space_env%active_orbitals(i1, ispin)
2702 DO i2 = 1, nmo_active
2703 m2 = active_space_env%active_orbitals(i2, ispin)
2704 mval = p_act_mo(i2 + (i1 - 1)*nmo_active)
2705 CALL cp_fm_set_element(p_active, m1, m2, mval)
2709 END SUBROUTINE update_active_density
2717 SUBROUTINE rsdft_embedding(qs_env, active_space_env, as_input)
2718 TYPE(qs_environment_type),
POINTER :: qs_env
2719 TYPE(active_space_type),
POINTER :: active_space_env
2720 TYPE(section_vals_type),
POINTER :: as_input
2722 CHARACTER(len=*),
PARAMETER :: routinen =
'rsdft_embedding'
2726 CALL timeset(routinen, handle)
2727 cpabort(
"CP2K was compiled with the __NO_SOCKETS option!")
2729 mark_used(active_space_env)
2733 INTEGER :: iw, client_fd, socket_fd, iter, max_iter
2734 LOGICAL :: converged, do_scf_embedding, ionode
2735 REAL(kind=dp) :: alpha, delta_e, energy_corr, energy_new, &
2736 energy_old, energy_scf, eps_iter, t1, &
2738 TYPE(cp_logger_type),
POINTER :: logger
2739 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: rho_ao
2740 TYPE(mo_set_type),
DIMENSION(:),
POINTER :: mos_active
2741 TYPE(mp_para_env_type),
POINTER :: para_env
2742 TYPE(qs_energy_type),
POINTER :: energy
2743 TYPE(qs_rho_type),
POINTER :: rho
2745 CALL timeset(routinen, handle)
2749 logger => cp_get_default_logger()
2750 iw = cp_logger_get_default_io_unit(logger)
2752 CALL get_qs_env(qs_env, para_env=para_env)
2753 ionode = para_env%is_source()
2756 CALL section_vals_val_get(as_input,
"SCF_EMBEDDING", l_val=do_scf_embedding)
2757 active_space_env%do_scf_embedding = do_scf_embedding
2758 CALL section_vals_val_get(as_input,
"MAX_ITER", i_val=max_iter)
2759 CALL section_vals_val_get(as_input,
"EPS_ITER", r_val=eps_iter)
2763 CALL initialize_socket(socket_fd, client_fd, as_input, ionode)
2764 CALL para_env%sync()
2767 CALL send_eri_to_client(client_fd, active_space_env, para_env)
2770 CALL get_qs_env(qs_env, rho=rho, energy=energy)
2771 CALL qs_rho_get(rho, rho_ao=rho_ao)
2774 WRITE (unit=iw, fmt=
"(/,T3,A,T11,A,T21,A,T34,A,T55,A,T75,A,/,T3,A)") &
2775 "Iter",
"Update",
"Time",
"Corr. energy",
"Total energy",
"Change", repeat(
"-", 78)
2782 energy_scf = active_space_env%energy_ref
2783 energy_new = energy_scf
2784 mos_active => active_space_env%mos_active
2788 DO WHILE (iter < max_iter)
2793 CALL send_fock_to_client(client_fd, active_space_env, para_env)
2796 energy_old = energy_new
2797 energy_new = active_space_env%energy_total
2798 energy_corr = energy_new - energy_scf
2799 delta_e = energy_new - energy_old
2806 fmt=
"(T3,I4,T11,A,T21,F4.2,T28,F18.10,T49,F18.10,T70,ES11.2)") &
2807 iter,
'P_Mix', t2 - t1, energy_corr, energy_new, delta_e
2811 CALL update_density_ao(active_space_env, rho_ao)
2814 qs_env%requires_matrix_vxc = .true.
2815 CALL qs_rho_update_rho(rho, qs_env)
2816 CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.true.)
2817 CALL qs_ks_update_qs_env(qs_env)
2820 active_space_env%energy_ref = energy%total
2823 CALL calculate_operators(mos_active, qs_env, active_space_env)
2826 CALL subspace_fock_matrix(active_space_env)
2829 IF (.NOT. active_space_env%do_scf_embedding)
THEN
2831 WRITE (unit=iw, fmt=
"(/,T3,A,I5,A/)") &
2832 "*** one-shot embedding correction finished ***"
2837 ELSEIF (abs(delta_e) <= eps_iter)
THEN
2839 WRITE (unit=iw, fmt=
"(/,T3,A,I5,A/)") &
2840 "*** rsDFT embedding run converged in ", iter,
" iteration(s) ***"
2849 IF (.NOT. converged)
THEN
2851 WRITE (unit=iw, fmt=
"(/,T3,A,I5,A/)") &
2852 "*** rsDFT embedding did not converged after ", iter,
" iteration(s) ***"
2857 energy%total = active_space_env%energy_total
2859 CALL finalize_socket(socket_fd, client_fd, as_input, ionode)
2860 CALL para_env%sync()
2863 CALL timestop(handle)
2865 END SUBROUTINE rsdft_embedding
2875 SUBROUTINE initialize_socket(socket_fd, client_fd, as_input, ionode)
2876 INTEGER,
INTENT(OUT) :: socket_fd, client_fd
2877 TYPE(section_vals_type),
INTENT(IN),
POINTER :: as_input
2878 LOGICAL,
INTENT(IN) :: ionode
2880 CHARACTER(len=*),
PARAMETER :: routinen =
'initialize_socket'
2881 INTEGER,
PARAMETER :: backlog = 10
2883 CHARACTER(len=default_path_length) :: hostname
2884 INTEGER :: handle, iw, port, protocol
2886 TYPE(cp_logger_type),
POINTER :: logger
2888 CALL timeset(routinen, handle)
2890 logger => cp_get_default_logger()
2891 iw = cp_logger_get_default_io_unit(logger)
2894 CALL section_vals_val_get(as_input,
"SOCKET%INET", l_val=inet)
2900 CALL section_vals_val_get(as_input,
"SOCKET%HOST", c_val=hostname)
2901 CALL section_vals_val_get(as_input,
"SOCKET%PORT", i_val=port)
2904 CALL open_bind_socket(socket_fd, protocol, port, trim(hostname)//c_null_char)
2905 WRITE (iw,
'(/,T3,A,A)')
"@SERVER: Created socket with address ", trim(hostname)
2909 WRITE (iw,
'(T3,A)')
"@SERVER: Waiting for requests..."
2911 WRITE (iw,
'(T3,A,I2)')
"@SERVER: Accepted socket with fd ", client_fd
2914 CALL timestop(handle)
2916 END SUBROUTINE initialize_socket
2925 SUBROUTINE finalize_socket(socket_fd, client_fd, as_input, ionode)
2926 INTEGER,
INTENT(IN) :: socket_fd, client_fd
2927 TYPE(section_vals_type),
INTENT(IN),
POINTER :: as_input
2928 LOGICAL,
INTENT(IN) :: ionode
2930 CHARACTER(len=*),
PARAMETER :: routinen =
'finalize_socket'
2931 INTEGER,
PARAMETER :: header_len = 12
2933 CHARACTER(len=default_path_length) :: hostname
2936 CALL timeset(routinen, handle)
2938 CALL section_vals_val_get(as_input,
"SOCKET%HOST", c_val=hostname)
2948 IF (file_exists(trim(hostname)))
THEN
2953 CALL timestop(handle)
2955 END SUBROUTINE finalize_socket
2963 SUBROUTINE send_eri_to_client(client_fd, active_space_env, para_env)
2964 INTEGER,
INTENT(IN) :: client_fd
2965 TYPE(active_space_type),
INTENT(IN),
POINTER :: active_space_env
2966 TYPE(mp_para_env_type),
INTENT(IN),
POINTER :: para_env
2968 CHARACTER(len=*),
PARAMETER :: routinen =
'send_eri_to_client'
2969 INTEGER,
PARAMETER :: header_len = 12
2971 CHARACTER(len=default_string_length) ::
header
2972 INTEGER :: handle, iw
2974 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:) :: eri_aa, eri_ab, eri_bb
2975 TYPE(cp_logger_type),
POINTER :: logger
2977 CALL timeset(routinen, handle)
2979 logger => cp_get_default_logger()
2980 iw = cp_logger_get_default_io_unit(logger)
2981 ionode = para_env%is_source()
2984 ALLOCATE (eri_aa(active_space_env%nmo_active**4))
2985 CALL eri_to_array(active_space_env%eri, eri_aa, active_space_env%active_orbitals, 1, 1)
2986 IF (active_space_env%nspins == 2)
THEN
2987 ALLOCATE (eri_ab(active_space_env%nmo_active**4))
2988 CALL eri_to_array(active_space_env%eri, eri_ab, active_space_env%active_orbitals, 1, 2)
2989 ALLOCATE (eri_bb(active_space_env%nmo_active**4))
2990 CALL eri_to_array(active_space_env%eri, eri_bb, active_space_env%active_orbitals, 2, 2)
2994 IF (ionode)
CALL writebuffer(client_fd,
"STATUS ", header_len)
2997 CALL para_env%sync()
3002 CALL para_env%bcast(
header, para_env%source)
3006 IF (trim(
header) ==
"READY")
THEN
3008 CALL para_env%sync()
3010 CALL writebuffer(client_fd,
"TWOBODY ", header_len)
3011 CALL writebuffer(client_fd, active_space_env%nspins)
3012 CALL writebuffer(client_fd, active_space_env%nmo_active)
3013 CALL writebuffer(client_fd, active_space_env%nelec_active)
3014 CALL writebuffer(client_fd, active_space_env%multiplicity)
3018 IF (active_space_env%nspins == 2)
THEN
3023 ELSE IF (trim(
header) ==
"RECEIVED")
THEN
3029 IF (active_space_env%nspins == 2)
THEN
3034 CALL para_env%sync()
3036 CALL timestop(handle)
3038 END SUBROUTINE send_eri_to_client
3046 SUBROUTINE send_fock_to_client(client_fd, active_space_env, para_env)
3047 INTEGER,
INTENT(IN) :: client_fd
3048 TYPE(active_space_type),
INTENT(IN),
POINTER :: active_space_env
3049 TYPE(mp_para_env_type),
INTENT(IN),
POINTER :: para_env
3051 CHARACTER(len=*),
PARAMETER :: routinen =
'send_fock_to_client'
3052 INTEGER,
PARAMETER :: header_len = 12
3054 CHARACTER(len=default_string_length) ::
header
3055 INTEGER :: handle, iw
3056 LOGICAL :: debug, ionode
3057 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:) :: fock_a, fock_b, p_act_mo_a, p_act_mo_b
3058 TYPE(cp_logger_type),
POINTER :: logger
3060 CALL timeset(routinen, handle)
3065 logger => cp_get_default_logger()
3066 iw = cp_logger_get_default_io_unit(logger)
3067 ionode = para_env%is_source()
3069 ALLOCATE (p_act_mo_a(active_space_env%nmo_active**2))
3070 ALLOCATE (fock_a(active_space_env%nmo_active**2))
3071 IF (active_space_env%nspins == 2)
THEN
3072 ALLOCATE (p_act_mo_b(active_space_env%nmo_active**2))
3073 ALLOCATE (fock_b(active_space_env%nmo_active**2))
3077 associate(act_indices => active_space_env%active_orbitals(:, 1))
3078 CALL subspace_matrix_to_array(active_space_env%fock_sub(1), fock_a, act_indices, act_indices)
3081 IF (active_space_env%nspins == 2)
THEN
3082 associate(act_indices => active_space_env%active_orbitals(:, 2))
3083 CALL subspace_matrix_to_array(active_space_env%fock_sub(2), fock_b, act_indices, act_indices)
3088 IF (ionode)
CALL writebuffer(client_fd,
"STATUS ", header_len)
3092 CALL para_env%sync()
3094 IF (debug .AND. iw > 0)
WRITE (iw, *)
"@SERVER: Waiting for messages..."
3097 CALL para_env%bcast(
header, para_env%source)
3099 IF (debug .AND. iw > 0)
WRITE (iw, *)
"@SERVER: Message from client: ", trim(
header)
3101 IF (trim(
header) ==
"READY")
THEN
3103 CALL para_env%sync()
3105 CALL writebuffer(client_fd,
"ONEBODY ", header_len)
3106 CALL writebuffer(client_fd, active_space_env%energy_inactive)
3110 IF (active_space_env%nspins == 2)
THEN
3115 ELSE IF (trim(
header) ==
"HAVEDATA")
THEN
3117 CALL para_env%sync()
3119 IF (debug .AND. iw > 0)
WRITE (iw, *)
"@SERVER: Qiskit has data to transfer"
3120 CALL writebuffer(client_fd,
"GETDENSITY ", header_len)
3123 CALL readbuffer(client_fd, active_space_env%energy_active)
3124 CALL readbuffer(client_fd, p_act_mo_a,
SIZE(p_act_mo_a))
3125 IF (active_space_env%nspins == 2)
THEN
3126 CALL readbuffer(client_fd, p_act_mo_b,
SIZE(p_act_mo_b))
3131 CALL para_env%bcast(active_space_env%energy_active, para_env%source)
3132 CALL para_env%bcast(p_act_mo_a, para_env%source)
3133 IF (active_space_env%nspins == 2)
THEN
3134 CALL para_env%bcast(p_act_mo_b, para_env%source)
3138 active_space_env%energy_total = active_space_env%energy_inactive + active_space_env%energy_active
3141 CALL update_active_density(p_act_mo_a, active_space_env, 1)
3142 IF (active_space_env%nspins == 2)
THEN
3143 CALL update_active_density(p_act_mo_b, active_space_env, 2)
3152 DEALLOCATE (p_act_mo_a)
3154 IF (active_space_env%nspins == 2)
THEN
3155 DEALLOCATE (p_act_mo_b)
3159 CALL para_env%sync()
3161 CALL timestop(handle)
3163 END SUBROUTINE send_fock_to_client
Define the atomic kind types and their sub types.
pure real(dp) function, public srules(z, ne, n, l)
...
subroutine, public deallocate_sto_basis_set(sto_basis_set)
...
subroutine, public allocate_sto_basis_set(sto_basis_set)
...
subroutine, public create_gto_from_sto_basis(sto_basis_set, gto_basis_set, ngauss, ortho)
...
subroutine, public set_sto_basis_set(sto_basis_set, name, nshell, symbol, nq, lq, zet)
...
subroutine, public init_orb_basis_set(gto_basis_set)
Initialise a Gaussian-type orbital (GTO) basis set data set.
Handles all functions related to the CELL.
methods related to the blacs parallel environment
integer, parameter, public blacs_grid_square
subroutine, public cp_blacs_env_release(blacs_env)
releases the given blacs_env
subroutine, public cp_blacs_env_create(blacs_env, para_env, blacs_grid_layout, blacs_repeatable, row_major, grid_2d)
allocates and initializes a type that represent a blacs context
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
subroutine, public dbcsr_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset)
...
logical function, public dbcsr_iterator_blocks_left(iterator)
...
subroutine, public dbcsr_iterator_stop(iterator)
...
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
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)
...
real(kind=dp) function, public dbcsr_get_occupation(matrix)
...
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
subroutine, public dbcsr_set(matrix, alpha)
...
subroutine, public dbcsr_release(matrix)
...
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
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
subroutine, public cp_dbcsr_plus_fm_fm_t(sparse_matrix, matrix_v, matrix_g, ncol, alpha, keep_sparsity, symmetry_mode)
performs the multiplication sparse_matrix+dense_mat*dens_mat^T if matrix_g is not explicitly given,...
subroutine, public cp_dbcsr_m_by_n_from_template(matrix, template, m, n, sym)
Utility function to create an arbitrary shaped dbcsr matrix with the same processor grid as the templ...
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.
logical function, public file_exists(file_name)
Checks if file exists, considering also the file discovery mechanism.
represent the structure of a full matrix
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
represent a full matrix distributed on many processors
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
subroutine, public cp_fm_get_element(matrix, irow_global, icol_global, alpha, local)
returns an element of a fm this value is valid on every cpu using this call is expensive
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
subroutine, public cp_fm_set_element(matrix, irow_global, icol_global, alpha)
sets an element of a matrix
subroutine, public cp_fm_init_random(matrix, ncol, start_col)
fills a matrix with random numbers
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 ...
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
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, parameter, public debug_print_level
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)
...
integer, parameter, public low_print_level
integer, parameter, public medium_print_level
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
integer, parameter, public cp_p_file
integer, parameter, public high_print_level
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
integer, parameter, public silent_print_level
A wrapper around pw_to_cube() which accepts particle_list_type.
subroutine, public cp_pw_to_cube(pw, unit_nr, title, particles, stride, zero_tails, silent, mpi_io)
...
Types to describe group distributions.
Defines the basic variable types.
integer, parameter, public int_8
integer, parameter, public dp
integer, parameter, public default_string_length
integer, parameter, public default_path_length
Machine interface based on Fortran 2003 and POSIX.
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
Definition of mathematical constants and functions.
real(kind=dp), parameter, public fourpi
Utility routines for the memory handling.
Interface to the message passing library MPI.
subroutine, public mp_para_env_release(para_env)
releases the para object (to be called when you don't want anymore the shared copy of this object)
Calls routines to get RI integrals and calculate total energies.
subroutine, public grep_rows_in_subgroups(para_env, para_env_sub, mo_coeff, gd_array, c)
...
subroutine, public build_dbcsr_from_rows(para_env_sub, mo_coeff_to_build, cread, mat_munu, gd_array, eps_filter)
Encapsulate the building of dbcsr_matrices mo_coeff_(v,o,all)
subroutine, public create_mat_munu(mat_munu, qs_env, eps_grid, blacs_env_sub, do_ri_aux_basis, do_mixed_basis, group_size_prim, do_alloc_blocks_from_nbl, do_kpoints, sab_orb_sub, dbcsr_sym_type)
Encapsulate the building of dbcsr_matrix mat_munu.
basic linear algebra operations for full matrixes
represent a simple array based list of the given type
Define the data structure for the particle information.
Periodic Table related data definitions.
type(atom), dimension(0:nelem), public ptable
computes preconditioners, and implements methods to apply them currently used in qs_ot
methods of pw_env that have dependence on qs_env
subroutine, public pw_env_rebuild(pw_env, qs_env, external_para_env)
rebuilds the pw_env data (necessary if cell or cutoffs change)
subroutine, public pw_env_create(pw_env)
creates a pw_env, if qs_env is given calls pw_env_rebuild
container for various plainwaves related things
subroutine, public pw_env_release(pw_env, para_env)
releases the given pw_env (see doc/ReferenceCounting.html)
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
functions related to the poisson solver on regular grids
integer, parameter, public pw_poisson_periodic
integer, parameter, public periodic3d
integer, parameter, public analytic0d
integer, parameter, public pw_poisson_analytic
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
The module to read/write QCSchema HDF5 files for interfacing CP2K with other programs.
subroutine, public qcschema_env_release(qcschema_env)
Releases the allocated memory of a qcschema environment.
subroutine, public qcschema_env_create(qcschema_env, qs_env)
Create and initialize a qcschema object from a quickstep environment.
subroutine, public qcschema_to_hdf5(qcschema_env, filename)
Writes a qcschema object to an hdf5 file.
Determine active space Hamiltonian.
subroutine eri_fcidump_set(this, bra_start, ket_start)
Sets the starting indices of the bra and ket.
subroutine, public active_space_main(qs_env)
Main method for determining the active space Hamiltonian.
The types needed for the calculation of active space Hamiltonians.
subroutine, public csr_idx_from_combined(ij, n, i, j)
extracts indices i and j from combined index ij
integer function, public csr_idx_to_combined(i, j, n)
calculates combined index (ij)
integer function, dimension(2), public get_irange_csr(nindex, mp_group)
calculates index range for processor in group mp_group
subroutine, public create_active_space_type(active_space_env)
Creates an active space environment type, nullifying all quantities.
Contains utility routines for the active space module.
subroutine, public eri_to_array(eri_env, array, active_orbitals, spin1, spin2)
Copy the eri tensor for spins isp1 and isp2 to a standard 1D Fortran array.
subroutine, public subspace_matrix_to_array(source_matrix, target_array, row_index, col_index)
Copy a (square portion) of a cp_fm_type matrix to a standard 1D Fortran array.
Calculate the plane wave density by collocating the primitive Gaussian functions (pgf).
subroutine, public calculate_wavefunction(mo_vectors, ivector, rho, rho_gspace, atomic_kind_set, qs_kind_set, cell, dft_control, particle_set, pw_env, basis_type)
maps a given wavefunction on the grid
collects routines that calculate density matrices
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs)
Get the QUICKSTEP environment.
subroutine, public set_qs_env(qs_env, super_cell, mos, qmmm, qmmm_periodic, ewald_env, ewald_pw, mpools, rho_external, external_vxc, mask, scf_control, rel_control, qs_charges, ks_env, ks_qmmm_env, wf_history, scf_env, active_space, input, oce, rho_atom_set, rho0_atom_set, rho0_mpole, run_rtp, rtp, rhoz_set, rhoz_tot, ecoul_1c, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, efield, linres_control, xas_env, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, ls_scf_env, do_transport, transport_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, harris_env, gcp_env, mp2_env, bs_env, kg_env, force, kpoints, wanniercentres, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs)
Set the QUICKSTEP environment.
Integrate single or product functions over a potential on a RS grid.
Define the quickstep kind type and their sub types.
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 qs_ks_did_change(ks_env, s_mstruct_changed, rho_changed, potential_changed, full_reset)
tells that some of the things relevant to the ks calculation did change. has to be called when change...
Definition and initialisation of the mo data type.
subroutine, public write_mo_set_to_output_unit(mo_set, atomic_kind_set, qs_kind_set, particle_set, dft_section, before, kpoint, final_mos, spin, solver_method, rtp, cpart, sim_step, umo_set)
Write MO information to output file (eigenvalues, occupation numbers, coefficients)
collects routines that perform operations directly related to MOs
Definition and initialisation of the mo data type.
subroutine, public allocate_mo_set(mo_set, nao, nmo, nelectron, n_el_f, maxocc, flexible_electron_count)
Allocates a mo set and partially initializes it (nao,nmo,nelectron, and flexible_electron_count are v...
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count)
Get the components of a MO set data structure.
subroutine, public init_mo_set(mo_set, fm_pool, fm_ref, fm_struct, name)
initializes an allocated mo_set. eigenvalues, mo_coeff, occupation_numbers are valid only after this ...
Define the neighbor list data types and the corresponding functionality.
subroutine, public release_neighbor_list_sets(nlists)
releases an array of neighbor_list_sets
an eigen-space solver for the generalised symmetric eigenvalue problem for sparse matrices,...
subroutine, public ot_eigensolver(matrix_h, matrix_s, matrix_orthogonal_space_fm, matrix_c_fm, preconditioner, eps_gradient, iter_max, size_ortho_space, silent, ot_settings)
...
methods of the rho structure (defined in qs_rho_types)
subroutine, public qs_rho_update_rho(rho_struct, qs_env, rho_xc_external, local_rho_set, task_list_external, task_list_external_soft, pw_env_external, para_env_external)
updates rho_r and rho_g to the rhorho_ao. if use_kinetic_energy_density also computes tau_r and tau_g...
superstucture that hold various representations of the density and keeps track of which ones are vali...
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
types that represent a quickstep subsys
subroutine, public qs_subsys_get(subsys, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell, cell_ref, use_ref_cell, energy, force, qs_kind_set, cp_subsys, nelectron_total, nelectron_spin)
...
parameters that control an scf iteration
Implements UNIX and INET sockets.
generate the tasks lists used by collocate and integrate routines
subroutine, public generate_qs_task_list(ks_env, task_list, reorder_rs_grid_ranks, skip_load_balance_distributed, soft_valid, basis_type, pw_env_external, sab_orb_external)
...
subroutine, public deallocate_task_list(task_list)
deallocates the components and the object itself
subroutine, public allocate_task_list(task_list)
allocates and initialised the components of the task_list_type
All kind of helpful little routines.
pure integer function, dimension(2), public get_limit(m, n, me)
divide m entries into n parts, return size of part me
void open_bind_socket(int *psockfd, int *inet, int *port, char *host)
Opens and binds a socket.
void writebuffer(int *psockfd, char *data, int *plen)
Writes to a socket.
void accept_socket(int *psockfd, int *pclientfd)
Listens to a socket.
void close_socket(int *psockfd)
Closes a socket.
void listen_socket(int *psockfd, int *backlog)
Listens to a socket.
void readbuffer(int *psockfd, char *data, int *plen)
Reads from a socket.
void remove_socket_file(char *host)
Removes a socket file.
Provides all information about an atomic kind.
Type defining parameters related to the simulation cell.
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
keeps the information about the structure of a full matrix
type of a logger, at the moment it contains just a print level starting at which level it should be l...
stores all the informations relevant to an mpi environment
represent a list of objects
contained for different pw related things
contains all the informations needed by the fft based poisson solvers
environment for the poisson solver
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
The full QCSchema output type. For more information refer to: https://molssi-qc-schema....
Abstract function object for the eri_type_eri_foreach method.
Provides all information about a quickstep kind.
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
keeps the density in various representations, keeping track of which ones are valid.