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

Add Z-matrix formalism for linear response (#3689)

Browse source
This commit is contained in:
Rangsiman Ketkaew 2024-09-23 15:53:39 +02:00 • committed by GitHub
parent 3206ad1121
commit 3f890bbe40
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
9 changed files with 545 additions and 65 deletions
Download patch file Download diff file Expand all files Collapse all files

7
src/input_cp2k_properties_dft.F
View file

@ -401,6 +401,13 @@ CONTAINS
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="Z_MATRIX_METHOD", &
description="Use Z_matrix method to solve the response equation", &
usage="Z_MATRIX_METHOD T", &
default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
NULLIFY (subsection)
CALL section_create(subsection, __LOCATION__, name="PRINT", &
description="print results of the magnetic dipole moment calculation", &

141
src/qs_dcdr.F
View file

@ -60,7 +60,9 @@ MODULE qs_dcdr
qs_kind_type
USE qs_linres_methods, ONLY: linres_solver
USE qs_linres_types, ONLY: dcdr_env_type,&
linres_control_type
get_polar_env,&
linres_control_type,&
polar_env_type
USE qs_mo_types, ONLY: get_mo_set,&
mo_set_type
USE qs_moments, ONLY: build_local_moment_matrix,&
@ -183,7 +185,10 @@ CONTAINS
!> \brief Build the operator for the position perturbation
!> \param dcdr_env ...
!> \param qs_env ...
!> \authors SL, ED
!> \authors Sandra Luber
!> Edward Ditler
!> Ravi Kumar
!> Rangsiman Ketkaew
! **************************************************************************************************
SUBROUTINE dcdr_build_op_dR(dcdr_env, qs_env)
@ -241,6 +246,11 @@ CONTAINS
! SL multiply by -1 for response solver (H-S<H> C + dR_coupled= - (op_dR)
CALL cp_fm_scale(-1.0_dp, dcdr_env%op_dR(ispin))
IF (dcdr_env%z_matrix_method) THEN
CALL cp_fm_to_fm(dcdr_env%op_dR(ispin), dcdr_env%matrix_m_alpha(dcdr_env%beta, ispin))
END IF
END DO
CALL dbcsr_deallocate_matrix_set(opdr_sym)
@ -357,7 +367,10 @@ CONTAINS
!> \brief Calculate atomic polar tensor
!> \param qs_env ...
!> \param dcdr_env ...
!> \author Edward Ditler
!> \authors Sandra Luber
!> Edward Ditler
!> Ravi Kumar
!> Rangsiman Ketkaew
! **************************************************************************************************
SUBROUTINE apt_dR(qs_env, dcdr_env)
TYPE(qs_environment_type), POINTER :: qs_env
@ -368,11 +381,14 @@ CONTAINS
INTEGER :: alpha, handle, ikind, ispin, nao, nmo
LOGICAL :: ghost
REAL(dp) :: apt_basis_derivative, &
apt_coeff_derivative, charge, f_spin
apt_coeff_derivative, charge, f_spin, &
temp1, temp2
REAL(dp), DIMENSION(:, :, :), POINTER :: apt_el, apt_nuc
TYPE(cp_fm_type) :: overlap1_MO, tmp_fm_like_mos
TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: dBerry_psi0, psi1_dBerry
TYPE(cp_fm_type), POINTER :: mo_coeff
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(polar_env_type), POINTER :: polar_env
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
apt_basis_derivative = 0._dp
@ -414,20 +430,48 @@ CONTAINS
CALL cp_fm_release(overlap1_MO)
DO alpha = 1, 3
! FIRST CONTRIBUTION: dCR * moments * mo
CALL cp_fm_set_all(tmp_fm_like_mos, 0._dp)
CALL dbcsr_desymmetrize(dcdr_env%matrix_s1(1)%matrix, dcdr_env%matrix_nosym_temp(1)%matrix)
CALL dbcsr_desymmetrize(dcdr_env%moments(alpha)%matrix, dcdr_env%matrix_nosym_temp(2)%matrix)
CALL dbcsr_add(dcdr_env%matrix_nosym_temp(1)%matrix, dcdr_env%matrix_nosym_temp(2)%matrix, &
-dcdr_env%ref_point(alpha), 1._dp)
IF (.NOT. dcdr_env%z_matrix_method) THEN
CALL cp_dbcsr_sm_fm_multiply(dcdr_env%matrix_nosym_temp(1)%matrix, dcdr_env%dCR_prime(ispin), &
tmp_fm_like_mos, ncol=nmo)
CALL cp_fm_trace(mo_coeff, tmp_fm_like_mos, apt_coeff_derivative)
! FIRST CONTRIBUTION: dCR * moments * mo
CALL cp_fm_set_all(tmp_fm_like_mos, 0._dp)
CALL dbcsr_desymmetrize(dcdr_env%matrix_s1(1)%matrix, dcdr_env%matrix_nosym_temp(1)%matrix)
CALL dbcsr_desymmetrize(dcdr_env%moments(alpha)%matrix, dcdr_env%matrix_nosym_temp(2)%matrix)
CALL dbcsr_add(dcdr_env%matrix_nosym_temp(1)%matrix, dcdr_env%matrix_nosym_temp(2)%matrix, &
-dcdr_env%ref_point(alpha), 1._dp)
apt_coeff_derivative = (-2._dp)*f_spin*apt_coeff_derivative
apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) &
= apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) + apt_coeff_derivative
CALL cp_dbcsr_sm_fm_multiply(dcdr_env%matrix_nosym_temp(1)%matrix, dcdr_env%dCR_prime(ispin), &
tmp_fm_like_mos, ncol=nmo)
CALL cp_fm_trace(mo_coeff, tmp_fm_like_mos, apt_coeff_derivative)
apt_coeff_derivative = (-2._dp)*f_spin*apt_coeff_derivative
apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) &
= apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) + apt_coeff_derivative
ELSE
CALL get_qs_env(qs_env=qs_env, polar_env=polar_env)
CALL get_polar_env(polar_env=polar_env, psi1_dBerry=psi1_dBerry, &
dBerry_psi0=dBerry_psi0)
! Note that here dcdr_env%dCR_prime contains only occ-occ block contribution,
! dcdr_env%dCR(ispin) is zero because we didn't run response calculation for dcdR.
CALL cp_fm_trace(dBerry_psi0(alpha, ispin), &
dcdr_env%dCR_prime(ispin), &
temp1)
CALL cp_fm_trace(dcdr_env%matrix_m_alpha(dcdr_env%beta, ispin), &
psi1_dBerry(alpha, ispin), &
temp2)
apt_coeff_derivative = temp1 - temp2
! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
! - apt_coeff_derivative , here the trace is negative to compensate the
! -ve sign in APTs= - 2 Z. M_alpha
apt_coeff_derivative = (-2._dp)*f_spin*apt_coeff_derivative
apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) &
= apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) + apt_coeff_derivative
END IF
! SECOND CONTRIBUTION: We assemble all combinations of r_i, d(chi)/d(idir)
! difdip contains derivatives with respect to atom dcdr_env%lambda
@ -442,6 +486,7 @@ CONTAINS
apt_basis_derivative = -f_spin*apt_basis_derivative
apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) = &
apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) + apt_basis_derivative
END DO ! alpha
CALL cp_fm_release(tmp_fm_like_mos)
@ -466,7 +511,9 @@ CONTAINS
!> \brief Calculate atomic polar tensor using the localized dipole operator
!> \param qs_env ...
!> \param dcdr_env ...
!> \author Edward Ditler
!> \authors Edward Ditler
!> Ravi Kumar
!> Rangsiman Ketkaew
! **************************************************************************************************
SUBROUTINE apt_dR_localization(qs_env, dcdr_env)
TYPE(qs_environment_type), POINTER :: qs_env
@ -485,16 +532,19 @@ CONTAINS
apt_coeff_derivative, charge, f_spin, &
smallest_r, this_factor, tmp_aptcontr, &
tmp_r
REAL(dp), ALLOCATABLE, DIMENSION(:) :: diagonal_elements
REAL(dp), ALLOCATABLE, DIMENSION(:) :: diagonal_elements, diagonal_elements2
REAL(dp), DIMENSION(3) :: distance, r_shifted
REAL(dp), DIMENSION(:, :, :), POINTER :: apt_el, apt_nuc
REAL(dp), DIMENSION(:, :, :, :), POINTER :: apt_center, apt_subset
TYPE(cell_type), POINTER :: cell
TYPE(cp_2d_r_p_type), DIMENSION(:), POINTER :: centers_set
TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: dBerry_psi0, psi1_dBerry
TYPE(cp_fm_type), POINTER :: mo_coeff, overlap1_MO, tmp_fm, &
tmp_fm_like_mos, tmp_fm_momo
tmp_fm_like_mos, tmp_fm_momo, &
tmp_fm_momo2
TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(polar_env_type), POINTER :: polar_env
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
CALL timeset(routineN, handle)
@ -581,30 +631,53 @@ CONTAINS
CALL cp_fm_release(overlap1_MO)
ALLOCATE (diagonal_elements(nmo))
ALLOCATE (diagonal_elements2(nmo))
! Allocate temporary matrices
ALLOCATE (tmp_fm)
ALLOCATE (tmp_fm_momo)
ALLOCATE (tmp_fm_momo2)
CALL cp_fm_create(tmp_fm, dcdr_env%likemos_fm_struct(ispin)%struct)
CALL cp_fm_create(tmp_fm_momo, dcdr_env%momo_fm_struct(ispin)%struct)
CALL cp_fm_create(tmp_fm_momo2, dcdr_env%momo_fm_struct(ispin)%struct)
! FIRST CONTRIBUTION: dCR * moments * mo
this_factor = -2._dp*f_spin
DO alpha = 1, 3
DO icenter = 1, dcdr_env%nbr_center(ispin)
CALL dbcsr_set(dcdr_env%moments(alpha)%matrix, 0.0_dp)
CALL build_local_moment_matrix(qs_env, dcdr_env%moments, 1, &
ref_point=centers_set(ispin)%array(1:3, icenter))
CALL multiply_localization(ao_matrix=dcdr_env%moments(alpha)%matrix, &
mo_coeff=dcdr_env%dCR_prime(ispin), work=tmp_fm, nmo=nmo, &
icenter=icenter, &
res=tmp_fm_like_mos)
END DO
IF (.NOT. dcdr_env%z_matrix_method) THEN
CALL parallel_gemm("T", "N", nmo, nmo, nao, &
1.0_dp, mo_coeff, tmp_fm_like_mos, &
0.0_dp, tmp_fm_momo)
CALL cp_fm_get_diag(tmp_fm_momo, diagonal_elements)
DO icenter = 1, dcdr_env%nbr_center(ispin)
CALL dbcsr_set(dcdr_env%moments(alpha)%matrix, 0.0_dp)
CALL build_local_moment_matrix(qs_env, dcdr_env%moments, 1, &
ref_point=centers_set(ispin)%array(1:3, icenter))
CALL multiply_localization(ao_matrix=dcdr_env%moments(alpha)%matrix, &
mo_coeff=dcdr_env%dCR_prime(ispin), work=tmp_fm, nmo=nmo, &
icenter=icenter, &
res=tmp_fm_like_mos)
END DO
CALL parallel_gemm("T", "N", nmo, nmo, nao, &
1.0_dp, mo_coeff, tmp_fm_like_mos, &
0.0_dp, tmp_fm_momo)
CALL cp_fm_get_diag(tmp_fm_momo, diagonal_elements)
ELSE
CALL get_qs_env(qs_env=qs_env, polar_env=polar_env)
CALL get_polar_env(polar_env=polar_env, psi1_dBerry=psi1_dBerry, &
dBerry_psi0=dBerry_psi0)
CALL parallel_gemm("T", "N", nmo, nmo, nao, &
1.0_dp, dcdr_env%dCR_prime(ispin), dBerry_psi0(alpha, ispin), &
0.0_dp, tmp_fm_momo)
CALL cp_fm_get_diag(tmp_fm_momo, diagonal_elements)
CALL parallel_gemm("T", "N", nmo, nmo, nao, &
1.0_dp, dcdr_env%matrix_m_alpha(dcdr_env%beta, ispin), &
psi1_dBerry(alpha, ispin), 0.0_dp, tmp_fm_momo2)
CALL cp_fm_get_diag(tmp_fm_momo2, diagonal_elements2)
diagonal_elements(:) = diagonal_elements(:) - diagonal_elements2(:)
END IF
DO icenter = 1, dcdr_env%nbr_center(ispin)
map_atom = mapping_wannier_atom(icenter, ispin)
@ -667,14 +740,17 @@ CONTAINS
END DO ! alpha
DEALLOCATE (diagonal_elements)
DEALLOCATE (diagonal_elements2)
CALL cp_fm_release(tmp_fm)
CALL cp_fm_release(tmp_fm_like_mos)
CALL cp_fm_release(tmp_fm_momo)
CALL cp_fm_release(tmp_fm_momo2)
DEALLOCATE (overlap1_MO)
DEALLOCATE (tmp_fm)
DEALLOCATE (tmp_fm_like_mos)
DEALLOCATE (tmp_fm_momo)
DEALLOCATE (tmp_fm_momo2)
END DO !ispin
! Finally the nuclear contribution: nuclear charge * Kronecker_delta_{dcdr_env%beta,i}
@ -695,3 +771,4 @@ CONTAINS
END SUBROUTINE apt_dR_localization
END MODULE qs_dcdr

16
src/qs_dcdr_utils.F
View file

@ -732,6 +732,8 @@ CONTAINS
CALL section_vals_val_get(dcdr_section, "DISTRIBUTED_ORIGIN", l_val=dcdr_env%distributed_origin)
CALL section_vals_val_get(loc_section, "_SECTION_PARAMETERS_", l_val=dcdr_env%localized_psi0)
CALL section_vals_val_get(lr_section, "RESTART", l_val=qs_env%linres_control%linres_restart)
CALL section_vals_val_get(dcdr_section, "Z_MATRIX_METHOD", l_val=dcdr_env%z_matrix_method)
dcdr_env%ref_point = 0._dp
! List of atoms
@ -844,6 +846,16 @@ CONTAINS
CALL cp_fm_to_fm(mo_coeff, dcdr_env%mo_coeff(ispin))
END DO
IF (dcdr_env%z_matrix_method) THEN
ALLOCATE (dcdr_env%matrix_m_alpha(3, nspins))
DO i = 1, 3
DO ispin = 1, nspins
CALL cp_fm_create(dcdr_env%matrix_m_alpha(i, ispin), dcdr_env%likemos_fm_struct(1)%struct)
CALL cp_fm_set_all(dcdr_env%matrix_m_alpha(i, ispin), 0.0_dp)
END DO
END DO
END IF
! DBCSR matrices
NULLIFY (dcdr_env%hamiltonian1)
NULLIFY (dcdr_env%moments)
@ -1038,6 +1050,10 @@ CONTAINS
CALL cp_fm_release(dcdr_env%chc)
CALL cp_fm_release(dcdr_env%op_dR)
IF (dcdr_env%z_matrix_method) THEN
CALL cp_fm_release(dcdr_env%matrix_m_alpha)
END IF
! DBCSR matrices
CALL dbcsr_deallocate_matrix_set(dcdr_env%perturbed_dm_correction)
CALL dbcsr_deallocate_matrix_set(dcdr_env%hamiltonian1)

90
src/qs_linres_module.F
View file

@ -80,18 +80,16 @@ MODULE qs_linres_module
nmr_env_init
USE qs_linres_op, ONLY: current_operators,&
issc_operators,&
polar_operators
polar_operators,&
polar_operators_local,&
polar_operators_local_wannier
USE qs_linres_polar_utils, ONLY: polar_env_init,&
polar_polar,&
polar_print,&
polar_response
USE qs_linres_types, ONLY: current_env_type,&
dcdr_env_type,&
epr_env_type,&
issc_env_type,&
linres_control_type,&
nmr_env_type,&
vcd_env_type
USE qs_linres_types, ONLY: &
current_env_type, dcdr_env_type, epr_env_type, get_polar_env, issc_env_type, &
linres_control_type, nmr_env_type, polar_env_type, vcd_env_type
USE qs_mfp, ONLY: mfp_aat,&
mfp_build_operator_gauge_dependent,&
mfp_build_operator_gauge_independent,&
@ -122,7 +120,6 @@ MODULE qs_linres_module
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_linres_module'
CONTAINS
! *****************************************************************************
!> \brief Calculates the derivatives of the MO coefficients dC/dV^lambda_beta
!> wrt to nuclear velocities. The derivative is indexed by `beta`, the
@ -231,31 +228,72 @@ CONTAINS
INTEGER :: beta, latom
TYPE(dcdr_env_type) :: dcdr_env
TYPE(polar_env_type), POINTER :: polar_env
CALL cite_reference(Ditler2021)
CALL dcdr_env_init(dcdr_env, qs_env)
DO latom = 1, SIZE(dcdr_env%list_of_atoms)
dcdr_env%lambda = dcdr_env%list_of_atoms(latom)
CALL prepare_per_atom(dcdr_env, qs_env)
DO beta = 1, 3 ! in every direction
dcdr_env%beta = beta
dcdr_env%deltaR(dcdr_env%beta, dcdr_env%lambda) = 1._dp
IF (.NOT. dcdr_env%z_matrix_method) THEN
CALL dcdr_build_op_dR(dcdr_env, qs_env)
CALL dcdr_response_dR(dcdr_env, p_env, qs_env)
DO latom = 1, SIZE(dcdr_env%list_of_atoms)
dcdr_env%lambda = dcdr_env%list_of_atoms(latom)
CALL prepare_per_atom(dcdr_env, qs_env)
IF (.NOT. dcdr_env%localized_psi0) THEN
CALL apt_dR(qs_env, dcdr_env)
ELSE IF (dcdr_env%localized_psi0) THEN
CALL apt_dR_localization(qs_env, dcdr_env)
END IF
DO beta = 1, 3 ! in every direction
dcdr_env%beta = beta
dcdr_env%deltaR(dcdr_env%beta, dcdr_env%lambda) = 1._dp
END DO !beta
CALL dcdr_build_op_dR(dcdr_env, qs_env)
CALL dcdr_response_dR(dcdr_env, p_env, qs_env)
dcdr_env%apt_total_dcdr(:, :, dcdr_env%lambda) = &
dcdr_env%apt_el_dcdr(:, :, dcdr_env%lambda) + dcdr_env%apt_nuc_dcdr(:, :, dcdr_env%lambda)
END DO !lambda
IF (.NOT. dcdr_env%localized_psi0) THEN
CALL apt_dR(qs_env, dcdr_env)
ELSE IF (dcdr_env%localized_psi0) THEN
CALL apt_dR_localization(qs_env, dcdr_env)
END IF
END DO !beta
dcdr_env%apt_total_dcdr(:, :, dcdr_env%lambda) = &
dcdr_env%apt_el_dcdr(:, :, dcdr_env%lambda) + dcdr_env%apt_nuc_dcdr(:, :, dcdr_env%lambda)
END DO !lambda
ELSE
CALL polar_env_init(qs_env)
CALL get_qs_env(qs_env=qs_env, polar_env=polar_env)
CALL get_polar_env(polar_env=polar_env)
IF (.NOT. dcdr_env%localized_psi0) THEN
CALL polar_operators_local(qs_env)
ELSE
CALL polar_operators_local_wannier(qs_env, dcdr_env)
END IF
polar_env%do_periodic = .FALSE.
CALL polar_response(p_env, qs_env)
DO latom = 1, SIZE(dcdr_env%list_of_atoms)
dcdr_env%lambda = dcdr_env%list_of_atoms(latom)
CALL prepare_per_atom(dcdr_env, qs_env)
DO beta = 1, 3 ! in every direction
dcdr_env%beta = beta
dcdr_env%deltaR(dcdr_env%beta, dcdr_env%lambda) = 1._dp
CALL dcdr_build_op_dR(dcdr_env, qs_env)
IF (.NOT. dcdr_env%localized_psi0) THEN
CALL apt_dR(qs_env, dcdr_env)
ELSE
CALL apt_dR_localization(qs_env, dcdr_env)
END IF
END DO !beta
dcdr_env%apt_total_dcdr(:, :, dcdr_env%lambda) = &
dcdr_env%apt_el_dcdr(:, :, dcdr_env%lambda) + dcdr_env%apt_nuc_dcdr(:, :, dcdr_env%lambda)
END DO !lambda
END IF
CALL dcdr_print(dcdr_env, qs_env)
CALL dcdr_env_cleanup(qs_env, dcdr_env)

163
src/qs_linres_op.F
View file

@ -59,10 +59,14 @@ MODULE qs_linres_op
USE kinds, ONLY: dp
USE mathconstants, ONLY: twopi
USE message_passing, ONLY: mp_para_env_type
USE molecule_types, ONLY: molecule_of_atom,&
molecule_type
USE orbital_pointers, ONLY: coset
USE parallel_gemm_api, ONLY: parallel_gemm
USE particle_methods, ONLY: get_particle_set
USE particle_types, ONLY: particle_type
USE qs_dcdr_utils, ONLY: multiply_localization,&
shift_wannier_into_cell
USE qs_elec_field, ONLY: build_efg_matrix
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
@ -70,6 +74,7 @@ MODULE qs_linres_op
USE qs_kind_types, ONLY: get_qs_kind_set,&
qs_kind_type
USE qs_linres_types, ONLY: current_env_type,&
dcdr_env_type,&
get_current_env,&
get_issc_env,&
get_polar_env,&
@ -91,7 +96,8 @@ MODULE qs_linres_op
PRIVATE
PUBLIC :: current_operators, issc_operators, fac_vecp, ind_m2, set_vecp, set_vecp_rev, &
fm_scale_by_pbc_AC, polar_operators
fm_scale_by_pbc_AC, polar_operators, polar_operators_local, &
polar_operators_local_wannier, polar_operators_berry
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_linres_op'
@ -1009,6 +1015,161 @@ CONTAINS
END SUBROUTINE polar_operators_local
! **************************************************************************************************
!> \brief Calculate the dipole operator referenced at the Wannier centers in the MO basis
!> \param qs_env ...
!> \param dcdr_env ...
!> \par History
!> 01.2013 created [SL]
!> 06.2018 polar_env integrated into qs_env (MK)
!> \authors Ravi Kumar
!> Rangsiman Ketkaew
! **************************************************************************************************
SUBROUTINE polar_operators_local_wannier(qs_env, dcdr_env)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(dcdr_env_type) :: dcdr_env
CHARACTER(LEN=*), PARAMETER :: routineN = 'polar_operators_local_wannier'
INTEGER :: alpha, handle, i, icenter, ispin, &
map_atom, map_molecule, &
max_nbr_center, nao, natom, nmo, &
nsubset
INTEGER, ALLOCATABLE, DIMENSION(:) :: mapping_atom_molecule
INTEGER, ALLOCATABLE, DIMENSION(:, :) :: mapping_wannier_atom
REAL(dp) :: f_spin, smallest_r, tmp_r
REAL(dp), DIMENSION(3) :: distance, r_shifted
REAL(dp), DIMENSION(:, :, :), POINTER :: apt_el, apt_nuc
REAL(dp), DIMENSION(:, :, :, :), POINTER :: apt_center, apt_subset
TYPE(cell_type), POINTER :: cell
TYPE(cp_2d_r_p_type), DIMENSION(:), POINTER :: centers_set
TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: dBerry_psi0
TYPE(cp_fm_type), POINTER :: mo_coeff, overlap1_MO, tmp_fm, &
tmp_fm_like_mos, tmp_fm_momo
TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(polar_env_type), POINTER :: polar_env
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
CALL timeset(routineN, handle)
NULLIFY (qs_kind_set, particle_set, molecule_set, cell)
CALL get_qs_env(qs_env=qs_env, &
qs_kind_set=qs_kind_set, &
particle_set=particle_set, &
molecule_set=molecule_set, &
polar_env=polar_env, &
cell=cell)
CALL get_polar_env(polar_env=polar_env, dBerry_psi0=dBerry_psi0)
nsubset = SIZE(molecule_set)
natom = SIZE(particle_set)
apt_el => dcdr_env%apt_el_dcdr
apt_nuc => dcdr_env%apt_nuc_dcdr
apt_subset => dcdr_env%apt_el_dcdr_per_subset
apt_center => dcdr_env%apt_el_dcdr_per_center
! Map wannier functions to atoms
IF (dcdr_env%nspins == 1) THEN
max_nbr_center = dcdr_env%nbr_center(1)
ELSE
max_nbr_center = MAX(dcdr_env%nbr_center(1), dcdr_env%nbr_center(2))
END IF
ALLOCATE (mapping_wannier_atom(max_nbr_center, dcdr_env%nspins))
ALLOCATE (mapping_atom_molecule(natom))
centers_set => dcdr_env%centers_set
DO ispin = 1, dcdr_env%nspins
DO icenter = 1, dcdr_env%nbr_center(ispin)
! For every center we check which atom is closest
CALL shift_wannier_into_cell(r=centers_set(ispin)%array(1:3, icenter), &
cell=cell, &
r_shifted=r_shifted)
smallest_r = HUGE(0._dp)
DO i = 1, natom
distance = pbc(r_shifted, particle_set(i)%r(1:3), cell)
tmp_r = SUM(distance**2)
IF (tmp_r < smallest_r) THEN
mapping_wannier_atom(icenter, ispin) = i
smallest_r = tmp_r
END IF
END DO
END DO
! Map atoms to molecules
CALL molecule_of_atom(molecule_set, atom_to_mol=mapping_atom_molecule)
IF (dcdr_env%lambda == 1 .AND. dcdr_env%beta == 1) THEN
DO icenter = 1, dcdr_env%nbr_center(ispin)
map_atom = mapping_wannier_atom(icenter, ispin)
map_molecule = mapping_atom_molecule(map_atom)
END DO
END IF
END DO !ispin
nao = dcdr_env%nao
f_spin = 2._dp/dcdr_env%nspins
DO ispin = 1, dcdr_env%nspins
! Compute S^(1,R)_(ij)
ALLOCATE (tmp_fm_like_mos)
ALLOCATE (overlap1_MO)
CALL cp_fm_create(tmp_fm_like_mos, dcdr_env%likemos_fm_struct(ispin)%struct)
CALL cp_fm_create(overlap1_MO, dcdr_env%momo_fm_struct(ispin)%struct)
nmo = dcdr_env%nmo(ispin)
mo_coeff => dcdr_env%mo_coeff(ispin)
CALL cp_fm_set_all(tmp_fm_like_mos, 0.0_dp)
CALL cp_fm_scale_and_add(0._dp, dcdr_env%dCR_prime(ispin), 1._dp, dcdr_env%dCR(ispin))
! CALL cp_dbcsr_sm_fm_multiply(dcdr_env%matrix_s1(dcdr_env%beta + 1)%matrix, mo_coeff, &
! tmp_fm_like_mos, ncol=nmo)
CALL parallel_gemm("T", "N", nmo, nmo, nao, &
1.0_dp, mo_coeff, tmp_fm_like_mos, &
0.0_dp, overlap1_MO)
! C^1 <- -dCR - 0.5 * mo_coeff @ S1_ij
! We get the negative of the coefficients out of the linres solver
! And apply the constant correction due to the overlap derivative.
CALL parallel_gemm("N", "N", nao, nmo, nmo, &
-0.5_dp, mo_coeff, overlap1_MO, &
-1.0_dp, dcdr_env%dCR_prime(ispin))
CALL cp_fm_release(overlap1_MO)
! Allocate temporary matrices
ALLOCATE (tmp_fm)
ALLOCATE (tmp_fm_momo)
CALL cp_fm_create(tmp_fm, dcdr_env%likemos_fm_struct(ispin)%struct)
CALL cp_fm_create(tmp_fm_momo, dcdr_env%momo_fm_struct(ispin)%struct)
! this_factor = -2._dp*f_spin
DO alpha = 1, 3
DO icenter = 1, dcdr_env%nbr_center(ispin)
CALL dbcsr_set(dcdr_env%moments(alpha)%matrix, 0.0_dp)
CALL build_local_moment_matrix(qs_env, dcdr_env%moments, 1, &
ref_point=centers_set(ispin)%array(1:3, icenter))
CALL multiply_localization(ao_matrix=dcdr_env%moments(alpha)%matrix, &
mo_coeff=mo_coeff, work=tmp_fm, nmo=nmo, &
icenter=icenter, &
res=dBerry_psi0(alpha, ispin))
END DO
END DO
CALL cp_fm_release(tmp_fm)
CALL cp_fm_release(tmp_fm_like_mos)
CALL cp_fm_release(tmp_fm_momo)
DEALLOCATE (overlap1_MO)
DEALLOCATE (tmp_fm)
DEALLOCATE (tmp_fm_like_mos)
DEALLOCATE (tmp_fm_momo)
END DO !ispin
! And deallocate all the things!
CALL timestop(handle)
END SUBROUTINE polar_operators_local_wannier
! **************************************************************************************************
!> \brief Calculate the local dipole operator in the AO basis
!> afterwards multiply with the ground state MO coefficients

3
src/qs_linres_types.F
View file

@ -277,6 +277,8 @@ MODULE qs_linres_types
TYPE(cp_fm_type), DIMENSION(:), POINTER :: dCR_prime => NULL()
TYPE(cp_fm_type), DIMENSION(:), POINTER :: op_dR => NULL()
TYPE(cp_fm_type), DIMENSION(:), POINTER :: chc => NULL()
TYPE(cp_fm_type), DIMENSION(:), POINTER :: ch1c => NULL()
TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: matrix_m_alpha => NULL()
CHARACTER(LEN=30) :: orb_center_name = ""
TYPE(cp_2d_i_p_type), DIMENSION(:), POINTER :: center_list => NULL()
TYPE(cp_2d_r_p_type), DIMENSION(:), POINTER :: centers_set => NULL()
@ -287,6 +289,7 @@ MODULE qs_linres_types
LOGICAL :: localized_psi0 = .FALSE.
INTEGER, POINTER :: list_of_atoms(:) => NULL()
LOGICAL :: distributed_origin = .FALSE.
LOGICAL :: z_matrix_method = .FALSE.
TYPE(cp_fm_struct_type), POINTER :: aoao_fm_struct => NULL()
TYPE(cp_fm_struct_type), POINTER :: homohomo_fm_struct => NULL()
TYPE(cp_fm_struct_p_type), DIMENSION(:), POINTER :: momo_fm_struct => NULL()

14
tests/QS/regtest-dcdr/TEST_FILES
View file

@ -3,10 +3,12 @@
# e.g. 0 means do not compare anything, running is enough
# 1 compares the last total energy in the file
# for details see cp2k/tools/do_regtest
h2o_apt.inp 96 1e-06 -0.879586
h2o_apt_uks.inp 96 1e-06 -0.879582
h2o_apt_loc.inp 96 1e-06 -0.879586
h2o_apt_pbc.inp 96 1e-06 -0.938707
h2o_apt_pbc_loc.inp 96 1e-06 -0.938770
h2o_aat.inp 100 1e-06 -0.857427
h2o_apt.inp 96 1e-06 -0.879586
h2o_apt_z-matrix.inp 96 1e-06 -0.882548
h2o_apt_uks.inp 96 1e-06 -0.879582
h2o_apt_uks_z-matrix.inp 96 1e-06 -0.882544
h2o_apt_loc.inp 96 1e-06 -0.879586
h2o_apt_pbc.inp 96 1e-06 -0.938707
h2o_apt_pbc_loc.inp 96 1e-06 -0.938770
h2o_aat.inp 100 1e-06 -0.857427
#EOF

89
tests/QS/regtest-dcdr/h2o_apt_uks_z-matrix.inp Normal file
View file

@ -0,0 +1,89 @@
###################################
@SET RUN_TYPE ENERGY_FORCE
@SET CUTOFF 200
@SET FUNCTIONAL LDA
@SET PRINT_LEVEL MEDIUM
@SET BASIS_SET_FILE_NAME GTH_BASIS_SETS
@SET BASIS_SET SZV-GTH
@SET EPS_SCF 1.08E-5
@SET EPS_LINRES 5.0E-5
###################################
&GLOBAL
PRINT_LEVEL $PRINT_LEVEL
PROJECT second
RUN_TYPE $RUN_TYPE
&END GLOBAL
&FORCE_EVAL
METHOD Quickstep
&DFT
BASIS_SET_FILE_NAME $BASIS_SET_FILE_NAME
CHARGE 0
MULTIPLICITY 1
POTENTIAL_FILE_NAME POTENTIAL
UKS T
&MGRID
CUTOFF $CUTOFF
NGRIDS 4
&END MGRID
&POISSON
PERIODIC NONE
POISSON_SOLVER ANALYTIC
&END POISSON
&PRINT
&MOMENTS
PERIODIC FALSE
&END MOMENTS
&END PRINT
&QS
EXTRAPOLATION ASPC
EXTRAPOLATION_ORDER 3
METHOD GPW
&END QS
&SCF
EPS_SCF $EPS_SCF
SCF_GUESS ATOMIC
&OT
PRECONDITIONER FULL_SINGLE_INVERSE
&END OT
&END SCF
&XC
&XC_FUNCTIONAL $FUNCTIONAL
&END XC_FUNCTIONAL
&END XC
&END DFT
&PROPERTIES
&LINRES
EPS $EPS_LINRES
MAX_ITER 1000
PRECONDITIONER FULL_SINGLE_INVERSE
&DCDR
Z_MATRIX_METHOD T
&PRINT
&APT
FILENAME __STD_OUT__
&END APT
&END PRINT
&END DCDR
&PRINT
&PROGRAM_RUN_INFO
&END PROGRAM_RUN_INFO
&END PRINT
&END LINRES
&END PROPERTIES
&SUBSYS
&CELL
ABC [angstrom] 5.0 5.0 5.0
PERIODIC NONE
&END CELL
&COORD
O 0.000000 0.000000 0.000000
H 0.000000 0.769665 -0.591648
H 0.000000 -0.769665 -0.591648
&END COORD
&KIND DEFAULT
BASIS_SET $BASIS_SET
POTENTIAL GTH-$FUNCTIONAL
&END KIND
&END SUBSYS
&END FORCE_EVAL

87
tests/QS/regtest-dcdr/h2o_apt_z-matrix.inp Normal file
View file

@ -0,0 +1,87 @@
###################################
@SET RUN_TYPE ENERGY_FORCE
@SET CUTOFF 200
@SET FUNCTIONAL LDA
@SET PRINT_LEVEL MEDIUM
@SET BASIS_SET_FILE_NAME GTH_BASIS_SETS
@SET BASIS_SET SZV-GTH
@SET EPS_SCF 1.08E-5
@SET EPS_LINRES 5.0E-5
###################################
&GLOBAL
PRINT_LEVEL $PRINT_LEVEL
PROJECT second
RUN_TYPE $RUN_TYPE
&END GLOBAL
&FORCE_EVAL
METHOD Quickstep
&DFT
BASIS_SET_FILE_NAME $BASIS_SET_FILE_NAME
CHARGE 0
POTENTIAL_FILE_NAME POTENTIAL
&MGRID
CUTOFF $CUTOFF
NGRIDS 1
&END MGRID
&POISSON
PERIODIC NONE
POISSON_SOLVER ANALYTIC
&END POISSON
&PRINT
&MOMENTS
PERIODIC FALSE
&END MOMENTS
&END PRINT
&QS
EXTRAPOLATION ASPC
EXTRAPOLATION_ORDER 3
METHOD GPW
&END QS
&SCF
EPS_SCF $EPS_SCF
SCF_GUESS ATOMIC
&OT
PRECONDITIONER FULL_SINGLE_INVERSE
&END OT
&END SCF
&XC
&XC_FUNCTIONAL $FUNCTIONAL
&END XC_FUNCTIONAL
&END XC
&END DFT
&PROPERTIES
&LINRES
EPS $EPS_LINRES
MAX_ITER 1000
PRECONDITIONER FULL_SINGLE_INVERSE
&DCDR
Z_MATRIX_METHOD T
&PRINT
&APT
FILENAME __STD_OUT__
&END APT
&END PRINT
&END DCDR
&PRINT
&PROGRAM_RUN_INFO
&END PROGRAM_RUN_INFO
&END PRINT
&END LINRES
&END PROPERTIES
&SUBSYS
&CELL
ABC [angstrom] 5.0 5.0 5.0
PERIODIC NONE
&END CELL
&COORD
O 0.000000 0.000000 0.000000
H 0.000000 0.769665 -0.591648
H 0.000000 -0.769665 -0.591648
&END COORD
&KIND DEFAULT
BASIS_SET $BASIS_SET
POTENTIAL GTH-$FUNCTIONAL
&END KIND
&END SUBSYS
&END FORCE_EVAL