comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-26 13:15:21 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

SF-TDDFT implementation (#4446)

Browse source 
The linear-response TDDFT module was extended to spin-flip excitations. This extension is
capable to make collinear and noncollinear SF-TDDFT calculations. Some key details of the
implementation are:

  - The noncollinear kernel used is the Wang & Ziegler kernel with a screening in the
    denominator to stabilize it
  - The noncollinear kernel derivative is calculated by finite differences
This commit is contained in:
Luis I Hernández Segura 2025-09-21 10:31:49 +02:00 • committed by GitHub
parent 916bf6ee98
commit b69419fcbe
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
39 changed files with 3072 additions and 663 deletions
Download patch file Download diff file Expand all files Collapse all files

2
src/cp_control_types.F
View file

@ -499,6 +499,8 @@ MODULE cp_control_types
INTEGER :: nprocs = 0
!> type of kernel function/approximation to use
INTEGER :: kernel = 0
!> type of spin excitations to compute
INTEGER :: spinflip = 0
!> for full kernel, do we have HFX/ADMM
LOGICAL :: do_hfx = .FALSE.
LOGICAL :: do_admm = .FALSE.

1
src/cp_control_utils.F
View file

@ -1674,6 +1674,7 @@ CONTAINS
CALL section_vals_val_get(t_section, "NLUMO", i_val=t_control%nlumo)
CALL section_vals_val_get(t_section, "NPROC_STATE", i_val=t_control%nprocs)
CALL section_vals_val_get(t_section, "KERNEL", i_val=t_control%kernel)
CALL section_vals_val_get(t_section, "SPINFLIP", i_val=t_control%spinflip)
CALL section_vals_val_get(t_section, "OE_CORR", i_val=t_control%oe_corr)
CALL section_vals_val_get(t_section, "EV_SHIFT", r_val=t_control%ev_shift)
CALL section_vals_val_get(t_section, "EOS_SHIFT", r_val=t_control%eos_shift)

3
src/excited_states.F
View file

@ -106,6 +106,7 @@ CONTAINS
CALL set_qs_env(qs_env, force=lr_force)
!
tdlr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%TDDFPT%LINRES")
! Solve the Z-vector linear equation system
CALL response_equation(qs_env, p_env, ex_env%cpmos, unit_nr, tdlr_section)
!
CALL get_qs_env(qs_env, dft_control=dft_control)
@ -120,7 +121,7 @@ CONTAINS
CALL response_force(qs_env=qs_env, vh_rspace=ex_env%vh_rspace, &
vxc_rspace=ex_env%vxc_rspace, vtau_rspace=ex_env%vtau_rspace, &
vadmm_rspace=ex_env%vadmm_rspace, matrix_hz=ex_env%matrix_hz, &
matrix_pz=ex_env%matrix_px1, matrix_pz_admm=p_env%p1_admm, &
matrix_pz=p_env%p1, matrix_pz_admm=p_env%p1_admm, &
matrix_wz=p_env%w1, &
p_env=p_env, ex_env=ex_env, &
debug=ex_env%debug_forces)

33
src/hfx_derivatives.F
View file

@ -113,7 +113,7 @@ CONTAINS
! **************************************************************************************************
SUBROUTINE derivatives_four_center(qs_env, rho_ao, rho_ao_resp, hfx_section, para_env, &
irep, use_virial, adiabatic_rescale_factor, resp_only, &
external_x_data)
external_x_data, nspins)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao
@ -125,6 +125,7 @@ CONTAINS
REAL(dp), INTENT(IN), OPTIONAL :: adiabatic_rescale_factor
LOGICAL, INTENT(IN), OPTIONAL :: resp_only
TYPE(hfx_type), DIMENSION(:, :), POINTER, OPTIONAL :: external_x_data
INTEGER, INTENT(IN), OPTIONAL :: nspins
CHARACTER(LEN=*), PARAMETER :: routineN = 'derivatives_four_center'
@ -138,7 +139,7 @@ CONTAINS
INTEGER :: kind_kind_idx, kkind, kset, l_atom, l_max, latom, latom_block, latom_end, &
latom_start, lkind, lset, max_am, max_pgf, max_set, my_bin_id, my_bin_size, my_thread_id, &
n_threads, natom, nbits, ncob, ncos_max, nints, nkimages, nkind, nneighbors, nseta, &
nsetb, nsgf_max, nspins, sgfb, shm_task_counter, shm_total_bins, sphi_a_u1, sphi_a_u2, &
nsetb, nsgf_max, my_nspins, sgfb, shm_task_counter, shm_total_bins, sphi_a_u1, sphi_a_u2, &
sphi_a_u3, sphi_b_u1, sphi_b_u2, sphi_b_u3, sphi_c_u1, sphi_c_u2, sphi_c_u3, sphi_d_u1, &
sphi_d_u2, sphi_d_u3, swap_id, tmp_i4, unit_id
INTEGER(int_8) :: atom_block, counter, estimate_to_store_int, max_val_memory, &
@ -239,7 +240,11 @@ CONTAINS
cell=cell, &
dft_control=dft_control)
nspins = dft_control%nspins
IF (PRESENT(nspins)) THEN
my_nspins = nspins
ELSE
my_nspins = dft_control%nspins
END IF
nkimages = dft_control%nimages
CPASSERT(nkimages == 1)
@ -337,7 +342,7 @@ CONTAINS
!$OMP shm_pmax_atom,&
!$OMP use_virial,&
!$OMP do_print_load_balance_info,&
!$OMP nkimages,nspins,my_x_data)&
!$OMP nkimages,my_nspins,my_x_data)&
!$OMP PRIVATE(actual_x_data,atomic_kind_set,atomic_offset_ac,atomic_offset_ad,atomic_offset_bc,&
!$OMP atomic_offset_bd,atom_of_kind,basis_info,&
!$OMP basis_parameter,bin,bins_left,bintime_start,bintime_stop,bits_max_val,buffer_left,buffer_overflow,buffer_size,&
@ -514,7 +519,7 @@ CONTAINS
antisymmetric=is_anti_symmetric)
END IF
IF (nspins == 2) THEN
IF (my_nspins == 2) THEN
ALLOCATE (full_density_beta(shm_block_offset(ncpu + 1), 1))
CALL get_full_density(para_env, full_density_beta(:, 1), rho_ao(2, 1)%matrix, shm_number_of_p_entries, &
shm_block_offset, &
@ -555,7 +560,7 @@ CONTAINS
! maybe in the future
IF (with_resp_density .AND. .NOT. my_resp_only) THEN
full_density_alpha(:, 1) = full_density_alpha(:, 1) - full_density_resp
IF (nspins == 2) THEN
IF (my_nspins == 2) THEN
full_density_beta(:, 1) = &
full_density_beta(:, 1) - full_density_resp_beta
END IF
@ -638,7 +643,7 @@ CONTAINS
ALLOCATE (pad_buf_resp(nsgf_max**2))
ALLOCATE (pac_buf_resp(nsgf_max**2))
END IF
IF (nspins == 2) THEN
IF (my_nspins == 2) THEN
ALLOCATE (pbd_buf_beta(nsgf_max**2))
ALLOCATE (pbc_buf_beta(nsgf_max**2))
ALLOCATE (pad_buf_beta(nsgf_max**2))
@ -1035,7 +1040,7 @@ CONTAINS
forces_map(4, 1) = lkind
forces_map(4, 2) = l_atom
IF (nspins == 1) THEN
IF (my_nspins == 1) THEN
fac = 0.25_dp*hf_fraction
ELSE
fac = 0.5_dp*hf_fraction
@ -1339,8 +1344,8 @@ CONTAINS
spherical_estimate = MAX(spherical_estimate, ABS(primitive_forces(i)))
END DO
spherical_estimate_virial = 0.0_dp
IF (use_virial) THEN
spherical_estimate_virial = 0.0_dp
DO i = 1, 3*nints
spherical_estimate_virial = MAX(spherical_estimate_virial, ABS(primitive_forces_virial(i)))
END DO
@ -1414,7 +1419,7 @@ CONTAINS
offset_ad_set, offset_ac_set, atomic_offset_bd, atomic_offset_bc, &
atomic_offset_ad, atomic_offset_ac, is_anti_symmetric)
END IF
IF (nspins == 2) THEN
IF (my_nspins == 2) THEN
IF (with_resp_density) THEN
CALL prefetch_density_matrix(nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), &
full_density_beta(:, 1), pbd_buf_beta, pbc_buf_beta, &
@ -1453,7 +1458,7 @@ CONTAINS
pbd_buf, pbc_buf, pad_buf, pac_buf, fac, &
T2, force, forces_map, coord)
END IF
IF (nspins == 2) THEN
IF (my_nspins == 2) THEN
IF (with_resp_density) THEN
CALL update_forces(nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), &
pbd_buf_beta, pbc_buf_beta, pad_buf_beta, pac_buf_beta, &
@ -1486,7 +1491,7 @@ CONTAINS
pbd_buf, pbc_buf, pad_buf, pac_buf, fac, &
T2, tmp_virial, coord, i)
END IF
IF (nspins == 2) THEN
IF (my_nspins == 2) THEN
IF (with_resp_density) THEN
CALL update_virial(nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), &
pbd_buf_beta, pbc_buf_beta, pad_buf_beta, pac_buf_beta, &
@ -1684,7 +1689,7 @@ CONTAINS
IF (with_resp_density) THEN
DEALLOCATE (full_density_resp)
END IF
IF (nspins == 2) THEN
IF (my_nspins == 2) THEN
DEALLOCATE (full_density_beta)
IF (with_resp_density) THEN
DEALLOCATE (full_density_resp_beta)
@ -1728,7 +1733,7 @@ CONTAINS
DEALLOCATE (nimages)
IF (nspins == 2) THEN
IF (my_nspins == 2) THEN
DEALLOCATE (pbd_buf_beta, pbc_buf_beta, pad_buf_beta, pac_buf_beta)
IF (with_resp_density) THEN
DEALLOCATE (pbd_buf_resp_beta)

46
src/hfx_energy_potential.F
View file

@ -133,7 +133,7 @@ CONTAINS
! **************************************************************************************************
SUBROUTINE integrate_four_center(qs_env, x_data, ks_matrix, ehfx, rho_ao, hfx_section, para_env, &
geometry_did_change, irep, distribute_fock_matrix, &
ispin)
ispin, nspins)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
@ -146,6 +146,7 @@ CONTAINS
INTEGER :: irep
LOGICAL, INTENT(IN) :: distribute_fock_matrix
INTEGER, INTENT(IN) :: ispin
INTEGER, INTENT(IN), OPTIONAL :: nspins
CHARACTER(LEN=*), PARAMETER :: routineN = 'integrate_four_center'
@ -161,7 +162,7 @@ CONTAINS
INTEGER :: katom_start, kind_kind_idx, kkind, kset, l_max, latom, latom_block, latom_end, &
latom_start, lkind, lset, ma, max_am, max_pgf, max_set, mb, my_bin_id, my_bin_size, &
my_thread_id, n_threads, natom, nbits, ncob, ncos_max, nints, nkimages, nkind, &
nneighbors, nseta, nsetb, nsgf_max, nspins, pa, sgfb, shm_task_counter, shm_total_bins, &
nneighbors, nseta, nsetb, nsgf_max, my_nspins, pa, sgfb, shm_task_counter, shm_total_bins, &
sphi_a_u1, sphi_a_u2, sphi_a_u3, sphi_b_u1, sphi_b_u2, sphi_b_u3, sphi_c_u1, sphi_c_u2, &
sphi_c_u3, sphi_d_u1, sphi_d_u2, sphi_d_u3, swap_id, tmp_i4, unit_id
INTEGER(int_8) :: atom_block, counter, estimate_to_store_int, max_val_memory, &
@ -313,6 +314,10 @@ CONTAINS
n_threads = 1
!$ n_threads = omp_get_max_threads()
! This initialization is needed to prevent a segmentation fault. The correct assigment is done below
my_nspins = 0
IF (PRESENT(nspins)) my_nspins = nspins
shm_neris_total = 0
shm_nprim_ints = 0
shm_neris_onthefly = 0
@ -374,7 +379,8 @@ CONTAINS
!$OMP shm_pmax_block,&
!$OMP shm_atomic_pair_list,&
!$OMP shm_mem_compression_counter,&
!$OMP do_print_load_balance_info) &
!$OMP do_print_load_balance_info,&
!$OMP my_nspins) &
!$OMP PRIVATE(ln_10,i_thread,actual_x_data,do_periodic,screening_parameter,potential_parameter,&
!$OMP general_parameter,load_balance_parameter,memory_parameter,cache_size,bits_max_val,&
!$OMP basis_parameter,basis_info,treat_lsd_in_core,ncpu,n_processes,neris_total,neris_incore,&
@ -382,7 +388,7 @@ CONTAINS
!$OMP compression_factor_disk,nprim_ints,neris_tmp,max_val_memory,max_am,do_p_screening,&
!$OMP max_set,particle_set,atomic_kind_set,natom,kind_of,ncos_max,nsgf_max,ikind,&
!$OMP nseta,npgfa,la_max,nsgfa,primitive_integrals,pbd_buf,pbc_buf,pad_buf,pac_buf,kbd_buf,kbc_buf,&
!$OMP kad_buf,kac_buf,ee_work,ee_work2,ee_buffer1,ee_buffer2,ee_primitives_tmp,nspins,max_contraction,&
!$OMP kad_buf,kac_buf,ee_work,ee_work2,ee_buffer1,ee_buffer2,ee_primitives_tmp,max_contraction,&
!$OMP max_pgf,jkind,lb_max,nsetb,npgfb,first_sgfb,sphi_b,nsgfb,ncob,sgfb,nneighbors,pgf_list_ij,pgf_list_kl,&
!$OMP pgf_product_list,nimages,ks_fully_occ,subtr_size_mb,use_disk_storage,counter,do_disk_storage,&
!$OMP maxval_container_disk,maxval_cache_disk,integral_containers_disk,integral_caches_disk,eps_schwarz,&
@ -517,7 +523,7 @@ CONTAINS
ALLOCATE (ee_buffer2(ncos_max**4))
ALLOCATE (ee_primitives_tmp(nsgf_max**4))
nspins = dft_control%nspins
IF (my_nspins .EQ. 0) my_nspins = dft_control%nspins
ALLOCATE (max_contraction(max_set, natom))
@ -607,7 +613,7 @@ CONTAINS
subtr_size_mb = 2_int_8*shm_block_offset(ncpu + 1)
! ** if non restricted ==> alpha/beta spin
IF (.NOT. treat_lsd_in_core) THEN
IF (nspins == 2) subtr_size_mb = subtr_size_mb*2_int_8
IF (my_nspins == 2) subtr_size_mb = subtr_size_mb*2_int_8
END IF
! ** Initial P only MAX(alpha,beta) (shared)
IF (do_p_screening .OR. screening_parameter%do_p_screening_forces) THEN
@ -719,7 +725,7 @@ CONTAINS
!! Get the full density from all the processors
NULLIFY (full_density_alpha, full_density_beta)
ALLOCATE (full_density_alpha(shm_block_offset(ncpu + 1), nkimages))
IF (.NOT. treat_lsd_in_core .OR. nspins == 1) THEN
IF (.NOT. treat_lsd_in_core .OR. my_nspins == 1) THEN
CALL timeset(routineN//"_getP", handle_getP)
DO img = 1, nkimages
CALL get_full_density(para_env, full_density_alpha(:, img), rho_ao(ispin, img)%matrix, shm_number_of_p_entries, &
@ -727,7 +733,7 @@ CONTAINS
kind_of, basis_parameter, get_max_vals_spin=.FALSE., antisymmetric=is_anti_symmetric)
END DO
IF (nspins == 2) THEN
IF (my_nspins == 2) THEN
ALLOCATE (full_density_beta(shm_block_offset(ncpu + 1), nkimages))
DO img = 1, nkimages
CALL get_full_density(para_env, full_density_beta(:, img), rho_ao(2, img)%matrix, shm_number_of_p_entries, &
@ -798,7 +804,7 @@ CONTAINS
full_ks_alpha = 0.0_dp
IF (.NOT. treat_lsd_in_core) THEN
IF (nspins == 2) THEN
IF (my_nspins == 2) THEN
ALLOCATE (shm_master_x_data%full_ks_beta(shm_block_offset(ncpu + 1), nkimages))
full_ks_beta => shm_master_x_data%full_ks_beta
full_ks_beta = 0.0_dp
@ -840,7 +846,7 @@ CONTAINS
!$OMP BARRIER
!! Initialize a prefactor depending on the fraction and the number of spins
IF (nspins == 1) THEN
IF (my_nspins == 1) THEN
fac = 0.5_dp*hf_fraction
ELSE
fac = 1.0_dp*hf_fraction
@ -1679,7 +1685,7 @@ CONTAINS
iset, jset, kset, lset, offset_bd_set, offset_bc_set, offset_ad_set, offset_ac_set, &
atomic_offset_bd, atomic_offset_bc, atomic_offset_ad, atomic_offset_ac)
IF (.NOT. treat_lsd_in_core) THEN
IF (nspins == 2) THEN
IF (my_nspins == 2) THEN
CALL update_fock_matrix( &
nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), &
fac, symm_fac, full_density_beta(:, 1), full_ks_beta(:, 1), &
@ -1699,7 +1705,7 @@ CONTAINS
iset, jset, kset, lset, offset_bd_set, offset_bc_set, offset_ad_set, offset_ac_set, &
atomic_offset_bd, atomic_offset_bc, atomic_offset_ad, atomic_offset_ac)
IF (.NOT. treat_lsd_in_core) THEN
IF (nspins == 2) THEN
IF (my_nspins == 2) THEN
CALL update_fock_matrix_as( &
nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), &
fac, symm_fac, full_density_beta(:, 1), full_ks_beta(:, 1), &
@ -1807,7 +1813,7 @@ CONTAINS
mb_size_p = shm_block_offset(ncpu + 1)/1024/128
mb_size_f = shm_block_offset(ncpu + 1)/1024/128
IF (.NOT. treat_lsd_in_core) THEN
IF (nspins == 2) THEN
IF (my_nspins == 2) THEN
mb_size_f = mb_size_f*2
mb_size_p = mb_size_p*2
END IF
@ -1872,7 +1878,7 @@ CONTAINS
IF (i > j) THEN
full_ks_alpha(i, img) = (full_ks_alpha(i, img) + full_ks_alpha(j, img)*afac)
full_ks_alpha(j, img) = full_ks_alpha(i, img)*afac
IF (.NOT. treat_lsd_in_core .AND. nspins == 2) THEN
IF (.NOT. treat_lsd_in_core .AND. my_nspins == 2) THEN
full_ks_beta(i, img) = (full_ks_beta(i, img) + full_ks_beta(j, img)*afac)
full_ks_beta(j, img) = full_ks_beta(i, img)*afac
END IF
@ -1880,7 +1886,7 @@ CONTAINS
! ** For adiabatically rescaled functionals we need the energy coming from the diagonal elements
IF (i == j) THEN
ene_x_aa_diag = ene_x_aa_diag + full_ks_alpha(i, img)*full_density_alpha(i, img)
IF (.NOT. treat_lsd_in_core .AND. nspins == 2) THEN
IF (.NOT. treat_lsd_in_core .AND. my_nspins == 2) THEN
ene_x_bb_diag = ene_x_bb_diag + full_ks_beta(i, img)*full_density_beta(i, img)
END IF
END IF
@ -1908,7 +1914,7 @@ CONTAINS
NULLIFY (full_ks_alpha)
DEALLOCATE (shm_master_x_data%full_ks_alpha)
IF (.NOT. treat_lsd_in_core) THEN
IF (nspins == 2) THEN
IF (my_nspins == 2) THEN
DO img = 1, nkimages
CALL distribute_ks_matrix(para_env, full_ks_beta(:, img), ks_matrix(2, img)%matrix, shm_number_of_p_entries, &
shm_block_offset, kind_of, basis_parameter, &
@ -1936,7 +1942,7 @@ CONTAINS
END IF
ene_x_bb = 0.0_dp
IF (nspins == 2 .AND. .NOT. treat_lsd_in_core) THEN
IF (my_nspins == 2 .AND. .NOT. treat_lsd_in_core) THEN
DO img = 1, nkimages
CALL dbcsr_dot_threadsafe(ks_matrix(2, img)%matrix, rho_ao(2, img)%matrix, etmp)
ene_x_bb = ene_x_bb + etmp
@ -1961,7 +1967,7 @@ CONTAINS
END DO
! ** Now correct
ene_x_aa = (ene_x_aa + ene_x_aa_diag)*0.5_dp
IF (nspins == 2) THEN
IF (my_nspins == 2) THEN
DO img = 1, nkimages
DO pa = 1, SIZE(full_ks_beta, 1)
ene_x_bb = ene_x_bb + full_ks_beta(pa, img)*full_density_beta(pa, img)
@ -1971,7 +1977,7 @@ CONTAINS
ene_x_bb = (ene_x_bb + ene_x_bb_diag)*0.5_dp
END IF
CALL para_env%sum(ene_x_aa)
IF (nspins == 2) CALL para_env%sum(ene_x_bb)
IF (my_nspins == 2) CALL para_env%sum(ene_x_bb)
ehfx = 0.5_dp*(ene_x_aa + ene_x_bb)
END IF
@ -2126,7 +2132,7 @@ CONTAINS
!$OMP MASTER
DEALLOCATE (full_density_alpha)
IF (.NOT. treat_lsd_in_core) THEN
IF (nspins == 2) THEN
IF (my_nspins == 2) THEN
DEALLOCATE (full_density_beta)
END IF
END IF

3
src/input_constants.F
View file

@ -583,6 +583,9 @@ MODULE input_constants
INTEGER, PARAMETER, PUBLIC :: tddfpt_kernel_none = 2, &
tddfpt_kernel_full = 1, &
tddfpt_kernel_stda = 0
INTEGER, PARAMETER, PUBLIC :: no_sf_tddfpt = 0, &
tddfpt_sf_col = 1, &
tddfpt_sf_noncol = 2
INTEGER, PARAMETER, PUBLIC :: oe_none = 0, &
oe_lb = 1, &
oe_gllb = 2, &

16
src/input_cp2k_properties_dft.F
View file

@ -38,8 +38,8 @@ MODULE input_cp2k_properties_dft
ot_precond_full_single, ot_precond_full_single_inverse, ot_precond_none, &
ot_precond_s_inverse, scan_x, scan_xy, scan_xyz, scan_xz, scan_y, scan_yz, scan_z, &
tddfpt_dipole_berry, tddfpt_dipole_length, tddfpt_dipole_velocity, tddfpt_kernel_full, &
tddfpt_kernel_none, tddfpt_kernel_stda, use_mom_ref_coac, use_mom_ref_com, &
use_mom_ref_user, use_mom_ref_zero
tddfpt_kernel_none, tddfpt_kernel_stda, no_sf_tddfpt, tddfpt_sf_col, tddfpt_sf_noncol, &
use_mom_ref_coac, use_mom_ref_com, use_mom_ref_user, use_mom_ref_zero
USE input_cp2k_atprop, ONLY: create_atprop_section
USE input_cp2k_dft, ONLY: create_interp_section, &
create_mgrid_section
@ -1508,7 +1508,7 @@ CONTAINS
"Density Functional Perturbation Theory. "// &
"Current implementation works for hybrid functionals. "// &
"Can be used with Gaussian and Plane Waves (GPW) method only.", &
n_keywords=12, n_subsections=4, repeats=.FALSE., &
n_keywords=14, n_subsections=4, repeats=.FALSE., &
citations=(/Iannuzzi2005/))
NULLIFY (keyword, print_key, subsection)
@ -1570,6 +1570,16 @@ CONTAINS
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
! spin-flip TDDFPT options
CALL keyword_create(keyword, __LOCATION__, name="SPINFLIP", &
description="Selects the type of spin-flip TDDFPT kernel", &
usage="SPINFLIP NONCOLLINEAR", &
enum_c_vals=s2a("NONE", "COLLINEAR", "NONCOLLINEAR"), &
enum_i_vals=(/no_sf_tddfpt, tddfpt_sf_col, tddfpt_sf_noncol/), &
default_i_val=no_sf_tddfpt)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="OE_CORR", &
description="Orbital energy correction potential.", &
enum_c_vals=s2a("NONE", "LB94", "GLLB", "SAOP", "SHIFT"), &

2
src/qs_energy.F
View file

@ -159,8 +159,10 @@ CONTAINS
END IF
IF (.NOT. loverlap_deltat) THEN
! Calculate energy, response vector and some contributions to the force
CALL qs_energies_properties(qs_env, calc_forces)
! Update total energy of the selected excited state
CALL excited_state_energy(qs_env, calculate_forces=.FALSE.)
END IF

5
src/qs_energy_utils.F
View file

@ -84,7 +84,8 @@ CONTAINS
! **************************************************************************************************
!> \brief Refactoring of qs_energies_scf. Moves computation of properties
!> into separate subroutine
!> into separate subroutine. Calculates energy and if calc_forces then the response vector
!> and some contributions to the force
!> \param qs_env ...
!> \param calc_forces ...
!> \par History
@ -191,10 +192,12 @@ CONTAINS
END IF
! Compute Linear Response properties as post-scf
! Initialization of TDDFPT calculation
IF (.NOT. qs_env%linres_run) THEN
CALL linres_calculation_low(qs_env)
END IF
! Calculates the energy, response vector and some contributions to the force
IF (dft_control%tddfpt2_control%enabled) THEN
CALL tddfpt(qs_env, calc_forces)
END IF

4
src/qs_force.F
View file

@ -197,7 +197,7 @@ CONTAINS
dft_control%qs_control%cdft_control%save_pot = .TRUE.
END IF
! recalculate energy with forces
! recalculate energy and the response vector for the Z-vector linear equation system if calc_force = .true.
CALL qs_energies(qs_env, calc_forces=.TRUE.)
NULLIFY (para_env)
@ -374,6 +374,8 @@ CONTAINS
END IF
! Excited state forces
! Solve the response linear equation system for the Z-vector method
! and calculate remaining terms of the force
CALL excited_state_energy(qs_env, calculate_forces=.TRUE.)
! replicate forces (get current pointer)

247
src/qs_fxc.F
View file

@ -27,7 +27,9 @@
MODULE qs_fxc
USE cp_control_types, ONLY: dft_control_type
USE input_section_types, ONLY: section_vals_get_subs_vals,&
USE input_section_types, ONLY: section_get_ival,&
section_get_rval,&
section_vals_get_subs_vals,&
section_vals_type
USE kinds, ONLY: dp
USE pw_env_types, ONLY: pw_env_get,&
@ -41,20 +43,26 @@ MODULE qs_fxc
USE qs_ks_types, ONLY: get_ks_env,&
qs_ks_env_type
USE qs_rho_methods, ONLY: qs_rho_copy,&
qs_rho_scale_and_add
qs_rho_scale_and_add,&
qs_rho_scale_and_add_b
USE qs_rho_types, ONLY: qs_rho_create,&
qs_rho_get,&
qs_rho_release,&
qs_rho_type
USE qs_vxc, ONLY: qs_vxc_create
USE xc, ONLY: xc_calc_2nd_deriv,&
xc_prep_2nd_deriv
xc_calc_2nd_deriv_analytical,&
xc_calc_3rd_deriv_analytical,&
xc_prep_2nd_deriv,&
xc_prep_3rd_deriv
USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
xc_dset_release
USE xc_derivatives, ONLY: xc_functionals_get_needs
USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
USE xc_rho_set_types, ONLY: xc_rho_set_release,&
xc_rho_set_type
USE xc_rho_set_types, ONLY: xc_rho_set_create,&
xc_rho_set_release,&
xc_rho_set_type,&
xc_rho_set_update
#include "./base/base_uses.f90"
IMPLICIT NONE
@ -62,7 +70,7 @@ MODULE qs_fxc
PRIVATE
! *** Public subroutines ***
PUBLIC :: qs_fxc_fdiff, qs_fxc_analytic, qs_fgxc_gdiff, qs_fgxc_create, qs_fgxc_release
PUBLIC :: qs_fxc_fdiff, qs_fxc_analytic, qs_fgxc_gdiff, qs_fgxc_analytic, qs_fgxc_create, qs_fgxc_release
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_fxc'
@ -80,8 +88,9 @@ CONTAINS
!> \param is_triplet ...
!> \param v_xc ...
!> \param v_xc_tau ...
!> \param spinflip ...
! **************************************************************************************************
SUBROUTINE qs_fxc_analytic(rho0, rho1_r, tau1_r, xc_section, auxbas_pw_pool, is_triplet, v_xc, v_xc_tau)
SUBROUTINE qs_fxc_analytic(rho0, rho1_r, tau1_r, xc_section, auxbas_pw_pool, is_triplet, v_xc, v_xc_tau, spinflip)
TYPE(qs_rho_type), POINTER :: rho0
TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho1_r, tau1_r
@ -89,12 +98,13 @@ CONTAINS
TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
LOGICAL, INTENT(IN) :: is_triplet
TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_xc, v_xc_tau
LOGICAL, OPTIONAL :: spinflip
CHARACTER(len=*), PARAMETER :: routineN = 'qs_fxc_analytic'
INTEGER :: handle, nspins
INTEGER, DIMENSION(2, 3) :: bo
LOGICAL :: lsd
LOGICAL :: do_sf, lsd
REAL(KIND=dp) :: fac
TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho0_g, rho1_g
TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho0_r, tau0_r
@ -108,6 +118,9 @@ CONTAINS
CPASSERT(.NOT. ASSOCIATED(v_xc))
CPASSERT(.NOT. ASSOCIATED(v_xc_tau))
do_sf = .FALSE.
IF (PRESENT(spinflip)) do_sf = spinflip
CALL qs_rho_get(rho0, rho_r=rho0_r, rho_g=rho0_g, tau_r=tau0_r)
nspins = SIZE(rho0_r)
@ -119,10 +132,12 @@ CONTAINS
bo = rho1_r(1)%pw_grid%bounds_local
xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
needs = xc_functionals_get_needs(xc_fun_section, lsd, .TRUE.)
! calculate the arguments needed by the functionals
! calculate the arguments needed by the functionals and the values of the functional on the grid
CALL xc_prep_2nd_deriv(deriv_set, rho0_set, rho0_r, auxbas_pw_pool, xc_section=xc_section, tau_r=tau0_r)
! Folds the density rho1 with the functional
CALL xc_calc_2nd_deriv(v_xc, v_xc_tau, deriv_set, rho0_set, rho1_r, rho1_g, tau1_r, &
auxbas_pw_pool, xc_section=xc_section, gapw=.FALSE., do_triplet=is_triplet)
auxbas_pw_pool, xc_section=xc_section, gapw=.FALSE., do_triplet=is_triplet, &
do_sf=do_sf)
CALL xc_dset_release(deriv_set)
CALL xc_rho_set_release(rho0_set)
@ -265,6 +280,161 @@ CONTAINS
END SUBROUTINE qs_fxc_fdiff
! **************************************************************************************************
!> \brief Calculates the values at the grid points in real space (r_i), of the second and third
!> functional derivatives of the exchange-correlation energy functional.
!> fxc_rho(r_i) = fxc[n](r_i)*n^(1)(r_i) ! Second functional derivative
!> gxc_rho(r_i) = n^(1)(r_i)*gxc[n](r_i)*n^(1)(r_i) ! Third functional derivative
!> \param rho0_struct Ground state density, n(r).
!> \param rho1_struct Density used to fold the functional derivatives, n^(1)(r).
!> \param xc_section ...
!> \param pw_pool ...
!> \param fxc_rho Second functional derivative with respect to the density, n(r).
!> \param fxc_tau mGGA contribution to the second functional derivative with respect to the density.
!> \param gxc_rho Third functional derivative with respect to the density, n(r).
!> \param gxc_tau mGGA contribution to the third functional derivative with respect to the density.
!> \param spinflip Flag used to activate the spin-flip noncollinear kernel and kernel derivatives.
!> \par History
!> * 07.2024 Created [LHS]
! **************************************************************************************************
SUBROUTINE qs_fgxc_analytic(rho0_struct, rho1_struct, xc_section, pw_pool, &
fxc_rho, fxc_tau, gxc_rho, gxc_tau, spinflip)
TYPE(qs_rho_type), POINTER :: rho0_struct, rho1_struct
TYPE(section_vals_type), POINTER :: xc_section
TYPE(pw_pool_type), POINTER :: pw_pool
TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: fxc_rho, fxc_tau, gxc_rho, gxc_tau
LOGICAL, INTENT(IN), OPTIONAL :: spinflip
CHARACTER(len=*), PARAMETER :: routineN = 'qs_fgxc_analytic'
INTEGER :: handle, ispin, nspins, spindim
INTEGER, DIMENSION(2, 3) :: bo
LOGICAL :: do_sf, lsd
REAL(KIND=dp) :: fac
TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho0_g, rho1_g
TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho0_r, rho1_r, tau0_r, tau1_r
TYPE(section_vals_type), POINTER :: xc_fun_section
TYPE(xc_derivative_set_type) :: deriv_set
TYPE(xc_rho_cflags_type) :: needs
TYPE(xc_rho_set_type) :: rho0_set, rho1_set
CALL timeset(routineN, handle)
! Only rho0 and rho1 should be associated
CPASSERT(.NOT. ASSOCIATED(fxc_rho))
CPASSERT(.NOT. ASSOCIATED(fxc_tau))
CPASSERT(.NOT. ASSOCIATED(gxc_rho))
CPASSERT(.NOT. ASSOCIATED(gxc_tau))
CPASSERT(ASSOCIATED(rho0_struct))
CPASSERT(ASSOCIATED(rho1_struct))
! Initialize parameters
do_sf = .FALSE.
IF (PRESENT(spinflip)) do_sf = spinflip
!
! Get the values on the gridpoints of the rho0 density
CALL qs_rho_get(rho0_struct, rho_r=rho0_r, rho_g=rho0_g, tau_r=tau0_r)
nspins = SIZE(rho0_r)
lsd = (nspins == 2)
!
IF (do_sf) THEN
spindim = 1
ELSE
spindim = nspins
END IF
!
fac = 0._dp
IF (nspins == 1) THEN
fac = 1.0_dp
END IF
! Read xc functional section and find out what the functional actually needs
xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
needs = xc_functionals_get_needs(xc_fun_section, lsd, .TRUE.)
! Create fields for the kernel and kernel derivative
ALLOCATE (fxc_rho(spindim), gxc_rho(nspins))
DO ispin = 1, spindim
CALL pw_pool%create_pw(fxc_rho(ispin))
CALL pw_zero(fxc_rho(ispin))
END DO
DO ispin = 1, nspins
CALL pw_pool%create_pw(gxc_rho(ispin))
CALL pw_zero(gxc_rho(ispin))
END DO
! Create fields for mGGA functionals. This implementation is not ready yet!
IF (needs%tau .OR. needs%tau_spin) THEN
IF (.NOT. ASSOCIATED(tau1_r)) &
CPABORT("Tau-dependent functionals requires allocated kinetic energy density grid")
ALLOCATE (fxc_tau(spindim), gxc_tau(nspins))
DO ispin = 1, spindim
CALL pw_pool%create_pw(fxc_tau(ispin))
CALL pw_zero(fxc_tau(ispin))
END DO
DO ispin = 1, nspins
CALL pw_pool%create_pw(gxc_tau(ispin))
CALL pw_zero(gxc_tau(ispin))
END DO
END IF
! Build rho0_set
! calculate the arguments needed by the functionals
! Needs
! deriv_set xc_derivative_set_type just declared
! rho0_set xc_rho_set_type just declared
! rho0_r pw_type calculated by qs_rho_get
! pw_pool given by the calling subroutine
! xc_section given by the calling subroutine
! tau0_r pw_type calculated by qs_rho_get
CALL xc_prep_3rd_deriv(deriv_set, rho0_set, rho0_r, pw_pool, xc_section, tau_r=tau0_r, &
do_sf=do_sf)
! Build rho1_set
! Get the values on the gridpoints of the rho1 density
CALL qs_rho_get(rho1_struct, rho_r=rho1_r, rho_g=rho1_g, tau_r=tau1_r)
bo = rho1_r(1)%pw_grid%bounds_local
! create the place where to store the argument for the functionals
! Needs
! rho1_set xc_rho_set_type just declared
! bo 2x3 integer matrix should have bounds_local or rho1_r
CALL xc_rho_set_create(rho1_set, bo, &
rho_cutoff=section_get_rval(xc_section, "DENSITY_CUTOFF"), &
drho_cutoff=section_get_rval(xc_section, "GRADIENT_CUTOFF"), &
tau_cutoff=section_get_rval(xc_section, "TAU_CUTOFF"))
! calculate the arguments needed by the functionals
! Needs
! rho1_set object created by xc_rho_set_create
! rho1_r,rho1_g,tau1_r pw_type values of rho1 in real space grid
! needs xc_rho_cflags_type defined through xc_functionals_get_needs
! pw_pool Given by the calling subroutine
CALL xc_rho_set_update(rho1_set, rho1_r, rho1_g, tau1_r, needs, &
section_get_ival(xc_section, "XC_GRID%XC_DERIV"), &
section_get_ival(xc_section, "XC_GRID%XC_SMOOTH_RHO"), &
pw_pool, spinflip=do_sf)
! Calculate exchange correlation kernel
! Needs
! fxc_rho, fxc_tau pw_type not associated
! deriv_set created and defined by xc_prep_3rd_deriv
! rho0_set xc_rho_set_type build by xc_prep_3rd_deriv
! rho1_set xc_rho_set_type build by xc_rho_set_create/update
! pw_pool needs to be given by the calling subroutine
! xc_section needs to be given by the calling subroutine
CALL xc_calc_2nd_deriv_analytical(fxc_rho, fxc_tau, deriv_set, rho0_set, rho1_set, pw_pool, &
xc_section, .FALSE., spinflip=do_sf, tddfpt_fac=fac)
! Calculate exchange correlation kernel derivative
CALL xc_calc_3rd_deriv_analytical(gxc_rho, gxc_tau, deriv_set, rho0_set, rho1_set, pw_pool, &
xc_section, spinflip=do_sf)
CALL xc_dset_release(deriv_set)
CALL xc_rho_set_release(rho0_set)
CALL xc_rho_set_release(rho1_set)
CALL timestop(handle)
END SUBROUTINE qs_fgxc_analytic
! **************************************************************************************************
!> \brief ...
!> \param ks_env ...
@ -278,9 +448,10 @@ CONTAINS
!> \param fxc_tau ...
!> \param gxc_rho ...
!> \param gxc_tau ...
!> \param spinflip ...
! **************************************************************************************************
SUBROUTINE qs_fgxc_gdiff(ks_env, rho0_struct, rho1_struct, xc_section, accuracy, epsrho, &
is_triplet, fxc_rho, fxc_tau, gxc_rho, gxc_tau)
is_triplet, fxc_rho, fxc_tau, gxc_rho, gxc_tau, spinflip)
TYPE(qs_ks_env_type), POINTER :: ks_env
TYPE(qs_rho_type), POINTER :: rho0_struct, rho1_struct
@ -289,16 +460,19 @@ CONTAINS
REAL(KIND=dp), INTENT(IN) :: epsrho
LOGICAL, INTENT(IN) :: is_triplet
TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: fxc_rho, fxc_tau, gxc_rho, gxc_tau
LOGICAL, OPTIONAL :: spinflip
CHARACTER(len=*), PARAMETER :: routineN = 'qs_fgxc_gdiff'
INTEGER :: handle, ispin, istep, nspins, nstep
LOGICAL :: do_sf
REAL(KIND=dp) :: alpha, beta, exc, oeps1
REAL(KIND=dp), DIMENSION(-4:4) :: ak
TYPE(dft_control_type), POINTER :: dft_control
TYPE(pw_env_type), POINTER :: pw_env
TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_tau_rspace, vxc00
TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho1_r, tau1_r, v_tau_rspace, vxc00, &
vxc00b
TYPE(qs_rho_type), POINTER :: rhoin
CALL timeset(routineN, handle)
@ -310,6 +484,9 @@ CONTAINS
CPASSERT(ASSOCIATED(rho0_struct))
CPASSERT(ASSOCIATED(rho1_struct))
do_sf = .FALSE.
IF (PRESENT(spinflip)) do_sf = spinflip
ak = 0.0_dp
SELECT CASE (accuracy)
CASE (:4)
@ -330,8 +507,14 @@ CONTAINS
nspins = dft_control%nspins
exc = 0.0_dp
CALL qs_fxc_fdiff(ks_env, rho0_struct, rho1_struct, xc_section, accuracy, is_triplet, &
fxc_rho, fxc_tau)
IF (do_sf) THEN
CALL qs_rho_get(rho1_struct, rho_r=rho1_r, tau_r=tau1_r)
CALL qs_fxc_analytic(rho0_struct, rho1_r, tau1_r, xc_section, &
auxbas_pw_pool, is_triplet, fxc_rho, fxc_tau, spinflip=do_sf)
ELSE
CALL qs_fxc_fdiff(ks_env, rho0_struct, rho1_struct, xc_section, accuracy, is_triplet, &
fxc_rho, fxc_tau)
END IF
DO istep = -nstep, nstep
@ -341,12 +524,25 @@ CONTAINS
NULLIFY (rhoin)
ALLOCATE (rhoin)
CALL qs_rho_create(rhoin)
NULLIFY (vxc00, v_tau_rspace)
NULLIFY (vxc00, vxc00b, v_tau_rspace)
CALL qs_rho_copy(rho0_struct, rhoin, auxbas_pw_pool, nspins)
CALL qs_rho_scale_and_add(rhoin, rho1_struct, alpha, beta)
CALL qs_fxc_fdiff(ks_env=ks_env, rho0_struct=rhoin, rho1_struct=rho1_struct, &
xc_section=xc_section, accuracy=accuracy, is_triplet=is_triplet, &
fxc_rho=vxc00, fxc_tau=v_tau_rspace)
IF (do_sf) THEN
! variation in alpha density
CALL qs_fxc_analytic(rhoin, rho1_r, tau1_r, &
xc_section, auxbas_pw_pool, is_triplet, &
vxc00, v_tau_rspace, spinflip=do_sf)
! variation in beta density
CALL qs_rho_copy(rho0_struct, rhoin, auxbas_pw_pool, nspins)
CALL qs_rho_scale_and_add_b(rhoin, rho1_struct, alpha, beta)
CALL qs_fxc_analytic(rhoin, rho1_r, tau1_r, &
xc_section, auxbas_pw_pool, is_triplet, &
vxc00b, v_tau_rspace, spinflip=do_sf)
ELSE
CALL qs_fxc_fdiff(ks_env=ks_env, rho0_struct=rhoin, rho1_struct=rho1_struct, &
xc_section=xc_section, accuracy=accuracy, is_triplet=is_triplet, &
fxc_rho=vxc00, fxc_tau=v_tau_rspace)
END IF
CALL qs_rho_release(rhoin)
DEALLOCATE (rhoin)
IF (.NOT. ASSOCIATED(gxc_rho)) THEN
@ -356,13 +552,22 @@ CONTAINS
CALL pw_zero(gxc_rho(ispin))
END DO
END IF
DO ispin = 1, nspins
CALL pw_axpy(vxc00(ispin), gxc_rho(ispin), ak(istep))
END DO
IF (do_sf) THEN
CALL pw_axpy(vxc00(1), gxc_rho(1), ak(istep))
CALL pw_axpy(vxc00b(1), gxc_rho(2), ak(istep))
ELSE
DO ispin = 1, nspins
CALL pw_axpy(vxc00(ispin), gxc_rho(ispin), ak(istep))
END DO
END IF
DO ispin = 1, SIZE(vxc00)
CALL auxbas_pw_pool%give_back_pw(vxc00(ispin))
END DO
DEALLOCATE (vxc00)
IF (ASSOCIATED(vxc00b)) THEN
CALL auxbas_pw_pool%give_back_pw(vxc00b(1))
DEALLOCATE (vxc00b)
END IF
IF (ASSOCIATED(v_tau_rspace)) THEN
IF (.NOT. ASSOCIATED(gxc_tau)) THEN
ALLOCATE (gxc_tau(nspins))

10
src/qs_integrate_potential_product.F
View file

@ -541,10 +541,14 @@ CONTAINS
! **************************************************************************************************
!> \brief computes matrix elements corresponding to a given potential
!> \param v_rspace ...
!> \param hmat ...
!> If compute_tau
!> hmat(\mu,\nu) = \int v_rspace(r)\nabla\mu(r)\cdot\nabla\nu(r) dr
!> else
!> hmat(\mu,\nu) = \int v_rspace(r)\mu(r)\nu(r) dr
!> \param v_rspace Potential in real space
!> \param hmat Matrix to write the result
!> \param hmat_kp ...
!> \param pmat ...
!> \param pmat Input density matrix, not used in the calculation but the screening. See note below.
!> \param pmat_kp ...
!> \param qs_env ...
!> \param calculate_forces ...

120
src/qs_rho_methods.F
View file

@ -70,7 +70,8 @@ MODULE qs_rho_methods
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_rho_methods'
PUBLIC :: qs_rho_update_rho, qs_rho_update_tddfpt, &
qs_rho_rebuild, qs_rho_copy, qs_rho_scale_and_add
qs_rho_rebuild, qs_rho_copy, qs_rho_scale_and_add, &
qs_rho_scale_and_add_b
PUBLIC :: duplicate_rho_type, allocate_rho_ao_imag_from_real
CONTAINS
@ -1057,6 +1058,123 @@ CONTAINS
END SUBROUTINE qs_rho_copy
! **************************************************************************************************
!> \brief rhoa(2) = alpha*rhoa(2)+beta*rhob(1)
!> \param rhoa ...
!> \param rhob ...
!> \param alpha ...
!> \param beta ...
! **************************************************************************************************
SUBROUTINE qs_rho_scale_and_add_b(rhoa, rhob, alpha, beta)
TYPE(qs_rho_type), INTENT(IN) :: rhoa, rhob
REAL(KIND=dp), INTENT(IN) :: alpha, beta
CHARACTER(len=*), PARAMETER :: routineN = 'qs_rho_scale_and_add_b'
INTEGER :: handle
REAL(KIND=dp), DIMENSION(:), POINTER :: tot_rho_g_a, tot_rho_g_b, tot_rho_r_a, &
tot_rho_r_b
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao_a, rho_ao_b, rho_ao_im_a, &
rho_ao_im_b
TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g_a, rho_g_b, tau_g_a, tau_g_b
TYPE(pw_c1d_gs_type), DIMENSION(:, :), POINTER :: drho_g_a, drho_g_b
TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r_a, rho_r_b, tau_r_a, tau_r_b
TYPE(pw_r3d_rs_type), DIMENSION(:, :), POINTER :: drho_r_a, drho_r_b
TYPE(pw_r3d_rs_type), POINTER :: rho_r_sccs_a, rho_r_sccs_b
CALL timeset(routineN, handle)
NULLIFY (rho_ao_a, rho_ao_im_a, rho_r_a, rho_g_a, drho_r_a, &
drho_g_a, tau_r_a, tau_g_a, tot_rho_r_a, tot_rho_g_a, rho_r_sccs_a)
CALL qs_rho_get(rhoa, &
rho_ao=rho_ao_a, &
rho_ao_im=rho_ao_im_a, &
rho_r=rho_r_a, &
rho_g=rho_g_a, &
drho_r=drho_r_a, &
drho_g=drho_g_a, &
tau_r=tau_r_a, &
tau_g=tau_g_a, &
tot_rho_r=tot_rho_r_a, &
tot_rho_g=tot_rho_g_a, &
rho_r_sccs=rho_r_sccs_a)
NULLIFY (rho_ao_b, rho_ao_im_b, rho_r_b, rho_g_b, drho_r_b, &
drho_g_b, tau_r_b, tau_g_b, tot_rho_r_b, tot_rho_g_b, rho_r_sccs_b)
CALL qs_rho_get(rhob, &
rho_ao=rho_ao_b, &
rho_ao_im=rho_ao_im_b, &
rho_r=rho_r_b, &
rho_g=rho_g_b, &
drho_r=drho_r_b, &
drho_g=drho_g_b, &
tau_r=tau_r_b, &
tau_g=tau_g_b, &
tot_rho_r=tot_rho_r_b, &
tot_rho_g=tot_rho_g_b, &
rho_r_sccs=rho_r_sccs_b)
! rho_ao
IF (ASSOCIATED(rho_ao_a) .AND. ASSOCIATED(rho_ao_b)) THEN
CALL dbcsr_add(rho_ao_a(2)%matrix, rho_ao_b(1)%matrix, alpha, beta)
END IF
! rho_ao_im
IF (ASSOCIATED(rho_ao_im_a) .AND. ASSOCIATED(rho_ao_im_b)) THEN
CALL dbcsr_add(rho_ao_im_a(2)%matrix, rho_ao_im_b(1)%matrix, alpha, beta)
END IF
! rho_r
IF (ASSOCIATED(rho_r_a) .AND. ASSOCIATED(rho_r_b)) THEN
CALL pw_axpy(rho_r_b(1), rho_r_a(2), beta, alpha)
END IF
! rho_g
IF (ASSOCIATED(rho_g_a) .AND. ASSOCIATED(rho_g_b)) THEN
CALL pw_axpy(rho_g_b(1), rho_g_a(2), beta, alpha)
END IF
! SCCS
IF (ASSOCIATED(rho_r_sccs_a) .AND. ASSOCIATED(rho_r_sccs_b)) THEN
CALL pw_axpy(rho_r_sccs_b, rho_r_sccs_a, beta, alpha)
END IF
! drho_r
IF (ASSOCIATED(drho_r_a) .AND. ASSOCIATED(drho_r_b)) THEN
CPASSERT(ASSOCIATED(drho_r_a) .AND. ASSOCIATED(drho_r_b)) ! not implemented
END IF
! drho_g
IF (ASSOCIATED(drho_g_a) .AND. ASSOCIATED(drho_g_b)) THEN
CPASSERT(ASSOCIATED(drho_g_a) .AND. ASSOCIATED(drho_g_b)) ! not implemented
END IF
! tau_r
IF (ASSOCIATED(tau_r_a) .AND. ASSOCIATED(tau_r_b)) THEN
CALL pw_axpy(tau_r_b(1), tau_r_a(2), beta, alpha)
END IF
! tau_g
IF (ASSOCIATED(tau_g_a) .AND. ASSOCIATED(tau_g_b)) THEN
CALL pw_axpy(tau_g_b(1), tau_g_a(2), beta, alpha)
END IF
! tot_rho_r
IF (ASSOCIATED(tot_rho_r_a) .AND. ASSOCIATED(tot_rho_r_b)) THEN
tot_rho_r_a(2) = alpha*tot_rho_r_a(2) + beta*tot_rho_r_b(1)
END IF
! tot_rho_g
IF (ASSOCIATED(tot_rho_g_a) .AND. ASSOCIATED(tot_rho_g_b)) THEN
tot_rho_g_a(2) = alpha*tot_rho_g_a(2) + beta*tot_rho_g_b(1)
END IF
CALL timestop(handle)
END SUBROUTINE qs_rho_scale_and_add_b
! **************************************************************************************************
!> \brief rhoa = alpha*rhoa+beta*rhob
!> \param rhoa ...

137
src/qs_tddfpt2_eigensolver.F
View file

@ -27,7 +27,8 @@ MODULE qs_tddfpt2_eigensolver
cp_fm_release, cp_fm_set_all, cp_fm_set_submatrix, cp_fm_to_fm, cp_fm_type
USE cp_log_handling, ONLY: cp_logger_type
USE cp_output_handling, ONLY: cp_iterate
USE input_constants, ONLY: tddfpt_kernel_full,&
USE input_constants, ONLY: no_sf_tddfpt,&
tddfpt_kernel_full,&
tddfpt_kernel_none,&
tddfpt_kernel_stda
USE input_section_types, ONLY: section_vals_type
@ -113,7 +114,7 @@ CONTAINS
CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_orthogonalize_psi1_psi0'
INTEGER :: handle, ispin, ivect, nactive, nao, &
nspins, nvects
nspins, nvects, spin
TYPE(cp_fm_struct_type), POINTER :: matrix_struct
TYPE(cp_fm_type) :: evortho
TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos
@ -127,12 +128,17 @@ CONTAINS
IF (nvects > 0) THEN
IF (.NOT. tddfpt_control%do_smearing) THEN
DO ispin = 1, nspins
IF (tddfpt_control%spinflip == no_sf_tddfpt) THEN
spin = ispin
ELSE
spin = 2
END IF
CALL cp_fm_get_info(matrix=evects(ispin, 1), matrix_struct=matrix_struct, &
nrow_global=nao, ncol_global=nactive)
CALL cp_fm_create(evortho, matrix_struct)
DO ivect = 1, nvects
! evortho: C0 * C0^T * S * C1 == (S * C0 * C0^T)^T * C1
CALL parallel_gemm('T', 'N', nao, nactive, nao, 1.0_dp, S_C0_C0T(ispin), &
CALL parallel_gemm('T', 'N', nao, nactive, nao, 1.0_dp, S_C0_C0T(spin), &
evects(ispin, ivect), 0.0_dp, evortho)
CALL cp_fm_scale_and_add(1.0_dp, evects(ispin, ivect), -1.0_dp, evortho)
END DO
@ -172,22 +178,24 @@ CONTAINS
!> for each spin in atomic basis set, and S is the corresponding overlap matrix
!> \param max_norm the largest possible overlap between the ground state and
!> excited state wave functions
!> \param spinflip ...
!> \return true if trial vectors are non-orthogonal to occupied molecular orbitals
!> \par History
!> * 07.2016 created [Sergey Chulkov]
!> * 05.2019 use temporary work matrices [JHU]
! **************************************************************************************************
FUNCTION tddfpt_is_nonorthogonal_psi1_psi0(evects, S_C0, max_norm) &
FUNCTION tddfpt_is_nonorthogonal_psi1_psi0(evects, S_C0, max_norm, spinflip) &
RESULT(is_nonortho)
TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in) :: evects
TYPE(cp_fm_type), DIMENSION(:), INTENT(in) :: S_C0
REAL(kind=dp), INTENT(in) :: max_norm
INTEGER, INTENT(in) :: spinflip
LOGICAL :: is_nonortho
CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_is_nonorthogonal_psi1_psi0'
INTEGER :: handle, ispin, ivect, nactive, nao, &
nocc, nspins, nvects
nocc, nspins, nvects, spin
REAL(kind=dp) :: maxabs_val
TYPE(cp_fm_struct_type), POINTER :: matrix_struct, matrix_struct_tmp
TYPE(cp_fm_type) :: aortho
@ -200,7 +208,14 @@ CONTAINS
is_nonortho = .FALSE.
loop: DO ispin = 1, nspins
CALL cp_fm_get_info(matrix=S_C0(ispin), ncol_global=nocc)
IF (spinflip == no_sf_tddfpt) THEN
spin = ispin
ELSE
spin = 2
END IF
CALL cp_fm_get_info(matrix=S_C0(spin), ncol_global=nocc)
CALL cp_fm_get_info(matrix=evects(ispin, 1), matrix_struct=matrix_struct, &
nrow_global=nao, ncol_global=nactive)
CALL cp_fm_struct_create(matrix_struct_tmp, nrow_global=nocc, &
@ -208,8 +223,8 @@ CONTAINS
CALL cp_fm_create(aortho, matrix_struct_tmp)
CALL cp_fm_struct_release(matrix_struct_tmp)
DO ivect = 1, nvects
! aortho = S_0^T * S * C_1
CALL parallel_gemm('T', 'N', nocc, nactive, nao, 1.0_dp, S_C0(ispin), &
! aortho = S_C0^T * S * C_1
CALL parallel_gemm('T', 'N', nocc, nactive, nao, 1.0_dp, S_C0(spin), &
evects(ispin, ivect), 0.0_dp, aortho)
CALL cp_fm_maxabsval(aortho, maxabs_val)
is_nonortho = maxabs_val > max_norm
@ -290,7 +305,7 @@ CONTAINS
END DO
DO jvect = nvects_old + 1, nvects_total
! <psi1_i | psi1_j>
! Orthogonalization <psi1_i | psi1_j>
DO ivect = 1, jvect - 1
CALL cp_fm_trace(evects(:, jvect), S_evects(:, ivect), weights(1:nspins), accurate=.FALSE.)
norm = SUM(weights(1:nspins))
@ -300,7 +315,7 @@ CONTAINS
END DO
END DO
! <psi1_j | psi1_j>
! Normalization <psi1_j | psi1_j> = 1
DO ispin = 1, nspins
CALL cp_dbcsr_sm_fm_multiply(matrix_s, evects(ispin, jvect), S_evects(ispin, jvect), &
ncol=nactive(ispin), alpha=1.0_dp, beta=0.0_dp)
@ -367,7 +382,7 @@ CONTAINS
CALL timeset(routineN, handle)
nspins = SIZE(evects, 1)
nspins = SIZE(gs_mos, 1)
nvects = SIZE(evects, 2)
do_hfx = tddfpt_control%do_hfx
do_admm = tddfpt_control%do_admm
@ -381,9 +396,9 @@ CONTAINS
IF (debug_this_module) THEN
CPASSERT(nspins > 0)
CPASSERT(SIZE(Aop_evects, 1) == nspins)
CPASSERT(SIZE(Aop_evects, 1) == SIZE(evects, 1))
CPASSERT(SIZE(S_evects, 1) == SIZE(evects, 1))
CPASSERT(SIZE(Aop_evects, 2) == nvects)
CPASSERT(SIZE(S_evects, 1) == nspins)
CPASSERT(SIZE(S_evects, 2) == nvects)
CPASSERT(SIZE(gs_mos) == nspins)
END IF
@ -396,19 +411,19 @@ CONTAINS
CALL cp_fm_get_info(evects(1, 1), para_env=para_env)
IF (ALLOCATED(work_matrices%evects_sub)) THEN
DO ivect = 1, nvects
DO ispin = 1, nspins
DO ispin = 1, SIZE(evects, 1)
ASSOCIATE (evect => evects(ispin, ivect), work_matrix => work_matrices%evects_sub(ispin, ivect))
IF (ASSOCIATED(evect%matrix_struct)) THEN
IF (ASSOCIATED(work_matrix%matrix_struct)) THEN
CALL cp_fm_copy_general(evect, work_matrix, para_env)
ELSE
CALL cp_fm_copy_general(evect, fm_dummy, para_env)
END IF
ELSE IF (ASSOCIATED(work_matrix%matrix_struct)) THEN
CALL cp_fm_copy_general(fm_dummy, work_matrix, para_env)
ELSE
CALL cp_fm_copy_general(fm_dummy, fm_dummy, para_env)
END IF
IF (ASSOCIATED(evect%matrix_struct)) THEN
IF (ASSOCIATED(work_matrix%matrix_struct)) THEN
CALL cp_fm_copy_general(evect, work_matrix, para_env)
ELSE
CALL cp_fm_copy_general(evect, fm_dummy, para_env)
END IF
ELSE IF (ASSOCIATED(work_matrix%matrix_struct)) THEN
CALL cp_fm_copy_general(fm_dummy, work_matrix, para_env)
ELSE
CALL cp_fm_copy_general(fm_dummy, fm_dummy, para_env)
END IF
END ASSOCIATE
END DO
END DO
@ -427,7 +442,7 @@ CONTAINS
ELSE IF (tddfpt_control%kernel == tddfpt_kernel_none) THEN
! No kernel
DO ivect = 1, nvects
DO ispin = 1, nspins
DO ispin = 1, SIZE(evects, 1)
CALL cp_fm_set_all(Aop_evects(ispin, ivect), 0.0_dp)
END DO
END DO
@ -437,20 +452,20 @@ CONTAINS
IF (ALLOCATED(work_matrices%evects_sub)) THEN
DO ivect = 1, nvects
DO ispin = 1, nspins
DO ispin = 1, SIZE(evects, 1)
ASSOCIATE (Aop_evect => Aop_evects(ispin, ivect), &
work_matrix => work_matrices%Aop_evects_sub(ispin, ivect))
IF (ASSOCIATED(Aop_evect%matrix_struct)) THEN
IF (ASSOCIATED(work_matrix%matrix_struct)) THEN
CALL cp_fm_copy_general(work_matrix, Aop_evect, para_env)
ELSE
CALL cp_fm_copy_general(fm_dummy, Aop_evect, para_env)
END IF
ELSE IF (ASSOCIATED(work_matrix%matrix_struct)) THEN
CALL cp_fm_copy_general(work_matrix, fm_dummy, para_env)
ELSE
CALL cp_fm_copy_general(fm_dummy, fm_dummy, para_env)
END IF
IF (ASSOCIATED(Aop_evect%matrix_struct)) THEN
IF (ASSOCIATED(work_matrix%matrix_struct)) THEN
CALL cp_fm_copy_general(work_matrix, Aop_evect, para_env)
ELSE
CALL cp_fm_copy_general(fm_dummy, Aop_evect, para_env)
END IF
ELSE IF (ASSOCIATED(work_matrix%matrix_struct)) THEN
CALL cp_fm_copy_general(work_matrix, fm_dummy, para_env)
ELSE
CALL cp_fm_copy_general(fm_dummy, fm_dummy, para_env)
END IF
END ASSOCIATE
END DO
END DO
@ -458,7 +473,7 @@ CONTAINS
! orbital energy difference term
CALL tddfpt_apply_energy_diff(Aop_evects=Aop_evects, evects=evects, S_evects=S_evects, &
gs_mos=gs_mos, matrix_ks=matrix_ks)
gs_mos=gs_mos, matrix_ks=matrix_ks, tddfpt_control=tddfpt_control)
! if smeared occupation, then add aCCSX here
IF (tddfpt_control%do_smearing) THEN
@ -547,7 +562,7 @@ CONTAINS
!> \param evals Ritz eigenvalues (initialised on exit)
!> \param krylov_vects Krylov's vectors
!> \param Aop_krylov action of TDDFPT operator on Krylov's vectors
!> \param Atilde ...
!> \param Atilde TDDFPT matrix projected into the Krylov's vectors subspace
!> \param nkvo ...
!> \param nkvn ...
!> \par History
@ -653,25 +668,27 @@ CONTAINS
!> \param Aop_ritz approximate action of TDDFPT operator on Ritz vectors
!> \param gs_mos molecular orbitals optimised for the ground state
!> \param matrix_s overlap matrix
!> \param tddfpt_control ...
!> \par History
!> * 06.2016 created [Sergey Chulkov]
!> * 03.2017 refactored to achieve significant performance gain [Sergey Chulkov]
! **************************************************************************************************
SUBROUTINE tddfpt_compute_residual_vects(residual_vects, evals, ritz_vects, Aop_ritz, gs_mos, &
matrix_s)
matrix_s, tddfpt_control)
TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in) :: residual_vects
REAL(kind=dp), DIMENSION(:), INTENT(in) :: evals
TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in) :: ritz_vects, Aop_ritz
TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
INTENT(in) :: gs_mos
TYPE(dbcsr_type), POINTER :: matrix_s
TYPE(tddfpt2_control_type), POINTER :: tddfpt_control
CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_compute_residual_vects'
REAL(kind=dp), PARAMETER :: eref_scale = 0.99_dp, threshold = 16.0_dp*EPSILON(1.0_dp)
INTEGER :: handle, icol_local, irow_local, irv, &
ispin, nao, ncols_local, nrows_local, &
nrvs, nspins
nrvs, nspins, spin2, spinflip
INTEGER, DIMENSION(:), POINTER :: col_indices_local, row_indices_local
INTEGER, DIMENSION(maxspins) :: nactive, nmo_virt
REAL(kind=dp) :: e_occ_plus_lambda, eref, lambda
@ -684,18 +701,24 @@ CONTAINS
nspins = SIZE(residual_vects, 1)
nrvs = SIZE(residual_vects, 2)
spinflip = tddfpt_control%spinflip
IF (nrvs > 0) THEN
CALL dbcsr_get_info(matrix_s, nfullrows_total=nao)
ALLOCATE (awork(nspins), vomat(nspins))
DO ispin = 1, nspins
nmo_virt(ispin) = SIZE(gs_mos(ispin)%evals_virt)
IF (spinflip == no_sf_tddfpt) THEN
spin2 = ispin
ELSE
spin2 = 2
END IF
nmo_virt(spin2) = SIZE(gs_mos(spin2)%evals_virt)
!
CALL cp_fm_get_info(matrix=ritz_vects(ispin, 1), matrix_struct=ao_mo_struct, &
ncol_global=nactive(ispin))
CALL cp_fm_create(awork(ispin), ao_mo_struct)
!
CALL cp_fm_struct_create(virt_mo_struct, nrow_global=nmo_virt(ispin), &
CALL cp_fm_struct_create(virt_mo_struct, nrow_global=nmo_virt(spin2), &
ncol_global=nactive(ispin), template_fmstruct=ao_mo_struct)
CALL cp_fm_create(vomat(ispin), virt_mo_struct)
CALL cp_fm_struct_release(virt_mo_struct)
@ -706,6 +729,11 @@ CONTAINS
lambda = evals(irv)
DO ispin = 1, nspins
IF (spinflip == no_sf_tddfpt) THEN
spin2 = ispin
ELSE
spin2 = 2
END IF
CALL cp_fm_get_info(vomat(ispin), nrow_local=nrows_local, &
ncol_local=ncols_local, row_indices=row_indices_local, &
col_indices=col_indices_local, local_data=weights_ldata)
@ -715,14 +743,15 @@ CONTAINS
ncol=nactive(ispin), alpha=-lambda, beta=0.0_dp)
CALL cp_fm_scale_and_add(1.0_dp, awork(ispin), 1.0_dp, Aop_ritz(ispin, irv))
!
CALL parallel_gemm('T', 'N', nmo_virt(ispin), nactive(ispin), nao, 1.0_dp, gs_mos(ispin)%mos_virt, &
CALL parallel_gemm('T', 'N', nmo_virt(spin2), nactive(ispin), nao, 1.0_dp, gs_mos(spin2)%mos_virt, &
awork(ispin), 0.0_dp, vomat(ispin))
! Apply Davidson preconditioner to the residue vectors vomat to obtain new directions
DO icol_local = 1, ncols_local
e_occ_plus_lambda = gs_mos(ispin)%evals_occ(col_indices_local(icol_local)) + lambda
DO irow_local = 1, nrows_local
eref = gs_mos(ispin)%evals_virt(row_indices_local(irow_local)) - e_occ_plus_lambda
eref = gs_mos(spin2)%evals_virt(row_indices_local(irow_local)) - e_occ_plus_lambda
! eref = e_virt - e_occ - lambda = e_virt - e_occ - (eref_scale*lambda + (1-eref_scale)*lambda);
! eref_new = e_virt - e_occ - eref_scale*lambda = eref + (1 - eref_scale)*lambda
@ -733,7 +762,7 @@ CONTAINS
END DO
END DO
CALL parallel_gemm('N', 'N', nao, nactive(ispin), nmo_virt(ispin), 1.0_dp, gs_mos(ispin)%mos_virt, &
CALL parallel_gemm('N', 'N', nao, nactive(ispin), nmo_virt(spin2), 1.0_dp, gs_mos(spin2)%mos_virt, &
vomat(ispin), 0.0_dp, residual_vects(ispin, irv))
END DO
END DO
@ -805,9 +834,9 @@ CONTAINS
CALL timeset(routineN, handle)
nspins = SIZE(gs_mos)
nspins = SIZE(evects, 1)
nstates = tddfpt_control%nstates
nstates_total = tddfpt_total_number_of_states(gs_mos)
nstates_total = tddfpt_total_number_of_states(tddfpt_control, gs_mos)
IF (debug_this_module) THEN
CPASSERT(SIZE(evects, 1) == nspins)
@ -818,9 +847,7 @@ CONTAINS
CALL get_qs_env(qs_env, matrix_s=matrix_s)
! adjust the number of Krylov vectors
max_krylov_vects = tddfpt_control%nkvs
IF (max_krylov_vects < nstates) max_krylov_vects = nstates
IF (INT(max_krylov_vects, kind=int_8) > nstates_total) max_krylov_vects = INT(nstates_total)
max_krylov_vects = MIN(MAX(tddfpt_control%nkvs, nstates), INT(nstates_total))
ALLOCATE (Aop_ritz(nspins, nstates))
DO istate = 1, nstates
@ -856,6 +883,7 @@ CONTAINS
! davidson iteration
CALL cp_iterate(logger%iter_info, iter_nr_out=iter)
! Matrix-vector operations
CALL tddfpt_compute_Aop_evects(Aop_evects=Aop_krylov(:, nvects_exist + 1:nvects_exist + nvects_new), &
evects=krylov_vects(:, nvects_exist + 1:nvects_exist + nvects_new), &
S_evects=S_krylov(:, nvects_exist + 1:nvects_exist + nvects_new), &
@ -898,7 +926,7 @@ CONTAINS
CALL tddfpt_compute_residual_vects(residual_vects=krylov_vects(:, nvects_exist + 1:nvects_exist + nvects_new), &
evals=evals_last(1:nvects_new), &
ritz_vects=evects(:, 1:nvects_new), Aop_ritz=Aop_ritz(:, 1:nvects_new), &
gs_mos=gs_mos, matrix_s=matrix_s(1)%matrix)
gs_mos=gs_mos, matrix_s=matrix_s(1)%matrix, tddfpt_control=tddfpt_control)
CALL tddfpt_orthogonalize_psi1_psi0(krylov_vects(:, nvects_exist + 1:nvects_exist + nvects_new), &
work_matrices%S_C0_C0T, qs_env, &
@ -908,7 +936,8 @@ CONTAINS
S_krylov(:, 1:nvects_exist + nvects_new), matrix_s(1)%matrix)
is_nonortho = tddfpt_is_nonorthogonal_psi1_psi0(krylov_vects(:, nvects_exist + 1:nvects_exist + nvects_new), &
work_matrices%S_C0, tddfpt_control%orthogonal_eps)
work_matrices%S_C0, tddfpt_control%orthogonal_eps, &
tddfpt_control%spinflip)
ELSE
! convergence or the maximum number of Krylov vectors have been achieved
nvects_new = 0

107
src/qs_tddfpt2_fhxc.F
View file

@ -20,7 +20,10 @@ MODULE qs_tddfpt2_fhxc
cp_fm_get_info,&
cp_fm_release,&
cp_fm_type
USE input_constants, ONLY: do_admm_aux_exch_func_none
USE input_constants, ONLY: do_admm_aux_exch_func_none,&
no_sf_tddfpt,&
tddfpt_sf_col,&
tddfpt_sf_noncol
USE kinds, ONLY: default_string_length,&
dp
USE lri_environment_types, ONLY: lri_kind_type
@ -111,10 +114,11 @@ CONTAINS
CHARACTER(LEN=default_string_length) :: basis_type
INTEGER :: handle, ikind, ispin, ivect, nao, &
nao_aux, nkind, nspins, nvects
nao_aux, nkind, nspins, nvects, &
spinflip
INTEGER, DIMENSION(:), POINTER :: blk_sizes
INTEGER, DIMENSION(maxspins) :: nactive
LOGICAL :: gapw, gapw_xc
LOGICAL :: do_noncol, gapw, gapw_xc
TYPE(admm_type), POINTER :: admm_env
TYPE(cp_fm_type) :: work_aux_orb, work_orb_orb
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: A_xc_munu_sub, rho_ia_ao, &
@ -142,6 +146,9 @@ CONTAINS
gapw = dft_control%qs_control%gapw
gapw_xc = dft_control%qs_control%gapw_xc
spinflip = dft_control%tddfpt2_control%spinflip
do_noncol = spinflip .EQ. tddfpt_sf_noncol
CALL cp_fm_get_info(evects(1, 1), nrow_global=nao)
DO ispin = 1, nspins
@ -158,6 +165,8 @@ CONTAINS
END IF
DO ivect = 1, nvects
! Transform TDDFT vectors to AO space and store them into rho_ia_ao
IF (ALLOCATED(work_matrices%evects_sub)) THEN
IF (ASSOCIATED(work_matrices%evects_sub(1, ivect)%matrix_struct)) THEN
DO ispin = 1, nspins
@ -230,46 +239,54 @@ CONTAINS
CALL pw_zero(work_matrices%A_ia_rspace_sub(ispin))
END DO
! C_x d^{2}E_{x}^{DFT}[\rho] / d\rho^2
! + C_{HF} d^{2}E_{x, ADMM}^{DFT}[\rho] / d\rho^2 in case of ADMM calculation
IF (gapw_xc) THEN
IF (kernel_env%do_exck) THEN
CPABORT("NYA")
! Skip kernel if collinear xc-kernel for spin-flip is requested
IF (spinflip /= tddfpt_sf_col) THEN
! C_x d^{2}E_{x}^{DFT}[\rho] / d\rho^2
! + C_{HF} d^{2}E_{x, ADMM}^{DFT}[\rho] / d\rho^2 in case of ADMM calculation
IF (gapw_xc) THEN
IF (kernel_env%do_exck) THEN
CPABORT("NYA")
ELSE
CALL tddfpt_apply_xc(A_ia_rspace=work_matrices%A_ia_rspace_sub, kernel_env=kernel_env, &
rho_ia_struct=work_matrices%rho_xc_struct_sub, &
is_rks_triplets=is_rks_triplets, pw_env=sub_env%pw_env, &
work_v_xc=work_matrices%wpw_rspace_sub, &
work_v_xc_tau=work_matrices%wpw_tau_rspace_sub, &
spinflip=spinflip)
END IF
DO ispin = 1, nspins
CALL pw_scale(work_matrices%A_ia_rspace_sub(ispin), &
work_matrices%A_ia_rspace_sub(ispin)%pw_grid%dvol)
CALL integrate_v_rspace(v_rspace=work_matrices%A_ia_rspace_sub(ispin), &
hmat=work_matrices%A_ia_munu_sub(ispin), &
qs_env=qs_env, calculate_forces=.FALSE., gapw=gapw_xc, &
pw_env_external=sub_env%pw_env, &
task_list_external=sub_env%task_list_orb_soft)
CALL pw_zero(work_matrices%A_ia_rspace_sub(ispin))
END DO
ELSE
CALL tddfpt_apply_xc(A_ia_rspace=work_matrices%A_ia_rspace_sub, kernel_env=kernel_env, &
rho_ia_struct=work_matrices%rho_xc_struct_sub, &
is_rks_triplets=is_rks_triplets, pw_env=sub_env%pw_env, &
work_v_xc=work_matrices%wpw_rspace_sub, &
work_v_xc_tau=work_matrices%wpw_tau_rspace_sub)
IF (kernel_env%do_exck) THEN
CALL tddfpt_apply_xc_potential(work_matrices%A_ia_rspace_sub, work_matrices%fxc_rspace_sub, &
work_matrices%rho_orb_struct_sub, is_rks_triplets)
ELSE
CALL tddfpt_apply_xc(A_ia_rspace=work_matrices%A_ia_rspace_sub, kernel_env=kernel_env, &
rho_ia_struct=work_matrices%rho_orb_struct_sub, &
is_rks_triplets=is_rks_triplets, pw_env=sub_env%pw_env, &
work_v_xc=work_matrices%wpw_rspace_sub, &
work_v_xc_tau=work_matrices%wpw_tau_rspace_sub, &
spinflip=spinflip)
END IF
END IF
DO ispin = 1, nspins
CALL pw_scale(work_matrices%A_ia_rspace_sub(ispin), &
work_matrices%A_ia_rspace_sub(ispin)%pw_grid%dvol)
CALL integrate_v_rspace(v_rspace=work_matrices%A_ia_rspace_sub(ispin), &
hmat=work_matrices%A_ia_munu_sub(ispin), &
qs_env=qs_env, calculate_forces=.FALSE., gapw=gapw_xc, &
pw_env_external=sub_env%pw_env, &
task_list_external=sub_env%task_list_orb_soft)
CALL pw_zero(work_matrices%A_ia_rspace_sub(ispin))
END DO
ELSE
IF (kernel_env%do_exck) THEN
CALL tddfpt_apply_xc_potential(work_matrices%A_ia_rspace_sub, work_matrices%fxc_rspace_sub, &
work_matrices%rho_orb_struct_sub, is_rks_triplets)
ELSE
CALL tddfpt_apply_xc(A_ia_rspace=work_matrices%A_ia_rspace_sub, kernel_env=kernel_env, &
rho_ia_struct=work_matrices%rho_orb_struct_sub, &
is_rks_triplets=is_rks_triplets, pw_env=sub_env%pw_env, &
work_v_xc=work_matrices%wpw_rspace_sub, &
work_v_xc_tau=work_matrices%wpw_tau_rspace_sub)
IF (gapw .OR. gapw_xc) THEN
rho_atom_set => sub_env%local_rho_set%rho_atom_set
rho1_atom_set => work_matrices%local_rho_set%rho_atom_set
CALL calculate_xc_2nd_deriv_atom(rho_atom_set, rho1_atom_set, qs_env, kernel_env%xc_section, &
sub_env%para_env, do_tddfpt2=.TRUE., do_triplet=is_rks_triplets, &
do_sf=do_noncol)
END IF
END IF
IF (gapw .OR. gapw_xc) THEN
rho_atom_set => sub_env%local_rho_set%rho_atom_set
rho1_atom_set => work_matrices%local_rho_set%rho_atom_set
CALL calculate_xc_2nd_deriv_atom(rho_atom_set, rho1_atom_set, qs_env, kernel_env%xc_section, &
sub_env%para_env, do_tddfpt2=.TRUE., do_triplet=is_rks_triplets)
END IF
END IF ! spin-flip
! ADMM correction
IF (do_admm .AND. admm_xc_correction) THEN
@ -314,7 +331,8 @@ CONTAINS
rho_ia_struct=work_matrices%rho_aux_fit_struct_sub, &
is_rks_triplets=is_rks_triplets, pw_env=sub_env%pw_env, &
work_v_xc=work_matrices%wpw_rspace_sub, &
work_v_xc_tau=work_matrices%wpw_tau_rspace_sub)
work_v_xc_tau=work_matrices%wpw_tau_rspace_sub, &
spinflip=spinflip)
DO ispin = 1, nspins
CALL pw_scale(V_rspace_sub(ispin), V_rspace_sub(ispin)%pw_grid%dvol)
CALL integrate_v_rspace(v_rspace=V_rspace_sub(ispin), &
@ -330,7 +348,7 @@ CONTAINS
CALL calculate_xc_2nd_deriv_atom(rho_atom_set, rho1_atom_set, qs_env, &
kernel_env_admm_aux%xc_section, &
sub_env%para_env, do_tddfpt2=.TRUE., &
do_triplet=is_rks_triplets, &
do_triplet=is_rks_triplets, do_sf=do_noncol, &
kind_set_external=admm_env%admm_gapw_env%admm_kind_set)
CALL update_ks_atom(qs_env, A_xc_munu_sub, rho_ia_ao_aux_fit, forces=.FALSE., tddft=.TRUE., &
rho_atom_external=rho1_atom_set, &
@ -373,7 +391,8 @@ CONTAINS
rho_ia_struct=work_matrices%rho_aux_fit_struct_sub, &
is_rks_triplets=is_rks_triplets, pw_env=sub_env%pw_env, &
work_v_xc=work_matrices%wpw_rspace_sub, &
work_v_xc_tau=work_matrices%wpw_tau_rspace_sub)
work_v_xc_tau=work_matrices%wpw_tau_rspace_sub, &
spinflip=spinflip)
IF (admm_env%do_gapw) THEN
CPWARN("GAPW/ADMM needs symmetric ADMM kernel")
CPABORT("GAPW/ADMM@TDDFT")
@ -383,7 +402,7 @@ CONTAINS
END IF
! electron-hole Coulomb interaction
IF (.NOT. is_rks_triplets) THEN
IF ((.NOT. is_rks_triplets) .AND. (spinflip .EQ. no_sf_tddfpt)) THEN
! a sum J_i{alpha}a{alpha}_munu + J_i{beta}a{beta}_munu can be computed by solving
! the Poisson equation for combined density (rho_{ia,alpha} + rho_{ia,beta}) .
! The following action will destroy reciprocal-space grid in spin-unrestricted case.

367
src/qs_tddfpt2_fhxc_forces.F
View file

@ -62,7 +62,9 @@ MODULE qs_tddfpt2_fhxc_forces
hfx_ri_update_ks
USE hfx_types, ONLY: hfx_type
USE input_constants, ONLY: do_admm_aux_exch_func_none,&
no_sf_tddfpt,&
tddfpt_kernel_full,&
tddfpt_sf_col,&
xc_kernel_method_analytic,&
xc_kernel_method_best,&
xc_kernel_method_numeric,&
@ -94,7 +96,8 @@ MODULE qs_tddfpt2_fhxc_forces
qs_environment_type,&
set_qs_env
USE qs_force_types, ONLY: qs_force_type
USE qs_fxc, ONLY: qs_fgxc_create,&
USE qs_fxc, ONLY: qs_fgxc_analytic,&
qs_fgxc_create,&
qs_fgxc_gdiff,&
qs_fgxc_release
USE qs_gapw_densities, ONLY: prepare_gapw_den
@ -148,10 +151,12 @@ MODULE qs_tddfpt2_fhxc_forces
CONTAINS
! **************************************************************************************************
!> \brief Calculate direct tddft forces
!> \param qs_env ...
!> \param ex_env ...
!> \param gs_mos ...
!> \brief Calculate direct tddft forces. Calculate the three last terms of the response vector
!> in equation 49 and the first term of \Lambda_munu in equation 51 in
!> J. Chem. Theory Comput. 2022, 18, 7, 4186–4202 (https://doi.org/10.1021/acs.jctc.2c00144)
!> \param qs_env Holds all system information relevant for the calculation.
!> \param ex_env Holds the response vector ex_env%cpmos and Lambda ex_env%matrix_wx1.
!> \param gs_mos MO coefficients of the ground state.
!> \param full_kernel ...
!> \param debug_forces ...
!> \par History
@ -171,10 +176,10 @@ CONTAINS
CHARACTER(LEN=default_string_length) :: basis_type
INTEGER :: handle, iounit, ispin, mspin, myfun, &
n_rep_hf, nao, nao_aux, natom, nkind, &
norb, nspins, order
norb(2), nspins, order, spin
LOGICAL :: distribute_fock_matrix, do_admm, do_analytic, do_hfx, do_hfxlr, do_hfxsr, &
do_numeric, gapw, gapw_xc, hfx_treat_lsd_in_core, is_rks_triplets, s_mstruct_changed, &
use_virial
do_numeric, do_sf, gapw, gapw_xc, hfx_treat_lsd_in_core, is_rks_triplets, &
s_mstruct_changed, use_virial
REAL(KIND=dp) :: eh1, eh1c, eps_delta, eps_fit, focc, &
fscal, fval, kval, xehartree
REAL(KIND=dp), DIMENSION(3) :: fodeb
@ -183,7 +188,7 @@ CONTAINS
TYPE(cp_fm_struct_type), POINTER :: fm_struct, fm_struct_mat
TYPE(cp_fm_type) :: cvcmat, vcvec
TYPE(cp_fm_type), DIMENSION(:), POINTER :: cpmos, evect
TYPE(cp_fm_type), POINTER :: mos
TYPE(cp_fm_type), POINTER :: mos, mos2
TYPE(cp_logger_type), POINTER :: logger
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_fx, matrix_gx, matrix_hfx, &
matrix_hfx_admm, matrix_hfx_admm_asymm, matrix_hfx_asymm, matrix_hx, matrix_p, &
@ -232,6 +237,11 @@ CONTAINS
tddfpt_control => dft_control%tddfpt2_control
nspins = dft_control%nspins
is_rks_triplets = tddfpt_control%rks_triplets .AND. (nspins == 1)
IF (tddfpt_control%spinflip == no_sf_tddfpt) THEN
do_sf = .FALSE.
ELSE
do_sf = .TRUE.
END IF
CPASSERT(tddfpt_control%kernel == tddfpt_kernel_full)
do_hfx = tddfpt_control%do_hfx
do_hfxsr = tddfpt_control%do_hfxsr
@ -248,20 +258,22 @@ CONTAINS
focc = 1.0_dp
IF (nspins == 2) focc = 0.5_dp
DO ispin = 1, nspins
DO ispin = 1, SIZE(evect, 1)
! Calculate (C*X^T + X*C^T)/2
CALL dbcsr_set(matrix_px1(ispin)%matrix, 0.0_dp)
CALL cp_fm_get_info(evect(ispin), ncol_global=norb)
CALL cp_fm_get_info(evect(ispin), ncol_global=norb(ispin))
CALL cp_dbcsr_plus_fm_fm_t(matrix_px1(ispin)%matrix, &
matrix_v=evect(ispin), &
matrix_g=gs_mos(ispin)%mos_occ, &
ncol=norb, alpha=2.0_dp*focc, symmetry_mode=1)
ncol=norb(ispin), alpha=2.0_dp*focc, symmetry_mode=1)
! Calculate (C*X^T - X*C^T)/2
CALL dbcsr_set(matrix_px1_asymm(ispin)%matrix, 0.0_dp)
CALL cp_fm_get_info(evect(ispin), ncol_global=norb)
CALL cp_fm_get_info(evect(ispin), ncol_global=norb(ispin))
CALL cp_dbcsr_plus_fm_fm_t(matrix_px1_asymm(ispin)%matrix, &
matrix_v=gs_mos(ispin)%mos_occ, &
matrix_g=evect(ispin), &
ncol=norb, alpha=2.0_dp*focc, &
ncol=norb(ispin), alpha=2.0_dp*focc, &
symmetry_mode=-1)
END DO
!
@ -270,7 +282,7 @@ 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
DO ispin = 1, SIZE(evect, 1)
CALL dbcsr_scale(matrix_px1(ispin)%matrix, 2.0_dp)
END DO
END IF
@ -298,7 +310,7 @@ CONTAINS
sap_oce, eps_fit)
CALL set_qs_env(qs_env, oce=oce)
mpga(1:nspins, 1:1) => matrix_px1(1:nspins)
mpga(1:SIZE(evect, 1), 1:1) => matrix_px1(1:SIZE(evect, 1))
CALL calculate_rho_atom_coeff(qs_env, mpga, local_rho_set%rho_atom_set, &
qs_kind_set, oce, sab, para_env)
CALL prepare_gapw_den(qs_env, local_rho_set, do_rho0=gapw)
@ -317,7 +329,7 @@ CONTAINS
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
DO ispin = 1, SIZE(evect, 1)
CALL dbcsr_scale(matrix_px1(ispin)%matrix, 0.5_dp)
END DO
END IF
@ -327,8 +339,9 @@ CONTAINS
CALL get_qs_env(qs_env, admm_env=admm_env)
nao_aux = admm_env%nao_aux_fit
nao = admm_env%nao_orb
!
DO ispin = 1, nspins
! Fit the symmetrized and antisymmetrized matrices
DO ispin = 1, SIZE(evect, 1)
CALL copy_dbcsr_to_fm(matrix_px1(ispin)%matrix, admm_env%work_orb_orb)
CALL parallel_gemm('N', 'N', nao_aux, nao, nao, &
1.0_dp, admm_env%A, admm_env%work_orb_orb, 0.0_dp, &
@ -359,7 +372,7 @@ CONTAINS
CALL local_rho_set_create(local_rho_set_admm)
CALL allocate_rho_atom_internals(local_rho_set_admm%rho_atom_set, atomic_kind_set, &
admm_env%admm_gapw_env%admm_kind_set, dft_control, para_env)
mpga(1:nspins, 1:1) => matrix_px1_admm(1:nspins)
mpga(1:SIZE(evect, 1), 1:1) => matrix_px1_admm(1:SIZE(evect, 1))
CALL get_admm_env(admm_env, sab_aux_fit=sab_aux_fit)
CALL calculate_rho_atom_coeff(qs_env, mpga, local_rho_set_admm%rho_atom_set, &
admm_env%admm_gapw_env%admm_kind_set, &
@ -392,30 +405,33 @@ CONTAINS
CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
poisson_env=poisson_env)
ALLOCATE (rhox_r(nspins), rhox_g(nspins))
DO ispin = 1, nspins
ALLOCATE (rhox_r(SIZE(evect, 1)), rhox_g(SIZE(evect, 1)))
DO ispin = 1, SIZE(evect, 1)
CALL auxbas_pw_pool%create_pw(rhox_r(ispin))
CALL auxbas_pw_pool%create_pw(rhox_g(ispin))
END DO
CALL auxbas_pw_pool%create_pw(rhox_tot_gspace)
CALL pw_zero(rhox_tot_gspace)
DO ispin = 1, nspins
DO ispin = 1, SIZE(evect, 1)
! Calculate gridpoint values of the density associated to 2*matrix_px1 = C*X^T + X*C^T
IF (nspins == 2) CALL dbcsr_scale(matrix_px1(ispin)%matrix, 2.0_dp)
CALL calculate_rho_elec(ks_env=ks_env, matrix_p=matrix_px1(ispin)%matrix, &
rho=rhox_r(ispin), rho_gspace=rhox_g(ispin), &
soft_valid=gapw)
! rhox_tot_gspace contains the values on the grid points of rhox = sum_munu 4D^X_munu*mu(r)*nu(r)
CALL pw_axpy(rhox_g(ispin), rhox_tot_gspace)
! Recover matrix_px1 = (C*X^T + X*C^T)/2
IF (nspins == 2) CALL dbcsr_scale(matrix_px1(ispin)%matrix, 0.5_dp)
END DO
IF (gapw_xc) THEN
ALLOCATE (rhoxx_r(nspins), rhoxx_g(nspins))
DO ispin = 1, nspins
ALLOCATE (rhoxx_r(SIZE(evect, 1)), rhoxx_g(SIZE(evect, 1)))
DO ispin = 1, SIZE(evect, 1)
CALL auxbas_pw_pool%create_pw(rhoxx_r(ispin))
CALL auxbas_pw_pool%create_pw(rhoxx_g(ispin))
END DO
DO ispin = 1, nspins
DO ispin = 1, SIZE(evect, 1)
IF (nspins == 2) CALL dbcsr_scale(matrix_px1(ispin)%matrix, 2.0_dp)
CALL calculate_rho_elec(ks_env=ks_env, matrix_p=matrix_px1(ispin)%matrix, &
rho=rhoxx_r(ispin), rho_gspace=rhoxx_g(ispin), &
@ -426,7 +442,7 @@ CONTAINS
CALL get_qs_env(qs_env, matrix_s=matrix_s, force=force)
IF (.NOT. is_rks_triplets) THEN
IF (.NOT. (is_rks_triplets .OR. do_sf)) THEN
CALL auxbas_pw_pool%create_pw(xv_hartree_rspace)
CALL auxbas_pw_pool%create_pw(xv_hartree_gspace)
! calculate associated hartree potential
@ -456,7 +472,7 @@ CONTAINS
IF (debug_forces) THEN
fodeb(1:3) = force(1)%rho_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] ", fodeb
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Dx*dKh[Dx] ", fodeb
END IF
IF (gapw) THEN
IF (debug_forces) fodeb(1:3) = force(1)%g0s_Vh_elec(1:3, 1)
@ -472,7 +488,7 @@ CONTAINS
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
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Dx*dKh[Dx]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., &
@ -480,7 +496,7 @@ CONTAINS
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*dKh[X]PAW ", fodeb
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Dx*dKh[Dx]PAW", fodeb
END IF
END IF
END IF
@ -493,7 +509,7 @@ CONTAINS
xc_section => full_kernel%xc_section
CALL section_vals_val_get(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", &
i_val=myfun)
IF (myfun /= xc_none) THEN
IF (.NOT. ((myfun == xc_none) .OR. (tddfpt_control%spinflip == tddfpt_sf_col))) THEN
SELECT CASE (ex_env%xc_kernel_method)
CASE (xc_kernel_method_best)
do_analytic = .TRUE.
@ -518,6 +534,8 @@ CONTAINS
CALL qs_rho_get(rho, rho_ao=matrix_p)
NULLIFY (rhox)
ALLOCATE (rhox)
! Create rhox object to collect all information on matrix_px1, including its values on the
! grid points
CALL qs_rho_create(rhox)
IF (gapw_xc) THEN
CALL qs_rho_set(rho_struct=rhox, rho_ao=matrix_px1, rho_r=rhoxx_r, rho_g=rhoxx_g, &
@ -526,12 +544,19 @@ CONTAINS
CALL qs_rho_set(rho_struct=rhox, rho_ao=matrix_px1, rho_r=rhox_r, rho_g=rhox_g, &
rho_r_valid=.TRUE., rho_g_valid=.TRUE.)
END IF
! Calculate the exchange-correlation kernel derivative contribution, notice that for open-shell
! rhox_r contains a factor of 2!
IF (do_analytic .AND. .NOT. do_numeric) THEN
CPABORT("Analytic 3rd EXC derivatives not available")
IF (.NOT. do_sf) THEN
CPABORT("Analytic 3rd EXC derivatives not available")
ELSE !TODO
CALL qs_fgxc_analytic(rho, rhox, xc_section, auxbas_pw_pool, &
fxc_rho, fxc_tau, gxc_rho, gxc_tau, spinflip=do_sf)
END IF
ELSEIF (do_numeric) THEN
IF (do_analytic) THEN
CALL qs_fgxc_gdiff(ks_env, rho, rhox, xc_section, order, eps_delta, is_rks_triplets, &
fxc_rho, fxc_tau, gxc_rho, gxc_tau)
fxc_rho, fxc_tau, gxc_rho, gxc_tau, spinflip=do_sf)
ELSE
CALL qs_fgxc_create(ks_env, rho, rhox, xc_section, order, is_rks_triplets, &
fxc_rho, fxc_tau, gxc_rho, gxc_tau)
@ -572,13 +597,15 @@ CONTAINS
IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
NULLIFY (matrix_fx)
CALL dbcsr_allocate_matrix_set(matrix_fx, nspins)
DO ispin = 1, nspins
CALL dbcsr_allocate_matrix_set(matrix_fx, SIZE(fxc_rho))
DO ispin = 1, SIZE(fxc_rho, 1)
ALLOCATE (matrix_fx(ispin)%matrix)
CALL dbcsr_create(matrix_fx(ispin)%matrix, template=matrix_s(1)%matrix)
CALL dbcsr_copy(matrix_fx(ispin)%matrix, matrix_s(1)%matrix)
CALL dbcsr_set(matrix_fx(ispin)%matrix, 0.0_dp)
CALL pw_scale(fxc_rho(ispin), fxc_rho(ispin)%pw_grid%dvol)
! Calculate 2sum_sigmatau<munu|fxc|sigmatau>D^X_sigmatau
! fxc_rho here containes a factor of 2
CALL integrate_v_rspace(qs_env=qs_env, v_rspace=fxc_rho(ispin), &
pmat=matrix_px1(ispin), hmat=matrix_fx(ispin), &
gapw=(gapw .OR. gapw_xc), calculate_forces=.TRUE.)
@ -586,12 +613,14 @@ CONTAINS
IF (debug_forces) THEN
fodeb(1:3) = force(1)%rho_elec(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] ", fodeb
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Dx*dfxc[Dx] ", fodeb
END IF
IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
NULLIFY (matrix_gx)
CALL dbcsr_allocate_matrix_set(matrix_gx, nspins)
! Calculate exchange-correlation kernel derivative 2<D^X D^X|gxc|mu nu>
! gxc comes with a factor of 4, so a factor of 1/2 is introduced
DO ispin = 1, nspins
ALLOCATE (matrix_gx(ispin)%matrix)
CALL dbcsr_create(matrix_gx(ispin)%matrix, template=matrix_s(1)%matrix)
@ -607,7 +636,7 @@ CONTAINS
IF (debug_forces) THEN
fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
CALL para_env%sum(fodeb)
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dKg[X] ", fodeb
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Dx*dgxc[Dx]", fodeb
END IF
CALL qs_fgxc_release(ks_env, fxc_rho, fxc_tau, gxc_rho, gxc_tau)
@ -619,7 +648,7 @@ CONTAINS
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
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Dx*dfxc[Dx]PAW ", fodeb
END IF
IF (debug_forces) fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1)
IF (nspins == 1) THEN
@ -634,13 +663,13 @@ CONTAINS
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*dKg[X]PAW ", fodeb
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Dx*dgxc[Dx]PAW ", fodeb
END IF
END IF
END IF
! ADMM XC correction Exc[rho_admm]
IF (do_admm .AND. tddfpt_control%admm_xc_correction) THEN
IF (do_admm .AND. tddfpt_control%admm_xc_correction .AND. (tddfpt_control%spinflip .NE. tddfpt_sf_col)) THEN
IF (admm_env%aux_exch_func == do_admm_aux_exch_func_none) THEN
! nothing to do
ELSE
@ -658,13 +687,16 @@ CONTAINS
END IF
!
NULLIFY (mfx, mgx)
CALL dbcsr_allocate_matrix_set(mfx, nspins)
CALL dbcsr_allocate_matrix_set(mfx, SIZE(evect, 1))
CALL dbcsr_allocate_matrix_set(mgx, nspins)
DO ispin = 1, nspins
ALLOCATE (mfx(ispin)%matrix, mgx(ispin)%matrix)
DO ispin = 1, SIZE(evect, 1)
ALLOCATE (mfx(ispin)%matrix)
CALL dbcsr_create(mfx(ispin)%matrix, template=matrix_s_aux_fit(1)%matrix)
CALL dbcsr_copy(mfx(ispin)%matrix, matrix_s_aux_fit(1)%matrix)
CALL dbcsr_set(mfx(ispin)%matrix, 0.0_dp)
END DO
DO ispin = 1, nspins
ALLOCATE (mgx(ispin)%matrix)
CALL dbcsr_create(mgx(ispin)%matrix, template=matrix_s_aux_fit(1)%matrix)
CALL dbcsr_copy(mgx(ispin)%matrix, matrix_s_aux_fit(1)%matrix)
CALL dbcsr_set(mgx(ispin)%matrix, 0.0_dp)
@ -675,12 +707,12 @@ CONTAINS
CALL qs_rho_get(rho_aux_fit, rho_r=rho_r_aux, rho_g=rho_g_aux)
CALL qs_rho_get(rho_aux_fit, rho_ao=matrix_p_admm)
! rhox_aux
ALLOCATE (rhox_r_aux(nspins), rhox_g_aux(nspins))
DO ispin = 1, nspins
ALLOCATE (rhox_r_aux(SIZE(evect, 1)), rhox_g_aux(SIZE(evect, 1)))
DO ispin = 1, SIZE(evect, 1)
CALL auxbas_pw_pool%create_pw(rhox_r_aux(ispin))
CALL auxbas_pw_pool%create_pw(rhox_g_aux(ispin))
END DO
DO ispin = 1, nspins
DO ispin = 1, SIZE(evect, 1)
CALL calculate_rho_elec(ks_env=ks_env, matrix_p=matrix_px1_admm(ispin)%matrix, &
rho=rhox_r_aux(ispin), rho_gspace=rhox_g_aux(ispin), &
basis_type=basis_type, &
@ -713,7 +745,7 @@ CONTAINS
! because this would release the arrays. Instead we're simply going to deallocate rhox_aux.
DEALLOCATE (rhox_aux)
DO ispin = 1, nspins
DO ispin = 1, SIZE(evect, 1)
CALL auxbas_pw_pool%give_back_pw(rhox_r_aux(ispin))
CALL auxbas_pw_pool%give_back_pw(rhox_g_aux(ispin))
END DO
@ -722,7 +754,7 @@ CONTAINS
IF (nspins == 2) fscal = 2.0_dp
!
IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
DO ispin = 1, nspins
DO ispin = 1, SIZE(evect, 1)
CALL pw_scale(fxc_rho(ispin), fxc_rho(ispin)%pw_grid%dvol)
CALL integrate_v_rspace(qs_env=qs_env, v_rspace=fxc_rho(ispin), &
hmat=mfx(ispin), &
@ -735,11 +767,11 @@ CONTAINS
IF (debug_forces) THEN
fodeb(1:3) = force(1)%rho_elec(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]ADMM", fodeb
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Dx*dfxc[Dx]ADMM", fodeb
END IF
IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
DO ispin = 1, nspins
DO ispin = 1, SIZE(evect, 1)
CALL pw_scale(gxc_rho(ispin), gxc_rho(ispin)%pw_grid%dvol)
CALL pw_scale(gxc_rho(ispin), 0.5_dp)
CALL integrate_v_rspace(qs_env=qs_env, v_rspace=gxc_rho(ispin), &
@ -754,7 +786,7 @@ CONTAINS
IF (debug_forces) THEN
fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
CALL para_env%sum(fodeb)
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dKg[X]ADMM", fodeb
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Dx*dgxc[Dx]ADMM", fodeb
END IF
CALL qs_fgxc_release(ks_env, fxc_rho, fxc_tau, gxc_rho, gxc_tau)
!
@ -799,7 +831,7 @@ CONTAINS
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]ADMM-PAW ", fodeb
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Dx*dfxc[Dx]ADMM-PAW ", fodeb
END IF
IF (debug_forces) fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1)
IF (nspins == 1) THEN
@ -818,7 +850,7 @@ CONTAINS
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:: P*dKg[X]ADMM-PAW ", fodeb
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Dx*dgxc[Dx]ADMM-PAW ", fodeb
END IF
END IF
!
@ -848,7 +880,7 @@ CONTAINS
nao_aux = admm_env%nao_aux_fit
ALLOCATE (dbwork)
CALL dbcsr_create(dbwork, template=matrix_fx(1)%matrix)
DO ispin = 1, nspins
DO ispin = 1, SIZE(evect, 1)
! fx
CALL cp_dbcsr_sm_fm_multiply(mfx(ispin)%matrix, admm_env%A, &
admm_env%work_aux_orb, nao)
@ -877,20 +909,20 @@ CONTAINS
END IF
END IF
DO ispin = 1, nspins
DO ispin = 1, SIZE(evect, 1)
CALL auxbas_pw_pool%give_back_pw(rhox_r(ispin))
CALL auxbas_pw_pool%give_back_pw(rhox_g(ispin))
END DO
DEALLOCATE (rhox_r, rhox_g)
CALL auxbas_pw_pool%give_back_pw(rhox_tot_gspace)
IF (gapw_xc) THEN
DO ispin = 1, nspins
DO ispin = 1, SIZE(evect, 1)
CALL auxbas_pw_pool%give_back_pw(rhoxx_r(ispin))
CALL auxbas_pw_pool%give_back_pw(rhoxx_g(ispin))
END DO
DEALLOCATE (rhoxx_r, rhoxx_g)
END IF
IF (.NOT. is_rks_triplets) THEN
IF (.NOT. (is_rks_triplets .OR. do_sf)) THEN
CALL auxbas_pw_pool%give_back_pw(xv_hartree_rspace)
CALL auxbas_pw_pool%give_back_pw(xv_hartree_gspace)
END IF
@ -898,9 +930,9 @@ CONTAINS
! HFX
IF (do_hfx) THEN
NULLIFY (matrix_hfx, matrix_hfx_asymm)
CALL dbcsr_allocate_matrix_set(matrix_hfx, nspins)
CALL dbcsr_allocate_matrix_set(matrix_hfx_asymm, nspins)
DO ispin = 1, nspins
CALL dbcsr_allocate_matrix_set(matrix_hfx, SIZE(evect, 1))
CALL dbcsr_allocate_matrix_set(matrix_hfx_asymm, SIZE(evect, 1))
DO ispin = 1, SIZE(evect, 1)
ALLOCATE (matrix_hfx(ispin)%matrix)
CALL dbcsr_create(matrix_hfx(ispin)%matrix, template=matrix_s(1)%matrix)
CALL dbcsr_copy(matrix_hfx(ispin)%matrix, matrix_s(1)%matrix)
@ -919,7 +951,7 @@ CONTAINS
CALL section_vals_val_get(hfx_section, "TREAT_LSD_IN_CORE", l_val=hfx_treat_lsd_in_core, &
i_rep_section=1)
mspin = 1
IF (hfx_treat_lsd_in_core) mspin = nspins
IF (hfx_treat_lsd_in_core) mspin = SIZE(evect, 1)
!
CALL get_qs_env(qs_env=qs_env, x_data=x_data, s_mstruct_changed=s_mstruct_changed)
distribute_fock_matrix = .TRUE.
@ -927,9 +959,9 @@ CONTAINS
IF (do_admm) THEN
CALL get_admm_env(qs_env%admm_env, matrix_s_aux_fit=matrix_s_aux_fit)
NULLIFY (matrix_hfx_admm, matrix_hfx_admm_asymm)
CALL dbcsr_allocate_matrix_set(matrix_hfx_admm, nspins)
CALL dbcsr_allocate_matrix_set(matrix_hfx_admm_asymm, nspins)
DO ispin = 1, nspins
CALL dbcsr_allocate_matrix_set(matrix_hfx_admm, SIZE(evect, 1))
CALL dbcsr_allocate_matrix_set(matrix_hfx_admm_asymm, SIZE(evect, 1))
DO ispin = 1, SIZE(evect, 1)
ALLOCATE (matrix_hfx_admm(ispin)%matrix)
CALL dbcsr_create(matrix_hfx_admm(ispin)%matrix, template=matrix_s_aux_fit(1)%matrix)
CALL dbcsr_copy(matrix_hfx_admm(ispin)%matrix, matrix_s_aux_fit(1)%matrix)
@ -942,8 +974,8 @@ CONTAINS
END DO
!
NULLIFY (mpe, mhe)
ALLOCATE (mpe(nspins, 1), mhe(nspins, 1))
DO ispin = 1, nspins
ALLOCATE (mpe(SIZE(evect, 1), 1), mhe(SIZE(evect, 1), 1))
DO ispin = 1, SIZE(evect)
mhe(ispin, 1)%matrix => matrix_hfx_admm(ispin)%matrix
mpe(ispin, 1)%matrix => matrix_px1_admm(ispin)%matrix
END DO
@ -957,11 +989,11 @@ CONTAINS
eh1 = 0.0
CALL integrate_four_center(qs_env, x_data, mhe, eh1, mpe, hfx_section, &
para_env, s_mstruct_changed, 1, distribute_fock_matrix, &
ispin=ispin)
ispin=ispin, nspins=SIZE(evect, 1))
END DO
END IF
!anti-symmetric density
DO ispin = 1, nspins
DO ispin = 1, SIZE(evect, 1)
mhe(ispin, 1)%matrix => matrix_hfx_admm_asymm(ispin)%matrix
mpe(ispin, 1)%matrix => matrix_px1_admm_asymm(ispin)%matrix
END DO
@ -975,7 +1007,7 @@ CONTAINS
eh1 = 0.0
CALL integrate_four_center(qs_env, x_data, mhe, eh1, mpe, hfx_section, &
para_env, s_mstruct_changed, 1, distribute_fock_matrix, &
ispin=ispin)
ispin=ispin, nspins=SIZE(evect, 1))
END DO
END IF
!
@ -984,7 +1016,7 @@ CONTAINS
ALLOCATE (dbwork, dbwork_asymm)
CALL dbcsr_create(dbwork, template=matrix_hfx(1)%matrix)
CALL dbcsr_create(dbwork_asymm, template=matrix_hfx(1)%matrix, matrix_type=dbcsr_type_antisymmetric)
DO ispin = 1, nspins
DO ispin = 1, SIZE(evect, 1)
CALL cp_dbcsr_sm_fm_multiply(matrix_hfx_admm(ispin)%matrix, admm_env%A, &
admm_env%work_aux_orb, nao)
CALL parallel_gemm('T', 'N', nao, nao, nao_aux, &
@ -1023,8 +1055,10 @@ CONTAINS
use_virial = .FALSE.
NULLIFY (mdum)
fval = 2.0_dp*REAL(nspins, KIND=dp)*0.5_dp !0.5 factor because of symemtry/anti-symmetry
! For SF TDDFT integrate_four_center and derivatives_four_center routines introduce a factor of 1/2
IF (do_sf) fval = fval*2.0_dp
IF (debug_forces) fodeb(1:3) = force(1)%fock_4c(1:3, 1)
DO ispin = 1, nspins
DO ispin = 1, SIZE(evect, 1)
mpe(ispin, 1)%matrix => matrix_px1_admm(ispin)%matrix
END DO
IF (x_data(1, 1)%do_hfx_ri) THEN
@ -1034,9 +1068,9 @@ CONTAINS
use_virial=use_virial, rescale_factor=fval)
ELSE
CALL derivatives_four_center(qs_env, mpe, mdum, hfx_section, para_env, 1, use_virial, &
adiabatic_rescale_factor=fval)
adiabatic_rescale_factor=fval, nspins=SIZE(mpe, 1))
END IF
DO ispin = 1, nspins
DO ispin = 1, SIZE(evect, 1)
mpe(ispin, 1)%matrix => matrix_px1_admm_asymm(ispin)%matrix
END DO
IF (x_data(1, 1)%do_hfx_ri) THEN
@ -1046,12 +1080,12 @@ CONTAINS
use_virial=use_virial, rescale_factor=fval)
ELSE
CALL derivatives_four_center(qs_env, mpe, mdum, hfx_section, para_env, 1, use_virial, &
adiabatic_rescale_factor=fval)
adiabatic_rescale_factor=fval, nspins=SIZE(mpe, 1))
END IF
IF (debug_forces) THEN
fodeb(1:3) = force(1)%fock_4c(1:3, 1) - fodeb(1:3)
CALL para_env%sum(fodeb)
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Px*hfx'*Px ", fodeb
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Dx*dhfx'*Dx ", fodeb
END IF
!
DEALLOCATE (mpe, mhe)
@ -1060,8 +1094,8 @@ CONTAINS
CALL dbcsr_deallocate_matrix_set(matrix_hfx_admm_asymm)
ELSE
NULLIFY (mpe, mhe)
ALLOCATE (mpe(nspins, 1), mhe(nspins, 1))
DO ispin = 1, nspins
ALLOCATE (mpe(SIZE(evect, 1), 1), mhe(SIZE(evect, 1), 1))
DO ispin = 1, SIZE(evect, 1)
mhe(ispin, 1)%matrix => matrix_hfx(ispin)%matrix
mpe(ispin, 1)%matrix => matrix_px1(ispin)%matrix
END DO
@ -1075,12 +1109,12 @@ CONTAINS
eh1 = 0.0
CALL integrate_four_center(qs_env, x_data, mhe, eh1, mpe, hfx_section, &
para_env, s_mstruct_changed, 1, distribute_fock_matrix, &
ispin=ispin)
ispin=ispin, nspins=SIZE(mpe, 1))
END DO
END IF
!anti-symmetric density matrix
DO ispin = 1, nspins
DO ispin = 1, SIZE(evect, 1)
mhe(ispin, 1)%matrix => matrix_hfx_asymm(ispin)%matrix
mpe(ispin, 1)%matrix => matrix_px1_asymm(ispin)%matrix
END DO
@ -1094,15 +1128,17 @@ CONTAINS
eh1 = 0.0
CALL integrate_four_center(qs_env, x_data, mhe, eh1, mpe, hfx_section, &
para_env, s_mstruct_changed, 1, distribute_fock_matrix, &
ispin=ispin)
ispin=ispin, nspins=SIZE(mpe, 1))
END DO
END IF
! forces
use_virial = .FALSE.
NULLIFY (mdum)
fval = 2.0_dp*REAL(nspins, KIND=dp)*0.5_dp !extra 0.5 factor because of symmetry/antisymemtry
! For SF TDDFT integrate_four_center and derivatives_four_center routines introduce a factor of 1/2
IF (do_sf) fval = fval*2.0_dp
IF (debug_forces) fodeb(1:3) = force(1)%fock_4c(1:3, 1)
DO ispin = 1, nspins
DO ispin = 1, SIZE(evect, 1)
mpe(ispin, 1)%matrix => matrix_px1(ispin)%matrix
END DO
IF (x_data(1, 1)%do_hfx_ri) THEN
@ -1112,9 +1148,9 @@ CONTAINS
use_virial=use_virial, rescale_factor=fval)
ELSE
CALL derivatives_four_center(qs_env, mpe, mdum, hfx_section, para_env, 1, use_virial, &
adiabatic_rescale_factor=fval)
adiabatic_rescale_factor=fval, nspins=SIZE(mpe, 1))
END IF
DO ispin = 1, nspins
DO ispin = 1, SIZE(evect, 1)
mpe(ispin, 1)%matrix => matrix_px1_asymm(ispin)%matrix
END DO
IF (x_data(1, 1)%do_hfx_ri) THEN
@ -1124,18 +1160,20 @@ CONTAINS
use_virial=use_virial, rescale_factor=fval)
ELSE
CALL derivatives_four_center(qs_env, mpe, mdum, hfx_section, para_env, 1, use_virial, &
adiabatic_rescale_factor=fval)
adiabatic_rescale_factor=fval, nspins=SIZE(mpe, 1))
END IF
IF (debug_forces) THEN
fodeb(1:3) = force(1)%fock_4c(1:3, 1) - fodeb(1:3)
CALL para_env%sum(fodeb)
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Px*hfx'*Px ", fodeb
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Dx*dhfx*Dx ", fodeb
END IF
!
DEALLOCATE (mpe, mhe)
END IF
fval = 2.0_dp*REAL(nspins, KIND=dp)*0.5_dp !extra 0.5 because of symm/antisymm
DO ispin = 1, nspins
! For SF TDDFT integrate_four_center and derivatives_four_center routines introduce a factor of 1/2
IF (do_sf) fval = fval*2.0_dp
DO ispin = 1, SIZE(evect, 1)
CALL dbcsr_scale(matrix_hfx(ispin)%matrix, fval)
CALL dbcsr_scale(matrix_hfx_asymm(ispin)%matrix, fval)
END DO
@ -1176,83 +1214,140 @@ CONTAINS
cpmos => ex_env%cpmos
focc = 2.0_dp
IF (nspins == 2) focc = 1.0_dp
! Initialize mos and dimensions of occupied space
! In the following comments mos is referred to as Ca and mos2 as Cb
spin = 1
mos => gs_mos(1)%mos_occ
CALL cp_fm_get_info(mos, ncol_global=norb(1))
IF (nspins == 2) THEN
mos2 => gs_mos(2)%mos_occ
CALL cp_fm_get_info(mos2, ncol_global=norb(2))
END IF
! Build response vector, Eq. 49, and the third term of \Lamda_munu, Eq. 51
DO ispin = 1, nspins
mos => gs_mos(ispin)%mos_occ
CALL cp_fm_get_info(evect(ispin), ncol_global=norb)
CALL cp_fm_create(vcvec, mos%matrix_struct, "vcvec")
! Initialize mos and dimensions of occupied space
IF (.NOT. do_sf) THEN
spin = ispin
mos => gs_mos(ispin)%mos_occ
mos2 => gs_mos(ispin)%mos_occ
END IF
! Create working fields for the response vector
CALL cp_fm_create(vcvec, gs_mos(ispin)%mos_occ%matrix_struct, "vcvec")
CALL cp_fm_get_info(vcvec, matrix_struct=fm_struct, nrow_global=nao)
CALL cp_fm_struct_create(fm_struct_mat, context=fm_struct%context, nrow_global=norb, &
ncol_global=norb, para_env=fm_struct%para_env)
CALL cp_fm_struct_create(fm_struct_mat, context=fm_struct%context, nrow_global=norb(spin), &
ncol_global=norb(ispin), para_env=fm_struct%para_env)
CALL cp_fm_create(cvcmat, fm_struct_mat)
CALL cp_fm_struct_release(fm_struct_mat)
!
! Allocate and initialize the Lambda matrix
ALLOCATE (matrix_wx1(ispin)%matrix)
CALL dbcsr_create(matrix=matrix_wx1(ispin)%matrix, template=matrix_s(1)%matrix)
CALL cp_dbcsr_alloc_block_from_nbl(matrix_wx1(ispin)%matrix, sab_orb)
CALL dbcsr_set(matrix_wx1(ispin)%matrix, 0.0_dp)
!
IF (.NOT. is_rks_triplets) THEN
! Add Hartree contributions to the perturbation vector
IF (.NOT. (is_rks_triplets .OR. do_sf)) THEN
CALL cp_dbcsr_sm_fm_multiply(matrix_hx(ispin)%matrix, evect(ispin), &
cpmos(ispin), norb, alpha=focc, beta=1.0_dp)
CALL cp_dbcsr_sm_fm_multiply(matrix_hx(ispin)%matrix, mos, vcvec, norb, alpha=1.0_dp, beta=0.0_dp)
CALL parallel_gemm("T", "N", norb, norb, nao, 1.0_dp, mos, vcvec, 0.0_dp, cvcmat)
CALL parallel_gemm("N", "N", nao, norb, norb, 1.0_dp, evect(ispin), cvcmat, 0.0_dp, vcvec)
cpmos(ispin), norb(ispin), alpha=focc, beta=1.0_dp)
CALL cp_dbcsr_sm_fm_multiply(matrix_hx(ispin)%matrix, mos, vcvec, norb(ispin), alpha=1.0_dp, beta=0.0_dp)
CALL parallel_gemm("T", "N", norb(ispin), norb(ispin), nao, 1.0_dp, mos, vcvec, 0.0_dp, cvcmat)
CALL parallel_gemm("N", "N", nao, norb(ispin), norb(ispin), 1.0_dp, evect(ispin), cvcmat, 0.0_dp, vcvec)
CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, vcvec, cpmos(ispin), &
norb, alpha=-focc, beta=1.0_dp)
norb(ispin), alpha=-focc, beta=1.0_dp)
!
CALL cp_dbcsr_plus_fm_fm_t(matrix_wx1(ispin)%matrix, matrix_v=mos, matrix_g=vcvec, &
ncol=norb, alpha=2.0_dp, symmetry_mode=1)
ncol=norb(ispin), alpha=2.0_dp, symmetry_mode=1)
END IF
!
IF (myfun /= xc_none) THEN
CALL cp_dbcsr_sm_fm_multiply(matrix_fx(ispin)%matrix, evect(ispin), &
cpmos(ispin), norb, alpha=focc, beta=1.0_dp)
CALL cp_dbcsr_sm_fm_multiply(matrix_fx(ispin)%matrix, mos, vcvec, norb, alpha=1.0_dp, beta=0.0_dp)
CALL parallel_gemm("T", "N", norb, norb, nao, 1.0_dp, mos, vcvec, 0.0_dp, cvcmat)
CALL parallel_gemm("N", "N", nao, norb, norb, 1.0_dp, evect(ispin), cvcmat, 0.0_dp, vcvec)
CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, vcvec, cpmos(ispin), &
norb, alpha=-focc, beta=1.0_dp)
! Add exchange-correlation kernel and exchange-correlation kernel derivative contributions to the response vector
IF ((myfun /= xc_none) .AND. (tddfpt_control%spinflip /= tddfpt_sf_col)) THEN
! XC Kernel contributions
! For spin-flip excitations this is the only contribution to the alpha response vector
IF (.NOT. (do_sf .AND. (ispin .EQ. 2))) THEN
! F*X
CALL cp_dbcsr_sm_fm_multiply(matrix_fx(spin)%matrix, evect(spin), &
cpmos(ispin), norb(ispin), alpha=focc, beta=1.0_dp)
END IF
! For spin-flip excitations this is the only contribution to the beta response vector
IF (.NOT. (do_sf .AND. (ispin .EQ. 1))) THEN
! F*Cb
CALL cp_dbcsr_sm_fm_multiply(matrix_fx(spin)%matrix, mos2, vcvec, norb(ispin), alpha=1.0_dp, beta=0.0_dp)
! Ca^T*F*Cb
CALL parallel_gemm("T", "N", norb(spin), norb(ispin), nao, 1.0_dp, mos, vcvec, 0.0_dp, cvcmat)
! X*Ca^T*F*Cb
CALL parallel_gemm("N", "N", nao, norb(ispin), norb(spin), 1.0_dp, evect(spin), cvcmat, 0.0_dp, vcvec)
! -S*X*Ca^T*F*Cb
CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, vcvec, cpmos(ispin), &
norb(ispin), alpha=-focc, beta=1.0_dp)
! Add contributions to the \Lambda_munu for the perturbed overlap matrix term, third term of Eq. 51
! 2X*Ca^T*F*Cb*Cb^T
CALL cp_dbcsr_plus_fm_fm_t(matrix_wx1(ispin)%matrix, matrix_v=vcvec, matrix_g=gs_mos(ispin)%mos_occ, &
ncol=norb(ispin), alpha=2.0_dp, symmetry_mode=1)
END IF
!
CALL cp_dbcsr_plus_fm_fm_t(matrix_wx1(ispin)%matrix, matrix_v=mos, matrix_g=vcvec, &
ncol=norb, alpha=2.0_dp, symmetry_mode=1)
! XC g (third functional derivative) contributions
! g*Ca*focc
CALL cp_dbcsr_sm_fm_multiply(matrix_gx(ispin)%matrix, gs_mos(ispin)%mos_occ, &
cpmos(ispin), norb(ispin), alpha=focc, beta=1.0_dp)
! Add contributions to the \Lambda_munu for the perturbed overlap matrix term, third term of Eq. 51
! g*Ca
CALL cp_dbcsr_sm_fm_multiply(matrix_gx(ispin)%matrix, gs_mos(ispin)%mos_occ, vcvec, norb(ispin), &
alpha=1.0_dp, beta=0.0_dp)
! Ca^T*g*Ca
CALL parallel_gemm("T", "N", norb(ispin), norb(ispin), nao, 1.0_dp, gs_mos(ispin)%mos_occ, vcvec, 0.0_dp, cvcmat)
! Ca*Ca^T*g*Ca
CALL parallel_gemm("N", "N", nao, norb(ispin), norb(ispin), 1.0_dp, gs_mos(ispin)%mos_occ, cvcmat, 0.0_dp, vcvec)
! Ca*Ca^T*g*Ca*Ca^T
CALL cp_dbcsr_plus_fm_fm_t(matrix_wx1(ispin)%matrix, matrix_v=vcvec, matrix_g=gs_mos(ispin)%mos_occ, &
ncol=norb(ispin), alpha=1.0_dp, symmetry_mode=1)
!
CALL cp_dbcsr_sm_fm_multiply(matrix_gx(ispin)%matrix, mos, &
cpmos(ispin), norb, alpha=focc, beta=1.0_dp)
!
CALL cp_dbcsr_sm_fm_multiply(matrix_gx(ispin)%matrix, mos, vcvec, norb, alpha=1.0_dp, beta=0.0_dp)
CALL parallel_gemm("T", "N", norb, norb, nao, 1.0_dp, mos, vcvec, 0.0_dp, cvcmat)
CALL parallel_gemm("N", "N", nao, norb, norb, 1.0_dp, mos, cvcmat, 0.0_dp, vcvec)
CALL cp_dbcsr_plus_fm_fm_t(matrix_wx1(ispin)%matrix, matrix_v=mos, matrix_g=vcvec, &
ncol=norb, alpha=1.0_dp, symmetry_mode=1)
END IF
!
! Add Fock contributions to the response vector
IF (do_hfx) THEN
CALL cp_dbcsr_sm_fm_multiply(matrix_hfx(ispin)%matrix, evect(ispin), &
cpmos(ispin), norb, alpha=focc, beta=1.0_dp)
CALL cp_dbcsr_sm_fm_multiply(matrix_hfx(ispin)%matrix, mos, vcvec, norb, alpha=1.0_dp, beta=0.0_dp)
CALL cp_dbcsr_sm_fm_multiply(matrix_hfx_asymm(ispin)%matrix, evect(ispin), &
cpmos(ispin), norb, alpha=focc, beta=1.0_dp)
CALL cp_dbcsr_sm_fm_multiply(matrix_hfx_asymm(ispin)%matrix, mos, vcvec, norb, alpha=1.0_dp, beta=1.0_dp)
!
CALL parallel_gemm("T", "N", norb, norb, nao, 1.0_dp, mos, vcvec, 0.0_dp, cvcmat)
CALL parallel_gemm("N", "N", nao, norb, norb, 1.0_dp, evect(ispin), cvcmat, 0.0_dp, vcvec)
CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, vcvec, cpmos(ispin), &
norb, alpha=-focc, beta=1.0_dp)
!
CALL cp_dbcsr_plus_fm_fm_t(matrix_wx1(ispin)%matrix, matrix_v=mos, matrix_g=vcvec, &
ncol=norb, alpha=2.0_dp, symmetry_mode=1)
! For spin-flip excitations this is the only contribution to the alpha response vector
IF (.NOT. (do_sf .AND. (ispin .EQ. 2))) THEN
! F^sym*X
CALL cp_dbcsr_sm_fm_multiply(matrix_hfx(spin)%matrix, evect(spin), &
cpmos(ispin), norb(ispin), alpha=focc, beta=1.0_dp)
! F^asym*X
CALL cp_dbcsr_sm_fm_multiply(matrix_hfx_asymm(spin)%matrix, evect(spin), &
cpmos(ispin), norb(ispin), alpha=focc, beta=1.0_dp)
END IF
! For spin-flip excitations this is the only contribution to the beta response vector
IF (.NOT. (do_sf .AND. (ispin .EQ. 1))) THEN
! F^sym*Cb
CALL cp_dbcsr_sm_fm_multiply(matrix_hfx(spin)%matrix, mos2, vcvec, norb(ispin), alpha=1.0_dp, beta=0.0_dp)
! -F^asym*Cb
CALL cp_dbcsr_sm_fm_multiply(matrix_hfx_asymm(spin)%matrix, mos2, vcvec, norb(ispin), &
alpha=1.0_dp, beta=1.0_dp)
! Ca^T*F*Cb
CALL parallel_gemm("T", "N", norb(spin), norb(ispin), nao, 1.0_dp, mos, vcvec, 0.0_dp, cvcmat)
! X*Ca^T*F*Cb
CALL parallel_gemm("N", "N", nao, norb(ispin), norb(spin), 1.0_dp, evect(spin), cvcmat, 0.0_dp, vcvec)
! -S*X*Ca^T*F*Cb
CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, vcvec, cpmos(ispin), &
norb(ispin), alpha=-focc, beta=1.0_dp)
! Add contributions to the \Lambda_munu for the perturbed overlap matrix term, third term of Eq. 51
! 2X*Ca^T*F*Cb*Cb^T
CALL cp_dbcsr_plus_fm_fm_t(matrix_wx1(ispin)%matrix, matrix_v=vcvec, matrix_g=mos2, &
ncol=norb(ispin), alpha=2.0_dp, symmetry_mode=1)
END IF
END IF
!
CALL cp_fm_release(vcvec)
CALL cp_fm_release(cvcmat)
END DO
IF (.NOT. is_rks_triplets) THEN
IF (.NOT. (is_rks_triplets .OR. do_sf)) THEN
CALL dbcsr_deallocate_matrix_set(matrix_hx)
END IF
IF (ASSOCIATED(ex_env%matrix_wx1)) CALL dbcsr_deallocate_matrix_set(ex_env%matrix_wx1)
ex_env%matrix_wx1 => matrix_wx1
IF (myfun /= xc_none) THEN
IF (.NOT. ((myfun == xc_none) .OR. (tddfpt_control%spinflip == tddfpt_sf_col))) THEN
CALL dbcsr_deallocate_matrix_set(matrix_fx)
CALL dbcsr_deallocate_matrix_set(matrix_gx)
END IF

265
src/qs_tddfpt2_forces.F
View file

@ -49,6 +49,7 @@ MODULE qs_tddfpt2_forces
USE hfx_ri, ONLY: hfx_ri_update_ks
USE hfx_types, ONLY: hfx_type
USE input_constants, ONLY: do_admm_aux_exch_func_none,&
no_sf_tddfpt,&
oe_shift,&
tddfpt_kernel_full,&
tddfpt_kernel_none,&
@ -143,10 +144,14 @@ MODULE qs_tddfpt2_forces
CONTAINS
! **************************************************************************************************
!> \brief Perform TDDFPT gradient calculation.
!> \brief Perform TDDFPT gradient calculation. This routine calculates the response vector R of Eq. 49
!> in J. Chem. Theory Comput. 2022, 18, 4186−4202 (https://doi.org/10.1021/acs.jctc.2c00144)
!> in ex_env%cpmos and a few contributions to the gradient.
!> \param qs_env Quickstep environment
!> \param gs_mos ...
!> \param ex_env ...
!> \param ex_env Holds: Response vector ex_env%cpmos = R
!> Difference density ex_env%matrix_pe = T
!> Matrix ex_env%matrix_hz = H_munu[T]
!> \param kernel_env ...
!> \param sub_env ...
!> \param work_matrices ...
@ -164,7 +169,8 @@ CONTAINS
CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_forces_main'
INTEGER :: handle, ispin, nspins
INTEGER :: handle, ispin, nspins, spin
LOGICAL :: do_sf
TYPE(admm_type), POINTER :: admm_env
TYPE(cp_fm_struct_type), POINTER :: matrix_struct
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_pe_asymm, matrix_pe_symm, &
@ -180,6 +186,11 @@ CONTAINS
nspins = dft_control%nspins
tddfpt_control => dft_control%tddfpt2_control
IF (tddfpt_control%spinflip .EQ. no_sf_tddfpt) THEN
do_sf = .FALSE.
ELSE
do_sf = .TRUE.
END IF
! rhs of linres equation
IF (ASSOCIATED(ex_env%cpmos)) THEN
DO ispin = 1, SIZE(ex_env%cpmos)
@ -188,17 +199,24 @@ CONTAINS
DEALLOCATE (ex_env%cpmos)
END IF
ALLOCATE (ex_env%cpmos(nspins))
! Create and initialize rectangular matrices of nao*occ dimension for alpha and beta R vectors
! for the Z-vector equation system: AZ=-R
DO ispin = 1, nspins
CALL cp_fm_get_info(matrix=ex_env%evect(ispin), matrix_struct=matrix_struct)
CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, matrix_struct=matrix_struct)
CALL cp_fm_create(ex_env%cpmos(ispin), matrix_struct)
CALL cp_fm_set_all(ex_env%cpmos(ispin), 0.0_dp)
END DO
CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s)
NULLIFY (matrix_pe_asymm, matrix_pe_symm)
! Build difference density matrix_pe = X*X^T - (C*X^T*S*X*C^T + (C*X^T*S*X*C^T)^T)/2
!
CALL dbcsr_allocate_matrix_set(ex_env%matrix_pe, nspins)
CALL dbcsr_allocate_matrix_set(matrix_pe_symm, nspins)
CALL dbcsr_allocate_matrix_set(matrix_pe_asymm, nspins)
DO ispin = 1, nspins
! Initialize matrix_pe as a sparse matrix with zeros
ALLOCATE (ex_env%matrix_pe(ispin)%matrix)
CALL dbcsr_create(ex_env%matrix_pe(ispin)%matrix, template=matrix_s(1)%matrix)
CALL dbcsr_copy(ex_env%matrix_pe(ispin)%matrix, matrix_s(1)%matrix)
@ -213,11 +231,17 @@ CONTAINS
matrix_type=dbcsr_type_antisymmetric)
CALL dbcsr_complete_redistribute(ex_env%matrix_pe(ispin)%matrix, matrix_pe_asymm(ispin)%matrix)
CALL tddfpt_resvec1(ex_env%evect(ispin), gs_mos(ispin)%mos_occ, &
matrix_s(1)%matrix, ex_env%matrix_pe(ispin)%matrix)
IF (do_sf) THEN
spin = 1
ELSE
spin = ispin
END IF
! Add difference density to matrix_pe
CALL tddfpt_resvec1(ex_env%evect(spin), gs_mos(spin)%mos_occ, &
matrix_s(1)%matrix, ex_env%matrix_pe(ispin)%matrix, ispin, do_sf)
END DO
!
! ground state ADMM!
! Project the difference density into auxiliary basis for ADMM
IF (dft_control%do_admm) THEN
CALL get_qs_env(qs_env, admm_env=admm_env)
CALL get_admm_env(admm_env, matrix_s_aux_fit=matrix_s_aux_fit)
@ -232,6 +256,7 @@ CONTAINS
END DO
END IF
!
CALL dbcsr_allocate_matrix_set(ex_env%matrix_hz, nspins)
DO ispin = 1, nspins
ALLOCATE (ex_env%matrix_hz(ispin)%matrix)
@ -239,6 +264,7 @@ CONTAINS
CALL dbcsr_copy(ex_env%matrix_hz(ispin)%matrix, matrix_s(1)%matrix)
CALL dbcsr_set(ex_env%matrix_hz(ispin)%matrix, 0.0_dp)
END DO
! Calculate first term of R vector: H_{\mu i\sigma}[T]
IF (dft_control%qs_control%xtb) THEN
CALL tddfpt_resvec2_xtb(qs_env, ex_env%matrix_pe, gs_mos, ex_env%matrix_hz, ex_env%cpmos)
ELSE
@ -246,9 +272,9 @@ CONTAINS
gs_mos, ex_env%matrix_hz, ex_env%cpmos)
END IF
!
CALL dbcsr_allocate_matrix_set(ex_env%matrix_px1, nspins)
CALL dbcsr_allocate_matrix_set(ex_env%matrix_px1_asymm, nspins)
DO ispin = 1, nspins
CALL dbcsr_allocate_matrix_set(ex_env%matrix_px1, SIZE(ex_env%evect, 1))
CALL dbcsr_allocate_matrix_set(ex_env%matrix_px1_asymm, SIZE(ex_env%evect, 1))
DO ispin = 1, SIZE(ex_env%evect, 1)
ALLOCATE (ex_env%matrix_px1(ispin)%matrix)
CALL dbcsr_create(ex_env%matrix_px1(ispin)%matrix, template=matrix_s(1)%matrix)
CALL dbcsr_copy(ex_env%matrix_px1(ispin)%matrix, matrix_s(1)%matrix)
@ -263,9 +289,9 @@ CONTAINS
IF (tddfpt_control%do_admm) THEN
CALL get_qs_env(qs_env, admm_env=admm_env)
CALL get_admm_env(admm_env, matrix_s_aux_fit=matrix_s_aux_fit)
CALL dbcsr_allocate_matrix_set(ex_env%matrix_px1_admm, nspins)
CALL dbcsr_allocate_matrix_set(ex_env%matrix_px1_admm_asymm, nspins)
DO ispin = 1, nspins
CALL dbcsr_allocate_matrix_set(ex_env%matrix_px1_admm, SIZE(ex_env%evect, 1))
CALL dbcsr_allocate_matrix_set(ex_env%matrix_px1_admm_asymm, SIZE(ex_env%evect, 1))
DO ispin = 1, SIZE(ex_env%evect, 1)
ALLOCATE (ex_env%matrix_px1_admm(ispin)%matrix)
CALL dbcsr_create(ex_env%matrix_px1_admm(ispin)%matrix, template=matrix_s_aux_fit(1)%matrix)
CALL dbcsr_copy(ex_env%matrix_px1_admm(ispin)%matrix, matrix_s_aux_fit(1)%matrix)
@ -278,9 +304,9 @@ CONTAINS
ex_env%matrix_px1_admm_asymm(ispin)%matrix)
END DO
END IF
! TDA forces
! TDA forces. Calculates and adds all missing terms for the response vector, Eq. 49.
CALL tddfpt_forces(qs_env, ex_env, gs_mos, kernel_env, sub_env, work_matrices)
! Rotate res vector cpmos into original frame of occupied orbitals
! Rotate res vector cpmos into original frame of occupied orbitals. Why do we need this?
CALL tddfpt_resvec3(qs_env, ex_env%cpmos, work_matrices)
CALL dbcsr_deallocate_matrix_set(matrix_pe_symm)
@ -341,6 +367,7 @@ CONTAINS
ELSE
!
CALL exstate_potential_release(ex_env)
! Build the values of hartree, fock and exchange-correlation potential on the grid
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, &
@ -360,9 +387,14 @@ CONTAINS
END SUBROUTINE tddfpt_forces
! **************************************************************************************************
!> \brief Calculate direct tddft forces
!> \brief Calculate direct tddft forces.
!> J. Chem. Theory Comput. 2022, 18, 7, 4186–4202 (https://doi.org/10.1021/acs.jctc.2c00144)
!> \param qs_env ...
!> \param ex_env ...
!> \param ex_env Holds on exit
!> cpmos = R, Response vector, Eq. 49.
!> matrix_pe = T, Difference density, Eq. 44.
!> matrix_wx1 = CK[D^X]X^T, Third term of Eq. 51.
!> matrix_wz = CX^T(\omegaS - K)XC^T, Last term of Eq. 51.
!> \param gs_mos ...
!> \param kernel_env ...
!> \param sub_env ...
@ -386,7 +418,8 @@ CONTAINS
CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_force_direct'
INTEGER :: handle, iounit, ispin, natom, norb, &
nspins
nspins, spin
LOGICAL :: do_sf
REAL(KIND=dp) :: evalue
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: ftot1, ftot2
REAL(KIND=dp), DIMENSION(3) :: fodeb
@ -418,6 +451,11 @@ CONTAINS
sab_orb=sab_orb, dft_control=dft_control, force=force)
NULLIFY (tddfpt_control)
tddfpt_control => dft_control%tddfpt2_control
IF (tddfpt_control%spinflip .EQ. no_sf_tddfpt) THEN
do_sf = .FALSE.
ELSE
do_sf = .TRUE.
END IF
nspins = dft_control%nspins
IF (debug_forces) THEN
@ -426,27 +464,43 @@ CONTAINS
CALL total_qs_force(ftot1, force, atomic_kind_set)
END IF
! Build last terms of the response vector, Eq. 49, and third term of Lambda_munu, Eq. 51.
! the response vector is in ex_env%cpmos and Lambda is in ex_env%matrix_wx1
CALL tddfpt_kernel_force(qs_env, ex_env, gs_mos, kernel_env, sub_env, work_matrices, debug_forces)
! Overlap matrix
! Overlap matrix, build the Lambda matrix, Eq. 51.
NULLIFY (matrix_wx1, matrix_wz)
CALL dbcsr_allocate_matrix_set(matrix_wz, nspins)
matrix_wx1 => ex_env%matrix_wx1
CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s, matrix_ks=matrix_ks)
NULLIFY (matrix_wz)
CALL dbcsr_allocate_matrix_set(matrix_wz, nspins)
DO ispin = 1, nspins
IF (do_sf) THEN
spin = 1
ELSE
spin = ispin
END IF
! Create and initialize the Lambda matrix as a sparse matrix
ALLOCATE (matrix_wz(ispin)%matrix)
CALL dbcsr_create(matrix=matrix_wz(ispin)%matrix, template=matrix_s(1)%matrix)
CALL cp_dbcsr_alloc_block_from_nbl(matrix_wz(ispin)%matrix, sab_orb)
CALL dbcsr_set(matrix_wz(ispin)%matrix, 0.0_dp)
CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, ncol_global=norb)
CALL cp_dbcsr_plus_fm_fm_t(matrix_wz(ispin)%matrix, matrix_v=evect(ispin), ncol=norb)
evalue = ex_env%evalue
IF (tddfpt_control%oe_corr == oe_shift) THEN
evalue = ex_env%evalue - tddfpt_control%ev_shift
! For spin-flip excitations only the beta component of the Lambda matrix
! contains the excitation energy term
IF (.NOT. (do_sf .AND. (ispin == 1))) THEN
CALL cp_fm_get_info(gs_mos(spin)%mos_occ, ncol_global=norb)
CALL cp_dbcsr_plus_fm_fm_t(matrix_wz(ispin)%matrix, matrix_v=evect(spin), ncol=norb)
evalue = ex_env%evalue
IF (tddfpt_control%oe_corr == oe_shift) THEN
evalue = ex_env%evalue - tddfpt_control%ev_shift
END IF
CALL dbcsr_scale(matrix_wz(ispin)%matrix, evalue)
END IF
! For spin-flip excitations only the alpha component of the Lambda matrix
! contains the occupied MO energy term
IF (.NOT. (do_sf .AND. (ispin == 2))) THEN
CALL calculate_wx_matrix(gs_mos(ispin)%mos_occ, evect(spin), matrix_ks(ispin)%matrix, &
matrix_wz(ispin)%matrix)
END IF
CALL dbcsr_scale(matrix_wz(ispin)%matrix, evalue)
CALL calculate_wx_matrix(gs_mos(ispin)%mos_occ, evect(ispin), matrix_ks(ispin)%matrix, &
matrix_wz(ispin)%matrix)
END DO
IF (nspins == 2) THEN
CALL dbcsr_add(matrix_wz(1)%matrix, matrix_wz(2)%matrix, &
@ -454,6 +508,8 @@ CONTAINS
END IF
NULLIFY (scrm)
IF (debug_forces) fodeb(1:3) = force(1)%overlap(1:3, 1)
! Calculate the force contribution of matrix_xz into the force in ks_env.
! force%overlap = Tr(dS*matrix_wz), last term of Eq. 51.
CALL build_overlap_matrix(ks_env, matrix_s=scrm, &
matrix_name="OVERLAP MATRIX", &
basis_type_a="ORB", basis_type_b="ORB", &
@ -467,7 +523,8 @@ CONTAINS
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Wx*dS ", fodeb
END IF
! Overlap matrix
! Overlap matrix. Build a part of the first term of Lamda, Eq. 51, corresponding to
! the second term of Eq. 41. matrix_wz = C*X^T*(omega*S - K)*X*C^T
CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s, matrix_ks=matrix_ks)
NULLIFY (matrix_wz)
CALL dbcsr_allocate_matrix_set(matrix_wz, nspins)
@ -476,13 +533,22 @@ CONTAINS
CALL dbcsr_create(matrix=matrix_wz(ispin)%matrix, template=matrix_s(1)%matrix)
CALL cp_dbcsr_alloc_block_from_nbl(matrix_wz(ispin)%matrix, sab_orb)
CALL dbcsr_set(matrix_wz(ispin)%matrix, 0.0_dp)
CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, ncol_global=norb)
evalue = ex_env%evalue
IF (tddfpt_control%oe_corr == oe_shift) THEN
evalue = ex_env%evalue - tddfpt_control%ev_shift
! For spin-flip excitations only the alpha component of Lambda has contributions
! of this term, so skip beta
IF (.NOT. (do_sf .AND. (ispin == 2))) THEN
CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, ncol_global=norb)
evalue = ex_env%evalue
IF (tddfpt_control%oe_corr == oe_shift) THEN
evalue = ex_env%evalue - tddfpt_control%ev_shift
END IF
IF (do_sf) THEN
spin = 2
ELSE
spin = ispin
END IF
CALL calculate_xwx_matrix(gs_mos(ispin)%mos_occ, evect(ispin), matrix_s(1)%matrix, &
matrix_ks(spin)%matrix, matrix_wz(ispin)%matrix, evalue)
END IF
CALL calculate_xwx_matrix(gs_mos(ispin)%mos_occ, evect(ispin), matrix_s(1)%matrix, &
matrix_ks(ispin)%matrix, matrix_wz(ispin)%matrix, evalue)
END DO
IF (nspins == 2) THEN
CALL dbcsr_add(matrix_wz(1)%matrix, matrix_wz(2)%matrix, &
@ -503,11 +569,16 @@ CONTAINS
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: xWx*dS ", fodeb
END IF
! Overlap matrix
! Compute force contribution of the first term of Eq. 41 in the first term of Eq. 51
! that was calculated in tddfpt_kernel_force,
! force%overlap = 0.5C*H[T]*C^T
IF (ASSOCIATED(matrix_wx1)) THEN
IF (nspins == 2) THEN
IF (nspins == 2 .AND. .NOT. do_sf) THEN
CALL dbcsr_add(matrix_wx1(1)%matrix, matrix_wx1(2)%matrix, &
alpha_scalar=0.5_dp, beta_scalar=0.5_dp)
ELSE IF (nspins == 2 .AND. do_sf) THEN
CALL dbcsr_add(matrix_wx1(1)%matrix, matrix_wx1(2)%matrix, &
alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
END IF
NULLIFY (scrm)
IF (debug_forces) fodeb(1:3) = force(1)%overlap(1:3, 1)
@ -520,7 +591,7 @@ CONTAINS
IF (debug_forces) THEN
fodeb(1:3) = force(1)%overlap(1:3, 1) - fodeb(1:3)
CALL para_env%sum(fodeb)
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: WK*dS ", fodeb
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: D^XKP*dS ", fodeb
END IF
END IF
@ -538,16 +609,21 @@ CONTAINS
END SUBROUTINE tddfpt_force_direct
! **************************************************************************************************
!> \brief ...
!> \brief Build the spin difference density,
!> matrix_pe = matrix_pe + X*X^T - (C*X^T*S*X*C^T + (C*X^T*S*X*C^T)^T)/2
!> \param evect ...
!> \param mos_occ ...
!> \param matrix_s ...
!> \param matrix_pe ...
!> \param spin ...
!> \param do_sf ...
! **************************************************************************************************
SUBROUTINE tddfpt_resvec1(evect, mos_occ, matrix_s, matrix_pe)
SUBROUTINE tddfpt_resvec1(evect, mos_occ, matrix_s, matrix_pe, spin, do_sf)
TYPE(cp_fm_type), INTENT(IN) :: evect, mos_occ
TYPE(dbcsr_type), POINTER :: matrix_s, matrix_pe
INTEGER, INTENT(IN) :: spin
LOGICAL, INTENT(IN) :: do_sf
CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_resvec1'
@ -558,39 +634,78 @@ CONTAINS
TYPE(cp_logger_type), POINTER :: logger
CALL timeset(routineN, handle)
! X*X^T
CALL cp_fm_get_info(mos_occ, nrow_global=nao, ncol_global=norb)
CALL cp_dbcsr_plus_fm_fm_t(matrix_pe, matrix_v=evect, ncol=norb)
! X^T*S*X
CALL cp_fm_get_info(evect, matrix_struct=fmstruct)
NULLIFY (fmstruct2)
CALL cp_fm_struct_create(fmstruct=fmstruct2, template_fmstruct=fmstruct, &
nrow_global=norb, ncol_global=norb)
CALL cp_fm_create(xxmat, matrix_struct=fmstruct2)
CALL cp_fm_struct_release(fmstruct2)
CALL cp_fm_create(cxmat, matrix_struct=fmstruct)
CALL cp_dbcsr_sm_fm_multiply(matrix_s, evect, cxmat, norb, alpha=1.0_dp, beta=0.0_dp)
CALL parallel_gemm('T', 'N', norb, norb, nao, 1.0_dp, cxmat, evect, 0.0_dp, xxmat)
CALL parallel_gemm('N', 'N', nao, norb, norb, 1.0_dp, mos_occ, xxmat, 0.0_dp, cxmat)
CALL cp_fm_release(xxmat)
! C*C^T*XX
CALL cp_dbcsr_plus_fm_fm_t(matrix_pe, matrix_v=mos_occ, matrix_g=cxmat, &
ncol=norb, alpha=-1.0_dp, symmetry_mode=1)
CALL cp_fm_release(cxmat)
! matrix_pe = X*X^T
IF (.NOT. do_sf .OR. (do_sf .AND. (spin .EQ. 2))) THEN
CALL cp_dbcsr_plus_fm_fm_t(matrix_pe, matrix_v=evect, ncol=norb)
END IF
! matrix_pe = matrix_pe - (C*X^T*S*X*C^T + (C*X^T*S*X*C^T)^T)/2
IF (.NOT. do_sf .OR. (do_sf .AND. (spin .EQ. 1))) THEN
CALL cp_fm_get_info(evect, matrix_struct=fmstruct)
NULLIFY (fmstruct2)
CALL cp_fm_struct_create(fmstruct=fmstruct2, template_fmstruct=fmstruct, &
nrow_global=norb, ncol_global=norb)
CALL cp_fm_create(xxmat, matrix_struct=fmstruct2)
CALL cp_fm_struct_release(fmstruct2)
CALL cp_fm_create(cxmat, matrix_struct=fmstruct)
! S*X
CALL cp_dbcsr_sm_fm_multiply(matrix_s, evect, cxmat, norb, alpha=1.0_dp, beta=0.0_dp)
! (S*X)^T*X
CALL parallel_gemm('T', 'N', norb, norb, nao, 1.0_dp, cxmat, evect, 0.0_dp, xxmat)
! C*X^T*S*X
CALL parallel_gemm('N', 'N', nao, norb, norb, 1.0_dp, mos_occ, xxmat, 0.0_dp, cxmat)
CALL cp_fm_release(xxmat)
! matrix_pe = matrix_pe - (C*(C^T*X^T*S*X)^T + C^T*(C^T*X^T*S*X))/2
CALL cp_dbcsr_plus_fm_fm_t(matrix_pe, matrix_v=mos_occ, matrix_g=cxmat, &
ncol=norb, alpha=-1.0_dp, symmetry_mode=1)
CALL cp_fm_release(cxmat)
END IF
!
! Test for Tr(Pe*S)=0
CALL dbcsr_dot(matrix_pe, matrix_s, tmp)
IF (ABS(tmp) > 1.e-08_dp) THEN
logger => cp_get_default_logger()
IF (logger%para_env%is_source()) THEN
iounit = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
ELSE
iounit = -1
IF (.NOT. do_sf) THEN
IF (ABS(tmp) > 1.e-08_dp) THEN
logger => cp_get_default_logger()
IF (logger%para_env%is_source()) THEN
iounit = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
ELSE
iounit = -1
END IF
CPWARN("Electron count of excitation density matrix is non-zero.")
IF (iounit > 0) THEN
WRITE (iounit, "(T2,A,T61,G20.10)") "Measured electron count is ", tmp
WRITE (iounit, "(T2,A,/)") REPEAT("*", 79)
END IF
END IF
CPWARN("Electron count of excitation density matrix is non-zero.")
IF (iounit > 0) THEN
WRITE (iounit, "(T2,A,T61,G20.10)") "Measured electron count is ", tmp
WRITE (iounit, "(T2,A,/)") REPEAT("*", 79)
ELSE IF (spin .EQ. 1) THEN
IF (ABS(tmp + 1) > 1.e-08_dp) THEN
logger => cp_get_default_logger()
IF (logger%para_env%is_source()) THEN
iounit = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
ELSE
iounit = -1
END IF
CPWARN("Count of occupied occupation number change is not -1.")
IF (iounit > 0) THEN
WRITE (iounit, "(T2,A,T61,G20.10)") "Measured electron count is ", tmp
WRITE (iounit, "(T2,A,/)") REPEAT("*", 79)
END IF
END IF
ELSE IF (spin .EQ. 2) THEN
IF (ABS(tmp - 1) > 1.e-08_dp) THEN
logger => cp_get_default_logger()
IF (logger%para_env%is_source()) THEN
iounit = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
ELSE
iounit = -1
END IF
CPWARN("Count of unoccupied occupation number change is not 1.")
IF (iounit > 0) THEN
WRITE (iounit, "(T2,A,T61,G20.10)") "Measured electron count is ", tmp
WRITE (iounit, "(T2,A,/)") REPEAT("*", 79)
END IF
END IF
END IF
!
@ -634,12 +749,14 @@ CONTAINS
END SUBROUTINE tddfpt_resvec1_admm
! **************************************************************************************************
!> \brief ...
!> \brief Calculates the action of the H operator as in the first term of equation 49 in
!> https://doi.org/10.1021/acs.jctc.2c00144 (J. Chem. Theory Comput. 2022, 18, 4186−4202)
!> cpmos = H_{\mu i\sigma}[matrix_pe]
!> \param qs_env ...
!> \param matrix_pe ...
!> \param matrix_pe Input square matrix with the size of the number of atomic orbitals squared nao^2
!> \param matrix_pe_admm ...
!> \param gs_mos ...
!> \param matrix_hz ...
!> \param matrix_hz Holds H_{\mu\nu\sigma}[matrix_pe] on exit
!> \param cpmos ...
! **************************************************************************************************
SUBROUTINE tddfpt_resvec2(qs_env, matrix_pe, matrix_pe_admm, gs_mos, matrix_hz, cpmos)
@ -907,9 +1024,7 @@ CONTAINS
DEALLOCATE (v_xc_tau)
END IF
IF (dft_control%do_admm) THEN
IF (qs_env%admm_env%aux_exch_func == do_admm_aux_exch_func_none) THEN
! nothing to do
ELSE
IF (qs_env%admm_env%aux_exch_func /= do_admm_aux_exch_func_none) THEN
! add ADMM xc_section_aux terms: f_x[rhoz_ADMM]
CALL get_qs_env(qs_env, admm_env=admm_env)
CALL get_admm_env(admm_env, rho_aux_fit=rho_aux_fit, matrix_s_aux_fit=msaux, &

6
src/qs_tddfpt2_fprint.F
View file

@ -205,7 +205,7 @@ CONTAINS
ex_env%evalue = evals(istate)
CALL cp_fm_release(ex_env%evect)
ALLOCATE (ex_env%evect(nspins))
DO ispin = 1, nspins
DO ispin = 1, SIZE(evects, 1)
CALL cp_fm_get_info(matrix=evects(ispin, 1), matrix_struct=matrix_struct)
CALL cp_fm_create(ex_env%evect(ispin), matrix_struct)
CALL cp_fm_to_fm(evects(ispin, istate), ex_env%evect(ispin))
@ -233,12 +233,14 @@ CONTAINS
CALL response_force(qs_env=qs_env, vh_rspace=ex_env%vh_rspace, &
vxc_rspace=ex_env%vxc_rspace, vtau_rspace=ex_env%vtau_rspace, &
vadmm_rspace=ex_env%vadmm_rspace, matrix_hz=ex_env%matrix_hz, &
matrix_pz=ex_env%matrix_px1, matrix_pz_admm=p_env%p1_admm, &
matrix_pz=p_env%p1, matrix_pz_admm=p_env%p1_admm, &
matrix_wz=p_env%w1, p_env=p_env, ex_env=ex_env, debug=debug_forces)
END IF
CALL p_env_release(p_env)
!
CALL replicate_qs_force(td_force, para_env)
CALL pforce(iwunit, gs_force, atomic_kind_set, natom)
CALL pforce(iwunit, td_force, atomic_kind_set, natom)
CALL sum_qs_force(td_force, gs_force)
CALL pforce(iwunit, td_force, atomic_kind_set, natom)
!

82
src/qs_tddfpt2_methods.F
View file

@ -44,8 +44,8 @@ MODULE qs_tddfpt2_methods
hfx_create
USE input_constants, ONLY: &
do_admm_aux_exch_func_none, do_admm_basis_projection, do_admm_exch_scaling_none, &
do_admm_purify_none, do_potential_truncated, tddfpt_dipole_velocity, tddfpt_kernel_full, &
tddfpt_kernel_none, tddfpt_kernel_stda
do_admm_purify_none, do_potential_truncated, no_sf_tddfpt, tddfpt_dipole_velocity, &
tddfpt_kernel_full, tddfpt_kernel_none, tddfpt_kernel_stda, tddfpt_sf_col, tddfpt_sf_noncol
USE input_section_types, ONLY: section_vals_get,&
section_vals_get_subs_vals,&
section_vals_type,&
@ -133,7 +133,8 @@ MODULE qs_tddfpt2_methods
CONTAINS
! **************************************************************************************************
!> \brief Perform TDDFPT calculation.
!> \brief Perform TDDFPT calculation. If calc_forces then it also builds the response vector for the
!> Z-vector method and calculates some contributions to the force
!> \param qs_env Quickstep environment
!> \param calc_forces ...
!> \param rixs_env ...
@ -155,8 +156,8 @@ CONTAINS
nstate_max, nstates, nvirt, old_state
INTEGER, DIMENSION(maxspins) :: nactive
LOGICAL :: do_admm, do_exck, do_hfx, do_hfxlr, &
do_hfxsr, do_rixs, do_soc, lmult_tmp, &
state_change
do_hfxsr, do_rixs, do_sf, do_soc, &
lmult_tmp, state_change
REAL(kind=dp) :: gsmin, gsval, xsval
REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: evals, ostrength
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
@ -234,6 +235,7 @@ CONTAINS
tddfpt_control%hfxsr_re_int = .TRUE.
tddfpt_control%do_hfxlr = do_hfxlr
tddfpt_control%do_exck = do_exck
do_sf = tddfpt_control%spinflip /= no_sf_tddfpt
IF (tddfpt_control%do_hfxlr) THEN
kernel_section => section_vals_get_subs_vals(tddfpt_section, "XC%HFX_KERNEL%HFXLR")
CALL section_vals_val_get(kernel_section, "RCUT", r_val=tddfpt_control%hfxlr_rcut)
@ -348,9 +350,13 @@ CONTAINS
nstates = tddfpt_control%nstates
!! Too many states can lead to Problems
!! You should be warned if there are more states
!! then occ-virt Combinations!!
!! than occ-virt Combinations!!
CALL cp_fm_get_info(gs_mos(1)%mos_occ, ncol_global=nocc)
CALL cp_fm_get_info(gs_mos(1)%mos_virt, ncol_global=nvirt)
IF (tddfpt_control%spinflip == no_sf_tddfpt) THEN
CALL cp_fm_get_info(gs_mos(1)%mos_virt, ncol_global=nvirt)
ELSE
CALL cp_fm_get_info(gs_mos(2)%mos_virt, ncol_global=nvirt)
END IF
nstate_max = nocc*nvirt
IF (nstates > nstate_max) THEN
CPWARN("NUMBER OF EXCITED STATES COULD LEAD TO PROBLEMS!")
@ -359,7 +365,7 @@ CONTAINS
tddfpt_control%nstates = nstate_max
END IF
CALL tddfpt_create_work_matrices(work_matrices, gs_mos, nstates, do_hfx, do_admm, &
do_hfxlr, do_exck, qs_env, sub_env)
do_hfxlr, do_exck, do_sf, qs_env, sub_env)
! create full_kernel and admm_kernel within tddfpt_energies
kernel_env%full_kernel => full_kernel_env
@ -405,12 +411,17 @@ CONTAINS
NULLIFY (kernel_env%stda_kernel)
END IF
ALLOCATE (evects(nspins, nstates))
IF (do_sf) THEN
! only alpha -> beta excitations are considered in spin-flip TDDFT
ALLOCATE (evects(1, nstates))
ELSE
ALLOCATE (evects(nspins, nstates))
END IF
ALLOCATE (evals(nstates))
ALLOCATE (S_evects(SIZE(evects, 1), nstates))
ALLOCATE (S_evects(nspins, nstates))
DO istate = 1, nstates
DO ispin = 1, nspins
DO ispin = 1, SIZE(evects, 1)
CALL fm_pool_create_fm( &
work_matrices%fm_pool_ao_mo_occ(ispin)%pool, &
S_evects(ispin, istate))
@ -419,7 +430,7 @@ CONTAINS
IF (.NOT. do_soc) THEN
! compute tddfpt excitation energies of multiplicity mult
CALL tddfpt_energies(qs_env, nstates, work_matrices, &
CALL tddfpt_energies(qs_env, nstates, nspins, work_matrices, &
tddfpt_control, logger, tddfpt_print_section, evects, evals, &
gs_mos, tddfpt_section, S_evects, matrix_s, kernel_env, matrix_ks, &
sub_env, ostrength, dipole_op_mos_occ, mult, xc_section, full_kernel_env, &
@ -483,8 +494,8 @@ CONTAINS
ex_env%evalue = evals(my_state)
! excitation vector
CALL cp_fm_release(ex_env%evect)
ALLOCATE (ex_env%evect(nspins))
DO ispin = 1, nspins
ALLOCATE (ex_env%evect(SIZE(evects, 1)))
DO ispin = 1, SIZE(evects, 1)
CALL cp_fm_get_info(matrix=evects(ispin, 1), &
matrix_struct=matrix_struct)
CALL cp_fm_create(ex_env%evect(ispin), matrix_struct)
@ -507,6 +518,7 @@ CONTAINS
my_state, " with excitation energy ", ex_env%evalue*evolt, " eV"
END IF
! Calculate response vector
IF (calc_forces) THEN
CALL tddfpt_forces_main(qs_env, gs_mos, ex_env, kernel_env, &
sub_env, work_matrices)
@ -867,7 +879,16 @@ CONTAINS
IF (tddfpt_control%rks_triplets) THEN
WRITE (log_unit, "(T2,A,T74,A7)") "KERNEL| Spin symmetry of excitations", "Triplet"
ELSE IF (lsd) THEN
WRITE (log_unit, "(T2,A,T69,A12)") "KERNEL| Spin symmetry of excitations", "Unrestricted"
! Spin-conserving excitations where requested
IF (tddfpt_control%spinflip == no_sf_tddfpt) THEN
WRITE (log_unit, "(T2,A,T69,A12)") "KERNEL| Spin symmetry of excitations", "Unrestricted"
! Spin-flip excitations with collinear exchange-correlation kernel requested
ELSE IF (tddfpt_control%spinflip == tddfpt_sf_col) THEN
WRITE (log_unit, "(T2,A,T72,A9)") "KERNEL| Spin flip", "Collinear"
! Spin-flip excitations with noncollinear exchange-correlation kernel requested
ELSE IF (tddfpt_control%spinflip == tddfpt_sf_noncol) THEN
WRITE (log_unit, "(T2,A,T69,A12)") "KERNEL| Spin flip", "Noncollinear"
END IF
ELSE
WRITE (log_unit, "(T2,A,T74,A7)") "KERNEL| Spin symmetry of excitations", "Singlet"
END IF
@ -883,6 +904,7 @@ CONTAINS
!> \brief The energy calculation has been moved to its own subroutine
!> \param qs_env ...
!> \param nstates ...
!> \param nspins ...
!> \param work_matrices ...
!> \param tddfpt_control ...
!> \param logger ...
@ -903,14 +925,14 @@ CONTAINS
!> \param full_kernel_env ...
!> \param kernel_env_admm_aux ...
! **************************************************************************************************
SUBROUTINE tddfpt_energies(qs_env, nstates, work_matrices, &
SUBROUTINE tddfpt_energies(qs_env, nstates, nspins, work_matrices, &
tddfpt_control, logger, tddfpt_print_section, evects, evals, &
gs_mos, tddfpt_section, S_evects, matrix_s, kernel_env, matrix_ks, &
sub_env, ostrength, dipole_op_mos_occ, mult, xc_section, full_kernel_env, &
kernel_env_admm_aux)
TYPE(qs_environment_type), POINTER :: qs_env
INTEGER :: nstates
INTEGER :: nstates, nspins
TYPE(tddfpt_work_matrices) :: work_matrices
TYPE(tddfpt2_control_type), POINTER :: tddfpt_control
TYPE(cp_logger_type), POINTER :: logger
@ -1020,12 +1042,21 @@ CONTAINS
tddfpt_print_section)
END IF
!! Too many states can lead to problems
!! Too many states can lead to Problems
!! You should be warned if there are more states
!! then occ-virt combinations!!
!! than occ-virt Combinations!!
CALL cp_fm_get_info(gs_mos(1)%mos_occ, ncol_global=nocc)
CALL cp_fm_get_info(gs_mos(1)%mos_virt, ncol_global=nvirt)
IF (tddfpt_control%spinflip == no_sf_tddfpt) THEN
CALL cp_fm_get_info(gs_mos(1)%mos_virt, ncol_global=nvirt)
ELSE
CALL cp_fm_get_info(gs_mos(2)%mos_virt, ncol_global=nvirt)
END IF
nstate_max = nocc*nvirt
IF ((SIZE(gs_mos, 1) == 2) .AND. (tddfpt_control%spinflip == no_sf_tddfpt)) THEN
CALL cp_fm_get_info(gs_mos(2)%mos_occ, ncol_global=nocc)
CALL cp_fm_get_info(gs_mos(2)%mos_virt, ncol_global=nvirt)
nstate_max = nocc*nvirt + nstate_max
END IF
IF (nstates > nstate_max) THEN
CPWARN("NUMBER OF EXCITED STATES COULD LEAD TO PROBLEMS!")
CPWARN("Experimental: CHANGED NSTATES TO ITS MAXIMUM VALUE!")
@ -1056,7 +1087,7 @@ CONTAINS
! according to their excitation energies
log_unit = cp_print_key_unit_nr(logger, tddfpt_print_section, &
"GUESS_VECTORS", extension=".tddfptLog")
CALL tddfpt_guess_vectors(evects=evects, evals=evals, &
CALL tddfpt_guess_vectors(tddfpt_control, nspins, evects=evects, evals=evals, &
gs_mos=gs_mos, log_unit=log_unit)
CALL cp_print_key_finished_output(log_unit, logger, &
tddfpt_print_section, "GUESS_VECTORS")
@ -1200,6 +1231,7 @@ CONTAINS
evals, &
gs_mos, &
matrix_s(1)%matrix, &
tddfpt_control%spinflip, &
min_amplitude=tddfpt_control%min_excitation_amplitude)
CALL tddfpt_print_nto_analysis(qs_env, &
evects, evals, &
@ -1287,6 +1319,7 @@ CONTAINS
INTEGER :: handle, ispin, istate, log_unit, mult, &
nspins
LOGICAL :: do_sf
REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: evals_mult, ostrength_mult
TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: evects_mult
@ -1300,6 +1333,7 @@ CONTAINS
nspins = SIZE(gs_mos)
ALLOCATE (evects_mult(nspins, nstates))
ALLOCATE (evals_mult(nstates))
do_sf = tddfpt_control%spinflip /= no_sf_tddfpt
! First multiplicity
IF (lmult_tmp) THEN
@ -1320,7 +1354,7 @@ CONTAINS
mult = 3
END IF
CALL tddfpt_energies(qs_env, nstates, work_matrices, tddfpt_control, logger, &
CALL tddfpt_energies(qs_env, nstates, nspins, work_matrices, tddfpt_control, logger, &
tddfpt_print_section, evects_mult, evals_mult, &
gs_mos, tddfpt_section, S_evects, matrix_s, &
kernel_env, matrix_ks, sub_env, ostrength_mult, &
@ -1334,7 +1368,7 @@ CONTAINS
CALL tddfpt_release_work_matrices(work_matrices, sub_env)
CALL tddfpt_create_work_matrices(work_matrices, gs_mos, nstates, tddfpt_control%do_hfx, &
tddfpt_control%do_admm, tddfpt_control%do_hfxlr, &
tddfpt_control%do_exck, qs_env, sub_env)
tddfpt_control%do_exck, do_sf, qs_env, sub_env)
END IF
DO istate = 1, nstates
@ -1376,7 +1410,7 @@ CONTAINS
mult = 1
END IF
CALL tddfpt_energies(qs_env, nstates, work_matrices, tddfpt_control, logger, &
CALL tddfpt_energies(qs_env, nstates, nspins, work_matrices, tddfpt_control, logger, &
tddfpt_print_section, evects, evals, &
gs_mos, tddfpt_section, S_evects, matrix_s, &
kernel_env, matrix_ks, sub_env, ostrength, &

51
src/qs_tddfpt2_operators.F
View file

@ -9,6 +9,7 @@ MODULE qs_tddfpt2_operators
USE admm_types, ONLY: admm_type
USE cell_types, ONLY: cell_type,&
pbc
USE cp_control_types, ONLY: tddfpt2_control_type
USE cp_dbcsr_api, ONLY: &
dbcsr_create, dbcsr_filter, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_p_type, &
@ -28,6 +29,8 @@ MODULE qs_tddfpt2_operators
USE hartree_local_types, ONLY: hartree_local_type
USE hfx_admm_utils, ONLY: tddft_hfx_matrix
USE hfx_types, ONLY: hfx_type
USE input_constants, ONLY: no_sf_tddfpt,&
tddfpt_sf_noncol
USE input_section_types, ONLY: section_vals_get,&
section_vals_get_subs_vals,&
section_vals_type
@ -95,23 +98,25 @@ CONTAINS
!> \param gs_mos molecular orbitals optimised for the ground state (only occupied orbital
!> energies [component %evals_occ] are needed)
!> \param matrix_ks Kohn-Sham matrix
!> \param tddfpt_control ...
!> \par History
!> * 05.2016 initialise all matrix elements in one go [Sergey Chulkov]
!> * 03.2017 renamed from tddfpt_init_energy_diff(), altered prototype [Sergey Chulkov]
!> \note Based on the subroutine p_op_l1() which was originally created by
!> Thomas Chassaing on 08.2002.
! **************************************************************************************************
SUBROUTINE tddfpt_apply_energy_diff(Aop_evects, evects, S_evects, gs_mos, matrix_ks)
SUBROUTINE tddfpt_apply_energy_diff(Aop_evects, evects, S_evects, gs_mos, matrix_ks, tddfpt_control)
TYPE(cp_fm_type), DIMENSION(:, :), INTENT(INOUT) :: Aop_evects
TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: evects, S_evects
TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
INTENT(in) :: gs_mos
TYPE(dbcsr_p_type), DIMENSION(:), INTENT(in) :: matrix_ks
TYPE(tddfpt2_control_type), INTENT(in), POINTER :: tddfpt_control
CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_apply_energy_diff'
INTEGER :: handle, ispin, ivect, nactive, nao, &
nspins, nvects
nspins, nvects, spin2
TYPE(cp_fm_struct_type), POINTER :: matrix_struct
TYPE(cp_fm_type) :: hevec
@ -120,13 +125,19 @@ CONTAINS
nspins = SIZE(evects, 1)
nvects = SIZE(evects, 2)
DO ispin = 1, nspins
DO ispin = 1, SIZE(evects, 1)
CALL cp_fm_get_info(matrix=evects(ispin, 1), matrix_struct=matrix_struct, &
nrow_global=nao, ncol_global=nactive)
CALL cp_fm_create(hevec, matrix_struct)
IF (tddfpt_control%spinflip == no_sf_tddfpt) THEN
spin2 = ispin
ELSE
spin2 = 2
END IF
DO ivect = 1, nvects
CALL cp_dbcsr_sm_fm_multiply(matrix_ks(ispin)%matrix, evects(ispin, ivect), &
CALL cp_dbcsr_sm_fm_multiply(matrix_ks(spin2)%matrix, evects(ispin, ivect), &
Aop_evects(ispin, ivect), ncol=nactive, &
alpha=1.0_dp, beta=1.0_dp)
@ -267,9 +278,10 @@ CONTAINS
!> \param work_v_xc work real-space grid to store the gradient of the exchange-correlation
!> potential with respect to the response density (modified on exit)
!> \param work_v_xc_tau ...
!> \param spinflip ...
! **************************************************************************************************
SUBROUTINE tddfpt_apply_xc(A_ia_rspace, kernel_env, rho_ia_struct, is_rks_triplets, &
pw_env, work_v_xc, work_v_xc_tau)
pw_env, work_v_xc, work_v_xc_tau, spinflip)
TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(INOUT) :: A_ia_rspace
TYPE(full_kernel_env_type), INTENT(IN) :: kernel_env
@ -277,14 +289,17 @@ CONTAINS
LOGICAL, INTENT(in) :: is_rks_triplets
TYPE(pw_env_type), POINTER :: pw_env
TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: work_v_xc, work_v_xc_tau
INTEGER, INTENT(in), OPTIONAL :: spinflip
INTEGER :: ispin, nspins
INTEGER :: ispin, my_spinflip, nspins
nspins = SIZE(A_ia_rspace)
my_spinflip = 0
IF (PRESENT(spinflip)) my_spinflip = spinflip
IF (kernel_env%deriv2_analytic) THEN
CALL tddfpt_apply_xc_analytic(kernel_env, rho_ia_struct, is_rks_triplets, nspins, &
pw_env, work_v_xc, work_v_xc_tau)
CALL tddfpt_apply_xc_analytic(kernel_env, rho_ia_struct, is_rks_triplets, my_spinflip, &
nspins, pw_env, work_v_xc, work_v_xc_tau)
ELSE
CALL tddfpt_apply_xc_fd(kernel_env, rho_ia_struct, is_rks_triplets, nspins, &
pw_env, work_v_xc, work_v_xc_tau)
@ -343,6 +358,7 @@ CONTAINS
!> \param rho_ia_struct response density for the given trial vector
!> \param is_rks_triplets indicates that the triplet excited states calculation using
!> spin-unpolarised molecular orbitals has been requested
!> \param spinflip ...
!> \param nspins ...
!> \param pw_env plain wave environment
!> \param work_v_xc work real-space grid to store the gradient of the exchange-correlation
@ -357,18 +373,19 @@ CONTAINS
!> \note Based on the subroutine kpp1_calc_k_p_p1() which was originally created by
!> Mohamed Fawzi on 10.2002.
! **************************************************************************************************
SUBROUTINE tddfpt_apply_xc_analytic(kernel_env, rho_ia_struct, is_rks_triplets, nspins, &
pw_env, work_v_xc, work_v_xc_tau)
SUBROUTINE tddfpt_apply_xc_analytic(kernel_env, rho_ia_struct, is_rks_triplets, spinflip, &
nspins, pw_env, work_v_xc, work_v_xc_tau)
TYPE(full_kernel_env_type), INTENT(in) :: kernel_env
TYPE(qs_rho_type), POINTER :: rho_ia_struct
LOGICAL, INTENT(in) :: is_rks_triplets
INTEGER, INTENT(in) :: nspins
INTEGER, INTENT(in) :: spinflip, nspins
TYPE(pw_env_type), POINTER :: pw_env
TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: work_v_xc, work_v_xc_tau
CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_apply_xc_analytic'
INTEGER :: handle, ispin
LOGICAL :: do_spinflip
TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_ia_g, rho_ia_g2
TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_ia_r, rho_ia_r2, tau_ia_r, tau_ia_r2
@ -383,6 +400,13 @@ CONTAINS
CPASSERT(SIZE(rho_ia_r) == nspins)
CPASSERT((.NOT. ASSOCIATED(tau_ia_r)) .OR. SIZE(tau_ia_r) == nspins)
CPASSERT((.NOT. is_rks_triplets) .OR. nspins == 1)
CPASSERT((spinflip .NE. no_sf_tddfpt) .EQV. (nspins == 2))
END IF
IF (spinflip == tddfpt_sf_noncol) THEN
do_spinflip = .TRUE.
ELSE
do_spinflip = .FALSE.
END IF
NULLIFY (tau_ia_r2)
@ -418,7 +442,8 @@ CONTAINS
CALL xc_calc_2nd_deriv_analytical(v_xc=work_v_xc, v_xc_tau=work_v_xc_tau, deriv_set=kernel_env%xc_deriv_set, &
rho_set=kernel_env%xc_rho_set, &
rho1_set=kernel_env%xc_rho1_set, pw_pool=auxbas_pw_pool, &
xc_section=kernel_env%xc_section, gapw=.FALSE., tddfpt_fac=kernel_env%beta)
xc_section=kernel_env%xc_section, gapw=.FALSE., tddfpt_fac=kernel_env%beta, &
spinflip=do_spinflip)
IF (is_rks_triplets) THEN
DEALLOCATE (rho_ia_r2)
@ -554,7 +579,7 @@ CONTAINS
nspins = SIZE(evects, 1)
nvects = SIZE(evects, 2)
IF (nspins > 1) THEN
IF (SIZE(gs_mos) > 1) THEN
alpha = 1.0_dp
ELSE
alpha = 2.0_dp

74
src/qs_tddfpt2_properties.F
View file

@ -57,7 +57,8 @@ MODULE qs_tddfpt2_properties
cp_print_key_should_output,&
cp_print_key_unit_nr
USE cp_realspace_grid_cube, ONLY: cp_pw_to_cube
USE input_constants, ONLY: tddfpt_dipole_berry,&
USE input_constants, ONLY: no_sf_tddfpt,&
tddfpt_dipole_berry,&
tddfpt_dipole_length,&
tddfpt_dipole_velocity
USE input_section_types, ONLY: section_vals_get_subs_vals,&
@ -632,26 +633,29 @@ CONTAINS
!> \param evals TDDFPT eigenvalues
!> \param gs_mos molecular orbitals optimised for the ground state
!> \param matrix_s overlap matrix
!> \param spinflip ...
!> \param min_amplitude the smallest excitation amplitude to print
!> \par History
!> * 05.2016 created as 'tddfpt_print_summary' [Sergey Chulkov]
!> * 08.2018 splited of from 'tddfpt_print_summary' [Sergey Chulkov]
! **************************************************************************************************
SUBROUTINE tddfpt_print_excitation_analysis(log_unit, evects, evals, gs_mos, matrix_s, min_amplitude)
SUBROUTINE tddfpt_print_excitation_analysis(log_unit, evects, evals, gs_mos, matrix_s, spinflip, min_amplitude)
INTEGER, INTENT(in) :: log_unit
TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in) :: evects
REAL(kind=dp), DIMENSION(:), INTENT(in) :: evals
TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
INTENT(in) :: gs_mos
TYPE(dbcsr_type), POINTER :: matrix_s
INTEGER :: spinflip
REAL(kind=dp), INTENT(in) :: min_amplitude
CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_print_excitation_analysis'
CHARACTER(len=5) :: spin_label
CHARACTER(len=5) :: spin_label, spin_label2
INTEGER :: handle, icol, iproc, irow, ispin, &
istate, nao, ncols_local, nrows_local, &
nspins, nstates, state_spin
nspins, nstates, spin2, state_spin, &
state_spin2
INTEGER(kind=int_8) :: iexc, imo_occ, imo_virt, ind, nexcs, &
nexcs_local, nexcs_max_local, &
nmo_virt_occ, nmo_virt_occ_alpha
@ -675,7 +679,7 @@ CONTAINS
CALL timeset(routineN, handle)
nspins = SIZE(evects, 1)
nspins = SIZE(gs_mos, 1)
nstates = SIZE(evects, 2)
do_exc_analysis = min_amplitude < 1.0_dp
@ -705,28 +709,40 @@ CONTAINS
IF (nspins == 1) THEN
state_spin = 1
state_spin2 = 2
spin_label = ' '
spin_label2 = ' '
ELSE IF (spinflip .NE. no_sf_tddfpt) THEN
state_spin = 1
state_spin2 = 2
spin_label = '(alp)'
spin_label2 = '(bet)'
END IF
END IF
ALLOCATE (S_mos_virt(nspins), weights_fm(nspins))
DO ispin = 1, nspins
CALL cp_fm_get_info(gs_mos(ispin)%mos_virt, matrix_struct=fm_struct)
ALLOCATE (S_mos_virt(SIZE(evects, 1)), weights_fm(SIZE(evects, 1)))
DO ispin = 1, SIZE(evects, 1)
IF (spinflip == no_sf_tddfpt) THEN
spin2 = ispin
ELSE
spin2 = 2
END IF
CALL cp_fm_get_info(gs_mos(spin2)%mos_virt, matrix_struct=fm_struct)
CALL cp_fm_create(S_mos_virt(ispin), fm_struct)
CALL cp_dbcsr_sm_fm_multiply(matrix_s, &
gs_mos(ispin)%mos_virt, &
gs_mos(spin2)%mos_virt, &
S_mos_virt(ispin), &
ncol=nmo_virt(ispin), alpha=1.0_dp, beta=0.0_dp)
ncol=nmo_virt(spin2), alpha=1.0_dp, beta=0.0_dp)
NULLIFY (fm_struct)
CALL cp_fm_struct_create(fm_struct, nrow_global=nmo_virt(ispin), ncol_global=nmo_occ(ispin), &
CALL cp_fm_struct_create(fm_struct, nrow_global=nmo_virt(spin2), ncol_global=nmo_occ(ispin), &
context=blacs_env)
CALL cp_fm_create(weights_fm(ispin), fm_struct)
CALL cp_fm_struct_release(fm_struct)
END DO
nexcs_max_local = 0
DO ispin = 1, nspins
DO ispin = 1, SIZE(evects, 1)
CALL cp_fm_get_info(weights_fm(ispin), nrow_local=nrows_local, ncol_local=ncols_local)
nexcs_max_local = nexcs_max_local + INT(nrows_local, int_8)*INT(ncols_local, int_8)
END DO
@ -739,9 +755,14 @@ CONTAINS
! analyse matrix elements locally and transfer only significant
! excitations to the master node for subsequent ordering
DO ispin = 1, nspins
DO ispin = 1, SIZE(evects, 1)
IF (spinflip == no_sf_tddfpt) THEN
spin2 = ispin
ELSE
spin2 = 2
END IF
! compute excitation amplitudes
CALL parallel_gemm('T', 'N', nmo_virt(ispin), nmo_occ(ispin), nao, 1.0_dp, S_mos_virt(ispin), &
CALL parallel_gemm('T', 'N', nmo_virt(spin2), nmo_occ(ispin), nao, 1.0_dp, S_mos_virt(ispin), &
evects(ispin, istate), 0.0_dp, weights_fm(ispin))
CALL cp_fm_get_info(weights_fm(ispin), nrow_local=nrows_local, ncol_local=ncols_local, &
@ -757,12 +778,12 @@ CONTAINS
weights_local(nexcs_local) = local_data(irow, icol)
! index of single excitation (ivirt, iocc, ispin) in compressed form
inds_local(nexcs_local) = nmo_virt_occ + INT(row_indices(irow), int_8) + &
INT(col_indices(icol) - 1, int_8)*nmo_virt8(ispin)
INT(col_indices(icol) - 1, int_8)*nmo_virt8(spin2)
END IF
END DO
END DO
nmo_virt_occ = nmo_virt_occ + nmo_virt8(ispin)*nmo_occ8(ispin)
nmo_virt_occ = nmo_virt_occ + nmo_virt8(spin2)*nmo_occ8(ispin)
END DO
IF (para_env%is_source()) THEN
@ -846,24 +867,33 @@ CONTAINS
WRITE (log_unit, '(T7,I8,F10.5,A)') istate, evals(istate)*evolt, " eV"
! This reinitialization is needed to prevent the intel fortran compiler from introduce
! a bug when using optimization level 3 flag
state_spin = 1
state_spin2 = 1
IF (spinflip .NE. no_sf_tddfpt) THEN
state_spin = 1
state_spin2 = 2
END IF
DO iexc = 1, nexcs
ind = inds_recv(inds(iexc)) - 1
IF (nspins > 1) THEN
IF ((nspins > 1) .AND. (spinflip == no_sf_tddfpt)) THEN
IF (ind < nmo_virt_occ_alpha) THEN
state_spin = 1
spin_label = '(alp)'
spin_label2 = '(alp)'
ELSE
state_spin = 2
state_spin2 = 2
ind = ind - nmo_virt_occ_alpha
spin_label = '(bet)'
spin_label2 = '(bet)'
END IF
END IF
imo_occ = ind/nmo_virt8(state_spin) + 1
imo_virt = MOD(ind, nmo_virt8(state_spin)) + 1
imo_occ = ind/nmo_virt8(state_spin2) + 1
imo_virt = MOD(ind, nmo_virt8(state_spin2)) + 1
WRITE (log_unit, '(T27,I8,1X,A5,T48,I8,1X,A5,T70,F9.6)') imo_occ, spin_label, &
nmo_occ8(state_spin) + imo_virt, spin_label, weights_recv(inds(iexc))
nmo_occ8(state_spin2) + imo_virt, spin_label2, weights_recv(inds(iexc))
END DO
END IF

27
src/qs_tddfpt2_types.F
View file

@ -214,25 +214,27 @@ CONTAINS
!> \param do_admm ...
!> \param do_hfxlr ...
!> \param do_exck ...
!> \param do_sf ...
!> \param qs_env Quickstep environment
!> \param sub_env parallel group environment
!> \par History
!> * 02.2017 created [Sergey Chulkov]
! **************************************************************************************************
SUBROUTINE tddfpt_create_work_matrices(work_matrices, gs_mos, nstates, do_hfx, do_admm, &
do_hfxlr, do_exck, qs_env, sub_env)
do_hfxlr, do_exck, do_sf, qs_env, sub_env)
TYPE(tddfpt_work_matrices), INTENT(out) :: work_matrices
TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
INTENT(in) :: gs_mos
INTEGER, INTENT(in) :: nstates
LOGICAL, INTENT(in) :: do_hfx, do_admm, do_hfxlr, do_exck
LOGICAL, INTENT(in) :: do_hfx, do_admm, do_hfxlr, do_exck, do_sf
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(tddfpt_subgroup_env_type), INTENT(in) :: sub_env
CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_create_work_matrices'
INTEGER :: handle, igroup, ispin, istate, nao, &
nao_aux, natom, ngroups, nspins
INTEGER :: evecs_dim, handle, igroup, ispin, &
istate, nao, nao_aux, natom, ngroups, &
nspins
INTEGER, DIMENSION(maxspins) :: nmo_occ, nmo_virt
TYPE(admm_type), POINTER :: admm_env
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
@ -271,6 +273,11 @@ CONTAINS
NULLIFY (work_matrices%rho_xc_struct_sub)
nspins = SIZE(gs_mos)
IF (do_sf) THEN
evecs_dim = 1
ELSE
evecs_dim = nspins
END IF
CALL get_qs_env(qs_env, blacs_env=blacs_env, matrix_s=matrix_s)
CALL dbcsr_get_info(matrix_s(1)%matrix, nfullrows_total=nao)
@ -313,25 +320,25 @@ CONTAINS
END DO
IF (sub_env%is_split) THEN
DO ispin = 1, nspins
DO ispin = 1, evecs_dim
CALL cp_fm_struct_create(fm_struct_evects(ispin)%struct, nrow_global=nao, &
ncol_global=nmo_occ(ispin), context=sub_env%blacs_env)
END DO
ALLOCATE (work_matrices%evects_sub(nspins, nstates), work_matrices%Aop_evects_sub(nspins, nstates))
ALLOCATE (work_matrices%evects_sub(evecs_dim, nstates), work_matrices%Aop_evects_sub(evecs_dim, nstates))
CALL blacs_env%get(para_env=para_env)
igroup = sub_env%group_distribution(para_env%mepos)
ngroups = sub_env%ngroups
DO istate = ngroups - igroup, nstates, ngroups
DO ispin = 1, nspins
DO ispin = 1, evecs_dim
CALL cp_fm_create(work_matrices%evects_sub(ispin, istate), fm_struct_evects(ispin)%struct)
CALL cp_fm_create(work_matrices%Aop_evects_sub(ispin, istate), fm_struct_evects(ispin)%struct)
END DO
END DO
DO ispin = nspins, 1, -1
DO ispin = evecs_dim, 1, -1
CALL cp_fm_struct_release(fm_struct_evects(ispin)%struct)
END DO
END IF
@ -356,8 +363,8 @@ CONTAINS
! group-specific dbcsr matrices
NULLIFY (work_matrices%A_ia_munu_sub)
CALL dbcsr_allocate_matrix_set(work_matrices%A_ia_munu_sub, nspins)
DO ispin = 1, nspins
CALL dbcsr_allocate_matrix_set(work_matrices%A_ia_munu_sub, evecs_dim)
DO ispin = 1, evecs_dim
CALL dbcsr_init_p(work_matrices%A_ia_munu_sub(ispin)%matrix)
CALL tddfpt_dbcsr_create_by_dist(work_matrices%A_ia_munu_sub(ispin)%matrix, template=matrix_s(1)%matrix, &
dbcsr_dist=sub_env%dbcsr_dist, sab=sub_env%sab_orb)

137
src/qs_tddfpt2_utils.F
View file

@ -38,8 +38,8 @@ MODULE qs_tddfpt2_utils
cp_logger_get_default_io_unit,&
cp_logger_type
USE input_constants, ONLY: &
cholesky_dbcsr, cholesky_inverse, cholesky_off, cholesky_restore, oe_gllb, oe_lb, oe_none, &
oe_saop, oe_shift
cholesky_dbcsr, cholesky_inverse, cholesky_off, cholesky_restore, no_sf_tddfpt, oe_gllb, &
oe_lb, oe_none, oe_saop, oe_shift
USE input_section_types, ONLY: section_vals_create,&
section_vals_get,&
section_vals_get_subs_vals,&
@ -726,12 +726,14 @@ CONTAINS
! **************************************************************************************************
!> \brief Compute the number of possible singly excited states (occ -> virt)
!> \param tddfpt_control ...
!> \param gs_mos occupied and virtual molecular orbitals optimised for the ground state
!> \return the number of possible single excitations
!> \par History
!> * 01.2017 created [Sergey Chulkov]
! **************************************************************************************************
PURE FUNCTION tddfpt_total_number_of_states(gs_mos) RESULT(nstates_total)
PURE FUNCTION tddfpt_total_number_of_states(tddfpt_control, gs_mos) RESULT(nstates_total)
TYPE(tddfpt2_control_type), POINTER :: tddfpt_control
TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
INTENT(in) :: gs_mos
INTEGER(kind=int_8) :: nstates_total
@ -741,11 +743,18 @@ CONTAINS
nstates_total = 0
nspins = SIZE(gs_mos)
DO ispin = 1, nspins
nstates_total = nstates_total + &
SIZE(gs_mos(ispin)%evals_occ, kind=int_8)* &
SIZE(gs_mos(ispin)%evals_virt, kind=int_8)
END DO
! Total number of possible excitations for spin-conserving TDDFT
IF (tddfpt_control%spinflip == no_sf_tddfpt) THEN
DO ispin = 1, nspins
nstates_total = nstates_total + &
SIZE(gs_mos(ispin)%evals_occ, kind=int_8)* &
SIZE(gs_mos(ispin)%evals_virt, kind=int_8)
END DO
! Total number of possible excitations for spin-flip TDDFT
ELSE
nstates_total = SIZE(gs_mos(1)%evals_occ, kind=int_8)* &
SIZE(gs_mos(2)%evals_virt, kind=int_8)
END IF
END FUNCTION tddfpt_total_number_of_states
! **************************************************************************************************
@ -909,6 +918,8 @@ CONTAINS
! **************************************************************************************************
!> \brief Generate missed guess vectors.
!> \param tddfpt_control ...
!> \param nspins ...
!> \param evects guess vectors distributed across all processors (initialised on exit)
!> \param evals guessed transition energies (initialised on exit)
!> \param gs_mos occupied and virtual molecular orbitals optimised for the ground state
@ -927,20 +938,22 @@ CONTAINS
!> a restart file).
!> \endparblock
! **************************************************************************************************
SUBROUTINE tddfpt_guess_vectors(evects, evals, gs_mos, log_unit)
SUBROUTINE tddfpt_guess_vectors(tddfpt_control, nspins, evects, evals, gs_mos, log_unit)
TYPE(tddfpt2_control_type), POINTER :: tddfpt_control
INTEGER :: nspins
TYPE(cp_fm_type), DIMENSION(:, :), INTENT(inout) :: evects
REAL(kind=dp), DIMENSION(:), INTENT(inout) :: evals
TYPE(tddfpt_ground_state_mos), &
DIMENSION(SIZE(evects, 1)), INTENT(in) :: gs_mos
TYPE(tddfpt_ground_state_mos), DIMENSION(nspins), &
INTENT(in) :: gs_mos
INTEGER, INTENT(in) :: log_unit
CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_guess_vectors'
CHARACTER(len=5) :: spin_label
CHARACTER(len=5) :: spin_label1, spin_label2
INTEGER :: handle, imo_occ, imo_virt, ind, ispin, &
istate, jspin, nspins, nstates, &
istate, jspin, nstates, &
nstates_occ_virt_alpha, &
nstates_selected
nstates_selected, spin1, spin2
INTEGER, ALLOCATABLE, DIMENSION(:) :: inds
INTEGER, DIMENSION(maxspins) :: nmo_occ_avail, nmo_occ_selected, &
nmo_virt_selected
@ -950,7 +963,6 @@ CONTAINS
CALL timeset(routineN, handle)
nspins = SIZE(evects, 1)
nstates = SIZE(evects, 2)
IF (debug_this_module) THEN
@ -971,23 +983,41 @@ CONTAINS
! TO DO: the variable 'nstates_selected' should probably be declared as INTEGER(kind=int_8),
! however we need a special version of the subroutine sort() in order to do so
nstates_selected = DOT_PRODUCT(nmo_occ_selected(1:nspins), nmo_virt_selected(1:nspins))
IF (tddfpt_control%spinflip == no_sf_tddfpt) THEN
nstates_selected = DOT_PRODUCT(nmo_occ_selected(1:nspins), nmo_virt_selected(1:nspins))
ELSE
nstates_selected = nmo_occ_selected(1)*nmo_virt_selected(2)
END IF
ALLOCATE (inds(nstates_selected))
ALLOCATE (e_virt_minus_occ(nstates_selected))
istate = 0
DO ispin = 1, nspins
DO imo_occ = 1, nmo_occ_selected(ispin)
! Here imo_occ enumerate Occupied orbitals in inverse order (from the last to the first element)
e_occ = gs_mos(ispin)%evals_occ(nmo_occ_avail(ispin) - imo_occ + 1)
IF (tddfpt_control%spinflip == no_sf_tddfpt) THEN
! Guess for spin-conserving TDDFT
DO ispin = 1, nspins
DO imo_occ = 1, nmo_occ_selected(ispin)
! Here imo_occ enumerate Occupied orbitals in inverse order (from the last to the first element)
e_occ = gs_mos(ispin)%evals_occ(nmo_occ_avail(ispin) - imo_occ + 1)
DO imo_virt = 1, nmo_virt_selected(ispin)
istate = istate + 1
e_virt_minus_occ(istate) = gs_mos(ispin)%evals_virt(imo_virt) - e_occ
DO imo_virt = 1, nmo_virt_selected(ispin)
istate = istate + 1
e_virt_minus_occ(istate) = gs_mos(ispin)%evals_virt(imo_virt) - e_occ
END DO
END DO
END DO
END DO
ELSE
! Guess for spin-flip TDDFT
DO imo_occ = 1, nmo_occ_selected(1)
! Here imo_occ enumerate alpha Occupied orbitals in inverse order (from the last to the first element)
e_occ = gs_mos(1)%evals_occ(nmo_occ_avail(1) - imo_occ + 1)
DO imo_virt = 1, nmo_virt_selected(2)
istate = istate + 1
e_virt_minus_occ(istate) = gs_mos(2)%evals_virt(imo_virt) - e_occ
END DO
END DO
END IF
IF (debug_this_module) THEN
CPASSERT(istate == nstates_selected)
@ -995,12 +1025,27 @@ CONTAINS
CALL sort(e_virt_minus_occ, nstates_selected, inds)
! Labels and spin component for closed-shell
IF (nspins == 1) THEN
ispin = 1
spin_label = ' '
spin1 = 1
spin2 = spin1
spin_label1 = ' '
spin_label2 = spin_label1
! Labels and spin component for spin-flip excitations
ELSE IF (tddfpt_control%spinflip .NE. no_sf_tddfpt) THEN
spin1 = 1
spin2 = 2
spin_label1 = '(alp)'
spin_label2 = '(bet)'
END IF
nstates_occ_virt_alpha = nmo_occ_selected(1)*nmo_virt_selected(1)
IF (tddfpt_control%spinflip == no_sf_tddfpt) THEN
! Calculate maximum number of alpha excitations
nstates_occ_virt_alpha = nmo_occ_selected(1)*nmo_virt_selected(1)
ELSE
! Calculate maximum number of spin-flip excitations
nstates_occ_virt_alpha = nmo_occ_selected(1)*nmo_virt_selected(2)
END IF
IF (log_unit > 0) THEN
WRITE (log_unit, "(1X,A)") "", &
"-------------------------------------------------------------------------------", &
@ -1012,45 +1057,59 @@ CONTAINS
END IF
DO istate = 1, nstates
! Initial guess vector read from restart file
IF (ASSOCIATED(evects(1, istate)%matrix_struct)) THEN
IF (log_unit > 0) &
WRITE (log_unit, '(T7,I8,T28,A19,T60,F14.5)') &
istate, "*** restarted ***", evals(istate)*evolt
! New initial guess vector
ELSE
! Index of excited state - 1
ind = inds(istate) - 1
IF (nspins > 1) THEN
! Labels and spin component for open-shell spin-conserving excitations
IF ((nspins > 1) .AND. (tddfpt_control%spinflip == no_sf_tddfpt)) THEN
IF (ind < nstates_occ_virt_alpha) THEN
ispin = 1
spin_label = '(alp)'
spin1 = 1
spin2 = 1
spin_label1 = '(alp)'
spin_label2 = '(alp)'
ELSE
ispin = 2
ind = ind - nstates_occ_virt_alpha
spin_label = '(bet)'
spin1 = 2
spin2 = 2
spin_label1 = '(bet)'
spin_label2 = '(bet)'
END IF
END IF
imo_occ = nmo_occ_avail(ispin) - ind/nmo_virt_selected(ispin)
imo_virt = MOD(ind, nmo_virt_selected(ispin)) + 1
! Recover index of occupied MO (imo_occ) and unoccupied MO (imo_virt)
! associated to the excited state index (ind+1)
imo_occ = nmo_occ_avail(spin1) - ind/nmo_virt_selected(spin2)
imo_virt = MOD(ind, nmo_virt_selected(spin2)) + 1
! Assign initial guess for excitation energy
evals(istate) = e_virt_minus_occ(istate)
IF (log_unit > 0) &
WRITE (log_unit, '(T7,I8,T24,I8,T37,A5,T45,I8,T54,A5,T60,F14.5)') &
istate, imo_occ, spin_label, nmo_occ_avail(ispin) + imo_virt, spin_label, e_virt_minus_occ(istate)*evolt
istate, imo_occ, spin_label1, nmo_occ_avail(spin2) + imo_virt, spin_label2, e_virt_minus_occ(istate)*evolt
DO jspin = 1, nspins
DO jspin = 1, SIZE(evects, 1)
! .NOT. ASSOCIATED(evects(jspin, istate)%matrix_struct))
CALL cp_fm_create(evects(jspin, istate), fm_struct_evects(jspin)%struct)
CALL cp_fm_set_all(evects(jspin, istate), 0.0_dp)
IF (jspin == ispin) &
CALL cp_fm_to_fm(gs_mos(ispin)%mos_virt, evects(ispin, istate), &
IF (jspin == spin1) &
! Half transform excitation vector to ao space:
! evects_mi = c_ma*X_ai
CALL cp_fm_to_fm(gs_mos(spin2)%mos_virt, evects(spin1, istate), &
ncol=1, source_start=imo_virt, target_start=imo_occ)
END DO
END IF
END DO
IF (log_unit > 0) THEN
WRITE (log_unit, '(/,T7,A,T50,I24)') 'Number of active states:', tddfpt_total_number_of_states(gs_mos)
WRITE (log_unit, '(/,T7,A,T50,I24)') 'Number of active states:', tddfpt_total_number_of_states(tddfpt_control, gs_mos)
WRITE (log_unit, "(1X,A)") &
"-------------------------------------------------------------------------------"
END IF

18
src/qs_vxc_atom.F
View file

@ -435,16 +435,17 @@ CONTAINS
!> \param para_env ...
!> \param do_tddfpt2 New implementation of TDDFT.
!> \param do_triplet ...
!> \param do_sf ...
!> \param kind_set_external ...
! **************************************************************************************************
SUBROUTINE calculate_xc_2nd_deriv_atom(rho_atom_set, rho1_atom_set, qs_env, xc_section, para_env, &
do_tddfpt2, do_triplet, kind_set_external)
do_tddfpt2, do_triplet, do_sf, kind_set_external)
TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set, rho1_atom_set
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(section_vals_type), POINTER :: xc_section
TYPE(mp_para_env_type), INTENT(IN) :: para_env
LOGICAL, INTENT(IN), OPTIONAL :: do_tddfpt2, do_triplet
LOGICAL, INTENT(IN), OPTIONAL :: do_tddfpt2, do_triplet, do_sf
TYPE(qs_kind_type), DIMENSION(:), OPTIONAL, &
POINTER :: kind_set_external
@ -455,8 +456,8 @@ CONTAINS
INTEGER, DIMENSION(2) :: local_loop_limit
INTEGER, DIMENSION(2, 3) :: bounds
INTEGER, DIMENSION(:), POINTER :: atom_list
LOGICAL :: gradient_functional, lsd, paw_atom, &
scale_rho, tau_f
LOGICAL :: gradient_functional, lsd, my_do_sf, &
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
@ -511,6 +512,9 @@ CONTAINS
CALL section_vals_val_get(xc_section, "TAU_CUTOFF", &
r_val=tau_cut)
my_do_sf = .FALSE.
IF (PRESENT(do_sf)) my_do_sf = do_sf
xc_fun_section => section_vals_get_subs_vals(xc_section, &
"XC_FUNCTIONAL")
IF (lsd) THEN
@ -682,11 +686,13 @@ CONTAINS
CALL xc_2nd_deriv_of_r(xc_section=xc_section, &
rho_set=rho_set_h, rho1_set=rho1_set_h, &
deriv_set=deriv_set, &
w=weight, vxc=vxc_h, vxg=vxg_h, do_triplet=do_triplet)
w=weight, vxc=vxc_h, vxg=vxg_h, do_triplet=do_triplet, &
do_sf=my_do_sf)
CALL xc_2nd_deriv_of_r(xc_section=xc_section, &
rho_set=rho_set_s, rho1_set=rho1_set_s, &
deriv_set=deriv_set, &
w=weight, vxc=vxc_s, vxg=vxg_s, do_triplet=do_triplet)
w=weight, vxc=vxc_s, vxg=vxg_s, do_triplet=do_triplet, &
do_sf=my_do_sf)
CALL get_rho_atom(rho_atom=rho1_atom, ga_Vlocal_gb_h=int_hh, ga_Vlocal_gb_s=int_ss)
IF (gradient_functional) THEN

10
src/response_solver.F
View file

@ -645,9 +645,11 @@ CONTAINS
! **************************************************************************************************
!> \brief Initializes vectors for MO-coefficient based linear response solver
!> and calculates response density, and energy-weighted response density matrix
!> J. Chem. Theory Comput. 2022, 18, 4186−4202 (https://doi.org/10.1021/acs.jctc.2c00144)
!>
!> \param qs_env ...
!> \param p_env ...
!> \param p_env Holds the two results of this routine, p_env%p1 = CZ^T + ZC^T,
!> p_env%w1 = 0.5\sum_i(C_i*\epsilon_i*Z_i^T + Z_i*\epsilon_i*C_i^T)
!> \param cpmos RHS of equation as Ax + b = 0 (sign of b)
!> \param iounit ...
!> \param lr_section ...
@ -780,7 +782,7 @@ CONTAINS
END DO
END IF
! Calculate Wz = 0.5*(psi1*eps*psi0^T + psi0*eps*psi1^T)
! Calculate the second term of Eq. 51 Wz = 0.5*(psi1*eps*psi0^T + psi0*eps*psi1^T)
CALL get_qs_env(qs_env, matrix_ks=matrix_ks)
DO ispin = 1, nspins
CALL calculate_wz_matrix(mos(ispin), psi1(ispin), matrix_ks(ispin)%matrix, &
@ -965,7 +967,7 @@ CONTAINS
mpa => matrix_pz
END IF
!
! START OF Tr(P+Z)Hcore
! START OF Tr[(P+Z)Hcore]
!
! Kinetic energy matrix
@ -1060,7 +1062,7 @@ CONTAINS
CALL total_qs_force(ftot2, force, atomic_kind_set)
fodeb(1:3) = ftot2(1:3, 1) - ftot1(1:3, 1)
CALL para_env%sum(fodeb)
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pz*dHcore", fodeb
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: (T+Dz)*dHcore", fodeb
END IF
IF (debug_stress .AND. use_virial) THEN
stdeb = fconv*(virial%pv_virial - sttot)

1072
src/xc/xc.F
View file

File diff suppressed because it is too large Load diff

12
src/xc/xc_atom.F
View file

@ -339,9 +339,10 @@ CONTAINS
!> \param vxc ...
!> \param vxg ...
!> \param do_triplet ...
!> \param do_sf ...
! **************************************************************************************************
SUBROUTINE xc_2nd_deriv_of_r(rho_set, rho1_set, xc_section, &
deriv_set, w, vxc, vxg, do_triplet)
deriv_set, w, vxc, vxg, do_triplet, do_sf)
! As input of this routine one gets rho and drho on a one dimensional grid.
! The grid is the angular grid corresponding to a given point ir on the radial grid.
@ -356,12 +357,12 @@ CONTAINS
REAL(dp), DIMENSION(:, :), POINTER :: w
REAL(dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: vxc
REAL(dp), DIMENSION(:, :, :, :), POINTER :: vxg
LOGICAL, INTENT(IN), OPTIONAL :: do_triplet
LOGICAL, INTENT(IN), OPTIONAL :: do_triplet, do_sf
CHARACTER(LEN=*), PARAMETER :: routineN = 'xc_2nd_deriv_of_r'
INTEGER :: handle, ispin, nspins
LOGICAL :: lsd
LOGICAL :: lsd, my_do_sf
REAL(dp) :: drho_cutoff, my_fac_triplet
TYPE(cp_sll_xc_deriv_type), POINTER :: pos
TYPE(pw_pool_type), POINTER :: pw_pool
@ -381,6 +382,9 @@ CONTAINS
IF (do_triplet) my_fac_triplet = -1.0_dp
END IF
my_do_sf = .FALSE.
IF (PRESENT(do_sf)) my_do_sf = do_sf
CALL xc_rho_set_get(rho_set, drho_cutoff=drho_cutoff)
xc_fun_section => section_vals_get_subs_vals(xc_section, &
"XC_FUNCTIONAL")
@ -409,7 +413,7 @@ CONTAINS
NULLIFY (vxc_tau_pw)
CALL xc_calc_2nd_deriv_analytical(vxc_pw, vxc_tau_pw, deriv_set, rho_set, rho1_set, pw_pool, &
xc_section, gapw=.TRUE., vxg=vxg, tddfpt_fac=my_fac_triplet)
xc_section, gapw=.TRUE., vxg=vxg, tddfpt_fac=my_fac_triplet, spinflip=do_sf)
DEALLOCATE (vxc_pw)

176
src/xc/xc_rho_set_types.F
View file

@ -593,10 +593,10 @@ CONTAINS
!> \param needs the components of rho that are needed
!> \param xc_deriv_method_id ...
!> \param xc_rho_smooth_id ...
!> \param pw_pool pool for the allocation of pw and cr3d
!> \param pw_pool pool for the allocation of pw and array
! **************************************************************************************************
SUBROUTINE xc_rho_set_update(rho_set, rho_r, rho_g, tau, needs, &
xc_deriv_method_id, xc_rho_smooth_id, pw_pool)
xc_deriv_method_id, xc_rho_smooth_id, pw_pool, spinflip)
TYPE(xc_rho_set_type), INTENT(INOUT) :: rho_set
TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(IN) :: rho_r
TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
@ -604,18 +604,22 @@ CONTAINS
TYPE(xc_rho_cflags_type), INTENT(in) :: needs
INTEGER, INTENT(IN) :: xc_deriv_method_id, xc_rho_smooth_id
TYPE(pw_pool_type), POINTER :: pw_pool
LOGICAL, OPTIONAL :: spinflip
REAL(KIND=dp), PARAMETER :: f13 = (1.0_dp/3.0_dp)
INTEGER :: i, idir, ispin, j, k, nspins
LOGICAL :: gradient_f, my_rho_g_local, &
needs_laplace, needs_rho_g
needs_laplace, needs_rho_g, do_sf
REAL(kind=dp) :: rho_cutoff
TYPE(pw_r3d_rs_type), DIMENSION(2) :: laplace_rho_r
TYPE(pw_r3d_rs_type), DIMENSION(3, 2) :: drho_r
TYPE(pw_c1d_gs_type) :: my_rho_g, tmp_g
TYPE(pw_r3d_rs_type), DIMENSION(2) :: my_rho_r
do_sf = .FALSE.
IF (PRESENT(spinflip)) do_sf = spinflip
IF (ANY(rho_set%local_bounds /= pw_pool%pw_grid%bounds_local)) &
CPABORT("pw_pool cr3d have different size than expected")
nspins = SIZE(rho_r)
@ -626,12 +630,20 @@ CONTAINS
! some checks
SELECT CASE (nspins)
CASE (1)
CPASSERT(.NOT. needs%rho_spin)
CPASSERT(.NOT. needs%drho_spin)
CPASSERT(.NOT. needs%norm_drho_spin)
CPASSERT(.NOT. needs%rho_spin_1_3)
CPASSERT(.NOT. needs%tau_spin)
CPASSERT(.NOT. needs%laplace_rho_spin)
IF (.NOT. do_sf) THEN
CPASSERT(.NOT. needs%rho_spin)
CPASSERT(.NOT. needs%drho_spin)
CPASSERT(.NOT. needs%norm_drho_spin)
CPASSERT(.NOT. needs%rho_spin_1_3)
CPASSERT(.NOT. needs%tau_spin)
CPASSERT(.NOT. needs%laplace_rho_spin)
ELSE
CPASSERT(.NOT. needs%rho)
CPASSERT(.NOT. needs%drho)
CPASSERT(.NOT. needs%rho_1_3)
CPASSERT(.NOT. needs%tau)
CPASSERT(.NOT. needs%laplace_rho)
END IF
CASE (2)
CPASSERT(.NOT. needs%rho)
CPASSERT(.NOT. needs%drho)
@ -724,55 +736,113 @@ CONTAINS
SELECT CASE (nspins)
CASE (1)
IF (needs%rho_1_3) THEN
CALL pw_pool%create_cr3d(rho_set%rho_1_3)
!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i,j,k) SHARED(rho_set,my_rho_r)
DO k = rho_set%local_bounds(1, 3), rho_set%local_bounds(2, 3)
DO j = rho_set%local_bounds(1, 2), rho_set%local_bounds(2, 2)
DO i = rho_set%local_bounds(1, 1), rho_set%local_bounds(2, 1)
rho_set%rho_1_3(i, j, k) = MAX(my_rho_r(1)%array(i, j, k), 0.0_dp)**f13
IF (.NOT. do_sf) THEN
IF (needs%rho_1_3) THEN
CALL pw_pool%create_cr3d(rho_set%rho_1_3)
!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i,j,k) SHARED(rho_set,my_rho_r)
DO k = rho_set%local_bounds(1, 3), rho_set%local_bounds(2, 3)
DO j = rho_set%local_bounds(1, 2), rho_set%local_bounds(2, 2)
DO i = rho_set%local_bounds(1, 1), rho_set%local_bounds(2, 1)
rho_set%rho_1_3(i, j, k) = MAX(my_rho_r(1)%array(i, j, k), 0.0_dp)**f13
END DO
END DO
END DO
END DO
rho_set%owns%rho_1_3 = .TRUE.
rho_set%has%rho_1_3 = .TRUE.
END IF
IF (needs%rho) THEN
rho_set%rho => my_rho_r(1)%array
NULLIFY (my_rho_r(1)%array)
rho_set%owns%rho = .TRUE.
rho_set%has%rho = .TRUE.
END IF
IF (needs%norm_drho) THEN
CALL pw_pool%create_cr3d(rho_set%norm_drho)
!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i,j,k) SHARED(rho_set,drho_r)
DO k = rho_set%local_bounds(1, 3), rho_set%local_bounds(2, 3)
DO j = rho_set%local_bounds(1, 2), rho_set%local_bounds(2, 2)
DO i = rho_set%local_bounds(1, 1), rho_set%local_bounds(2, 1)
rho_set%norm_drho(i, j, k) = SQRT( &
drho_r(1, 1)%array(i, j, k)**2 + &
drho_r(2, 1)%array(i, j, k)**2 + &
drho_r(3, 1)%array(i, j, k)**2)
rho_set%owns%rho_1_3 = .TRUE.
rho_set%has%rho_1_3 = .TRUE.
END IF
IF (needs%rho) THEN
rho_set%rho => my_rho_r(1)%array
NULLIFY (my_rho_r(1)%array)
rho_set%owns%rho = .TRUE.
rho_set%has%rho = .TRUE.
END IF
IF (needs%norm_drho) THEN
CALL pw_pool%create_cr3d(rho_set%norm_drho)
!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i,j,k) SHARED(rho_set,drho_r)
DO k = rho_set%local_bounds(1, 3), rho_set%local_bounds(2, 3)
DO j = rho_set%local_bounds(1, 2), rho_set%local_bounds(2, 2)
DO i = rho_set%local_bounds(1, 1), rho_set%local_bounds(2, 1)
rho_set%norm_drho(i, j, k) = SQRT( &
drho_r(1, 1)%array(i, j, k)**2 + &
drho_r(2, 1)%array(i, j, k)**2 + &
drho_r(3, 1)%array(i, j, k)**2)
END DO
END DO
END DO
END DO
rho_set%owns%norm_drho = .TRUE.
rho_set%has%norm_drho = .TRUE.
END IF
IF (needs%laplace_rho) THEN
rho_set%laplace_rho => laplace_rho_r(1)%array
NULLIFY (laplace_rho_r(1)%array)
rho_set%owns%laplace_rho = .TRUE.
rho_set%has%laplace_rho = .TRUE.
END IF
rho_set%owns%norm_drho = .TRUE.
rho_set%has%norm_drho = .TRUE.
END IF
IF (needs%laplace_rho) THEN
rho_set%laplace_rho => laplace_rho_r(1)%array
NULLIFY (laplace_rho_r(1)%array)
rho_set%owns%laplace_rho = .TRUE.
rho_set%has%laplace_rho = .TRUE.
END IF
IF (needs%drho) THEN
DO idir = 1, 3
rho_set%drho(idir)%array => drho_r(idir, 1)%array
NULLIFY (drho_r(idir, 1)%array)
END DO
rho_set%owns%drho = .TRUE.
rho_set%has%drho = .TRUE.
IF (needs%drho) THEN
DO idir = 1, 3
rho_set%drho(idir)%array => drho_r(idir, 1)%array
NULLIFY (drho_r(idir, 1)%array)
END DO
rho_set%owns%drho = .TRUE.
rho_set%has%drho = .TRUE.
END IF
ELSE
IF (needs%norm_drho) THEN
CALL pw_pool%create_cr3d(rho_set%norm_drho)
!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i,j,k) SHARED(rho_set,drho_r)
DO k = rho_set%local_bounds(1, 3), rho_set%local_bounds(2, 3)
DO j = rho_set%local_bounds(1, 2), rho_set%local_bounds(2, 2)
DO i = rho_set%local_bounds(1, 1), rho_set%local_bounds(2, 1)
rho_set%norm_drho(i, j, k) = SQRT( &
drho_r(1, 1)%array(i, j, k)**2 + &
drho_r(2, 1)%array(i, j, k)**2 + &
drho_r(3, 1)%array(i, j, k)**2)
END DO
END DO
END DO
rho_set%owns%norm_drho = .TRUE.
rho_set%has%norm_drho = .TRUE.
END IF
IF (needs%rho_spin) THEN
rho_set%rhoa => my_rho_r(1)%array
NULLIFY (my_rho_r(1)%array)
rho_set%owns%rho_spin = .TRUE.
rho_set%has%rho_spin = .TRUE.
END IF
IF (needs%norm_drho_spin) THEN
CALL pw_pool%create_cr3d(rho_set%norm_drhoa)
!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i,j,k) SHARED(rho_set,drho_r)
DO k = rho_set%local_bounds(1, 3), rho_set%local_bounds(2, 3)
DO j = rho_set%local_bounds(1, 2), rho_set%local_bounds(2, 2)
DO i = rho_set%local_bounds(1, 1), rho_set%local_bounds(2, 1)
rho_set%norm_drhoa(i, j, k) = SQRT( &
drho_r(1, 1)%array(i, j, k)**2 + &
drho_r(2, 1)%array(i, j, k)**2 + &
drho_r(3, 1)%array(i, j, k)**2)
END DO
END DO
END DO
rho_set%owns%norm_drho_spin = .TRUE.
rho_set%has%norm_drho_spin = .TRUE.
END IF
IF (needs%laplace_rho_spin) THEN
rho_set%laplace_rhoa => laplace_rho_r(1)%array
NULLIFY (laplace_rho_r(1)%array)
rho_set%owns%laplace_rho_spin = .TRUE.
rho_set%has%laplace_rho_spin = .TRUE.
END IF
IF (needs%drho_spin) THEN
DO idir = 1, 3
rho_set%drhoa(idir)%array => drho_r(idir, 1)%array
NULLIFY (drho_r(idir, 1)%array)
END DO
rho_set%owns%drho_spin = .TRUE.
rho_set%has%drho_spin = .TRUE.
END IF
END IF
CASE (2)
IF (needs%rho_spin_1_3) THEN

11
tests/QS/regtest-tddfpt-sf-force/TEST_FILES.toml Normal file
View file

@ -0,0 +1,11 @@
# 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
# Check noncollinear kernel for optimizations using LDA, GGA and global Hybrid functional
"h2o_pade_opt.inp" = [{matcher="E_total", tol=1.0E-11, ref=-16.83108497678968}]
"h2o_pbe_opt.inp" = [{matcher="E_total", tol=1.0E-06, ref=-16.93971229185894}]
"h2o_pbe0_opt.inp" = [{matcher="E_total", tol=1.0E-11, ref=-16.92940602177875}]
#
#EOF

92
tests/QS/regtest-tddfpt-sf-force/h2o_pade_opt.inp Normal file
View file

@ -0,0 +1,92 @@
&GLOBAL
PRINT_LEVEL LOW
PROJECT ftest
RUN_TYPE GEO_OPT
&END GLOBAL
&MOTION
&GEO_OPT
MAX_ITER 1
&END GEO_OPT
&END MOTION
&FORCE_EVAL
METHOD Quickstep
&DFT
MULTIPLICITY 3
UKS
&EXCITED_STATES T
STATE 1
&END EXCITED_STATES
&MGRID
CUTOFF 280
REL_CUTOFF 60
&END MGRID
&POISSON
PERIODIC NONE
POISSON_SOLVER MT
&END POISSON
&QS
METHOD GPW
&END QS
&SCF
EPS_SCF 1.0E-8
MAX_SCF 50
SCF_GUESS ATOMIC
&OT
MINIMIZER DIIS
PRECONDITIONER FULL_SINGLE_INVERSE
&END OT
&OUTER_SCF
EPS_SCF 1.0E-8
MAX_SCF 10
&END OUTER_SCF
&END SCF
&XC
&XC_FUNCTIONAL PADE
&END XC_FUNCTIONAL
&END XC
&END DFT
&PRINT
&FORCES
&END FORCES
&END PRINT
&PROPERTIES
&TDDFPT
CONVERGENCE [eV] 1.0e-7
KERNEL FULL
MAX_ITER 50
NSTATES 5
SPINFLIP NONCOLLINEAR
&LINRES
EPS 1.0E-10
EPS_FILTER 1.0E-15
MAX_ITER 20
PRECONDITIONER FULL_SINGLE_INVERSE
&END LINRES
&END TDDFPT
&END PROPERTIES
&SUBSYS
&CELL
ABC [angstrom] 4.0 4.0 4.0
PERIODIC NONE
&END CELL
&COORD
O 0.000000 0.000000 0.000000
H 0.000000 -0.757136 0.580545
H 0.000000 0.757136 0.580545
&END COORD
&KIND H
BASIS_SET DZV-GTH-PADE
POTENTIAL GTH-PADE-q1
&END KIND
&KIND O
BASIS_SET DZVP-GTH-PADE
POTENTIAL GTH-PADE-q6
&END KIND
&TOPOLOGY
&CENTER_COORDINATES
&END CENTER_COORDINATES
&END TOPOLOGY
&END SUBSYS
&END FORCE_EVAL

107
tests/QS/regtest-tddfpt-sf-force/h2o_pbe0_opt.inp Normal file
View file

@ -0,0 +1,107 @@
&GLOBAL
PRINT_LEVEL LOW
PROJECT ftest
RUN_TYPE GEO_OPT
&END GLOBAL
&MOTION
&GEO_OPT
MAX_ITER 1
&END GEO_OPT
&END MOTION
&FORCE_EVAL
METHOD Quickstep
&DFT
MULTIPLICITY 3
UKS
&EXCITED_STATES T
STATE 1
&END EXCITED_STATES
&MGRID
CUTOFF 280
REL_CUTOFF 60
&END MGRID
&POISSON
PERIODIC NONE
POISSON_SOLVER MT
&END POISSON
&QS
METHOD GPW
&END QS
&SCF
EPS_SCF 1.0E-8
MAX_SCF 50
SCF_GUESS ATOMIC
&OT
MINIMIZER DIIS
PRECONDITIONER FULL_SINGLE_INVERSE
&END OT
&OUTER_SCF
EPS_SCF 1.0E-8
MAX_SCF 10
&END OUTER_SCF
&END SCF
&XC
&HF
FRACTION 0.25
&INTERACTION_POTENTIAL
CUTOFF_RADIUS 10.0
POTENTIAL_TYPE TRUNCATED
T_C_G_DATA ./t_c_g.dat
&END INTERACTION_POTENTIAL
&SCREENING
EPS_SCHWARZ 1.0E-12
&END SCREENING
&END HF
&XC_FUNCTIONAL
&PBE
SCALE_C 1.0
SCALE_X 0.75
&END PBE
&END XC_FUNCTIONAL
&END XC
&END DFT
&PRINT
&FORCES
&END FORCES
&END PRINT
&PROPERTIES
&TDDFPT
CONVERGENCE [eV] 1.0e-7
KERNEL FULL
MAX_ITER 50
NSTATES 5
SPINFLIP COLLINEAR
&LINRES
EPS 1.0E-10
EPS_FILTER 1.0E-15
MAX_ITER 20
PRECONDITIONER FULL_SINGLE_INVERSE
&END LINRES
&END TDDFPT
&END PROPERTIES
&SUBSYS
&CELL
ABC [angstrom] 4.0 4.0 4.0
PERIODIC NONE
&END CELL
&COORD
O 0.000000 0.000000 0.000000
H 0.000000 -0.757136 0.580545
H 0.000000 0.757136 0.580545
&END COORD
&KIND H
BASIS_SET DZV-GTH-PADE
POTENTIAL GTH-PADE-q1
&END KIND
&KIND O
BASIS_SET DZVP-GTH-PADE
POTENTIAL GTH-PADE-q6
&END KIND
&TOPOLOGY
&CENTER_COORDINATES
&END CENTER_COORDINATES
&END TOPOLOGY
&END SUBSYS
&END FORCE_EVAL

92
tests/QS/regtest-tddfpt-sf-force/h2o_pbe_opt.inp Normal file
View file

@ -0,0 +1,92 @@
&GLOBAL
PRINT_LEVEL LOW
PROJECT ftest
RUN_TYPE GEO_OPT
&END GLOBAL
&MOTION
&GEO_OPT
MAX_ITER 1
&END GEO_OPT
&END MOTION
&FORCE_EVAL
METHOD Quickstep
&DFT
MULTIPLICITY 3
UKS
&EXCITED_STATES T
STATE 1
&END EXCITED_STATES
&MGRID
CUTOFF 280
REL_CUTOFF 60
&END MGRID
&POISSON
PERIODIC NONE
POISSON_SOLVER MT
&END POISSON
&QS
METHOD GPW
&END QS
&SCF
EPS_SCF 1.0E-8
MAX_SCF 50
SCF_GUESS ATOMIC
&OT
MINIMIZER DIIS
PRECONDITIONER FULL_SINGLE_INVERSE
&END OT
&OUTER_SCF
EPS_SCF 1.0E-8
MAX_SCF 10
&END OUTER_SCF
&END SCF
&XC
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&END XC
&END DFT
&PRINT
&FORCES
&END FORCES
&END PRINT
&PROPERTIES
&TDDFPT
CONVERGENCE [eV] 1.0e-7
KERNEL FULL
MAX_ITER 50
NSTATES 5
SPINFLIP NONCOLLINEAR
&LINRES
EPS 1.0E-10
EPS_FILTER 1.0E-15
MAX_ITER 20
PRECONDITIONER FULL_SINGLE_INVERSE
&END LINRES
&END TDDFPT
&END PROPERTIES
&SUBSYS
&CELL
ABC [angstrom] 4.0 4.0 4.0
PERIODIC NONE
&END CELL
&COORD
O 0.000000 0.000000 0.000000
H 0.000000 -0.757136 0.580545
H 0.000000 0.757136 0.580545
&END COORD
&KIND H
BASIS_SET DZV-GTH-PADE
POTENTIAL GTH-PADE-q1
&END KIND
&KIND O
BASIS_SET DZVP-GTH-PADE
POTENTIAL GTH-PADE-q6
&END KIND
&TOPOLOGY
&CENTER_COORDINATES
&END CENTER_COORDINATES
&END TOPOLOGY
&END SUBSYS
&END FORCE_EVAL

18
tests/QS/regtest-tddfpt-sf/TEST_FILES.toml Normal file
View file

@ -0,0 +1,18 @@
# 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_pade.inp" = [{matcher="M068", tol=1.0E-03, ref=-6.92175}]
"h2o_pbe.inp" = [{matcher="M068", tol=1.0E-03, ref=-6.95678}]
"h2o_pbe0_col.inp" = [{matcher="M068", tol=1.0E-03, ref=11.73433}]
"h2o_pbe0_admm_pbex.inp" = [{matcher="M068", tol=1.0E-03, ref=9.59331}]
"h2o_pbe0_admm_none.inp" = [{matcher="M068", tol=1.0E-03, ref=9.74160}]
#
#EOF

63
tests/QS/regtest-tddfpt-sf/h2o_pade.inp Normal file
View file

@ -0,0 +1,63 @@
&GLOBAL
PRINT_LEVEL LOW
PROJECT h2o_pade
RUN_TYPE ENERGY
&END GLOBAL
&FORCE_EVAL
METHOD Quickstep
&DFT
BASIS_SET_FILE_NAME BASIS_SET
MULTIPLICITY 3
POTENTIAL_FILE_NAME POTENTIAL
UKS
&MGRID
CUTOFF 100
&END MGRID
&QS
&END QS
&SCF
IGNORE_CONVERGENCE_FAILURE
MAX_SCF 10
SCF_GUESS atomic
&END SCF
&XC
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&END XC
&END DFT
&PROPERTIES
&TDDFPT
CONVERGENCE 1.0e-5
MAX_ITER 10
MAX_KV 10
NSTATES 3
SPINFLIP NONCOLLINEAR
&XC
&XC_FUNCTIONAL PADE
&END XC_FUNCTIONAL
&XC_GRID
XC_DERIV SPLINE2_SMOOTH
&END XC_GRID
&END XC
&END TDDFPT
&END PROPERTIES
&SUBSYS
&CELL
ABC 5.0 5.0 5.0
&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
&END COORD
&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

63
tests/QS/regtest-tddfpt-sf/h2o_pbe.inp Normal file
View file

@ -0,0 +1,63 @@
&GLOBAL
PRINT_LEVEL LOW
PROJECT h2o_pbe
RUN_TYPE ENERGY
&END GLOBAL
&FORCE_EVAL
METHOD Quickstep
&DFT
BASIS_SET_FILE_NAME BASIS_SET
MULTIPLICITY 3
POTENTIAL_FILE_NAME POTENTIAL
UKS
&MGRID
CUTOFF 100
&END MGRID
&QS
&END QS
&SCF
IGNORE_CONVERGENCE_FAILURE
MAX_SCF 10
SCF_GUESS atomic
&END SCF
&XC
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&END XC
&END DFT
&PROPERTIES
&TDDFPT
CONVERGENCE 1.0e-5
MAX_ITER 10
MAX_KV 10
NSTATES 3
SPINFLIP NONCOLLINEAR
&XC
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&XC_GRID
XC_DERIV SPLINE2_SMOOTH
&END XC_GRID
&END XC
&END TDDFPT
&END PROPERTIES
&SUBSYS
&CELL
ABC 5.0 5.0 5.0
&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
&END COORD
&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

88
tests/QS/regtest-tddfpt-sf/h2o_pbe0_admm_none.inp Normal file
View file

@ -0,0 +1,88 @@
&GLOBAL
PRINT_LEVEL LOW
PROJECT h2o_pbe0_admm_none
RUN_TYPE ENERGY
&END GLOBAL
&FORCE_EVAL
METHOD Quickstep
&DFT
BASIS_SET_FILE_NAME BASIS_SET
BASIS_SET_FILE_NAME BASIS_ADMM
MULTIPLICITY 3
POTENTIAL_FILE_NAME POTENTIAL
UKS
&AUXILIARY_DENSITY_MATRIX_METHOD
ADMM_PURIFICATION_METHOD NONE
ADMM_TYPE NONE
EXCH_CORRECTION_FUNC NONE
EXCH_SCALING_MODEL NONE
METHOD BASIS_PROJECTION
&END AUXILIARY_DENSITY_MATRIX_METHOD
&MGRID
CUTOFF 100
&END MGRID
&QS
&END QS
&SCF
IGNORE_CONVERGENCE_FAILURE
MAX_SCF 10
SCF_GUESS atomic
&END SCF
&XC
&HF
FRACTION 0.25
&INTERACTION_POTENTIAL
CUTOFF_RADIUS 10.0
POTENTIAL_TYPE TRUNCATED
T_C_G_DATA ./t_c_g.dat
&END INTERACTION_POTENTIAL
&SCREENING
EPS_SCHWARZ 1.0E-12
&END SCREENING
&END HF
&XC_FUNCTIONAL
&PBE
SCALE_C 1.0
SCALE_X 0.75
&END PBE
&END XC_FUNCTIONAL
&END XC
&END DFT
&PROPERTIES
&TDDFPT
CONVERGENCE 1.0e-5
MAX_ITER 10
MAX_KV 10
NSTATES 3
SPINFLIP COLLINEAR
&XC
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&XC_GRID
XC_DERIV SPLINE2_SMOOTH
&END XC_GRID
&END XC
&END TDDFPT
&END PROPERTIES
&SUBSYS
&CELL
ABC 5.0 5.0 5.0
&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
&END COORD
&KIND H
BASIS_SET DZVP-GTH-PADE
BASIS_SET AUX_FIT cFIT3
POTENTIAL GTH-PADE-q1
&END KIND
&KIND O
BASIS_SET DZVP-GTH-PADE
BASIS_SET AUX_FIT cFIT3
POTENTIAL GTH-PADE-q6
&END KIND
&END SUBSYS
&END FORCE_EVAL

88
tests/QS/regtest-tddfpt-sf/h2o_pbe0_admm_pbex.inp Normal file
View file

@ -0,0 +1,88 @@
&GLOBAL
PRINT_LEVEL LOW
PROJECT h2o_pbe0_admm_pbex
RUN_TYPE ENERGY
&END GLOBAL
&FORCE_EVAL
METHOD Quickstep
&DFT
BASIS_SET_FILE_NAME BASIS_SET
BASIS_SET_FILE_NAME BASIS_ADMM
MULTIPLICITY 3
POTENTIAL_FILE_NAME POTENTIAL
UKS
&AUXILIARY_DENSITY_MATRIX_METHOD
ADMM_PURIFICATION_METHOD NONE
ADMM_TYPE NONE
EXCH_CORRECTION_FUNC PBEx
EXCH_SCALING_MODEL NONE
METHOD BASIS_PROJECTION
&END AUXILIARY_DENSITY_MATRIX_METHOD
&MGRID
CUTOFF 100
&END MGRID
&QS
&END QS
&SCF
IGNORE_CONVERGENCE_FAILURE
MAX_SCF 10
SCF_GUESS atomic
&END SCF
&XC
&HF
FRACTION 0.25
&INTERACTION_POTENTIAL
CUTOFF_RADIUS 10.0
POTENTIAL_TYPE TRUNCATED
T_C_G_DATA ./t_c_g.dat
&END INTERACTION_POTENTIAL
&SCREENING
EPS_SCHWARZ 1.0E-12
&END SCREENING
&END HF
&XC_FUNCTIONAL
&PBE
SCALE_C 1.0
SCALE_X 0.75
&END PBE
&END XC_FUNCTIONAL
&END XC
&END DFT
&PROPERTIES
&TDDFPT
CONVERGENCE 1.0e-5
MAX_ITER 10
MAX_KV 10
NSTATES 3
SPINFLIP COLLINEAR
&XC
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&XC_GRID
XC_DERIV SPLINE2_SMOOTH
&END XC_GRID
&END XC
&END TDDFPT
&END PROPERTIES
&SUBSYS
&CELL
ABC 5.0 5.0 5.0
&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
&END COORD
&KIND H
BASIS_SET DZVP-GTH-PADE
BASIS_SET AUX_FIT cFIT3
POTENTIAL GTH-PADE-q1
&END KIND
&KIND O
BASIS_SET DZVP-GTH-PADE
BASIS_SET AUX_FIT cFIT3
POTENTIAL GTH-PADE-q6
&END KIND
&END SUBSYS
&END FORCE_EVAL

78
tests/QS/regtest-tddfpt-sf/h2o_pbe0_col.inp Normal file
View file

@ -0,0 +1,78 @@
&GLOBAL
PRINT_LEVEL LOW
PROJECT h2o_pbe0_col
RUN_TYPE ENERGY
&END GLOBAL
&FORCE_EVAL
METHOD Quickstep
&DFT
BASIS_SET_FILE_NAME BASIS_SET
MULTIPLICITY 3
POTENTIAL_FILE_NAME POTENTIAL
UKS
&MGRID
CUTOFF 100
&END MGRID
&QS
&END QS
&SCF
IGNORE_CONVERGENCE_FAILURE
MAX_SCF 10
SCF_GUESS atomic
&END SCF
&XC
&HF
FRACTION 0.25
&INTERACTION_POTENTIAL
CUTOFF_RADIUS 10.0
POTENTIAL_TYPE TRUNCATED
T_C_G_DATA ./t_c_g.dat
&END INTERACTION_POTENTIAL
&SCREENING
EPS_SCHWARZ 1.0E-12
&END SCREENING
&END HF
&XC_FUNCTIONAL
&PBE
SCALE_C 1.0
SCALE_X 0.75
&END PBE
&END XC_FUNCTIONAL
&END XC
&END DFT
&PROPERTIES
&TDDFPT
CONVERGENCE 1.0e-5
MAX_ITER 10
MAX_KV 10
NSTATES 3
SPINFLIP COLLINEAR
&XC
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&XC_GRID
XC_DERIV SPLINE2_SMOOTH
&END XC_GRID
&END XC
&END TDDFPT
&END PROPERTIES
&SUBSYS
&CELL
ABC 5.0 5.0 5.0
&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
&END COORD
&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

2
tests/TEST_DIRS
View file

@ -3,6 +3,8 @@
# Directories have been reordered according the execution time needed for a gfortran pdbg run using 2 MPI tasks
# in case a new directory is added just add it at the top of the list..
# the order will be regularly checked and modified...
QS/regtest-tddfpt-sf-force libint
QS/regtest-tddfpt-sf libint
xTB/regtest-tblite-gfn1-grad tblite
QS/regtest-spgr-id spglib
HP-DFT