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

GAPW TDDFPT (#2932)

Browse source
This commit is contained in:
Juerg Hutter 2023-08-18 15:18:36 +02:00 • committed by GitHub
parent 627b89bbdc
commit e709bc9b96
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
37 changed files with 648 additions and 386 deletions
Download patch file Download diff file Expand all files Collapse all files

2
src/qs_environment.F
View file

@ -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

2
src/qs_fxc.F
View file

@ -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.)

23
src/qs_ks_atom.F
View file

@ -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

21
src/qs_ks_reference.F
View file

@ -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

3
src/qs_p_env_methods.F
View file

@ -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)

23
src/qs_rho0_ggrid.F
View file

@ -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

8
src/qs_rho0_methods.F
View file

@ -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, &

19
src/qs_tddfpt2_fhxc.F
View file

@ -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 !

42
src/qs_tddfpt2_fhxc_forces.F
View file

@ -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)

13
src/qs_tddfpt2_forces.F
View file

@ -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)

3
src/qs_tddfpt2_subgroups.F
View file

@ -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)

3
src/qs_tddfpt2_types.F
View file

@ -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)

162
src/qs_vxc_atom.F
View file

@ -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

242
src/response_solver.F
View file

@ -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

3
src/xc/xc_atom.F
View file

@ -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

23
tests/QS/regtest-tddfpt-force-gapw/TEST_FILES
View file

@ -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

22
tests/QS/regtest-tddfpt-force-gapw/ch2o_f01_pbe.inp → tests/QS/regtest-tddfpt-force-gapw/h2o_f01_pbe.inp
View file

@ -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

0
tests/QS/regtest-tddfpt-force-gapw/ch2o_f01_pbe_gapwxc.inp → tests/QS/regtest-tddfpt-force-gapw/h2o_f01_pbe_gapwxc.inp
View file

0
tests/QS/regtest-tddfpt-force-gapw/ch2o_f02_coulomb_only.inp → tests/QS/regtest-tddfpt-force-gapw/h2o_f02_coulomb_only.inp
View file

2
tests/QS/regtest-tddfpt-force-gapw/h2o_t06.inp
View file

@ -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

90
tests/QS/regtest-tddfpt-force-gapw/h2o_t09.inp Normal file
View file

@ -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

90
tests/QS/regtest-tddfpt-force-gapw/h2o_t10.inp Normal file
View file

@ -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

90
tests/QS/regtest-tddfpt-force-gapw/h2o_t11.inp Normal file
View file

@ -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

90
tests/QS/regtest-tddfpt-force-gapw/h2o_t12.inp Normal file
View file

@ -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

0
tests/QS/regtest-tddfpt/H2O_GAPW_1.inp → tests/QS/regtest-tddfpt-gapw/H2O_GAPW_1.inp
View file

16
tests/QS/regtest-tddfpt/H2O_GAPW_1_triplet.inp → tests/QS/regtest-tddfpt-gapw/H2O_GAPW_1_triplet.inp
View file

@ -1,3 +1,13 @@
#
# PBE ENERGY : -76.3600550759
# TDA Triplets
# Excited State 1: Triplet-B1 6.5020 eV 190.69 nm f=0.0000 <S**2>=2.000
# 5 -> 6 0.70683
# Excited State 2: Triplet-A1 8.5315 eV 145.32 nm f=0.0000 <S**2>=2.000
# 4 -> 6 0.70572
# Excited State 3: Triplet-A2 8.5956 eV 144.24 nm f=0.0000 <S**2>=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

0
tests/QS/regtest-tddfpt/H2O_GAPW_2.inp → tests/QS/regtest-tddfpt-gapw/H2O_GAPW_2.inp
View file

0
tests/QS/regtest-tddfpt/H2O_GAPW_3.inp → tests/QS/regtest-tddfpt-gapw/H2O_GAPW_3.inp
View file

0
tests/QS/regtest-tddfpt/H2O_GAPW_4.inp → tests/QS/regtest-tddfpt-gapw/H2O_GAPW_4.inp
View file

0
tests/QS/regtest-tddfpt/H2O_GAPW_XC_1.inp → tests/QS/regtest-tddfpt-gapw/H2O_GAPW_XC_1.inp
View file

0
tests/QS/regtest-tddfpt/H2O_GAPW_XC_2.inp → tests/QS/regtest-tddfpt-gapw/H2O_GAPW_XC_2.inp
View file

0
tests/QS/regtest-tddfpt/H2O_GAPW_XC_3.inp → tests/QS/regtest-tddfpt-gapw/H2O_GAPW_XC_3.inp
View file

0
tests/QS/regtest-tddfpt/Ne_GAPW_triplet.inp → tests/QS/regtest-tddfpt-gapw/Ne_GAPW_triplet.inp
View file

16
tests/QS/regtest-tddfpt-gapw/TEST_FILES Normal file
View file

@ -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

14
tests/QS/regtest-tddfpt-soc/TEST_FILES
View file

@ -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

9
tests/QS/regtest-tddfpt/TEST_FILES
View file

@ -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

3
tests/TEST_DIRS
View file

@ -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