62 USE dbcsr_api,
ONLY: &
63 dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_desymmetrize, dbcsr_distribution_type, &
64 dbcsr_dot, dbcsr_filter, dbcsr_get_info, dbcsr_init_p, dbcsr_multiply, dbcsr_p_type, &
65 dbcsr_release, dbcsr_scale, dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry, &
130 USE pw_types,
ONLY: pw_c1d_gs_type,&
185 #include "./base/base_uses.f90"
193 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'energy_corrections'
212 TYPE(qs_environment_type),
POINTER :: qs_env
213 LOGICAL,
INTENT(IN),
OPTIONAL :: ec_init, calculate_forces
215 CHARACTER(len=*),
PARAMETER :: routinen =
'energy_correction'
217 INTEGER :: handle, unit_nr
218 LOGICAL :: my_calc_forces
219 TYPE(cp_logger_type),
POINTER :: logger
220 TYPE(energy_correction_type),
POINTER :: ec_env
221 TYPE(qs_energy_type),
POINTER :: energy
222 TYPE(qs_force_type),
DIMENSION(:),
POINTER :: ks_force
224 CALL timeset(routinen, handle)
227 IF (logger%para_env%is_source())
THEN
239 IF (.NOT. ec_env%do_skip)
THEN
241 ec_env%should_update = .true.
242 IF (
PRESENT(ec_init)) ec_env%should_update = ec_init
244 my_calc_forces = .false.
245 IF (
PRESENT(calculate_forces)) my_calc_forces = calculate_forces
247 IF (ec_env%should_update)
THEN
248 ec_env%old_etotal = 0.0_dp
249 ec_env%etotal = 0.0_dp
250 ec_env%eband = 0.0_dp
251 ec_env%ehartree = 0.0_dp
255 ec_env%edispersion = 0.0_dp
256 ec_env%exc_aux_fit = 0.0_dp
260 ec_env%old_etotal = energy%total
264 IF (my_calc_forces)
THEN
265 IF (unit_nr > 0)
THEN
266 WRITE (unit_nr,
'(T2,A,A,A,A,A)')
"!", repeat(
"-", 25), &
267 " Energy Correction Forces ", repeat(
"-", 26),
"!"
272 IF (unit_nr > 0)
THEN
273 WRITE (unit_nr,
'(T2,A,A,A,A,A)')
"!", repeat(
"-", 29), &
274 " Energy Correction ", repeat(
"-", 29),
"!"
279 CALL energy_correction_low(qs_env, ec_env, my_calc_forces, unit_nr)
282 IF (ec_env%should_update)
THEN
283 energy%nonscf_correction = ec_env%etotal - ec_env%old_etotal
284 energy%total = ec_env%etotal
287 IF (.NOT. my_calc_forces .AND. unit_nr > 0)
THEN
288 WRITE (unit_nr,
'(T3,A,T56,F25.15)')
"Energy Correction ", energy%nonscf_correction
290 IF (unit_nr > 0)
THEN
291 WRITE (unit_nr,
'(T2,A,A,A)')
"!", repeat(
"-", 77),
"!"
298 IF (unit_nr > 0)
THEN
299 WRITE (unit_nr,
'(T2,A,A,A)')
"!", repeat(
"-", 77),
"!"
300 WRITE (unit_nr,
'(T2,A,A,A,A,A)')
"!", repeat(
"-", 26), &
301 " Skip Energy Correction ", repeat(
"-", 27),
"!"
302 WRITE (unit_nr,
'(T2,A,A,A)')
"!", repeat(
"-", 77),
"!"
307 CALL timestop(handle)
322 SUBROUTINE energy_correction_low(qs_env, ec_env, calculate_forces, unit_nr)
323 TYPE(qs_environment_type),
POINTER :: qs_env
324 TYPE(energy_correction_type),
POINTER :: ec_env
325 LOGICAL,
INTENT(IN) :: calculate_forces
326 INTEGER,
INTENT(IN) :: unit_nr
328 INTEGER :: ispin, nspins
331 TYPE(dft_control_type),
POINTER :: dft_control
332 TYPE(pw_env_type),
POINTER :: pw_env
333 TYPE(pw_pool_type),
POINTER :: auxbas_pw_pool
335 IF (ec_env%should_update)
THEN
336 CALL ec_build_neighborlist(qs_env, ec_env)
340 ec_env%vtau_rspace, &
341 ec_env%vadmm_rspace, &
342 ec_env%ehartree, exc)
344 SELECT CASE (ec_env%energy_functional)
347 CALL ec_build_core_hamiltonian(qs_env, ec_env)
348 CALL ec_build_ks_matrix(qs_env, ec_env)
353 NULLIFY (ec_env%mao_coef)
354 CALL mao_generate_basis(qs_env, ec_env%mao_coef, ref_basis_set=
"HARRIS", molecular=.true., &
355 max_iter=ec_env%mao_max_iter, eps_grad=ec_env%mao_eps_grad, &
356 eps1_mao=ec_env%mao_eps1, iolevel=ec_env%mao_iolevel, unit_nr=unit_nr)
359 CALL ec_ks_solver(qs_env, ec_env)
361 CALL evaluate_ec_core_matrix_traces(qs_env, ec_env)
366 CALL ec_dc_energy(qs_env, ec_env, calculate_forces=.false.)
371 CALL ec_build_ks_matrix(qs_env, ec_env)
374 cpabort(
"unknown energy correction")
378 CALL ec_disp(qs_env, ec_env, calculate_forces=.false.)
381 CALL ec_energy(ec_env, unit_nr)
385 IF (calculate_forces)
THEN
387 debug_f = ec_env%debug_forces .OR. ec_env%debug_stress
389 CALL ec_disp(qs_env, ec_env, calculate_forces=.true.)
391 SELECT CASE (ec_env%energy_functional)
394 CALL ec_build_core_hamiltonian_force(qs_env, ec_env, &
398 CALL ec_build_ks_matrix_force(qs_env, ec_env)
403 CALL ec_dc_energy(qs_env, ec_env, calculate_forces=.true.)
405 CALL ec_build_core_hamiltonian_force(qs_env, ec_env, &
409 CALL ec_dc_build_ks_matrix_force(qs_env, ec_env)
412 cpabort(
"unknown energy correction")
419 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, pw_env=pw_env)
420 nspins = dft_control%nspins
422 cpassert(
ASSOCIATED(pw_env))
423 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
424 ALLOCATE (ec_env%rhoz_r(nspins))
426 CALL auxbas_pw_pool%create_pw(ec_env%rhoz_r(ispin))
430 vh_rspace=ec_env%vh_rspace, &
431 vxc_rspace=ec_env%vxc_rspace, &
432 vtau_rspace=ec_env%vtau_rspace, &
433 vadmm_rspace=ec_env%vadmm_rspace, &
434 matrix_hz=ec_env%matrix_hz, &
435 matrix_pz=ec_env%matrix_z, &
436 matrix_pz_admm=ec_env%z_admm, &
437 matrix_wz=ec_env%matrix_wz, &
438 rhopz_r=ec_env%rhoz_r, &
439 zehartree=ec_env%ehartree, &
441 zexc_aux_fit=ec_env%exc_aux_fit, &
442 p_env=ec_env%p_env, &
445 CALL ec_properties(qs_env, ec_env)
449 CALL auxbas_pw_pool%give_back_pw(ec_env%rhoout_r(ispin))
450 CALL auxbas_pw_pool%give_back_pw(ec_env%rhoz_r(ispin))
452 DEALLOCATE (ec_env%rhoout_r, ec_env%rhoz_r)
467 END SUBROUTINE energy_correction_low
476 SUBROUTINE evaluate_ec_core_matrix_traces(qs_env, ec_env)
477 TYPE(qs_environment_type),
POINTER :: qs_env
478 TYPE(energy_correction_type),
POINTER :: ec_env
480 CHARACTER(LEN=*),
PARAMETER :: routinen =
'evaluate_ec_core_matrix_traces'
483 TYPE(dft_control_type),
POINTER :: dft_control
484 TYPE(qs_energy_type),
POINTER :: energy
486 CALL timeset(routinen, handle)
489 CALL get_qs_env(qs_env, dft_control=dft_control, energy=energy)
492 CALL calculate_ptrace(ec_env%matrix_h, ec_env%matrix_p, energy%core, dft_control%nspins)
495 CALL calculate_ptrace(ec_env%matrix_t, ec_env%matrix_p, energy%kinetic, dft_control%nspins)
497 CALL timestop(handle)
499 END SUBROUTINE evaluate_ec_core_matrix_traces
512 SUBROUTINE ec_dc_energy(qs_env, ec_env, calculate_forces)
513 TYPE(qs_environment_type),
POINTER :: qs_env
514 TYPE(energy_correction_type),
POINTER :: ec_env
515 LOGICAL,
INTENT(IN) :: calculate_forces
517 CHARACTER(LEN=*),
PARAMETER :: routinen =
'ec_dc_energy'
519 CHARACTER(LEN=default_string_length) :: headline
520 INTEGER :: handle, ispin, nspins
521 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_h, matrix_p, matrix_s, matrix_w
522 TYPE(dft_control_type),
POINTER :: dft_control
523 TYPE(qs_ks_env_type),
POINTER :: ks_env
524 TYPE(qs_rho_type),
POINTER :: rho
526 CALL timeset(routinen, handle)
528 NULLIFY (dft_control, ks_env, matrix_h, matrix_p, matrix_s, matrix_w, rho)
530 dft_control=dft_control, &
532 matrix_h_kp=matrix_h, &
533 matrix_s_kp=matrix_s, &
534 matrix_w_kp=matrix_w, &
537 nspins = dft_control%nspins
542 matrix_name=
"OVERLAP MATRIX", &
543 basis_type_a=
"HARRIS", &
544 basis_type_b=
"HARRIS", &
545 sab_nl=ec_env%sab_orb)
550 headline =
"CORE HAMILTONIAN MATRIX"
551 ALLOCATE (ec_env%matrix_h(1, 1)%matrix)
552 CALL dbcsr_create(ec_env%matrix_h(1, 1)%matrix, name=trim(headline), &
553 template=matrix_s(1, 1)%matrix, matrix_type=dbcsr_type_symmetric)
555 CALL dbcsr_copy(ec_env%matrix_h(1, 1)%matrix, matrix_h(1, 1)%matrix)
560 headline =
"DENSITY MATRIX"
562 ALLOCATE (ec_env%matrix_p(ispin, 1)%matrix)
563 CALL dbcsr_create(ec_env%matrix_p(ispin, 1)%matrix, name=trim(headline), &
564 template=matrix_s(1, 1)%matrix, matrix_type=dbcsr_type_symmetric)
566 CALL dbcsr_copy(ec_env%matrix_p(ispin, 1)%matrix, matrix_p(ispin, 1)%matrix)
569 IF (calculate_forces)
THEN
574 headline =
"ENERGY-WEIGHTED DENSITY MATRIX"
576 ALLOCATE (ec_env%matrix_w(ispin, 1)%matrix)
577 CALL dbcsr_create(ec_env%matrix_w(ispin, 1)%matrix, name=trim(headline), &
578 template=matrix_s(1, 1)%matrix, matrix_type=dbcsr_type_symmetric)
580 CALL dbcsr_copy(ec_env%matrix_w(ispin, 1)%matrix, matrix_w(ispin, 1)%matrix)
586 ec_env%efield_nuclear = 0.0_dp
587 ec_env%efield_elec = 0.0_dp
590 CALL timestop(handle)
592 END SUBROUTINE ec_dc_energy
603 SUBROUTINE ec_dc_build_ks_matrix_force(qs_env, ec_env)
604 TYPE(qs_environment_type),
POINTER :: qs_env
605 TYPE(energy_correction_type),
POINTER :: ec_env
607 CHARACTER(LEN=*),
PARAMETER :: routinen =
'ec_dc_build_ks_matrix_force'
609 INTEGER :: handle, i, iounit, ispin, natom, nspins
610 LOGICAL :: debug_forces, debug_stress, do_ec_hfx, &
612 REAL(
dp) :: dummy_real, dummy_real2(2), ehartree, &
614 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: ftot
615 REAL(
dp),
DIMENSION(3) :: fodeb, fodeb2
616 REAL(kind=
dp),
DIMENSION(3, 3) :: h_stress, pv_loc, stdeb, sttot
617 TYPE(atomic_kind_type),
DIMENSION(:),
POINTER :: atomic_kind_set
618 TYPE(cell_type),
POINTER :: cell
619 TYPE(cp_logger_type),
POINTER :: logger
620 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_ks, matrix_s, scrm
621 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_p
622 TYPE(dft_control_type),
POINTER :: dft_control
623 TYPE(mp_para_env_type),
POINTER :: para_env
624 TYPE(neighbor_list_set_p_type),
DIMENSION(:), &
626 TYPE(pw_c1d_gs_type) :: rho_tot_gspace, v_hartree_gspace
627 TYPE(pw_env_type),
POINTER :: pw_env
628 TYPE(pw_grid_type),
POINTER :: pw_grid
629 TYPE(pw_poisson_type),
POINTER :: poisson_env
630 TYPE(pw_pool_type),
POINTER :: auxbas_pw_pool
631 TYPE(pw_r3d_rs_type) :: v_hartree_rspace
632 TYPE(pw_r3d_rs_type),
DIMENSION(:),
POINTER :: rho_r, v_rspace, v_rspace_in, &
634 TYPE(qs_force_type),
DIMENSION(:),
POINTER :: force
635 TYPE(qs_ks_env_type),
POINTER :: ks_env
636 TYPE(qs_rho_type),
POINTER :: rho
637 TYPE(section_vals_type),
POINTER :: ec_hfx_sections
638 TYPE(virial_type),
POINTER :: virial
640 CALL timeset(routinen, handle)
642 debug_forces = ec_env%debug_forces
643 debug_stress = ec_env%debug_stress
646 IF (logger%para_env%is_source())
THEN
652 NULLIFY (atomic_kind_set, cell, dft_control, force, ks_env, matrix_ks, &
653 matrix_p, matrix_s, para_env, pw_env, rho, sab_orb, virial)
656 dft_control=dft_control, &
659 matrix_ks=matrix_ks, &
666 cpassert(
ASSOCIATED(pw_env))
668 nspins = dft_control%nspins
669 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
671 fconv = 1.0e-9_dp*
pascal/cell%deth
672 IF (debug_stress .AND. use_virial)
THEN
673 sttot = virial%pv_virial
679 NULLIFY (auxbas_pw_pool, poisson_env)
681 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
682 poisson_env=poisson_env)
685 CALL auxbas_pw_pool%create_pw(v_hartree_gspace)
686 CALL auxbas_pw_pool%create_pw(rho_tot_gspace)
687 CALL auxbas_pw_pool%create_pw(v_hartree_rspace)
696 h_stress(:, :) = 0.0_dp
697 CALL pw_poisson_solve(poisson_env, &
698 density=rho_tot_gspace, &
700 vhartree=v_hartree_gspace, &
703 virial%pv_ehartree = virial%pv_ehartree + h_stress/real(para_env%num_pe,
dp)
704 virial%pv_virial = virial%pv_virial + h_stress/real(para_env%num_pe,
dp)
706 IF (debug_stress)
THEN
707 stdeb = fconv*(h_stress/real(para_env%num_pe,
dp))
708 CALL para_env%sum(stdeb)
709 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
714 CALL pw_poisson_solve(poisson_env, rho_tot_gspace, ehartree, &
718 CALL pw_transfer(v_hartree_gspace, v_hartree_rspace)
719 CALL pw_scale(v_hartree_rspace, v_hartree_rspace%pw_grid%dvol)
723 ALLOCATE (ec_env%rhoout_r(nspins))
725 CALL auxbas_pw_pool%create_pw(ec_env%rhoout_r(ispin))
726 CALL pw_copy(rho_r(ispin), ec_env%rhoout_r(ispin))
731 IF (debug_forces) fodeb(1:3) = force(1)%rho_core(1:3, 1)
732 IF (debug_stress .AND. use_virial) stdeb = virial%pv_ehartree
733 CALL integrate_v_core_rspace(v_hartree_rspace, qs_env)
734 IF (debug_forces)
THEN
735 fodeb(1:3) = force(1)%rho_core(1:3, 1) - fodeb(1:3)
736 CALL para_env%sum(fodeb)
737 IF (iounit > 0)
WRITE (iounit,
"(T3,A,T33,3F16.8)")
"DEBUG:: Vtot*dncore", fodeb
739 IF (debug_stress .AND. use_virial)
THEN
740 stdeb = fconv*(virial%pv_ehartree - stdeb)
741 CALL para_env%sum(stdeb)
742 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
748 NULLIFY (v_rspace, v_tau_rspace)
751 IF (use_virial) virial%pv_calculate = .true.
754 CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho, xc_section=ec_env%xc_section, &
755 vxc_rho=v_rspace, vxc_tau=v_tau_rspace, exc=exc, just_energy=.false.)
757 IF (.NOT.
ASSOCIATED(v_rspace))
THEN
758 ALLOCATE (v_rspace(nspins))
760 CALL auxbas_pw_pool%create_pw(v_rspace(ispin))
761 CALL pw_zero(v_rspace(ispin))
766 virial%pv_exc = virial%pv_exc - virial%pv_xc
767 virial%pv_virial = virial%pv_virial - virial%pv_xc
775 ALLOCATE (scrm(ispin)%matrix)
776 CALL dbcsr_create(scrm(ispin)%matrix, template=ec_env%matrix_ks(ispin, 1)%matrix)
777 CALL dbcsr_copy(scrm(ispin)%matrix, ec_env%matrix_ks(ispin, 1)%matrix)
778 CALL dbcsr_set(scrm(ispin)%matrix, 0.0_dp)
781 pw_grid => v_hartree_rspace%pw_grid
782 ALLOCATE (v_rspace_in(nspins))
784 CALL v_rspace_in(ispin)%create(pw_grid)
790 CALL pw_transfer(ec_env%vxc_rspace(ispin), v_rspace_in(ispin))
792 CALL pw_axpy(ec_env%vh_rspace, v_rspace_in(ispin))
803 IF (debug_forces) fodeb(1:3) = force(1)%fock_4c(1:3, 1)
804 IF (debug_forces) fodeb2(1:3) = force(1)%overlap_admm(1:3, 1)
811 hfx_sections=ec_hfx_sections, &
812 x_data=ec_env%x_data, &
814 do_admm=ec_env%do_ec_admm, &
815 calc_forces=.true., &
816 reuse_hfx=ec_env%reuse_hfx, &
817 do_im_time=.false., &
818 e_ex_from_gw=dummy_real, &
819 e_admm_from_gw=dummy_real2, &
822 IF (debug_forces)
THEN
823 fodeb(1:3) = force(1)%fock_4c(1:3, 1) - fodeb(1:3)
824 CALL para_env%sum(fodeb)
825 IF (iounit > 0)
WRITE (iounit,
"(T3,A,T33,3F16.8)")
"DEBUG:: P*hfx_DC ", fodeb
827 fodeb2(1:3) = force(1)%overlap_admm(1:3, 1) - fodeb2(1:3)
828 CALL para_env%sum(fodeb2)
829 IF (iounit > 0)
WRITE (iounit,
"(T3,A,T33,3F16.8)")
"DEBUG:: P*hfx_DC*S ", fodeb2
831 IF (debug_stress .AND. use_virial)
THEN
832 stdeb = -1.0_dp*fconv*virial%pv_fock_4c
833 CALL para_env%sum(stdeb)
834 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
845 pv_loc = virial%pv_virial
848 IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
849 IF (debug_stress .AND. use_virial) stdeb = virial%pv_virial
853 CALL pw_scale(v_rspace(ispin), v_rspace(ispin)%pw_grid%dvol)
854 CALL pw_axpy(v_hartree_rspace, v_rspace(ispin))
856 CALL integrate_v_rspace(v_rspace=v_rspace(ispin), &
858 pmat=matrix_p(ispin, 1), &
860 calculate_forces=.true., &
861 basis_type=
"HARRIS", &
862 task_list_external=ec_env%task_list)
865 IF (debug_forces)
THEN
866 fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
867 CALL para_env%sum(fodeb)
868 IF (iounit > 0)
WRITE (iounit,
"(T3,A,T33,3F16.8)")
"DEBUG:: Pout*dVhxc ", fodeb
870 IF (debug_stress .AND. use_virial)
THEN
871 stdeb = fconv*(virial%pv_virial - stdeb)
872 CALL para_env%sum(stdeb)
873 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
877 IF (
ASSOCIATED(v_tau_rspace))
THEN
878 IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
879 IF (debug_stress .AND. use_virial) stdeb = virial%pv_virial
881 CALL pw_scale(v_tau_rspace(ispin), v_tau_rspace(ispin)%pw_grid%dvol)
883 CALL integrate_v_rspace(v_rspace=v_tau_rspace(ispin), &
885 pmat=matrix_p(ispin, 1), &
887 calculate_forces=.true., &
888 compute_tau=.true., &
889 basis_type=
"HARRIS", &
890 task_list_external=ec_env%task_list)
893 IF (debug_forces)
THEN
894 fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
895 CALL para_env%sum(fodeb)
896 IF (iounit > 0)
WRITE (iounit,
"(T3,A,T33,3F16.8)")
"DEBUG:: Pout*dVhxc_tau ", fodeb
898 IF (debug_stress .AND. use_virial)
THEN
899 stdeb = fconv*(virial%pv_virial - stdeb)
900 CALL para_env%sum(stdeb)
901 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
908 virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
927 NULLIFY (ec_env%matrix_hz)
930 ALLOCATE (ec_env%matrix_hz(ispin)%matrix)
931 CALL dbcsr_create(ec_env%matrix_hz(ispin)%matrix, template=matrix_s(1)%matrix)
932 CALL dbcsr_copy(ec_env%matrix_hz(ispin)%matrix, matrix_s(1)%matrix)
933 CALL dbcsr_set(ec_env%matrix_hz(ispin)%matrix, 0.0_dp)
939 CALL pw_axpy(v_rspace_in(ispin), v_rspace(ispin), -1.0_dp)
943 CALL integrate_v_rspace(v_rspace=v_rspace(ispin), &
944 hmat=ec_env%matrix_hz(ispin), &
945 pmat=matrix_p(ispin, 1), &
947 calculate_forces=.false., &
948 basis_type=
"HARRIS", &
949 task_list_external=ec_env%task_list)
953 IF (dft_control%use_kinetic_energy_density)
THEN
956 IF (.NOT.
ASSOCIATED(v_tau_rspace))
THEN
957 ALLOCATE (v_tau_rspace(nspins))
959 CALL auxbas_pw_pool%create_pw(v_tau_rspace(ispin))
960 CALL pw_zero(v_tau_rspace(ispin))
966 IF (
ASSOCIATED(ec_env%vtau_rspace))
THEN
967 CALL pw_axpy(ec_env%vtau_rspace(ispin), v_tau_rspace(ispin), -1.0_dp)
970 CALL integrate_v_rspace(v_rspace=v_tau_rspace(ispin), &
971 hmat=ec_env%matrix_hz(ispin), &
972 pmat=matrix_p(ispin, 1), &
974 calculate_forces=.false., compute_tau=.true., &
975 basis_type=
"HARRIS", &
976 task_list_external=ec_env%task_list)
984 ext_hfx_section=ec_hfx_sections, &
985 x_data=ec_env%x_data, &
986 recalc_integrals=.false., &
987 do_admm=ec_env%do_ec_admm, &
990 reuse_hfx=ec_env%reuse_hfx)
993 IF (debug_forces) fodeb(1:3) = force(1)%core_overlap(1:3, 1)
994 IF (debug_stress .AND. use_virial) stdeb = virial%pv_ecore_overlap
996 IF (debug_forces)
THEN
997 fodeb(1:3) = force(1)%core_overlap(1:3, 1) - fodeb(1:3)
998 CALL para_env%sum(fodeb)
999 IF (iounit > 0)
WRITE (iounit,
"(T3,A,T33,3F16.8)")
"DEBUG:: CoreOverlap", fodeb
1001 IF (debug_stress .AND. use_virial)
THEN
1002 stdeb = fconv*(stdeb - virial%pv_ecore_overlap)
1003 CALL para_env%sum(stdeb)
1004 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
1008 IF (debug_forces)
THEN
1009 CALL get_qs_env(qs_env, natom=natom, atomic_kind_set=atomic_kind_set)
1010 ALLOCATE (ftot(3, natom))
1012 fodeb(1:3) = ftot(1:3, 1)
1014 CALL para_env%sum(fodeb)
1015 IF (iounit > 0)
WRITE (iounit,
"(T3,A,T33,3F16.8)")
"DEBUG:: Force Explicit", fodeb
1019 DO ispin = 1, nspins
1020 CALL auxbas_pw_pool%give_back_pw(v_rspace(ispin))
1021 CALL auxbas_pw_pool%give_back_pw(v_rspace_in(ispin))
1022 IF (
ASSOCIATED(v_tau_rspace))
THEN
1023 CALL auxbas_pw_pool%give_back_pw(v_tau_rspace(ispin))
1027 DEALLOCATE (v_rspace, v_rspace_in)
1028 IF (
ASSOCIATED(v_tau_rspace))
DEALLOCATE (v_tau_rspace)
1030 CALL auxbas_pw_pool%give_back_pw(rho_tot_gspace)
1031 CALL auxbas_pw_pool%give_back_pw(v_hartree_gspace)
1032 CALL auxbas_pw_pool%give_back_pw(v_hartree_rspace)
1036 IF (use_virial)
THEN
1037 IF (qs_env%energy_correction)
THEN
1038 ec_env%ehartree = ehartree
1043 IF (debug_stress .AND. use_virial)
THEN
1045 stdeb = -1.0_dp*fconv*ehartree
1046 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
1049 stdeb = -1.0_dp*fconv*exc
1050 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
1056 TYPE(virial_type) :: virdeb
1059 CALL para_env%sum(virdeb%pv_overlap)
1060 CALL para_env%sum(virdeb%pv_ekinetic)
1061 CALL para_env%sum(virdeb%pv_ppl)
1062 CALL para_env%sum(virdeb%pv_ppnl)
1063 CALL para_env%sum(virdeb%pv_ecore_overlap)
1064 CALL para_env%sum(virdeb%pv_ehartree)
1065 CALL para_env%sum(virdeb%pv_exc)
1066 CALL para_env%sum(virdeb%pv_exx)
1067 CALL para_env%sum(virdeb%pv_vdw)
1068 CALL para_env%sum(virdeb%pv_mp2)
1069 CALL para_env%sum(virdeb%pv_nlcc)
1070 CALL para_env%sum(virdeb%pv_gapw)
1071 CALL para_env%sum(virdeb%pv_lrigpw)
1072 CALL para_env%sum(virdeb%pv_virial)
1077 virdeb%pv_ehartree(i, i) = virdeb%pv_ehartree(i, i) - 2.0_dp*ehartree
1078 virdeb%pv_virial(i, i) = virdeb%pv_virial(i, i) - exc &
1080 virdeb%pv_exc(i, i) = virdeb%pv_exc(i, i) - exc
1086 CALL para_env%sum(sttot)
1087 stdeb = fconv*(virdeb%pv_virial - sttot)
1088 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
1091 stdeb = fconv*(virdeb%pv_virial)
1092 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
1101 CALL timestop(handle)
1103 END SUBROUTINE ec_dc_build_ks_matrix_force
1111 SUBROUTINE ec_disp(qs_env, ec_env, calculate_forces)
1112 TYPE(qs_environment_type),
POINTER :: qs_env
1113 TYPE(energy_correction_type),
POINTER :: ec_env
1114 LOGICAL,
INTENT(IN) :: calculate_forces
1116 REAL(kind=dp) :: edisp
1118 CALL calculate_dispersion_pairpot(qs_env, ec_env%dispersion_env, edisp, calculate_forces)
1119 IF (.NOT. calculate_forces) ec_env%edispersion = ec_env%edispersion + edisp
1121 END SUBROUTINE ec_disp
1130 SUBROUTINE ec_build_core_hamiltonian(qs_env, ec_env)
1131 TYPE(qs_environment_type),
POINTER :: qs_env
1132 TYPE(energy_correction_type),
POINTER :: ec_env
1134 CHARACTER(LEN=*),
PARAMETER :: routinen =
'ec_build_core_hamiltonian'
1136 INTEGER :: handle, nder, nimages
1137 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
1138 LOGICAL :: calculate_forces, use_virial
1139 REAL(kind=dp) :: eps_ppnl
1140 TYPE(atomic_kind_type),
DIMENSION(:),
POINTER :: atomic_kind_set
1141 TYPE(dft_control_type),
POINTER :: dft_control
1142 TYPE(neighbor_list_set_p_type),
DIMENSION(:), &
1143 POINTER :: sab_orb, sac_ae, sac_ppl, sap_ppnl
1144 TYPE(particle_type),
DIMENSION(:),
POINTER :: particle_set
1145 TYPE(qs_force_type),
DIMENSION(:),
POINTER :: force
1146 TYPE(qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
1147 TYPE(qs_ks_env_type),
POINTER :: ks_env
1148 TYPE(virial_type),
POINTER :: virial
1150 CALL timeset(routinen, handle)
1152 NULLIFY (atomic_kind_set, cell_to_index, dft_control, ks_env, particle_set, qs_kind_set, virial)
1154 CALL get_qs_env(qs_env=qs_env, &
1155 atomic_kind_set=atomic_kind_set, &
1156 dft_control=dft_control, &
1157 particle_set=particle_set, &
1158 qs_kind_set=qs_kind_set, &
1162 nimages = dft_control%nimages
1163 IF (nimages /= 1)
THEN
1164 cpabort(
"K-points for Harris functional not implemented")
1168 IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc)
THEN
1169 cpabort(
"Harris functional for GAPW not implemented")
1173 use_virial = .false.
1174 calculate_forces = .false.
1177 NULLIFY (sab_orb, sac_ae, sac_ppl, sap_ppnl)
1178 sab_orb => ec_env%sab_orb
1179 sac_ppl => ec_env%sac_ppl
1180 sap_ppnl => ec_env%sap_ppnl
1184 CALL build_overlap_matrix(ks_env, matrixkp_s=ec_env%matrix_s, &
1185 matrix_name=
"OVERLAP MATRIX", &
1186 basis_type_a=
"HARRIS", &
1187 basis_type_b=
"HARRIS", &
1189 CALL build_kinetic_matrix(ks_env, matrixkp_t=ec_env%matrix_t, &
1190 matrix_name=
"KINETIC ENERGY MATRIX", &
1191 basis_type=
"HARRIS", &
1195 CALL dbcsr_allocate_matrix_set(ec_env%matrix_h, 1, 1)
1196 ALLOCATE (ec_env%matrix_h(1, 1)%matrix)
1197 CALL dbcsr_create(ec_env%matrix_h(1, 1)%matrix, template=ec_env%matrix_s(1, 1)%matrix)
1198 CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_h(1, 1)%matrix, sab_orb)
1201 CALL dbcsr_copy(ec_env%matrix_h(1, 1)%matrix, ec_env%matrix_t(1, 1)%matrix, &
1202 keep_sparsity=.true., name=
"CORE HAMILTONIAN MATRIX")
1205 IF (
ASSOCIATED(sac_ae))
THEN
1206 cpabort(
"ECP/AE not available for energy corrections")
1209 CALL build_core_ae(ec_env%matrix_h, ec_env%matrix_p, force, &
1210 virial, calculate_forces, use_virial, nder, &
1211 qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ae, &
1212 nimages, cell_to_index)
1215 IF (
ASSOCIATED(sac_ppl))
THEN
1216 CALL build_core_ppl(ec_env%matrix_h, ec_env%matrix_p, force, &
1217 virial, calculate_forces, use_virial, nder, &
1218 qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ppl, &
1219 nimages, cell_to_index,
"HARRIS")
1223 eps_ppnl = dft_control%qs_control%eps_ppnl
1224 IF (
ASSOCIATED(sap_ppnl))
THEN
1225 CALL build_core_ppnl(ec_env%matrix_h, ec_env%matrix_p, force, &
1226 virial, calculate_forces, use_virial, nder, &
1227 qs_kind_set, atomic_kind_set, particle_set, sab_orb, sap_ppnl, &
1228 eps_ppnl, nimages, cell_to_index,
"HARRIS")
1232 ec_env%efield_nuclear = 0.0_dp
1233 CALL ec_efield_local_operator(qs_env, ec_env, calculate_forces)
1235 CALL timestop(handle)
1237 END SUBROUTINE ec_build_core_hamiltonian
1248 SUBROUTINE ec_build_ks_matrix(qs_env, ec_env)
1249 TYPE(qs_environment_type),
POINTER :: qs_env
1250 TYPE(energy_correction_type),
POINTER :: ec_env
1252 CHARACTER(LEN=*),
PARAMETER :: routinen =
'ec_build_ks_matrix'
1254 CHARACTER(LEN=default_string_length) :: headline
1255 INTEGER :: handle, iounit, ispin, nspins
1256 LOGICAL :: calculate_forces, &
1257 do_adiabatic_rescaling, do_ec_hfx, &
1258 hfx_treat_lsd_in_core, use_virial
1259 REAL(dp) :: dummy_real, dummy_real2(2), eexc, evhxc, &
1261 TYPE(cp_logger_type),
POINTER :: logger
1262 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: ks_mat
1263 TYPE(dft_control_type),
POINTER :: dft_control
1264 TYPE(pw_env_type),
POINTER :: pw_env
1265 TYPE(pw_pool_type),
POINTER :: auxbas_pw_pool
1266 TYPE(pw_r3d_rs_type),
DIMENSION(:),
POINTER :: rho_r, tau_r, v_rspace, v_tau_rspace
1267 TYPE(qs_energy_type),
POINTER :: energy
1268 TYPE(qs_ks_env_type),
POINTER :: ks_env
1269 TYPE(qs_rho_type),
POINTER :: rho
1270 TYPE(section_vals_type),
POINTER :: adiabatic_rescaling_section, &
1271 ec_hfx_sections, ec_section
1273 CALL timeset(routinen, handle)
1275 logger => cp_get_default_logger()
1276 IF (logger%para_env%is_source())
THEN
1277 iounit = cp_logger_get_default_unit_nr(logger, local=.true.)
1283 NULLIFY (auxbas_pw_pool, dft_control, energy, ks_env, rho, rho_r, tau_r)
1284 CALL get_qs_env(qs_env=qs_env, &
1285 dft_control=dft_control, &
1288 nspins = dft_control%nspins
1289 calculate_forces = .false.
1290 use_virial = .false.
1293 IF (
ASSOCIATED(ec_env%matrix_ks))
CALL dbcsr_deallocate_matrix_set(ec_env%matrix_ks)
1294 CALL dbcsr_allocate_matrix_set(ec_env%matrix_ks, nspins, 1)
1295 DO ispin = 1, nspins
1296 headline =
"KOHN-SHAM MATRIX"
1297 ALLOCATE (ec_env%matrix_ks(ispin, 1)%matrix)
1298 CALL dbcsr_create(ec_env%matrix_ks(ispin, 1)%matrix, name=trim(headline), &
1299 template=ec_env%matrix_s(1, 1)%matrix, matrix_type=dbcsr_type_symmetric)
1300 CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_ks(ispin, 1)%matrix, ec_env%sab_orb)
1301 CALL dbcsr_set(ec_env%matrix_ks(ispin, 1)%matrix, 0.0_dp)
1305 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
1306 cpassert(
ASSOCIATED(pw_env))
1309 ec_section => section_vals_get_subs_vals(qs_env%input,
"DFT%ENERGY_CORRECTION")
1310 ec_hfx_sections => section_vals_get_subs_vals(ec_section,
"XC%HF")
1311 CALL section_vals_get(ec_hfx_sections, explicit=do_ec_hfx)
1316 adiabatic_rescaling_section => section_vals_get_subs_vals(ec_section,
"XC%ADIABATIC_RESCALING")
1317 CALL section_vals_get(adiabatic_rescaling_section, explicit=do_adiabatic_rescaling)
1318 IF (do_adiabatic_rescaling)
THEN
1319 CALL cp_abort(__location__,
"Adiabatic rescaling NYI for energy correction")
1321 CALL section_vals_val_get(ec_hfx_sections,
"TREAT_LSD_IN_CORE", l_val=hfx_treat_lsd_in_core)
1322 IF (hfx_treat_lsd_in_core)
THEN
1323 CALL cp_abort(__location__,
"HFX_TREAT_LSD_IN_CORE NYI for energy correction")
1327 IF (dft_control%do_admm)
THEN
1329 IF (dft_control%do_admm_mo)
THEN
1330 IF (qs_env%run_rtp)
THEN
1331 CALL rtp_admm_calc_rho_aux(qs_env)
1333 CALL admm_mo_calc_rho_aux(qs_env)
1335 ELSEIF (dft_control%do_admm_dm)
THEN
1336 CALL admm_dm_calc_rho_aux(qs_env)
1343 CALL get_qs_env(qs_env, energy=energy)
1344 CALL calculate_exx(qs_env=qs_env, &
1346 hfx_sections=ec_hfx_sections, &
1347 x_data=ec_env%x_data, &
1349 do_admm=ec_env%do_ec_admm, &
1350 calc_forces=.false., &
1351 reuse_hfx=ec_env%reuse_hfx, &
1352 do_im_time=.false., &
1353 e_ex_from_gw=dummy_real, &
1354 e_admm_from_gw=dummy_real2, &
1358 ec_env%ex = energy%ex
1360 IF (ec_env%do_ec_admm)
THEN
1361 ec_env%exc_aux_fit = energy%exc_aux_fit + energy%exc
1367 ks_mat => ec_env%matrix_ks(:, 1)
1368 CALL add_exx_to_rhs(rhs=ks_mat, &
1370 ext_hfx_section=ec_hfx_sections, &
1371 x_data=ec_env%x_data, &
1372 recalc_integrals=.false., &
1373 do_admm=ec_env%do_ec_admm, &
1376 reuse_hfx=ec_env%reuse_hfx)
1381 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
1382 NULLIFY (v_rspace, v_tau_rspace)
1383 CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho, xc_section=ec_env%xc_section, &
1384 vxc_rho=v_rspace, vxc_tau=v_tau_rspace, exc=eexc, just_energy=.false.)
1386 IF (.NOT.
ASSOCIATED(v_rspace))
THEN
1387 ALLOCATE (v_rspace(nspins))
1388 DO ispin = 1, nspins
1389 CALL auxbas_pw_pool%create_pw(v_rspace(ispin))
1390 CALL pw_zero(v_rspace(ispin))
1395 CALL qs_rho_get(rho, rho_r=rho_r)
1396 IF (
ASSOCIATED(v_tau_rspace))
THEN
1397 CALL qs_rho_get(rho, tau_r=tau_r)
1399 DO ispin = 1, nspins
1401 CALL pw_scale(v_rspace(ispin), v_rspace(ispin)%pw_grid%dvol)
1402 CALL pw_axpy(ec_env%vh_rspace, v_rspace(ispin))
1404 CALL integrate_v_rspace(v_rspace=v_rspace(ispin), &
1405 hmat=ec_env%matrix_ks(ispin, 1), &
1407 calculate_forces=.false., &
1408 basis_type=
"HARRIS", &
1409 task_list_external=ec_env%task_list)
1411 IF (
ASSOCIATED(v_tau_rspace))
THEN
1413 CALL pw_scale(v_tau_rspace(ispin), v_tau_rspace(ispin)%pw_grid%dvol)
1414 CALL integrate_v_rspace(v_rspace=v_tau_rspace(ispin), &
1415 hmat=ec_env%matrix_ks(ispin, 1), &
1417 calculate_forces=.false., &
1418 compute_tau=.true., &
1419 basis_type=
"HARRIS", &
1420 task_list_external=ec_env%task_list)
1424 evhxc = evhxc + pw_integral_ab(rho_r(ispin), v_rspace(ispin))/ &
1425 v_rspace(1)%pw_grid%dvol
1426 IF (
ASSOCIATED(v_tau_rspace))
THEN
1427 evhxc = evhxc + pw_integral_ab(tau_r(ispin), v_tau_rspace(ispin))/ &
1428 v_tau_rspace(ispin)%pw_grid%dvol
1434 DO ispin = 1, nspins
1435 CALL auxbas_pw_pool%give_back_pw(v_rspace(ispin))
1436 IF (
ASSOCIATED(v_tau_rspace))
THEN
1437 CALL auxbas_pw_pool%give_back_pw(v_tau_rspace(ispin))
1440 DEALLOCATE (v_rspace)
1441 IF (
ASSOCIATED(v_tau_rspace))
DEALLOCATE (v_tau_rspace)
1448 DO ispin = 1, nspins
1449 CALL dbcsr_add(ec_env%matrix_ks(ispin, 1)%matrix, ec_env%matrix_h(1, 1)%matrix, &
1450 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
1451 CALL dbcsr_filter(ec_env%matrix_ks(ispin, 1)%matrix, &
1452 dft_control%qs_control%eps_filter_matrix)
1455 CALL timestop(handle)
1457 END SUBROUTINE ec_build_ks_matrix
1469 SUBROUTINE ec_build_core_hamiltonian_force(qs_env, ec_env, matrix_p, matrix_s, matrix_w)
1470 TYPE(qs_environment_type),
POINTER :: qs_env
1471 TYPE(energy_correction_type),
POINTER :: ec_env
1472 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_p, matrix_s, matrix_w
1474 CHARACTER(LEN=*),
PARAMETER :: routinen =
'ec_build_core_hamiltonian_force'
1476 INTEGER :: handle, iounit, nder, nimages
1477 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
1478 LOGICAL :: calculate_forces, debug_forces, &
1479 debug_stress, use_virial
1480 REAL(kind=dp) :: eps_ppnl, fconv
1481 REAL(kind=dp),
DIMENSION(3) :: fodeb
1482 REAL(kind=dp),
DIMENSION(3, 3) :: stdeb, sttot
1483 TYPE(atomic_kind_type),
DIMENSION(:),
POINTER :: atomic_kind_set
1484 TYPE(cell_type),
POINTER :: cell
1485 TYPE(cp_logger_type),
POINTER :: logger
1486 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: scrm
1487 TYPE(dft_control_type),
POINTER :: dft_control
1488 TYPE(mp_para_env_type),
POINTER :: para_env
1489 TYPE(neighbor_list_set_p_type),
DIMENSION(:), &
1490 POINTER :: sab_orb, sac_ppl, sap_ppnl
1491 TYPE(particle_type),
DIMENSION(:),
POINTER :: particle_set
1492 TYPE(qs_force_type),
DIMENSION(:),
POINTER :: force
1493 TYPE(qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
1494 TYPE(qs_ks_env_type),
POINTER :: ks_env
1495 TYPE(virial_type),
POINTER :: virial
1497 CALL timeset(routinen, handle)
1499 debug_forces = ec_env%debug_forces
1500 debug_stress = ec_env%debug_stress
1502 logger => cp_get_default_logger()
1503 IF (logger%para_env%is_source())
THEN
1504 iounit = cp_logger_get_default_unit_nr(logger, local=.true.)
1509 calculate_forces = .true.
1512 NULLIFY (cell, dft_control, force, ks_env, para_env, virial)
1513 CALL get_qs_env(qs_env=qs_env, &
1515 dft_control=dft_control, &
1518 para_env=para_env, &
1520 nimages = dft_control%nimages
1521 IF (nimages /= 1)
THEN
1522 cpabort(
"K-points for Harris functional not implemented")
1526 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
1528 fconv = 1.0e-9_dp*pascal/cell%deth
1529 IF (debug_stress .AND. use_virial)
THEN
1530 sttot = virial%pv_virial
1534 IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc)
THEN
1535 cpabort(
"Harris functional for GAPW not implemented")
1539 NULLIFY (sab_orb, sac_ppl, sap_ppnl)
1540 sab_orb => ec_env%sab_orb
1541 sac_ppl => ec_env%sac_ppl
1542 sap_ppnl => ec_env%sap_ppnl
1546 CALL dbcsr_allocate_matrix_set(scrm, 1, 1)
1547 ALLOCATE (scrm(1, 1)%matrix)
1548 CALL dbcsr_create(scrm(1, 1)%matrix, template=matrix_s(1, 1)%matrix)
1549 CALL cp_dbcsr_alloc_block_from_nbl(scrm(1, 1)%matrix, sab_orb)
1552 IF (
SIZE(matrix_p, 1) == 2)
THEN
1553 CALL dbcsr_add(matrix_p(1, 1)%matrix, matrix_p(2, 1)%matrix, &
1554 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
1555 CALL dbcsr_add(matrix_w(1, 1)%matrix, matrix_w(2, 1)%matrix, &
1556 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
1560 IF (debug_forces) fodeb(1:3) = force(1)%overlap(1:3, 1)
1561 IF (debug_stress .AND. use_virial) stdeb = virial%pv_overlap
1562 CALL build_overlap_matrix(ks_env, matrixkp_s=scrm, &
1563 matrix_name=
"OVERLAP MATRIX", &
1564 basis_type_a=
"HARRIS", &
1565 basis_type_b=
"HARRIS", &
1566 sab_nl=sab_orb, calculate_forces=.true., &
1567 matrixkp_p=matrix_w)
1569 IF (debug_forces)
THEN
1570 fodeb(1:3) = force(1)%overlap(1:3, 1) - fodeb(1:3)
1571 CALL para_env%sum(fodeb)
1572 IF (iounit > 0)
WRITE (iounit,
"(T3,A,T33,3F16.8)")
"DEBUG:: Wout*dS ", fodeb
1573 fodeb(1:3) = force(1)%kinetic(1:3, 1)
1575 IF (debug_stress .AND. use_virial)
THEN
1576 stdeb = fconv*(virial%pv_overlap - stdeb)
1577 CALL para_env%sum(stdeb)
1578 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
1579 'STRESS| Wout*dS', one_third_sum_diag(stdeb),
det_3x3(stdeb)
1580 stdeb = virial%pv_ekinetic
1582 CALL build_kinetic_matrix(ks_env, matrixkp_t=scrm, &
1583 matrix_name=
"KINETIC ENERGY MATRIX", &
1584 basis_type=
"HARRIS", &
1585 sab_nl=sab_orb, calculate_forces=.true., &
1586 matrixkp_p=matrix_p)
1587 IF (debug_forces)
THEN
1588 fodeb(1:3) = force(1)%kinetic(1:3, 1) - fodeb(1:3)
1589 CALL para_env%sum(fodeb)
1590 IF (iounit > 0)
WRITE (iounit,
"(T3,A,T33,3F16.8)")
"DEBUG:: Pout*dT ", fodeb
1592 IF (debug_stress .AND. use_virial)
THEN
1593 stdeb = fconv*(virial%pv_ekinetic - stdeb)
1594 CALL para_env%sum(stdeb)
1595 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
1596 'STRESS| Pout*dT', one_third_sum_diag(stdeb),
det_3x3(stdeb)
1598 IF (
SIZE(matrix_p, 1) == 2)
THEN
1599 CALL dbcsr_add(matrix_p(1, 1)%matrix, matrix_p(2, 1)%matrix, &
1600 alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
1604 NULLIFY (atomic_kind_set, particle_set, qs_kind_set)
1605 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, particle_set=particle_set, &
1606 atomic_kind_set=atomic_kind_set)
1608 IF (
ASSOCIATED(sac_ppl))
THEN
1609 IF (calculate_forces .AND. debug_forces) fodeb(1:3) = force(1)%gth_ppl(1:3, 1)
1610 IF (debug_stress .AND. use_virial) stdeb = virial%pv_ppl
1611 CALL build_core_ppl(scrm, matrix_p, force, &
1612 virial, calculate_forces, use_virial, nder, &
1613 qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ppl, &
1614 nimages, cell_to_index,
"HARRIS")
1615 IF (calculate_forces .AND. debug_forces)
THEN
1616 fodeb(1:3) = force(1)%gth_ppl(1:3, 1) - fodeb(1:3)
1617 CALL para_env%sum(fodeb)
1618 IF (iounit > 0)
WRITE (iounit,
"(T3,A,T33,3F16.8)")
"DEBUG:: Pout*dH_PPL ", fodeb
1620 IF (debug_stress .AND. use_virial)
THEN
1621 stdeb = fconv*(virial%pv_ppl - stdeb)
1622 CALL para_env%sum(stdeb)
1623 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
1624 'STRESS| Pout*dH_PPL', one_third_sum_diag(stdeb),
det_3x3(stdeb)
1629 eps_ppnl = dft_control%qs_control%eps_ppnl
1630 IF (
ASSOCIATED(sap_ppnl))
THEN
1631 IF (calculate_forces .AND. debug_forces) fodeb(1:3) = force(1)%gth_ppnl(1:3, 1)
1632 IF (debug_stress .AND. use_virial) stdeb = virial%pv_ppnl
1633 CALL build_core_ppnl(scrm, matrix_p, force, &
1634 virial, calculate_forces, use_virial, nder, &
1635 qs_kind_set, atomic_kind_set, particle_set, &
1636 sab_orb, sap_ppnl, eps_ppnl, &
1637 nimages, cell_to_index,
"HARRIS")
1638 IF (calculate_forces .AND. debug_forces)
THEN
1639 fodeb(1:3) = force(1)%gth_ppnl(1:3, 1) - fodeb(1:3)
1640 CALL para_env%sum(fodeb)
1641 IF (iounit > 0)
WRITE (iounit,
"(T3,A,T33,3F16.8)")
"DEBUG:: Pout*dH_PPNL", fodeb
1643 IF (debug_stress .AND. use_virial)
THEN
1644 stdeb = fconv*(virial%pv_ppnl - stdeb)
1645 CALL para_env%sum(stdeb)
1646 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
1647 'STRESS| Pout*dH_PPNL', one_third_sum_diag(stdeb),
det_3x3(stdeb)
1652 ec_env%efield_nuclear = 0.0_dp
1653 IF (calculate_forces .AND. debug_forces) fodeb(1:3) = force(1)%efield(1:3, 1)
1654 CALL ec_efield_local_operator(qs_env, ec_env, calculate_forces)
1655 IF (calculate_forces .AND. debug_forces)
THEN
1656 fodeb(1:3) = force(1)%efield(1:3, 1) - fodeb(1:3)
1657 CALL para_env%sum(fodeb)
1658 IF (iounit > 0)
WRITE (iounit,
"(T3,A,T33,3F16.8)")
"DEBUG:: Pout*dEfield", fodeb
1660 IF (debug_stress .AND. use_virial)
THEN
1661 stdeb = fconv*(virial%pv_virial - sttot)
1662 CALL para_env%sum(stdeb)
1663 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
1664 'STRESS| Stress Pout*dHcore ', one_third_sum_diag(stdeb),
det_3x3(stdeb)
1665 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))")
' '
1669 CALL dbcsr_deallocate_matrix_set(scrm)
1671 CALL timestop(handle)
1673 END SUBROUTINE ec_build_core_hamiltonian_force
1684 SUBROUTINE ec_build_ks_matrix_force(qs_env, ec_env)
1685 TYPE(qs_environment_type),
POINTER :: qs_env
1686 TYPE(energy_correction_type),
POINTER :: ec_env
1688 CHARACTER(LEN=*),
PARAMETER :: routinen =
'ec_build_ks_matrix_force'
1690 INTEGER :: handle, i, iounit, ispin, natom, nspins
1691 LOGICAL :: debug_forces, debug_stress, do_ec_hfx, &
1693 REAL(dp) :: dehartree, dummy_real, dummy_real2(2), &
1694 eexc, ehartree, eovrl, exc, fconv
1695 REAL(dp),
ALLOCATABLE,
DIMENSION(:, :) :: ftot
1696 REAL(dp),
DIMENSION(3) :: fodeb
1697 REAL(kind=dp),
DIMENSION(3, 3) :: h_stress, pv_loc, stdeb, sttot
1698 TYPE(atomic_kind_type),
DIMENSION(:),
POINTER :: atomic_kind_set
1699 TYPE(cell_type),
POINTER :: cell
1700 TYPE(cp_logger_type),
POINTER :: logger
1701 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_ks, scrm
1702 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_p, matrix_s
1703 TYPE(dft_control_type),
POINTER :: dft_control
1704 TYPE(mp_para_env_type),
POINTER :: para_env
1705 TYPE(neighbor_list_set_p_type),
DIMENSION(:), &
1707 TYPE(pw_c1d_gs_type) :: rho_tot_gspace, rhodn_tot_gspace, &
1709 TYPE(pw_c1d_gs_type),
DIMENSION(:),
POINTER :: rho_g, rhoout_g
1710 TYPE(pw_c1d_gs_type),
POINTER :: rho_core
1711 TYPE(pw_env_type),
POINTER :: pw_env
1712 TYPE(pw_poisson_type),
POINTER :: poisson_env
1713 TYPE(pw_pool_type),
POINTER :: auxbas_pw_pool
1714 TYPE(pw_r3d_rs_type) :: dv_hartree_rspace, v_hartree_rspace, &
1716 TYPE(pw_r3d_rs_type),
DIMENSION(:),
POINTER :: rho_r, rhoout_r, tau_r, tauout_r, &
1717 v_rspace, v_tau_rspace, v_xc, v_xc_tau
1718 TYPE(qs_force_type),
DIMENSION(:),
POINTER :: force
1719 TYPE(qs_ks_env_type),
POINTER :: ks_env
1720 TYPE(qs_rho_type),
POINTER :: rho
1721 TYPE(section_vals_type),
POINTER :: ec_hfx_sections, xc_section
1722 TYPE(virial_type),
POINTER :: virial
1724 CALL timeset(routinen, handle)
1726 debug_forces = ec_env%debug_forces
1727 debug_stress = ec_env%debug_stress
1729 logger => cp_get_default_logger()
1730 IF (logger%para_env%is_source())
THEN
1731 iounit = cp_logger_get_default_unit_nr(logger, local=.true.)
1737 NULLIFY (atomic_kind_set, cell, dft_control, force, ks_env, &
1738 matrix_ks, matrix_p, matrix_s, para_env, rho, rho_core, &
1739 rho_g, rho_r, sab_orb, tau_r, virial)
1740 CALL get_qs_env(qs_env=qs_env, &
1742 dft_control=dft_control, &
1745 matrix_ks=matrix_ks, &
1746 para_env=para_env, &
1751 nspins = dft_control%nspins
1752 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
1754 fconv = 1.0e-9_dp*pascal/cell%deth
1755 IF (debug_stress .AND. use_virial)
THEN
1756 sttot = virial%pv_virial
1760 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
1761 cpassert(
ASSOCIATED(pw_env))
1763 NULLIFY (auxbas_pw_pool, poisson_env)
1765 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
1766 poisson_env=poisson_env)
1769 CALL auxbas_pw_pool%create_pw(v_hartree_gspace)
1770 CALL auxbas_pw_pool%create_pw(rhodn_tot_gspace)
1771 CALL auxbas_pw_pool%create_pw(v_hartree_rspace)
1773 CALL pw_transfer(ec_env%vh_rspace, v_hartree_rspace)
1778 CALL qs_rho_get(rho, rho_r=rho_r, rho_g=rho_g, tau_r=tau_r)
1779 NULLIFY (rhoout_r, rhoout_g)
1780 ALLOCATE (rhoout_r(nspins), rhoout_g(nspins))
1781 DO ispin = 1, nspins
1782 CALL auxbas_pw_pool%create_pw(rhoout_r(ispin))
1783 CALL auxbas_pw_pool%create_pw(rhoout_g(ispin))
1785 CALL auxbas_pw_pool%create_pw(dv_hartree_rspace)
1786 CALL auxbas_pw_pool%create_pw(vtot_rspace)
1788 CALL pw_zero(rhodn_tot_gspace)
1789 DO ispin = 1, nspins
1790 CALL calculate_rho_elec(ks_env=ks_env, matrix_p=ec_env%matrix_p(ispin, 1)%matrix, &
1791 rho=rhoout_r(ispin), &
1792 rho_gspace=rhoout_g(ispin), &
1793 basis_type=
"HARRIS", &
1794 task_list_external=ec_env%task_list)
1798 ALLOCATE (ec_env%rhoout_r(nspins))
1799 DO ispin = 1, nspins
1800 CALL auxbas_pw_pool%create_pw(ec_env%rhoout_r(ispin))
1801 CALL pw_copy(rhoout_r(ispin), ec_env%rhoout_r(ispin))
1805 IF (dft_control%use_kinetic_energy_density)
THEN
1807 TYPE(pw_c1d_gs_type) :: tauout_g
1808 ALLOCATE (tauout_r(nspins))
1809 DO ispin = 1, nspins
1810 CALL auxbas_pw_pool%create_pw(tauout_r(ispin))
1812 CALL auxbas_pw_pool%create_pw(tauout_g)
1814 DO ispin = 1, nspins
1815 CALL calculate_rho_elec(ks_env=ks_env, matrix_p=ec_env%matrix_p(ispin, 1)%matrix, &
1816 rho=tauout_r(ispin), &
1817 rho_gspace=tauout_g, &
1818 compute_tau=.true., &
1819 basis_type=
"HARRIS", &
1820 task_list_external=ec_env%task_list)
1823 CALL auxbas_pw_pool%give_back_pw(tauout_g)
1827 IF (use_virial)
THEN
1830 CALL auxbas_pw_pool%create_pw(rho_tot_gspace)
1833 CALL calc_rho_tot_gspace(rho_tot_gspace, qs_env, rho)
1836 CALL get_qs_env(qs_env=qs_env, rho_core=rho_core)
1837 CALL pw_copy(rho_core, rhodn_tot_gspace)
1838 DO ispin = 1, dft_control%nspins
1839 CALL pw_axpy(rhoout_g(ispin), rhodn_tot_gspace)
1843 h_stress(:, :) = 0.0_dp
1844 CALL pw_poisson_solve(poisson_env, &
1845 density=rho_tot_gspace, &
1846 ehartree=ehartree, &
1847 vhartree=v_hartree_gspace, &
1848 h_stress=h_stress, &
1849 aux_density=rhodn_tot_gspace)
1851 virial%pv_ehartree = virial%pv_ehartree + h_stress/real(para_env%num_pe, dp)
1852 virial%pv_virial = virial%pv_virial + h_stress/real(para_env%num_pe, dp)
1854 IF (debug_stress)
THEN
1855 stdeb = fconv*(h_stress/real(para_env%num_pe, dp))
1856 CALL para_env%sum(stdeb)
1857 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
1858 'STRESS| GREEN 1st v_H[n_in]*n_out ', one_third_sum_diag(stdeb),
det_3x3(stdeb)
1862 virial%pv_calculate = .true.
1864 NULLIFY (v_rspace, v_tau_rspace)
1865 CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho, xc_section=ec_env%xc_section, &
1866 vxc_rho=v_rspace, vxc_tau=v_tau_rspace, exc=exc, just_energy=.false.)
1869 virial%pv_exc = virial%pv_exc - virial%pv_xc
1870 virial%pv_virial = virial%pv_virial - virial%pv_xc
1872 IF (debug_stress)
THEN
1873 stdeb = -1.0_dp*fconv*virial%pv_xc
1874 CALL para_env%sum(stdeb)
1875 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
1876 'STRESS| GGA 1st E_xc[Pin] ', one_third_sum_diag(stdeb),
det_3x3(stdeb)
1879 IF (
ASSOCIATED(v_rspace))
THEN
1880 DO ispin = 1, nspins
1881 CALL auxbas_pw_pool%give_back_pw(v_rspace(ispin))
1883 DEALLOCATE (v_rspace)
1885 IF (
ASSOCIATED(v_tau_rspace))
THEN
1886 DO ispin = 1, nspins
1887 CALL auxbas_pw_pool%give_back_pw(v_tau_rspace(ispin))
1889 DEALLOCATE (v_tau_rspace)
1891 CALL pw_zero(rhodn_tot_gspace)
1896 DO ispin = 1, nspins
1897 CALL pw_axpy(rho_r(ispin), rhoout_r(ispin), -1.0_dp)
1898 CALL pw_axpy(rho_g(ispin), rhoout_g(ispin), -1.0_dp)
1899 CALL pw_axpy(rhoout_g(ispin), rhodn_tot_gspace)
1900 IF (dft_control%use_kinetic_energy_density)
CALL pw_axpy(tau_r(ispin), tauout_r(ispin), -1.0_dp)
1904 IF (use_virial)
THEN
1907 h_stress(:, :) = 0.0_dp
1908 CALL pw_poisson_solve(poisson_env, &
1909 density=rhodn_tot_gspace, &
1910 ehartree=dehartree, &
1911 vhartree=v_hartree_gspace, &
1912 h_stress=h_stress, &
1913 aux_density=rho_tot_gspace)
1915 CALL auxbas_pw_pool%give_back_pw(rho_tot_gspace)
1917 virial%pv_ehartree = virial%pv_ehartree + h_stress/real(para_env%num_pe, dp)
1918 virial%pv_virial = virial%pv_virial + h_stress/real(para_env%num_pe, dp)
1920 IF (debug_stress)
THEN
1921 stdeb = fconv*(h_stress/real(para_env%num_pe, dp))
1922 CALL para_env%sum(stdeb)
1923 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
1924 'STRESS| GREEN 2nd V_H[dP]*n_in ', one_third_sum_diag(stdeb),
det_3x3(stdeb)
1929 CALL pw_poisson_solve(poisson_env, rhodn_tot_gspace, dehartree, &
1933 CALL pw_transfer(v_hartree_gspace, dv_hartree_rspace)
1934 CALL pw_scale(dv_hartree_rspace, dv_hartree_rspace%pw_grid%dvol)
1937 CALL pw_transfer(v_hartree_rspace, vtot_rspace)
1938 CALL pw_axpy(dv_hartree_rspace, vtot_rspace)
1939 IF (debug_forces) fodeb(1:3) = force(1)%rho_core(1:3, 1)
1940 IF (debug_stress .AND. use_virial) stdeb = virial%pv_ehartree
1941 CALL integrate_v_core_rspace(vtot_rspace, qs_env)
1942 IF (debug_forces)
THEN
1943 fodeb(1:3) = force(1)%rho_core(1:3, 1) - fodeb(1:3)
1944 CALL para_env%sum(fodeb)
1945 IF (iounit > 0)
WRITE (iounit,
"(T3,A,T33,3F16.8)")
"DEBUG:: Vtot*dncore", fodeb
1947 IF (debug_stress .AND. use_virial)
THEN
1948 stdeb = fconv*(virial%pv_ehartree - stdeb)
1949 CALL para_env%sum(stdeb)
1950 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
1951 'STRESS| Vtot*dncore', one_third_sum_diag(stdeb),
det_3x3(stdeb)
1957 xc_section => ec_env%xc_section
1959 IF (use_virial) virial%pv_xc = 0.0_dp
1960 NULLIFY (v_xc, v_xc_tau)
1961 CALL create_kernel(qs_env, &
1968 xc_section=xc_section, &
1969 compute_virial=use_virial, &
1970 virial_xc=virial%pv_xc)
1972 IF (use_virial)
THEN
1974 virial%pv_exc = virial%pv_exc + virial%pv_xc
1975 virial%pv_virial = virial%pv_virial + virial%pv_xc
1977 IF (debug_stress .AND. use_virial)
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 2nd f_Hxc[dP]*Pin ', one_third_sum_diag(stdeb),
det_3x3(stdeb)
1984 CALL get_qs_env(qs_env=qs_env, rho=rho, matrix_s_kp=matrix_s)
1985 NULLIFY (ec_env%matrix_hz)
1986 CALL dbcsr_allocate_matrix_set(ec_env%matrix_hz, nspins)
1987 DO ispin = 1, nspins
1988 ALLOCATE (ec_env%matrix_hz(ispin)%matrix)
1989 CALL dbcsr_create(ec_env%matrix_hz(ispin)%matrix, template=matrix_s(1, 1)%matrix)
1990 CALL dbcsr_copy(ec_env%matrix_hz(ispin)%matrix, matrix_s(1, 1)%matrix)
1991 CALL dbcsr_set(ec_env%matrix_hz(ispin)%matrix, 0.0_dp)
1993 CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
1995 IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
1996 IF (debug_stress .AND. use_virial) stdeb = virial%pv_virial
1999 IF (use_virial)
THEN
2000 pv_loc = virial%pv_virial
2003 DO ispin = 1, nspins
2004 CALL pw_scale(v_xc(ispin), v_xc(ispin)%pw_grid%dvol)
2005 CALL pw_axpy(dv_hartree_rspace, v_xc(ispin))
2006 CALL integrate_v_rspace(v_rspace=v_xc(ispin), &
2007 hmat=ec_env%matrix_hz(ispin), &
2008 pmat=matrix_p(ispin, 1), &
2010 calculate_forces=.true.)
2013 IF (debug_forces)
THEN
2014 fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
2015 CALL para_env%sum(fodeb)
2016 IF (iounit > 0)
WRITE (iounit,
"(T3,A,T33,3F16.8)")
"DEBUG:: Pin*dKdrho", fodeb
2018 IF (debug_stress .AND. use_virial)
THEN
2019 stdeb = fconv*(virial%pv_virial - stdeb)
2020 CALL para_env%sum(stdeb)
2021 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
2022 'STRESS| INT 2nd f_Hxc[dP]*Pin ', one_third_sum_diag(stdeb),
det_3x3(stdeb)
2025 IF (
ASSOCIATED(v_xc_tau))
THEN
2026 IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
2027 IF (debug_stress .AND. use_virial) stdeb = virial%pv_virial
2029 DO ispin = 1, nspins
2030 CALL pw_scale(v_xc_tau(ispin), v_xc_tau(ispin)%pw_grid%dvol)
2031 CALL integrate_v_rspace(v_rspace=v_xc_tau(ispin), &
2032 hmat=ec_env%matrix_hz(ispin), &
2033 pmat=matrix_p(ispin, 1), &
2035 compute_tau=.true., &
2036 calculate_forces=.true.)
2039 IF (debug_forces)
THEN
2040 fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
2041 CALL para_env%sum(fodeb)
2042 IF (iounit > 0)
WRITE (iounit,
"(T3,A,T33,3F16.8)")
"DEBUG:: Pin*dKtaudtau", fodeb
2044 IF (debug_stress .AND. use_virial)
THEN
2045 stdeb = fconv*(virial%pv_virial - stdeb)
2046 CALL para_env%sum(stdeb)
2047 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
2048 'STRESS| INT 2nd f_xctau[dP]*Pin ', one_third_sum_diag(stdeb),
det_3x3(stdeb)
2052 IF (use_virial)
THEN
2053 virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
2058 CALL dbcsr_allocate_matrix_set(scrm, nspins)
2059 DO ispin = 1, nspins
2060 ALLOCATE (scrm(ispin)%matrix)
2061 CALL dbcsr_create(scrm(ispin)%matrix, template=ec_env%matrix_ks(ispin, 1)%matrix)
2062 CALL dbcsr_copy(scrm(ispin)%matrix, ec_env%matrix_ks(ispin, 1)%matrix)
2063 CALL dbcsr_set(scrm(ispin)%matrix, 0.0_dp)
2067 NULLIFY (v_rspace, v_tau_rspace)
2069 CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho, xc_section=ec_env%xc_section, &
2070 vxc_rho=v_rspace, vxc_tau=v_tau_rspace, exc=eexc, just_energy=.false.)
2072 IF (use_virial)
THEN
2074 IF (
ASSOCIATED(v_rspace))
THEN
2075 DO ispin = 1, nspins
2077 eexc = eexc + pw_integral_ab(rhoout_r(ispin), v_rspace(ispin))
2080 IF (
ASSOCIATED(v_tau_rspace))
THEN
2081 DO ispin = 1, nspins
2083 eexc = eexc + pw_integral_ab(tauout_r(ispin), v_tau_rspace(ispin))
2088 IF (.NOT.
ASSOCIATED(v_rspace))
THEN
2089 ALLOCATE (v_rspace(nspins))
2090 DO ispin = 1, nspins
2091 CALL auxbas_pw_pool%create_pw(v_rspace(ispin))
2092 CALL pw_zero(v_rspace(ispin))
2098 IF (use_virial)
THEN
2099 pv_loc = virial%pv_virial
2102 IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
2103 IF (debug_stress .AND. use_virial) stdeb = virial%pv_virial
2104 DO ispin = 1, nspins
2106 CALL pw_scale(v_rspace(ispin), v_rspace(ispin)%pw_grid%dvol)
2107 CALL pw_axpy(v_hartree_rspace, v_rspace(ispin))
2109 CALL integrate_v_rspace(v_rspace=v_rspace(ispin), &
2111 pmat=ec_env%matrix_p(ispin, 1), &
2113 calculate_forces=.true., &
2114 basis_type=
"HARRIS", &
2115 task_list_external=ec_env%task_list)
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:: Pout*dVhxc ", 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 Pout*dVhxc ', one_third_sum_diag(stdeb),
det_3x3(stdeb)
2131 IF (use_virial)
THEN
2132 virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
2135 IF (
ASSOCIATED(v_tau_rspace))
THEN
2136 IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
2137 DO ispin = 1, nspins
2139 CALL pw_scale(v_tau_rspace(ispin), v_tau_rspace(ispin)%pw_grid%dvol)
2140 CALL integrate_v_rspace(v_rspace=v_tau_rspace(ispin), &
2142 pmat=ec_env%matrix_p(ispin, 1), &
2144 calculate_forces=.true., &
2145 compute_tau=.true., &
2146 basis_type=
"HARRIS", &
2147 task_list_external=ec_env%task_list)
2149 IF (debug_forces)
THEN
2150 fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
2151 CALL para_env%sum(fodeb)
2152 IF (iounit > 0)
WRITE (iounit,
"(T3,A,T33,3F16.8)")
"DEBUG:: Pout*dVhxc_tau ", fodeb
2161 ec_hfx_sections => section_vals_get_subs_vals(qs_env%input,
"DFT%ENERGY_CORRECTION%XC%HF")
2162 CALL section_vals_get(ec_hfx_sections, explicit=do_ec_hfx)
2166 IF (debug_forces) fodeb(1:3) = force(1)%fock_4c(1:3, 1)
2167 IF (use_virial) virial%pv_fock_4c = 0.0_dp
2169 CALL calculate_exx(qs_env=qs_env, &
2171 hfx_sections=ec_hfx_sections, &
2172 x_data=ec_env%x_data, &
2174 do_admm=ec_env%do_ec_admm, &
2175 calc_forces=.true., &
2176 reuse_hfx=ec_env%reuse_hfx, &
2177 do_im_time=.false., &
2178 e_ex_from_gw=dummy_real, &
2179 e_admm_from_gw=dummy_real2, &
2182 IF (use_virial)
THEN
2183 virial%pv_exx = virial%pv_exx - virial%pv_fock_4c
2184 virial%pv_virial = virial%pv_virial - virial%pv_fock_4c
2185 virial%pv_calculate = .false.
2187 IF (debug_forces)
THEN
2188 fodeb(1:3) = force(1)%fock_4c(1:3, 1) - fodeb(1:3)
2189 CALL para_env%sum(fodeb)
2190 IF (iounit > 0)
WRITE (iounit,
"(T3,A,T33,3F16.8)")
"DEBUG:: Pout*hfx ", fodeb
2192 IF (debug_stress .AND. use_virial)
THEN
2193 stdeb = -1.0_dp*fconv*virial%pv_fock_4c
2194 CALL para_env%sum(stdeb)
2195 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
2196 'STRESS| Pout*hfx ', one_third_sum_diag(stdeb),
det_3x3(stdeb)
2204 CALL dbcsr_deallocate_matrix_set(scrm)
2207 CALL auxbas_pw_pool%give_back_pw(v_hartree_rspace)
2208 DO ispin = 1, nspins
2209 CALL auxbas_pw_pool%give_back_pw(v_rspace(ispin))
2210 IF (
ASSOCIATED(v_tau_rspace))
THEN
2211 CALL auxbas_pw_pool%give_back_pw(v_tau_rspace(ispin))
2214 IF (
ASSOCIATED(v_tau_rspace))
DEALLOCATE (v_tau_rspace)
2217 IF (debug_forces) fodeb(1:3) = force(1)%core_overlap(1:3, 1)
2218 IF (debug_stress .AND. use_virial) stdeb = virial%pv_ecore_overlap
2219 CALL calculate_ecore_overlap(qs_env, para_env, .true., e_overlap_core=eovrl)
2220 IF (debug_forces)
THEN
2221 fodeb(1:3) = force(1)%core_overlap(1:3, 1) - fodeb(1:3)
2222 CALL para_env%sum(fodeb)
2223 IF (iounit > 0)
WRITE (iounit,
"(T3,A,T33,3F16.8)")
"DEBUG:: CoreOverlap", fodeb
2225 IF (debug_stress .AND. use_virial)
THEN
2226 stdeb = fconv*(stdeb - virial%pv_ecore_overlap)
2227 CALL para_env%sum(stdeb)
2228 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
2229 'STRESS| CoreOverlap ', one_third_sum_diag(stdeb),
det_3x3(stdeb)
2232 IF (debug_forces)
THEN
2233 CALL get_qs_env(qs_env, natom=natom, atomic_kind_set=atomic_kind_set)
2234 ALLOCATE (ftot(3, natom))
2235 CALL total_qs_force(ftot, force, atomic_kind_set)
2236 fodeb(1:3) = ftot(1:3, 1)
2238 CALL para_env%sum(fodeb)
2239 IF (iounit > 0)
WRITE (iounit,
"(T3,A,T33,3F16.8)")
"DEBUG:: Force Explicit", fodeb
2242 DEALLOCATE (v_rspace)
2244 CALL auxbas_pw_pool%give_back_pw(dv_hartree_rspace)
2245 CALL auxbas_pw_pool%give_back_pw(vtot_rspace)
2246 DO ispin = 1, nspins
2247 CALL auxbas_pw_pool%give_back_pw(rhoout_r(ispin))
2248 CALL auxbas_pw_pool%give_back_pw(rhoout_g(ispin))
2249 CALL auxbas_pw_pool%give_back_pw(v_xc(ispin))
2251 DEALLOCATE (rhoout_r, rhoout_g, v_xc)
2252 IF (
ASSOCIATED(tauout_r))
THEN
2253 DO ispin = 1, nspins
2254 CALL auxbas_pw_pool%give_back_pw(tauout_r(ispin))
2256 DEALLOCATE (tauout_r)
2258 IF (
ASSOCIATED(v_xc_tau))
THEN
2259 DO ispin = 1, nspins
2260 CALL auxbas_pw_pool%give_back_pw(v_xc_tau(ispin))
2262 DEALLOCATE (v_xc_tau)
2264 CALL auxbas_pw_pool%give_back_pw(v_hartree_gspace)
2265 CALL auxbas_pw_pool%give_back_pw(rhodn_tot_gspace)
2269 IF (use_virial)
THEN
2270 IF (qs_env%energy_correction)
THEN
2271 ec_env%ehartree = ehartree + dehartree
2272 ec_env%exc = exc + eexc
2276 IF (debug_stress .AND. use_virial)
THEN
2278 stdeb = -1.0_dp*fconv*ehartree
2279 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
2280 'STRESS| VOL 1st v_H[n_in]*n_out', one_third_sum_diag(stdeb),
det_3x3(stdeb)
2282 stdeb = -1.0_dp*fconv*exc
2283 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
2284 'STRESS| VOL 1st E_XC[n_in]', one_third_sum_diag(stdeb),
det_3x3(stdeb)
2286 stdeb = -1.0_dp*fconv*dehartree
2287 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
2288 'STRESS| VOL 2nd v_H[dP]*n_in', one_third_sum_diag(stdeb),
det_3x3(stdeb)
2290 stdeb = -1.0_dp*fconv*eexc
2291 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
2292 'STRESS| VOL 2nd v_XC[n_in]*dP', one_third_sum_diag(stdeb),
det_3x3(stdeb)
2297 TYPE(virial_type) :: virdeb
2300 CALL para_env%sum(virdeb%pv_overlap)
2301 CALL para_env%sum(virdeb%pv_ekinetic)
2302 CALL para_env%sum(virdeb%pv_ppl)
2303 CALL para_env%sum(virdeb%pv_ppnl)
2304 CALL para_env%sum(virdeb%pv_ecore_overlap)
2305 CALL para_env%sum(virdeb%pv_ehartree)
2306 CALL para_env%sum(virdeb%pv_exc)
2307 CALL para_env%sum(virdeb%pv_exx)
2308 CALL para_env%sum(virdeb%pv_vdw)
2309 CALL para_env%sum(virdeb%pv_mp2)
2310 CALL para_env%sum(virdeb%pv_nlcc)
2311 CALL para_env%sum(virdeb%pv_gapw)
2312 CALL para_env%sum(virdeb%pv_lrigpw)
2313 CALL para_env%sum(virdeb%pv_virial)
2314 CALL symmetrize_virial(virdeb)
2318 virdeb%pv_ehartree(i, i) = virdeb%pv_ehartree(i, i) - 2.0_dp*(ehartree + dehartree)
2319 virdeb%pv_virial(i, i) = virdeb%pv_virial(i, i) - exc - eexc &
2320 - 2.0_dp*(ehartree + dehartree)
2321 virdeb%pv_exc(i, i) = virdeb%pv_exc(i, i) - exc - eexc
2327 CALL para_env%sum(sttot)
2328 stdeb = fconv*(virdeb%pv_virial - sttot)
2329 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
2330 'STRESS| Explicit electronic stress ', one_third_sum_diag(stdeb),
det_3x3(stdeb)
2332 stdeb = fconv*(virdeb%pv_virial)
2333 IF (iounit > 0)
WRITE (unit=iounit, fmt=
"(T2,A,T41,2(1X,ES19.11))") &
2334 'STRESS| Explicit total stress ', one_third_sum_diag(stdeb),
det_3x3(stdeb)
2336 CALL write_stress_tensor_components(virdeb, iounit, cell)
2337 CALL write_stress_tensor(virdeb%pv_virial, iounit, cell, .false.)
2342 CALL timestop(handle)
2344 END SUBROUTINE ec_build_ks_matrix_force
2354 SUBROUTINE ec_ks_solver(qs_env, ec_env)
2356 TYPE(qs_environment_type),
POINTER :: qs_env
2357 TYPE(energy_correction_type),
POINTER :: ec_env
2359 CHARACTER(LEN=*),
PARAMETER :: routinen =
'ec_ks_solver'
2361 CHARACTER(LEN=default_string_length) :: headline
2362 INTEGER :: handle, ispin, nspins
2363 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: ksmat, pmat, smat, wmat
2364 TYPE(dft_control_type),
POINTER :: dft_control
2366 CALL timeset(routinen, handle)
2368 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
2369 nspins = dft_control%nspins
2372 IF (.NOT.
ASSOCIATED(ec_env%matrix_p))
THEN
2373 headline =
"DENSITY MATRIX"
2374 CALL dbcsr_allocate_matrix_set(ec_env%matrix_p, nspins, 1)
2375 DO ispin = 1, nspins
2376 ALLOCATE (ec_env%matrix_p(ispin, 1)%matrix)
2377 CALL dbcsr_create(ec_env%matrix_p(ispin, 1)%matrix, name=trim(headline), &
2378 template=ec_env%matrix_s(1, 1)%matrix)
2379 CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_p(ispin, 1)%matrix, ec_env%sab_orb)
2383 IF (.NOT.
ASSOCIATED(ec_env%matrix_w))
THEN
2384 headline =
"ENERGY WEIGHTED DENSITY MATRIX"
2385 CALL dbcsr_allocate_matrix_set(ec_env%matrix_w, nspins, 1)
2386 DO ispin = 1, nspins
2387 ALLOCATE (ec_env%matrix_w(ispin, 1)%matrix)
2388 CALL dbcsr_create(ec_env%matrix_w(ispin, 1)%matrix, name=trim(headline), &
2389 template=ec_env%matrix_s(1, 1)%matrix)
2390 CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_w(ispin, 1)%matrix, ec_env%sab_orb)
2394 IF (ec_env%mao)
THEN
2395 CALL mao_create_matrices(ec_env, ksmat, smat, pmat, wmat)
2397 ksmat => ec_env%matrix_ks
2398 smat => ec_env%matrix_s
2399 pmat => ec_env%matrix_p
2400 wmat => ec_env%matrix_w
2403 SELECT CASE (ec_env%ks_solver)
2404 CASE (ec_diagonalization)
2405 CALL ec_diag_solver(qs_env, ksmat, smat, pmat, wmat)
2407 CALL ec_ot_diag_solver(qs_env, ec_env, ksmat, smat, pmat, wmat)
2408 CASE (ec_matrix_sign, ec_matrix_trs4, ec_matrix_tc2)
2409 CALL ec_ls_init(qs_env, ksmat, smat)
2410 CALL ec_ls_solver(qs_env, pmat, wmat, ec_ls_method=ec_env%ks_solver)
2466 IF (ec_env%mao)
THEN
2467 CALL mao_release_matrices(ec_env, ksmat, smat, pmat, wmat)
2470 CALL timestop(handle)
2472 END SUBROUTINE ec_ks_solver
2485 SUBROUTINE mao_create_matrices(ec_env, ksmat, smat, pmat, wmat)
2487 TYPE(energy_correction_type),
POINTER :: ec_env
2488 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: ksmat, smat, pmat, wmat
2490 CHARACTER(LEN=*),
PARAMETER :: routinen =
'mao_create_matrices'
2492 INTEGER :: handle, ispin, nspins
2493 INTEGER,
DIMENSION(:),
POINTER :: col_blk_sizes
2494 TYPE(dbcsr_distribution_type) :: dbcsr_dist
2495 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: mao_coef
2496 TYPE(dbcsr_type) :: cgmat
2498 CALL timeset(routinen, handle)
2500 mao_coef => ec_env%mao_coef
2502 NULLIFY (ksmat, smat, pmat, wmat)
2503 nspins =
SIZE(ec_env%matrix_ks, 1)
2504 CALL dbcsr_get_info(mao_coef(1)%matrix, col_blk_size=col_blk_sizes, distribution=dbcsr_dist)
2505 CALL dbcsr_allocate_matrix_set(ksmat, nspins, 1)
2506 CALL dbcsr_allocate_matrix_set(smat, nspins, 1)
2507 DO ispin = 1, nspins
2508 ALLOCATE (ksmat(ispin, 1)%matrix)
2509 CALL dbcsr_create(ksmat(ispin, 1)%matrix, dist=dbcsr_dist, name=
"MAO KS mat", &
2510 matrix_type=dbcsr_type_symmetric, row_blk_size=col_blk_sizes, &
2511 col_blk_size=col_blk_sizes, nze=0)
2512 ALLOCATE (smat(ispin, 1)%matrix)
2513 CALL dbcsr_create(smat(ispin, 1)%matrix, dist=dbcsr_dist, name=
"MAO S mat", &
2514 matrix_type=dbcsr_type_symmetric, row_blk_size=col_blk_sizes, &
2515 col_blk_size=col_blk_sizes, nze=0)
2518 CALL dbcsr_create(cgmat, name=
"TEMP matrix", template=mao_coef(1)%matrix)
2519 DO ispin = 1, nspins
2520 CALL dbcsr_multiply(
"N",
"N", 1.0_dp, ec_env%matrix_s(1, 1)%matrix, mao_coef(ispin)%matrix, &
2522 CALL dbcsr_multiply(
"T",
"N", 1.0_dp, mao_coef(ispin)%matrix, cgmat, 0.0_dp, smat(ispin, 1)%matrix)
2523 CALL dbcsr_multiply(
"N",
"N", 1.0_dp, ec_env%matrix_ks(1, 1)%matrix, mao_coef(ispin)%matrix, &
2525 CALL dbcsr_multiply(
"T",
"N", 1.0_dp, mao_coef(ispin)%matrix, cgmat, 0.0_dp, ksmat(ispin, 1)%matrix)
2527 CALL dbcsr_release(cgmat)
2529 CALL dbcsr_allocate_matrix_set(pmat, nspins, 1)
2530 DO ispin = 1, nspins
2531 ALLOCATE (pmat(ispin, 1)%matrix)
2532 CALL dbcsr_create(pmat(ispin, 1)%matrix, template=smat(1, 1)%matrix, name=
"MAO P mat")
2533 CALL cp_dbcsr_alloc_block_from_nbl(pmat(ispin, 1)%matrix, ec_env%sab_orb)
2536 CALL dbcsr_allocate_matrix_set(wmat, nspins, 1)
2537 DO ispin = 1, nspins
2538 ALLOCATE (wmat(ispin, 1)%matrix)
2539 CALL dbcsr_create(wmat(ispin, 1)%matrix, template=smat(1, 1)%matrix, name=
"MAO W mat")
2540 CALL cp_dbcsr_alloc_block_from_nbl(wmat(ispin, 1)%matrix, ec_env%sab_orb)
2543 CALL timestop(handle)
2545 END SUBROUTINE mao_create_matrices
2558 SUBROUTINE mao_release_matrices(ec_env, ksmat, smat, pmat, wmat)
2560 TYPE(energy_correction_type),
POINTER :: ec_env
2561 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: ksmat, smat, pmat, wmat
2563 CHARACTER(LEN=*),
PARAMETER :: routinen =
'mao_release_matrices'
2565 INTEGER :: handle, ispin, nspins
2566 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: mao_coef
2567 TYPE(dbcsr_type) :: cgmat
2569 CALL timeset(routinen, handle)
2571 mao_coef => ec_env%mao_coef
2572 nspins =
SIZE(mao_coef, 1)
2575 CALL dbcsr_create(cgmat, name=
"TEMP matrix", template=mao_coef(1)%matrix)
2576 DO ispin = 1, nspins
2577 CALL dbcsr_multiply(
"N",
"N", 1.0_dp, mao_coef(ispin)%matrix, pmat(ispin, 1)%matrix, 0.0_dp, cgmat)
2578 CALL dbcsr_multiply(
"N",
"T", 1.0_dp, mao_coef(ispin)%matrix, cgmat, 0.0_dp, &
2579 ec_env%matrix_p(ispin, 1)%matrix, retain_sparsity=.true.)
2580 CALL dbcsr_multiply(
"N",
"N", 1.0_dp, mao_coef(ispin)%matrix, wmat(ispin, 1)%matrix, 0.0_dp, cgmat)
2581 CALL dbcsr_multiply(
"N",
"T", 1.0_dp, mao_coef(ispin)%matrix, cgmat, 0.0_dp, &
2582 ec_env%matrix_w(ispin, 1)%matrix, retain_sparsity=.true.)
2584 CALL dbcsr_release(cgmat)
2586 CALL dbcsr_deallocate_matrix_set(ksmat)
2587 CALL dbcsr_deallocate_matrix_set(smat)
2588 CALL dbcsr_deallocate_matrix_set(pmat)
2589 CALL dbcsr_deallocate_matrix_set(wmat)
2591 CALL timestop(handle)
2593 END SUBROUTINE mao_release_matrices
2606 SUBROUTINE ec_diag_solver(qs_env, matrix_ks, matrix_s, matrix_p, matrix_w)
2608 TYPE(qs_environment_type),
POINTER :: qs_env
2609 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_ks, matrix_s, matrix_p, matrix_w
2611 CHARACTER(LEN=*),
PARAMETER :: routinen =
'ec_diag_solver'
2613 INTEGER :: handle, ispin, nmo(2), nsize, nspins
2614 REAL(kind=dp) :: eps_filter, focc(2)
2615 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:) :: eigenvalues
2616 TYPE(cp_blacs_env_type),
POINTER :: blacs_env
2617 TYPE(cp_fm_struct_type),
POINTER :: fm_struct
2618 TYPE(cp_fm_type) :: fm_ks, fm_mo, fm_ortho
2619 TYPE(dbcsr_type),
POINTER :: buf1_dbcsr, buf2_dbcsr, buf3_dbcsr, &
2620 ortho_dbcsr, ref_matrix
2621 TYPE(dft_control_type),
POINTER :: dft_control
2622 TYPE(mp_para_env_type),
POINTER :: para_env
2624 CALL timeset(routinen, handle)
2626 NULLIFY (blacs_env, para_env)
2627 CALL get_qs_env(qs_env=qs_env, blacs_env=blacs_env, para_env=para_env)
2629 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
2630 eps_filter = dft_control%qs_control%eps_filter_matrix
2631 nspins = dft_control%nspins
2634 CALL get_qs_env(qs_env=qs_env, nelectron_spin=nmo)
2636 IF (nspins == 1)
THEN
2641 CALL dbcsr_get_info(matrix_ks(1, 1)%matrix, nfullrows_total=nsize)
2642 ALLOCATE (eigenvalues(nsize))
2644 NULLIFY (fm_struct, ref_matrix)
2645 CALL cp_fm_struct_create(fm_struct, context=blacs_env, nrow_global=nsize, &
2646 ncol_global=nsize, para_env=para_env)
2647 CALL cp_fm_create(fm_ortho, fm_struct)
2648 CALL cp_fm_create(fm_ks, fm_struct)
2649 CALL cp_fm_create(fm_mo, fm_struct)
2650 CALL cp_fm_struct_release(fm_struct)
2653 ref_matrix => matrix_s(1, 1)%matrix
2654 NULLIFY (ortho_dbcsr, buf1_dbcsr, buf2_dbcsr, buf3_dbcsr)
2655 CALL dbcsr_init_p(ortho_dbcsr)
2656 CALL dbcsr_create(ortho_dbcsr, template=ref_matrix, &
2657 matrix_type=dbcsr_type_no_symmetry)
2658 CALL dbcsr_init_p(buf1_dbcsr)
2659 CALL dbcsr_create(buf1_dbcsr, template=ref_matrix, &
2660 matrix_type=dbcsr_type_no_symmetry)
2661 CALL dbcsr_init_p(buf2_dbcsr)
2662 CALL dbcsr_create(buf2_dbcsr, template=ref_matrix, &
2663 matrix_type=dbcsr_type_no_symmetry)
2664 CALL dbcsr_init_p(buf3_dbcsr)
2665 CALL dbcsr_create(buf3_dbcsr, template=ref_matrix, &
2666 matrix_type=dbcsr_type_no_symmetry)
2668 ref_matrix => matrix_s(1, 1)%matrix
2669 CALL copy_dbcsr_to_fm(ref_matrix, fm_ortho)
2670 CALL cp_fm_cholesky_decompose(fm_ortho)
2671 CALL cp_fm_triangular_invert(fm_ortho)
2672 CALL cp_fm_set_all(fm_ks, 0.0_dp)
2673 CALL cp_fm_to_fm_triangular(fm_ortho, fm_ks,
"U")
2674 CALL copy_fm_to_dbcsr(fm_ks, ortho_dbcsr)
2675 DO ispin = 1, nspins
2677 CALL dbcsr_desymmetrize(matrix_ks(ispin, 1)%matrix, buf1_dbcsr)
2678 CALL dbcsr_multiply(
"N",
"N", 1.0_dp, buf1_dbcsr, ortho_dbcsr, &
2679 0.0_dp, buf2_dbcsr, filter_eps=eps_filter)
2680 CALL dbcsr_multiply(
"T",
"N", 1.0_dp, ortho_dbcsr, buf2_dbcsr, &
2681 0.0_dp, buf1_dbcsr, filter_eps=eps_filter)
2683 CALL copy_dbcsr_to_fm(buf1_dbcsr, fm_ks)
2684 CALL choose_eigv_solver(fm_ks, fm_mo, eigenvalues)
2686 CALL copy_fm_to_dbcsr(fm_mo, buf1_dbcsr)
2687 CALL dbcsr_multiply(
"N",
"N", 1.0_dp, ortho_dbcsr, buf1_dbcsr, &
2688 0.0_dp, buf2_dbcsr, filter_eps=eps_filter)
2690 CALL dbcsr_set(matrix_p(ispin, 1)%matrix, 0.0_dp)
2691 CALL dbcsr_multiply(
"N",
"T", focc(ispin), buf2_dbcsr, buf2_dbcsr, &
2692 1.0_dp, matrix_p(ispin, 1)%matrix, &
2693 retain_sparsity=.true., last_k=nmo(ispin))
2696 CALL dbcsr_set(matrix_w(ispin, 1)%matrix, 0.0_dp)
2697 CALL cp_fm_column_scale(fm_mo, eigenvalues)
2698 CALL copy_fm_to_dbcsr(fm_mo, buf1_dbcsr)
2699 CALL dbcsr_multiply(
"N",
"N", 1.0_dp, ortho_dbcsr, buf1_dbcsr, &
2700 0.0_dp, buf3_dbcsr, filter_eps=eps_filter)
2701 CALL dbcsr_multiply(
"N",
"T", focc(ispin), buf2_dbcsr, buf3_dbcsr, &
2702 1.0_dp, matrix_w(ispin, 1)%matrix, &
2703 retain_sparsity=.true., last_k=nmo(ispin))
2706 CALL cp_fm_release(fm_ks)
2707 CALL cp_fm_release(fm_mo)
2708 CALL cp_fm_release(fm_ortho)
2709 CALL dbcsr_release(ortho_dbcsr)
2710 CALL dbcsr_release(buf1_dbcsr)
2711 CALL dbcsr_release(buf2_dbcsr)
2712 CALL dbcsr_release(buf3_dbcsr)
2713 DEALLOCATE (ortho_dbcsr, buf1_dbcsr, buf2_dbcsr, buf3_dbcsr)
2714 DEALLOCATE (eigenvalues)
2716 CALL timestop(handle)
2718 END SUBROUTINE ec_diag_solver
2726 SUBROUTINE ec_energy(ec_env, unit_nr)
2727 TYPE(energy_correction_type) :: ec_env
2728 INTEGER,
INTENT(IN) :: unit_nr
2730 CHARACTER(LEN=*),
PARAMETER :: routinen =
'ec_energy'
2732 INTEGER :: handle, ispin, nspins
2733 REAL(kind=dp) :: eband, trace
2735 CALL timeset(routinen, handle)
2737 nspins =
SIZE(ec_env%matrix_ks, 1)
2738 DO ispin = 1, nspins
2739 CALL dbcsr_dot(ec_env%matrix_p(ispin, 1)%matrix, ec_env%matrix_s(1, 1)%matrix, trace)
2740 IF (unit_nr > 0)
WRITE (unit_nr,
'(T3,A,T65,F16.10)')
'Tr[PS] ', trace
2744 SELECT CASE (ec_env%energy_functional)
2745 CASE (ec_functional_harris)
2749 DO ispin = 1, nspins
2750 CALL dbcsr_dot(ec_env%matrix_ks(ispin, 1)%matrix, ec_env%matrix_p(ispin, 1)%matrix, trace)
2751 eband = eband + trace
2753 ec_env%eband = eband + ec_env%efield_nuclear
2756 ec_env%etotal = ec_env%eband + ec_env%ehartree + ec_env%exc - ec_env%vhxc + ec_env%edispersion - ec_env%ex
2757 IF (unit_nr > 0)
THEN
2758 WRITE (unit_nr,
'(T3,A,T56,F25.15)')
"Eband ", ec_env%eband
2759 WRITE (unit_nr,
'(T3,A,T56,F25.15)')
"Ehartree ", ec_env%ehartree
2760 WRITE (unit_nr,
'(T3,A,T56,F25.15)')
"Exc ", ec_env%exc
2761 WRITE (unit_nr,
'(T3,A,T56,F25.15)')
"Ex ", ec_env%ex
2762 WRITE (unit_nr,
'(T3,A,T56,F25.15)')
"Evhxc ", ec_env%vhxc
2763 WRITE (unit_nr,
'(T3,A,T56,F25.15)')
"Edisp ", ec_env%edispersion
2764 WRITE (unit_nr,
'(T3,A,T56,F25.15)')
"Etotal Harris Functional ", ec_env%etotal
2767 CASE (ec_functional_dc)
2770 CALL calculate_ptrace(ec_env%matrix_h, ec_env%matrix_p, ec_env%ecore,
SIZE(ec_env%matrix_p, 1))
2772 ec_env%ecore = ec_env%ecore + ec_env%efield_nuclear
2773 ec_env%etotal = ec_env%ecore + ec_env%ehartree + ec_env%exc + ec_env%edispersion &
2774 + ec_env%ex + ec_env%exc_aux_fit
2776 IF (unit_nr > 0)
THEN
2777 WRITE (unit_nr,
'(T3,A,T56,F25.15)')
"Ecore ", ec_env%ecore
2778 WRITE (unit_nr,
'(T3,A,T56,F25.15)')
"Ehartree ", ec_env%ehartree
2779 WRITE (unit_nr,
'(T3,A,T56,F25.15)')
"Exc ", ec_env%exc
2780 WRITE (unit_nr,
'(T3,A,T56,F25.15)')
"Ex ", ec_env%ex
2781 WRITE (unit_nr,
'(T3,A,T56,F25.15)')
"Exc_aux_fit", ec_env%exc_aux_fit
2782 WRITE (unit_nr,
'(T3,A,T56,F25.15)')
"Edisp ", ec_env%edispersion
2783 WRITE (unit_nr,
'(T3,A,T56,F25.15)')
"Etotal Energy Functional ", ec_env%etotal
2792 CALL timestop(handle)
2794 END SUBROUTINE ec_energy
2806 SUBROUTINE ec_build_neighborlist(qs_env, ec_env)
2807 TYPE(qs_environment_type),
POINTER :: qs_env
2808 TYPE(energy_correction_type),
POINTER :: ec_env
2810 CHARACTER(LEN=*),
PARAMETER :: routinen =
'ec_build_neighborlist'
2812 INTEGER :: handle, ikind, nkind, zat
2813 LOGICAL :: gth_potential_present, &
2814 sgp_potential_present, &
2815 skip_load_balance_distributed
2816 LOGICAL,
ALLOCATABLE,
DIMENSION(:) :: default_present, orb_present, &
2817 ppl_present, ppnl_present
2818 REAL(dp) :: subcells
2819 REAL(dp),
ALLOCATABLE,
DIMENSION(:) :: c_radius, orb_radius, ppl_radius, &
2821 REAL(dp),
ALLOCATABLE,
DIMENSION(:, :) :: pair_radius
2822 TYPE(atomic_kind_type),
DIMENSION(:),
POINTER :: atomic_kind_set
2823 TYPE(cell_type),
POINTER :: cell
2824 TYPE(dft_control_type),
POINTER :: dft_control
2825 TYPE(distribution_1d_type),
POINTER :: distribution_1d
2826 TYPE(distribution_2d_type),
POINTER :: distribution_2d
2827 TYPE(gth_potential_type),
POINTER :: gth_potential
2828 TYPE(gto_basis_set_type),
POINTER :: basis_set
2829 TYPE(local_atoms_type),
ALLOCATABLE,
DIMENSION(:) :: atom2d
2830 TYPE(molecule_type),
DIMENSION(:),
POINTER :: molecule_set
2831 TYPE(neighbor_list_set_p_type),
DIMENSION(:), &
2832 POINTER :: sab_cn, sab_vdw
2833 TYPE(particle_type),
DIMENSION(:),
POINTER :: particle_set
2834 TYPE(qs_dispersion_type),
POINTER :: dispersion_env
2835 TYPE(qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
2836 TYPE(qs_kind_type),
POINTER :: qs_kind
2837 TYPE(qs_ks_env_type),
POINTER :: ks_env
2838 TYPE(sgp_potential_type),
POINTER :: sgp_potential
2840 CALL timeset(routinen, handle)
2842 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set)
2843 CALL get_qs_kind_set(qs_kind_set, gth_potential_present=gth_potential_present, &
2844 sgp_potential_present=sgp_potential_present)
2845 nkind =
SIZE(qs_kind_set)
2846 ALLOCATE (c_radius(nkind), default_present(nkind))
2847 ALLOCATE (orb_radius(nkind), ppl_radius(nkind), ppnl_radius(nkind))
2848 ALLOCATE (orb_present(nkind), ppl_present(nkind), ppnl_present(nkind))
2849 ALLOCATE (pair_radius(nkind, nkind))
2850 ALLOCATE (atom2d(nkind))
2852 CALL get_qs_env(qs_env, &
2853 atomic_kind_set=atomic_kind_set, &
2855 distribution_2d=distribution_2d, &
2856 local_particles=distribution_1d, &
2857 particle_set=particle_set, &
2858 molecule_set=molecule_set)
2860 CALL atom2d_build(atom2d, distribution_1d, distribution_2d, atomic_kind_set, &
2861 molecule_set, .false., particle_set)
2864 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom2d(ikind)%list)
2865 qs_kind => qs_kind_set(ikind)
2866 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type=
"HARRIS")
2867 IF (
ASSOCIATED(basis_set))
THEN
2868 orb_present(ikind) = .true.
2869 CALL get_gto_basis_set(gto_basis_set=basis_set, kind_radius=orb_radius(ikind))
2871 orb_present(ikind) = .false.
2872 orb_radius(ikind) = 0.0_dp
2874 CALL get_qs_kind(qs_kind, gth_potential=gth_potential, sgp_potential=sgp_potential)
2875 IF (gth_potential_present .OR. sgp_potential_present)
THEN
2876 IF (
ASSOCIATED(gth_potential))
THEN
2877 CALL get_potential(potential=gth_potential, &
2878 ppl_present=ppl_present(ikind), &
2879 ppl_radius=ppl_radius(ikind), &
2880 ppnl_present=ppnl_present(ikind), &
2881 ppnl_radius=ppnl_radius(ikind))
2882 ELSE IF (
ASSOCIATED(sgp_potential))
THEN
2883 CALL get_potential(potential=sgp_potential, &
2884 ppl_present=ppl_present(ikind), &
2885 ppl_radius=ppl_radius(ikind), &
2886 ppnl_present=ppnl_present(ikind), &
2887 ppnl_radius=ppnl_radius(ikind))
2889 ppl_present(ikind) = .false.
2890 ppl_radius(ikind) = 0.0_dp
2891 ppnl_present(ikind) = .false.
2892 ppnl_radius(ikind) = 0.0_dp
2897 CALL section_vals_val_get(qs_env%input,
"DFT%SUBCELLS", r_val=subcells)
2900 CALL pair_radius_setup(orb_present, orb_present, orb_radius, orb_radius, pair_radius)
2901 CALL build_neighbor_lists(ec_env%sab_orb, particle_set, atom2d, cell, pair_radius, &
2902 subcells=subcells, nlname=
"sab_orb")
2904 IF (gth_potential_present .OR. sgp_potential_present)
THEN
2905 IF (any(ppl_present))
THEN
2906 CALL pair_radius_setup(orb_present, ppl_present, orb_radius, ppl_radius, pair_radius)
2907 CALL build_neighbor_lists(ec_env%sac_ppl, particle_set, atom2d, cell, pair_radius, &
2908 subcells=subcells, operator_type=
"ABC", nlname=
"sac_ppl")
2911 IF (any(ppnl_present))
THEN
2912 CALL pair_radius_setup(orb_present, ppnl_present, orb_radius, ppnl_radius, pair_radius)
2913 CALL build_neighbor_lists(ec_env%sap_ppnl, particle_set, atom2d, cell, pair_radius, &
2914 subcells=subcells, operator_type=
"ABBA", nlname=
"sap_ppnl")
2919 c_radius(:) = 0.0_dp
2920 dispersion_env => ec_env%dispersion_env
2921 sab_vdw => dispersion_env%sab_vdw
2922 sab_cn => dispersion_env%sab_cn
2923 IF (dispersion_env%type == xc_vdw_fun_pairpot)
THEN
2924 c_radius(:) = dispersion_env%rc_disp
2925 default_present = .true.
2926 CALL pair_radius_setup(default_present, default_present, c_radius, c_radius, pair_radius)
2927 CALL build_neighbor_lists(sab_vdw, particle_set, atom2d, cell, pair_radius, &
2928 subcells=subcells, operator_type=
"PP", nlname=
"sab_vdw")
2929 dispersion_env%sab_vdw => sab_vdw
2930 IF (dispersion_env%pp_type == vdw_pairpot_dftd3 .OR. &
2931 dispersion_env%pp_type == vdw_pairpot_dftd3bj)
THEN
2934 CALL get_atomic_kind(atomic_kind_set(ikind), z=zat)
2935 c_radius(ikind) = 4._dp*ptable(zat)%covalent_radius*bohr
2937 CALL pair_radius_setup(default_present, default_present, c_radius, c_radius, pair_radius)
2938 CALL build_neighbor_lists(sab_cn, particle_set, atom2d, cell, pair_radius, &
2939 subcells=subcells, operator_type=
"PP", nlname=
"sab_cn")
2940 dispersion_env%sab_cn => sab_cn
2945 CALL atom2d_cleanup(atom2d)
2947 DEALLOCATE (orb_present, default_present, ppl_present, ppnl_present)
2948 DEALLOCATE (orb_radius, ppl_radius, ppnl_radius, c_radius)
2949 DEALLOCATE (pair_radius)
2952 CALL get_qs_env(qs_env, ks_env=ks_env, dft_control=dft_control)
2953 skip_load_balance_distributed = dft_control%qs_control%skip_load_balance_distributed
2954 IF (
ASSOCIATED(ec_env%task_list))
CALL deallocate_task_list(ec_env%task_list)
2955 CALL allocate_task_list(ec_env%task_list)
2956 CALL generate_qs_task_list(ks_env, ec_env%task_list, &
2957 reorder_rs_grid_ranks=.false., soft_valid=.false., &
2958 skip_load_balance_distributed=skip_load_balance_distributed, &
2959 basis_type=
"HARRIS", sab_orb_external=ec_env%sab_orb)
2961 CALL timestop(handle)
2963 END SUBROUTINE ec_build_neighborlist
2970 SUBROUTINE ec_properties(qs_env, ec_env)
2971 TYPE(qs_environment_type),
POINTER :: qs_env
2972 TYPE(energy_correction_type),
POINTER :: ec_env
2974 CHARACTER(LEN=*),
PARAMETER :: routinen =
'ec_properties'
2976 CHARACTER(LEN=8),
DIMENSION(3) :: rlab
2977 CHARACTER(LEN=default_path_length) :: filename, my_pos_voro
2978 CHARACTER(LEN=default_string_length) :: description
2979 INTEGER :: akind, handle, i, ia, iatom, idir, ikind, iounit, ispin, maxmom, nspins, &
2980 reference, should_print_bqb, should_print_voro, unit_nr, unit_nr_voro
2981 LOGICAL :: append_voro, magnetic, periodic, &
2983 REAL(kind=dp) :: charge, dd, focc, tmp
2984 REAL(kind=dp),
DIMENSION(3) :: cdip, pdip, rcc, rdip, ria, tdip
2985 REAL(kind=dp),
DIMENSION(:),
POINTER :: ref_point
2986 TYPE(atomic_kind_type),
POINTER :: atomic_kind
2987 TYPE(cell_type),
POINTER :: cell
2988 TYPE(cp_logger_type),
POINTER :: logger
2989 TYPE(cp_result_type),
POINTER :: results
2990 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_s, moments
2991 TYPE(dft_control_type),
POINTER :: dft_control
2992 TYPE(distribution_1d_type),
POINTER :: local_particles
2993 TYPE(mp_para_env_type),
POINTER :: para_env
2994 TYPE(particle_type),
DIMENSION(:),
POINTER :: particle_set
2995 TYPE(pw_env_type),
POINTER :: pw_env
2996 TYPE(pw_pool_p_type),
DIMENSION(:),
POINTER :: pw_pools
2997 TYPE(pw_pool_type),
POINTER :: auxbas_pw_pool
2998 TYPE(pw_r3d_rs_type) :: rho_elec_rspace
2999 TYPE(qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
3000 TYPE(section_vals_type),
POINTER :: ec_section, print_key, print_key_bqb, &
3003 CALL timeset(routinen, handle)
3009 logger => cp_get_default_logger()
3010 IF (logger%para_env%is_source())
THEN
3011 iounit = cp_logger_get_default_unit_nr(logger, local=.true.)
3016 NULLIFY (dft_control)
3017 CALL get_qs_env(qs_env, dft_control=dft_control)
3018 nspins = dft_control%nspins
3020 ec_section => section_vals_get_subs_vals(qs_env%input,
"DFT%ENERGY_CORRECTION")
3021 print_key => section_vals_get_subs_vals(section_vals=ec_section, &
3022 subsection_name=
"PRINT%MOMENTS")
3024 IF (btest(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file))
THEN
3026 maxmom = section_get_ival(section_vals=ec_section, &
3027 keyword_name=
"PRINT%MOMENTS%MAX_MOMENT")
3028 periodic = section_get_lval(section_vals=ec_section, &
3029 keyword_name=
"PRINT%MOMENTS%PERIODIC")
3030 reference = section_get_ival(section_vals=ec_section, &
3031 keyword_name=
"PRINT%MOMENTS%REFERENCE")
3032 magnetic = section_get_lval(section_vals=ec_section, &
3033 keyword_name=
"PRINT%MOMENTS%MAGNETIC")
3035 CALL section_vals_val_get(ec_section,
"PRINT%MOMENTS%REF_POINT", r_vals=ref_point)
3036 unit_nr = cp_print_key_unit_nr(logger=logger, basis_section=ec_section, &
3037 print_key_path=
"PRINT%MOMENTS", extension=
".dat", &
3038 middle_name=
"moments", log_filename=.false.)
3040 IF (iounit > 0)
THEN
3041 IF (unit_nr /= iounit .AND. unit_nr > 0)
THEN
3042 INQUIRE (unit=unit_nr, name=filename)
3043 WRITE (unit=iounit, fmt=
"(/,T2,A,2(/,T3,A),/)") &
3044 "MOMENTS",
"The electric/magnetic moments are written to file:", &
3047 WRITE (unit=iounit, fmt=
"(/,T2,A)")
"ELECTRIC/MAGNETIC MOMENTS"
3052 cpabort(
"Periodic moments not implemented with EC")
3054 cpassert(maxmom < 2)
3055 cpassert(.NOT. magnetic)
3056 IF (maxmom == 1)
THEN
3057 CALL get_qs_env(qs_env=qs_env, cell=cell, para_env=para_env)
3059 CALL get_reference_point(rcc, qs_env=qs_env, reference=reference, ref_point=ref_point)
3062 CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, &
3063 qs_kind_set=qs_kind_set, local_particles=local_particles)
3064 DO ikind = 1,
SIZE(local_particles%n_el)
3065 DO ia = 1, local_particles%n_el(ikind)
3066 iatom = local_particles%list(ikind)%array(ia)
3068 ria =
pbc(particle_set(iatom)%r - rcc, cell) + rcc
3070 atomic_kind => particle_set(iatom)%atomic_kind
3071 CALL get_atomic_kind(atomic_kind, kind_number=akind)
3072 CALL get_qs_kind(qs_kind_set(akind), core_charge=charge)
3073 cdip(1:3) = cdip(1:3) - charge*ria(1:3)
3076 CALL para_env%sum(cdip)
3079 CALL ec_efield_integrals(qs_env, ec_env, rcc)
3082 DO ispin = 1, nspins
3084 CALL dbcsr_dot(ec_env%matrix_p(ispin, 1)%matrix, &
3085 ec_env%efield%dipmat(idir)%matrix, tmp)
3086 pdip(idir) = pdip(idir) + tmp
3091 CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s)
3093 CALL dbcsr_allocate_matrix_set(moments, 4)
3095 ALLOCATE (moments(i)%matrix)
3096 CALL dbcsr_copy(moments(i)%matrix, matrix_s(1)%matrix,
"Moments")
3097 CALL dbcsr_set(moments(i)%matrix, 0.0_dp)
3099 CALL build_local_moment_matrix(qs_env, moments, 1, ref_point=rcc)
3102 IF (nspins == 2) focc = 1.0_dp
3104 DO ispin = 1, nspins
3106 CALL dbcsr_dot(ec_env%matrix_z(ispin)%matrix, moments(idir)%matrix, tmp)
3107 rdip(idir) = rdip(idir) + tmp
3110 CALL dbcsr_deallocate_matrix_set(moments)
3112 tdip = -(rdip + pdip + cdip)
3113 IF (unit_nr > 0)
THEN
3114 WRITE (unit_nr,
"(T3,A)")
"Dipoles are based on the traditional operator."
3115 dd = sqrt(sum(tdip(1:3)**2))*debye
3116 WRITE (unit_nr,
"(T3,A)")
"Dipole moment [Debye]"
3117 WRITE (unit_nr,
"(T5,3(A,A,F14.8,1X),T60,A,T67,F14.8)") &
3118 (trim(rlab(i)),
"=", tdip(i)*debye, i=1, 3),
"Total=", dd
3123 CALL cp_print_key_finished_output(unit_nr=unit_nr, logger=logger, &
3124 basis_section=ec_section, print_key_path=
"PRINT%MOMENTS")
3125 CALL get_qs_env(qs_env=qs_env, results=results)
3126 description =
"[DIPOLE]"
3127 CALL cp_results_erase(results=results, description=description)
3128 CALL put_results(results=results, description=description, values=tdip(1:3))
3132 print_key_voro => section_vals_get_subs_vals(ec_section,
"PRINT%VORONOI")
3133 print_key_bqb => section_vals_get_subs_vals(ec_section,
"PRINT%E_DENSITY_BQB")
3134 IF (btest(cp_print_key_should_output(logger%iter_info, print_key_voro), cp_p_file))
THEN
3135 should_print_voro = 1
3137 should_print_voro = 0
3139 IF (btest(cp_print_key_should_output(logger%iter_info, print_key_bqb), cp_p_file))
THEN
3140 should_print_bqb = 1
3142 should_print_bqb = 0
3144 IF ((should_print_voro /= 0) .OR. (should_print_bqb /= 0))
THEN
3146 CALL get_qs_env(qs_env=qs_env, &
3148 CALL pw_env_get(pw_env=pw_env, &
3149 auxbas_pw_pool=auxbas_pw_pool, &
3151 CALL auxbas_pw_pool%create_pw(pw=rho_elec_rspace)
3153 IF (dft_control%nspins > 1)
THEN
3156 CALL pw_copy(ec_env%rhoout_r(1), rho_elec_rspace)
3157 CALL pw_axpy(ec_env%rhoout_r(2), rho_elec_rspace)
3159 CALL pw_axpy(ec_env%rhoz_r(1), rho_elec_rspace)
3160 CALL pw_axpy(ec_env%rhoz_r(2), rho_elec_rspace)
3164 CALL pw_copy(ec_env%rhoout_r(1), rho_elec_rspace)
3165 CALL pw_axpy(ec_env%rhoz_r(1), rho_elec_rspace)
3168 IF (should_print_voro /= 0)
THEN
3169 CALL section_vals_val_get(print_key_voro,
"OUTPUT_TEXT", l_val=voro_print_txt)
3170 IF (voro_print_txt)
THEN
3171 append_voro = section_get_lval(ec_section,
"PRINT%VORONOI%APPEND")
3172 my_pos_voro =
"REWIND"
3173 IF (append_voro)
THEN
3174 my_pos_voro =
"APPEND"
3176 unit_nr_voro = cp_print_key_unit_nr(logger, ec_section,
"PRINT%VORONOI", extension=
".voronoi", &
3177 file_position=my_pos_voro, log_filename=.false.)
3185 CALL entry_voronoi_or_bqb(should_print_voro, should_print_bqb, print_key_voro, print_key_bqb, &
3186 unit_nr_voro, qs_env, rho_elec_rspace)
3188 CALL auxbas_pw_pool%give_back_pw(rho_elec_rspace)
3190 IF (unit_nr_voro > 0)
THEN
3191 CALL cp_print_key_finished_output(unit_nr_voro, logger, ec_section,
"PRINT%VORONOI")
3196 CALL timestop(handle)
3198 END SUBROUTINE ec_properties
3211 SUBROUTINE ec_ls_init(qs_env, matrix_ks, matrix_s)
3212 TYPE(qs_environment_type),
POINTER :: qs_env
3213 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_ks, matrix_s
3215 CHARACTER(len=*),
PARAMETER :: routinen =
'ec_ls_init'
3217 INTEGER :: handle, ispin, nspins
3218 TYPE(dft_control_type),
POINTER :: dft_control
3219 TYPE(energy_correction_type),
POINTER :: ec_env
3220 TYPE(ls_scf_env_type),
POINTER :: ls_env
3222 CALL timeset(routinen, handle)
3224 CALL get_qs_env(qs_env=qs_env, &
3225 dft_control=dft_control, &
3227 nspins = dft_control%nspins
3228 ls_env => ec_env%ls_env
3231 CALL matrix_ls_create(matrix_ls=ls_env%matrix_s, matrix_qs=matrix_s(1, 1)%matrix, &
3232 ls_mstruct=ls_env%ls_mstruct)
3234 IF (
ALLOCATED(ls_env%matrix_p))
THEN
3235 DO ispin = 1,
SIZE(ls_env%matrix_p)
3236 CALL dbcsr_release(ls_env%matrix_p(ispin))
3239 ALLOCATE (ls_env%matrix_p(nspins))
3242 DO ispin = 1, nspins
3243 CALL dbcsr_create(ls_env%matrix_p(ispin), template=ls_env%matrix_s, &
3244 matrix_type=dbcsr_type_no_symmetry)
3247 ALLOCATE (ls_env%matrix_ks(nspins))
3248 DO ispin = 1, nspins
3249 CALL dbcsr_create(ls_env%matrix_ks(ispin), template=ls_env%matrix_s, &
3250 matrix_type=dbcsr_type_no_symmetry)
3254 CALL ls_scf_init_matrix_s(matrix_s(1, 1)%matrix, ls_env)
3258 DO ispin = 1, nspins
3259 CALL matrix_qs_to_ls(matrix_ls=ls_env%matrix_ks(ispin), &
3260 matrix_qs=matrix_ks(ispin, 1)%matrix, &
3261 ls_mstruct=ls_env%ls_mstruct, &
3265 CALL timestop(handle)
3267 END SUBROUTINE ec_ls_init
3283 SUBROUTINE ec_ls_solver(qs_env, matrix_p, matrix_w, ec_ls_method)
3284 TYPE(qs_environment_type),
POINTER :: qs_env
3285 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_p, matrix_w
3286 INTEGER,
INTENT(IN) :: ec_ls_method
3288 CHARACTER(LEN=*),
PARAMETER :: routinen =
'ec_ls_solver'
3290 INTEGER :: handle, ispin, nelectron_spin_real, &
3292 INTEGER,
DIMENSION(2) :: nmo
3293 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: wmat
3294 TYPE(dbcsr_type),
ALLOCATABLE,
DIMENSION(:) :: matrix_ks_deviation
3295 TYPE(dft_control_type),
POINTER :: dft_control
3296 TYPE(energy_correction_type),
POINTER :: ec_env
3297 TYPE(ls_scf_env_type),
POINTER :: ls_env
3298 TYPE(mp_para_env_type),
POINTER :: para_env
3300 CALL timeset(routinen, handle)
3303 CALL get_qs_env(qs_env, &
3304 dft_control=dft_control, &
3306 nspins = dft_control%nspins
3307 ec_env => qs_env%ec_env
3308 ls_env => ec_env%ls_env
3311 CALL get_qs_env(qs_env=qs_env, nelectron_spin=nmo)
3312 IF (nspins == 1) nmo(1) = nmo(1)/2
3313 ls_env%homo_spin(:) = 0.0_dp
3314 ls_env%lumo_spin(:) = 0.0_dp
3316 ALLOCATE (matrix_ks_deviation(nspins))
3317 DO ispin = 1, nspins
3318 CALL dbcsr_create(matrix_ks_deviation(ispin), template=ls_env%matrix_ks(ispin))
3319 CALL dbcsr_set(matrix_ks_deviation(ispin), 0.0_dp)
3323 IF (ls_env%has_s_preconditioner)
THEN
3324 DO ispin = 1, nspins
3325 CALL apply_matrix_preconditioner(ls_env%matrix_ks(ispin),
"forward", &
3326 ls_env%matrix_bs_sqrt, ls_env%matrix_bs_sqrt_inv)
3328 CALL dbcsr_filter(ls_env%matrix_ks(ispin), ls_env%eps_filter)
3332 DO ispin = 1, nspins
3333 nelectron_spin_real = ls_env%nelectron_spin(ispin)
3334 IF (ls_env%nspins == 1) nelectron_spin_real = nelectron_spin_real/2
3336 SELECT CASE (ec_ls_method)
3337 CASE (ec_matrix_sign)
3338 CALL density_matrix_sign(ls_env%matrix_p(ispin), &
3339 ls_env%mu_spin(ispin), &
3341 ls_env%sign_method, &
3342 ls_env%sign_order, &
3343 ls_env%matrix_ks(ispin), &
3345 ls_env%matrix_s_inv, &
3346 nelectron_spin_real, &
3349 CASE (ec_matrix_trs4)
3350 CALL density_matrix_trs4( &
3351 ls_env%matrix_p(ispin), &
3352 ls_env%matrix_ks(ispin), &
3353 ls_env%matrix_s_sqrt_inv, &
3354 nelectron_spin_real, &
3355 ec_env%eps_default, &
3356 ls_env%homo_spin(ispin), &
3357 ls_env%lumo_spin(ispin), &
3358 ls_env%mu_spin(ispin), &
3359 matrix_ks_deviation=matrix_ks_deviation(ispin), &
3360 dynamic_threshold=ls_env%dynamic_threshold, &
3361 eps_lanczos=ls_env%eps_lanczos, &
3362 max_iter_lanczos=ls_env%max_iter_lanczos)
3364 CASE (ec_matrix_tc2)
3365 CALL density_matrix_tc2( &
3366 ls_env%matrix_p(ispin), &
3367 ls_env%matrix_ks(ispin), &
3368 ls_env%matrix_s_sqrt_inv, &
3369 nelectron_spin_real, &
3370 ec_env%eps_default, &
3371 ls_env%homo_spin(ispin), &
3372 ls_env%lumo_spin(ispin), &
3373 non_monotonic=ls_env%non_monotonic, &
3374 eps_lanczos=ls_env%eps_lanczos, &
3375 max_iter_lanczos=ls_env%max_iter_lanczos)
3382 IF (ls_env%has_s_preconditioner)
THEN
3383 DO ispin = 1, nspins
3385 CALL apply_matrix_preconditioner(ls_env%matrix_p(ispin),
"forward", &
3386 ls_env%matrix_bs_sqrt, ls_env%matrix_bs_sqrt_inv)
3388 CALL dbcsr_filter(ls_env%matrix_p(ispin), ls_env%eps_filter)
3393 IF (nspins == 1)
CALL dbcsr_scale(ls_env%matrix_p(1), 2.0_dp)
3395 IF (ls_env%report_all_sparsities)
CALL post_scf_sparsities(ls_env)
3399 DO ispin = 1, nspins
3400 CALL matrix_ls_to_qs(matrix_qs=matrix_p(ispin, 1)%matrix, &
3401 matrix_ls=ls_env%matrix_p(ispin), &
3402 ls_mstruct=ls_env%ls_mstruct, &
3406 wmat => matrix_w(:, 1)
3407 CALL calculate_w_matrix_ls(wmat, ec_env%ls_env)
3410 CALL dbcsr_release(ls_env%matrix_s)
3411 IF (ls_env%has_s_preconditioner)
THEN
3412 CALL dbcsr_release(ls_env%matrix_bs_sqrt)
3413 CALL dbcsr_release(ls_env%matrix_bs_sqrt_inv)
3415 IF (ls_env%needs_s_inv)
THEN
3416 CALL dbcsr_release(ls_env%matrix_s_inv)
3418 IF (ls_env%use_s_sqrt)
THEN
3419 CALL dbcsr_release(ls_env%matrix_s_sqrt)
3420 CALL dbcsr_release(ls_env%matrix_s_sqrt_inv)
3423 DO ispin = 1,
SIZE(ls_env%matrix_ks)
3424 CALL dbcsr_release(ls_env%matrix_ks(ispin))
3426 DEALLOCATE (ls_env%matrix_ks)
3428 DO ispin = 1, nspins
3429 CALL dbcsr_release(matrix_ks_deviation(ispin))
3431 DEALLOCATE (matrix_ks_deviation)
3433 CALL timestop(handle)
3435 END SUBROUTINE ec_ls_solver
3453 SUBROUTINE ec_ot_diag_solver(qs_env, ec_env, matrix_ks, matrix_s, matrix_p, matrix_w)
3454 TYPE(qs_environment_type),
POINTER :: qs_env
3455 TYPE(energy_correction_type),
POINTER :: ec_env
3456 TYPE(dbcsr_p_type),
DIMENSION(:, :),
INTENT(IN), &
3457 POINTER :: matrix_ks, matrix_s
3458 TYPE(dbcsr_p_type),
DIMENSION(:, :), &
3459 INTENT(INOUT),
POINTER :: matrix_p, matrix_w
3461 CHARACTER(len=*),
PARAMETER :: routinen =
'ec_ot_diag_solver'
3463 INTEGER :: handle, homo, ikind, iounit, ispin, &
3464 max_iter, nao, nkind, nmo, nspins
3465 INTEGER,
DIMENSION(2) :: nelectron_spin
3466 REAL(kind=dp),
DIMENSION(:),
POINTER :: eigenvalues
3467 TYPE(atomic_kind_type),
DIMENSION(:),
POINTER :: atomic_kind_set
3468 TYPE(cp_blacs_env_type),
POINTER :: blacs_env
3469 TYPE(cp_fm_type) :: sv
3470 TYPE(cp_fm_type),
POINTER :: mo_coeff
3471 TYPE(cp_logger_type),
POINTER :: logger
3472 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: p_rmpv
3473 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: rho_ao
3474 TYPE(dft_control_type),
POINTER :: dft_control
3475 TYPE(gto_basis_set_type),
POINTER :: basis_set, harris_basis
3476 TYPE(mo_set_type),
DIMENSION(:),
POINTER :: mos
3477 TYPE(mp_para_env_type),
POINTER :: para_env
3478 TYPE(particle_type),
DIMENSION(:),
POINTER :: particle_set
3479 TYPE(preconditioner_type),
POINTER :: local_preconditioner
3480 TYPE(qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
3481 TYPE(qs_kind_type),
POINTER :: qs_kind
3482 TYPE(qs_rho_type),
POINTER :: rho
3484 CALL timeset(routinen, handle)
3486 logger => cp_get_default_logger()
3487 iounit = cp_logger_get_default_unit_nr(logger)
3489 cpassert(
ASSOCIATED(qs_env))
3490 cpassert(
ASSOCIATED(ec_env))
3491 cpassert(
ASSOCIATED(matrix_ks))
3492 cpassert(
ASSOCIATED(matrix_s))
3493 cpassert(
ASSOCIATED(matrix_p))
3494 cpassert(
ASSOCIATED(matrix_w))
3496 CALL get_qs_env(qs_env=qs_env, &
3497 atomic_kind_set=atomic_kind_set, &
3498 blacs_env=blacs_env, &
3499 dft_control=dft_control, &
3500 nelectron_spin=nelectron_spin, &
3501 para_env=para_env, &
3502 particle_set=particle_set, &
3503 qs_kind_set=qs_kind_set)
3504 nspins = dft_control%nspins
3511 IF (dft_control%qs_control%do_ls_scf)
THEN
3512 CALL ec_mos_init(qs_env, matrix_s(1, 1)%matrix)
3515 CALL get_qs_env(qs_env, mos=mos)
3522 CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, nkind=nkind)
3527 IF (ec_env%basis ==
"HARRIS")
THEN
3529 qs_kind => qs_kind_set(ikind)
3531 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type=
"ORB")
3533 CALL get_qs_kind(qs_kind=qs_kind, basis_set=harris_basis, basis_type=
"HARRIS")
3535 IF (basis_set%name .NE. harris_basis%name)
THEN
3536 cpabort(
"OT-Diag initial guess: Harris and ground state need to use the same basis")
3541 IF (ec_env%mao)
THEN
3542 cpabort(
"OT-Diag initial guess: not implemented for MAOs")
3546 SELECT CASE (ec_env%ec_initial_guess)
3549 p_rmpv => matrix_p(:, 1)
3551 CALL calculate_atomic_block_dm(p_rmpv, matrix_s(1, 1)%matrix, atomic_kind_set, qs_kind_set, &
3552 nspins, nelectron_spin, iounit, para_env)
3556 CALL get_qs_env(qs_env, rho=rho)
3557 CALL qs_rho_get(rho, rho_ao_kp=rho_ao)
3558 p_rmpv => rho_ao(:, 1)
3561 cpabort(
"Unknown inital guess for OT-Diagonalization (Harris functional)")
3564 DO ispin = 1, nspins
3565 CALL get_mo_set(mo_set=mos(ispin), &
3566 mo_coeff=mo_coeff, &
3572 CALL cp_fm_set_all(mo_coeff, 0.0_dp)
3573 CALL cp_fm_init_random(mo_coeff, nmo)
3575 CALL cp_fm_create(sv, mo_coeff%matrix_struct,
"SV")
3578 CALL cp_dbcsr_sm_fm_multiply(matrix_s(1, 1)%matrix, mo_coeff, sv, nmo)
3579 CALL cp_dbcsr_sm_fm_multiply(p_rmpv(ispin)%matrix, sv, mo_coeff, homo)
3580 CALL cp_fm_release(sv)
3583 CALL make_basis_sm(mo_coeff, nmo, matrix_s(1, 1)%matrix)
3587 NULLIFY (local_preconditioner)
3588 ALLOCATE (local_preconditioner)
3589 CALL init_preconditioner(local_preconditioner, para_env=para_env, &
3590 blacs_env=blacs_env)
3591 DO ispin = 1, nspins
3592 CALL make_preconditioner(local_preconditioner, &
3593 precon_type=ot_precond_full_single_inverse, &
3594 solver_type=ot_precond_solver_default, &
3595 matrix_h=matrix_ks(ispin, 1)%matrix, &
3596 matrix_s=matrix_s(ispin, 1)%matrix, &
3597 convert_precond_to_dbcsr=.true., &
3598 mo_set=mos(ispin), energy_gap=0.2_dp)
3600 CALL get_mo_set(mos(ispin), &
3601 mo_coeff=mo_coeff, &
3602 eigenvalues=eigenvalues, &
3605 CALL ot_eigensolver(matrix_h=matrix_ks(ispin, 1)%matrix, &
3606 matrix_s=matrix_s(1, 1)%matrix, &
3607 matrix_c_fm=mo_coeff, &
3609 eps_gradient=ec_env%eps_default, &
3610 iter_max=max_iter, &
3612 CALL calculate_subspace_eigenvalues(mo_coeff, matrix_ks(ispin, 1)%matrix, &
3613 evals_arg=eigenvalues, do_rotation=.true.)
3616 CALL destroy_preconditioner(local_preconditioner)
3617 DEALLOCATE (local_preconditioner)
3620 CALL copy_fm_to_dbcsr(mos(ispin)%mo_coeff, &
3621 mos(ispin)%mo_coeff_b)
3625 DO ispin = 1, nspins
3626 CALL calculate_density_matrix(mos(ispin), matrix_p(ispin, 1)%matrix)
3628 CALL calculate_w_matrix(mos(ispin), matrix_w(ispin, 1)%matrix)
3632 IF (dft_control%qs_control%do_ls_scf)
THEN
3633 DO ispin = 1, nspins
3634 CALL deallocate_mo_set(mos(ispin))
3636 IF (
ASSOCIATED(qs_env%mos))
THEN
3637 DO ispin = 1,
SIZE(qs_env%mos)
3638 CALL deallocate_mo_set(qs_env%mos(ispin))
3640 DEALLOCATE (qs_env%mos)
3644 CALL timestop(handle)
3646 END SUBROUTINE ec_ot_diag_solver
real(kind=dp) function det_3x3(a)
...
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
subroutine uppercase(string)
...
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)
...
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.
Calculation of the nuclear attraction contribution to the core Hamiltonian <a|erfc|b> :we only calcul...
subroutine, public build_core_ae(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ae, nimages, cell_to_index)
...
Calculation of the local pseudopotential contribution to the core Hamiltonian <a|V(local)|b> = <a|Sum...
subroutine, public build_core_ppl(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ppl, nimages, cell_to_index, basis_type, deltaR)
...
Calculation of the non-local pseudopotential contribution to the core Hamiltonian <a|V(non-local)|b> ...
subroutine, public build_core_ppnl(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, qs_kind_set, atomic_kind_set, particle_set, sab_orb, sap_ppnl, eps_ppnl, nimages, cell_to_index, basis_type, deltaR, matrix_l)
...
methods related to the blacs parallel environment
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
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.
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
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 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_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)
compute the density matrix using a trace-resetting algorithm
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)
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 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.
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_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
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, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, 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_r3d_rs_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_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)
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.
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
subroutine, public write_stress_tensor_components(virial, iw, cell)
...
subroutine, public write_stress_tensor(pv_virial, iw, cell, numerical)
Print stress tensor to output file.
pure real(kind=dp) function, public one_third_sum_diag(a)
...
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.