From e709bc9b963055c87db2d7fb03b6a1948d4b7f58 Mon Sep 17 00:00:00 2001 From: Juerg Hutter Date: Fri, 18 Aug 2023 15:18:36 +0200 Subject: [PATCH] GAPW TDDFPT (#2932) --- src/qs_environment.F | 2 +- src/qs_fxc.F | 2 + src/qs_ks_atom.F | 23 +- src/qs_ks_reference.F | 21 +- src/qs_p_env_methods.F | 3 +- src/qs_rho0_ggrid.F | 23 +- src/qs_rho0_methods.F | 8 +- src/qs_tddfpt2_fhxc.F | 19 +- src/qs_tddfpt2_fhxc_forces.F | 42 ++- src/qs_tddfpt2_forces.F | 13 +- src/qs_tddfpt2_subgroups.F | 3 +- src/qs_tddfpt2_types.F | 3 +- src/qs_vxc_atom.F | 162 ++++++++---- src/response_solver.F | 242 +----------------- src/xc/xc_atom.F | 3 +- tests/QS/regtest-tddfpt-force-gapw/TEST_FILES | 23 +- .../{ch2o_f01_pbe.inp => h2o_f01_pbe.inp} | 22 +- ..._pbe_gapwxc.inp => h2o_f01_pbe_gapwxc.inp} | 0 ...lomb_only.inp => h2o_f02_coulomb_only.inp} | 0 .../QS/regtest-tddfpt-force-gapw/h2o_t06.inp | 2 +- .../QS/regtest-tddfpt-force-gapw/h2o_t09.inp | 90 +++++++ .../QS/regtest-tddfpt-force-gapw/h2o_t10.inp | 90 +++++++ .../QS/regtest-tddfpt-force-gapw/h2o_t11.inp | 90 +++++++ .../QS/regtest-tddfpt-force-gapw/h2o_t12.inp | 90 +++++++ .../H2O_GAPW_1.inp | 0 .../H2O_GAPW_1_triplet.inp | 16 +- .../H2O_GAPW_2.inp | 0 .../H2O_GAPW_3.inp | 0 .../H2O_GAPW_4.inp | 0 .../H2O_GAPW_XC_1.inp | 0 .../H2O_GAPW_XC_2.inp | 0 .../H2O_GAPW_XC_3.inp | 0 .../Ne_GAPW_triplet.inp | 0 tests/QS/regtest-tddfpt-gapw/TEST_FILES | 16 ++ tests/QS/regtest-tddfpt-soc/TEST_FILES | 14 +- tests/QS/regtest-tddfpt/TEST_FILES | 9 - tests/TEST_DIRS | 3 +- 37 files changed, 648 insertions(+), 386 deletions(-) rename tests/QS/regtest-tddfpt-force-gapw/{ch2o_f01_pbe.inp => h2o_f01_pbe.inp} (82%) rename tests/QS/regtest-tddfpt-force-gapw/{ch2o_f01_pbe_gapwxc.inp => h2o_f01_pbe_gapwxc.inp} (100%) rename tests/QS/regtest-tddfpt-force-gapw/{ch2o_f02_coulomb_only.inp => h2o_f02_coulomb_only.inp} (100%) create mode 100644 tests/QS/regtest-tddfpt-force-gapw/h2o_t09.inp create mode 100644 tests/QS/regtest-tddfpt-force-gapw/h2o_t10.inp create mode 100644 tests/QS/regtest-tddfpt-force-gapw/h2o_t11.inp create mode 100644 tests/QS/regtest-tddfpt-force-gapw/h2o_t12.inp rename tests/QS/{regtest-tddfpt => regtest-tddfpt-gapw}/H2O_GAPW_1.inp (100%) rename tests/QS/{regtest-tddfpt => regtest-tddfpt-gapw}/H2O_GAPW_1_triplet.inp (66%) rename tests/QS/{regtest-tddfpt => regtest-tddfpt-gapw}/H2O_GAPW_2.inp (100%) rename tests/QS/{regtest-tddfpt => regtest-tddfpt-gapw}/H2O_GAPW_3.inp (100%) rename tests/QS/{regtest-tddfpt => regtest-tddfpt-gapw}/H2O_GAPW_4.inp (100%) rename tests/QS/{regtest-tddfpt => regtest-tddfpt-gapw}/H2O_GAPW_XC_1.inp (100%) rename tests/QS/{regtest-tddfpt => regtest-tddfpt-gapw}/H2O_GAPW_XC_2.inp (100%) rename tests/QS/{regtest-tddfpt => regtest-tddfpt-gapw}/H2O_GAPW_XC_3.inp (100%) rename tests/QS/{regtest-tddfpt => regtest-tddfpt-gapw}/Ne_GAPW_triplet.inp (100%) create mode 100644 tests/QS/regtest-tddfpt-gapw/TEST_FILES diff --git a/src/qs_environment.F b/src/qs_environment.F index bfe4fe4cc9..a82129eb77 100644 --- a/src/qs_environment.F +++ b/src/qs_environment.F @@ -1146,7 +1146,7 @@ CONTAINS IF (dft_control%qs_control%method_id /= do_method_gapw_xc) THEN CALL get_qs_env(qs_env=qs_env, local_rho_set=local_rho_set, natom=natom) ! *** Allocate and initialize the compensation density rho0 *** - CALL init_rho0(local_rho_set, qs_env, gapw_control, .FALSE.) + CALL init_rho0(local_rho_set, qs_env, gapw_control) ! *** Allocate and Initialize the local coulomb term *** CALL init_coulomb_local(qs_env%hartree_local, natom) END IF diff --git a/src/qs_fxc.F b/src/qs_fxc.F index 50da18bf38..870856f4e8 100644 --- a/src/qs_fxc.F +++ b/src/qs_fxc.F @@ -203,7 +203,9 @@ CONTAINS NULLIFY (vxc00, v_tau_rspace) IF (is_triplet) THEN CPASSERT(nspins == 1) + ! rhoin = (0.5 rho0, 0.5 rho0) CALL qs_rho_copy(rho0_struct, rhoin, auxbas_pw_pool, 2) + ! rhoin = (0.5 rho0 + 0.5 rho1, 0.5 rho0) CALL qs_rho_scale_and_add(rhoin, rho1_struct, alpha, 0.5_dp*beta) CALL qs_vxc_create(ks_env=ks_env, rho_struct=rhoin, xc_section=xc_section, & vxc_rho=vxc00, vxc_tau=v_tau_rspace, exc=exc, just_energy=.FALSE.) diff --git a/src/qs_ks_atom.F b/src/qs_ks_atom.F index 8132ff7fbb..f351f5f614 100644 --- a/src/qs_ks_atom.F +++ b/src/qs_ks_atom.F @@ -94,6 +94,8 @@ CONTAINS !> \param oce_external ... !> \param sab_external ... !> \param kscale ... +!> \param kintegral ... +!> \param kforce ... !> \param fscale ... !> \par History !> created [MI] @@ -103,7 +105,8 @@ CONTAINS !> Allow for external kind_set, rho_atom_set, oce, sab 12.2019 (A. Bussy) ! ************************************************************************************************** SUBROUTINE update_ks_atom(qs_env, ksmat, pmat, forces, tddft, rho_atom_external, & - kind_set_external, oce_external, sab_external, kscale, fscale) + kind_set_external, oce_external, sab_external, kscale, & + kintegral, kforce, fscale) TYPE(qs_environment_type), POINTER :: qs_env TYPE(dbcsr_p_type), DIMENSION(*), INTENT(INOUT) :: ksmat, pmat @@ -116,7 +119,8 @@ CONTAINS TYPE(oce_matrix_type), OPTIONAL, POINTER :: oce_external TYPE(neighbor_list_set_p_type), DIMENSION(:), & OPTIONAL, POINTER :: sab_external - REAL(KIND=dp), INTENT(IN), OPTIONAL :: kscale, fscale(2) + REAL(KIND=dp), INTENT(IN), OPTIONAL :: kscale, kintegral, kforce + REAL(KIND=dp), DIMENSION(2), INTENT(IN), OPTIONAL :: fscale CHARACTER(len=*), PARAMETER :: routineN = 'update_ks_atom' @@ -134,10 +138,11 @@ CONTAINS REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: a_matrix, p_matrix REAL(dp), DIMENSION(3) :: rac, rbc REAL(dp), DIMENSION(3, 3) :: force_tmp - REAL(kind=dp) :: eps_cpc, factor1, factor2, force_fac(2) + REAL(kind=dp) :: eps_cpc, factor1, factor2 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: C_int_h, C_int_s, coc REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dCPC_h, dCPC_s REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: PC_h, PC_s + REAL(KIND=dp), DIMENSION(2) :: force_fac REAL(KIND=dp), DIMENSION(3, 3) :: pv_virial_thread REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: C_coeff_hh_a, C_coeff_hh_b, & C_coeff_ss_a, C_coeff_ss_b @@ -185,11 +190,12 @@ CONTAINS nspins = dft_control%nspins nimages = dft_control%nimages + factor1 = 1.0_dp + factor2 = 1.0_dp + !deal with externals my_tddft = .FALSE. IF (PRESENT(tddft)) my_tddft = tddft - factor1 = 1.0_dp - factor2 = 1.0_dp IF (my_tddft) THEN IF (nspins == 1) factor1 = 2.0_dp CPASSERT(nimages == 1) @@ -198,7 +204,8 @@ CONTAINS factor1 = factor1*kscale factor2 = factor2*kscale END IF - + IF (PRESENT(kintegral)) factor1 = kintegral + IF (PRESENT(kforce)) factor2 = kforce force_fac = 1.0_dp IF (PRESENT(fscale)) force_fac(:) = fscale(:) @@ -376,6 +383,7 @@ CONTAINS DO ispin = 1, nspins NULLIFY (mat_h(ispin)%array, mat_p(ispin)%array) + found = .FALSE. IF (iatom <= jatom) THEN CALL dbcsr_get_block_p(matrix=ksmat(nspins*(img - 1) + ispin)%matrix, & row=iatom, col=jatom, & @@ -385,8 +393,10 @@ CONTAINS row=jatom, col=iatom, & BLOCK=mat_h(ispin)%array, found=found) END IF + CPASSERT(found) IF (forces) THEN + found = .FALSE. IF (iatom <= jatom) THEN CALL dbcsr_get_block_p(matrix=pmat(nspins*(img - 1) + ispin)%matrix, & row=iatom, col=jatom, & @@ -396,6 +406,7 @@ CONTAINS row=jatom, col=iatom, & BLOCK=mat_p(ispin)%array, found=found) END IF + CPASSERT(found) END IF END DO diff --git a/src/qs_ks_reference.F b/src/qs_ks_reference.F index fd74cc141a..24434ff38c 100644 --- a/src/qs_ks_reference.F +++ b/src/qs_ks_reference.F @@ -19,8 +19,7 @@ MODULE qs_ks_reference USE dbcsr_api, ONLY: dbcsr_p_type USE hartree_local_methods, ONLY: Vh_1c_gg_integrals,& init_coulomb_local - USE hartree_local_types, ONLY: ecoul_1center_type,& - hartree_local_create,& + USE hartree_local_types, ONLY: hartree_local_create,& hartree_local_release,& hartree_local_type USE input_constants, ONLY: do_admm_aux_exch_func_none @@ -58,7 +57,8 @@ MODULE qs_ks_reference local_rho_type USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type USE qs_oce_types, ONLY: oce_matrix_type - USE qs_rho0_ggrid, ONLY: rho0_s_grid_create + USE qs_rho0_ggrid, ONLY: integrate_vhg0_rspace,& + rho0_s_grid_create USE qs_rho0_methods, ONLY: init_rho0 USE qs_rho_atom_methods, ONLY: allocate_rho_atom_internals,& calculate_rho_atom_coeff @@ -301,13 +301,15 @@ CONTAINS !> \param qs_env ... !> \param local_rho_set ... !> \param local_rho_set_admm ... +!> \param v_hartree_rspace ... !> \par History !> 07.2022 created [JGH] !> \author JGH ! ************************************************************************************************** - SUBROUTINE ks_ref_potential_atom(qs_env, local_rho_set, local_rho_set_admm) + SUBROUTINE ks_ref_potential_atom(qs_env, local_rho_set, local_rho_set_admm, v_hartree_rspace) TYPE(qs_environment_type), POINTER :: qs_env TYPE(local_rho_type), POINTER :: local_rho_set, local_rho_set_admm + TYPE(pw_type), INTENT(IN) :: v_hartree_rspace CHARACTER(LEN=*), PARAMETER :: routineN = 'ks_ref_potential_atom' @@ -318,7 +320,6 @@ CONTAINS TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_aux, rho_ao_kp TYPE(dft_control_type), POINTER :: dft_control - TYPE(ecoul_1center_type), DIMENSION(:), POINTER :: ecoul_1c TYPE(hartree_local_type), POINTER :: hartree_local TYPE(mp_para_env_type), POINTER :: para_env TYPE(neighbor_list_set_p_type), DIMENSION(:), & @@ -349,7 +350,7 @@ CONTAINS qs_kind_set, dft_control, para_env) IF (gapw) THEN CALL get_qs_env(qs_env, natom=natom) - CALL init_rho0(local_rho_set, qs_env, dft_control%qs_control%gapw_control, .FALSE.) + CALL init_rho0(local_rho_set, qs_env, dft_control%qs_control%gapw_control) CALL rho0_s_grid_create(pw_env, local_rho_set%rho0_mpole) CALL hartree_local_create(hartree_local) CALL init_coulomb_local(hartree_local, natom) @@ -362,9 +363,9 @@ CONTAINS CALL prepare_gapw_den(qs_env, local_rho_set, do_rho0=gapw) IF (gapw) THEN - CALL get_qs_env(qs_env, ecoul_1c=ecoul_1c) - CALL Vh_1c_gg_integrals(qs_env, eh1c, ecoul_1c, local_rho_set, para_env, tddft=.FALSE., & - core_2nd=.FALSE.) + CALL Vh_1c_gg_integrals(qs_env, eh1c, hartree_local%ecoul_1c, local_rho_set, para_env, .FALSE.) + CALL integrate_vhg0_rspace(qs_env, v_hartree_rspace, para_env, calculate_forces=.FALSE., & + local_rho_set=local_rho_set) END IF IF (dft_control%do_admm) THEN CALL get_qs_env(qs_env, admm_env=admm_env) @@ -386,7 +387,7 @@ CONTAINS CALL qs_rho_get(rho, rho_ao_kp=rho_ao_aux) CALL calculate_rho_atom_coeff(qs_env, rho_ao_aux, local_rho_set_admm%rho_atom_set, & admm_env%admm_gapw_env%admm_kind_set, oce, sab, para_env) - CALL prepare_gapw_den(qs_env, local_rho_set=admm_env%admm_gapw_env%local_rho_set, & + CALL prepare_gapw_den(qs_env, local_rho_set=local_rho_set_admm, & do_rho0=.FALSE., kind_set_external=admm_env%admm_gapw_env%admm_kind_set) !compute the potential due to atomic densities xc_section => admm_env%xc_section_aux diff --git a/src/qs_p_env_methods.F b/src/qs_p_env_methods.F index c35355e50d..376b099efc 100644 --- a/src/qs_p_env_methods.F +++ b/src/qs_p_env_methods.F @@ -262,7 +262,8 @@ CONTAINS CALL allocate_rho_atom_internals(p_env%local_rho_set%rho_atom_set, atomic_kind_set, & qs_kind_set, dft_control, para_env) - CALL init_rho0(p_env%local_rho_set, qs_env, dft_control%qs_control%gapw_control, .TRUE.) + CALL init_rho0(p_env%local_rho_set, qs_env, dft_control%qs_control%gapw_control, & + zcore=0.0_dp) CALL rho0_s_grid_create(pw_env, p_env%local_rho_set%rho0_mpole) CALL hartree_local_create(p_env%hartree_local) CALL init_coulomb_local(p_env%hartree_local, natom) diff --git a/src/qs_rho0_ggrid.F b/src/qs_rho0_ggrid.F index 2b705604ce..5b2a7afd71 100644 --- a/src/qs_rho0_ggrid.F +++ b/src/qs_rho0_ggrid.F @@ -319,9 +319,10 @@ CONTAINS !> \param local_rho_set ... !> \param local_rho_set_2nd ... !> \param atener ... +!> \param kforce ... ! ************************************************************************************************** SUBROUTINE integrate_vhg0_rspace(qs_env, v_rspace, para_env, calculate_forces, local_rho_set, & - local_rho_set_2nd, atener) + local_rho_set_2nd, atener, kforce) TYPE(qs_environment_type), POINTER :: qs_env TYPE(pw_type), INTENT(IN) :: v_rspace @@ -329,6 +330,7 @@ CONTAINS LOGICAL, INTENT(IN) :: calculate_forces TYPE(local_rho_type), OPTIONAL, POINTER :: local_rho_set, local_rho_set_2nd REAL(KIND=dp), DIMENSION(:), OPTIONAL :: atener + REAL(KIND=dp), INTENT(IN), OPTIONAL :: kforce CHARACTER(LEN=*), PARAMETER :: routineN = 'integrate_vhg0_rspace' @@ -340,8 +342,8 @@ CONTAINS INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cg_list INTEGER, DIMENSION(:), POINTER :: atom_list, lmax, lmin, npgf LOGICAL :: grid_distributed, paw_atom, use_virial - REAL(KIND=dp) :: eps_rho_rspace, force_tmp(3), ra(3), & - rpgf0, scale, zet0 + REAL(KIND=dp) :: eps_rho_rspace, force_tmp(3), fscale, & + ra(3), rpgf0, zet0 REAL(KIND=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b REAL(KIND=dp), DIMENSION(:), POINTER :: hab_sph, norm_l, Qlm REAL(KIND=dp), DIMENSION(:, :), POINTER :: hab, hdab_sph, intloc, pab @@ -455,8 +457,8 @@ CONTAINS CALL pw_pool_give_back_pw(pw_aux, coeff_gaux) CALL pw_pool_create_pw(pw_aux, coeff_raux, use_data=REALDATA3D, & in_space=REALSPACE) - scale = coeff_rspace%pw_grid%dvol/coeff_raux%pw_grid%dvol - coeff_rspace%cr3d = scale*coeff_rspace%cr3d + fscale = coeff_rspace%pw_grid%dvol/coeff_raux%pw_grid%dvol + coeff_rspace%cr3d = fscale*coeff_rspace%cr3d CALL pw_pool_give_back_pw(pw_aux, coeff_raux) ELSE @@ -496,6 +498,11 @@ CONTAINS grid_distributed = rs_v%desc%distributed + fscale = 1.0_dp + IF (PRESENT(kforce)) THEN + fscale = kforce + END IF + DO ikind = 1, SIZE(atomic_kind_set, 1) NULLIFY (basis_1c_set, atom_list, harmonics) CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=nat) @@ -655,7 +662,7 @@ CONTAINS force_tmp(2) = force_tmp(2) + Qlm(iso)*hdab_sph(2, iso) force_tmp(3) = force_tmp(3) + Qlm(iso)*hdab_sph(3, iso) END DO - force(ikind)%g0s_Vh_elec(1:3, iat) = force(ikind)%g0s_Vh_elec(1:3, iat) + force_tmp(1:3) + force(ikind)%g0s_Vh_elec(1:3, iat) = force(ikind)%g0s_Vh_elec(1:3, iat) + fscale*force_tmp(1:3) END IF IF (use_virial) THEN my_virial_a = 0.0_dp @@ -663,8 +670,8 @@ CONTAINS DO ii = 1, 3 DO i = 1, 3 ! Q from local_rho_set - virial%pv_gapw(i, ii) = virial%pv_gapw(i, ii) + Qlm(iso)*a_hdab_sph(i, ii, iso) - virial%pv_virial(i, ii) = virial%pv_virial(i, ii) + Qlm(iso)*a_hdab_sph(i, ii, iso) + virial%pv_gapw(i, ii) = virial%pv_gapw(i, ii) + fscale*Qlm(iso)*a_hdab_sph(i, ii, iso) + virial%pv_virial(i, ii) = virial%pv_virial(i, ii) + fscale*Qlm(iso)*a_hdab_sph(i, ii, iso) END DO END DO END DO diff --git a/src/qs_rho0_methods.F b/src/qs_rho0_methods.F index e3c2b6832c..4dc94f6d84 100644 --- a/src/qs_rho0_methods.F +++ b/src/qs_rho0_methods.F @@ -342,14 +342,14 @@ CONTAINS !> \param local_rho_set ... !> \param qs_env ... !> \param gapw_control ... -!> \param tddft ... +!> \param zcore ... ! ************************************************************************************************** - SUBROUTINE init_rho0(local_rho_set, qs_env, gapw_control, tddft) + SUBROUTINE init_rho0(local_rho_set, qs_env, gapw_control, zcore) TYPE(local_rho_type), POINTER :: local_rho_set TYPE(qs_environment_type), POINTER :: qs_env TYPE(gapw_control_type), POINTER :: gapw_control - LOGICAL, INTENT(in) :: tddft + REAL(KIND=dp), INTENT(IN), OPTIONAL :: zcore CHARACTER(len=*), PARAMETER :: routineN = 'init_rho0' @@ -430,7 +430,7 @@ CONTAINS basis_set=basis_1c, basis_type="GAPW_1C") ! Set charge distribution of ionic cores to zero when computing the response-density - IF (tddft) zeff = 0.0_dp + IF (PRESENT(zcore)) zeff = zcore CALL get_gto_basis_set(gto_basis_set=basis_1c, & maxl=maxl, & diff --git a/src/qs_tddfpt2_fhxc.F b/src/qs_tddfpt2_fhxc.F index f626245dc8..b6a5144b97 100644 --- a/src/qs_tddfpt2_fhxc.F +++ b/src/qs_tddfpt2_fhxc.F @@ -200,7 +200,8 @@ CONTAINS pw_env_external=sub_env%pw_env, & task_list_external=sub_env%task_list_orb_soft, & para_env_external=sub_env%para_env) - CALL prepare_gapw_den(qs_env, work_matrices%local_rho_set) + CALL prepare_gapw_den(qs_env, work_matrices%local_rho_set, & + do_rho0=(.NOT. is_rks_triplets)) ELSEIF (gapw_xc) THEN CALL qs_rho_update_rho(work_matrices%rho_orb_struct_sub, qs_env, & rho_xc_external=work_matrices%rho_xc_struct_sub, & @@ -408,9 +409,6 @@ CONTAINS qs_env=qs_env, calculate_forces=.FALSE., gapw=gapw, & pw_env_external=sub_env%pw_env, & task_list_external=sub_env%task_list_orb_soft) - ! rho_ia_ao will not be touched - CALL update_ks_atom(qs_env, work_matrices%A_ia_munu_sub, rho_ia_ao, forces=.FALSE., tddft=.TRUE., & - rho_atom_external=work_matrices%local_rho_set%rho_atom_set) ELSEIF (gapw_xc) THEN IF (.NOT. is_rks_triplets) THEN CALL integrate_v_rspace(v_rspace=work_matrices%A_ia_rspace_sub(ispin), & @@ -418,9 +416,6 @@ CONTAINS qs_env=qs_env, calculate_forces=.FALSE., gapw=.FALSE., & pw_env_external=sub_env%pw_env, task_list_external=sub_env%task_list_orb) END IF - ! rho_ia_ao will not be touched - CALL update_ks_atom(qs_env, work_matrices%A_ia_munu_sub, rho_ia_ao, forces=.FALSE., tddft=.TRUE., & - rho_atom_external=work_matrices%local_rho_set%rho_atom_set) ELSE CALL integrate_v_rspace(v_rspace=work_matrices%A_ia_rspace_sub(ispin), & hmat=work_matrices%A_ia_munu_sub(ispin), & @@ -443,6 +438,16 @@ CONTAINS END IF ! for full kernel using lri END DO + ! local atom contributions + IF (.NOT. do_lrigpw) THEN + IF (gapw .OR. gapw_xc) THEN + ! rho_ia_ao will not be touched + CALL update_ks_atom(qs_env, work_matrices%A_ia_munu_sub, rho_ia_ao, forces=.FALSE., & + rho_atom_external=work_matrices%local_rho_set%rho_atom_set, & + tddft=.TRUE.) + END IF + END IF + ! calculate Coulomb contribution to response vector for lrigpw ! ! this is restricting lri to Coulomb only at the moment ! IF (do_lrigpw .AND. (.NOT. is_rks_triplets)) THEN ! diff --git a/src/qs_tddfpt2_fhxc_forces.F b/src/qs_tddfpt2_fhxc_forces.F index 2a981ca15d..86dedae617 100644 --- a/src/qs_tddfpt2_fhxc_forces.F +++ b/src/qs_tddfpt2_fhxc_forces.F @@ -174,7 +174,8 @@ CONTAINS norb, nspins LOGICAL :: distribute_fock_matrix, do_admm, do_analytic, do_hfx, do_numeric, gapw, gapw_xc, & hfx_treat_lsd_in_core, is_rks_triplets, s_mstruct_changed, use_virial - REAL(KIND=dp) :: eh1, eh1c, eps_fit, focc, fval, xehartree + REAL(KIND=dp) :: eh1, eh1c, eps_fit, focc, fval, kval, & + xehartree REAL(KIND=dp), DIMENSION(3) :: fodeb TYPE(admm_type), POINTER :: admm_env TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set @@ -264,6 +265,11 @@ CONTAINS ! NULLIFY (hartree_local, local_rho_set, local_rho_set_admm) IF (gapw .OR. gapw_xc) THEN + IF (nspins == 2) THEN + DO ispin = 1, nspins + CALL dbcsr_scale(matrix_px1(ispin)%matrix, 2.0_dp) + END DO + END IF CALL get_qs_env(qs_env, & atomic_kind_set=atomic_kind_set, & qs_kind_set=qs_kind_set) @@ -272,7 +278,8 @@ CONTAINS qs_kind_set, dft_control, para_env) IF (gapw) THEN CALL get_qs_env(qs_env, natom=natom) - CALL init_rho0(local_rho_set, qs_env, dft_control%qs_control%gapw_control, .TRUE.) + CALL init_rho0(local_rho_set, qs_env, dft_control%qs_control%gapw_control, & + zcore=0.0_dp) CALL rho0_s_grid_create(pw_env, local_rho_set%rho0_mpole) CALL hartree_local_create(hartree_local) CALL init_coulomb_local(hartree_local, natom) @@ -305,6 +312,11 @@ CONTAINS CALL calculate_rho_atom_coeff(qs_env, mpga, local_rho_set_g%rho_atom_set, & qs_kind_set, oce, sab, para_env) CALL prepare_gapw_den(qs_env, local_rho_set_g, do_rho0=.FALSE.) + IF (nspins == 2) THEN + DO ispin = 1, nspins + CALL dbcsr_scale(matrix_px1(ispin)%matrix, 0.5_dp) + END DO + END IF END IF ! IF (do_admm) THEN @@ -426,16 +438,21 @@ CONTAINS IF (debug_forces) fodeb(1:3) = force(1)%g0s_Vh_elec(1:3, 1) CALL Vh_1c_gg_integrals(qs_env, eh1c, hartree_local%ecoul_1c, local_rho_set, para_env, tddft=.TRUE., & core_2nd=.TRUE.) + IF (nspins == 1) THEN + kval = 1.0_dp + ELSE + kval = 0.5_dp + END IF CALL integrate_vhg0_rspace(qs_env, xv_hartree_rspace, para_env, calculate_forces=.TRUE., & - local_rho_set=local_rho_set) + local_rho_set=local_rho_set, kforce=kval) IF (debug_forces) THEN fodeb(1:3) = force(1)%g0s_Vh_elec(1:3, 1) - fodeb(1:3) CALL para_env%sum(fodeb) IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Px*dKh[X]PAWg0", fodeb END IF IF (debug_forces) fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1) - CALL update_ks_atom(qs_env, matrix_hx, matrix_px1, forces=.TRUE., tddft=.TRUE., & - rho_atom_external=local_rho_set%rho_atom_set, kscale=0.5_dp) + CALL update_ks_atom(qs_env, matrix_hx, matrix_px1, forces=.TRUE., & + rho_atom_external=local_rho_set%rho_atom_set) IF (debug_forces) THEN fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1) - fodeb(1:3) CALL para_env%sum(fodeb) @@ -552,16 +569,23 @@ CONTAINS IF (gapw .OR. gapw_xc) THEN IF (debug_forces) fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1) CALL update_ks_atom(qs_env, matrix_fx, matrix_px1, forces=.TRUE., tddft=.TRUE., & - rho_atom_external=local_rho_set_f%rho_atom_set, kscale=0.5_dp) + rho_atom_external=local_rho_set_f%rho_atom_set, & + kintegral=1.0_dp, kforce=0.5_dp) IF (debug_forces) THEN fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1) - fodeb(1:3) CALL para_env%sum(fodeb) IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Px*dKf[X]PAW ", fodeb END IF IF (debug_forces) fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1) - CALL update_ks_atom(qs_env, matrix_gx, matrix_p, forces=.TRUE., tddft=.TRUE., & - rho_atom_external=local_rho_set_g%rho_atom_set, & - kscale=0.5_dp) + IF (nspins == 1) THEN + CALL update_ks_atom(qs_env, matrix_gx, matrix_p, forces=.TRUE., tddft=.TRUE., & + rho_atom_external=local_rho_set_g%rho_atom_set, & + kscale=0.5_dp) + ELSE + CALL update_ks_atom(qs_env, matrix_gx, matrix_p, forces=.TRUE., & + rho_atom_external=local_rho_set_g%rho_atom_set, & + kintegral=0.5_dp, kforce=0.25_dp) + END IF IF (debug_forces) THEN fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1) - fodeb(1:3) CALL para_env%sum(fodeb) diff --git a/src/qs_tddfpt2_forces.F b/src/qs_tddfpt2_forces.F index 7cc0fe05ce..e5f933ea2d 100644 --- a/src/qs_tddfpt2_forces.F +++ b/src/qs_tddfpt2_forces.F @@ -341,7 +341,8 @@ CONTAINS CALL exstate_potential_release(ex_env) CALL ks_ref_potential(qs_env, ex_env%vh_rspace, ex_env%vxc_rspace, & ex_env%vtau_rspace, ex_env%vadmm_rspace, ehartree, exc) - CALL ks_ref_potential_atom(qs_env, ex_env%local_rho_set, ex_env%local_rho_set_admm) + CALL ks_ref_potential_atom(qs_env, ex_env%local_rho_set, ex_env%local_rho_set_admm, & + ex_env%vh_rspace) CALL tddfpt_force_direct(qs_env, ex_env, gs_mos, kernel_env, sub_env, & work_matrices, debug_forces) END IF @@ -739,7 +740,8 @@ CONTAINS CALL local_rho_set_create(local_rho_set) CALL allocate_rho_atom_internals(local_rho_set%rho_atom_set, atomic_kind_set, & qs_kind_set, dft_control, para_env) - CALL init_rho0(local_rho_set, qs_env, dft_control%qs_control%gapw_control, .TRUE.) + CALL init_rho0(local_rho_set, qs_env, dft_control%qs_control%gapw_control, & + zcore=0.0_dp) CALL rho0_s_grid_create(pw_env, local_rho_set%rho0_mpole) CALL hartree_local_create(hartree_local) CALL init_coulomb_local(hartree_local, natom) @@ -805,8 +807,7 @@ CONTAINS CALL pw_scale(v_hartree_rspace, v_hartree_rspace%pw_grid%dvol) IF (gapw) THEN CALL Vh_1c_gg_integrals(qs_env, thartree, hartree_local%ecoul_1c, & - local_rho_set, & - para_env, tddft=.TRUE.) + local_rho_set, para_env, tddft=.TRUE.) CALL integrate_vhg0_rspace(qs_env, v_hartree_rspace, para_env, & calculate_forces=.FALSE., & local_rho_set=local_rho_set) @@ -876,8 +877,8 @@ CONTAINS IF (gapw .OR. gapw_xc) THEN mhz(1:nspins, 1:1) => matrix_hz(1:nspins) mpe(1:nspins, 1:1) => matrix_pe(1:nspins) - CALL update_ks_atom(qs_env, mhz, mpe, forces=.FALSE., tddft=.TRUE., & - rho_atom_external=local_rho_set%rho_atom_set, kscale=0.5_dp) + CALL update_ks_atom(qs_env, mhz, mpe, forces=.FALSE., & + rho_atom_external=local_rho_set%rho_atom_set) END IF CALL pw_pool_give_back_pw(auxbas_pw_pool, v_hartree_gspace) diff --git a/src/qs_tddfpt2_subgroups.F b/src/qs_tddfpt2_subgroups.F index 723627fe71..8c00d93f38 100644 --- a/src/qs_tddfpt2_subgroups.F +++ b/src/qs_tddfpt2_subgroups.F @@ -346,7 +346,8 @@ CONTAINS CALL allocate_rho_atom_internals(sub_env%local_rho_set%rho_atom_set, atomic_kind_set, & qs_kind_set, dft_control, sub_env%para_env) - CALL init_rho0(sub_env%local_rho_set, qs_env, dft_control%qs_control%gapw_control, .TRUE.) + CALL init_rho0(sub_env%local_rho_set, qs_env, dft_control%qs_control%gapw_control, & + zcore=0.0_dp) CALL rho0_s_grid_create(sub_env%pw_env, sub_env%local_rho_set%rho0_mpole) CALL hartree_local_create(sub_env%hartree_local) CALL init_coulomb_local(sub_env%hartree_local, natom) diff --git a/src/qs_tddfpt2_types.F b/src/qs_tddfpt2_types.F index 3e5dfb26c0..c48b45ed1f 100644 --- a/src/qs_tddfpt2_types.F +++ b/src/qs_tddfpt2_types.F @@ -446,7 +446,8 @@ CONTAINS CALL local_rho_set_create(work_matrices%local_rho_set) CALL allocate_rho_atom_internals(work_matrices%local_rho_set%rho_atom_set, atomic_kind_set, & qs_kind_set, dft_control, sub_env%para_env) - CALL init_rho0(work_matrices%local_rho_set, qs_env, dft_control%qs_control%gapw_control, .TRUE.) + CALL init_rho0(work_matrices%local_rho_set, qs_env, dft_control%qs_control%gapw_control, & + zcore=0.0_dp) CALL rho0_s_grid_create(sub_env%pw_env, work_matrices%local_rho_set%rho0_mpole) CALL hartree_local_create(work_matrices%hartree_local) CALL init_coulomb_local(work_matrices%hartree_local, natom) diff --git a/src/qs_vxc_atom.F b/src/qs_vxc_atom.F index 81b4eac990..c22656dfd2 100644 --- a/src/qs_vxc_atom.F +++ b/src/qs_vxc_atom.F @@ -377,7 +377,6 @@ CONTAINS IF (.NOT. energy_only) THEN NULLIFY (int_hh, int_ss) - rho_atom => my_rho_atom_set(iatom) CALL get_rho_atom(rho_atom=rho_atom, ga_Vlocal_gb_h=int_hh, ga_Vlocal_gb_s=int_ss) IF (gradient_f) THEN CALL gaVxcgb_GC(vxc_h, vxc_s, vxg_h, vxg_s, int_hh, int_ss, & @@ -460,7 +459,7 @@ CONTAINS INTEGER, DIMENSION(2, 3) :: bounds INTEGER, DIMENSION(:), POINTER :: atom_list LOGICAL :: gradient_functional, lsd, lsd_singlets, & - my_tddft, paw_atom, tau_f + my_tddft, paw_atom, scale_rho, tau_f REAL(KIND=dp) :: density_cut, gradient_cut, rtot, tau_cut REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), & POINTER :: vxc_h, vxc_s @@ -532,6 +531,23 @@ CONTAINS ELSE nspins = 1 END IF + + scale_rho = .FALSE. + IF (my_tddft) THEN + IF (excitations == tddfpt_excitations) THEN + IF (nspins == 1 .AND. (lsd_singlets .OR. res_etype == tddfpt_triplet)) THEN + lsd = .TRUE. + END IF + END IF + ELSEIF (PRESENT(do_tddfpt2) .AND. PRESENT(do_triplet)) THEN + IF (nspins == 1 .AND. do_triplet) THEN + lsd = .TRUE. + scale_rho = .TRUE. + END IF + ELSEIF (PRESENT(do_triplet)) THEN + IF (nspins == 1 .AND. do_triplet) lsd = .TRUE. + END IF + needs = xc_functionals_get_needs(xc_fun_section, lsd=lsd, & calc_potential=.TRUE.) gradient_functional = needs%drho .OR. needs%drho_spin @@ -542,18 +558,6 @@ CONTAINS rtau = 0.0_dp END IF - IF (my_tddft) THEN - IF (excitations == tddfpt_excitations) THEN - IF (nspins == 1 .AND. (lsd_singlets .OR. res_etype == tddfpt_triplet)) THEN - lsd = .TRUE. - END IF - END IF - ELSEIF (PRESENT(do_tddfpt2) .AND. PRESENT(do_triplet)) THEN - IF (nspins == 1 .AND. do_triplet) lsd = .TRUE. - ELSEIF (PRESENT(do_triplet)) THEN - IF (nspins == 1 .AND. do_triplet) lsd = .TRUE. - END IF - ! Here starts the loop over all the atoms DO ikind = 1, SIZE(atomic_kind_set) @@ -666,6 +670,14 @@ CONTAINS ir, r1_h, r1_s, rho1_h, rho1_s, dr1_h, dr1_s, r1_h_d, r1_s_d, & drho1_h, drho1_s) END DO + IF (scale_rho) THEN + rho_h = 2.0_dp*rho_h + rho_s = 2.0_dp*rho_s + IF (gradient_functional) THEN + drho_h = 2.0_dp*drho_h + drho_s = 2.0_dp*drho_s + END IF + END IF DO ir = 1, nr IF (tau_f) THEN @@ -752,8 +764,8 @@ CONTAINS REAL(KIND=dp), PARAMETER :: epsrho = 5.e-4_dp INTEGER :: bo(2), handle, iat, iatom, ikind, ir, & - istep, myfun, na, natom, nf, nr, ns, & - nspins, nstep, num_pe + istep, mspins, myfun, na, natom, nf, & + nr, ns, nspins, nstep, num_pe, nx INTEGER, DIMENSION(2, 3) :: bounds INTEGER, DIMENSION(:), POINTER :: atom_list LOGICAL :: donlcc, gradient_f, lsd, nlcc, paw_atom, & @@ -826,13 +838,15 @@ CONTAINS nlcc = has_nlcc(kind_set) lsd = dft_control%lsd nspins = dft_control%nspins + mspins = nspins + IF (is_triplet) THEN + CPASSERT(nspins == 1) + lsd = .TRUE. + mspins = 2 + END IF needs = xc_functionals_get_needs(xc_fun_section, lsd=lsd, calc_potential=.TRUE.) gradient_f = (needs%drho .OR. needs%drho_spin) tau_f = (needs%tau .OR. needs%tau_spin) - IF (is_triplet) THEN - CPASSERT(nspins == 1) - CPABORT("Missing Code") - END IF ! Here starts the loop over all the atoms DO ikind = 1, SIZE(atomic_kind_set) @@ -864,26 +878,26 @@ CONTAINS drho_cutoff=gradient_cut, tau_cutoff=tau_cut) ! allocate the required 3d arrays where to store rho and drho - CALL xc_rho_set_atom_update(rho_set_h, needs, nspins, bounds) - CALL xc_rho_set_atom_update(rho_set_s, needs, nspins, bounds) + CALL xc_rho_set_atom_update(rho_set_h, needs, mspins, bounds) + CALL xc_rho_set_atom_update(rho_set_s, needs, mspins, bounds) weight => grid_atom%weight - ALLOCATE (rho_h(na, nr, nspins), rho_s(na, nr, nspins), & + ALLOCATE (rho_h(na, nr, mspins), rho_s(na, nr, mspins), & rho0_h(na, nr, nspins), rho0_s(na, nr, nspins), & rho1_h(na, nr, nspins), rho1_s(na, nr, nspins)) - ALLOCATE (vxc_h(na, nr, nspins), vxc_s(na, nr, nspins)) + ALLOCATE (vxc_h(na, nr, mspins), vxc_s(na, nr, mspins)) IF (gradient_f) THEN - ALLOCATE (drho_h(4, na, nr, nspins), drho_s(4, na, nr, nspins), & + ALLOCATE (drho_h(4, na, nr, mspins), drho_s(4, na, nr, mspins), & drho0_h(4, na, nr, nspins), drho0_s(4, na, nr, nspins), & drho1_h(4, na, nr, nspins), drho1_s(4, na, nr, nspins)) - ALLOCATE (vxg_h(3, na, nr, nspins), vxg_s(3, na, nr, nspins)) + ALLOCATE (vxg_h(3, na, nr, mspins), vxg_s(3, na, nr, mspins)) END IF IF (tau_f) THEN - ALLOCATE (tau_h(na, nr, nspins), tau_s(na, nr, nspins), & + ALLOCATE (tau_h(na, nr, mspins), tau_s(na, nr, mspins), & tau0_h(na, nr, nspins), tau0_s(na, nr, nspins), & tau1_h(na, nr, nspins), tau1_s(na, nr, nspins)) - ALLOCATE (vtau_h(na, nr, nspins), vtau_s(na, nr, nspins)) + ALLOCATE (vtau_h(na, nr, mspins), vtau_s(na, nr, mspins)) END IF ! ! NLCC: prepare rho and drho of the core charge for this KIND @@ -904,11 +918,12 @@ CONTAINS NULLIFY (int_hh, int_ss) rho0_atom => rho0_atom_set(iatom) CALL get_rho_atom(rho_atom=rho0_atom, ga_Vlocal_gb_h=int_hh, ga_Vlocal_gb_s=int_ss) - ALLOCATE (fint_ss(nspins), fint_hh(nspins)) - DO ns = 1, nspins - nf = SIZE(int_ss(ns)%r_coef, 1) + ALLOCATE (fint_ss(mspins), fint_hh(mspins)) + DO ns = 1, mspins + nx = MIN(nspins, ns) + nf = SIZE(int_ss(nx)%r_coef, 1) ALLOCATE (fint_ss(ns)%r_coef(nf, nf)) - nf = SIZE(int_hh(ns)%r_coef, 1) + nf = SIZE(int_hh(nx)%r_coef, 1) ALLOCATE (fint_hh(ns)%r_coef(nf, nf)) END DO @@ -972,27 +987,52 @@ CONTAINS beta = REAL(istep, KIND=dp)*epsrho - rho_h = rho0_h + beta*rho1_h - rho_s = rho0_s + beta*rho1_s - IF (gradient_f) THEN - drho_h = drho0_h + beta*drho1_h - drho_s = drho0_s + beta*drho1_s - END IF - IF (tau_f) THEN - tau_h = tau0_h + beta*tau1_h - tau_s = tau0_s + beta*tau1_s + IF (is_triplet) THEN + rho_h(:, :, 1) = rho0_h(:, :, 1) + beta*rho1_h(:, :, 1) + rho_h(:, :, 2) = rho0_h(:, :, 1) + rho_h = 0.5_dp*rho_h + rho_s(:, :, 1) = rho0_s(:, :, 1) + beta*rho1_s(:, :, 1) + rho_s(:, :, 2) = rho0_s(:, :, 1) + rho_s = 0.5_dp*rho_s + IF (gradient_f) THEN + drho_h(:, :, :, 1) = drho0_h(:, :, :, 1) + beta*drho1_h(:, :, :, 1) + drho_h(:, :, :, 2) = drho0_h(:, :, :, 1) + drho_h = 0.5_dp*drho_h + drho_s(:, :, :, 1) = drho0_s(:, :, :, 1) + beta*drho1_s(:, :, :, 1) + drho_s(:, :, :, 2) = drho0_s(:, :, :, 1) + drho_s = 0.5_dp*drho_s + END IF + IF (tau_f) THEN + tau_h(:, :, 1) = tau0_h(:, :, 1) + beta*tau1_h(:, :, 1) + tau_h(:, :, 2) = tau0_h(:, :, 1) + tau_h = 0.5_dp*tau0_h + tau_s(:, :, 1) = tau0_s(:, :, 1) + beta*tau1_s(:, :, 1) + tau_s(:, :, 2) = tau0_s(:, :, 1) + tau_s = 0.5_dp*tau0_s + END IF + ELSE + rho_h = rho0_h + beta*rho1_h + rho_s = rho0_s + beta*rho1_s + IF (gradient_f) THEN + drho_h = drho0_h + beta*drho1_h + drho_s = drho0_s + beta*drho1_s + END IF + IF (tau_f) THEN + tau_h = tau0_h + beta*tau1_h + tau_s = tau0_s + beta*tau1_s + END IF END IF DO ir = 1, nr IF (tau_f) THEN - CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, drho_h, tau_h, na, ir) - CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, drho_s, tau_s, na, ir) + CALL fill_rho_set(rho_set_h, lsd, mspins, needs, rho_h, drho_h, tau_h, na, ir) + CALL fill_rho_set(rho_set_s, lsd, mspins, needs, rho_s, drho_s, tau_s, na, ir) ELSE IF (gradient_f) THEN - CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, drho_h, tau_d, na, ir) - CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, drho_s, tau_d, na, ir) + CALL fill_rho_set(rho_set_h, lsd, mspins, needs, rho_h, drho_h, tau_d, na, ir) + CALL fill_rho_set(rho_set_s, lsd, mspins, needs, rho_s, drho_s, tau_d, na, ir) ELSE - CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, rho_d, tau_d, na, ir) - CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, rho_d, tau_d, na, ir) + CALL fill_rho_set(rho_set_h, lsd, mspins, needs, rho_h, rho_d, tau_d, na, ir) + CALL fill_rho_set(rho_set_s, lsd, mspins, needs, rho_s, rho_d, tau_d, na, ir) END IF END DO @@ -1000,25 +1040,43 @@ CONTAINS CALL xc_dset_zero_all(deriv_set) CALL vxc_of_r_new(xc_fun_section, rho_set_h, deriv_set, 1, needs, weight, & lsd, na, nr, exc_h, vxc_h, vxg_h, vtau_h) + IF (is_triplet) THEN + vxc_h(:, :, 1) = vxc_h(:, :, 1) - vxc_h(:, :, 2) + IF (gradient_f) THEN + vxg_h(:, :, :, 1) = vxg_h(:, :, :, 1) - vxg_h(:, :, :, 2) + END IF + IF (tau_f) THEN + vtau_h(:, :, 1) = vtau_h(:, :, 1) - vtau_h(:, :, 2) + END IF + END IF ! soft atom density ! CALL xc_dset_zero_all(deriv_set) CALL vxc_of_r_new(xc_fun_section, rho_set_s, deriv_set, 1, needs, weight, & lsd, na, nr, exc_s, vxc_s, vxg_s, vtau_s) + IF (is_triplet) THEN + vxc_s(:, :, 1) = vxc_s(:, :, 1) - vxc_s(:, :, 2) + IF (gradient_f) THEN + vxg_s(:, :, :, 1) = vxg_s(:, :, :, 1) - vxg_s(:, :, :, 2) + END IF + IF (tau_f) THEN + vtau_s(:, :, 1) = vtau_s(:, :, 1) - vtau_s(:, :, 2) + END IF + END IF ! potentials - DO ns = 1, nspins + DO ns = 1, mspins fint_ss(ns)%r_coef(:, :) = 0.0_dp fint_hh(ns)%r_coef(:, :) = 0.0_dp END DO IF (gradient_f) THEN CALL gaVxcgb_GC(vxc_h, vxc_s, vxg_h, vxg_s, fint_hh, fint_ss, & - grid_atom, basis_1c, harmonics, nspins) + grid_atom, basis_1c, harmonics, mspins) ELSE CALL gaVxcgb_noGC(vxc_h, vxc_s, fint_hh, fint_ss, & - grid_atom, basis_1c, harmonics, nspins) + grid_atom, basis_1c, harmonics, mspins) END IF IF (tau_f) THEN CALL dgaVtaudgb(vtau_h, vtau_s, fint_hh, fint_ss, & - grid_atom, basis_1c, harmonics, nspins) + grid_atom, basis_1c, harmonics, mspins) END IF ! first derivative fxc NULLIFY (int_hh, int_ss) @@ -1036,7 +1094,7 @@ CONTAINS END DO END DO ! - DO ns = 1, nspins + DO ns = 1, mspins DEALLOCATE (fint_ss(ns)%r_coef) DEALLOCATE (fint_hh(ns)%r_coef) END DO diff --git a/src/response_solver.F b/src/response_solver.F index a4d5637504..98e7d400b6 100644 --- a/src/response_solver.F +++ b/src/response_solver.F @@ -84,7 +84,6 @@ MODULE response_solver USE physcon, ONLY: pascal USE pw_env_types, ONLY: pw_env_get,& pw_env_type - USE pw_grid_types, ONLY: pw_grid_type USE pw_methods, ONLY: pw_axpy,& pw_copy,& pw_integral_ab,& @@ -100,13 +99,9 @@ MODULE response_solver REALDATA3D,& REALSPACE,& RECIPROCALSPACE,& - pw_create,& - pw_release,& pw_type USE qs_2nd_kernel_ao, ONLY: build_dm_response USE qs_collocate_density, ONLY: calculate_rho_elec - USE qs_core_energies, ONLY: calculate_ecore_overlap,& - calculate_ecore_self USE qs_density_matrices, ONLY: calculate_whz_matrix,& calculate_wz_matrix USE qs_energy_types, ONLY: qs_energy_type @@ -150,7 +145,6 @@ MODULE response_solver calculate_rho_atom_coeff USE qs_rho_types, ONLY: qs_rho_get,& qs_rho_type - USE qs_vxc, ONLY: qs_vxc_create USE qs_vxc_atom, ONLY: calculate_vxc_atom,& calculate_xc_2nd_deriv_atom USE task_list_types, ONLY: task_list_type @@ -1171,7 +1165,7 @@ CONTAINS CALL local_rho_set_create(local_rho_set_gs) CALL allocate_rho_atom_internals(local_rho_set_gs%rho_atom_set, atomic_kind_set, & qs_kind_set, dft_control, para_env) - CALL init_rho0(local_rho_set_gs, qs_env, dft_control%qs_control%gapw_control, .FALSE.) + CALL init_rho0(local_rho_set_gs, qs_env, dft_control%qs_control%gapw_control) CALL rho0_s_grid_create(pw_env, local_rho_set_gs%rho0_mpole) CALL calculate_rho_atom_coeff(qs_env, matrix_p(:, 1), local_rho_set_gs%rho_atom_set, & qs_kind_set, oce, sab_orb, para_env) @@ -1181,7 +1175,8 @@ CONTAINS CALL local_rho_set_create(local_rho_set_t) CALL allocate_rho_atom_internals(local_rho_set_t%rho_atom_set, atomic_kind_set, & qs_kind_set, dft_control, para_env) - CALL init_rho0(local_rho_set_t, qs_env, dft_control%qs_control%gapw_control, .TRUE.) + CALL init_rho0(local_rho_set_t, qs_env, dft_control%qs_control%gapw_control, & + zcore=0.0_dp) CALL rho0_s_grid_create(pw_env, local_rho_set_t%rho0_mpole) CALL calculate_rho_atom_coeff(qs_env, mpa(:), local_rho_set_t%rho_atom_set, & qs_kind_set, oce, sab_orb, para_env) @@ -1286,8 +1281,8 @@ CONTAINS calculate_forces=.TRUE.) END DO IF (myfun /= xc_none) THEN - CALL pw_zero(vhxc_rspace) DO ispin = 1, nspins + CALL pw_zero(vhxc_rspace) CALL pw_axpy(vxc_rspace(ispin), vhxc_rspace) CALL integrate_v_rspace(v_rspace=vhxc_rspace, & hmat=scrm(ispin), pmat=mpa(ispin), & @@ -1322,11 +1317,9 @@ CONTAINS ! ! HXC term IF (debug_forces) fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1) IF (gapw) CALL update_ks_atom(qs_env, scrm, mpa, forces=.TRUE., tddft=.FALSE., & - rho_atom_external=local_rho_set_gs%rho_atom_set, & - kscale=1.0_dp) + rho_atom_external=local_rho_set_gs%rho_atom_set) IF (myfun /= xc_none) CALL update_ks_atom(qs_env, scrm, mpa, forces=.TRUE., tddft=.FALSE., & - rho_atom_external=local_rho_set_vxc%rho_atom_set, & - kscale=1.0_dp) + rho_atom_external=local_rho_set_vxc%rho_atom_set) IF (debug_forces) THEN fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1) - fodeb(1:3) CALL para_env%sum(fodeb) @@ -1622,7 +1615,8 @@ CONTAINS CALL local_rho_set_create(local_rho_set_t) CALL allocate_rho_atom_internals(local_rho_set_t%rho_atom_set, atomic_kind_set, & qs_kind_set, dft_control, para_env) - CALL init_rho0(local_rho_set_t, qs_env, dft_control%qs_control%gapw_control, .TRUE.) + CALL init_rho0(local_rho_set_t, qs_env, dft_control%qs_control%gapw_control, & + zcore=0.0_dp) CALL rho0_s_grid_create(pw_env, local_rho_set_t%rho0_mpole) CALL calculate_rho_atom_coeff(qs_env, mpa(:), local_rho_set_t%rho_atom_set, & qs_kind_set, oce, sab_orb, para_env) @@ -1631,7 +1625,7 @@ CONTAINS CALL local_rho_set_create(local_rho_set_gs) CALL allocate_rho_atom_internals(local_rho_set_gs%rho_atom_set, atomic_kind_set, & qs_kind_set, dft_control, para_env) - CALL init_rho0(local_rho_set_gs, qs_env, dft_control%qs_control%gapw_control, .FALSE.) + CALL init_rho0(local_rho_set_gs, qs_env, dft_control%qs_control%gapw_control) CALL rho0_s_grid_create(pw_env, local_rho_set_gs%rho0_mpole) CALL calculate_rho_atom_coeff(qs_env, matrix_p(:, 1), local_rho_set_gs%rho_atom_set, & qs_kind_set, oce, sab_orb, para_env) @@ -1778,8 +1772,7 @@ CONTAINS IF (myfun /= xc_none) THEN IF (debug_forces) fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1) CALL update_ks_atom(qs_env, matrix_hz, matrix_p, forces=.TRUE., tddft=.FALSE., & - rho_atom_external=local_rho_set_f%rho_atom_set, & - kscale=1.0_dp) + rho_atom_external=local_rho_set_f%rho_atom_set) IF (debug_forces) THEN fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1) - fodeb(1:3) CALL para_env%sum(fodeb) @@ -1790,8 +1783,7 @@ CONTAINS IF (gapw) THEN IF (debug_forces) fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1) CALL update_ks_atom(qs_env, matrix_ht, matrix_p, forces=.TRUE., tddft=.FALSE., & - rho_atom_external=local_rho_set_t%rho_atom_set, & - kscale=1.0_dp) + rho_atom_external=local_rho_set_t%rho_atom_set) IF (debug_forces) THEN fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1) - fodeb(1:3) CALL para_env%sum(fodeb) @@ -2642,218 +2634,6 @@ CONTAINS END SUBROUTINE response_force_xtb -! ************************************************************************************************** -!> \brief calculate the Kohn-Sham reference potential -!> \param qs_env ... -!> \param vh_rspace ... -!> \param vxc_rspace ... -!> \param vtau_rspace ... -!> \param vadmm_rspace ... -!> \param ehartree ... -!> \param exc ... -!> \param h_stress container for the stress tensor of the Hartree term -!> \par History -!> 10.2019 created [JGH] -!> \author JGH -! ************************************************************************************************** - SUBROUTINE ks_ref_potential(qs_env, vh_rspace, vxc_rspace, vtau_rspace, vadmm_rspace, ehartree, exc, h_stress) - TYPE(qs_environment_type), POINTER :: qs_env - TYPE(pw_type), INTENT(INOUT) :: vh_rspace - TYPE(pw_type), DIMENSION(:), POINTER :: vxc_rspace, vtau_rspace, vadmm_rspace - REAL(KIND=dp), INTENT(OUT) :: ehartree, exc - REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT), & - OPTIONAL :: h_stress - - CHARACTER(LEN=*), PARAMETER :: routineN = 'ks_ref_potential' - - INTEGER :: handle, iab, ispin, nspins - REAL(dp) :: eadmm, eovrl, eself - REAL(KIND=dp), DIMENSION(3, 3) :: virial_xc - TYPE(admm_type), POINTER :: admm_env - TYPE(dft_control_type), POINTER :: dft_control - TYPE(mp_para_env_type), POINTER :: para_env - TYPE(pw_env_type), POINTER :: pw_env - TYPE(pw_grid_type), POINTER :: pw_grid - TYPE(pw_poisson_type), POINTER :: poisson_env - TYPE(pw_pool_type), POINTER :: auxbas_pw_pool - TYPE(pw_type) :: rho_tot_gspace, v_hartree_gspace, & - v_hartree_rspace - TYPE(pw_type), DIMENSION(:), POINTER :: v_admm_rspace, v_admm_tau_rspace, & - v_rspace, v_tau_rspace - TYPE(pw_type), POINTER :: rho_core - TYPE(qs_ks_env_type), POINTER :: ks_env - TYPE(qs_rho_type), POINTER :: rho - TYPE(section_vals_type), POINTER :: xc_section - TYPE(virial_type), POINTER :: virial - - CALL timeset(routineN, handle) - - ! get all information on the electronic density - NULLIFY (rho, ks_env) - CALL get_qs_env(qs_env=qs_env, rho=rho, dft_control=dft_control, & - para_env=para_env, ks_env=ks_env, rho_core=rho_core) - - nspins = dft_control%nspins - - NULLIFY (pw_env) - CALL get_qs_env(qs_env=qs_env, pw_env=pw_env) - CPASSERT(ASSOCIATED(pw_env)) - - NULLIFY (auxbas_pw_pool, poisson_env) - ! gets the tmp grids - CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, & - poisson_env=poisson_env) - - ! Calculate the Hartree potential - CALL pw_pool_create_pw(auxbas_pw_pool, v_hartree_gspace, & - use_data=COMPLEXDATA1D, in_space=RECIPROCALSPACE) - CALL pw_pool_create_pw(auxbas_pw_pool, v_hartree_rspace, & - use_data=REALDATA3D, in_space=REALSPACE) - CALL pw_pool_create_pw(auxbas_pw_pool, rho_tot_gspace, & - use_data=COMPLEXDATA1D, in_space=RECIPROCALSPACE) - - ! Get the total density in g-space [ions + electrons] - CALL calc_rho_tot_gspace(rho_tot_gspace, qs_env, rho) - - CALL pw_poisson_solve(poisson_env, rho_tot_gspace, ehartree, & - v_hartree_gspace, h_stress=h_stress, rho_core=rho_core) - CALL pw_transfer(v_hartree_gspace, v_hartree_rspace) - CALL pw_scale(v_hartree_rspace, v_hartree_rspace%pw_grid%dvol) - - CALL pw_pool_give_back_pw(auxbas_pw_pool, v_hartree_gspace) - CALL pw_pool_give_back_pw(auxbas_pw_pool, rho_tot_gspace) - ! - CALL calculate_ecore_self(qs_env, E_self_core=eself) - CALL calculate_ecore_overlap(qs_env, para_env, PRESENT(h_stress), E_overlap_core=eovrl) - ehartree = ehartree + eovrl + eself - - ! v_rspace and v_tau_rspace are generated from the auxbas pool - IF (dft_control%do_admm) THEN - CALL get_qs_env(qs_env, admm_env=admm_env) - xc_section => admm_env%xc_section_primary - ELSE - xc_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC") - END IF - NULLIFY (v_rspace, v_tau_rspace) - CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho, xc_section=xc_section, & - vxc_rho=v_rspace, vxc_tau=v_tau_rspace, exc=exc, just_energy=.FALSE.) - - NULLIFY (v_admm_rspace, v_admm_tau_rspace) - IF (dft_control%do_admm) THEN - IF (qs_env%admm_env%aux_exch_func /= do_admm_aux_exch_func_none) THEN - ! For the virial, we have to save the pv_xc component because it will be reset in qs_vxc_create - IF (PRESENT(h_stress)) THEN - CALL get_qs_env(qs_env, virial=virial) - virial_xc = virial%pv_xc - END IF - CALL get_admm_env(admm_env, rho_aux_fit=rho) - xc_section => admm_env%xc_section_aux - CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho, xc_section=xc_section, & - vxc_rho=v_admm_rspace, vxc_tau=v_admm_tau_rspace, exc=eadmm, just_energy=.FALSE.) - IF (PRESENT(h_stress)) virial%pv_xc = virial%pv_xc + virial_xc - END IF - END IF - - ! allocate potentials - IF (ASSOCIATED(vh_rspace%pw_grid)) THEN - CALL pw_release(vh_rspace) - END IF - IF (ASSOCIATED(vxc_rspace)) THEN - DO iab = 1, SIZE(vxc_rspace) - CALL pw_release(vxc_rspace(iab)) - END DO - ELSE - ALLOCATE (vxc_rspace(nspins)) - END IF - IF (ASSOCIATED(v_tau_rspace)) THEN - IF (ASSOCIATED(vtau_rspace)) THEN - DO iab = 1, SIZE(vtau_rspace) - CALL pw_release(vtau_rspace(iab)) - END DO - ELSE - ALLOCATE (vtau_rspace(nspins)) - END IF - ELSE - NULLIFY (vtau_rspace) - END IF - IF (ASSOCIATED(v_admm_rspace)) THEN - IF (ASSOCIATED(vadmm_rspace)) THEN - DO iab = 1, SIZE(vadmm_rspace) - CALL pw_release(vadmm_rspace(iab)) - END DO - ELSE - ALLOCATE (vadmm_rspace(nspins)) - END IF - ELSE - NULLIFY (vadmm_rspace) - END IF - - pw_grid => v_hartree_rspace%pw_grid - CALL pw_create(vh_rspace, pw_grid, use_data=REALDATA3D, in_space=REALSPACE) - DO ispin = 1, nspins - CALL pw_create(vxc_rspace(ispin), pw_grid, & - use_data=REALDATA3D, in_space=REALSPACE) - IF (ASSOCIATED(vtau_rspace)) THEN - CALL pw_create(vtau_rspace(ispin), pw_grid, & - use_data=REALDATA3D, in_space=REALSPACE) - END IF - IF (ASSOCIATED(vadmm_rspace)) THEN - CALL pw_create(vadmm_rspace(ispin), pw_grid, & - use_data=REALDATA3D, in_space=REALSPACE) - END IF - END DO - ! - CALL pw_transfer(v_hartree_rspace, vh_rspace) - IF (ASSOCIATED(v_rspace)) THEN - DO ispin = 1, nspins - CALL pw_transfer(v_rspace(ispin), vxc_rspace(ispin)) - CALL pw_scale(vxc_rspace(ispin), v_rspace(ispin)%pw_grid%dvol) - IF (ASSOCIATED(v_tau_rspace)) THEN - CALL pw_transfer(v_tau_rspace(ispin), vtau_rspace(ispin)) - CALL pw_scale(vtau_rspace(ispin), v_tau_rspace(ispin)%pw_grid%dvol) - END IF - END DO - ELSE - DO ispin = 1, nspins - CALL pw_zero(vxc_rspace(ispin)) - END DO - END IF - IF (ASSOCIATED(v_admm_rspace)) THEN - DO ispin = 1, nspins - CALL pw_transfer(v_admm_rspace(ispin), vadmm_rspace(ispin)) - CALL pw_scale(vadmm_rspace(ispin), vadmm_rspace(ispin)%pw_grid%dvol) - END DO - END IF - - ! return pw grids - CALL pw_pool_give_back_pw(auxbas_pw_pool, v_hartree_rspace) - IF (ASSOCIATED(v_rspace)) THEN - DO ispin = 1, nspins - CALL pw_pool_give_back_pw(auxbas_pw_pool, v_rspace(ispin)) - IF (ASSOCIATED(v_tau_rspace)) THEN - CALL pw_pool_give_back_pw(auxbas_pw_pool, v_tau_rspace(ispin)) - END IF - END DO - DEALLOCATE (v_rspace) - END IF - IF (ASSOCIATED(v_tau_rspace)) DEALLOCATE (v_tau_rspace) - IF (ASSOCIATED(v_admm_rspace)) THEN - DO ispin = 1, nspins - CALL pw_pool_give_back_pw(auxbas_pw_pool, v_admm_rspace(ispin)) - END DO - DEALLOCATE (v_admm_rspace) - END IF - IF (ASSOCIATED(v_admm_tau_rspace)) THEN - DO ispin = 1, nspins - CALL pw_pool_give_back_pw(auxbas_pw_pool, v_admm_tau_rspace(ispin)) - END DO - DEALLOCATE (v_admm_tau_rspace) - END IF - - CALL timestop(handle) - - END SUBROUTINE ks_ref_potential - ! ************************************************************************************************** !> \brief Win = focc*(P*(H[P_out - P_in] + H[Z] )*P) !> Langrange multiplier matrix with response and perturbation (Harris) kernel matrices diff --git a/src/xc/xc_atom.F b/src/xc/xc_atom.F index 0fe2dc5d88..467e94d681 100644 --- a/src/xc/xc_atom.F +++ b/src/xc/xc_atom.F @@ -392,10 +392,9 @@ CONTAINS CALL divide_by_norm_drho(deriv_set, rho_set, lsd) - ! multiply by -w + ! multiply by w pos => deriv_set%derivs DO WHILE (cp_sll_xc_deriv_next(pos, el_att=deriv_att)) - !deriv_att%deriv_data(:,:,1) = -w(:,:)*deriv_att%deriv_data(:,:,1) deriv_att%deriv_data(:, :, 1) = w(:, :)*deriv_att%deriv_data(:, :, 1) END DO diff --git a/tests/QS/regtest-tddfpt-force-gapw/TEST_FILES b/tests/QS/regtest-tddfpt-force-gapw/TEST_FILES index d37d77c0f7..f56adac1d8 100644 --- a/tests/QS/regtest-tddfpt-force-gapw/TEST_FILES +++ b/tests/QS/regtest-tddfpt-force-gapw/TEST_FILES @@ -3,21 +3,20 @@ # 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_f01_coulomb_only.inp 1 1.0E-11 -17.13968888100846 -ch2o_f01_pbe_gapwxc.inp 1 1.0E-11 -16.99417766936085 -ch2o_f01_pbe.inp 1 1.0E-11 -16.98005256434762 -ch2o_f02_coulomb_only.inp 1 1.0E-11 -12.74156401757076 h2o_f01_coulomb_only.inp 37 1.0E-4 0.101128E+01 -ch2o_f01_pbe_gapwxc.inp 37 1.0E-4 0.138361E+00 -ch2o_f01_pbe.inp 37 1.0E-4 0.145756E+00 -ch2o_f02_coulomb_only.inp 37 1.0E-4 0.202091E+00 -#ch2o_f01_pbe.inp 8 1.0E-07 +h2o_f01_pbe_gapwxc.inp 37 1.0E-4 0.138361E+00 +##h2o_f01_pbe.inp 37 1.0E-4 0.145756E+00 +h2o_f02_coulomb_only.inp 37 1.0E-4 0.202091E+00 h2o_t01.inp 0 h2o_t02.inp 0 h2o_t03.inp 0 h2o_t04.inp 0 -#h2o_t05.inp 0 -#h2o_t06.inp 0 -#h2o_t07.inp 0 -#h2o_t08.inp 0 +h2o_t05.inp 0 +h2o_t06.inp 0 +h2o_t07.inp 0 +h2o_t08.inp 0 +h2o_t09.inp 0 +h2o_t10.inp 0 +##h2o_t11.inp 0 +h2o_t12.inp 0 #EOF diff --git a/tests/QS/regtest-tddfpt-force-gapw/ch2o_f01_pbe.inp b/tests/QS/regtest-tddfpt-force-gapw/h2o_f01_pbe.inp similarity index 82% rename from tests/QS/regtest-tddfpt-force-gapw/ch2o_f01_pbe.inp rename to tests/QS/regtest-tddfpt-force-gapw/h2o_f01_pbe.inp index b061c6fb16..bd479dcdb6 100644 --- a/tests/QS/regtest-tddfpt-force-gapw/ch2o_f01_pbe.inp +++ b/tests/QS/regtest-tddfpt-force-gapw/h2o_f01_pbe.inp @@ -9,14 +9,14 @@ &END XC NSTATES 1 MAX_ITER 50 - CONVERGENCE [eV] 1.0e-5 + CONVERGENCE [eV] 1.0e-9 &END TDDFPT &END PROPERTIES &DFT BASIS_SET_FILE_NAME BASIS_SET POTENTIAL_FILE_NAME GTH_POTENTIALS &MGRID - CUTOFF 200 + CUTOFF 280 REL_CUTOFF 60 &END MGRID &QS @@ -47,16 +47,6 @@ &XC_FUNCTIONAL PBE &END XC_FUNCTIONAL &END XC - &PRINT - # &MOMENTS ON - # PERIODIC .FALSE. - # REFERENCE COM - # &END - # &AO_MATRICES - # NDIGITS 12 - # W_MATRIX - # &END AO_MATRICES - &END &POISSON PERIODIC NONE POISSON_SOLVER MT @@ -64,7 +54,7 @@ &END DFT &SUBSYS &CELL - ABC [angstrom] 4.0 4.0 4.0 + ABC [angstrom] 6.0 6.0 6.0 PERIODIC NONE &END &COORD @@ -93,15 +83,13 @@ &GLOBAL PRINT_LEVEL LOW PROJECT td_dipole - RUN_TYPE ENERGY_FORCE -# RUN_TYPE DEBUG + RUN_TYPE DEBUG &END GLOBAL &DEBUG DEBUG_FORCES .TRUE. DEBUG_STRESS_TENSOR .FALSE. DEBUG_DIPOLE .FALSE. -### CHECK_DIPOLE_DIRS Z DEBUG_POLARIZABILITY .FALSE. - DE 0.0001 + DE 0.002 &END diff --git a/tests/QS/regtest-tddfpt-force-gapw/ch2o_f01_pbe_gapwxc.inp b/tests/QS/regtest-tddfpt-force-gapw/h2o_f01_pbe_gapwxc.inp similarity index 100% rename from tests/QS/regtest-tddfpt-force-gapw/ch2o_f01_pbe_gapwxc.inp rename to tests/QS/regtest-tddfpt-force-gapw/h2o_f01_pbe_gapwxc.inp diff --git a/tests/QS/regtest-tddfpt-force-gapw/ch2o_f02_coulomb_only.inp b/tests/QS/regtest-tddfpt-force-gapw/h2o_f02_coulomb_only.inp similarity index 100% rename from tests/QS/regtest-tddfpt-force-gapw/ch2o_f02_coulomb_only.inp rename to tests/QS/regtest-tddfpt-force-gapw/h2o_f02_coulomb_only.inp diff --git a/tests/QS/regtest-tddfpt-force-gapw/h2o_t06.inp b/tests/QS/regtest-tddfpt-force-gapw/h2o_t06.inp index 92358b0e48..a185342590 100644 --- a/tests/QS/regtest-tddfpt-force-gapw/h2o_t06.inp +++ b/tests/QS/regtest-tddfpt-force-gapw/h2o_t06.inp @@ -82,7 +82,7 @@ &DEBUG DEBUG_FORCES .TRUE. DEBUG_STRESS_TENSOR .FALSE. - CHECK_ATOM_FORCE 1 z + CHECK_ATOM_FORCE 1 y STOP_ON_MISMATCH T &END diff --git a/tests/QS/regtest-tddfpt-force-gapw/h2o_t09.inp b/tests/QS/regtest-tddfpt-force-gapw/h2o_t09.inp new file mode 100644 index 0000000000..497828ca71 --- /dev/null +++ b/tests/QS/regtest-tddfpt-force-gapw/h2o_t09.inp @@ -0,0 +1,90 @@ +&FORCE_EVAL + METHOD Quickstep + &PROPERTIES + &TDDFPT + KERNEL FULL + NSTATES 1 + MAX_ITER 50 + CONVERGENCE [eV] 1.0e-7 + RKS_TRIPLETS T + &XC + &XC_FUNCTIONAL NONE + &END XC_FUNCTIONAL + &END XC + &END TDDFPT + &END PROPERTIES + &DFT + BASIS_SET_FILE_NAME BASIS_SET + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 200 + &END MGRID + &QS + METHOD GAPW + EPS_DEFAULT 1.E-10 + &END QS + &EXCITED_STATES T + STATE 1 + DEBUG_FORCES T + &END EXCITED_STATES + &SCF + SCF_GUESS ATOMIC + &OT + PRECONDITIONER FULL_ALL + MINIMIZER DIIS + STEPSIZE 0.1 + &END + &OUTER_SCF + MAX_SCF 20 + EPS_SCF 1.0E-6 + &END + MAX_SCF 10 + EPS_SCF 1.0E-6 + &END SCF + &XC + &XC_FUNCTIONAL PBE + &END XC_FUNCTIONAL + &END XC + &POISSON + PERIODIC NONE + POISSON_SOLVER MT + &END + &END DFT + &SUBSYS + &CELL + ABC [angstrom] 6.0 6.0 6.0 + PERIODIC NONE + &END + &COORD + O 0.000000 0.000000 0.000000 + H 0.000000 -0.757136 0.500545 + H 0.000000 0.757136 0.500545 + &END COORD + &TOPOLOGY + &CENTER_COORDINATES + &END + &END + &KIND H + BASIS_SET DZVP-GTH-PADE + POTENTIAL GTH-PADE-q1 + &END KIND + &KIND O + BASIS_SET DZVP-GTH-PADE + POTENTIAL GTH-PADE-q6 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PRINT_LEVEL LOW + PROJECT td_force + RUN_TYPE DEBUG +&END GLOBAL + +&DEBUG + DEBUG_FORCES .TRUE. + DEBUG_STRESS_TENSOR .FALSE. + CHECK_ATOM_FORCE 1 z + STOP_ON_MISMATCH T +&END + + diff --git a/tests/QS/regtest-tddfpt-force-gapw/h2o_t10.inp b/tests/QS/regtest-tddfpt-force-gapw/h2o_t10.inp new file mode 100644 index 0000000000..3289482efd --- /dev/null +++ b/tests/QS/regtest-tddfpt-force-gapw/h2o_t10.inp @@ -0,0 +1,90 @@ +&FORCE_EVAL + METHOD Quickstep + &PROPERTIES + &TDDFPT + KERNEL FULL + NSTATES 1 + MAX_ITER 50 + CONVERGENCE [eV] 1.0e-7 + RKS_TRIPLETS T + &XC + &XC_FUNCTIONAL PADE + &END XC_FUNCTIONAL + &END XC + &END TDDFPT + &END PROPERTIES + &DFT + BASIS_SET_FILE_NAME BASIS_SET + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 200 + &END MGRID + &QS + METHOD GAPW + EPS_DEFAULT 1.E-10 + &END QS + &EXCITED_STATES T + STATE 1 + DEBUG_FORCES T + &END EXCITED_STATES + &SCF + SCF_GUESS ATOMIC + &OT + PRECONDITIONER FULL_ALL + MINIMIZER DIIS + STEPSIZE 0.1 + &END + &OUTER_SCF + MAX_SCF 20 + EPS_SCF 1.0E-6 + &END + MAX_SCF 10 + EPS_SCF 1.0E-6 + &END SCF + &XC + &XC_FUNCTIONAL PBE + &END XC_FUNCTIONAL + &END XC + &POISSON + PERIODIC NONE + POISSON_SOLVER MT + &END + &END DFT + &SUBSYS + &CELL + ABC [angstrom] 6.0 6.0 6.0 + PERIODIC NONE + &END + &COORD + O 0.000000 0.000000 0.000000 + H 0.000000 -0.757136 0.500545 + H 0.000000 0.757136 0.500545 + &END COORD + &TOPOLOGY + &CENTER_COORDINATES + &END + &END + &KIND H + BASIS_SET DZVP-GTH-PADE + POTENTIAL GTH-PADE-q1 + &END KIND + &KIND O + BASIS_SET DZVP-GTH-PADE + POTENTIAL GTH-PADE-q6 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PRINT_LEVEL LOW + PROJECT td_force + RUN_TYPE DEBUG +&END GLOBAL + +&DEBUG + DEBUG_FORCES .TRUE. + DEBUG_STRESS_TENSOR .FALSE. + CHECK_ATOM_FORCE 1 z + STOP_ON_MISMATCH T +&END + + diff --git a/tests/QS/regtest-tddfpt-force-gapw/h2o_t11.inp b/tests/QS/regtest-tddfpt-force-gapw/h2o_t11.inp new file mode 100644 index 0000000000..986a71057a --- /dev/null +++ b/tests/QS/regtest-tddfpt-force-gapw/h2o_t11.inp @@ -0,0 +1,90 @@ +&FORCE_EVAL + METHOD Quickstep + &PROPERTIES + &TDDFPT + KERNEL FULL + NSTATES 1 + MAX_ITER 50 + CONVERGENCE [eV] 1.0e-9 + RKS_TRIPLETS T + &XC + &XC_FUNCTIONAL PBE + &END XC_FUNCTIONAL + &END XC + &END TDDFPT + &END PROPERTIES + &DFT + BASIS_SET_FILE_NAME ALL_BASIS_SETS + POTENTIAL_FILE_NAME ALL_POTENTIALS + &MGRID + CUTOFF 200 + &END MGRID + &QS + METHOD GAPW + EPS_DEFAULT 1.E-10 + &END QS + &EXCITED_STATES T + STATE 1 + DEBUG_FORCES T + &END EXCITED_STATES + &SCF + SCF_GUESS ATOMIC + &OT + PRECONDITIONER FULL_ALL + MINIMIZER DIIS + STEPSIZE 0.1 + &END + &OUTER_SCF + MAX_SCF 20 + EPS_SCF 1.0E-7 + &END + MAX_SCF 10 + EPS_SCF 1.0E-7 + &END SCF + &XC + &XC_FUNCTIONAL PBE + &END XC_FUNCTIONAL + &END XC + &POISSON + PERIODIC NONE + POISSON_SOLVER MT + &END + &END DFT + &SUBSYS + &CELL + ABC [angstrom] 6.0 6.0 6.0 + PERIODIC NONE + &END + &COORD + O 0.000000 0.000000 0.000000 + H 0.000000 -0.757136 0.500545 + H 0.000000 0.757136 0.500545 + &END COORD + &TOPOLOGY + &CENTER_COORDINATES + &END + &END + &KIND H + BASIS_SET DZVP-ALL + POTENTIAL ALL + &END KIND + &KIND O + BASIS_SET DZVP-ALL + POTENTIAL ALL + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PRINT_LEVEL LOW + PROJECT td_force + RUN_TYPE DEBUG +&END GLOBAL + +&DEBUG + DEBUG_FORCES .TRUE. + DEBUG_STRESS_TENSOR .FALSE. + CHECK_ATOM_FORCE 1 z + STOP_ON_MISMATCH T +&END + + diff --git a/tests/QS/regtest-tddfpt-force-gapw/h2o_t12.inp b/tests/QS/regtest-tddfpt-force-gapw/h2o_t12.inp new file mode 100644 index 0000000000..f44bb0637b --- /dev/null +++ b/tests/QS/regtest-tddfpt-force-gapw/h2o_t12.inp @@ -0,0 +1,90 @@ +&FORCE_EVAL + METHOD Quickstep + &PROPERTIES + &TDDFPT + KERNEL FULL + NSTATES 1 + MAX_ITER 50 + CONVERGENCE [eV] 1.0e-7 + RKS_TRIPLETS T + &XC + &XC_FUNCTIONAL PBE0 + &END XC_FUNCTIONAL + &END XC + &END TDDFPT + &END PROPERTIES + &DFT + BASIS_SET_FILE_NAME BASIS_SET + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 200 + &END MGRID + &QS + METHOD GAPW + EPS_DEFAULT 1.E-10 + &END QS + &EXCITED_STATES T + STATE 1 + DEBUG_FORCES T + &END EXCITED_STATES + &SCF + SCF_GUESS ATOMIC + &OT + PRECONDITIONER FULL_ALL + MINIMIZER DIIS + STEPSIZE 0.02 + &END + &OUTER_SCF + MAX_SCF 10 + EPS_SCF 1.0E-6 + &END + MAX_SCF 10 + EPS_SCF 1.0E-6 + &END SCF + &XC + &XC_FUNCTIONAL PBE0 + &END XC_FUNCTIONAL + &END XC + &POISSON + PERIODIC NONE + POISSON_SOLVER MT + &END + &END DFT + &SUBSYS + &CELL + ABC [angstrom] 6.0 6.0 6.0 + PERIODIC NONE + &END + &COORD + O 0.000000 0.000000 0.000000 + H 0.000000 -0.757136 0.504545 + H 0.000000 0.757136 0.504545 + &END COORD + &TOPOLOGY + &CENTER_COORDINATES + &END + &END + &KIND H + BASIS_SET DZVP-GTH-PADE + POTENTIAL GTH-PADE-q1 + &END KIND + &KIND O + BASIS_SET DZVP-GTH-PADE + POTENTIAL GTH-PADE-q6 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PRINT_LEVEL LOW + PROJECT td_force + RUN_TYPE DEBUG +&END GLOBAL + +&DEBUG + DEBUG_FORCES .TRUE. + DEBUG_STRESS_TENSOR .FALSE. + CHECK_ATOM_FORCE 1 z + STOP_ON_MISMATCH T +&END + + diff --git a/tests/QS/regtest-tddfpt/H2O_GAPW_1.inp b/tests/QS/regtest-tddfpt-gapw/H2O_GAPW_1.inp similarity index 100% rename from tests/QS/regtest-tddfpt/H2O_GAPW_1.inp rename to tests/QS/regtest-tddfpt-gapw/H2O_GAPW_1.inp diff --git a/tests/QS/regtest-tddfpt/H2O_GAPW_1_triplet.inp b/tests/QS/regtest-tddfpt-gapw/H2O_GAPW_1_triplet.inp similarity index 66% rename from tests/QS/regtest-tddfpt/H2O_GAPW_1_triplet.inp rename to tests/QS/regtest-tddfpt-gapw/H2O_GAPW_1_triplet.inp index 6444292ce3..07ea79e7d0 100644 --- a/tests/QS/regtest-tddfpt/H2O_GAPW_1_triplet.inp +++ b/tests/QS/regtest-tddfpt-gapw/H2O_GAPW_1_triplet.inp @@ -1,3 +1,13 @@ +# +# PBE ENERGY : -76.3600550759 +# TDA Triplets +# Excited State 1: Triplet-B1 6.5020 eV 190.69 nm f=0.0000 =2.000 +# 5 -> 6 0.70683 +# Excited State 2: Triplet-A1 8.5315 eV 145.32 nm f=0.0000 =2.000 +# 4 -> 6 0.70572 +# Excited State 3: Triplet-A2 8.5956 eV 144.24 nm f=0.0000 =2.000 +# 5 -> 7 0.70664 +# &GLOBAL PROJECT H2O_GAPW RUN_TYPE ENERGY @@ -47,9 +57,9 @@ PERIODIC NONE &END CELL &COORD - O 0.000000 0.000000 -0.065587 H2O - H 0.000000 -0.757136 0.520545 H2O - H 0.000000 0.757136 0.520545 H2O + O 0.000000 0.000000 0.117226 H2O + H 0.000000 -0.757136 -0.468906 H2O + H 0.000000 0.757136 -0.468906 H2O &END COORD &TOPOLOGY &CENTER_COORDINATES diff --git a/tests/QS/regtest-tddfpt/H2O_GAPW_2.inp b/tests/QS/regtest-tddfpt-gapw/H2O_GAPW_2.inp similarity index 100% rename from tests/QS/regtest-tddfpt/H2O_GAPW_2.inp rename to tests/QS/regtest-tddfpt-gapw/H2O_GAPW_2.inp diff --git a/tests/QS/regtest-tddfpt/H2O_GAPW_3.inp b/tests/QS/regtest-tddfpt-gapw/H2O_GAPW_3.inp similarity index 100% rename from tests/QS/regtest-tddfpt/H2O_GAPW_3.inp rename to tests/QS/regtest-tddfpt-gapw/H2O_GAPW_3.inp diff --git a/tests/QS/regtest-tddfpt/H2O_GAPW_4.inp b/tests/QS/regtest-tddfpt-gapw/H2O_GAPW_4.inp similarity index 100% rename from tests/QS/regtest-tddfpt/H2O_GAPW_4.inp rename to tests/QS/regtest-tddfpt-gapw/H2O_GAPW_4.inp diff --git a/tests/QS/regtest-tddfpt/H2O_GAPW_XC_1.inp b/tests/QS/regtest-tddfpt-gapw/H2O_GAPW_XC_1.inp similarity index 100% rename from tests/QS/regtest-tddfpt/H2O_GAPW_XC_1.inp rename to tests/QS/regtest-tddfpt-gapw/H2O_GAPW_XC_1.inp diff --git a/tests/QS/regtest-tddfpt/H2O_GAPW_XC_2.inp b/tests/QS/regtest-tddfpt-gapw/H2O_GAPW_XC_2.inp similarity index 100% rename from tests/QS/regtest-tddfpt/H2O_GAPW_XC_2.inp rename to tests/QS/regtest-tddfpt-gapw/H2O_GAPW_XC_2.inp diff --git a/tests/QS/regtest-tddfpt/H2O_GAPW_XC_3.inp b/tests/QS/regtest-tddfpt-gapw/H2O_GAPW_XC_3.inp similarity index 100% rename from tests/QS/regtest-tddfpt/H2O_GAPW_XC_3.inp rename to tests/QS/regtest-tddfpt-gapw/H2O_GAPW_XC_3.inp diff --git a/tests/QS/regtest-tddfpt/Ne_GAPW_triplet.inp b/tests/QS/regtest-tddfpt-gapw/Ne_GAPW_triplet.inp similarity index 100% rename from tests/QS/regtest-tddfpt/Ne_GAPW_triplet.inp rename to tests/QS/regtest-tddfpt-gapw/Ne_GAPW_triplet.inp diff --git a/tests/QS/regtest-tddfpt-gapw/TEST_FILES b/tests/QS/regtest-tddfpt-gapw/TEST_FILES new file mode 100644 index 0000000000..783bbb16f0 --- /dev/null +++ b/tests/QS/regtest-tddfpt-gapw/TEST_FILES @@ -0,0 +1,16 @@ +# runs are executed in the same order as in this file +# the second field tells which test should be run in order to compare with the last available output +# 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_GAPW_1.inp 37 4.0E-06 0.543793E+00 +H2O_GAPW_2.inp 37 4.0E-06 0.586179E+00 +H2O_GAPW_3.inp 37 4.0E-06 0.590838E+00 +H2O_GAPW_4.inp 37 4.0E-06 0.619386E+00 +H2O_GAPW_XC_1.inp 37 4.0E-06 0.619451E+00 +H2O_GAPW_XC_2.inp 37 4.0E-06 0.812642E+00 +H2O_GAPW_XC_3.inp 37 4.0E-06 0.836577E+00 +H2O_GAPW_1_triplet.inp 37 4.0E-06 0.504851E+00 +Ne_GAPW_triplet.inp 37 4.0E-06 0.166970E+01 +#EOF diff --git a/tests/QS/regtest-tddfpt-soc/TEST_FILES b/tests/QS/regtest-tddfpt-soc/TEST_FILES index 504e49dfd4..41f88080d7 100644 --- a/tests/QS/regtest-tddfpt-soc/TEST_FILES +++ b/tests/QS/regtest-tddfpt-soc/TEST_FILES @@ -3,9 +3,15 @@ # 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 -regtest-Ar-S.inp 104 4.0E-06 20.8357607201 -regtest-Ar-T.inp 104 4.0E-06 20.8357607201 -regtest-Ne-S.inp 104 4.0E-06 45.3245164440 -regtest-Ne-PPs-GAPW.inp 105 1.0E-05 148.205584 +# OLD TEST RESULTS +## regtest-Ar-S.inp 104 4.0E-06 20.8357607201 +## regtest-Ar-T.inp 104 4.0E-06 20.8357607201 +## regtest-Ne-S.inp 104 4.0E-06 45.3245164440 +## regtest-Ne-PPs-GAPW.inp 105 1.0E-05 148.205584 +# NEW TEST RESULTS please recheck +regtest-Ar-S.inp 104 4.0E-06 20.8315290764 +regtest-Ar-T.inp 104 4.0E-06 20.8315290764 +regtest-Ne-S.inp 104 4.0E-06 45.2284637224 +regtest-Ne-PPs-GAPW.inp 105 1.0E-05 148.045212 #regtest-Ne-PPs-GPW.inp 105 1.0E-01 92.832116 #EOF diff --git a/tests/QS/regtest-tddfpt/TEST_FILES b/tests/QS/regtest-tddfpt/TEST_FILES index 24889f7a4f..2ae412dece 100644 --- a/tests/QS/regtest-tddfpt/TEST_FILES +++ b/tests/QS/regtest-tddfpt/TEST_FILES @@ -17,13 +17,4 @@ NO_tddfpt-t-3.inp 37 4.0E-06 H2O_tddfpt_NTO.inp 37 4.0E-06 0.542432E+00 H2O_tddfpt_NTO_slist.inp 37 4.0E-06 0.542432E+00 # -H2O_GAPW_1.inp 37 4.0E-06 0.543793E+00 -H2O_GAPW_2.inp 37 4.0E-06 0.586179E+00 -H2O_GAPW_3.inp 37 4.0E-06 0.590838E+00 -H2O_GAPW_4.inp 37 4.0E-06 0.619386E+00 -H2O_GAPW_XC_1.inp 37 4.0E-06 0.619451E+00 -H2O_GAPW_XC_2.inp 37 4.0E-06 0.812642E+00 -H2O_GAPW_XC_3.inp 37 4.0E-06 0.836577E+00 -H2O_GAPW_1_triplet.inp 37 4.0E-06 0.505991E+00 -Ne_GAPW_triplet.inp 37 4.0E-06 0.168075E+01 #EOF diff --git a/tests/TEST_DIRS b/tests/TEST_DIRS index 2e5e720a78..511b8d9ab3 100644 --- a/tests/TEST_DIRS +++ b/tests/TEST_DIRS @@ -65,6 +65,8 @@ SIRIUS/regtest-1 sirius QS/regtest-embed libint QS/regtest-pod QS/regtest-tddfpt libint +QS/regtest-tddfpt-force-gapw libint +QS/regtest-tddfpt-gapw libint QS/regtest-tddfpt-admm libint QS/regtest-debug-1 QS/regtest-debug-2 libint @@ -151,7 +153,6 @@ QMMM/SE/regtest QS/regtest-hfx-periodic libint QS/regtest-tddfpt-stda libint QS/regtest-tddfpt-lri libint -QS/regtest-tddfpt-force-gapw libint Fist/regtest-opt QS/regtest-nmr-6 QS/regtest-gpw-1 libint libvori