200#include "./base/base_uses.f90"
208 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'energy_corrections'
229 LOGICAL,
INTENT(IN),
OPTIONAL :: ec_init, calculate_forces
231 CHARACTER(len=*),
PARAMETER :: routinen =
'energy_correction'
233 INTEGER :: handle, unit_nr
234 LOGICAL :: my_calc_forces
241 CALL timeset(routinen, handle)
244 IF (logger%para_env%is_source())
THEN
256 IF (.NOT. ec_env%do_skip)
THEN
258 ec_env%should_update = .true.
259 IF (
PRESENT(ec_init)) ec_env%should_update = ec_init
261 my_calc_forces = .false.
262 IF (
PRESENT(calculate_forces)) my_calc_forces = calculate_forces
264 IF (ec_env%should_update)
THEN
265 ec_env%old_etotal = 0.0_dp
266 ec_env%etotal = 0.0_dp
267 ec_env%eband = 0.0_dp
268 ec_env%ehartree = 0.0_dp
272 ec_env%edispersion = 0.0_dp
273 ec_env%exc_aux_fit = 0.0_dp
277 ec_env%old_etotal = energy%total
281 IF (my_calc_forces)
THEN
282 IF (unit_nr > 0)
THEN
283 WRITE (unit_nr,
'(T2,A,A,A,A,A)')
"!", repeat(
"-", 25), &
284 " Energy Correction Forces ", repeat(
"-", 26),
"!"
286 CALL get_qs_env(qs_env, force=ks_force, virial=virial)
290 IF (unit_nr > 0)
THEN
291 WRITE (unit_nr,
'(T2,A,A,A,A,A)')
"!", repeat(
"-", 29), &
292 " Energy Correction ", repeat(
"-", 29),
"!"
297 CALL energy_correction_low(qs_env, ec_env, my_calc_forces, unit_nr)
300 IF (ec_env%should_update)
THEN
301 energy%nonscf_correction = ec_env%etotal - ec_env%old_etotal
302 energy%total = ec_env%etotal
305 IF (.NOT. my_calc_forces .AND. unit_nr > 0)
THEN
306 WRITE (unit_nr,
'(T3,A,T56,F25.15)')
"Energy Correction ", energy%nonscf_correction
308 IF (unit_nr > 0)
THEN
309 WRITE (unit_nr,
'(T2,A,A,A)')
"!", repeat(
"-", 77),
"!"
316 IF (unit_nr > 0)
THEN
317 WRITE (unit_nr,
'(T2,A,A,A)')
"!", repeat(
"-", 77),
"!"
318 WRITE (unit_nr,
'(T2,A,A,A,A,A)')
"!", repeat(
"-", 26), &
319 " Skip Energy Correction ", repeat(
"-", 27),
"!"
320 WRITE (unit_nr,
'(T2,A,A,A)')
"!", repeat(
"-", 77),
"!"
325 CALL timestop(handle)
340 SUBROUTINE energy_correction_low(qs_env, ec_env, calculate_forces, unit_nr)
343 LOGICAL,
INTENT(IN) :: calculate_forces
344 INTEGER,
INTENT(IN) :: unit_nr
346 INTEGER :: ispin, nkind, nspins
347 LOGICAL :: debug_f, gapw, gapw_xc
348 REAL(kind=
dp) :: eps_fit, exc
356 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
358 IF (ec_env%should_update)
THEN
359 CALL ec_build_neighborlist(qs_env, ec_env)
363 ec_env%vtau_rspace, &
364 ec_env%vadmm_rspace, &
365 ec_env%ehartree, exc)
367 SELECT CASE (ec_env%energy_functional)
370 CALL ec_build_core_hamiltonian(qs_env, ec_env)
371 CALL ec_build_ks_matrix(qs_env, ec_env)
376 NULLIFY (ec_env%mao_coef)
377 CALL mao_generate_basis(qs_env, ec_env%mao_coef, ref_basis_set=
"HARRIS", molecular=.true., &
378 max_iter=ec_env%mao_max_iter, eps_grad=ec_env%mao_eps_grad, &
379 eps1_mao=ec_env%mao_eps1, iolevel=ec_env%mao_iolevel, unit_nr=unit_nr)
382 CALL ec_ks_solver(qs_env, ec_env)
384 CALL evaluate_ec_core_matrix_traces(qs_env, ec_env)
389 CALL ec_dc_energy(qs_env, ec_env, calculate_forces=.false.)
394 CALL ec_build_ks_matrix(qs_env, ec_env)
401 cpabort(
"unknown energy correction")
405 CALL ec_disp(qs_env, ec_env, calculate_forces=.false.)
408 CALL ec_energy(ec_env, unit_nr)
412 IF (calculate_forces)
THEN
414 debug_f = ec_env%debug_forces .OR. ec_env%debug_stress
416 CALL ec_disp(qs_env, ec_env, calculate_forces=.true.)
418 SELECT CASE (ec_env%energy_functional)
421 CALL ec_build_core_hamiltonian_force(qs_env, ec_env, &
425 CALL ec_build_ks_matrix_force(qs_env, ec_env)
426 IF (ec_env%debug_external)
THEN
427 CALL write_response_interface(qs_env, ec_env)
428 CALL init_response_deriv(qs_env, ec_env)
435 CALL ec_dc_energy(qs_env, ec_env, calculate_forces=.true.)
437 CALL ec_build_core_hamiltonian_force(qs_env, ec_env, &
441 CALL ec_dc_build_ks_matrix_force(qs_env, ec_env)
442 IF (ec_env%debug_external)
THEN
443 CALL write_response_interface(qs_env, ec_env)
444 CALL init_response_deriv(qs_env, ec_env)
449 CALL init_response_deriv(qs_env, ec_env)
453 cpabort(
"unknown energy correction")
456 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
457 nspins = dft_control%nspins
458 gapw = dft_control%qs_control%gapw
459 gapw_xc = dft_control%qs_control%gapw_xc
460 IF (gapw .OR. gapw_xc)
THEN
462 qs_kind_set=qs_kind_set, particle_set=particle_set)
463 CALL get_qs_env(qs_env=qs_env, oce=oce, sap_oce=sap_oce)
466 eps_fit = dft_control%qs_control%gapw_control%eps_fit
478 cpassert(
ASSOCIATED(pw_env))
479 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
480 ALLOCATE (ec_env%rhoz_r(nspins))
482 CALL auxbas_pw_pool%create_pw(ec_env%rhoz_r(ispin))
486 vh_rspace=ec_env%vh_rspace, &
487 vxc_rspace=ec_env%vxc_rspace, &
488 vtau_rspace=ec_env%vtau_rspace, &
489 vadmm_rspace=ec_env%vadmm_rspace, &
490 matrix_hz=ec_env%matrix_hz, &
491 matrix_pz=ec_env%matrix_z, &
492 matrix_pz_admm=ec_env%z_admm, &
493 matrix_wz=ec_env%matrix_wz, &
494 rhopz_r=ec_env%rhoz_r, &
495 zehartree=ec_env%ehartree, &
497 zexc_aux_fit=ec_env%exc_aux_fit, &
498 p_env=ec_env%p_env, &
501 CALL output_response_deriv(qs_env, ec_env, unit_nr)
503 CALL ec_properties(qs_env, ec_env)
506 IF (
ASSOCIATED(ec_env%rhoout_r))
THEN
508 CALL auxbas_pw_pool%give_back_pw(ec_env%rhoout_r(ispin))
510 DEALLOCATE (ec_env%rhoout_r)
512 IF (
ASSOCIATED(ec_env%rhoz_r))
THEN
514 CALL auxbas_pw_pool%give_back_pw(ec_env%rhoz_r(ispin))
516 DEALLOCATE (ec_env%rhoz_r)
532 END SUBROUTINE energy_correction_low
540 SUBROUTINE write_response_interface(qs_env, ec_env)
547 NULLIFY (trexio_section)
551 CALL write_trexio(qs_env, trexio_section, ec_env%matrix_hz)
553 END SUBROUTINE write_response_interface
561 SUBROUTINE init_response_deriv(qs_env, ec_env)
571 ALLOCATE (ec_env%rf(3, natom))
574 CALL get_qs_env(qs_env, force=force, virial=virial)
576 CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set)
579 IF (virial%pv_availability .AND. (.NOT. virial%pv_numer))
THEN
580 ec_env%rpv = virial%pv_virial
583 END SUBROUTINE init_response_deriv
592 SUBROUTINE output_response_deriv(qs_env, ec_env, unit_nr)
595 INTEGER,
INTENT(IN) :: unit_nr
597 CHARACTER(LEN=default_string_length) :: unit_string
598 INTEGER :: funit, ia, natom
599 REAL(kind=
dp) :: evol, fconv
600 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: ftot
608 IF (
ASSOCIATED(ec_env%rf))
THEN
610 ALLOCATE (ftot(3, natom))
612 CALL get_qs_env(qs_env, force=force, virial=virial, para_env=para_env)
614 CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set)
616 ec_env%rf(1:3, 1:natom) = ftot(1:3, 1:natom) - ec_env%rf(1:3, 1:natom)
617 CALL para_env%sum(ec_env%rf)
620 IF (virial%pv_availability .AND. (.NOT. virial%pv_numer))
THEN
621 ec_env%rpv = virial%pv_virial - ec_env%rpv
622 CALL para_env%sum(ec_env%rpv)
624 evol = ec_env%exc + ec_env%exc_aux_fit + 2.0_dp*ec_env%ehartree
625 ec_env%rpv(1, 1) = ec_env%rpv(1, 1) - evol
626 ec_env%rpv(2, 2) = ec_env%rpv(2, 2) - evol
627 ec_env%rpv(3, 3) = ec_env%rpv(3, 3) - evol
630 CALL get_qs_env(qs_env, particle_set=particle_set, cell=cell)
634 IF (unit_nr > 0)
THEN
635 WRITE (unit_nr,
'(T2,A)')
"Write EXTERNAL Response Derivative: "//trim(ec_env%exresult_fn)
637 CALL open_file(ec_env%exresult_fn, file_status=
"REPLACE", file_form=
"FORMATTED", &
638 file_action=
"WRITE", unit_number=funit)
639 WRITE (funit,
"(T8,A,T58,A)")
"COORDINATES [Bohr]",
"RESPONSE FORCES [Hartree/Bohr]"
641 WRITE (funit,
"(2(3F15.8,5x))") particle_set(ia)%r(1:3), ec_env%rf(1:3, ia)
644 WRITE (funit,
"(T8,A,T58,A)")
"CELL [Bohr]",
"RESPONSE PRESSURE [GPa]"
646 WRITE (funit,
"(3F15.8,5x,3F15.8)") cell%hmat(ia, 1:3), -fconv*ec_env%rpv(ia, 1:3)
653 END SUBROUTINE output_response_deriv
662 SUBROUTINE evaluate_ec_core_matrix_traces(qs_env, ec_env)
666 CHARACTER(LEN=*),
PARAMETER :: routinen =
'evaluate_ec_core_matrix_traces'
672 CALL timeset(routinen, handle)
675 CALL get_qs_env(qs_env, dft_control=dft_control, energy=energy)
678 CALL calculate_ptrace(ec_env%matrix_h, ec_env%matrix_p, energy%core, dft_control%nspins)
681 CALL calculate_ptrace(ec_env%matrix_t, ec_env%matrix_p, energy%kinetic, dft_control%nspins)
683 CALL timestop(handle)
685 END SUBROUTINE evaluate_ec_core_matrix_traces
698 SUBROUTINE ec_dc_energy(qs_env, ec_env, calculate_forces)
701 LOGICAL,
INTENT(IN) :: calculate_forces
703 CHARACTER(LEN=*),
PARAMETER :: routinen =
'ec_dc_energy'
705 CHARACTER(LEN=default_string_length) :: headline
706 INTEGER :: handle, ispin, nspins
707 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_h, matrix_p, matrix_s, matrix_w
712 CALL timeset(routinen, handle)
714 NULLIFY (dft_control, ks_env, matrix_h, matrix_p, matrix_s, matrix_w, rho)
716 dft_control=dft_control, &
718 matrix_h_kp=matrix_h, &
719 matrix_s_kp=matrix_s, &
720 matrix_w_kp=matrix_w, &
723 nspins = dft_control%nspins
728 matrix_name=
"OVERLAP MATRIX", &
729 basis_type_a=
"HARRIS", &
730 basis_type_b=
"HARRIS", &
731 sab_nl=ec_env%sab_orb)
736 headline =
"CORE HAMILTONIAN MATRIX"
737 ALLOCATE (ec_env%matrix_h(1, 1)%matrix)
738 CALL dbcsr_create(ec_env%matrix_h(1, 1)%matrix, name=trim(headline), &
739 template=matrix_s(1, 1)%matrix, matrix_type=dbcsr_type_symmetric)
741 CALL dbcsr_copy(ec_env%matrix_h(1, 1)%matrix, matrix_h(1, 1)%matrix)
746 headline =
"DENSITY MATRIX"
748 ALLOCATE (ec_env%matrix_p(ispin, 1)%matrix)
749 CALL dbcsr_create(ec_env%matrix_p(ispin, 1)%matrix, name=trim(headline), &
750 template=matrix_s(1, 1)%matrix, matrix_type=dbcsr_type_symmetric)
752 CALL dbcsr_copy(ec_env%matrix_p(ispin, 1)%matrix, matrix_p(ispin, 1)%matrix)
755 IF (calculate_forces)
THEN
760 headline =
"ENERGY-WEIGHTED DENSITY MATRIX"
762 ALLOCATE (ec_env%matrix_w(ispin, 1)%matrix)
763 CALL dbcsr_create(ec_env%matrix_w(ispin, 1)%matrix, name=trim(headline), &
764 template=matrix_s(1, 1)%matrix, matrix_type=dbcsr_type_symmetric)
766 CALL dbcsr_copy(ec_env%matrix_w(ispin, 1)%matrix, matrix_w(ispin, 1)%matrix)
772 ec_env%efield_nuclear = 0.0_dp
773 ec_env%efield_elec = 0.0_dp
776 CALL timestop(handle)
778 END SUBROUTINE ec_dc_energy
789 SUBROUTINE ec_dc_build_ks_matrix_force(qs_env, ec_env)
793 CHARACTER(LEN=*),
PARAMETER :: routinen =
'ec_dc_build_ks_matrix_force'
795 CHARACTER(LEN=default_string_length) :: unit_string
796 INTEGER :: handle, i, iounit, ispin, natom, nspins
797 LOGICAL :: debug_forces, debug_stress, do_ec_hfx, &
799 REAL(
dp) :: dummy_real, dummy_real2(2), ehartree, &
801 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: ftot
802 REAL(
dp),
DIMENSION(3) :: fodeb, fodeb2
803 REAL(kind=
dp),
DIMENSION(3, 3) :: h_stress, pv_loc, stdeb, sttot
807 TYPE(
dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_ks, matrix_s, scrm
808 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_p
819 TYPE(
pw_r3d_rs_type),
DIMENSION(:),
POINTER :: rho_r, v_rspace, v_rspace_in, &
827 CALL timeset(routinen, handle)
829 debug_forces = ec_env%debug_forces
830 debug_stress = ec_env%debug_stress
833 IF (logger%para_env%is_source())
THEN
839 NULLIFY (atomic_kind_set, cell, dft_control, force, ks_env, matrix_ks, &
840 matrix_p, matrix_s, para_env, pw_env, rho, sab_orb, virial)
843 dft_control=dft_control, &
846 matrix_ks=matrix_ks, &
853 cpassert(
ASSOCIATED(pw_env))
855 nspins = dft_control%nspins
856 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
858 fconv = 1.0e-9_dp*
pascal/cell%deth
859 IF (debug_stress .AND. use_virial)
THEN
860 sttot = virial%pv_virial
866 NULLIFY (auxbas_pw_pool, poisson_env)
868 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
869 poisson_env=poisson_env)
872 CALL auxbas_pw_pool%create_pw(v_hartree_gspace)
873 CALL auxbas_pw_pool%create_pw(rho_tot_gspace)
874 CALL auxbas_pw_pool%create_pw(v_hartree_rspace)
883 h_stress(:, :) = 0.0_dp
885 density=rho_tot_gspace, &
887 vhartree=v_hartree_gspace, &
890 virial%pv_ehartree = virial%pv_ehartree + h_stress/real(para_env%num_pe,
dp)
891 virial%pv_virial = virial%pv_virial + h_stress/real(para_env%num_pe,
dp)
893 IF (debug_stress)
THEN
894 stdeb = fconv*(h_stress/real(para_env%num_pe,
dp))
895 CALL para_env%sum(stdeb)
896 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
905 CALL pw_transfer(v_hartree_gspace, v_hartree_rspace)
906 CALL pw_scale(v_hartree_rspace, v_hartree_rspace%pw_grid%dvol)
910 ALLOCATE (ec_env%rhoout_r(nspins))
912 CALL auxbas_pw_pool%create_pw(ec_env%rhoout_r(ispin))
913 CALL pw_copy(rho_r(ispin), ec_env%rhoout_r(ispin))
918 IF (debug_forces) fodeb(1:3) = force(1)%rho_core(1:3, 1)
919 IF (debug_stress .AND. use_virial) stdeb = virial%pv_ehartree
920 CALL integrate_v_core_rspace(v_hartree_rspace, qs_env)
921 IF (debug_forces)
THEN
922 fodeb(1:3) = force(1)%rho_core(1:3, 1) - fodeb(1:3)
923 CALL para_env%sum(fodeb)
924 IF (iounit > 0)
WRITE (iounit,
"(T3,A,T33,3F16.8)")
"DEBUG:: Vtot*dncore", fodeb
926 IF (debug_stress .AND. use_virial)
THEN
927 stdeb = fconv*(virial%pv_ehartree - stdeb)
928 CALL para_env%sum(stdeb)
929 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
935 NULLIFY (v_rspace, v_tau_rspace)
938 IF (use_virial) virial%pv_calculate = .true.
941 CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho, xc_section=ec_env%xc_section, &
942 vxc_rho=v_rspace, vxc_tau=v_tau_rspace, exc=exc, just_energy=.false.)
944 IF (.NOT.
ASSOCIATED(v_rspace))
THEN
945 ALLOCATE (v_rspace(nspins))
947 CALL auxbas_pw_pool%create_pw(v_rspace(ispin))
953 virial%pv_exc = virial%pv_exc - virial%pv_xc
954 virial%pv_virial = virial%pv_virial - virial%pv_xc
962 ALLOCATE (scrm(ispin)%matrix)
963 CALL dbcsr_create(scrm(ispin)%matrix, template=ec_env%matrix_ks(ispin, 1)%matrix)
964 CALL dbcsr_copy(scrm(ispin)%matrix, ec_env%matrix_ks(ispin, 1)%matrix)
965 CALL dbcsr_set(scrm(ispin)%matrix, 0.0_dp)
968 pw_grid => v_hartree_rspace%pw_grid
969 ALLOCATE (v_rspace_in(nspins))
971 CALL v_rspace_in(ispin)%create(pw_grid)
977 CALL pw_transfer(ec_env%vxc_rspace(ispin), v_rspace_in(ispin))
979 CALL pw_axpy(ec_env%vh_rspace, v_rspace_in(ispin))
990 IF (debug_forces) fodeb(1:3) = force(1)%fock_4c(1:3, 1)
991 IF (debug_forces) fodeb2(1:3) = force(1)%overlap_admm(1:3, 1)
998 hfx_sections=ec_hfx_sections, &
999 x_data=ec_env%x_data, &
1001 do_admm=ec_env%do_ec_admm, &
1002 calc_forces=.true., &
1003 reuse_hfx=ec_env%reuse_hfx, &
1004 do_im_time=.false., &
1005 e_ex_from_gw=dummy_real, &
1006 e_admm_from_gw=dummy_real2, &
1009 IF (debug_forces)
THEN
1010 fodeb(1:3) = force(1)%fock_4c(1:3, 1) - fodeb(1:3)
1011 CALL para_env%sum(fodeb)
1012 IF (iounit > 0)
WRITE (iounit,
"(T3,A,T33,3F16.8)")
"DEBUG:: P*hfx_DC ", fodeb
1014 fodeb2(1:3) = force(1)%overlap_admm(1:3, 1) - fodeb2(1:3)
1015 CALL para_env%sum(fodeb2)
1016 IF (iounit > 0)
WRITE (iounit,
"(T3,A,T33,3F16.8)")
"DEBUG:: P*hfx_DC*S ", fodeb2
1018 IF (debug_stress .AND. use_virial)
THEN
1019 stdeb = -1.0_dp*fconv*virial%pv_fock_4c
1020 CALL para_env%sum(stdeb)
1021 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
1031 IF (use_virial)
THEN
1032 pv_loc = virial%pv_virial
1035 IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
1036 IF (debug_stress .AND. use_virial) stdeb = virial%pv_virial
1038 DO ispin = 1, nspins
1040 CALL pw_scale(v_rspace(ispin), v_rspace(ispin)%pw_grid%dvol)
1041 CALL pw_axpy(v_hartree_rspace, v_rspace(ispin))
1043 CALL integrate_v_rspace(v_rspace=v_rspace(ispin), &
1045 pmat=matrix_p(ispin, 1), &
1047 calculate_forces=.true., &
1048 basis_type=
"HARRIS", &
1049 task_list_external=ec_env%task_list)
1052 IF (debug_forces)
THEN
1053 fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
1054 CALL para_env%sum(fodeb)
1055 IF (iounit > 0)
WRITE (iounit,
"(T3,A,T33,3F16.8)")
"DEBUG:: Pout*dVhxc ", fodeb
1057 IF (debug_stress .AND. use_virial)
THEN
1058 stdeb = fconv*(virial%pv_virial - stdeb)
1059 CALL para_env%sum(stdeb)
1060 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
1064 IF (
ASSOCIATED(v_tau_rspace))
THEN
1065 IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
1066 IF (debug_stress .AND. use_virial) stdeb = virial%pv_virial
1067 DO ispin = 1, nspins
1068 CALL pw_scale(v_tau_rspace(ispin), v_tau_rspace(ispin)%pw_grid%dvol)
1070 CALL integrate_v_rspace(v_rspace=v_tau_rspace(ispin), &
1072 pmat=matrix_p(ispin, 1), &
1074 calculate_forces=.true., &
1075 compute_tau=.true., &
1076 basis_type=
"HARRIS", &
1077 task_list_external=ec_env%task_list)
1080 IF (debug_forces)
THEN
1081 fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
1082 CALL para_env%sum(fodeb)
1083 IF (iounit > 0)
WRITE (iounit,
"(T3,A,T33,3F16.8)")
"DEBUG:: Pout*dVhxc_tau ", fodeb
1085 IF (debug_stress .AND. use_virial)
THEN
1086 stdeb = fconv*(virial%pv_virial - stdeb)
1087 CALL para_env%sum(stdeb)
1088 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
1094 IF (use_virial)
THEN
1095 virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
1114 NULLIFY (ec_env%matrix_hz)
1116 DO ispin = 1, nspins
1117 ALLOCATE (ec_env%matrix_hz(ispin)%matrix)
1118 CALL dbcsr_create(ec_env%matrix_hz(ispin)%matrix, template=matrix_s(1)%matrix)
1119 CALL dbcsr_copy(ec_env%matrix_hz(ispin)%matrix, matrix_s(1)%matrix)
1120 CALL dbcsr_set(ec_env%matrix_hz(ispin)%matrix, 0.0_dp)
1123 DO ispin = 1, nspins
1126 CALL pw_axpy(v_rspace_in(ispin), v_rspace(ispin), -1.0_dp)
1129 DO ispin = 1, nspins
1130 CALL integrate_v_rspace(v_rspace=v_rspace(ispin), &
1131 hmat=ec_env%matrix_hz(ispin), &
1132 pmat=matrix_p(ispin, 1), &
1134 calculate_forces=.false., &
1135 basis_type=
"HARRIS", &
1136 task_list_external=ec_env%task_list)
1140 IF (dft_control%use_kinetic_energy_density)
THEN
1143 IF (.NOT.
ASSOCIATED(v_tau_rspace))
THEN
1144 ALLOCATE (v_tau_rspace(nspins))
1145 DO ispin = 1, nspins
1146 CALL auxbas_pw_pool%create_pw(v_tau_rspace(ispin))
1147 CALL pw_zero(v_tau_rspace(ispin))
1151 DO ispin = 1, nspins
1153 IF (
ASSOCIATED(ec_env%vtau_rspace))
THEN
1154 CALL pw_axpy(ec_env%vtau_rspace(ispin), v_tau_rspace(ispin), -1.0_dp)
1157 CALL integrate_v_rspace(v_rspace=v_tau_rspace(ispin), &
1158 hmat=ec_env%matrix_hz(ispin), &
1159 pmat=matrix_p(ispin, 1), &
1161 calculate_forces=.false., compute_tau=.true., &
1162 basis_type=
"HARRIS", &
1163 task_list_external=ec_env%task_list)
1171 ext_hfx_section=ec_hfx_sections, &
1172 x_data=ec_env%x_data, &
1173 recalc_integrals=.false., &
1174 do_admm=ec_env%do_ec_admm, &
1177 reuse_hfx=ec_env%reuse_hfx)
1180 IF (debug_forces) fodeb(1:3) = force(1)%core_overlap(1:3, 1)
1181 IF (debug_stress .AND. use_virial) stdeb = virial%pv_ecore_overlap
1183 IF (debug_forces)
THEN
1184 fodeb(1:3) = force(1)%core_overlap(1:3, 1) - fodeb(1:3)
1185 CALL para_env%sum(fodeb)
1186 IF (iounit > 0)
WRITE (iounit,
"(T3,A,T33,3F16.8)")
"DEBUG:: CoreOverlap", fodeb
1188 IF (debug_stress .AND. use_virial)
THEN
1189 stdeb = fconv*(stdeb - virial%pv_ecore_overlap)
1190 CALL para_env%sum(stdeb)
1191 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
1195 IF (debug_forces)
THEN
1196 CALL get_qs_env(qs_env, natom=natom, atomic_kind_set=atomic_kind_set)
1197 ALLOCATE (ftot(3, natom))
1199 fodeb(1:3) = ftot(1:3, 1)
1201 CALL para_env%sum(fodeb)
1202 IF (iounit > 0)
WRITE (iounit,
"(T3,A,T33,3F16.8)")
"DEBUG:: Force Explicit", fodeb
1206 DO ispin = 1, nspins
1207 CALL auxbas_pw_pool%give_back_pw(v_rspace(ispin))
1208 CALL auxbas_pw_pool%give_back_pw(v_rspace_in(ispin))
1209 IF (
ASSOCIATED(v_tau_rspace))
THEN
1210 CALL auxbas_pw_pool%give_back_pw(v_tau_rspace(ispin))
1214 DEALLOCATE (v_rspace, v_rspace_in)
1215 IF (
ASSOCIATED(v_tau_rspace))
DEALLOCATE (v_tau_rspace)
1217 CALL auxbas_pw_pool%give_back_pw(rho_tot_gspace)
1218 CALL auxbas_pw_pool%give_back_pw(v_hartree_gspace)
1219 CALL auxbas_pw_pool%give_back_pw(v_hartree_rspace)
1223 IF (use_virial)
THEN
1224 IF (qs_env%energy_correction)
THEN
1225 ec_env%ehartree = ehartree
1230 IF (debug_stress .AND. use_virial)
THEN
1232 stdeb = -1.0_dp*fconv*ehartree
1233 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
1236 stdeb = -1.0_dp*fconv*exc
1237 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
1246 CALL para_env%sum(virdeb%pv_overlap)
1247 CALL para_env%sum(virdeb%pv_ekinetic)
1248 CALL para_env%sum(virdeb%pv_ppl)
1249 CALL para_env%sum(virdeb%pv_ppnl)
1250 CALL para_env%sum(virdeb%pv_ecore_overlap)
1251 CALL para_env%sum(virdeb%pv_ehartree)
1252 CALL para_env%sum(virdeb%pv_exc)
1253 CALL para_env%sum(virdeb%pv_exx)
1254 CALL para_env%sum(virdeb%pv_vdw)
1255 CALL para_env%sum(virdeb%pv_mp2)
1256 CALL para_env%sum(virdeb%pv_nlcc)
1257 CALL para_env%sum(virdeb%pv_gapw)
1258 CALL para_env%sum(virdeb%pv_lrigpw)
1259 CALL para_env%sum(virdeb%pv_virial)
1264 virdeb%pv_ehartree(i, i) = virdeb%pv_ehartree(i, i) - 2.0_dp*ehartree
1265 virdeb%pv_virial(i, i) = virdeb%pv_virial(i, i) - exc &
1267 virdeb%pv_exc(i, i) = virdeb%pv_exc(i, i) - exc
1273 CALL para_env%sum(sttot)
1274 stdeb = fconv*(virdeb%pv_virial - sttot)
1275 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
1278 stdeb = fconv*(virdeb%pv_virial)
1279 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
1289 CALL timestop(handle)
1291 END SUBROUTINE ec_dc_build_ks_matrix_force
1299 SUBROUTINE ec_disp(qs_env, ec_env, calculate_forces)
1300 TYPE(qs_environment_type),
POINTER :: qs_env
1301 TYPE(energy_correction_type),
POINTER :: ec_env
1302 LOGICAL,
INTENT(IN) :: calculate_forces
1304 REAL(kind=dp) :: edisp
1306 CALL calculate_dispersion_pairpot(qs_env, ec_env%dispersion_env, edisp, calculate_forces)
1307 IF (.NOT. calculate_forces) ec_env%edispersion = ec_env%edispersion + edisp
1309 END SUBROUTINE ec_disp
1318 SUBROUTINE ec_build_core_hamiltonian(qs_env, ec_env)
1319 TYPE(qs_environment_type),
POINTER :: qs_env
1320 TYPE(energy_correction_type),
POINTER :: ec_env
1322 CHARACTER(LEN=*),
PARAMETER :: routinen =
'ec_build_core_hamiltonian'
1324 CHARACTER(LEN=default_string_length) :: basis_type
1325 INTEGER :: handle, nder, nimages
1326 LOGICAL :: calculate_forces, use_virial
1327 TYPE(atomic_kind_type),
DIMENSION(:),
POINTER :: atomic_kind_set
1328 TYPE(dft_control_type),
POINTER :: dft_control
1329 TYPE(neighbor_list_set_p_type),
DIMENSION(:), &
1330 POINTER :: sab_orb, sac_ae, sac_ppl, sap_ppnl
1331 TYPE(particle_type),
DIMENSION(:),
POINTER :: particle_set
1332 TYPE(qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
1333 TYPE(qs_ks_env_type),
POINTER :: ks_env
1335 CALL timeset(routinen, handle)
1337 NULLIFY (atomic_kind_set, dft_control, ks_env, particle_set, &
1340 CALL get_qs_env(qs_env=qs_env, &
1341 atomic_kind_set=atomic_kind_set, &
1342 dft_control=dft_control, &
1343 particle_set=particle_set, &
1344 qs_kind_set=qs_kind_set, &
1348 nimages = dft_control%nimages
1349 IF (nimages /= 1)
THEN
1350 cpabort(
"K-points for Harris functional not implemented")
1354 IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc)
THEN
1355 cpabort(
"Harris functional for GAPW not implemented")
1359 use_virial = .false.
1360 calculate_forces = .false.
1363 NULLIFY (sab_orb, sac_ae, sac_ppl, sap_ppnl)
1364 sab_orb => ec_env%sab_orb
1365 sac_ae => ec_env%sac_ae
1366 sac_ppl => ec_env%sac_ppl
1367 sap_ppnl => ec_env%sap_ppnl
1369 basis_type =
"HARRIS"
1373 CALL build_overlap_matrix(ks_env, matrixkp_s=ec_env%matrix_s, &
1374 matrix_name=
"OVERLAP MATRIX", &
1375 basis_type_a=basis_type, &
1376 basis_type_b=basis_type, &
1378 CALL build_kinetic_matrix(ks_env, matrixkp_t=ec_env%matrix_t, &
1379 matrix_name=
"KINETIC ENERGY MATRIX", &
1380 basis_type=basis_type, &
1384 CALL dbcsr_allocate_matrix_set(ec_env%matrix_h, 1, 1)
1385 ALLOCATE (ec_env%matrix_h(1, 1)%matrix)
1386 CALL dbcsr_create(ec_env%matrix_h(1, 1)%matrix, template=ec_env%matrix_s(1, 1)%matrix)
1387 CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_h(1, 1)%matrix, sab_orb)
1390 CALL dbcsr_copy(ec_env%matrix_h(1, 1)%matrix, ec_env%matrix_t(1, 1)%matrix, &
1391 keep_sparsity=.true., name=
"CORE HAMILTONIAN MATRIX")
1393 CALL core_matrices(qs_env, ec_env%matrix_h, ec_env%matrix_p, calculate_forces, nder, &
1394 ec_env=ec_env, ec_env_matrices=.true., basis_type=basis_type)
1397 ec_env%efield_nuclear = 0.0_dp
1398 CALL ec_efield_local_operator(qs_env, ec_env, calculate_forces)
1400 CALL timestop(handle)
1402 END SUBROUTINE ec_build_core_hamiltonian
1413 SUBROUTINE ec_build_ks_matrix(qs_env, ec_env)
1414 TYPE(qs_environment_type),
POINTER :: qs_env
1415 TYPE(energy_correction_type),
POINTER :: ec_env
1417 CHARACTER(LEN=*),
PARAMETER :: routinen =
'ec_build_ks_matrix'
1419 CHARACTER(LEN=default_string_length) :: headline
1420 INTEGER :: handle, iounit, ispin, nspins
1421 LOGICAL :: calculate_forces, &
1422 do_adiabatic_rescaling, do_ec_hfx, &
1423 hfx_treat_lsd_in_core, use_virial
1424 REAL(dp) :: dummy_real, dummy_real2(2), eexc, evhxc, &
1426 TYPE(cp_logger_type),
POINTER :: logger
1427 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: ks_mat
1428 TYPE(dft_control_type),
POINTER :: dft_control
1429 TYPE(pw_env_type),
POINTER :: pw_env
1430 TYPE(pw_pool_type),
POINTER :: auxbas_pw_pool
1431 TYPE(pw_r3d_rs_type),
DIMENSION(:),
POINTER :: rho_r, tau_r, v_rspace, v_tau_rspace
1432 TYPE(qs_energy_type),
POINTER :: energy
1433 TYPE(qs_ks_env_type),
POINTER :: ks_env
1434 TYPE(qs_rho_type),
POINTER :: rho
1435 TYPE(section_vals_type),
POINTER :: adiabatic_rescaling_section, &
1436 ec_hfx_sections, ec_section
1438 CALL timeset(routinen, handle)
1440 logger => cp_get_default_logger()
1441 IF (logger%para_env%is_source())
THEN
1442 iounit = cp_logger_get_default_unit_nr(logger, local=.true.)
1448 NULLIFY (auxbas_pw_pool, dft_control, energy, ks_env, rho, rho_r, tau_r)
1449 CALL get_qs_env(qs_env=qs_env, &
1450 dft_control=dft_control, &
1453 nspins = dft_control%nspins
1454 calculate_forces = .false.
1455 use_virial = .false.
1458 IF (
ASSOCIATED(ec_env%matrix_ks))
CALL dbcsr_deallocate_matrix_set(ec_env%matrix_ks)
1459 CALL dbcsr_allocate_matrix_set(ec_env%matrix_ks, nspins, 1)
1460 DO ispin = 1, nspins
1461 headline =
"KOHN-SHAM MATRIX"
1462 ALLOCATE (ec_env%matrix_ks(ispin, 1)%matrix)
1463 CALL dbcsr_create(ec_env%matrix_ks(ispin, 1)%matrix, name=trim(headline), &
1464 template=ec_env%matrix_s(1, 1)%matrix, matrix_type=dbcsr_type_symmetric)
1465 CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_ks(ispin, 1)%matrix, ec_env%sab_orb)
1466 CALL dbcsr_set(ec_env%matrix_ks(ispin, 1)%matrix, 0.0_dp)
1470 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
1471 cpassert(
ASSOCIATED(pw_env))
1474 ec_section => section_vals_get_subs_vals(qs_env%input,
"DFT%ENERGY_CORRECTION")
1475 ec_hfx_sections => section_vals_get_subs_vals(ec_section,
"XC%HF")
1476 CALL section_vals_get(ec_hfx_sections, explicit=do_ec_hfx)
1481 adiabatic_rescaling_section => section_vals_get_subs_vals(ec_section,
"XC%ADIABATIC_RESCALING")
1482 CALL section_vals_get(adiabatic_rescaling_section, explicit=do_adiabatic_rescaling)
1483 IF (do_adiabatic_rescaling)
THEN
1484 CALL cp_abort(__location__,
"Adiabatic rescaling NYI for energy correction")
1486 CALL section_vals_val_get(ec_hfx_sections,
"TREAT_LSD_IN_CORE", l_val=hfx_treat_lsd_in_core)
1487 IF (hfx_treat_lsd_in_core)
THEN
1488 CALL cp_abort(__location__,
"HFX_TREAT_LSD_IN_CORE NYI for energy correction")
1492 IF (dft_control%do_admm)
THEN
1494 IF (dft_control%do_admm_mo)
THEN
1495 IF (qs_env%run_rtp)
THEN
1496 CALL rtp_admm_calc_rho_aux(qs_env)
1498 CALL admm_mo_calc_rho_aux(qs_env)
1500 ELSEIF (dft_control%do_admm_dm)
THEN
1501 CALL admm_dm_calc_rho_aux(qs_env)
1508 CALL get_qs_env(qs_env, energy=energy)
1509 CALL calculate_exx(qs_env=qs_env, &
1511 hfx_sections=ec_hfx_sections, &
1512 x_data=ec_env%x_data, &
1514 do_admm=ec_env%do_ec_admm, &
1515 calc_forces=.false., &
1516 reuse_hfx=ec_env%reuse_hfx, &
1517 do_im_time=.false., &
1518 e_ex_from_gw=dummy_real, &
1519 e_admm_from_gw=dummy_real2, &
1523 ec_env%ex = energy%ex
1525 IF (ec_env%do_ec_admm)
THEN
1526 ec_env%exc_aux_fit = energy%exc_aux_fit + energy%exc
1532 ks_mat => ec_env%matrix_ks(:, 1)
1533 CALL add_exx_to_rhs(rhs=ks_mat, &
1535 ext_hfx_section=ec_hfx_sections, &
1536 x_data=ec_env%x_data, &
1537 recalc_integrals=.false., &
1538 do_admm=ec_env%do_ec_admm, &
1541 reuse_hfx=ec_env%reuse_hfx)
1546 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
1547 NULLIFY (v_rspace, v_tau_rspace)
1548 CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho, xc_section=ec_env%xc_section, &
1549 vxc_rho=v_rspace, vxc_tau=v_tau_rspace, exc=eexc, just_energy=.false.)
1551 IF (.NOT.
ASSOCIATED(v_rspace))
THEN
1552 ALLOCATE (v_rspace(nspins))
1553 DO ispin = 1, nspins
1554 CALL auxbas_pw_pool%create_pw(v_rspace(ispin))
1555 CALL pw_zero(v_rspace(ispin))
1560 CALL qs_rho_get(rho, rho_r=rho_r)
1561 IF (
ASSOCIATED(v_tau_rspace))
THEN
1562 CALL qs_rho_get(rho, tau_r=tau_r)
1564 DO ispin = 1, nspins
1566 CALL pw_scale(v_rspace(ispin), v_rspace(ispin)%pw_grid%dvol)
1567 CALL pw_axpy(ec_env%vh_rspace, v_rspace(ispin))
1569 CALL integrate_v_rspace(v_rspace=v_rspace(ispin), &
1570 hmat=ec_env%matrix_ks(ispin, 1), &
1572 calculate_forces=.false., &
1573 basis_type=
"HARRIS", &
1574 task_list_external=ec_env%task_list)
1576 IF (
ASSOCIATED(v_tau_rspace))
THEN
1578 CALL pw_scale(v_tau_rspace(ispin), v_tau_rspace(ispin)%pw_grid%dvol)
1579 CALL integrate_v_rspace(v_rspace=v_tau_rspace(ispin), &
1580 hmat=ec_env%matrix_ks(ispin, 1), &
1582 calculate_forces=.false., &
1583 compute_tau=.true., &
1584 basis_type=
"HARRIS", &
1585 task_list_external=ec_env%task_list)
1589 evhxc = evhxc + pw_integral_ab(rho_r(ispin), v_rspace(ispin))/ &
1590 v_rspace(1)%pw_grid%dvol
1591 IF (
ASSOCIATED(v_tau_rspace))
THEN
1592 evhxc = evhxc + pw_integral_ab(tau_r(ispin), v_tau_rspace(ispin))/ &
1593 v_tau_rspace(ispin)%pw_grid%dvol
1599 DO ispin = 1, nspins
1600 CALL auxbas_pw_pool%give_back_pw(v_rspace(ispin))
1601 IF (
ASSOCIATED(v_tau_rspace))
THEN
1602 CALL auxbas_pw_pool%give_back_pw(v_tau_rspace(ispin))
1605 DEALLOCATE (v_rspace)
1606 IF (
ASSOCIATED(v_tau_rspace))
DEALLOCATE (v_tau_rspace)
1613 DO ispin = 1, nspins
1614 CALL dbcsr_add(ec_env%matrix_ks(ispin, 1)%matrix, ec_env%matrix_h(1, 1)%matrix, &
1615 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
1616 CALL dbcsr_filter(ec_env%matrix_ks(ispin, 1)%matrix, &
1617 dft_control%qs_control%eps_filter_matrix)
1620 CALL timestop(handle)
1622 END SUBROUTINE ec_build_ks_matrix
1634 SUBROUTINE ec_build_core_hamiltonian_force(qs_env, ec_env, matrix_p, matrix_s, matrix_w)
1635 TYPE(qs_environment_type),
POINTER :: qs_env
1636 TYPE(energy_correction_type),
POINTER :: ec_env
1637 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_p, matrix_s, matrix_w
1639 CHARACTER(LEN=*),
PARAMETER :: routinen =
'ec_build_core_hamiltonian_force'
1641 CHARACTER(LEN=default_string_length) :: basis_type
1642 INTEGER :: handle, iounit, nder, nimages
1643 LOGICAL :: calculate_forces, debug_forces, &
1644 debug_stress, use_virial
1645 REAL(kind=dp) :: fconv
1646 REAL(kind=dp),
DIMENSION(3) :: fodeb
1647 REAL(kind=dp),
DIMENSION(3, 3) :: stdeb, sttot
1648 TYPE(cell_type),
POINTER :: cell
1649 TYPE(cp_logger_type),
POINTER :: logger
1650 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: scrm
1651 TYPE(dft_control_type),
POINTER :: dft_control
1652 TYPE(mp_para_env_type),
POINTER :: para_env
1653 TYPE(neighbor_list_set_p_type),
DIMENSION(:), &
1655 TYPE(qs_force_type),
DIMENSION(:),
POINTER :: force
1656 TYPE(qs_ks_env_type),
POINTER :: ks_env
1657 TYPE(virial_type),
POINTER :: virial
1659 CALL timeset(routinen, handle)
1661 debug_forces = ec_env%debug_forces
1662 debug_stress = ec_env%debug_stress
1664 logger => cp_get_default_logger()
1665 IF (logger%para_env%is_source())
THEN
1666 iounit = cp_logger_get_default_unit_nr(logger, local=.true.)
1671 calculate_forces = .true.
1673 basis_type =
"HARRIS"
1676 NULLIFY (cell, dft_control, force, ks_env, para_env, virial)
1677 CALL get_qs_env(qs_env=qs_env, &
1679 dft_control=dft_control, &
1682 para_env=para_env, &
1684 nimages = dft_control%nimages
1685 IF (nimages /= 1)
THEN
1686 cpabort(
"K-points for Harris functional not implemented")
1690 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
1692 fconv = 1.0e-9_dp*pascal/cell%deth
1693 IF (debug_stress .AND. use_virial)
THEN
1694 sttot = virial%pv_virial
1698 IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc)
THEN
1699 IF (ec_env%energy_functional == ec_functional_harris)
THEN
1700 cpabort(
"Harris functional for GAPW not implemented")
1705 sab_orb => ec_env%sab_orb
1709 CALL dbcsr_allocate_matrix_set(scrm, 1, 1)
1710 ALLOCATE (scrm(1, 1)%matrix)
1711 CALL dbcsr_create(scrm(1, 1)%matrix, template=matrix_s(1, 1)%matrix)
1712 CALL cp_dbcsr_alloc_block_from_nbl(scrm(1, 1)%matrix, sab_orb)
1715 IF (
SIZE(matrix_p, 1) == 2)
THEN
1716 CALL dbcsr_add(matrix_w(1, 1)%matrix, matrix_w(2, 1)%matrix, &
1717 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
1721 IF (debug_forces) fodeb(1:3) = force(1)%overlap(1:3, 1)
1722 IF (debug_stress .AND. use_virial) stdeb = virial%pv_overlap
1723 CALL build_overlap_matrix(ks_env, matrixkp_s=scrm, &
1724 matrix_name=
"OVERLAP MATRIX", &
1725 basis_type_a=basis_type, &
1726 basis_type_b=basis_type, &
1727 sab_nl=sab_orb, calculate_forces=.true., &
1728 matrixkp_p=matrix_w)
1730 IF (debug_forces)
THEN
1731 fodeb(1:3) = force(1)%overlap(1:3, 1) - fodeb(1:3)
1732 CALL para_env%sum(fodeb)
1733 IF (iounit > 0)
WRITE (iounit,
"(T3,A,T33,3F16.8)")
"DEBUG:: Wout*dS ", fodeb
1735 IF (debug_stress .AND. use_virial)
THEN
1736 stdeb = fconv*(virial%pv_overlap - stdeb)
1737 CALL para_env%sum(stdeb)
1738 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
1739 'STRESS| Wout*dS', one_third_sum_diag(stdeb), det_3x3(stdeb)
1742 CALL kinetic_energy_matrix(qs_env, matrixkp_t=scrm, matrix_p=matrix_p, &
1743 calculate_forces=.true., &
1745 basis_type=basis_type, &
1746 debug_forces=debug_forces, debug_stress=debug_stress)
1748 CALL core_matrices(qs_env, scrm, matrix_p, calculate_forces, nder, &
1749 ec_env=ec_env, ec_env_matrices=.false., basis_type=basis_type, &
1750 debug_forces=debug_forces, debug_stress=debug_stress)
1753 ec_env%efield_nuclear = 0.0_dp
1754 IF (calculate_forces .AND. debug_forces) fodeb(1:3) = force(1)%efield(1:3, 1)
1755 CALL ec_efield_local_operator(qs_env, ec_env, calculate_forces)
1756 IF (calculate_forces .AND. debug_forces)
THEN
1757 fodeb(1:3) = force(1)%efield(1:3, 1) - fodeb(1:3)
1758 CALL para_env%sum(fodeb)
1759 IF (iounit > 0)
WRITE (iounit,
"(T3,A,T33,3F16.8)")
"DEBUG:: Pout*dEfield", fodeb
1761 IF (debug_stress .AND. use_virial)
THEN
1762 stdeb = fconv*(virial%pv_virial - sttot)
1763 CALL para_env%sum(stdeb)
1764 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
1765 'STRESS| Stress Pout*dHcore ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1766 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))")
' '
1770 CALL dbcsr_deallocate_matrix_set(scrm)
1772 CALL timestop(handle)
1774 END SUBROUTINE ec_build_core_hamiltonian_force
1785 SUBROUTINE ec_build_ks_matrix_force(qs_env, ec_env)
1786 TYPE(qs_environment_type),
POINTER :: qs_env
1787 TYPE(energy_correction_type),
POINTER :: ec_env
1789 CHARACTER(LEN=*),
PARAMETER :: routinen =
'ec_build_ks_matrix_force'
1791 CHARACTER(LEN=default_string_length) :: unit_string
1792 INTEGER :: handle, i, iounit, ispin, natom, nspins
1793 LOGICAL :: debug_forces, debug_stress, do_ec_hfx, &
1795 REAL(dp) :: dehartree, dummy_real, dummy_real2(2), &
1796 eexc, ehartree, eovrl, exc, fconv
1797 REAL(dp),
ALLOCATABLE,
DIMENSION(:, :) :: ftot
1798 REAL(dp),
DIMENSION(3) :: fodeb
1799 REAL(kind=dp),
DIMENSION(3, 3) :: h_stress, pv_loc, stdeb, sttot
1800 TYPE(atomic_kind_type),
DIMENSION(:),
POINTER :: atomic_kind_set
1801 TYPE(cell_type),
POINTER :: cell
1802 TYPE(cp_logger_type),
POINTER :: logger
1803 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_ks, scrm
1804 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_p, matrix_s
1805 TYPE(dft_control_type),
POINTER :: dft_control
1806 TYPE(mp_para_env_type),
POINTER :: para_env
1807 TYPE(neighbor_list_set_p_type),
DIMENSION(:), &
1809 TYPE(pw_c1d_gs_type) :: rho_tot_gspace, rhodn_tot_gspace, &
1811 TYPE(pw_c1d_gs_type),
DIMENSION(:),
POINTER :: rho_g, rhoout_g
1812 TYPE(pw_c1d_gs_type),
POINTER :: rho_core
1813 TYPE(pw_env_type),
POINTER :: pw_env
1814 TYPE(pw_poisson_type),
POINTER :: poisson_env
1815 TYPE(pw_pool_type),
POINTER :: auxbas_pw_pool
1816 TYPE(pw_r3d_rs_type) :: dv_hartree_rspace, v_hartree_rspace, &
1818 TYPE(pw_r3d_rs_type),
DIMENSION(:),
POINTER :: rho_r, rhoout_r, tau_r, tauout_r, &
1819 v_rspace, v_tau_rspace, v_xc, v_xc_tau
1820 TYPE(qs_force_type),
DIMENSION(:),
POINTER :: force
1821 TYPE(qs_ks_env_type),
POINTER :: ks_env
1822 TYPE(qs_rho_type),
POINTER :: rho
1823 TYPE(section_vals_type),
POINTER :: ec_hfx_sections, xc_section
1824 TYPE(virial_type),
POINTER :: virial
1826 CALL timeset(routinen, handle)
1828 debug_forces = ec_env%debug_forces
1829 debug_stress = ec_env%debug_stress
1831 logger => cp_get_default_logger()
1832 IF (logger%para_env%is_source())
THEN
1833 iounit = cp_logger_get_default_unit_nr(logger, local=.true.)
1839 NULLIFY (atomic_kind_set, cell, dft_control, force, ks_env, &
1840 matrix_ks, matrix_p, matrix_s, para_env, rho, rho_core, &
1841 rho_g, rho_r, sab_orb, tau_r, virial)
1842 CALL get_qs_env(qs_env=qs_env, &
1844 dft_control=dft_control, &
1847 matrix_ks=matrix_ks, &
1848 para_env=para_env, &
1853 nspins = dft_control%nspins
1854 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
1858 fconv = cp_unit_from_cp2k(1.0_dp/cell%deth, trim(unit_string))
1860 IF (debug_stress .AND. use_virial)
THEN
1861 sttot = virial%pv_virial
1865 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
1866 cpassert(
ASSOCIATED(pw_env))
1868 NULLIFY (auxbas_pw_pool, poisson_env)
1870 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
1871 poisson_env=poisson_env)
1874 CALL auxbas_pw_pool%create_pw(v_hartree_gspace)
1875 CALL auxbas_pw_pool%create_pw(rhodn_tot_gspace)
1876 CALL auxbas_pw_pool%create_pw(v_hartree_rspace)
1878 CALL pw_transfer(ec_env%vh_rspace, v_hartree_rspace)
1883 CALL qs_rho_get(rho, rho_r=rho_r, rho_g=rho_g, tau_r=tau_r)
1884 NULLIFY (rhoout_r, rhoout_g)
1885 ALLOCATE (rhoout_r(nspins), rhoout_g(nspins))
1886 DO ispin = 1, nspins
1887 CALL auxbas_pw_pool%create_pw(rhoout_r(ispin))
1888 CALL auxbas_pw_pool%create_pw(rhoout_g(ispin))
1890 CALL auxbas_pw_pool%create_pw(dv_hartree_rspace)
1891 CALL auxbas_pw_pool%create_pw(vtot_rspace)
1893 CALL pw_zero(rhodn_tot_gspace)
1894 DO ispin = 1, nspins
1895 CALL calculate_rho_elec(ks_env=ks_env, matrix_p=ec_env%matrix_p(ispin, 1)%matrix, &
1896 rho=rhoout_r(ispin), &
1897 rho_gspace=rhoout_g(ispin), &
1898 basis_type=
"HARRIS", &
1899 task_list_external=ec_env%task_list)
1903 ALLOCATE (ec_env%rhoout_r(nspins))
1904 DO ispin = 1, nspins
1905 CALL auxbas_pw_pool%create_pw(ec_env%rhoout_r(ispin))
1906 CALL pw_copy(rhoout_r(ispin), ec_env%rhoout_r(ispin))
1910 IF (dft_control%use_kinetic_energy_density)
THEN
1912 TYPE(pw_c1d_gs_type) :: tauout_g
1913 ALLOCATE (tauout_r(nspins))
1914 DO ispin = 1, nspins
1915 CALL auxbas_pw_pool%create_pw(tauout_r(ispin))
1917 CALL auxbas_pw_pool%create_pw(tauout_g)
1919 DO ispin = 1, nspins
1920 CALL calculate_rho_elec(ks_env=ks_env, matrix_p=ec_env%matrix_p(ispin, 1)%matrix, &
1921 rho=tauout_r(ispin), &
1922 rho_gspace=tauout_g, &
1923 compute_tau=.true., &
1924 basis_type=
"HARRIS", &
1925 task_list_external=ec_env%task_list)
1928 CALL auxbas_pw_pool%give_back_pw(tauout_g)
1932 IF (use_virial)
THEN
1935 CALL auxbas_pw_pool%create_pw(rho_tot_gspace)
1938 CALL calc_rho_tot_gspace(rho_tot_gspace, qs_env, rho)
1941 CALL get_qs_env(qs_env=qs_env, rho_core=rho_core)
1942 CALL pw_copy(rho_core, rhodn_tot_gspace)
1943 DO ispin = 1, dft_control%nspins
1944 CALL pw_axpy(rhoout_g(ispin), rhodn_tot_gspace)
1948 h_stress(:, :) = 0.0_dp
1949 CALL pw_poisson_solve(poisson_env, &
1950 density=rho_tot_gspace, &
1951 ehartree=ehartree, &
1952 vhartree=v_hartree_gspace, &
1953 h_stress=h_stress, &
1954 aux_density=rhodn_tot_gspace)
1956 virial%pv_ehartree = virial%pv_ehartree + h_stress/real(para_env%num_pe, dp)
1957 virial%pv_virial = virial%pv_virial + h_stress/real(para_env%num_pe, dp)
1959 IF (debug_stress)
THEN
1960 stdeb = fconv*(h_stress/real(para_env%num_pe, dp))
1961 CALL para_env%sum(stdeb)
1962 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
1963 'STRESS| GREEN 1st v_H[n_in]*n_out ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1967 virial%pv_calculate = .true.
1969 NULLIFY (v_rspace, v_tau_rspace)
1970 CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho, xc_section=ec_env%xc_section, &
1971 vxc_rho=v_rspace, vxc_tau=v_tau_rspace, exc=exc, just_energy=.false.)
1974 virial%pv_exc = virial%pv_exc - virial%pv_xc
1975 virial%pv_virial = virial%pv_virial - virial%pv_xc
1977 IF (debug_stress)
THEN
1978 stdeb = -1.0_dp*fconv*virial%pv_xc
1979 CALL para_env%sum(stdeb)
1980 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
1981 'STRESS| GGA 1st E_xc[Pin] ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1984 IF (
ASSOCIATED(v_rspace))
THEN
1985 DO ispin = 1, nspins
1986 CALL auxbas_pw_pool%give_back_pw(v_rspace(ispin))
1988 DEALLOCATE (v_rspace)
1990 IF (
ASSOCIATED(v_tau_rspace))
THEN
1991 DO ispin = 1, nspins
1992 CALL auxbas_pw_pool%give_back_pw(v_tau_rspace(ispin))
1994 DEALLOCATE (v_tau_rspace)
1996 CALL pw_zero(rhodn_tot_gspace)
2001 DO ispin = 1, nspins
2002 CALL pw_axpy(rho_r(ispin), rhoout_r(ispin), -1.0_dp)
2003 CALL pw_axpy(rho_g(ispin), rhoout_g(ispin), -1.0_dp)
2004 CALL pw_axpy(rhoout_g(ispin), rhodn_tot_gspace)
2005 IF (dft_control%use_kinetic_energy_density)
CALL pw_axpy(tau_r(ispin), tauout_r(ispin), -1.0_dp)
2009 IF (use_virial)
THEN
2012 h_stress(:, :) = 0.0_dp
2013 CALL pw_poisson_solve(poisson_env, &
2014 density=rhodn_tot_gspace, &
2015 ehartree=dehartree, &
2016 vhartree=v_hartree_gspace, &
2017 h_stress=h_stress, &
2018 aux_density=rho_tot_gspace)
2020 CALL auxbas_pw_pool%give_back_pw(rho_tot_gspace)
2022 virial%pv_ehartree = virial%pv_ehartree + h_stress/real(para_env%num_pe, dp)
2023 virial%pv_virial = virial%pv_virial + h_stress/real(para_env%num_pe, dp)
2025 IF (debug_stress)
THEN
2026 stdeb = fconv*(h_stress/real(para_env%num_pe, dp))
2027 CALL para_env%sum(stdeb)
2028 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
2029 'STRESS| GREEN 2nd V_H[dP]*n_in ', one_third_sum_diag(stdeb), det_3x3(stdeb)
2034 CALL pw_poisson_solve(poisson_env, rhodn_tot_gspace, dehartree, &
2038 CALL pw_transfer(v_hartree_gspace, dv_hartree_rspace)
2039 CALL pw_scale(dv_hartree_rspace, dv_hartree_rspace%pw_grid%dvol)
2042 CALL pw_transfer(v_hartree_rspace, vtot_rspace)
2043 CALL pw_axpy(dv_hartree_rspace, vtot_rspace)
2044 IF (debug_forces) fodeb(1:3) = force(1)%rho_core(1:3, 1)
2045 IF (debug_stress .AND. use_virial) stdeb = virial%pv_ehartree
2046 CALL integrate_v_core_rspace(vtot_rspace, qs_env)
2047 IF (debug_forces)
THEN
2048 fodeb(1:3) = force(1)%rho_core(1:3, 1) - fodeb(1:3)
2049 CALL para_env%sum(fodeb)
2050 IF (iounit > 0)
WRITE (iounit,
"(T3,A,T33,3F16.8)")
"DEBUG:: Vtot*dncore", fodeb
2052 IF (debug_stress .AND. use_virial)
THEN
2053 stdeb = fconv*(virial%pv_ehartree - stdeb)
2054 CALL para_env%sum(stdeb)
2055 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
2056 'STRESS| Vtot*dncore', one_third_sum_diag(stdeb), det_3x3(stdeb)
2062 xc_section => ec_env%xc_section
2064 IF (use_virial) virial%pv_xc = 0.0_dp
2065 NULLIFY (v_xc, v_xc_tau)
2066 CALL create_kernel(qs_env, &
2073 xc_section=xc_section, &
2074 compute_virial=use_virial, &
2075 virial_xc=virial%pv_xc)
2077 IF (use_virial)
THEN
2079 virial%pv_exc = virial%pv_exc + virial%pv_xc
2080 virial%pv_virial = virial%pv_virial + virial%pv_xc
2082 IF (debug_stress .AND. use_virial)
THEN
2083 stdeb = 1.0_dp*fconv*virial%pv_xc
2084 CALL para_env%sum(stdeb)
2085 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
2086 'STRESS| GGA 2nd f_Hxc[dP]*Pin ', one_third_sum_diag(stdeb), det_3x3(stdeb)
2089 CALL get_qs_env(qs_env=qs_env, rho=rho, matrix_s_kp=matrix_s)
2090 NULLIFY (ec_env%matrix_hz)
2091 CALL dbcsr_allocate_matrix_set(ec_env%matrix_hz, nspins)
2092 DO ispin = 1, nspins
2093 ALLOCATE (ec_env%matrix_hz(ispin)%matrix)
2094 CALL dbcsr_create(ec_env%matrix_hz(ispin)%matrix, template=matrix_s(1, 1)%matrix)
2095 CALL dbcsr_copy(ec_env%matrix_hz(ispin)%matrix, matrix_s(1, 1)%matrix)
2096 CALL dbcsr_set(ec_env%matrix_hz(ispin)%matrix, 0.0_dp)
2098 CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
2100 IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
2101 IF (debug_stress .AND. use_virial) stdeb = virial%pv_virial
2104 IF (use_virial)
THEN
2105 pv_loc = virial%pv_virial
2108 DO ispin = 1, nspins
2109 CALL pw_scale(v_xc(ispin), v_xc(ispin)%pw_grid%dvol)
2110 CALL pw_axpy(dv_hartree_rspace, v_xc(ispin))
2111 CALL integrate_v_rspace(v_rspace=v_xc(ispin), &
2112 hmat=ec_env%matrix_hz(ispin), &
2113 pmat=matrix_p(ispin, 1), &
2115 calculate_forces=.true.)
2118 IF (debug_forces)
THEN
2119 fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
2120 CALL para_env%sum(fodeb)
2121 IF (iounit > 0)
WRITE (iounit,
"(T3,A,T33,3F16.8)")
"DEBUG:: Pin*dKdrho", fodeb
2123 IF (debug_stress .AND. use_virial)
THEN
2124 stdeb = fconv*(virial%pv_virial - stdeb)
2125 CALL para_env%sum(stdeb)
2126 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
2127 'STRESS| INT 2nd f_Hxc[dP]*Pin ', one_third_sum_diag(stdeb), det_3x3(stdeb)
2130 IF (
ASSOCIATED(v_xc_tau))
THEN
2131 IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
2132 IF (debug_stress .AND. use_virial) stdeb = virial%pv_virial
2134 DO ispin = 1, nspins
2135 CALL pw_scale(v_xc_tau(ispin), v_xc_tau(ispin)%pw_grid%dvol)
2136 CALL integrate_v_rspace(v_rspace=v_xc_tau(ispin), &
2137 hmat=ec_env%matrix_hz(ispin), &
2138 pmat=matrix_p(ispin, 1), &
2140 compute_tau=.true., &
2141 calculate_forces=.true.)
2144 IF (debug_forces)
THEN
2145 fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
2146 CALL para_env%sum(fodeb)
2147 IF (iounit > 0)
WRITE (iounit,
"(T3,A,T33,3F16.8)")
"DEBUG:: Pin*dKtaudtau", fodeb
2149 IF (debug_stress .AND. use_virial)
THEN
2150 stdeb = fconv*(virial%pv_virial - stdeb)
2151 CALL para_env%sum(stdeb)
2152 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
2153 'STRESS| INT 2nd f_xctau[dP]*Pin ', one_third_sum_diag(stdeb), det_3x3(stdeb)
2157 IF (use_virial)
THEN
2158 virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
2163 CALL dbcsr_allocate_matrix_set(scrm, nspins)
2164 DO ispin = 1, nspins
2165 ALLOCATE (scrm(ispin)%matrix)
2166 CALL dbcsr_create(scrm(ispin)%matrix, template=ec_env%matrix_ks(ispin, 1)%matrix)
2167 CALL dbcsr_copy(scrm(ispin)%matrix, ec_env%matrix_ks(ispin, 1)%matrix)
2168 CALL dbcsr_set(scrm(ispin)%matrix, 0.0_dp)
2172 NULLIFY (v_rspace, v_tau_rspace)
2174 CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho, xc_section=ec_env%xc_section, &
2175 vxc_rho=v_rspace, vxc_tau=v_tau_rspace, exc=eexc, just_energy=.false.)
2177 IF (use_virial)
THEN
2179 IF (
ASSOCIATED(v_rspace))
THEN
2180 DO ispin = 1, nspins
2182 eexc = eexc + pw_integral_ab(rhoout_r(ispin), v_rspace(ispin))
2185 IF (
ASSOCIATED(v_tau_rspace))
THEN
2186 DO ispin = 1, nspins
2188 eexc = eexc + pw_integral_ab(tauout_r(ispin), v_tau_rspace(ispin))
2193 IF (.NOT.
ASSOCIATED(v_rspace))
THEN
2194 ALLOCATE (v_rspace(nspins))
2195 DO ispin = 1, nspins
2196 CALL auxbas_pw_pool%create_pw(v_rspace(ispin))
2197 CALL pw_zero(v_rspace(ispin))
2203 IF (use_virial)
THEN
2204 pv_loc = virial%pv_virial
2207 IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
2208 IF (debug_stress .AND. use_virial) stdeb = virial%pv_virial
2209 DO ispin = 1, nspins
2211 CALL pw_scale(v_rspace(ispin), v_rspace(ispin)%pw_grid%dvol)
2212 CALL pw_axpy(v_hartree_rspace, v_rspace(ispin))
2214 CALL integrate_v_rspace(v_rspace=v_rspace(ispin), &
2216 pmat=ec_env%matrix_p(ispin, 1), &
2218 calculate_forces=.true., &
2219 basis_type=
"HARRIS", &
2220 task_list_external=ec_env%task_list)
2223 IF (debug_forces)
THEN
2224 fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
2225 CALL para_env%sum(fodeb)
2226 IF (iounit > 0)
WRITE (iounit,
"(T3,A,T33,3F16.8)")
"DEBUG:: Pout*dVhxc ", fodeb
2228 IF (debug_stress .AND. use_virial)
THEN
2229 stdeb = fconv*(virial%pv_virial - stdeb)
2230 CALL para_env%sum(stdeb)
2231 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
2232 'STRESS| INT Pout*dVhxc ', one_third_sum_diag(stdeb), det_3x3(stdeb)
2236 IF (use_virial)
THEN
2237 virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
2240 IF (
ASSOCIATED(v_tau_rspace))
THEN
2241 IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
2242 DO ispin = 1, nspins
2244 CALL pw_scale(v_tau_rspace(ispin), v_tau_rspace(ispin)%pw_grid%dvol)
2245 CALL integrate_v_rspace(v_rspace=v_tau_rspace(ispin), &
2247 pmat=ec_env%matrix_p(ispin, 1), &
2249 calculate_forces=.true., &
2250 compute_tau=.true., &
2251 basis_type=
"HARRIS", &
2252 task_list_external=ec_env%task_list)
2254 IF (debug_forces)
THEN
2255 fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
2256 CALL para_env%sum(fodeb)
2257 IF (iounit > 0)
WRITE (iounit,
"(T3,A,T33,3F16.8)")
"DEBUG:: Pout*dVhxc_tau ", fodeb
2266 ec_hfx_sections => section_vals_get_subs_vals(qs_env%input,
"DFT%ENERGY_CORRECTION%XC%HF")
2267 CALL section_vals_get(ec_hfx_sections, explicit=do_ec_hfx)
2271 IF (debug_forces) fodeb(1:3) = force(1)%fock_4c(1:3, 1)
2272 IF (use_virial) virial%pv_fock_4c = 0.0_dp
2274 CALL calculate_exx(qs_env=qs_env, &
2276 hfx_sections=ec_hfx_sections, &
2277 x_data=ec_env%x_data, &
2279 do_admm=ec_env%do_ec_admm, &
2280 calc_forces=.true., &
2281 reuse_hfx=ec_env%reuse_hfx, &
2282 do_im_time=.false., &
2283 e_ex_from_gw=dummy_real, &
2284 e_admm_from_gw=dummy_real2, &
2287 IF (use_virial)
THEN
2288 virial%pv_exx = virial%pv_exx - virial%pv_fock_4c
2289 virial%pv_virial = virial%pv_virial - virial%pv_fock_4c
2290 virial%pv_calculate = .false.
2292 IF (debug_forces)
THEN
2293 fodeb(1:3) = force(1)%fock_4c(1:3, 1) - fodeb(1:3)
2294 CALL para_env%sum(fodeb)
2295 IF (iounit > 0)
WRITE (iounit,
"(T3,A,T33,3F16.8)")
"DEBUG:: Pout*hfx ", fodeb
2297 IF (debug_stress .AND. use_virial)
THEN
2298 stdeb = -1.0_dp*fconv*virial%pv_fock_4c
2299 CALL para_env%sum(stdeb)
2300 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
2301 'STRESS| Pout*hfx ', one_third_sum_diag(stdeb), det_3x3(stdeb)
2309 CALL dbcsr_deallocate_matrix_set(scrm)
2312 CALL auxbas_pw_pool%give_back_pw(v_hartree_rspace)
2313 DO ispin = 1, nspins
2314 CALL auxbas_pw_pool%give_back_pw(v_rspace(ispin))
2315 IF (
ASSOCIATED(v_tau_rspace))
THEN
2316 CALL auxbas_pw_pool%give_back_pw(v_tau_rspace(ispin))
2319 IF (
ASSOCIATED(v_tau_rspace))
DEALLOCATE (v_tau_rspace)
2322 IF (debug_forces) fodeb(1:3) = force(1)%core_overlap(1:3, 1)
2323 IF (debug_stress .AND. use_virial) stdeb = virial%pv_ecore_overlap
2324 CALL calculate_ecore_overlap(qs_env, para_env, .true., e_overlap_core=eovrl)
2325 IF (debug_forces)
THEN
2326 fodeb(1:3) = force(1)%core_overlap(1:3, 1) - fodeb(1:3)
2327 CALL para_env%sum(fodeb)
2328 IF (iounit > 0)
WRITE (iounit,
"(T3,A,T33,3F16.8)")
"DEBUG:: CoreOverlap", fodeb
2330 IF (debug_stress .AND. use_virial)
THEN
2331 stdeb = fconv*(stdeb - virial%pv_ecore_overlap)
2332 CALL para_env%sum(stdeb)
2333 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
2334 'STRESS| CoreOverlap ', one_third_sum_diag(stdeb), det_3x3(stdeb)
2337 IF (debug_forces)
THEN
2338 CALL get_qs_env(qs_env, natom=natom, atomic_kind_set=atomic_kind_set)
2339 ALLOCATE (ftot(3, natom))
2340 CALL total_qs_force(ftot, force, atomic_kind_set)
2341 fodeb(1:3) = ftot(1:3, 1)
2343 CALL para_env%sum(fodeb)
2344 IF (iounit > 0)
WRITE (iounit,
"(T3,A,T33,3F16.8)")
"DEBUG:: Force Explicit", fodeb
2347 DEALLOCATE (v_rspace)
2349 CALL auxbas_pw_pool%give_back_pw(dv_hartree_rspace)
2350 CALL auxbas_pw_pool%give_back_pw(vtot_rspace)
2351 DO ispin = 1, nspins
2352 CALL auxbas_pw_pool%give_back_pw(rhoout_r(ispin))
2353 CALL auxbas_pw_pool%give_back_pw(rhoout_g(ispin))
2354 CALL auxbas_pw_pool%give_back_pw(v_xc(ispin))
2356 DEALLOCATE (rhoout_r, rhoout_g, v_xc)
2357 IF (
ASSOCIATED(tauout_r))
THEN
2358 DO ispin = 1, nspins
2359 CALL auxbas_pw_pool%give_back_pw(tauout_r(ispin))
2361 DEALLOCATE (tauout_r)
2363 IF (
ASSOCIATED(v_xc_tau))
THEN
2364 DO ispin = 1, nspins
2365 CALL auxbas_pw_pool%give_back_pw(v_xc_tau(ispin))
2367 DEALLOCATE (v_xc_tau)
2369 CALL auxbas_pw_pool%give_back_pw(v_hartree_gspace)
2370 CALL auxbas_pw_pool%give_back_pw(rhodn_tot_gspace)
2374 IF (use_virial)
THEN
2375 IF (qs_env%energy_correction)
THEN
2376 ec_env%ehartree = ehartree + dehartree
2377 ec_env%exc = exc + eexc
2381 IF (debug_stress .AND. use_virial)
THEN
2383 stdeb = -1.0_dp*fconv*ehartree
2384 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
2385 'STRESS| VOL 1st v_H[n_in]*n_out', one_third_sum_diag(stdeb), det_3x3(stdeb)
2387 stdeb = -1.0_dp*fconv*exc
2388 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
2389 'STRESS| VOL 1st E_XC[n_in]', one_third_sum_diag(stdeb), det_3x3(stdeb)
2391 stdeb = -1.0_dp*fconv*dehartree
2392 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
2393 'STRESS| VOL 2nd v_H[dP]*n_in', one_third_sum_diag(stdeb), det_3x3(stdeb)
2395 stdeb = -1.0_dp*fconv*eexc
2396 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
2397 'STRESS| VOL 2nd v_XC[n_in]*dP', one_third_sum_diag(stdeb), det_3x3(stdeb)
2402 TYPE(virial_type) :: virdeb
2405 CALL para_env%sum(virdeb%pv_overlap)
2406 CALL para_env%sum(virdeb%pv_ekinetic)
2407 CALL para_env%sum(virdeb%pv_ppl)
2408 CALL para_env%sum(virdeb%pv_ppnl)
2409 CALL para_env%sum(virdeb%pv_ecore_overlap)
2410 CALL para_env%sum(virdeb%pv_ehartree)
2411 CALL para_env%sum(virdeb%pv_exc)
2412 CALL para_env%sum(virdeb%pv_exx)
2413 CALL para_env%sum(virdeb%pv_vdw)
2414 CALL para_env%sum(virdeb%pv_mp2)
2415 CALL para_env%sum(virdeb%pv_nlcc)
2416 CALL para_env%sum(virdeb%pv_gapw)
2417 CALL para_env%sum(virdeb%pv_lrigpw)
2418 CALL para_env%sum(virdeb%pv_virial)
2419 CALL symmetrize_virial(virdeb)
2423 virdeb%pv_ehartree(i, i) = virdeb%pv_ehartree(i, i) - 2.0_dp*(ehartree + dehartree)
2424 virdeb%pv_virial(i, i) = virdeb%pv_virial(i, i) - exc - eexc &
2425 - 2.0_dp*(ehartree + dehartree)
2426 virdeb%pv_exc(i, i) = virdeb%pv_exc(i, i) - exc - eexc
2432 CALL para_env%sum(sttot)
2433 stdeb = fconv*(virdeb%pv_virial - sttot)
2434 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
2435 'STRESS| Explicit electronic stress ', one_third_sum_diag(stdeb), det_3x3(stdeb)
2437 stdeb = fconv*(virdeb%pv_virial)
2438 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
2439 'STRESS| Explicit total stress ', one_third_sum_diag(stdeb), det_3x3(stdeb)
2441 CALL write_stress_tensor_components(virdeb, iounit, cell, unit_string)
2442 CALL write_stress_tensor(virdeb%pv_virial, iounit, cell, unit_string, .false.)
2447 CALL timestop(handle)
2449 END SUBROUTINE ec_build_ks_matrix_force
2459 SUBROUTINE ec_ks_solver(qs_env, ec_env)
2461 TYPE(qs_environment_type),
POINTER :: qs_env
2462 TYPE(energy_correction_type),
POINTER :: ec_env
2464 CHARACTER(LEN=*),
PARAMETER :: routinen =
'ec_ks_solver'
2466 CHARACTER(LEN=default_string_length) :: headline
2467 INTEGER :: handle, ispin, nspins
2468 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: ksmat, pmat, smat, wmat
2469 TYPE(dft_control_type),
POINTER :: dft_control
2471 CALL timeset(routinen, handle)
2473 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
2474 nspins = dft_control%nspins
2477 IF (.NOT.
ASSOCIATED(ec_env%matrix_p))
THEN
2478 headline =
"DENSITY MATRIX"
2479 CALL dbcsr_allocate_matrix_set(ec_env%matrix_p, nspins, 1)
2480 DO ispin = 1, nspins
2481 ALLOCATE (ec_env%matrix_p(ispin, 1)%matrix)
2482 CALL dbcsr_create(ec_env%matrix_p(ispin, 1)%matrix, name=trim(headline), &
2483 template=ec_env%matrix_s(1, 1)%matrix)
2484 CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_p(ispin, 1)%matrix, ec_env%sab_orb)
2488 IF (.NOT.
ASSOCIATED(ec_env%matrix_w))
THEN
2489 headline =
"ENERGY WEIGHTED DENSITY MATRIX"
2490 CALL dbcsr_allocate_matrix_set(ec_env%matrix_w, nspins, 1)
2491 DO ispin = 1, nspins
2492 ALLOCATE (ec_env%matrix_w(ispin, 1)%matrix)
2493 CALL dbcsr_create(ec_env%matrix_w(ispin, 1)%matrix, name=trim(headline), &
2494 template=ec_env%matrix_s(1, 1)%matrix)
2495 CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_w(ispin, 1)%matrix, ec_env%sab_orb)
2499 IF (ec_env%mao)
THEN
2500 CALL mao_create_matrices(ec_env, ksmat, smat, pmat, wmat)
2502 ksmat => ec_env%matrix_ks
2503 smat => ec_env%matrix_s
2504 pmat => ec_env%matrix_p
2505 wmat => ec_env%matrix_w
2508 SELECT CASE (ec_env%ks_solver)
2509 CASE (ec_diagonalization)
2510 CALL ec_diag_solver(qs_env, ksmat, smat, pmat, wmat)
2512 CALL ec_ot_diag_solver(qs_env, ec_env, ksmat, smat, pmat, wmat)
2513 CASE (ec_matrix_sign, ec_matrix_trs4, ec_matrix_tc2)
2514 CALL ec_ls_init(qs_env, ksmat, smat)
2515 CALL ec_ls_solver(qs_env, pmat, wmat, ec_ls_method=ec_env%ks_solver)
2520 IF (ec_env%mao)
THEN
2521 CALL mao_release_matrices(ec_env, ksmat, smat, pmat, wmat)
2524 CALL timestop(handle)
2526 END SUBROUTINE ec_ks_solver
2539 SUBROUTINE mao_create_matrices(ec_env, ksmat, smat, pmat, wmat)
2541 TYPE(energy_correction_type),
POINTER :: ec_env
2542 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: ksmat, smat, pmat, wmat
2544 CHARACTER(LEN=*),
PARAMETER :: routinen =
'mao_create_matrices'
2546 INTEGER :: handle, ispin, nspins
2547 INTEGER,
DIMENSION(:),
POINTER :: col_blk_sizes
2548 TYPE(dbcsr_distribution_type) :: dbcsr_dist
2549 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: mao_coef
2550 TYPE(dbcsr_type) :: cgmat
2552 CALL timeset(routinen, handle)
2554 mao_coef => ec_env%mao_coef
2556 NULLIFY (ksmat, smat, pmat, wmat)
2557 nspins =
SIZE(ec_env%matrix_ks, 1)
2558 CALL dbcsr_get_info(mao_coef(1)%matrix, col_blk_size=col_blk_sizes, distribution=dbcsr_dist)
2559 CALL dbcsr_allocate_matrix_set(ksmat, nspins, 1)
2560 CALL dbcsr_allocate_matrix_set(smat, nspins, 1)
2561 DO ispin = 1, nspins
2562 ALLOCATE (ksmat(ispin, 1)%matrix)
2563 CALL dbcsr_create(ksmat(ispin, 1)%matrix, dist=dbcsr_dist, name=
"MAO KS mat", &
2564 matrix_type=dbcsr_type_symmetric, row_blk_size=col_blk_sizes, &
2565 col_blk_size=col_blk_sizes)
2566 ALLOCATE (smat(ispin, 1)%matrix)
2567 CALL dbcsr_create(smat(ispin, 1)%matrix, dist=dbcsr_dist, name=
"MAO S mat", &
2568 matrix_type=dbcsr_type_symmetric, row_blk_size=col_blk_sizes, &
2569 col_blk_size=col_blk_sizes)
2572 CALL dbcsr_create(cgmat, name=
"TEMP matrix", template=mao_coef(1)%matrix)
2573 DO ispin = 1, nspins
2574 CALL dbcsr_multiply(
"N",
"N", 1.0_dp, ec_env%matrix_s(1, 1)%matrix, mao_coef(ispin)%matrix, &
2576 CALL dbcsr_multiply(
"T",
"N", 1.0_dp, mao_coef(ispin)%matrix, cgmat, 0.0_dp, smat(ispin, 1)%matrix)
2577 CALL dbcsr_multiply(
"N",
"N", 1.0_dp, ec_env%matrix_ks(1, 1)%matrix, mao_coef(ispin)%matrix, &
2579 CALL dbcsr_multiply(
"T",
"N", 1.0_dp, mao_coef(ispin)%matrix, cgmat, 0.0_dp, ksmat(ispin, 1)%matrix)
2581 CALL dbcsr_release(cgmat)
2583 CALL dbcsr_allocate_matrix_set(pmat, nspins, 1)
2584 DO ispin = 1, nspins
2585 ALLOCATE (pmat(ispin, 1)%matrix)
2586 CALL dbcsr_create(pmat(ispin, 1)%matrix, template=smat(1, 1)%matrix, name=
"MAO P mat")
2587 CALL cp_dbcsr_alloc_block_from_nbl(pmat(ispin, 1)%matrix, ec_env%sab_orb)
2590 CALL dbcsr_allocate_matrix_set(wmat, nspins, 1)
2591 DO ispin = 1, nspins
2592 ALLOCATE (wmat(ispin, 1)%matrix)
2593 CALL dbcsr_create(wmat(ispin, 1)%matrix, template=smat(1, 1)%matrix, name=
"MAO W mat")
2594 CALL cp_dbcsr_alloc_block_from_nbl(wmat(ispin, 1)%matrix, ec_env%sab_orb)
2597 CALL timestop(handle)
2599 END SUBROUTINE mao_create_matrices
2612 SUBROUTINE mao_release_matrices(ec_env, ksmat, smat, pmat, wmat)
2614 TYPE(energy_correction_type),
POINTER :: ec_env
2615 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: ksmat, smat, pmat, wmat
2617 CHARACTER(LEN=*),
PARAMETER :: routinen =
'mao_release_matrices'
2619 INTEGER :: handle, ispin, nspins
2620 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: mao_coef
2621 TYPE(dbcsr_type) :: cgmat
2623 CALL timeset(routinen, handle)
2625 mao_coef => ec_env%mao_coef
2626 nspins =
SIZE(mao_coef, 1)
2629 CALL dbcsr_create(cgmat, name=
"TEMP matrix", template=mao_coef(1)%matrix)
2630 DO ispin = 1, nspins
2631 CALL dbcsr_multiply(
"N",
"N", 1.0_dp, mao_coef(ispin)%matrix, pmat(ispin, 1)%matrix, 0.0_dp, cgmat)
2632 CALL dbcsr_multiply(
"N",
"T", 1.0_dp, mao_coef(ispin)%matrix, cgmat, 0.0_dp, &
2633 ec_env%matrix_p(ispin, 1)%matrix, retain_sparsity=.true.)
2634 CALL dbcsr_multiply(
"N",
"N", 1.0_dp, mao_coef(ispin)%matrix, wmat(ispin, 1)%matrix, 0.0_dp, cgmat)
2635 CALL dbcsr_multiply(
"N",
"T", 1.0_dp, mao_coef(ispin)%matrix, cgmat, 0.0_dp, &
2636 ec_env%matrix_w(ispin, 1)%matrix, retain_sparsity=.true.)
2638 CALL dbcsr_release(cgmat)
2640 CALL dbcsr_deallocate_matrix_set(ksmat)
2641 CALL dbcsr_deallocate_matrix_set(smat)
2642 CALL dbcsr_deallocate_matrix_set(pmat)
2643 CALL dbcsr_deallocate_matrix_set(wmat)
2645 CALL timestop(handle)
2647 END SUBROUTINE mao_release_matrices
2660 SUBROUTINE ec_diag_solver(qs_env, matrix_ks, matrix_s, matrix_p, matrix_w)
2662 TYPE(qs_environment_type),
POINTER :: qs_env
2663 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_ks, matrix_s, matrix_p, matrix_w
2665 CHARACTER(LEN=*),
PARAMETER :: routinen =
'ec_diag_solver'
2667 INTEGER :: handle, ispin, nmo(2), nsize, nspins
2668 REAL(kind=dp) :: eps_filter, focc(2)
2669 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:) :: eigenvalues
2670 TYPE(cp_blacs_env_type),
POINTER :: blacs_env
2671 TYPE(cp_fm_struct_type),
POINTER :: fm_struct
2672 TYPE(cp_fm_type) :: fm_ks, fm_mo, fm_ortho
2673 TYPE(dbcsr_type),
POINTER :: buf1_dbcsr, buf2_dbcsr, buf3_dbcsr, &
2674 ortho_dbcsr, ref_matrix
2675 TYPE(dft_control_type),
POINTER :: dft_control
2676 TYPE(mp_para_env_type),
POINTER :: para_env
2678 CALL timeset(routinen, handle)
2680 NULLIFY (blacs_env, para_env)
2681 CALL get_qs_env(qs_env=qs_env, blacs_env=blacs_env, para_env=para_env)
2683 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
2684 eps_filter = dft_control%qs_control%eps_filter_matrix
2685 nspins = dft_control%nspins
2688 CALL get_qs_env(qs_env=qs_env, nelectron_spin=nmo)
2690 IF (nspins == 1)
THEN
2695 CALL dbcsr_get_info(matrix_ks(1, 1)%matrix, nfullrows_total=nsize)
2696 ALLOCATE (eigenvalues(nsize))
2698 NULLIFY (fm_struct, ref_matrix)
2699 CALL cp_fm_struct_create(fm_struct, context=blacs_env, nrow_global=nsize, &
2700 ncol_global=nsize, para_env=para_env)
2701 CALL cp_fm_create(fm_ortho, fm_struct)
2702 CALL cp_fm_create(fm_ks, fm_struct)
2703 CALL cp_fm_create(fm_mo, fm_struct)
2704 CALL cp_fm_struct_release(fm_struct)
2707 ref_matrix => matrix_s(1, 1)%matrix
2708 NULLIFY (ortho_dbcsr, buf1_dbcsr, buf2_dbcsr, buf3_dbcsr)
2709 CALL dbcsr_init_p(ortho_dbcsr)
2710 CALL dbcsr_create(ortho_dbcsr, template=ref_matrix, &
2711 matrix_type=dbcsr_type_no_symmetry)
2712 CALL dbcsr_init_p(buf1_dbcsr)
2713 CALL dbcsr_create(buf1_dbcsr, template=ref_matrix, &
2714 matrix_type=dbcsr_type_no_symmetry)
2715 CALL dbcsr_init_p(buf2_dbcsr)
2716 CALL dbcsr_create(buf2_dbcsr, template=ref_matrix, &
2717 matrix_type=dbcsr_type_no_symmetry)
2718 CALL dbcsr_init_p(buf3_dbcsr)
2719 CALL dbcsr_create(buf3_dbcsr, template=ref_matrix, &
2720 matrix_type=dbcsr_type_no_symmetry)
2722 ref_matrix => matrix_s(1, 1)%matrix
2723 CALL copy_dbcsr_to_fm(ref_matrix, fm_ortho)
2724 CALL cp_fm_cholesky_decompose(fm_ortho)
2725 CALL cp_fm_triangular_invert(fm_ortho)
2726 CALL cp_fm_set_all(fm_ks, 0.0_dp)
2727 CALL cp_fm_to_fm_triangular(fm_ortho, fm_ks,
"U")
2728 CALL copy_fm_to_dbcsr(fm_ks, ortho_dbcsr)
2729 DO ispin = 1, nspins
2731 CALL dbcsr_desymmetrize(matrix_ks(ispin, 1)%matrix, buf1_dbcsr)
2732 CALL dbcsr_multiply(
"N",
"N", 1.0_dp, buf1_dbcsr, ortho_dbcsr, &
2733 0.0_dp, buf2_dbcsr, filter_eps=eps_filter)
2734 CALL dbcsr_multiply(
"T",
"N", 1.0_dp, ortho_dbcsr, buf2_dbcsr, &
2735 0.0_dp, buf1_dbcsr, filter_eps=eps_filter)
2737 CALL copy_dbcsr_to_fm(buf1_dbcsr, fm_ks)
2738 CALL choose_eigv_solver(fm_ks, fm_mo, eigenvalues)
2740 CALL copy_fm_to_dbcsr(fm_mo, buf1_dbcsr)
2741 CALL dbcsr_multiply(
"N",
"N", 1.0_dp, ortho_dbcsr, buf1_dbcsr, &
2742 0.0_dp, buf2_dbcsr, filter_eps=eps_filter)
2744 CALL dbcsr_set(matrix_p(ispin, 1)%matrix, 0.0_dp)
2745 CALL dbcsr_multiply(
"N",
"T", focc(ispin), buf2_dbcsr, buf2_dbcsr, &
2746 1.0_dp, matrix_p(ispin, 1)%matrix, &
2747 retain_sparsity=.true., last_k=nmo(ispin))
2750 CALL dbcsr_set(matrix_w(ispin, 1)%matrix, 0.0_dp)
2751 CALL cp_fm_column_scale(fm_mo, eigenvalues)
2752 CALL copy_fm_to_dbcsr(fm_mo, buf1_dbcsr)
2753 CALL dbcsr_multiply(
"N",
"N", 1.0_dp, ortho_dbcsr, buf1_dbcsr, &
2754 0.0_dp, buf3_dbcsr, filter_eps=eps_filter)
2755 CALL dbcsr_multiply(
"N",
"T", focc(ispin), buf2_dbcsr, buf3_dbcsr, &
2756 1.0_dp, matrix_w(ispin, 1)%matrix, &
2757 retain_sparsity=.true., last_k=nmo(ispin))
2760 CALL cp_fm_release(fm_ks)
2761 CALL cp_fm_release(fm_mo)
2762 CALL cp_fm_release(fm_ortho)
2763 CALL dbcsr_release(ortho_dbcsr)
2764 CALL dbcsr_release(buf1_dbcsr)
2765 CALL dbcsr_release(buf2_dbcsr)
2766 CALL dbcsr_release(buf3_dbcsr)
2767 DEALLOCATE (ortho_dbcsr, buf1_dbcsr, buf2_dbcsr, buf3_dbcsr)
2768 DEALLOCATE (eigenvalues)
2770 CALL timestop(handle)
2772 END SUBROUTINE ec_diag_solver
2780 SUBROUTINE ec_energy(ec_env, unit_nr)
2781 TYPE(energy_correction_type) :: ec_env
2782 INTEGER,
INTENT(IN) :: unit_nr
2784 CHARACTER(LEN=*),
PARAMETER :: routinen =
'ec_energy'
2786 INTEGER :: handle, ispin, nspins
2787 REAL(kind=dp) :: eband, trace
2789 CALL timeset(routinen, handle)
2791 nspins =
SIZE(ec_env%matrix_p, 1)
2792 DO ispin = 1, nspins
2793 CALL dbcsr_dot(ec_env%matrix_p(ispin, 1)%matrix, ec_env%matrix_s(1, 1)%matrix, trace)
2794 IF (unit_nr > 0)
WRITE (unit_nr,
'(T3,A,T65,F16.10)')
'Tr[PS] ', trace
2798 SELECT CASE (ec_env%energy_functional)
2799 CASE (ec_functional_harris)
2803 DO ispin = 1, nspins
2804 CALL dbcsr_dot(ec_env%matrix_ks(ispin, 1)%matrix, ec_env%matrix_p(ispin, 1)%matrix, trace)
2805 eband = eband + trace
2807 ec_env%eband = eband + ec_env%efield_nuclear
2810 ec_env%etotal = ec_env%eband + ec_env%ehartree + ec_env%exc - ec_env%vhxc + ec_env%edispersion - ec_env%ex
2811 IF (unit_nr > 0)
THEN
2812 WRITE (unit_nr,
'(T3,A,T56,F25.15)')
"Eband ", ec_env%eband
2813 WRITE (unit_nr,
'(T3,A,T56,F25.15)')
"Ehartree ", ec_env%ehartree
2814 WRITE (unit_nr,
'(T3,A,T56,F25.15)')
"Exc ", ec_env%exc
2815 WRITE (unit_nr,
'(T3,A,T56,F25.15)')
"Ex ", ec_env%ex
2816 WRITE (unit_nr,
'(T3,A,T56,F25.15)')
"Evhxc ", ec_env%vhxc
2817 WRITE (unit_nr,
'(T3,A,T56,F25.15)')
"Edisp ", ec_env%edispersion
2818 WRITE (unit_nr,
'(T3,A,T56,F25.15)')
"Etotal Harris Functional ", ec_env%etotal
2821 CASE (ec_functional_dc)
2824 CALL calculate_ptrace(ec_env%matrix_h, ec_env%matrix_p, ec_env%ecore,
SIZE(ec_env%matrix_p, 1))
2826 ec_env%ecore = ec_env%ecore + ec_env%efield_nuclear
2827 ec_env%etotal = ec_env%ecore + ec_env%ehartree + ec_env%exc + ec_env%edispersion &
2828 + ec_env%ex + ec_env%exc_aux_fit
2830 IF (unit_nr > 0)
THEN
2831 WRITE (unit_nr,
'(T3,A,T56,F25.15)')
"Ecore ", ec_env%ecore
2832 WRITE (unit_nr,
'(T3,A,T56,F25.15)')
"Ehartree ", ec_env%ehartree
2833 WRITE (unit_nr,
'(T3,A,T56,F25.15)')
"Exc ", ec_env%exc
2834 WRITE (unit_nr,
'(T3,A,T56,F25.15)')
"Ex ", ec_env%ex
2835 WRITE (unit_nr,
'(T3,A,T56,F25.15)')
"Exc_aux_fit", ec_env%exc_aux_fit
2836 WRITE (unit_nr,
'(T3,A,T56,F25.15)')
"Edisp ", ec_env%edispersion
2837 WRITE (unit_nr,
'(T3,A,T56,F25.15)')
"Etotal Energy Functional ", ec_env%etotal
2840 CASE (ec_functional_ext)
2842 ec_env%etotal = ec_env%ex
2843 IF (unit_nr > 0)
THEN
2844 WRITE (unit_nr,
'(T3,A,T56,F25.15)')
"Etotal Energy Functional ", ec_env%etotal
2853 CALL timestop(handle)
2855 END SUBROUTINE ec_energy
2867 SUBROUTINE ec_build_neighborlist(qs_env, ec_env)
2868 TYPE(qs_environment_type),
POINTER :: qs_env
2869 TYPE(energy_correction_type),
POINTER :: ec_env
2871 CHARACTER(LEN=*),
PARAMETER :: routinen =
'ec_build_neighborlist'
2873 INTEGER :: handle, ikind, nkind, zat
2874 LOGICAL :: all_potential_present, &
2875 gth_potential_present, &
2876 sgp_potential_present, &
2877 skip_load_balance_distributed
2878 LOGICAL,
ALLOCATABLE,
DIMENSION(:) :: all_present, default_present, &
2879 orb_present, ppl_present, ppnl_present
2880 REAL(dp) :: subcells
2881 REAL(dp),
ALLOCATABLE,
DIMENSION(:) :: all_radius, c_radius, orb_radius, &
2882 ppl_radius, ppnl_radius
2883 REAL(dp),
ALLOCATABLE,
DIMENSION(:, :) :: pair_radius
2884 TYPE(all_potential_type),
POINTER :: all_potential
2885 TYPE(atomic_kind_type),
DIMENSION(:),
POINTER :: atomic_kind_set
2886 TYPE(cell_type),
POINTER :: cell
2887 TYPE(dft_control_type),
POINTER :: dft_control
2888 TYPE(distribution_1d_type),
POINTER :: distribution_1d
2889 TYPE(distribution_2d_type),
POINTER :: distribution_2d
2890 TYPE(gth_potential_type),
POINTER :: gth_potential
2891 TYPE(gto_basis_set_type),
POINTER :: basis_set
2892 TYPE(local_atoms_type),
ALLOCATABLE,
DIMENSION(:) :: atom2d
2893 TYPE(molecule_type),
DIMENSION(:),
POINTER :: molecule_set
2894 TYPE(neighbor_list_set_p_type),
DIMENSION(:), &
2895 POINTER :: sab_cn, sab_vdw
2896 TYPE(particle_type),
DIMENSION(:),
POINTER :: particle_set
2897 TYPE(qs_dispersion_type),
POINTER :: dispersion_env
2898 TYPE(qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
2899 TYPE(qs_kind_type),
POINTER :: qs_kind
2900 TYPE(qs_ks_env_type),
POINTER :: ks_env
2901 TYPE(sgp_potential_type),
POINTER :: sgp_potential
2903 CALL timeset(routinen, handle)
2905 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set)
2906 CALL get_qs_kind_set(qs_kind_set, all_potential_present=all_potential_present, &
2907 gth_potential_present=gth_potential_present, &
2908 sgp_potential_present=sgp_potential_present)
2909 nkind =
SIZE(qs_kind_set)
2910 ALLOCATE (c_radius(nkind), default_present(nkind))
2911 ALLOCATE (orb_radius(nkind), all_radius(nkind), ppl_radius(nkind), ppnl_radius(nkind))
2912 ALLOCATE (orb_present(nkind), all_present(nkind), ppl_present(nkind), ppnl_present(nkind))
2913 ALLOCATE (pair_radius(nkind, nkind))
2914 ALLOCATE (atom2d(nkind))
2916 CALL get_qs_env(qs_env, &
2917 atomic_kind_set=atomic_kind_set, &
2919 distribution_2d=distribution_2d, &
2920 local_particles=distribution_1d, &
2921 particle_set=particle_set, &
2922 molecule_set=molecule_set)
2924 CALL atom2d_build(atom2d, distribution_1d, distribution_2d, atomic_kind_set, &
2925 molecule_set, .false., particle_set)
2928 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom2d(ikind)%list)
2929 qs_kind => qs_kind_set(ikind)
2930 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type=
"HARRIS")
2931 IF (
ASSOCIATED(basis_set))
THEN
2932 orb_present(ikind) = .true.
2933 CALL get_gto_basis_set(gto_basis_set=basis_set, kind_radius=orb_radius(ikind))
2935 orb_present(ikind) = .false.
2936 orb_radius(ikind) = 0.0_dp
2938 CALL get_qs_kind(qs_kind, all_potential=all_potential, &
2939 gth_potential=gth_potential, sgp_potential=sgp_potential)
2940 IF (gth_potential_present .OR. sgp_potential_present)
THEN
2941 IF (
ASSOCIATED(gth_potential))
THEN
2942 CALL get_potential(potential=gth_potential, &
2943 ppl_present=ppl_present(ikind), &
2944 ppl_radius=ppl_radius(ikind), &
2945 ppnl_present=ppnl_present(ikind), &
2946 ppnl_radius=ppnl_radius(ikind))
2947 ELSE IF (
ASSOCIATED(sgp_potential))
THEN
2948 CALL get_potential(potential=sgp_potential, &
2949 ppl_present=ppl_present(ikind), &
2950 ppl_radius=ppl_radius(ikind), &
2951 ppnl_present=ppnl_present(ikind), &
2952 ppnl_radius=ppnl_radius(ikind))
2954 ppl_present(ikind) = .false.
2955 ppl_radius(ikind) = 0.0_dp
2956 ppnl_present(ikind) = .false.
2957 ppnl_radius(ikind) = 0.0_dp
2961 IF (all_potential_present .OR. sgp_potential_present)
THEN
2962 all_present(ikind) = .false.
2963 all_radius(ikind) = 0.0_dp
2964 IF (
ASSOCIATED(all_potential))
THEN
2965 all_present(ikind) = .true.
2966 CALL get_potential(potential=all_potential, core_charge_radius=all_radius(ikind))
2967 ELSE IF (
ASSOCIATED(sgp_potential))
THEN
2968 IF (sgp_potential%ecp_local)
THEN
2969 all_present(ikind) = .true.
2970 CALL get_potential(potential=sgp_potential, core_charge_radius=all_radius(ikind))
2976 CALL section_vals_val_get(qs_env%input,
"DFT%SUBCELLS", r_val=subcells)
2979 CALL pair_radius_setup(orb_present, orb_present, orb_radius, orb_radius, pair_radius)
2980 CALL build_neighbor_lists(ec_env%sab_orb, particle_set, atom2d, cell, pair_radius, &
2981 subcells=subcells, nlname=
"sab_orb")
2983 IF (all_potential_present .OR. sgp_potential_present)
THEN
2984 IF (any(all_present))
THEN
2985 CALL pair_radius_setup(orb_present, all_present, orb_radius, all_radius, pair_radius)
2986 CALL build_neighbor_lists(ec_env%sac_ae, particle_set, atom2d, cell, pair_radius, &
2987 subcells=subcells, operator_type=
"ABC", nlname=
"sac_ae")
2991 IF (gth_potential_present .OR. sgp_potential_present)
THEN
2992 IF (any(ppl_present))
THEN
2993 CALL pair_radius_setup(orb_present, ppl_present, orb_radius, ppl_radius, pair_radius)
2994 CALL build_neighbor_lists(ec_env%sac_ppl, particle_set, atom2d, cell, pair_radius, &
2995 subcells=subcells, operator_type=
"ABC", nlname=
"sac_ppl")
2998 IF (any(ppnl_present))
THEN
2999 CALL pair_radius_setup(orb_present, ppnl_present, orb_radius, ppnl_radius, pair_radius)
3000 CALL build_neighbor_lists(ec_env%sap_ppnl, particle_set, atom2d, cell, pair_radius, &
3001 subcells=subcells, operator_type=
"ABBA", nlname=
"sap_ppnl")
3006 c_radius(:) = 0.0_dp
3007 dispersion_env => ec_env%dispersion_env
3008 sab_vdw => dispersion_env%sab_vdw
3009 sab_cn => dispersion_env%sab_cn
3010 IF (dispersion_env%type == xc_vdw_fun_pairpot)
THEN
3011 c_radius(:) = dispersion_env%rc_disp
3012 default_present = .true.
3013 CALL pair_radius_setup(default_present, default_present, c_radius, c_radius, pair_radius)
3014 CALL build_neighbor_lists(sab_vdw, particle_set, atom2d, cell, pair_radius, &
3015 subcells=subcells, operator_type=
"PP", nlname=
"sab_vdw")
3016 dispersion_env%sab_vdw => sab_vdw
3017 IF (dispersion_env%pp_type == vdw_pairpot_dftd3 .OR. &
3018 dispersion_env%pp_type == vdw_pairpot_dftd3bj)
THEN
3021 CALL get_atomic_kind(atomic_kind_set(ikind), z=zat)
3022 c_radius(ikind) = 4._dp*ptable(zat)%covalent_radius*bohr
3024 CALL pair_radius_setup(default_present, default_present, c_radius, c_radius, pair_radius)
3025 CALL build_neighbor_lists(sab_cn, particle_set, atom2d, cell, pair_radius, &
3026 subcells=subcells, operator_type=
"PP", nlname=
"sab_cn")
3027 dispersion_env%sab_cn => sab_cn
3032 CALL atom2d_cleanup(atom2d)
3034 DEALLOCATE (orb_present, default_present, all_present, ppl_present, ppnl_present)
3035 DEALLOCATE (orb_radius, all_radius, ppl_radius, ppnl_radius, c_radius)
3036 DEALLOCATE (pair_radius)
3039 CALL get_qs_env(qs_env, ks_env=ks_env, dft_control=dft_control)
3040 skip_load_balance_distributed = dft_control%qs_control%skip_load_balance_distributed
3041 IF (
ASSOCIATED(ec_env%task_list))
CALL deallocate_task_list(ec_env%task_list)
3042 CALL allocate_task_list(ec_env%task_list)
3043 CALL generate_qs_task_list(ks_env, ec_env%task_list, &
3044 reorder_rs_grid_ranks=.false., soft_valid=.false., &
3045 skip_load_balance_distributed=skip_load_balance_distributed, &
3046 basis_type=
"HARRIS", sab_orb_external=ec_env%sab_orb)
3048 CALL timestop(handle)
3050 END SUBROUTINE ec_build_neighborlist
3057 SUBROUTINE ec_properties(qs_env, ec_env)
3058 TYPE(qs_environment_type),
POINTER :: qs_env
3059 TYPE(energy_correction_type),
POINTER :: ec_env
3061 CHARACTER(LEN=*),
PARAMETER :: routinen =
'ec_properties'
3063 CHARACTER(LEN=8),
DIMENSION(3) :: rlab
3064 CHARACTER(LEN=default_path_length) :: filename, my_pos_voro
3065 CHARACTER(LEN=default_string_length) :: description
3066 INTEGER :: akind, handle, i, ia, iatom, idir, ikind, iounit, ispin, maxmom, nspins, &
3067 reference, should_print_bqb, should_print_voro, unit_nr, unit_nr_voro
3068 LOGICAL :: append_voro, magnetic, periodic, &
3070 REAL(kind=dp) :: charge, dd, focc, tmp
3071 REAL(kind=dp),
DIMENSION(3) :: cdip, pdip, rcc, rdip, ria, tdip
3072 REAL(kind=dp),
DIMENSION(:),
POINTER :: ref_point
3073 TYPE(atomic_kind_type),
POINTER :: atomic_kind
3074 TYPE(cell_type),
POINTER :: cell
3075 TYPE(cp_logger_type),
POINTER :: logger
3076 TYPE(cp_result_type),
POINTER :: results
3077 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_s, moments
3078 TYPE(dft_control_type),
POINTER :: dft_control
3079 TYPE(distribution_1d_type),
POINTER :: local_particles
3080 TYPE(mp_para_env_type),
POINTER :: para_env
3081 TYPE(particle_type),
DIMENSION(:),
POINTER :: particle_set
3082 TYPE(pw_env_type),
POINTER :: pw_env
3083 TYPE(pw_pool_p_type),
DIMENSION(:),
POINTER :: pw_pools
3084 TYPE(pw_pool_type),
POINTER :: auxbas_pw_pool
3085 TYPE(pw_r3d_rs_type) :: rho_elec_rspace
3086 TYPE(qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
3087 TYPE(section_vals_type),
POINTER :: ec_section, print_key, print_key_bqb, &
3090 CALL timeset(routinen, handle)
3096 logger => cp_get_default_logger()
3097 IF (logger%para_env%is_source())
THEN
3098 iounit = cp_logger_get_default_unit_nr(logger, local=.true.)
3103 NULLIFY (dft_control)
3104 CALL get_qs_env(qs_env, dft_control=dft_control)
3105 nspins = dft_control%nspins
3107 ec_section => section_vals_get_subs_vals(qs_env%input,
"DFT%ENERGY_CORRECTION")
3108 print_key => section_vals_get_subs_vals(section_vals=ec_section, &
3109 subsection_name=
"PRINT%MOMENTS")
3111 IF (btest(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file))
THEN
3113 maxmom = section_get_ival(section_vals=ec_section, &
3114 keyword_name=
"PRINT%MOMENTS%MAX_MOMENT")
3115 periodic = section_get_lval(section_vals=ec_section, &
3116 keyword_name=
"PRINT%MOMENTS%PERIODIC")
3117 reference = section_get_ival(section_vals=ec_section, &
3118 keyword_name=
"PRINT%MOMENTS%REFERENCE")
3119 magnetic = section_get_lval(section_vals=ec_section, &
3120 keyword_name=
"PRINT%MOMENTS%MAGNETIC")
3122 CALL section_vals_val_get(ec_section,
"PRINT%MOMENTS%REF_POINT", r_vals=ref_point)
3123 unit_nr = cp_print_key_unit_nr(logger=logger, basis_section=ec_section, &
3124 print_key_path=
"PRINT%MOMENTS", extension=
".dat", &
3125 middle_name=
"moments", log_filename=.false.)
3127 IF (iounit > 0)
THEN
3128 IF (unit_nr /= iounit .AND. unit_nr > 0)
THEN
3129 INQUIRE (unit=unit_nr, name=filename)
3130 WRITE (unit=iounit, fmt=
"(/,T2,A,2(/,T3,A),/)") &
3131 "MOMENTS",
"The electric/magnetic moments are written to file:", &
3134 WRITE (unit=iounit, fmt=
"(/,T2,A)")
"ELECTRIC/MAGNETIC MOMENTS"
3139 cpabort(
"Periodic moments not implemented with EC")
3141 cpassert(maxmom < 2)
3142 cpassert(.NOT. magnetic)
3143 IF (maxmom == 1)
THEN
3144 CALL get_qs_env(qs_env=qs_env, cell=cell, para_env=para_env)
3146 CALL get_reference_point(rcc, qs_env=qs_env, reference=reference, ref_point=ref_point)
3149 CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, &
3150 qs_kind_set=qs_kind_set, local_particles=local_particles)
3151 DO ikind = 1,
SIZE(local_particles%n_el)
3152 DO ia = 1, local_particles%n_el(ikind)
3153 iatom = local_particles%list(ikind)%array(ia)
3155 ria = pbc(particle_set(iatom)%r - rcc, cell) + rcc
3157 atomic_kind => particle_set(iatom)%atomic_kind
3158 CALL get_atomic_kind(atomic_kind, kind_number=akind)
3159 CALL get_qs_kind(qs_kind_set(akind), core_charge=charge)
3160 cdip(1:3) = cdip(1:3) - charge*ria(1:3)
3163 CALL para_env%sum(cdip)
3166 CALL ec_efield_integrals(qs_env, ec_env, rcc)
3169 DO ispin = 1, nspins
3171 CALL dbcsr_dot(ec_env%matrix_p(ispin, 1)%matrix, &
3172 ec_env%efield%dipmat(idir)%matrix, tmp)
3173 pdip(idir) = pdip(idir) + tmp
3178 CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s)
3180 CALL dbcsr_allocate_matrix_set(moments, 4)
3182 ALLOCATE (moments(i)%matrix)
3183 CALL dbcsr_copy(moments(i)%matrix, matrix_s(1)%matrix,
"Moments")
3184 CALL dbcsr_set(moments(i)%matrix, 0.0_dp)
3186 CALL build_local_moment_matrix(qs_env, moments, 1, ref_point=rcc)
3189 IF (nspins == 2) focc = 1.0_dp
3191 DO ispin = 1, nspins
3193 CALL dbcsr_dot(ec_env%matrix_z(ispin)%matrix, moments(idir)%matrix, tmp)
3194 rdip(idir) = rdip(idir) + tmp
3197 CALL dbcsr_deallocate_matrix_set(moments)
3199 tdip = -(rdip + pdip + cdip)
3200 IF (unit_nr > 0)
THEN
3201 WRITE (unit_nr,
"(T3,A)")
"Dipoles are based on the traditional operator."
3202 dd = sqrt(sum(tdip(1:3)**2))*debye
3203 WRITE (unit_nr,
"(T3,A)")
"Dipole moment [Debye]"
3204 WRITE (unit_nr,
"(T5,3(A,A,F14.8,1X),T60,A,T67,F14.8)") &
3205 (trim(rlab(i)),
"=", tdip(i)*debye, i=1, 3),
"Total=", dd
3210 CALL cp_print_key_finished_output(unit_nr=unit_nr, logger=logger, &
3211 basis_section=ec_section, print_key_path=
"PRINT%MOMENTS")
3212 CALL get_qs_env(qs_env=qs_env, results=results)
3213 description =
"[DIPOLE]"
3214 CALL cp_results_erase(results=results, description=description)
3215 CALL put_results(results=results, description=description, values=tdip(1:3))
3219 print_key_voro => section_vals_get_subs_vals(ec_section,
"PRINT%VORONOI")
3220 print_key_bqb => section_vals_get_subs_vals(ec_section,
"PRINT%E_DENSITY_BQB")
3221 IF (btest(cp_print_key_should_output(logger%iter_info, print_key_voro), cp_p_file))
THEN
3222 should_print_voro = 1
3224 should_print_voro = 0
3226 IF (btest(cp_print_key_should_output(logger%iter_info, print_key_bqb), cp_p_file))
THEN
3227 should_print_bqb = 1
3229 should_print_bqb = 0
3231 IF ((should_print_voro /= 0) .OR. (should_print_bqb /= 0))
THEN
3233 CALL get_qs_env(qs_env=qs_env, &
3235 CALL pw_env_get(pw_env=pw_env, &
3236 auxbas_pw_pool=auxbas_pw_pool, &
3238 CALL auxbas_pw_pool%create_pw(pw=rho_elec_rspace)
3240 IF (dft_control%nspins > 1)
THEN
3243 CALL pw_copy(ec_env%rhoout_r(1), rho_elec_rspace)
3244 CALL pw_axpy(ec_env%rhoout_r(2), rho_elec_rspace)
3246 CALL pw_axpy(ec_env%rhoz_r(1), rho_elec_rspace)
3247 CALL pw_axpy(ec_env%rhoz_r(2), rho_elec_rspace)
3251 CALL pw_copy(ec_env%rhoout_r(1), rho_elec_rspace)
3252 CALL pw_axpy(ec_env%rhoz_r(1), rho_elec_rspace)
3255 IF (should_print_voro /= 0)
THEN
3256 CALL section_vals_val_get(print_key_voro,
"OUTPUT_TEXT", l_val=voro_print_txt)
3257 IF (voro_print_txt)
THEN
3258 append_voro = section_get_lval(ec_section,
"PRINT%VORONOI%APPEND")
3259 my_pos_voro =
"REWIND"
3260 IF (append_voro)
THEN
3261 my_pos_voro =
"APPEND"
3263 unit_nr_voro = cp_print_key_unit_nr(logger, ec_section,
"PRINT%VORONOI", extension=
".voronoi", &
3264 file_position=my_pos_voro, log_filename=.false.)
3272 CALL entry_voronoi_or_bqb(should_print_voro, should_print_bqb, print_key_voro, print_key_bqb, &
3273 unit_nr_voro, qs_env, rho_elec_rspace)
3275 CALL auxbas_pw_pool%give_back_pw(rho_elec_rspace)
3277 IF (unit_nr_voro > 0)
THEN
3278 CALL cp_print_key_finished_output(unit_nr_voro, logger, ec_section,
"PRINT%VORONOI")
3283 CALL timestop(handle)
3285 END SUBROUTINE ec_properties
3298 SUBROUTINE ec_ls_init(qs_env, matrix_ks, matrix_s)
3299 TYPE(qs_environment_type),
POINTER :: qs_env
3300 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_ks, matrix_s
3302 CHARACTER(len=*),
PARAMETER :: routinen =
'ec_ls_init'
3304 INTEGER :: handle, ispin, nspins
3305 TYPE(dft_control_type),
POINTER :: dft_control
3306 TYPE(energy_correction_type),
POINTER :: ec_env
3307 TYPE(ls_scf_env_type),
POINTER :: ls_env
3309 CALL timeset(routinen, handle)
3311 CALL get_qs_env(qs_env=qs_env, &
3312 dft_control=dft_control, &
3314 nspins = dft_control%nspins
3315 ls_env => ec_env%ls_env
3318 CALL matrix_ls_create(matrix_ls=ls_env%matrix_s, matrix_qs=matrix_s(1, 1)%matrix, &
3319 ls_mstruct=ls_env%ls_mstruct)
3321 IF (
ALLOCATED(ls_env%matrix_p))
THEN
3322 DO ispin = 1,
SIZE(ls_env%matrix_p)
3323 CALL dbcsr_release(ls_env%matrix_p(ispin))
3326 ALLOCATE (ls_env%matrix_p(nspins))
3329 DO ispin = 1, nspins
3330 CALL dbcsr_create(ls_env%matrix_p(ispin), template=ls_env%matrix_s, &
3331 matrix_type=dbcsr_type_no_symmetry)
3334 ALLOCATE (ls_env%matrix_ks(nspins))
3335 DO ispin = 1, nspins
3336 CALL dbcsr_create(ls_env%matrix_ks(ispin), template=ls_env%matrix_s, &
3337 matrix_type=dbcsr_type_no_symmetry)
3341 CALL ls_scf_init_matrix_s(matrix_s(1, 1)%matrix, ls_env)
3345 DO ispin = 1, nspins
3346 CALL matrix_qs_to_ls(matrix_ls=ls_env%matrix_ks(ispin), &
3347 matrix_qs=matrix_ks(ispin, 1)%matrix, &
3348 ls_mstruct=ls_env%ls_mstruct, &
3352 CALL timestop(handle)
3354 END SUBROUTINE ec_ls_init
3370 SUBROUTINE ec_ls_solver(qs_env, matrix_p, matrix_w, ec_ls_method)
3371 TYPE(qs_environment_type),
POINTER :: qs_env
3372 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_p, matrix_w
3373 INTEGER,
INTENT(IN) :: ec_ls_method
3375 CHARACTER(LEN=*),
PARAMETER :: routinen =
'ec_ls_solver'
3377 INTEGER :: handle, ispin, nelectron_spin_real, &
3379 INTEGER,
DIMENSION(2) :: nmo
3380 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: wmat
3381 TYPE(dbcsr_type),
ALLOCATABLE,
DIMENSION(:) :: matrix_ks_deviation
3382 TYPE(dft_control_type),
POINTER :: dft_control
3383 TYPE(energy_correction_type),
POINTER :: ec_env
3384 TYPE(ls_scf_env_type),
POINTER :: ls_env
3385 TYPE(mp_para_env_type),
POINTER :: para_env
3387 CALL timeset(routinen, handle)
3390 CALL get_qs_env(qs_env, &
3391 dft_control=dft_control, &
3393 nspins = dft_control%nspins
3394 ec_env => qs_env%ec_env
3395 ls_env => ec_env%ls_env
3398 CALL get_qs_env(qs_env=qs_env, nelectron_spin=nmo)
3399 IF (nspins == 1) nmo(1) = nmo(1)/2
3400 ls_env%homo_spin(:) = 0.0_dp
3401 ls_env%lumo_spin(:) = 0.0_dp
3403 ALLOCATE (matrix_ks_deviation(nspins))
3404 DO ispin = 1, nspins
3405 CALL dbcsr_create(matrix_ks_deviation(ispin), template=ls_env%matrix_ks(ispin))
3406 CALL dbcsr_set(matrix_ks_deviation(ispin), 0.0_dp)
3410 IF (ls_env%has_s_preconditioner)
THEN
3411 DO ispin = 1, nspins
3412 CALL apply_matrix_preconditioner(ls_env%matrix_ks(ispin),
"forward", &
3413 ls_env%matrix_bs_sqrt, ls_env%matrix_bs_sqrt_inv)
3415 CALL dbcsr_filter(ls_env%matrix_ks(ispin), ls_env%eps_filter)
3419 DO ispin = 1, nspins
3420 nelectron_spin_real = ls_env%nelectron_spin(ispin)
3421 IF (ls_env%nspins == 1) nelectron_spin_real = nelectron_spin_real/2
3423 SELECT CASE (ec_ls_method)
3424 CASE (ec_matrix_sign)
3425 CALL density_matrix_sign(ls_env%matrix_p(ispin), &
3426 ls_env%mu_spin(ispin), &
3428 ls_env%sign_method, &
3429 ls_env%sign_order, &
3430 ls_env%matrix_ks(ispin), &
3432 ls_env%matrix_s_inv, &
3433 nelectron_spin_real, &
3436 CASE (ec_matrix_trs4)
3437 CALL density_matrix_trs4( &
3438 ls_env%matrix_p(ispin), &
3439 ls_env%matrix_ks(ispin), &
3440 ls_env%matrix_s_sqrt_inv, &
3441 nelectron_spin_real, &
3442 ec_env%eps_default, &
3443 ls_env%homo_spin(ispin), &
3444 ls_env%lumo_spin(ispin), &
3445 ls_env%mu_spin(ispin), &
3446 matrix_ks_deviation=matrix_ks_deviation(ispin), &
3447 dynamic_threshold=ls_env%dynamic_threshold, &
3448 eps_lanczos=ls_env%eps_lanczos, &
3449 max_iter_lanczos=ls_env%max_iter_lanczos)
3451 CASE (ec_matrix_tc2)
3452 CALL density_matrix_tc2( &
3453 ls_env%matrix_p(ispin), &
3454 ls_env%matrix_ks(ispin), &
3455 ls_env%matrix_s_sqrt_inv, &
3456 nelectron_spin_real, &
3457 ec_env%eps_default, &
3458 ls_env%homo_spin(ispin), &
3459 ls_env%lumo_spin(ispin), &
3460 non_monotonic=ls_env%non_monotonic, &
3461 eps_lanczos=ls_env%eps_lanczos, &
3462 max_iter_lanczos=ls_env%max_iter_lanczos)
3469 IF (ls_env%has_s_preconditioner)
THEN
3470 DO ispin = 1, nspins
3472 CALL apply_matrix_preconditioner(ls_env%matrix_p(ispin),
"forward", &
3473 ls_env%matrix_bs_sqrt, ls_env%matrix_bs_sqrt_inv)
3475 CALL dbcsr_filter(ls_env%matrix_p(ispin), ls_env%eps_filter)
3480 IF (nspins == 1)
CALL dbcsr_scale(ls_env%matrix_p(1), 2.0_dp)
3482 IF (ls_env%report_all_sparsities)
CALL post_scf_sparsities(ls_env)
3486 DO ispin = 1, nspins
3487 CALL matrix_ls_to_qs(matrix_qs=matrix_p(ispin, 1)%matrix, &
3488 matrix_ls=ls_env%matrix_p(ispin), &
3489 ls_mstruct=ls_env%ls_mstruct, &
3493 wmat => matrix_w(:, 1)
3494 CALL calculate_w_matrix_ls(wmat, ec_env%ls_env)
3497 CALL dbcsr_release(ls_env%matrix_s)
3498 IF (ls_env%has_s_preconditioner)
THEN
3499 CALL dbcsr_release(ls_env%matrix_bs_sqrt)
3500 CALL dbcsr_release(ls_env%matrix_bs_sqrt_inv)
3502 IF (ls_env%needs_s_inv)
THEN
3503 CALL dbcsr_release(ls_env%matrix_s_inv)
3505 IF (ls_env%use_s_sqrt)
THEN
3506 CALL dbcsr_release(ls_env%matrix_s_sqrt)
3507 CALL dbcsr_release(ls_env%matrix_s_sqrt_inv)
3510 DO ispin = 1,
SIZE(ls_env%matrix_ks)
3511 CALL dbcsr_release(ls_env%matrix_ks(ispin))
3513 DEALLOCATE (ls_env%matrix_ks)
3515 DO ispin = 1, nspins
3516 CALL dbcsr_release(matrix_ks_deviation(ispin))
3518 DEALLOCATE (matrix_ks_deviation)
3520 CALL timestop(handle)
3522 END SUBROUTINE ec_ls_solver
3540 SUBROUTINE ec_ot_diag_solver(qs_env, ec_env, matrix_ks, matrix_s, matrix_p, matrix_w)
3541 TYPE(qs_environment_type),
POINTER :: qs_env
3542 TYPE(energy_correction_type),
POINTER :: ec_env
3543 TYPE(dbcsr_p_type),
DIMENSION(:, :),
INTENT(IN), &
3544 POINTER :: matrix_ks, matrix_s
3545 TYPE(dbcsr_p_type),
DIMENSION(:, :), &
3546 INTENT(INOUT),
POINTER :: matrix_p, matrix_w
3548 CHARACTER(len=*),
PARAMETER :: routinen =
'ec_ot_diag_solver'
3550 INTEGER :: handle, homo, ikind, iounit, ispin, &
3551 max_iter, nao, nkind, nmo, nspins
3552 INTEGER,
DIMENSION(2) :: nelectron_spin
3553 REAL(kind=dp),
DIMENSION(:),
POINTER :: eigenvalues
3554 TYPE(atomic_kind_type),
DIMENSION(:),
POINTER :: atomic_kind_set
3555 TYPE(cp_blacs_env_type),
POINTER :: blacs_env
3556 TYPE(cp_fm_type) :: sv
3557 TYPE(cp_fm_type),
POINTER :: mo_coeff
3558 TYPE(cp_logger_type),
POINTER :: logger
3559 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: p_rmpv
3560 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: rho_ao
3561 TYPE(dft_control_type),
POINTER :: dft_control
3562 TYPE(gto_basis_set_type),
POINTER :: basis_set, harris_basis
3563 TYPE(mo_set_type),
DIMENSION(:),
POINTER :: mos
3564 TYPE(mp_para_env_type),
POINTER :: para_env
3565 TYPE(particle_type),
DIMENSION(:),
POINTER :: particle_set
3566 TYPE(preconditioner_type),
POINTER :: local_preconditioner
3567 TYPE(qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
3568 TYPE(qs_kind_type),
POINTER :: qs_kind
3569 TYPE(qs_rho_type),
POINTER :: rho
3571 CALL timeset(routinen, handle)
3573 logger => cp_get_default_logger()
3574 iounit = cp_logger_get_default_unit_nr(logger)
3576 cpassert(
ASSOCIATED(qs_env))
3577 cpassert(
ASSOCIATED(ec_env))
3578 cpassert(
ASSOCIATED(matrix_ks))
3579 cpassert(
ASSOCIATED(matrix_s))
3580 cpassert(
ASSOCIATED(matrix_p))
3581 cpassert(
ASSOCIATED(matrix_w))
3583 CALL get_qs_env(qs_env=qs_env, &
3584 atomic_kind_set=atomic_kind_set, &
3585 blacs_env=blacs_env, &
3586 dft_control=dft_control, &
3587 nelectron_spin=nelectron_spin, &
3588 para_env=para_env, &
3589 particle_set=particle_set, &
3590 qs_kind_set=qs_kind_set)
3591 nspins = dft_control%nspins
3598 IF (dft_control%qs_control%do_ls_scf)
THEN
3599 CALL ec_mos_init(qs_env, matrix_s(1, 1)%matrix)
3602 CALL get_qs_env(qs_env, mos=mos)
3609 CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, nkind=nkind)
3610 CALL uppercase(ec_env%basis)
3614 IF (ec_env%basis ==
"HARRIS")
THEN
3616 qs_kind => qs_kind_set(ikind)
3618 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type=
"ORB")
3620 CALL get_qs_kind(qs_kind=qs_kind, basis_set=harris_basis, basis_type=
"HARRIS")
3622 IF (basis_set%name .NE. harris_basis%name)
THEN
3623 cpabort(
"OT-Diag initial guess: Harris and ground state need to use the same basis")
3628 IF (ec_env%mao)
THEN
3629 cpabort(
"OT-Diag initial guess: not implemented for MAOs")
3633 SELECT CASE (ec_env%ec_initial_guess)
3636 p_rmpv => matrix_p(:, 1)
3638 CALL calculate_atomic_block_dm(p_rmpv, matrix_s(1, 1)%matrix, atomic_kind_set, qs_kind_set, &
3639 nspins, nelectron_spin, iounit, para_env)
3643 CALL get_qs_env(qs_env, rho=rho)
3644 CALL qs_rho_get(rho, rho_ao_kp=rho_ao)
3645 p_rmpv => rho_ao(:, 1)
3648 cpabort(
"Unknown inital guess for OT-Diagonalization (Harris functional)")
3651 DO ispin = 1, nspins
3652 CALL get_mo_set(mo_set=mos(ispin), &
3653 mo_coeff=mo_coeff, &
3659 CALL cp_fm_set_all(mo_coeff, 0.0_dp)
3660 CALL cp_fm_init_random(mo_coeff, nmo)
3662 CALL cp_fm_create(sv, mo_coeff%matrix_struct,
"SV")
3665 CALL cp_dbcsr_sm_fm_multiply(matrix_s(1, 1)%matrix, mo_coeff, sv, nmo)
3666 CALL cp_dbcsr_sm_fm_multiply(p_rmpv(ispin)%matrix, sv, mo_coeff, homo)
3667 CALL cp_fm_release(sv)
3670 CALL make_basis_sm(mo_coeff, nmo, matrix_s(1, 1)%matrix)
3674 NULLIFY (local_preconditioner)
3675 ALLOCATE (local_preconditioner)
3676 CALL init_preconditioner(local_preconditioner, para_env=para_env, &
3677 blacs_env=blacs_env)
3678 DO ispin = 1, nspins
3679 CALL make_preconditioner(local_preconditioner, &
3680 precon_type=ot_precond_full_single_inverse, &
3681 solver_type=ot_precond_solver_default, &
3682 matrix_h=matrix_ks(ispin, 1)%matrix, &
3683 matrix_s=matrix_s(ispin, 1)%matrix, &
3684 convert_precond_to_dbcsr=.true., &
3685 mo_set=mos(ispin), energy_gap=0.2_dp)
3687 CALL get_mo_set(mos(ispin), &
3688 mo_coeff=mo_coeff, &
3689 eigenvalues=eigenvalues, &
3692 CALL ot_eigensolver(matrix_h=matrix_ks(ispin, 1)%matrix, &
3693 matrix_s=matrix_s(1, 1)%matrix, &
3694 matrix_c_fm=mo_coeff, &
3696 eps_gradient=ec_env%eps_default, &
3697 iter_max=max_iter, &
3699 CALL calculate_subspace_eigenvalues(mo_coeff, matrix_ks(ispin, 1)%matrix, &
3700 evals_arg=eigenvalues, do_rotation=.true.)
3703 CALL destroy_preconditioner(local_preconditioner)
3704 DEALLOCATE (local_preconditioner)
3707 CALL copy_fm_to_dbcsr(mos(ispin)%mo_coeff, &
3708 mos(ispin)%mo_coeff_b)
3712 DO ispin = 1, nspins
3713 CALL calculate_density_matrix(mos(ispin), matrix_p(ispin, 1)%matrix)
3715 CALL calculate_w_matrix(mos(ispin), matrix_w(ispin, 1)%matrix)
3719 IF (dft_control%qs_control%do_ls_scf)
THEN
3720 DO ispin = 1, nspins
3721 CALL deallocate_mo_set(mos(ispin))
3723 IF (
ASSOCIATED(qs_env%mos))
THEN
3724 DO ispin = 1,
SIZE(qs_env%mos)
3725 CALL deallocate_mo_set(qs_env%mos(ispin))
3727 DEALLOCATE (qs_env%mos)
3731 CALL timestop(handle)
3733 END SUBROUTINE ec_ot_diag_solver
Contains ADMM methods which only require the density matrix.
subroutine, public admm_dm_calc_rho_aux(qs_env)
Entry methods: Calculates auxiliary density matrix from primary one.
Contains ADMM methods which require molecular orbitals.
subroutine, public admm_mo_calc_rho_aux(qs_env)
...
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
subroutine, public get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, npgf_sum, nshell_sum, maxder, short_kind_radius, npgf_seg_sum)
...
collects all references to literature in CP2K as new algorithms / method are included from literature...
integer, save, public belleflamme2023
Handles all functions related to the CELL.
methods related to the blacs parallel environment
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
subroutine, public dbcsr_scale(matrix, alpha_scalar)
...
subroutine, public dbcsr_desymmetrize(matrix_a, matrix_b)
...
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)
...
subroutine, public dbcsr_init_p(matrix)
...
subroutine, public dbcsr_filter(matrix, eps)
...
subroutine, public dbcsr_set(matrix, alpha)
...
subroutine, public dbcsr_release(matrix)
...
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
subroutine, public dbcsr_dot(matrix_a, matrix_b, trace)
Computes the dot product of two matrices, also known as the trace of their matrix product.
Routines that link DBCSR and CP2K concepts together.
subroutine, public cp_dbcsr_alloc_block_from_nbl(matrix, sab_orb, desymmetrize)
allocate the blocks of a dbcsr based on the neighbor list
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 copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
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.
Basic linear algebra operations for full matrices.
subroutine, public cp_fm_column_scale(matrixa, scaling)
scales column i of matrix a with scaling(i)
subroutine, public cp_fm_triangular_invert(matrix_a, uplo_tr)
inverts a triangular matrix
various cholesky decomposition related routines
subroutine, public cp_fm_cholesky_decompose(matrix, n, info_out)
used to replace a symmetric positive def. matrix M with its cholesky decomposition U: M = U^T * U,...
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
subroutine, public choose_eigv_solver(matrix, eigenvectors, eigenvalues, info)
Choose the Eigensolver depending on which library is available ELPA seems to be unstable for small sy...
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_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_init_random(matrix, ncol, start_col)
fills a matrix with random numbers
subroutine, public cp_fm_to_fm_triangular(msource, mtarget, uplo)
copy just a triangular 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
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
integer, parameter, public cp_p_file
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...
set of type/routines to handle the storage of results in force_envs
subroutine, public cp_results_erase(results, description, nval)
erase a part of result_list
set of type/routines to handle the storage of results in force_envs
real(kind=dp) function, public cp_unit_from_cp2k(value, unit_str, defaults, power)
converts from the internal cp2k units to the given unit
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
stores a mapping of 2D info (e.g. matrix) on a 2D processor distribution (i.e. blacs grid) where cpus...
lower level routines for linear scaling SCF
subroutine, public density_matrix_trs4(matrix_p, matrix_ks, matrix_s_sqrt_inv, nelectron, threshold, e_homo, e_lumo, e_mu, dynamic_threshold, matrix_ks_deviation, max_iter_lanczos, eps_lanczos, converged, iounit)
compute the density matrix using a trace-resetting algorithm
subroutine, public ls_scf_init_matrix_s(matrix_s, ls_scf_env)
initialize S matrix related properties (sqrt, inverse...) Might be factored-out since this seems comm...
subroutine, public apply_matrix_preconditioner(matrix, direction, matrix_bs_sqrt, matrix_bs_sqrt_inv)
apply a preconditioner either forward (precondition) inv(sqrt(bs)) * A * inv(sqrt(bs)) backward (rest...
subroutine, public density_matrix_sign(matrix_p, mu, fixed_mu, sign_method, sign_order, matrix_ks, matrix_s, matrix_s_inv, nelectron, threshold, sign_symmetric, submatrix_sign_method, matrix_s_sqrt_inv)
compute the density matrix with a trace that is close to nelectron. take a mu as input,...
subroutine, public density_matrix_tc2(matrix_p, matrix_ks, matrix_s_sqrt_inv, nelectron, threshold, e_homo, e_lumo, non_monotonic, eps_lanczos, max_iter_lanczos, iounit)
compute the density matrix using a non monotonic trace conserving algorithm based on SIAM DOI....
Routines for a linear scaling quickstep SCF run based on the density matrix, with a focus on the inte...
subroutine, public matrix_ls_to_qs(matrix_qs, matrix_ls, ls_mstruct, covariant, keep_sparsity)
second link to QS, copy a LS matrix to QS matrix used to isolate QS style matrices from LS style will...
subroutine, public matrix_ls_create(matrix_ls, matrix_qs, ls_mstruct)
create a matrix for use (and as a template) in ls based on a qs template
subroutine, public matrix_qs_to_ls(matrix_ls, matrix_qs, ls_mstruct, covariant)
first link to QS, copy a QS matrix to LS matrix used to isolate QS style matrices from LS style will ...
Types needed for a linear scaling quickstep SCF run based on the density matrix.
Routines for a linear scaling quickstep SCF run based on the density matrix.
subroutine, public post_scf_sparsities(ls_scf_env)
Report on the sparsities of various interesting matrices.
subroutine, public calculate_w_matrix_ls(matrix_w, ls_scf_env)
Compute matrix_w as needed for the forces.
Calculates the energy contribution and the mo_derivative of a static electric field (nonperiodic)
subroutine, public ec_efield_local_operator(qs_env, ec_env, calculate_forces)
...
subroutine, public ec_efield_integrals(qs_env, ec_env, rpoint)
...
Types needed for a for a Energy Correction.
Routines for an external energy correction on top of a Kohn-Sham calculation.
subroutine, public ec_ext_energy(qs_env, ec_env, calculate_forces)
External energy method.
Routines used for Harris functional Kohn-Sham calculation.
subroutine, public ec_mos_init(qs_env, matrix_s)
Allocate and initiate molecular orbitals environment.
subroutine, public create_kernel(qs_env, vxc, vxc_tau, rho, rho1_r, rho1_g, tau1_r, xc_section, compute_virial, virial_xc)
Creation of second derivative xc-potential.
Routines for an energy correction on top of a Kohn-Sham calculation.
subroutine, public energy_correction(qs_env, ec_init, calculate_forces)
Energy Correction to a Kohn-Sham simulation Available energy corrections: (1) Harris energy functiona...
Definition of the atomic potential types.
Routines to calculate EXX in RPA and energy correction methods.
subroutine, public calculate_exx(qs_env, unit_nr, hfx_sections, x_data, do_gw, do_admm, calc_forces, reuse_hfx, do_im_time, e_ex_from_gw, e_admm_from_gw, t3)
...
subroutine, public add_exx_to_rhs(rhs, qs_env, ext_hfx_section, x_data, recalc_integrals, do_admm, do_ec, do_exx, reuse_hfx)
Add the EXX contribution to the RHS of the Z-vector equation, namely the HF Hamiltonian.
Defines the basic variable types.
integer, parameter, public dp
integer, parameter, public default_string_length
integer, parameter, public default_path_length
Calculate MAO's and analyze wavefunctions.
subroutine, public mao_generate_basis(qs_env, mao_coef, ref_basis_set, pmat_external, smat_external, molecular, max_iter, eps_grad, nmao_external, eps1_mao, iolevel, unit_nr)
...
Collection of simple mathematical functions and subroutines.
Interface to the message passing library MPI.
Define the data structure for the molecule information.
Calculates the moment integrals <a|r^m|b>
subroutine, public get_reference_point(rpoint, drpoint, qs_env, fist_env, reference, ref_point, ifirst, ilast)
...
Define the data structure for the particle information.
Periodic Table related data definitions.
type(atom), dimension(0:nelem), public ptable
Definition of physical constants:
real(kind=dp), parameter, public pascal
real(kind=dp), parameter, public bohr
real(kind=dp), parameter, public debye
subroutine, public init_preconditioner(preconditioner_env, para_env, blacs_env)
...
subroutine, public destroy_preconditioner(preconditioner_env)
...
computes preconditioners, and implements methods to apply them currently used in qs_ot
subroutine, public make_preconditioner(preconditioner_env, precon_type, solver_type, matrix_h, matrix_s, matrix_t, mo_set, energy_gap, convert_precond_to_dbcsr, chol_type)
...
container for various plainwaves related things
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
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Routine to return block diagonal density matrix. Blocks correspond to the atomic densities.
subroutine, public calculate_atomic_block_dm(pmatrix, matrix_s, atomic_kind_set, qs_kind_set, nspin, nelectron_spin, ounit, para_env)
returns a block diagonal density matrix. Blocks correspond to the atomic densities.
Calculate the plane wave density by collocating the primitive Gaussian functions (pgf).
subroutine, public calculate_rho_elec(matrix_p, matrix_p_kp, rho, rho_gspace, total_rho, ks_env, soft_valid, compute_tau, compute_grad, basis_type, der_type, idir, task_list_external, pw_env_external)
computes the density corresponding to a given density matrix on the grid
Calculation of the energies concerning the core charge distribution.
subroutine, public calculate_ecore_overlap(qs_env, para_env, calculate_forces, molecular, e_overlap_core, atecc)
Calculate the overlap energy of the core charge distribution.
Calculation of the core Hamiltonian integral matrix <a|H|b> over Cartesian Gaussian-type functions.
subroutine, public kinetic_energy_matrix(qs_env, matrixkp_t, matrix_t, matrix_p, matrix_name, calculate_forces, nderivative, sab_orb, eps_filter, basis_type, debug_forces, debug_stress)
Calculate kinetic energy matrix and possible relativistic correction.
subroutine, public core_matrices(qs_env, matrix_h, matrix_p, calculate_forces, nder, ec_env, dcdr_env, ec_env_matrices, basis_type, debug_forces, debug_stress, atcore)
...
collects routines that calculate density matrices
Calculation of dispersion using pair potentials.
subroutine, public calculate_dispersion_pairpot(qs_env, dispersion_env, energy, calculate_forces, atevdw)
...
Definition of disperson types for DFT calculations.
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, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, 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, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
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, rhoz_cneo_set, 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, do_rixs, tb_tblite)
Set the QUICKSTEP environment.
subroutine, public zero_qs_force(qs_force)
Initialize a Quickstep force data structure.
subroutine, public total_qs_force(force, qs_force, atomic_kind_set)
Get current total force.
Integrate single or product functions over a potential on a RS grid.
Define the quickstep kind type and their sub types.
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, basis_type, total_zeff_corr, npgf_seg, cneo_potential_present, nkind_q, natom_q)
Get attributes of an atomic kind set.
Calculation of kinetic energy matrix and forces.
subroutine, public build_kinetic_matrix(ks_env, matrix_t, matrixkp_t, matrix_name, basis_type, sab_nl, calculate_forces, matrix_p, matrixkp_p, eps_filter, nderivative)
Calculation of the kinetic energy matrix over Cartesian Gaussian functions.
routines that build the Kohn-Sham matrix (i.e calculate the coulomb and xc parts
subroutine, public calc_rho_tot_gspace(rho_tot_gspace, qs_env, rho, skip_nuclear_density)
...
Calculate the KS reference potentials.
subroutine, public ks_ref_potential(qs_env, vh_rspace, vxc_rspace, vtau_rspace, vadmm_rspace, ehartree, exc, h_stress)
calculate the Kohn-Sham reference potential
collects routines that perform operations directly related to MOs
subroutine, public make_basis_sm(vmatrix, ncol, matrix_s)
returns an S-orthonormal basis v (v^T S v ==1)
Definition and initialisation of the mo data type.
subroutine, public deallocate_mo_set(mo_set)
Deallocate a wavefunction data structure.
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.
Calculates the moment integrals <a|r^m|b> and <a|r x d/dr|b>
subroutine, public build_local_moment_matrix(qs_env, moments, nmoments, ref_point, ref_points, basis_type)
...
Define the neighbor list data types and the corresponding functionality.
Generate the atomic neighbor lists.
subroutine, public atom2d_cleanup(atom2d)
free the internals of atom2d
subroutine, public pair_radius_setup(present_a, present_b, radius_a, radius_b, pair_radius, prmin)
...
subroutine, public build_neighbor_lists(ab_list, particle_set, atom, cell, pair_radius, subcells, mic, symmetric, molecular, subset_of_mol, current_subset, operator_type, nlname, atomb_to_keep)
Build simple pair neighbor lists.
subroutine, public atom2d_build(atom2d, distribution_1d, distribution_2d, atomic_kind_set, molecule_set, molecule_only, particle_set)
Build some distribution structure of atoms, refactored from build_qs_neighbor_lists.
Routines for the construction of the coefficients for the expansion of the atomic densities rho1_hard...
subroutine, public build_oce_matrices(intac, calculate_forces, nder, qs_kind_set, particle_set, sap_oce, eps_fit)
Set up the sparse matrix for the coefficients of one center expansions This routine uses the same log...
subroutine, public allocate_oce_set(oce_set, nkind)
Allocate and initialize the matrix set of oce coefficients.
subroutine, public create_oce_set(oce_set)
...
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)
...
Calculation of overlap matrix, its derivatives and forces.
subroutine, public build_overlap_matrix(ks_env, matrix_s, matrixkp_s, matrix_name, nderivative, basis_type_a, basis_type_b, sab_nl, calculate_forces, matrix_p, matrixkp_p)
Calculation of the overlap matrix over Cartesian Gaussian functions.
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...
subroutine, public qs_vxc_create(ks_env, rho_struct, xc_section, vxc_rho, vxc_tau, exc, just_energy, edisp, dispersion_env, adiabatic_rescale_factor, pw_env_external)
calculates and allocates the xc potential, already reducing it to the dependence on rho and the one o...
Calculate the CPKS equation and the resulting forces.
subroutine, public response_force(qs_env, vh_rspace, vxc_rspace, vtau_rspace, vadmm_rspace, matrix_hz, matrix_pz, matrix_pz_admm, matrix_wz, zehartree, zexc, zexc_aux_fit, rhopz_r, p_env, ex_env, debug)
...
subroutine, public response_calculation(qs_env, ec_env)
Initializes solver of linear response equation for energy correction.
Utilities for rtp in combination with admm methods adapted routines from admm_method (author Manuel G...
subroutine, public rtp_admm_calc_rho_aux(qs_env)
Compute the ADMM density matrix in case of rtp (complex MO's)
Utilities for string manipulations.
elemental subroutine, public uppercase(string)
Convert all lower case characters in a string to upper case.
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
The module to read/write TREX IO files for interfacing CP2K with other programs.
subroutine, public write_trexio(qs_env, trexio_section, energy_derivative)
Write a trexio file.
subroutine, public write_stress_tensor(pv_virial, iw, cell, unit_string, numerical)
Print stress tensor to output file.
subroutine, public write_stress_tensor_components(virial, iw, cell, unit_string)
...
pure real(kind=dp) function, public one_third_sum_diag(a)
...
subroutine, public zero_virial(virial, reset)
...
subroutine, public symmetrize_virial(virial)
Symmetrize the virial components.
Interface for Voronoi Integration and output of BQB files.
subroutine, public entry_voronoi_or_bqb(do_voro, do_bqb, input_voro, input_bqb, unit_voro, qs_env, rspace_pw)
Does a Voronoi integration of density or stores the density to compressed BQB format.
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...
contains arbitrary information which need to be stored
structure to store local (to a processor) ordered lists of integers.
distributes pairs on a 2d grid of processors
Contains information on the energy correction functional for KG.
stores all the informations relevant to an mpi environment
contained for different pw related things
environment for the poisson solver
to create arrays of pools
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
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.