From 3f890bbe40bbf237a88c26c2c87ca15893cb8c80 Mon Sep 17 00:00:00 2001 From: Rangsiman Ketkaew Date: Mon, 23 Sep 2024 15:53:39 +0200 Subject: [PATCH] Add Z-matrix formalism for linear response (#3689) --- src/input_cp2k_properties_dft.F | 7 + src/qs_dcdr.F | 141 +++++++++++---- src/qs_dcdr_utils.F | 16 ++ src/qs_linres_module.F | 90 +++++++--- src/qs_linres_op.F | 163 +++++++++++++++++- src/qs_linres_types.F | 3 + tests/QS/regtest-dcdr/TEST_FILES | 14 +- .../QS/regtest-dcdr/h2o_apt_uks_z-matrix.inp | 89 ++++++++++ tests/QS/regtest-dcdr/h2o_apt_z-matrix.inp | 87 ++++++++++ 9 files changed, 545 insertions(+), 65 deletions(-) create mode 100644 tests/QS/regtest-dcdr/h2o_apt_uks_z-matrix.inp create mode 100644 tests/QS/regtest-dcdr/h2o_apt_z-matrix.inp diff --git a/src/input_cp2k_properties_dft.F b/src/input_cp2k_properties_dft.F index 91c0f4ad77..3a05994b27 100644 --- a/src/input_cp2k_properties_dft.F +++ b/src/input_cp2k_properties_dft.F @@ -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", & diff --git a/src/qs_dcdr.F b/src/qs_dcdr.F index 0b3a30aec9..527fc0318d 100644 --- a/src/qs_dcdr.F +++ b/src/qs_dcdr.F @@ -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 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 + diff --git a/src/qs_dcdr_utils.F b/src/qs_dcdr_utils.F index 1ece6d9ba4..b87a4bb216 100644 --- a/src/qs_dcdr_utils.F +++ b/src/qs_dcdr_utils.F @@ -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) diff --git a/src/qs_linres_module.F b/src/qs_linres_module.F index 470c24d9f9..ab84afb35a 100644 --- a/src/qs_linres_module.F +++ b/src/qs_linres_module.F @@ -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) diff --git a/src/qs_linres_op.F b/src/qs_linres_op.F index 00b28f349b..24cb343174 100644 --- a/src/qs_linres_op.F +++ b/src/qs_linres_op.F @@ -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 diff --git a/src/qs_linres_types.F b/src/qs_linres_types.F index 8ecc4c8bc7..3d89220b17 100644 --- a/src/qs_linres_types.F +++ b/src/qs_linres_types.F @@ -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() diff --git a/tests/QS/regtest-dcdr/TEST_FILES b/tests/QS/regtest-dcdr/TEST_FILES index eca8733e08..c1b3425c51 100644 --- a/tests/QS/regtest-dcdr/TEST_FILES +++ b/tests/QS/regtest-dcdr/TEST_FILES @@ -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 diff --git a/tests/QS/regtest-dcdr/h2o_apt_uks_z-matrix.inp b/tests/QS/regtest-dcdr/h2o_apt_uks_z-matrix.inp new file mode 100644 index 0000000000..aa7276f224 --- /dev/null +++ b/tests/QS/regtest-dcdr/h2o_apt_uks_z-matrix.inp @@ -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 diff --git a/tests/QS/regtest-dcdr/h2o_apt_z-matrix.inp b/tests/QS/regtest-dcdr/h2o_apt_z-matrix.inp new file mode 100644 index 0000000000..626a15c39b --- /dev/null +++ b/tests/QS/regtest-dcdr/h2o_apt_z-matrix.inp @@ -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