comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-28 22:25:32 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Refactioring of kinetic energy and core matrices calls. (#4393)

Browse source
This commit is contained in:
Juerg Hutter 2025-08-26 11:22:41 +02:00 • committed by GitHub
parent 448620dc73
commit e3076134ed
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
15 changed files with 752 additions and 829 deletions
Download patch file Download diff file Expand all files Collapse all files

1
src/CMakeLists.txt
View file

@ -530,6 +530,7 @@ list(
qs_condnum.F
qs_core_energies.F
qs_core_hamiltonian.F
qs_core_matrices.F
qs_dcdr_ao.F
qs_dcdr.F
qs_dcdr_utils.F

4
src/core_ppl.F
View file

@ -168,6 +168,10 @@ CONTAINS
do_dR = PRESENT(deltaR)
doat = PRESENT(atcore)
IF ((calculate_forces .OR. doat) .AND. do_dR) THEN
CPABORT("core_ppl: incompatible options")
END IF
MARK_USED(int_8)
! Use internal integral routine for local ECP terms or use libgrrp

3
src/core_ppnl.F
View file

@ -162,6 +162,9 @@ CONTAINS
IF (PRESENT(deltaR)) do_dR = .TRUE.
doat = .FALSE.
IF (PRESENT(atcore)) doat = .TRUE.
IF ((calculate_forces .OR. doat) .AND. do_dR) THEN
CPABORT("core_ppl: incompatible options")
END IF
IF (calculate_forces) THEN
CALL timeset(routineN//"_forces", handle)

131
src/ec_external.F
View file

@ -12,11 +12,7 @@
!> \author JGH
! **************************************************************************************************
MODULE ec_external
USE atomic_kind_types, ONLY: atomic_kind_type
USE cell_types, ONLY: cell_type
USE core_ae, ONLY: build_core_ae
USE core_ppl, ONLY: build_core_ppl
USE core_ppnl, ONLY: build_core_ppnl
USE cp_control_types, ONLY: dft_control_type
USE cp_dbcsr_api, ONLY: dbcsr_add,&
dbcsr_copy,&
@ -50,14 +46,13 @@ MODULE ec_external
USE mathlib, ONLY: det_3x3
USE message_passing, ONLY: mp_para_env_type
USE parallel_gemm_api, ONLY: parallel_gemm
USE particle_types, ONLY: particle_type
USE physcon, ONLY: pascal
USE qs_core_energies, ONLY: calculate_ptrace
USE qs_core_matrices, ONLY: core_matrices,&
kinetic_energy_matrix
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
USE qs_force_types, ONLY: qs_force_type
USE qs_kind_types, ONLY: qs_kind_type
USE qs_kinetic, ONLY: build_kinetic_matrix
USE qs_ks_types, ONLY: qs_ks_env_type
USE qs_mo_types, ONLY: deallocate_mo_set,&
get_mo_set,&
@ -361,23 +356,16 @@ CONTAINS
CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_debug_force'
INTEGER :: handle, iounit, nder, nimages
INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
LOGICAL :: calculate_forces, debug_forces, &
debug_stress, use_virial
REAL(KIND=dp) :: eps_ppnl, fconv
REAL(KIND=dp), DIMENSION(3) :: fodeb
REAL(KIND=dp) :: fconv
REAL(KIND=dp), DIMENSION(3, 3) :: stdeb, sttot
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(cell_type), POINTER :: cell
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: scrm
TYPE(dft_control_type), POINTER :: dft_control
TYPE(mp_para_env_type), POINTER :: para_env
TYPE(neighbor_list_set_p_type), DIMENSION(:), &
POINTER :: sab_orb, sac_ae, sac_ppl, sap_ppnl
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(qs_force_type), DIMENSION(:), POINTER :: force
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(qs_ks_env_type), POINTER :: ks_env
POINTER :: sab_orb
TYPE(virial_type), POINTER :: virial
CALL timeset(routineN, handle)
@ -389,12 +377,9 @@ CONTAINS
calculate_forces = .TRUE.
! no k-points possible
NULLIFY (cell, dft_control, force, ks_env, para_env, virial)
CALL get_qs_env(qs_env=qs_env, &
cell=cell, &
dft_control=dft_control, &
force=force, &
ks_env=ks_env, &
para_env=para_env, &
virial=virial)
nimages = dft_control%nimages
@ -415,116 +400,22 @@ CONTAINS
CPABORT("GAPW not implemented")
END IF
! get neighbor lists
NULLIFY (sab_orb, sac_ae, sac_ppl, sap_ppnl)
CALL get_qs_env(qs_env=qs_env, &
sab_orb=sab_orb, sac_ae=sac_ae, &
sac_ppl=sac_ppl, sap_ppnl=sap_ppnl)
! initialize src matrix
NULLIFY (scrm)
CALL dbcsr_allocate_matrix_set(scrm, 1, 1)
ALLOCATE (scrm(1, 1)%matrix)
CALL dbcsr_create(scrm(1, 1)%matrix, template=matrix_p(1, 1)%matrix)
CALL get_qs_env(qs_env=qs_env, sab_orb=sab_orb)
CALL cp_dbcsr_alloc_block_from_nbl(scrm(1, 1)%matrix, sab_orb)
nder = 1
IF (SIZE(matrix_p, 1) == 2) THEN
CALL dbcsr_add(matrix_p(1, 1)%matrix, matrix_p(2, 1)%matrix, &
alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
END IF
! kinetic energy
IF (debug_forces) fodeb(1:3) = force(1)%kinetic(1:3, 1)
IF (debug_stress .AND. use_virial) stdeb = virial%pv_ekinetic
CALL build_kinetic_matrix(ks_env, matrixkp_t=scrm, &
matrix_name="KINETIC ENERGY MATRIX", &
basis_type="ORB", &
sab_nl=sab_orb, calculate_forces=.TRUE., &
matrixkp_p=matrix_p)
IF (debug_forces) THEN
fodeb(1:3) = force(1)%kinetic(1:3, 1) - fodeb(1:3)
CALL para_env%sum(fodeb)
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pout*dT ", fodeb
END IF
IF (debug_stress .AND. use_virial) THEN
stdeb = fconv*(virial%pv_ekinetic - stdeb)
CALL para_env%sum(stdeb)
IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
'STRESS| Pout*dT', one_third_sum_diag(stdeb), det_3x3(stdeb)
END IF
IF (SIZE(matrix_p, 1) == 2) THEN
CALL dbcsr_add(matrix_p(1, 1)%matrix, matrix_p(2, 1)%matrix, &
alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
END IF
CALL kinetic_energy_matrix(qs_env, matrixkp_t=scrm, matrix_p=matrix_p, &
calculate_forces=calculate_forces, &
debug_forces=debug_forces, debug_stress=debug_stress)
! compute the ppl contribution to the core hamiltonian
NULLIFY (atomic_kind_set, particle_set, qs_kind_set)
CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, particle_set=particle_set, &
atomic_kind_set=atomic_kind_set)
IF (ASSOCIATED(sac_ae)) THEN
IF (calculate_forces .AND. debug_forces) fodeb(1:3) = force(1)%all_potential(1:3, 1)
IF (debug_stress .AND. use_virial) stdeb = virial%pv_ppl
CALL build_core_ae(scrm, 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, "ORB")
IF (calculate_forces .AND. debug_forces) THEN
fodeb(1:3) = force(1)%all_potential(1:3, 1) - fodeb(1:3)
CALL para_env%sum(fodeb)
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pout*dH_AE ", fodeb
END IF
IF (debug_stress .AND. use_virial) THEN
stdeb = fconv*(virial%pv_ppl - stdeb)
CALL para_env%sum(stdeb)
IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
'STRESS| Pout*dH_AE ', one_third_sum_diag(stdeb), det_3x3(stdeb)
END IF
END IF
IF (ASSOCIATED(sac_ppl)) THEN
IF (calculate_forces .AND. debug_forces) fodeb(1:3) = force(1)%gth_ppl(1:3, 1)
IF (debug_stress .AND. use_virial) stdeb = virial%pv_ppl
CALL build_core_ppl(scrm, 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, "ORB")
IF (calculate_forces .AND. debug_forces) THEN
fodeb(1:3) = force(1)%gth_ppl(1:3, 1) - fodeb(1:3)
CALL para_env%sum(fodeb)
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pout*dH_PPL ", fodeb
END IF
IF (debug_stress .AND. use_virial) THEN
stdeb = fconv*(virial%pv_ppl - stdeb)
CALL para_env%sum(stdeb)
IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
'STRESS| Pout*dH_PPL', one_third_sum_diag(stdeb), det_3x3(stdeb)
END IF
END IF
! compute the ppnl contribution to the core hamiltonian ***
eps_ppnl = dft_control%qs_control%eps_ppnl
IF (ASSOCIATED(sap_ppnl)) THEN
IF (calculate_forces .AND. debug_forces) fodeb(1:3) = force(1)%gth_ppnl(1:3, 1)
IF (debug_stress .AND. use_virial) stdeb = virial%pv_ppnl
CALL build_core_ppnl(scrm, 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, "ORB")
IF (calculate_forces .AND. debug_forces) THEN
fodeb(1:3) = force(1)%gth_ppnl(1:3, 1) - fodeb(1:3)
CALL para_env%sum(fodeb)
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pout*dH_PPNL", fodeb
END IF
IF (debug_stress .AND. use_virial) THEN
stdeb = fconv*(virial%pv_ppnl - stdeb)
CALL para_env%sum(stdeb)
IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
'STRESS| Pout*dH_PPNL', one_third_sum_diag(stdeb), det_3x3(stdeb)
END IF
END IF
nder = 1
CALL core_matrices(qs_env, scrm, matrix_p, calculate_forces, nder, &
debug_forces=debug_forces, debug_stress=debug_stress)
IF (debug_stress .AND. use_virial) THEN
stdeb = fconv*(virial%pv_virial - sttot)

2
src/ed_analysis.F
View file

@ -72,7 +72,7 @@ MODULE ed_analysis
USE qs_core_energies, ONLY: calculate_ecore_alpha,&
calculate_ecore_overlap,&
calculate_ecore_self
USE qs_core_hamiltonian, ONLY: core_matrices
USE qs_core_matrices, ONLY: core_matrices
USE qs_dispersion_pairpot, ONLY: calculate_dispersion_pairpot
USE qs_dispersion_types, ONLY: qs_dispersion_type
USE qs_energy_types, ONLY: qs_energy_type

163
src/energy_corrections.F
View file

@ -26,9 +26,6 @@ MODULE energy_corrections
cite_reference
USE cell_types, ONLY: cell_type,&
pbc
USE core_ae, ONLY: build_core_ae
USE core_ppl, ONLY: build_core_ppl
USE core_ppnl, ONLY: build_core_ppnl
USE cp_blacs_env, ONLY: cp_blacs_env_type
USE cp_control_types, ONLY: dft_control_type
USE cp_dbcsr_api, ONLY: &
@ -142,6 +139,8 @@ MODULE energy_corrections
USE qs_collocate_density, ONLY: calculate_rho_elec
USE qs_core_energies, ONLY: calculate_ecore_overlap,&
calculate_ptrace
USE qs_core_matrices, ONLY: core_matrices,&
kinetic_energy_matrix
USE qs_density_matrices, ONLY: calculate_density_matrix,&
calculate_w_matrix
USE qs_dispersion_pairpot, ONLY: calculate_dispersion_pairpot
@ -1322,23 +1321,21 @@ CONTAINS
CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_build_core_hamiltonian'
CHARACTER(LEN=default_string_length) :: basis_type
INTEGER :: handle, nder, nimages
INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
LOGICAL :: calculate_forces, use_virial
REAL(KIND=dp) :: eps_ppnl
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(dft_control_type), POINTER :: dft_control
TYPE(neighbor_list_set_p_type), DIMENSION(:), &
POINTER :: sab_orb, sac_ae, sac_ppl, sap_ppnl
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(qs_force_type), DIMENSION(:), POINTER :: force
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(qs_ks_env_type), POINTER :: ks_env
TYPE(virial_type), POINTER :: virial
CALL timeset(routineN, handle)
NULLIFY (atomic_kind_set, cell_to_index, dft_control, ks_env, particle_set, qs_kind_set, virial)
NULLIFY (atomic_kind_set, dft_control, ks_env, particle_set, &
qs_kind_set)
CALL get_qs_env(qs_env=qs_env, &
atomic_kind_set=atomic_kind_set, &
@ -1369,16 +1366,18 @@ CONTAINS
sac_ppl => ec_env%sac_ppl
sap_ppnl => ec_env%sap_ppnl
basis_type = "HARRIS"
nder = 0
! Overlap and kinetic energy matrices
CALL build_overlap_matrix(ks_env, matrixkp_s=ec_env%matrix_s, &
matrix_name="OVERLAP MATRIX", &
basis_type_a="HARRIS", &
basis_type_b="HARRIS", &
basis_type_a=basis_type, &
basis_type_b=basis_type, &
sab_nl=sab_orb)
CALL build_kinetic_matrix(ks_env, matrixkp_t=ec_env%matrix_t, &
matrix_name="KINETIC ENERGY MATRIX", &
basis_type="HARRIS", &
basis_type=basis_type, &
sab_nl=sab_orb)
! initialize H matrix
@ -1391,29 +1390,8 @@ CONTAINS
CALL dbcsr_copy(ec_env%matrix_h(1, 1)%matrix, ec_env%matrix_t(1, 1)%matrix, &
keep_sparsity=.TRUE., name="CORE HAMILTONIAN MATRIX")
! Possible AE contributions (also ECP)
IF (ASSOCIATED(sac_ae)) THEN
CALL build_core_ae(ec_env%matrix_h, ec_env%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, "HARRIS")
END IF
! compute the ppl contribution to the core hamiltonian
IF (ASSOCIATED(sac_ppl)) THEN
CALL build_core_ppl(ec_env%matrix_h, ec_env%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, "HARRIS")
END IF
! compute the ppnl contribution to the core hamiltonian ***
eps_ppnl = dft_control%qs_control%eps_ppnl
IF (ASSOCIATED(sap_ppnl)) THEN
CALL build_core_ppnl(ec_env%matrix_h, ec_env%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, "HARRIS")
END IF
CALL core_matrices(qs_env, ec_env%matrix_h, ec_env%matrix_p, calculate_forces, nder, &
ec_env=ec_env, ec_env_matrices=.TRUE., basis_type=basis_type)
! External field (nonperiodic case)
ec_env%efield_nuclear = 0.0_dp
@ -1660,24 +1638,21 @@ CONTAINS
CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_build_core_hamiltonian_force'
CHARACTER(LEN=default_string_length) :: basis_type
INTEGER :: handle, iounit, nder, nimages
INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
LOGICAL :: calculate_forces, debug_forces, &
debug_stress, use_virial
REAL(KIND=dp) :: eps_ppnl, fconv
REAL(KIND=dp) :: fconv
REAL(KIND=dp), DIMENSION(3) :: fodeb
REAL(KIND=dp), DIMENSION(3, 3) :: stdeb, sttot
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(cell_type), POINTER :: cell
TYPE(cp_logger_type), POINTER :: logger
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: scrm
TYPE(dft_control_type), POINTER :: dft_control
TYPE(mp_para_env_type), POINTER :: para_env
TYPE(neighbor_list_set_p_type), DIMENSION(:), &
POINTER :: sab_orb, sac_ae, sac_ppl, sap_ppnl
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
POINTER :: sab_orb
TYPE(qs_force_type), DIMENSION(:), POINTER :: force
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(qs_ks_env_type), POINTER :: ks_env
TYPE(virial_type), POINTER :: virial
@ -1695,6 +1670,8 @@ CONTAINS
calculate_forces = .TRUE.
basis_type = "HARRIS"
! no k-points possible
NULLIFY (cell, dft_control, force, ks_env, para_env, virial)
CALL get_qs_env(qs_env=qs_env, &
@ -1725,11 +1702,7 @@ CONTAINS
END IF
! get neighbor lists, we need the full sab_orb list from the ec_env
NULLIFY (sab_orb, sac_ppl, sap_ppnl)
sab_orb => ec_env%sab_orb
sac_ae => ec_env%sac_ae
sac_ppl => ec_env%sac_ppl
sap_ppnl => ec_env%sap_ppnl
! initialize src matrix
NULLIFY (scrm)
@ -1740,8 +1713,6 @@ CONTAINS
nder = 1
IF (SIZE(matrix_p, 1) == 2) THEN
CALL dbcsr_add(matrix_p(1, 1)%matrix, matrix_p(2, 1)%matrix, &
alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
CALL dbcsr_add(matrix_w(1, 1)%matrix, matrix_w(2, 1)%matrix, &
alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
END IF
@ -1751,8 +1722,8 @@ CONTAINS
IF (debug_stress .AND. use_virial) stdeb = virial%pv_overlap
CALL build_overlap_matrix(ks_env, matrixkp_s=scrm, &
matrix_name="OVERLAP MATRIX", &
basis_type_a="HARRIS", &
basis_type_b="HARRIS", &
basis_type_a=basis_type, &
basis_type_b=basis_type, &
sab_nl=sab_orb, calculate_forces=.TRUE., &
matrixkp_p=matrix_w)
@ -1760,103 +1731,23 @@ CONTAINS
fodeb(1:3) = force(1)%overlap(1:3, 1) - fodeb(1:3)
CALL para_env%sum(fodeb)
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Wout*dS ", fodeb
fodeb(1:3) = force(1)%kinetic(1:3, 1)
END IF
IF (debug_stress .AND. use_virial) THEN
stdeb = fconv*(virial%pv_overlap - stdeb)
CALL para_env%sum(stdeb)
IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
'STRESS| Wout*dS', one_third_sum_diag(stdeb), det_3x3(stdeb)
stdeb = virial%pv_ekinetic
END IF
CALL build_kinetic_matrix(ks_env, matrixkp_t=scrm, &
matrix_name="KINETIC ENERGY MATRIX", &
basis_type="HARRIS", &
sab_nl=sab_orb, calculate_forces=.TRUE., &
matrixkp_p=matrix_p)
IF (debug_forces) THEN
fodeb(1:3) = force(1)%kinetic(1:3, 1) - fodeb(1:3)
CALL para_env%sum(fodeb)
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pout*dT ", fodeb
END IF
IF (debug_stress .AND. use_virial) THEN
stdeb = fconv*(virial%pv_ekinetic - stdeb)
CALL para_env%sum(stdeb)
IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
'STRESS| Pout*dT', one_third_sum_diag(stdeb), det_3x3(stdeb)
END IF
IF (SIZE(matrix_p, 1) == 2) THEN
CALL dbcsr_add(matrix_p(1, 1)%matrix, matrix_p(2, 1)%matrix, &
alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
END IF
! compute the ppl contribution to the core hamiltonian
NULLIFY (atomic_kind_set, particle_set, qs_kind_set)
CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, particle_set=particle_set, &
atomic_kind_set=atomic_kind_set)
CALL kinetic_energy_matrix(qs_env, matrixkp_t=scrm, matrix_p=matrix_p, &
calculate_forces=.TRUE., &
sab_orb=sab_orb, &
basis_type=basis_type, &
debug_forces=debug_forces, debug_stress=debug_stress)
IF (ASSOCIATED(sac_ae)) THEN
IF (calculate_forces .AND. debug_forces) fodeb(1:3) = force(1)%all_potential(1:3, 1)
IF (debug_stress .AND. use_virial) stdeb = virial%pv_ppl
CALL build_core_ae(scrm, 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, "HARRIS")
IF (calculate_forces .AND. debug_forces) THEN
fodeb(1:3) = force(1)%all_potential(1:3, 1) - fodeb(1:3)
CALL para_env%sum(fodeb)
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pout*dHae ", fodeb
END IF
IF (debug_stress .AND. use_virial) THEN
stdeb = fconv*(virial%pv_ppl - stdeb)
CALL para_env%sum(stdeb)
IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
'STRESS| Pout*dHae ', one_third_sum_diag(stdeb), det_3x3(stdeb)
END IF
END IF
IF (ASSOCIATED(sac_ppl)) THEN
IF (calculate_forces .AND. debug_forces) fodeb(1:3) = force(1)%gth_ppl(1:3, 1)
IF (debug_stress .AND. use_virial) stdeb = virial%pv_ppl
CALL build_core_ppl(scrm, 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, "HARRIS")
IF (calculate_forces .AND. debug_forces) THEN
fodeb(1:3) = force(1)%gth_ppl(1:3, 1) - fodeb(1:3)
CALL para_env%sum(fodeb)
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pout*dH_PPL ", fodeb
END IF
IF (debug_stress .AND. use_virial) THEN
stdeb = fconv*(virial%pv_ppl - stdeb)
CALL para_env%sum(stdeb)
IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
'STRESS| Pout*dH_PPL', one_third_sum_diag(stdeb), det_3x3(stdeb)
END IF
END IF
! compute the ppnl contribution to the core hamiltonian ***
eps_ppnl = dft_control%qs_control%eps_ppnl
IF (ASSOCIATED(sap_ppnl)) THEN
IF (calculate_forces .AND. debug_forces) fodeb(1:3) = force(1)%gth_ppnl(1:3, 1)
IF (debug_stress .AND. use_virial) stdeb = virial%pv_ppnl
CALL build_core_ppnl(scrm, 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, "HARRIS")
IF (calculate_forces .AND. debug_forces) THEN
fodeb(1:3) = force(1)%gth_ppnl(1:3, 1) - fodeb(1:3)
CALL para_env%sum(fodeb)
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pout*dH_PPNL", fodeb
END IF
IF (debug_stress .AND. use_virial) THEN
stdeb = fconv*(virial%pv_ppnl - stdeb)
CALL para_env%sum(stdeb)
IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
'STRESS| Pout*dH_PPNL', one_third_sum_diag(stdeb), det_3x3(stdeb)
END IF
END IF
CALL core_matrices(qs_env, scrm, matrix_p, calculate_forces, nder, &
ec_env=ec_env, ec_env_matrices=.FALSE., basis_type=basis_type, &
debug_forces=debug_forces, debug_stress=debug_stress)
! External field (nonperiodic case)
ec_env%efield_nuclear = 0.0_dp

120
src/mp2_cphf.F
View file

@ -15,11 +15,7 @@ MODULE mp2_cphf
USE admm_methods, ONLY: admm_projection_derivative
USE admm_types, ONLY: admm_type,&
get_admm_env
USE atomic_kind_types, ONLY: atomic_kind_type
USE cell_types, ONLY: cell_type
USE core_ae, ONLY: build_core_ae
USE core_ppl, ONLY: build_core_ppl
USE core_ppnl, ONLY: build_core_ppnl
USE cp_blacs_env, ONLY: cp_blacs_env_type
USE cp_control_types, ONLY: dft_control_type
USE cp_dbcsr_api, ONLY: dbcsr_add,&
@ -74,7 +70,6 @@ MODULE mp2_cphf
USE mp2_types, ONLY: mp2_type,&
ri_rpa_method_gpw
USE parallel_gemm_api, ONLY: parallel_gemm
USE particle_types, ONLY: particle_type
USE pw_env_types, ONLY: pw_env_get,&
pw_env_type
USE pw_methods, ONLY: pw_axpy,&
@ -90,6 +85,8 @@ MODULE mp2_cphf
USE pw_types, ONLY: pw_c1d_gs_type,&
pw_r3d_rs_type
USE qs_2nd_kernel_ao, ONLY: apply_2nd_order_kernel
USE qs_core_matrices, ONLY: core_matrices,&
kinetic_energy_matrix
USE qs_density_matrices, ONLY: calculate_whz_matrix
USE qs_dispersion_pairpot, ONLY: calculate_dispersion_pairpot
USE qs_dispersion_types, ONLY: qs_dispersion_type
@ -103,8 +100,6 @@ MODULE mp2_cphf
zero_qs_force
USE qs_integrate_potential, ONLY: integrate_v_core_rspace,&
integrate_v_rspace
USE qs_kind_types, ONLY: qs_kind_type
USE qs_kinetic, ONLY: build_kinetic_matrix
USE qs_ks_reference, ONLY: ks_ref_potential
USE qs_ks_types, ONLY: qs_ks_env_type,&
set_ks_env
@ -1291,13 +1286,11 @@ CONTAINS
INTEGER :: alpha, beta, handle, idir, iounit, &
ispin, nimages, nocc, nspins
INTEGER, DIMENSION(3) :: comp
INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
LOGICAL :: do_exx, do_hfx, use_virial
REAL(KIND=dp) :: e_dummy, e_hartree, e_xc, ehartree, exc
REAL(KIND=dp), DIMENSION(3) :: deb
REAL(KIND=dp), DIMENSION(3, 3) :: h_stress, pv_virial
TYPE(admm_type), POINTER :: admm_env
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(cell_type), POINTER :: cell
TYPE(cp_logger_type), POINTER :: logger
TYPE(dbcsr_p_type), ALLOCATABLE, DIMENSION(:), &
@ -1305,15 +1298,14 @@ CONTAINS
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_p_mp2, &
matrix_p_mp2_admm, matrix_s, rho1, &
rho_ao, rho_ao_aux, scrm
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp, scrm_kp
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, rho_ao_kp, scrm_kp
TYPE(dft_control_type), POINTER :: dft_control
TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
TYPE(linres_control_type), POINTER :: linres_control
TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
TYPE(mp_para_env_type), POINTER :: para_env
TYPE(neighbor_list_set_p_type), DIMENSION(:), &
POINTER :: sab_orb, sac_ae, sac_ppl, sap_ppnl
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
POINTER :: sab_orb
TYPE(pw_c1d_gs_type) :: pot_g, rho_tot_g, temp_pw_g
TYPE(pw_c1d_gs_type), ALLOCATABLE, DIMENSION(:) :: dvg
TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g, rho_mp2_g
@ -1328,7 +1320,6 @@ CONTAINS
TYPE(qs_dispersion_type), POINTER :: dispersion_env
TYPE(qs_energy_type), POINTER :: energy
TYPE(qs_force_type), DIMENSION(:), POINTER :: force
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(qs_ks_env_type), POINTER :: ks_env
TYPE(qs_p_env_type), POINTER :: p_env
TYPE(qs_rho_type), POINTER :: rho, rho_aux
@ -1373,7 +1364,6 @@ CONTAINS
nimages = dft_control%nimages
CPASSERT(nimages == 1)
NULLIFY (cell_to_index)
p_env => qs_env%mp2_env%ri_grad%p_env
@ -1391,101 +1381,17 @@ CONTAINS
CALL dbcsr_add(rho_ao(ispin)%matrix, matrix_p_mp2(ispin)%matrix, 1.0_dp, 1.0_dp)
END DO
IF (nspins == 2) THEN
CALL dbcsr_add(rho_ao(1)%matrix, rho_ao(2)%matrix, 1.0_dp, 1.0_dp)
END IF
! kinetic energy
NULLIFY (scrm_kp)
matrix_p(1:nspins, 1:1) => rho_ao(1:nspins)
CALL kinetic_energy_matrix(qs_env, matrixkp_t=scrm_kp, matrix_p=matrix_p, &
calculate_forces=.TRUE., &
debug_forces=debug_forces, debug_stress=debug_forces)
CALL dbcsr_deallocate_matrix_set(scrm_kp)
! Kinetic energy matrix
NULLIFY (scrm)
IF (debug_forces) THEN
deb(1:3) = force(1)%kinetic(1:3, 1)
IF (use_virial) e_dummy = third_tr(virial%pv_virial)
END IF
CALL build_kinetic_matrix(ks_env, matrix_t=scrm, &
matrix_name="KINETIC ENERGY MATRIX", &
basis_type="ORB", &
sab_nl=sab_orb, calculate_forces=.TRUE., &
matrix_p=rho_ao(1)%matrix)
IF (debug_forces) THEN
deb(1:3) = force(1)%kinetic(1:3, 1) - deb(1:3)
CALL para_env%sum(deb)
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dT ", deb
IF (use_virial) THEN
e_dummy = third_tr(virial%pv_virial) - e_dummy
CALL para_env%sum(e_dummy)
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: T ", e_dummy
END IF
END IF
CALL dbcsr_deallocate_matrix_set(scrm)
IF (nspins == 2) THEN
CALL dbcsr_add(rho_ao(1)%matrix, rho_ao(2)%matrix, 1.0_dp, -1.0_dp)
END IF
! Add pseudo potential terms
scrm_kp(1:nspins, 1:1) => matrix_ks(1:nspins)
CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, particle_set=particle_set, &
atomic_kind_set=atomic_kind_set, sac_ae=sac_ae, sac_ppl=sac_ppl, sap_ppnl=sap_ppnl)
IF (ASSOCIATED(sac_ae)) THEN
IF (debug_forces) THEN
deb(1:3) = force(1)%all_potential(1:3, 1)
IF (use_virial) e_dummy = third_tr(virial%pv_virial)
END IF
CALL build_core_ae(scrm_kp, rho_ao_kp, force, &
virial, .TRUE., use_virial, 1, &
qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ae, &
nimages, cell_to_index, basis_type="ORB")
IF (debug_forces) THEN
deb(1:3) = force(1)%all_potential(1:3, 1) - deb(1:3)
CALL para_env%sum(deb)
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dHae ", deb
IF (use_virial) THEN
e_dummy = third_tr(virial%pv_virial) - e_dummy
CALL para_env%sum(e_dummy)
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: Hae ", e_dummy
END IF
END IF
END IF
IF (ASSOCIATED(sac_ppl)) THEN
IF (debug_forces) THEN
deb(1:3) = force(1)%gth_ppl(1:3, 1)
IF (use_virial) e_dummy = third_tr(virial%pv_virial)
END IF
CALL build_core_ppl(scrm_kp, rho_ao_kp, force, &
virial, .TRUE., use_virial, 1, &
qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ppl, &
nimages=1, basis_type="ORB")
IF (debug_forces) THEN
deb(1:3) = force(1)%gth_ppl(1:3, 1) - deb(1:3)
CALL para_env%sum(deb)
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dHppl ", deb
IF (use_virial) THEN
e_dummy = third_tr(virial%pv_virial) - e_dummy
CALL para_env%sum(e_dummy)
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: Hppl ", e_dummy
END IF
END IF
END IF
IF (ASSOCIATED(sap_ppnl)) THEN
IF (debug_forces) THEN
deb(1:3) = force(1)%gth_ppnl(1:3, 1)
IF (use_virial) e_dummy = third_tr(virial%pv_virial)
END IF
CALL build_core_ppnl(scrm_kp, rho_ao_kp, force, &
virial, .TRUE., use_virial, 1, &
qs_kind_set, atomic_kind_set, particle_set, sab_orb, sap_ppnl, dft_control%qs_control%eps_ppnl, &
nimages=1, basis_type="ORB")
IF (debug_forces) THEN
deb(1:3) = force(1)%gth_ppnl(1:3, 1) - deb(1:3)
CALL para_env%sum(deb)
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dHppnl ", deb
IF (use_virial) THEN
e_dummy = third_tr(virial%pv_virial) - e_dummy
CALL para_env%sum(e_dummy)
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: Hppnl ", e_dummy
END IF
END IF
END IF
CALL core_matrices(qs_env, scrm_kp, rho_ao_kp, .TRUE., 1, &
debug_forces=debug_forces, debug_stress=debug_forces)
! Get the different components of the KS potential
NULLIFY (vxc_rspace, vtau_rspace, vadmm_rspace)

264
src/qs_core_hamiltonian.F
View file

@ -75,17 +75,17 @@
!> \author Matthias Krack (14.09.2000,21.03.02)
! **************************************************************************************************
MODULE qs_core_hamiltonian
USE atomic_kind_types, ONLY: atomic_kind_type,&
get_atomic_kind_set
USE core_ae, ONLY: build_core_ae
USE core_ppl, ONLY: build_core_ppl
USE core_ppnl, ONLY: build_core_ppnl
USE atomic_kind_types, ONLY: atomic_kind_type
USE cp_blacs_env, ONLY: cp_blacs_env_type
USE cp_control_types, ONLY: dft_control_type
USE cp_dbcsr_api, ONLY: &
dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_distribution_type, dbcsr_iterator_blocks_left, &
dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
dbcsr_p_type, dbcsr_set, dbcsr_type, dbcsr_type_antisymmetric
USE cp_dbcsr_api, ONLY: dbcsr_add,&
dbcsr_copy,&
dbcsr_create,&
dbcsr_distribution_type,&
dbcsr_p_type,&
dbcsr_set,&
dbcsr_type,&
dbcsr_type_antisymmetric
USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set,&
dbcsr_deallocate_matrix_set
@ -98,28 +98,22 @@ MODULE qs_core_hamiltonian
cp_print_key_should_output,&
cp_print_key_unit_nr
USE input_constants, ONLY: do_admm_purify_none,&
do_ppl_analytic,&
kg_tnadd_atomic,&
rel_none,&
rel_trans_atom
kg_tnadd_atomic
USE input_section_types, ONLY: section_vals_val_get
USE kg_environment_types, ONLY: kg_environment_type
USE kg_tnadd_mat, ONLY: build_tnadd_mat
USE kinds, ONLY: default_string_length,&
dp
USE kpoint_types, ONLY: get_kpoint_info,&
kpoint_type
USE lri_environment_types, ONLY: lri_environment_type
USE message_passing, ONLY: mp_para_env_type
USE particle_types, ONLY: particle_type
USE qs_condnum, ONLY: overlap_condnum
USE qs_core_matrices, ONLY: core_matrices,&
kinetic_energy_matrix
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type,&
set_qs_env
USE qs_force_types, ONLY: qs_force_type
USE qs_kind_types, ONLY: get_qs_kind,&
qs_kind_type
USE qs_kinetic, ONLY: build_kinetic_matrix
USE qs_kind_types, ONLY: qs_kind_type
USE qs_ks_types, ONLY: get_ks_env,&
qs_ks_env_type,&
set_ks_env
@ -140,7 +134,7 @@ MODULE qs_core_hamiltonian
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_core_hamiltonian'
PUBLIC :: build_core_hamiltonian_matrix, core_matrices
PUBLIC :: build_core_hamiltonian_matrix
PUBLIC :: dump_info_core_hamiltonian, qs_matrix_h_allocate_imag_from_real
CONTAINS
@ -160,9 +154,8 @@ CONTAINS
CHARACTER(LEN=*), PARAMETER :: routineN = 'build_core_hamiltonian_matrix'
INTEGER :: handle, ic, img, iw, nder, nders, &
nimages, nkind
INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
INTEGER :: handle, img, iw, nder, nders, nimages, &
nkind
LOGICAL :: h_is_complex, norml1, norml2, ofdft, &
use_arnoldi, use_virial
REAL(KIND=dp) :: eps_filter, eps_fit
@ -175,7 +168,6 @@ CONTAINS
matrix_w
TYPE(dft_control_type), POINTER :: dft_control
TYPE(kg_environment_type), POINTER :: kg_env
TYPE(kpoint_type), POINTER :: kpoints
TYPE(neighbor_list_set_p_type), DIMENSION(:), &
POINTER :: sab_orb, sap_oce
TYPE(oce_matrix_type), POINTER :: oce
@ -273,21 +265,19 @@ CONTAINS
alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
END DO
END IF
! S matrix
CALL build_overlap_matrix(ks_env, nderivative=nders, matrixkp_s=matrix_s, &
matrix_name="OVERLAP MATRIX", &
basis_type_a="ORB", &
basis_type_b="ORB", &
sab_nl=sab_orb, calculate_forces=.TRUE., &
matrixkp_p=matrix_w)
! T matrix
IF (.NOT. ofdft) &
CALL build_kinetic_matrix(ks_env, matrixkp_t=matrix_t, &
matrix_name="KINETIC ENERGY MATRIX", &
basis_type="ORB", &
sab_nl=sab_orb, calculate_forces=.TRUE., &
matrixkp_p=matrix_p, &
eps_filter=eps_filter)
ELSE
NULLIFY (matrix_p, matrix_w)
END IF
! S matrix
CALL build_overlap_matrix(ks_env, nderivative=nders, matrixkp_s=matrix_s, &
matrix_name="OVERLAP MATRIX", &
basis_type_a="ORB", &
basis_type_b="ORB", &
sab_nl=sab_orb, calculate_forces=calculate_forces, &
matrixkp_p=matrix_w)
IF (calculate_forces) THEN
! *** If LSD, then recover alpha density and beta density ***
! *** from the total density (1) and the spin density (2) ***
! *** The W matrix is neglected, since it will be destroyed ***
@ -300,22 +290,17 @@ CONTAINS
alpha_scalar=-1.0_dp, beta_scalar=1.0_dp)
END DO
END IF
ELSE
! S matrix
CALL build_overlap_matrix(ks_env, nderivative=nders, matrixkp_s=matrix_s, &
matrix_name="OVERLAP MATRIX", &
basis_type_a="ORB", &
basis_type_b="ORB", &
sab_nl=sab_orb)
! T matrix
IF (.NOT. ofdft) &
CALL build_kinetic_matrix(ks_env, matrixkp_t=matrix_t, &
matrix_name="KINETIC ENERGY MATRIX", &
basis_type="ORB", &
sab_nl=sab_orb, &
eps_filter=eps_filter)
END IF
! T matrix
CALL kinetic_energy_matrix(qs_env, matrixkp_t=matrix_t, &
matrix_p=matrix_p, &
matrix_name="KINETIC ENERGY MATRIX", &
basis_type="ORB", &
sab_orb=sab_orb, &
calculate_forces=calculate_forces, &
eps_filter=eps_filter)
! (Re-)allocate H matrix based on overlap matrix
CALL get_ks_env(ks_env, complex_ks=h_is_complex)
CALL qs_matrix_h_allocate(qs_env, matrix_s(1, 1)%matrix, is_complex=h_is_complex)
@ -334,24 +319,6 @@ CONTAINS
CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, atomic_kind_set=atomic_kind_set, &
particle_set=particle_set)
IF (.NOT. ofdft) THEN
! relativistic atomic correction to kinetic energy
IF (qs_env%rel_control%rel_method /= rel_none) THEN
IF (qs_env%rel_control%rel_transformation == rel_trans_atom) THEN
IF (nimages == 1) THEN
ic = 1
ELSE
CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
ic = cell_to_index(0, 0, 0)
END IF
CALL build_atomic_relmat(matrix_h(1, 1)%matrix, atomic_kind_set, qs_kind_set)
ELSE
CPABORT("Relativistic corrections of this type are currently not implemented")
END IF
END IF
END IF
! *** core and pseudopotentials
CALL core_matrices(qs_env, matrix_h, matrix_p, calculate_forces, nder)
@ -405,161 +372,6 @@ CONTAINS
END SUBROUTINE build_core_hamiltonian_matrix
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param matrix_h ...
!> \param matrix_p ...
!> \param calculate_forces ...
!> \param nder ...
!> \param atcore ...
! **************************************************************************************************
SUBROUTINE core_matrices(qs_env, matrix_h, matrix_p, calculate_forces, nder, atcore)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p
LOGICAL, INTENT(IN) :: calculate_forces
INTEGER, INTENT(IN) :: nder
REAL(KIND=dp), DIMENSION(:), OPTIONAL :: atcore
INTEGER :: natom, nimages
INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
LOGICAL :: all_present, my_gt_nl, ppl_present, &
ppnl_present, use_virial
REAL(KIND=dp) :: eps_ppnl
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(dft_control_type), POINTER :: dft_control
TYPE(kpoint_type), POINTER :: kpoints
TYPE(lri_environment_type), POINTER :: lri_env
TYPE(neighbor_list_set_p_type), DIMENSION(:), &
POINTER :: sab_orb, sac_ae, sac_ppl, sap_ppnl
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(qs_force_type), DIMENSION(:), POINTER :: force
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(qs_ks_env_type), POINTER :: ks_env
TYPE(virial_type), POINTER :: virial
NULLIFY (dft_control)
CALL get_qs_env(qs_env=qs_env, ks_env=ks_env, dft_control=dft_control, natom=natom)
nimages = dft_control%nimages
IF (PRESENT(atcore)) THEN
CPASSERT(SIZE(atcore) >= natom)
END IF
! check whether a gauge transformed version of the non-local potential part has to be used
my_gt_nl = .FALSE.
IF (qs_env%run_rtp) THEN
CPASSERT(ASSOCIATED(dft_control%rtp_control))
IF (dft_control%rtp_control%velocity_gauge) THEN
my_gt_nl = dft_control%rtp_control%nl_gauge_transform
END IF
END IF
! prepare for k-points
NULLIFY (cell_to_index)
IF (nimages > 1) THEN
CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
dft_control%qs_control%do_ppl_method = do_ppl_analytic
END IF
! force analytic ppl calculation for GAPW methods
IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) THEN
dft_control%qs_control%do_ppl_method = do_ppl_analytic
END IF
! force
NULLIFY (force)
IF (calculate_forces) CALL get_qs_env(qs_env=qs_env, force=force)
! virial
CALL get_qs_env(qs_env=qs_env, virial=virial)
use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
NULLIFY (qs_kind_set, atomic_kind_set, particle_set)
CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, atomic_kind_set=atomic_kind_set, &
particle_set=particle_set)
NULLIFY (sab_orb, sac_ae, sac_ppl, sap_ppnl)
CALL get_qs_env(qs_env=qs_env, &
sab_orb=sab_orb, &
sac_ae=sac_ae, &
sac_ppl=sac_ppl, &
sap_ppnl=sap_ppnl)
! *** compute the nuclear attraction contribution to the core hamiltonian ***
all_present = ASSOCIATED(sac_ae)
IF (all_present) THEN
CALL 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, "ORB", atcore=atcore)
END IF
! *** compute the ppl contribution to the core hamiltonian ***
ppl_present = ASSOCIATED(sac_ppl)
IF (ppl_present) THEN
IF (dft_control%qs_control%do_ppl_method == do_ppl_analytic) THEN
IF (dft_control%qs_control%lrigpw) THEN
CALL get_qs_env(qs_env, lri_env=lri_env)
IF (lri_env%ppl_ri) THEN
IF (lri_env%exact_1c_terms) THEN
CPABORT("not implemented")
END IF
ELSE
CALL 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, "ORB", atcore=atcore)
END IF
ELSE
CALL 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, "ORB", atcore=atcore)
END IF
END IF
END IF
! *** compute the ppnl contribution to the core hamiltonian ***
eps_ppnl = dft_control%qs_control%eps_ppnl
ppnl_present = ASSOCIATED(sap_ppnl)
IF (ppnl_present) THEN
IF (.NOT. my_gt_nl) THEN
CALL 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, "ORB", atcore=atcore)
END IF
END IF
END SUBROUTINE core_matrices
! **************************************************************************************************
!> \brief Adds atomic blocks of relativistic correction for the kinetic energy
!> \param matrix_h ...
!> \param atomic_kind_set ...
!> \param qs_kind_set ...
! **************************************************************************************************
SUBROUTINE build_atomic_relmat(matrix_h, atomic_kind_set, qs_kind_set)
TYPE(dbcsr_type), POINTER :: matrix_h
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
INTEGER :: iatom, ikind, jatom
INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
REAL(KIND=dp), DIMENSION(:, :), POINTER :: hblock, reltmat
TYPE(dbcsr_iterator_type) :: iter
CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
CALL dbcsr_iterator_start(iter, matrix_h)
DO WHILE (dbcsr_iterator_blocks_left(iter))
CALL dbcsr_iterator_next_block(iter, iatom, jatom, hblock)
IF (iatom == jatom) THEN
ikind = kind_of(iatom)
CALL get_qs_kind(qs_kind_set(ikind), reltmat=reltmat)
IF (ASSOCIATED(reltmat)) hblock = hblock + reltmat
END IF
END DO
CALL dbcsr_iterator_stop(iter)
END SUBROUTINE build_atomic_relmat
! **************************************************************************************************
!> \brief Possibly prints matrices after the construction of the Core
!> Hamiltonian Matrix

577
src/qs_core_matrices.F Normal file
View file

@ -0,0 +1,577 @@
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
! !
! SPDX-License-Identifier: GPL-2.0-or-later !
!--------------------------------------------------------------------------------------------------!
! **************************************************************************************************
!> \brief Calculation of the core Hamiltonian integral matrix <a|H|b> over
!> Cartesian Gaussian-type functions.
!>
!> Nuclear potential energy:
!>
!> a) Allelectron calculation:
!>
!> erfc(r)
!> <a|V|b> = -Z*<a|---------|b>
!> r
!>
!> 1 - erf(r)
!> = -Z*<a|------------|b>
!> r
!>
!> 1 erf(r)
!> = -Z*(<a|---|b> - <a|--------|b>)
!> r r
!>
!> 1
!> = -Z*(<a|---|b> - N*<ab||c>)
!> r
!>
!> -Z
!> = <a|---|b> + Z*N*<ab||c>
!> r
!> \_______/ \_____/
!> | |
!> nuclear coulomb
!>
!> b) Pseudopotential calculation (Goedecker, Teter and Hutter; GTH):
!>
!> <a|V|b> = <a|(V(local) + V(non-local))|b>
!>
!> = <a|(V(local)|b> + <a|V(non-local))|b>
!>
!> <a|V(local)|b> = <a|-Z(eff)*erf(SQRT(2)*alpha*r)/r +
!> (C1 + C2*(alpha*r)**2 + C3*(alpha*r)**4 +
!> C4*(alpha*r)**6)*exp(-(alpha*r)**2/2))|b>
!>
!> <a|V(non-local)|b> = <a|p(l,i)>*h(i,j)*<p(l,j)|b>
!> \par Literature
!> S. Goedecker, M. Teter and J. Hutter, Phys. Rev. B 54, 1703 (1996)
!> C. Hartwigsen, S. Goedecker and J. Hutter, Phys. Rev. B 58, 3641 (1998)
!> M. Krack and M. Parrinello, Phys. Chem. Chem. Phys. 2, 2105 (2000)
!> S. Obara and A. Saika, J. Chem. Phys. 84, 3963 (1986)
!> \par History
!> - Joost VandeVondele (April 2003) : added LSD forces
!> - Non-redundant calculation of the non-local part of the GTH PP
!> (22.05.2003,MK)
!> - New parallelization scheme (27.06.2003,MK)
!> - OpenMP version (07.12.2003,JGH)
!> - Binary search loop for VPPNL operators (09.01.2004,JGH,MK)
!> - Refactoring of pseudopotential and nuclear attraction integrals (25.02.2009,JGH)
!> - General refactoring (01.10.2010,JGH)
!> - Refactoring related to the new kinetic energy and overlap routines (07.2014,JGH)
!> - k-point functionality (07.2015,JGH)
!> - Refactor for PP functionality (08.2025,JGH)
!> \author Matthias Krack (14.09.2000,21.03.02)
! **************************************************************************************************
MODULE qs_core_matrices
USE atomic_kind_types, ONLY: atomic_kind_type,&
get_atomic_kind_set
USE cell_types, ONLY: cell_type
USE core_ae, ONLY: build_core_ae
USE core_ppl, ONLY: build_core_ppl
USE core_ppnl, ONLY: build_core_ppnl
USE cp_control_types, ONLY: dft_control_type
USE cp_dbcsr_api, ONLY: dbcsr_add,&
dbcsr_iterator_blocks_left,&
dbcsr_iterator_next_block,&
dbcsr_iterator_start,&
dbcsr_iterator_stop,&
dbcsr_iterator_type,&
dbcsr_p_type,&
dbcsr_type
USE cp_log_handling, ONLY: cp_get_default_logger,&
cp_logger_get_default_unit_nr,&
cp_logger_type
USE ec_env_types, ONLY: energy_correction_type
USE input_constants, ONLY: do_ppl_analytic,&
rel_none,&
rel_trans_atom
USE kinds, ONLY: default_string_length,&
dp
USE kpoint_types, ONLY: get_kpoint_info,&
kpoint_type
USE lri_environment_types, ONLY: lri_environment_type
USE mathlib, ONLY: det_3x3
USE message_passing, ONLY: mp_para_env_type
USE particle_types, ONLY: particle_type
USE physcon, ONLY: pascal
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
USE qs_force_types, ONLY: qs_force_type
USE qs_kind_types, ONLY: get_qs_kind,&
qs_kind_type
USE qs_kinetic, ONLY: build_kinetic_matrix
USE qs_ks_types, ONLY: get_ks_env,&
qs_ks_env_type
USE qs_linres_types, ONLY: dcdr_env_type
USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
USE virial_methods, ONLY: one_third_sum_diag
USE virial_types, ONLY: virial_type
#include "./base/base_uses.f90"
IMPLICIT NONE
PRIVATE
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_core_matrices'
PUBLIC :: core_matrices, kinetic_energy_matrix
CONTAINS
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param matrix_h ...
!> \param matrix_p ...
!> \param calculate_forces ...
!> \param nder ...
!> \param ec_env ...
!> \param dcdr_env ...
!> \param ec_env_matrices ...
!> \param basis_type ...
!> \param debug_forces ...
!> \param debug_stress ...
!> \param atcore ...
! **************************************************************************************************
SUBROUTINE core_matrices(qs_env, matrix_h, matrix_p, calculate_forces, nder, ec_env, dcdr_env, &
ec_env_matrices, basis_type, debug_forces, debug_stress, atcore)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p
LOGICAL, INTENT(IN) :: calculate_forces
INTEGER, INTENT(IN) :: nder
TYPE(energy_correction_type), OPTIONAL, POINTER :: ec_env
TYPE(dcdr_env_type), OPTIONAL :: dcdr_env
LOGICAL, INTENT(IN), OPTIONAL :: ec_env_matrices
CHARACTER(LEN=*), OPTIONAL :: basis_type
LOGICAL, INTENT(IN), OPTIONAL :: debug_forces, debug_stress
REAL(KIND=dp), DIMENSION(:), OPTIONAL :: atcore
CHARACTER(LEN=default_string_length) :: my_basis_type
INTEGER :: iounit, natom, nimages
INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
LOGICAL :: all_present, debfor, debstr, my_gt_nl, &
ppl_present, ppnl_present, use_virial
REAL(KIND=dp) :: eps_ppnl, fconv
REAL(KIND=dp), DIMENSION(3) :: fodeb
REAL(KIND=dp), DIMENSION(3, 3) :: stdeb
REAL(KIND=dp), DIMENSION(:, :), POINTER :: deltaR
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(cell_type), POINTER :: cell
TYPE(cp_logger_type), POINTER :: logger
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_hloc, matrix_ploc
TYPE(dft_control_type), POINTER :: dft_control
TYPE(kpoint_type), POINTER :: kpoints
TYPE(lri_environment_type), POINTER :: lri_env
TYPE(mp_para_env_type), POINTER :: para_env
TYPE(neighbor_list_set_p_type), DIMENSION(:), &
POINTER :: sab_orb, sac_ae, sac_ppl, sap_ppnl
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(qs_force_type), DIMENSION(:), POINTER :: force
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(qs_ks_env_type), POINTER :: ks_env
TYPE(virial_type), POINTER :: virial
logger => cp_get_default_logger()
IF (logger%para_env%is_source()) THEN
iounit = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
ELSE
iounit = -1
END IF
NULLIFY (dft_control)
CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, para_env=para_env)
! basis type
IF (PRESENT(basis_type)) THEN
my_basis_type = basis_type
ELSE
my_basis_type = "ORB"
END IF
IF (PRESENT(dcdr_env) .AND. PRESENT(ec_env)) THEN
CPABORT("core_matrices: conflicting options")
END IF
! check size of atcore array
IF (PRESENT(atcore)) THEN
CALL get_qs_env(qs_env=qs_env, natom=natom)
CPASSERT(SIZE(atcore) >= natom)
END IF
! check whether a gauge transformed version of the non-local potential part has to be used
my_gt_nl = .FALSE.
IF (qs_env%run_rtp) THEN
CPASSERT(ASSOCIATED(dft_control%rtp_control))
IF (dft_control%rtp_control%velocity_gauge) THEN
my_gt_nl = dft_control%rtp_control%nl_gauge_transform
END IF
END IF
! prepare for k-points
nimages = dft_control%nimages
NULLIFY (cell_to_index)
IF (nimages > 1) THEN
CALL get_qs_env(qs_env=qs_env, ks_env=ks_env)
CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
dft_control%qs_control%do_ppl_method = do_ppl_analytic
END IF
! Possible dc/dr terms
IF (PRESENT(dcdr_env)) THEN
deltaR => dcdr_env%delta_basis_function
matrix_hloc(1:3, 1:1) => dcdr_env%matrix_hc(1:3)
ELSE
NULLIFY (deltaR)
matrix_hloc => matrix_h
END IF
matrix_ploc => matrix_p
! force analytic ppl calculation for GAPW methods
IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) THEN
dft_control%qs_control%do_ppl_method = do_ppl_analytic
END IF
! force
NULLIFY (force)
IF (calculate_forces) CALL get_qs_env(qs_env=qs_env, force=force)
! virial
CALL get_qs_env(qs_env=qs_env, virial=virial)
use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
!debug
debfor = .FALSE.
IF (PRESENT(debug_forces)) debfor = debug_forces
debfor = debfor .AND. calculate_forces
debstr = .FALSE.
IF (PRESENT(debug_stress)) debstr = debug_stress
debstr = debstr .AND. use_virial
IF (debstr) THEN
CALL get_qs_env(qs_env=qs_env, cell=cell)
fconv = 1.0E-9_dp*pascal/cell%deth
END IF
NULLIFY (qs_kind_set, atomic_kind_set, particle_set)
CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, atomic_kind_set=atomic_kind_set, &
particle_set=particle_set)
NULLIFY (sab_orb, sac_ae, sac_ppl, sap_ppnl)
IF (PRESENT(ec_env)) THEN
sab_orb => ec_env%sab_orb
sac_ae => ec_env%sac_ae
sac_ppl => ec_env%sac_ppl
sap_ppnl => ec_env%sap_ppnl
ELSE
CALL get_qs_env(qs_env=qs_env, &
sab_orb=sab_orb, &
sac_ae=sac_ae, &
sac_ppl=sac_ppl, &
sap_ppnl=sap_ppnl)
END IF
IF (PRESENT(ec_env) .AND. PRESENT(ec_env_matrices)) THEN
IF (ec_env_matrices) THEN
matrix_hloc => ec_env%matrix_h
matrix_ploc => ec_env%matrix_p
END IF
END IF
! *** compute the nuclear attraction contribution to the core hamiltonian ***
all_present = ASSOCIATED(sac_ae)
IF (all_present) THEN
IF (PRESENT(dcdr_env)) THEN
CPABORT("ECP/AE functionality for dcdr missing")
END IF
IF (debfor) fodeb(1:3) = force(1)%all_potential(1:3, 1)
IF (debstr) stdeb = virial%pv_ppl
CALL build_core_ae(matrix_hloc, matrix_ploc, force, virial, calculate_forces, use_virial, nder, &
qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ae, &
nimages, cell_to_index, my_basis_type, atcore=atcore)
IF (debfor) THEN
fodeb(1:3) = force(1)%all_potential(1:3, 1) - fodeb(1:3)
CALL para_env%sum(fodeb)
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dHae ", fodeb
END IF
IF (debstr) THEN
stdeb = fconv*(virial%pv_ppl - stdeb)
CALL para_env%sum(stdeb)
IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
'STRESS| P*dHae ', one_third_sum_diag(stdeb), det_3x3(stdeb)
END IF
END IF
! *** compute the ppl contribution to the core hamiltonian ***
ppl_present = ASSOCIATED(sac_ppl)
IF (ppl_present) THEN
IF (dft_control%qs_control%do_ppl_method == do_ppl_analytic) THEN
CALL get_qs_env(qs_env, lri_env=lri_env)
IF (dft_control%qs_control%lrigpw .AND. lri_env%ppl_ri) THEN
IF (lri_env%exact_1c_terms) THEN
CPABORT("not implemented")
END IF
ELSE
IF (debfor) fodeb(1:3) = force(1)%gth_ppl(1:3, 1)
IF (debstr) stdeb = virial%pv_ppl
CALL build_core_ppl(matrix_hloc, matrix_ploc, force, &
virial, calculate_forces, use_virial, nder, &
qs_kind_set, atomic_kind_set, particle_set, &
sab_orb, sac_ppl, nimages, cell_to_index, my_basis_type, &
deltaR=deltaR, atcore=atcore)
IF (debfor) THEN
fodeb(1:3) = force(1)%gth_ppl(1:3, 1) - fodeb(1:3)
CALL para_env%sum(fodeb)
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dHppl ", fodeb
END IF
IF (debstr) THEN
stdeb = fconv*(virial%pv_ppl - stdeb)
CALL para_env%sum(stdeb)
IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
'STRESS| P*dHppl ', one_third_sum_diag(stdeb), det_3x3(stdeb)
END IF
END IF
END IF
END IF
! *** compute the ppnl contribution to the core hamiltonian ***
eps_ppnl = dft_control%qs_control%eps_ppnl
ppnl_present = ASSOCIATED(sap_ppnl)
IF (ppnl_present) THEN
IF (PRESENT(dcdr_env)) THEN
matrix_hloc(1:3, 1:1) => dcdr_env%matrix_ppnl_1(1:3)
END IF
IF (.NOT. my_gt_nl) THEN
IF (debfor) fodeb(1:3) = force(1)%gth_ppnl(1:3, 1)
IF (debstr) stdeb = virial%pv_ppnl
CALL build_core_ppnl(matrix_hloc, matrix_ploc, 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, my_basis_type, &
deltaR=deltaR, atcore=atcore)
IF (debfor) THEN
fodeb(1:3) = force(1)%gth_ppnl(1:3, 1) - fodeb(1:3)
CALL para_env%sum(fodeb)
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dHppnl ", fodeb
END IF
IF (debstr) THEN
stdeb = fconv*(virial%pv_ppnl - stdeb)
CALL para_env%sum(stdeb)
IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
'STRESS| P*dHppnl ', one_third_sum_diag(stdeb), det_3x3(stdeb)
END IF
END IF
END IF
END SUBROUTINE core_matrices
! **************************************************************************************************
!> \brief Calculate kinetic energy matrix and possible relativistic correction
!> \param qs_env ...
!> \param matrixkp_t ...
!> \param matrix_t ...
!> \param matrix_p ...
!> \param matrix_name ...
!> \param calculate_forces ...
!> \param nderivative ...
!> \param sab_orb ...
!> \param eps_filter ...
!> \param basis_type ...
!> \param debug_forces ...
!> \param debug_stress ...
!> \author Creation (21.08.2025,JGH)
! **************************************************************************************************
SUBROUTINE kinetic_energy_matrix(qs_env, matrixkp_t, matrix_t, matrix_p, &
matrix_name, calculate_forces, nderivative, &
sab_orb, eps_filter, basis_type, debug_forces, debug_stress)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
POINTER :: matrixkp_t
TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
POINTER :: matrix_t
TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
POINTER :: matrix_p
CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: matrix_name
LOGICAL, INTENT(IN) :: calculate_forces
INTEGER, INTENT(IN), OPTIONAL :: nderivative
TYPE(neighbor_list_set_p_type), DIMENSION(:), &
OPTIONAL, POINTER :: sab_orb
REAL(KIND=dp), INTENT(IN), OPTIONAL :: eps_filter
CHARACTER(LEN=*), OPTIONAL :: basis_type
LOGICAL, INTENT(IN), OPTIONAL :: debug_forces, debug_stress
CHARACTER(LEN=default_string_length) :: my_basis_type
INTEGER :: ic, img, iounit, nimages
INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
LOGICAL :: debfor, debstr, use_virial
REAL(KIND=dp) :: fconv
REAL(KIND=dp), DIMENSION(3) :: fodeb
REAL(KIND=dp), DIMENSION(3, 3) :: stdeb
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(cell_type), POINTER :: cell
TYPE(cp_logger_type), POINTER :: logger
TYPE(dft_control_type), POINTER :: dft_control
TYPE(kpoint_type), POINTER :: kpoints
TYPE(mp_para_env_type), POINTER :: para_env
TYPE(neighbor_list_set_p_type), DIMENSION(:), &
POINTER :: sab_nl
TYPE(qs_force_type), DIMENSION(:), POINTER :: force
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(qs_ks_env_type), POINTER :: ks_env
TYPE(virial_type), POINTER :: virial
logger => cp_get_default_logger()
IF (logger%para_env%is_source()) THEN
iounit = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
ELSE
iounit = -1
END IF
CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, para_env=para_env)
! is this a orbital-free method calculation
IF (dft_control%qs_control%ofgpw) RETURN
! Matrix images (kp)
nimages = dft_control%nimages
! basis type
IF (PRESENT(basis_type)) THEN
my_basis_type = basis_type
ELSE
my_basis_type = "ORB"
END IF
debfor = .FALSE.
IF (PRESENT(debug_forces)) debfor = debug_forces
debfor = debfor .AND. calculate_forces
CALL get_qs_env(qs_env=qs_env, virial=virial)
use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
debstr = .FALSE.
IF (PRESENT(debug_stress)) debstr = debug_stress
debstr = debstr .AND. use_virial
IF (debstr) THEN
CALL get_qs_env(qs_env=qs_env, cell=cell)
fconv = 1.0E-9_dp*pascal/cell%deth
END IF
NULLIFY (ks_env)
CALL get_qs_env(qs_env=qs_env, ks_env=ks_env)
NULLIFY (sab_nl)
IF (PRESENT(sab_orb)) THEN
sab_nl => sab_orb
ELSE
CALL get_qs_env(qs_env=qs_env, sab_orb=sab_nl)
END IF
IF (calculate_forces) THEN
IF (SIZE(matrix_p, 1) == 2) THEN
DO img = 1, nimages
CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
END DO
END IF
END IF
IF (debfor) THEN
CALL get_qs_env(qs_env=qs_env, force=force)
fodeb(1:3) = force(1)%kinetic(1:3, 1)
END IF
IF (debstr) THEN
stdeb = virial%pv_ekinetic
END IF
! T matrix
CALL build_kinetic_matrix(ks_env, matrixkp_t=matrixkp_t, matrix_t=matrix_t, &
matrix_name=matrix_name, &
basis_type=my_basis_type, &
sab_nl=sab_nl, &
calculate_forces=calculate_forces, &
matrixkp_p=matrix_p, &
nderivative=nderivative, &
eps_filter=eps_filter)
IF (debfor) THEN
fodeb(1:3) = force(1)%kinetic(1:3, 1) - fodeb(1:3)
CALL para_env%sum(fodeb)
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dT ", fodeb
END IF
IF (debstr) THEN
stdeb = fconv*(virial%pv_ekinetic - stdeb)
CALL para_env%sum(stdeb)
IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
'STRESS| P*dT', one_third_sum_diag(stdeb), det_3x3(stdeb)
END IF
IF (calculate_forces) THEN
IF (SIZE(matrix_p, 1) == 2) THEN
DO img = 1, nimages
CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
END DO
END IF
END IF
! relativistic atomic correction to kinetic energy
IF (qs_env%rel_control%rel_method /= rel_none) THEN
IF (qs_env%rel_control%rel_transformation == rel_trans_atom) THEN
CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, atomic_kind_set=atomic_kind_set)
IF (nimages == 1) THEN
ic = 1
ELSE
CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
ic = cell_to_index(0, 0, 0)
END IF
IF (my_basis_type /= "ORB") THEN
CPABORT("Basis mismatch for relativistic correction")
END IF
IF (PRESENT(matrixkp_t)) THEN
CALL build_atomic_relmat(matrixkp_t(1, ic)%matrix, atomic_kind_set, qs_kind_set)
ELSEIF (PRESENT(matrix_t)) THEN
CALL build_atomic_relmat(matrix_t(1)%matrix, atomic_kind_set, qs_kind_set)
END IF
ELSE
CPABORT("Relativistic corrections of this type are currently not implemented")
END IF
END IF
END SUBROUTINE kinetic_energy_matrix
! **************************************************************************************************
!> \brief Adds atomic blocks of relativistic correction for the kinetic energy
!> \param matrix_t ...
!> \param atomic_kind_set ...
!> \param qs_kind_set ...
! **************************************************************************************************
SUBROUTINE build_atomic_relmat(matrix_t, atomic_kind_set, qs_kind_set)
TYPE(dbcsr_type), POINTER :: matrix_t
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
INTEGER :: iatom, ikind, jatom
INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
REAL(KIND=dp), DIMENSION(:, :), POINTER :: reltmat, tblock
TYPE(dbcsr_iterator_type) :: iter
CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
CALL dbcsr_iterator_start(iter, matrix_t)
DO WHILE (dbcsr_iterator_blocks_left(iter))
CALL dbcsr_iterator_next_block(iter, iatom, jatom, tblock)
IF (iatom == jatom) THEN
ikind = kind_of(iatom)
CALL get_qs_kind(qs_kind_set(ikind), reltmat=reltmat)
IF (ASSOCIATED(reltmat)) tblock = tblock + reltmat
END IF
END DO
CALL dbcsr_iterator_stop(iter)
END SUBROUTINE build_atomic_relmat
END MODULE qs_core_matrices

90
src/qs_dcdr_ao.F
View file

@ -12,12 +12,8 @@
MODULE qs_dcdr_ao
USE atomic_kind_types, ONLY: atomic_kind_type
USE basis_set_types, ONLY: gto_basis_set_p_type,&
gto_basis_set_type
USE core_ae, ONLY: build_core_ae
USE core_ppl, ONLY: build_core_ppl
USE core_ppnl, ONLY: build_core_ppnl
USE cp_control_types, ONLY: dft_control_type
USE cp_dbcsr_api, ONLY: dbcsr_copy,&
dbcsr_get_block_p,&
@ -31,14 +27,12 @@ MODULE qs_dcdr_ao
cp_fm_type
USE cp_log_handling, ONLY: cp_get_default_logger,&
cp_logger_type
USE input_constants, ONLY: do_ppl_analytic
USE input_section_types, ONLY: section_vals_get_subs_vals,&
section_vals_type
USE kinds, ONLY: default_string_length,&
dp
USE orbital_pointers, ONLY: ncoset
USE parallel_gemm_api, ONLY: parallel_gemm
USE particle_types, ONLY: particle_type
USE pw_env_types, ONLY: pw_env_get,&
pw_env_type
USE pw_methods, ONLY: pw_axpy,&
@ -54,10 +48,10 @@ MODULE qs_dcdr_ao
pw_r3d_rs_type
USE qs_collocate_density, ONLY: calculate_drho_core,&
calculate_drho_elec_dR
USE qs_core_matrices, ONLY: core_matrices
USE qs_energy_types, ONLY: qs_energy_type
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
USE qs_force_types, ONLY: qs_force_type
USE qs_integral_utils, ONLY: basis_set_list_setup,&
get_memory_usage
USE qs_integrate_potential, ONLY: integrate_v_dbasis,&
@ -80,7 +74,6 @@ MODULE qs_dcdr_ao
qs_rho_release,&
qs_rho_type
USE qs_vxc, ONLY: qs_vxc_create
USE virial_types, ONLY: virial_type
USE xc, ONLY: xc_calc_2nd_deriv,&
xc_prep_2nd_deriv
USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
@ -403,95 +396,30 @@ CONTAINS
CHARACTER(LEN=default_string_length) :: my_basis_type
INTEGER :: handle, nder
INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
LOGICAL :: calculate_forces, failure, ppl_present, &
ppnl_present, use_virial
REAL(KIND=dp) :: eps_ppnl
REAL(KIND=dp), DIMENSION(:, :), POINTER :: deltaR
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
LOGICAL :: calculate_forces
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_hc_pass, matrix_p_pass, &
matrix_ppnl_1_pass
TYPE(dft_control_type), POINTER :: dft_control
TYPE(neighbor_list_set_p_type), DIMENSION(:), &
POINTER :: sab_orb, sac_ae, sac_ppl, sap_ppnl
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(qs_force_type), DIMENSION(:), POINTER :: force
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p_pass
TYPE(qs_ks_env_type), POINTER :: ks_env
TYPE(qs_rho_type), POINTER :: rho
TYPE(virial_type), POINTER :: virial
CALL timeset(routineN, handle)
failure = .FALSE.
NULLIFY (atomic_kind_set, qs_kind_set, ks_env, dft_control, particle_set, &
sab_orb, sac_ae, sac_ppl, sap_ppnl, virial, rho, rho_ao)
CALL get_qs_env(qs_env=qs_env, &
atomic_kind_set=atomic_kind_set, &
qs_kind_set=qs_kind_set, &
ks_env=ks_env, &
dft_control=dft_control, &
particle_set=particle_set, &
sab_orb=sab_orb, &
sac_ae=sac_ae, &
sac_ppl=sac_ppl, &
sap_ppnl=sap_ppnl, &
virial=virial)
CALL get_qs_env(qs_env=qs_env, ks_env=ks_env)
CALL get_ks_env(ks_env=ks_env, rho=rho)
CALL qs_rho_get(rho, rho_ao=rho_ao)
deltaR => dcdr_env%delta_basis_function
nder = 1
calculate_forces = .FALSE.
my_basis_type = "ORB"
! ECP/AE contribution to the core hamiltonian
IF (ASSOCIATED(sac_ae)) THEN
CPABORT("ECP/AE functionality in qs_dcdr_ao missing")
! Missing feature: deltaR weighting factors of the derivatives wrt. nuclear positions
matrix_hc_pass(1:3, 1:1) => dcdr_env%matrix_hc(1:3)
matrix_p_pass(1:1, 1:1) => rho_ao(1:1)
CALL build_core_ae(matrix_h=matrix_hc_pass, matrix_p=matrix_p_pass, &
force=force, virial=virial, calculate_forces=calculate_forces, &
use_virial=use_virial, nder=nder, qs_kind_set=qs_kind_set, &
atomic_kind_set=atomic_kind_set, particle_set=particle_set, &
sab_orb=sab_orb, sac_ae=sac_ae, nimages=1, &
basis_type=my_basis_type, cell_to_index=cell_to_index)
END IF
! *** compute the ppl contribution to the core hamiltonian ***
ppl_present = ASSOCIATED(sac_ppl)
IF (ppl_present) THEN
IF (dft_control%qs_control%do_ppl_method == do_ppl_analytic) THEN
matrix_hc_pass(1:3, 1:1) => dcdr_env%matrix_hc(1:3)
matrix_p_pass(1:1, 1:1) => rho_ao(1:1)
CALL build_core_ppl(matrix_h=matrix_hc_pass, matrix_p=matrix_p_pass, &
force=force, virial=virial, calculate_forces=calculate_forces, &
use_virial=use_virial, nder=nder, qs_kind_set=qs_kind_set, &
atomic_kind_set=atomic_kind_set, particle_set=particle_set, &
sab_orb=sab_orb, sac_ppl=sac_ppl, basis_type=my_basis_type, &
nimages=1, cell_to_index=cell_to_index, deltaR=deltaR)
END IF ! ppl_analytic
END IF ! ppl_present
! *** compute the ppnl contribution to the core hamiltonian ***
eps_ppnl = dft_control%qs_control%eps_ppnl
ppnl_present = ASSOCIATED(sap_ppnl)
IF (ppnl_present) THEN
matrix_ppnl_1_pass(1:3, 1:1) => dcdr_env%matrix_ppnl_1(1:3)
CALL build_core_ppnl(matrix_h=matrix_ppnl_1_pass, matrix_p=matrix_p_pass, force=force, virial=virial, &
calculate_forces=calculate_forces, use_virial=use_virial, nder=nder, &
qs_kind_set=qs_kind_set, atomic_kind_set=atomic_kind_set, &
particle_set=particle_set, sab_orb=sab_orb, sap_ppnl=sap_ppnl, &
eps_ppnl=eps_ppnl, nimages=1, cell_to_index=cell_to_index, &
basis_type=my_basis_type, deltaR=deltaR)
END IF
NULLIFY (matrix_h)
matrix_p_pass(1:1, 1:1) => rho_ao(1:1)
CALL core_matrices(qs_env, matrix_h, matrix_p_pass, calculate_forces, nder, &
dcdr_env=dcdr_env)
CALL timestop(handle)
END SUBROUTINE core_dR
! **************************************************************************************************

22
src/qs_dcdr_utils.F
View file

@ -64,9 +64,9 @@ MODULE qs_dcdr_utils
USE moments_utils, ONLY: get_reference_point
USE parallel_gemm_api, ONLY: parallel_gemm
USE particle_types, ONLY: particle_type
USE qs_core_matrices, ONLY: kinetic_energy_matrix
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
USE qs_kinetic, ONLY: build_kinetic_matrix
USE qs_ks_types, ONLY: qs_ks_env_type
USE qs_linres_types, ONLY: dcdr_env_type,&
linres_control_type
@ -666,6 +666,7 @@ CONTAINS
TYPE(cp_fm_type), POINTER :: mo_coeff
TYPE(cp_logger_type), POINTER :: logger
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p
TYPE(dft_control_type), POINTER :: dft_control
TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
@ -723,11 +724,20 @@ CONTAINS
sab_nl=sab_orb)
NULLIFY (dcdr_env%matrix_t)
CALL build_kinetic_matrix(ks_env, matrix_t=dcdr_env%matrix_t, &
matrix_name="KINETIC ENERGY MATRIX", &
basis_type="ORB", &
sab_nl=sab_orb, nderivative=1, &
eps_filter=dft_control%qs_control%eps_filter_matrix)
NULLIFY (matrix_p)
CALL kinetic_energy_matrix(qs_env, matrix_t=dcdr_env%matrix_t, matrix_p=matrix_p, &
matrix_name="KINETIC ENERGY MATRIX", &
calculate_forces=.FALSE., &
basis_type="ORB", &
sab_orb=sab_orb, &
nderivative=1, &
eps_filter=dft_control%qs_control%eps_filter_matrix)
! CALL build_kinetic_matrix(ks_env, matrix_t=dcdr_env%matrix_t, &
! matrix_name="KINETIC ENERGY MATRIX", &
! basis_type="ORB", &
! sab_nl=sab_orb, nderivative=1, &
! eps_filter=dft_control%qs_control%eps_filter_matrix)
! Get inputs
CALL section_vals_val_get(dcdr_section, "DISTRIBUTED_ORIGIN", l_val=dcdr_env%distributed_origin)

2
src/qs_energy_utils.F
View file

@ -40,7 +40,7 @@ MODULE qs_energy_utils
pw_scale
USE pw_pool_types, ONLY: pw_pool_type
USE pw_types, ONLY: pw_r3d_rs_type
USE qs_core_hamiltonian, ONLY: core_matrices
USE qs_core_matrices, ONLY: core_matrices
USE qs_dispersion_types, ONLY: qs_dispersion_type
USE qs_energy_types, ONLY: qs_energy_type
USE qs_environment_types, ONLY: get_qs_env,&

15
src/qs_kinetic.F
View file

@ -186,6 +186,7 @@ CONTAINS
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: dab
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: k_block
TYPE(dbcsr_type), POINTER :: matp
TYPE(dft_control_type), POINTER :: dft_control
TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
@ -260,7 +261,13 @@ CONTAINS
IF (dokp) THEN
CPASSERT(PRESENT(matrixkp_p))
ELSE
CPASSERT(PRESENT(matrix_p))
IF (PRESENT(matrix_p)) THEN
matp => matrix_p
ELSEIF (PRESENT(matrixkp_p)) THEN
matp => matrixkp_p(1, 1)%matrix
ELSE
CPABORT("Missing density matrix")
END IF
END IF
END IF
@ -273,7 +280,7 @@ CONTAINS
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP SHARED (do_forces, ldsab, use_cell_mapping, do_symmetric, dokp,&
!$OMP sab_nl, ncoset, maxder, nder, ndod, use_virial, matrix_t, matrixkp_t,&
!$OMP matrix_p, basis_set_list, atom_of_kind, cell_to_index, matrixkp_p, locks, natom) &
!$OMP matp, matrix_p, basis_set_list, atom_of_kind, cell_to_index, matrixkp_p, locks, natom) &
!$OMP PRIVATE (k_block, kab, qab, pab, ikind, jkind, iatom, jatom, rab, cell, basis_set_a, basis_set_b, f, &
!$OMP first_sgfa, la_max, la_min, npgfa, nsgfa, nseta, rpgfa, set_radius_a, ncoa, ncob, force_a, &
!$OMP zeta, first_sgfb, lb_max, lb_min, npgfb, nsetb, rpgfb, set_radius_b, nsgfb, p_block, dab, tab, &
@ -386,7 +393,9 @@ CONTAINS
row=irow, col=icol, block=p_block, found=found)
CPASSERT(found)
ELSE
CALL dbcsr_get_block_p(matrix=matrix_p, row=irow, col=icol, &
!deb CALL dbcsr_get_block_p(matrix=matrix_p, row=irow, col=icol, &
!deb block=p_block, found=found)
CALL dbcsr_get_block_p(matrix=matp, row=irow, col=icol, &
block=p_block, found=found)
CPASSERT(found)
END IF

126
src/response_solver.F
View file

@ -21,9 +21,6 @@ MODULE response_solver
USE atomic_kind_types, ONLY: atomic_kind_type,&
get_atomic_kind
USE cell_types, ONLY: cell_type
USE core_ae, ONLY: build_core_ae
USE core_ppl, ONLY: build_core_ppl
USE core_ppnl, ONLY: build_core_ppnl
USE cp_blacs_env, ONLY: cp_blacs_env_type
USE cp_control_types, ONLY: dft_control_type
USE cp_dbcsr_api, ONLY: &
@ -98,6 +95,8 @@ MODULE response_solver
pw_r3d_rs_type
USE qs_2nd_kernel_ao, ONLY: build_dm_response
USE qs_collocate_density, ONLY: calculate_rho_elec
USE qs_core_matrices, ONLY: core_matrices,&
kinetic_energy_matrix
USE qs_density_matrices, ONLY: calculate_whz_matrix,&
calculate_wz_matrix
USE qs_energy_types, ONLY: qs_energy_type
@ -112,7 +111,6 @@ MODULE response_solver
USE qs_kind_types, ONLY: get_qs_kind,&
get_qs_kind_set,&
qs_kind_type
USE qs_kinetic, ONLY: build_kinetic_matrix
USE qs_ks_atom, ONLY: update_ks_atom
USE qs_ks_methods, ONLY: calc_rho_tot_gspace
USE qs_ks_types, ONLY: qs_ks_env_type
@ -834,13 +832,12 @@ CONTAINS
CHARACTER(LEN=default_string_length) :: basis_type, unitstr
INTEGER :: handle, iounit, ispin, mspin, myfun, &
n_rep_hf, nao, nao_aux, natom, nder, &
nimages, nocc, nspins
INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
nocc, nspins
LOGICAL :: debug_forces, debug_stress, distribute_fock_matrix, do_ex, do_hfx, gapw, gapw_xc, &
hfx_treat_lsd_in_core, resp_only, s_mstruct_changed, use_virial
REAL(KIND=dp) :: eh1, ehartree, ekin_mol, eps_filter, &
eps_ppnl, exc, exc_aux_fit, fconv, &
focc, hartree_gs, hartree_t
exc, exc_aux_fit, fconv, focc, &
hartree_gs, hartree_t
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: ftot1, ftot2, ftot3
REAL(KIND=dp), DIMENSION(2) :: total_rho_gs, total_rho_t
REAL(KIND=dp), DIMENSION(3) :: fodeb
@ -853,7 +850,7 @@ CONTAINS
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ht, matrix_pd, matrix_pza, &
matrix_s, mpa, scrm
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p, mhd, mhx, mhy, mhz, &
mpd, mpz
mpa2, mpd, mpz, scrm2
TYPE(dbcsr_type), POINTER :: dbwork
TYPE(dft_control_type), POINTER :: dft_control
TYPE(hartree_local_type), POINTER :: hartree_local_gs, hartree_local_t
@ -865,8 +862,7 @@ CONTAINS
TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
TYPE(mp_para_env_type), POINTER :: para_env
TYPE(neighbor_list_set_p_type), DIMENSION(:), &
POINTER :: sab_aux_fit, sab_orb, sac_ae, sac_ppl, &
sap_ppnl
POINTER :: sab_aux_fit, sab_orb
TYPE(oce_matrix_type), POINTER :: oce
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(pw_c1d_gs_type) :: rho_tot_gspace, rho_tot_gspace_gs, rho_tot_gspace_t, &
@ -900,7 +896,6 @@ CONTAINS
debug_stress = .FALSE.
END IF
logger => cp_get_default_logger()
logger => cp_get_default_logger()
IF (logger%para_env%is_source()) THEN
iounit = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
@ -914,7 +909,7 @@ CONTAINS
CPASSERT(PRESENT(p_env))
END IF
NULLIFY (ks_env, sab_orb, sac_ae, sac_ppl, sap_ppnl, virial)
NULLIFY (ks_env, sab_orb, virial)
CALL get_qs_env(qs_env=qs_env, &
cell=cell, &
force=force, &
@ -922,9 +917,6 @@ CONTAINS
dft_control=dft_control, &
para_env=para_env, &
sab_orb=sab_orb, &
sac_ae=sac_ae, &
sac_ppl=sac_ppl, &
sap_ppnl=sap_ppnl, &
virial=virial)
nspins = dft_control%nspins
gapw = dft_control%qs_control%gapw
@ -975,39 +967,20 @@ CONTAINS
!
! START OF Tr(P+Z)Hcore
!
IF (nspins == 2) THEN
CALL dbcsr_add(mpa(1)%matrix, mpa(2)%matrix, 1.0_dp, 1.0_dp)
END IF
nimages = 1
! Kinetic energy matrix
NULLIFY (scrm)
IF (debug_forces) fodeb(1:3) = force(1)%kinetic(1:3, 1)
IF (debug_stress .AND. use_virial) stdeb = virial%pv_ekinetic
CALL build_kinetic_matrix(ks_env, matrix_t=scrm, &
matrix_name="KINETIC ENERGY MATRIX", &
basis_type="ORB", &
sab_nl=sab_orb, calculate_forces=.TRUE., &
matrix_p=mpa(1)%matrix)
IF (debug_forces) THEN
fodeb(1:3) = force(1)%kinetic(1:3, 1) - fodeb(1:3)
CALL para_env%sum(fodeb)
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pz*dT ", fodeb
END IF
IF (debug_stress .AND. use_virial) THEN
stdeb = fconv*(virial%pv_ekinetic - stdeb)
CALL para_env%sum(stdeb)
IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
'STRESS| Kinetic energy', one_third_sum_diag(stdeb), det_3x3(stdeb)
END IF
IF (nspins == 2) THEN
CALL dbcsr_add(mpa(1)%matrix, mpa(2)%matrix, 1.0_dp, -1.0_dp)
END IF
CALL dbcsr_deallocate_matrix_set(scrm)
NULLIFY (scrm2)
mpa2(1:nspins, 1:1) => mpa(1:nspins)
CALL kinetic_energy_matrix(qs_env, matrixkp_t=scrm2, matrix_p=mpa2, &
matrix_name="KINETIC ENERGY MATRIX", &
basis_type="ORB", &
sab_orb=sab_orb, calculate_forces=.TRUE., &
debug_forces=debug_forces, debug_stress=debug_stress)
CALL dbcsr_deallocate_matrix_set(scrm2)
! Initialize a matrix scrm, later used for scratch purposes
CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s)
NULLIFY (scrm)
CALL dbcsr_allocate_matrix_set(scrm, nspins)
DO ispin = 1, nspins
ALLOCATE (scrm(ispin)%matrix)
@ -1019,7 +992,6 @@ CONTAINS
CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, particle_set=particle_set, &
atomic_kind_set=atomic_kind_set)
NULLIFY (cell_to_index)
ALLOCATE (matrix_p(nspins, 1), matrix_h(nspins, 1))
DO ispin = 1, nspins
matrix_p(ispin, 1)%matrix => mpa(ispin)%matrix
@ -1027,68 +999,10 @@ CONTAINS
END DO
matrix_h(1, 1)%matrix => scrm(1)%matrix
IF (ASSOCIATED(sac_ae)) THEN
nder = 1
IF (debug_forces) fodeb(1:3) = force(1)%all_potential(1:3, 1)
IF (debug_stress .AND. use_virial) stdeb = virial%pv_ppl
CALL build_core_ae(matrix_h, matrix_p, force, virial, .TRUE., use_virial, nder, &
qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ae, &
nimages, cell_to_index, "ORB")
IF (debug_forces) THEN
fodeb(1:3) = force(1)%all_potential(1:3, 1) - fodeb(1:3)
CALL para_env%sum(fodeb)
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pz*dHae ", fodeb
END IF
IF (debug_stress .AND. use_virial) THEN
stdeb = fconv*(virial%pv_ppl - stdeb)
CALL para_env%sum(stdeb)
IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
'STRESS| Pz*dHae ', one_third_sum_diag(stdeb), det_3x3(stdeb)
END IF
END IF
nder = 1
CALL core_matrices(qs_env, matrix_h, matrix_p, .TRUE., nder, &
debug_forces=debug_forces, debug_stress=debug_stress)
IF (ASSOCIATED(sac_ppl)) THEN
nder = 1
IF (debug_forces) fodeb(1:3) = force(1)%gth_ppl(1:3, 1)
IF (debug_stress .AND. use_virial) stdeb = virial%pv_ppl
CALL build_core_ppl(matrix_h, matrix_p, force, &
virial, .TRUE., use_virial, nder, &
qs_kind_set, atomic_kind_set, particle_set, &
sab_orb, sac_ppl, nimages, cell_to_index, "ORB")
IF (debug_forces) THEN
fodeb(1:3) = force(1)%gth_ppl(1:3, 1) - fodeb(1:3)
CALL para_env%sum(fodeb)
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pz*dHppl ", fodeb
END IF
IF (debug_stress .AND. use_virial) THEN
stdeb = fconv*(virial%pv_ppl - stdeb)
CALL para_env%sum(stdeb)
IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
'STRESS| Pz*dHppl ', one_third_sum_diag(stdeb), det_3x3(stdeb)
END IF
END IF
eps_ppnl = dft_control%qs_control%eps_ppnl
IF (ASSOCIATED(sap_ppnl)) THEN
nder = 1
IF (debug_forces) fodeb(1:3) = force(1)%gth_ppnl(1:3, 1)
IF (debug_stress .AND. use_virial) stdeb = virial%pv_ppnl
CALL build_core_ppnl(matrix_h, matrix_p, force, &
virial, .TRUE., use_virial, nder, &
qs_kind_set, atomic_kind_set, particle_set, &
sab_orb, sap_ppnl, eps_ppnl, nimages, cell_to_index, "ORB")
IF (debug_forces) THEN
fodeb(1:3) = force(1)%gth_ppnl(1:3, 1) - fodeb(1:3)
CALL para_env%sum(fodeb)
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pz*dHppnl ", fodeb
END IF
IF (debug_stress .AND. use_virial) THEN
stdeb = fconv*(virial%pv_ppnl - stdeb)
CALL para_env%sum(stdeb)
IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
'STRESS| Pz*dHppnl ', one_third_sum_diag(stdeb), det_3x3(stdeb)
END IF
END IF
! Kim-Gordon subsystem DFT
! Atomic potential for nonadditive kinetic energy contribution
IF (dft_control%qs_control%do_kg) THEN

61
src/rpa_im_time_force_methods.F
View file

@ -22,9 +22,6 @@ MODULE rpa_im_time_force_methods
cite_reference
USE cell_types, ONLY: cell_type,&
pbc
USE core_ae, ONLY: build_core_ae
USE core_ppl, ONLY: build_core_ppl
USE core_ppnl, ONLY: build_core_ppnl
USE cp_blacs_env, ONLY: cp_blacs_env_type
USE cp_control_types, ONLY: dft_control_type
USE cp_dbcsr_api, ONLY: &
@ -118,6 +115,8 @@ MODULE rpa_im_time_force_methods
pw_r3d_rs_type
USE qs_collocate_density, ONLY: calculate_rho_elec,&
collocate_function
USE qs_core_matrices, ONLY: core_matrices,&
kinetic_energy_matrix
USE qs_density_matrices, ONLY: calculate_whz_matrix
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type,&
@ -129,7 +128,6 @@ MODULE rpa_im_time_force_methods
integrate_v_rspace
USE qs_interactions, ONLY: init_interaction_radii_orb_basis
USE qs_kind_types, ONLY: qs_kind_type
USE qs_kinetic, ONLY: build_kinetic_matrix
USE qs_ks_methods, ONLY: calc_rho_tot_gspace
USE qs_ks_reference, ONLY: ks_ref_potential
USE qs_ks_types, ONLY: set_ks_env
@ -3136,18 +3134,17 @@ CONTAINS
INTEGER :: handle, i, idens, ispin, n_rep_hf, nao, &
nao_aux, nder, nimages, nocc, nspins
INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
LOGICAL :: do_exx, do_hfx, do_tau, do_tau_admm, &
use_virial
REAL(dp) :: dummy_real1, dummy_real2, ehartree, &
eps_ppnl, exc, focc
REAL(dp) :: dummy_real1, dummy_real2, ehartree, exc, &
focc
REAL(dp), DIMENSION(3, 3) :: h_stress, pv_loc
TYPE(admm_type), POINTER :: admm_env
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: current_density, current_density_admm, &
current_mat_h, matrix_p_mp2, matrix_p_mp2_admm, matrix_s, matrix_s_aux_fit, matrix_w, &
rho_ao, rho_ao_aux, scrm, scrm_admm
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: dbcsr_work_h, dbcsr_work_p
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: dbcsr_work_h, dbcsr_work_p, mpa2
TYPE(dbcsr_type) :: dbcsr_work
TYPE(dft_control_type), POINTER :: dft_control
TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
@ -3172,12 +3169,14 @@ CONTAINS
TYPE(task_list_type), POINTER :: task_list_aux_fit
TYPE(virial_type), POINTER :: virial
NULLIFY (scrm, rho, dft_control, matrix_p_mp2, matrix_s, matrix_p_mp2_admm, admm_env, sab_orb, &
cell_to_index, dbcsr_work_p, dbcsr_work_h, sac_ae, sac_ppl, sap_ppnl, force, virial, &
qs_kind_set, atomic_kind_set, particle_set, pw_env, poisson_env, auxbas_pw_pool, &
task_list_aux_fit, matrix_s_aux_fit, scrm_admm, rho_aux_fit, rho_ao_aux, x_data, &
hfx_section, xc_section, para_env, rhoz_g, rhoz_r, tauz_r, v_xc, v_xc_tau, &
vxc_rspace, vtau_rspace, vadmm_rspace, rho_ao, matrix_w)
NULLIFY (scrm, rho, dft_control, matrix_p_mp2, matrix_s, &
matrix_p_mp2_admm, admm_env, sab_orb, dbcsr_work_p, &
dbcsr_work_h, sac_ae, sac_ppl, sap_ppnl, force, virial, &
qs_kind_set, atomic_kind_set, particle_set, pw_env, poisson_env, &
auxbas_pw_pool, task_list_aux_fit, matrix_s_aux_fit, scrm_admm, &
rho_aux_fit, rho_ao_aux, x_data, hfx_section, xc_section, &
para_env, rhoz_g, rhoz_r, tauz_r, v_xc, v_xc_tau, vxc_rspace, &
vtau_rspace, vadmm_rspace, rho_ao, matrix_w)
CALL timeset(routineN, handle)
@ -3202,13 +3201,11 @@ CONTAINS
!The kinetic term (only response density)
NULLIFY (scrm)
IF (nspins == 2) CALL dbcsr_add(matrix_p_mp2(1)%matrix, matrix_p_mp2(2)%matrix, 1.0_dp, 1.0_dp)
CALL build_kinetic_matrix(qs_env%ks_env, matrix_t=scrm, &
matrix_name="KINETIC ENERGY MATRIX", &
basis_type="ORB", &
sab_nl=sab_orb, calculate_forces=.TRUE., &
matrix_p=matrix_p_mp2(1)%matrix)
IF (nspins == 2) CALL dbcsr_add(matrix_p_mp2(1)%matrix, matrix_p_mp2(2)%matrix, 1.0_dp, -1.0_dp)
mpa2(1:nspins, 1:1) => matrix_p_mp2(1:nspins)
CALL kinetic_energy_matrix(qs_env, matrix_t=scrm, matrix_p=mpa2, &
matrix_name="KINETIC ENERGY MATRIX", &
basis_type="ORB", &
sab_orb=sab_orb, calculate_forces=.TRUE.)
CALL dbcsr_deallocate_matrix_set(scrm)
!The pseudo-potential terms (only reponse density)
@ -3222,34 +3219,14 @@ CONTAINS
nder = 1
nimages = 1
NULLIFY (cell_to_index)
ALLOCATE (dbcsr_work_p(nspins, 1), dbcsr_work_h(nspins, 1))
DO ispin = 1, nspins
dbcsr_work_p(ispin, 1)%matrix => matrix_p_mp2(ispin)%matrix
dbcsr_work_h(ispin, 1)%matrix => scrm(ispin)%matrix
END DO
IF (ASSOCIATED(sac_ae)) THEN
CALL build_core_ae(dbcsr_work_h, dbcsr_work_p, force, &
virial, .TRUE., use_virial, nder, &
qs_kind_set, atomic_kind_set, particle_set, &
sab_orb, sac_ae, nimages, cell_to_index, "ORB")
END IF
CALL core_matrices(qs_env, dbcsr_work_h, dbcsr_work_p, .TRUE., nder)
IF (ASSOCIATED(sac_ppl)) THEN
CALL build_core_ppl(dbcsr_work_h, dbcsr_work_p, force, &
virial, .TRUE., use_virial, nder, &
qs_kind_set, atomic_kind_set, particle_set, &
sab_orb, sac_ppl, nimages, cell_to_index, "ORB")
END IF
IF (ASSOCIATED(sap_ppnl)) THEN
eps_ppnl = dft_control%qs_control%eps_ppnl
CALL build_core_ppnl(dbcsr_work_h, dbcsr_work_p, force, &
virial, .TRUE., use_virial, nder, &
qs_kind_set, atomic_kind_set, particle_set, &
sab_orb, sap_ppnl, eps_ppnl, nimages, cell_to_index, "ORB")
END IF
DEALLOCATE (dbcsr_work_p, dbcsr_work_h)
IF (use_virial) THEN