diff --git a/src/ec_efield_local.F b/src/ec_efield_local.F index acb9c818b5..57bf4c1eae 100644 --- a/src/ec_efield_local.F +++ b/src/ec_efield_local.F @@ -25,10 +25,13 @@ MODULE ec_efield_local USE cp_para_types, ONLY: cp_para_env_type USE dbcsr_api, ONLY: dbcsr_add,& dbcsr_copy,& + dbcsr_dot,& dbcsr_get_block_p,& dbcsr_p_type,& dbcsr_set USE ec_env_types, ONLY: energy_correction_type + USE input_constants, ONLY: ec_functional_dc,& + ec_functional_harris USE kinds, ONLY: dp USE orbital_pointers, ONLY: ncoset USE particle_types, ONLY: particle_type @@ -48,6 +51,8 @@ MODULE ec_efield_local USE qs_period_efield_types, ONLY: efield_berry_type,& init_efield_matrices,& set_efield_matrices + USE qs_rho_types, ONLY: qs_rho_get,& + qs_rho_type #include "./base/base_uses.f90" IMPLICIT NONE @@ -160,7 +165,7 @@ CONTAINS npgfb, nsgfa, nsgfb INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb LOGICAL :: found, trans - REAL(dp) :: charge, dab, fdir + REAL(dp) :: charge, dab, ener_field, fdir, tmp REAL(dp), DIMENSION(3) :: ci, fieldpol, ra, rab, rac, rbc, ria REAL(dp), DIMENSION(3, 3) :: forcea, forceb REAL(dp), DIMENSION(:, :), POINTER :: p_block_a, p_block_b, pblock, pmat, work @@ -169,7 +174,7 @@ CONTAINS TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set TYPE(cell_type), POINTER :: cell TYPE(cp_para_env_type), POINTER :: para_env - TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: dipmat, matrix_ks + TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: dipmat, matrix_ks, matrix_p TYPE(dft_control_type), POINTER :: dft_control TYPE(efield_berry_type), POINTER :: efield TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list @@ -183,6 +188,9 @@ CONTAINS TYPE(qs_force_type), DIMENSION(:), POINTER :: force TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set TYPE(qs_kind_type), POINTER :: qs_kind + TYPE(qs_rho_type), POINTER :: rho + +! CALL timeset(routineN, handle) @@ -230,6 +238,27 @@ CONTAINS alpha_scalar=1.0_dp, beta_scalar=fieldpol(idir)) END DO END DO + + ! Handling of electronic efield energy contribution + ! Harris functional: Part of band structure energy T[P_out*F_KS] + ! DC-DFT : needs to be added here + SELECT CASE (ec_env%energy_functional) + CASE (ec_functional_harris) + CASE (ec_functional_dc) + ! Energy + NULLIFY (rho, matrix_p) + CALL get_qs_env(qs_env=qs_env, rho=rho) + CALL qs_rho_get(rho, rho_ao=matrix_p) + DO ispin = 1, SIZE(matrix_p) + DO idir = 1, 3 + CALL dbcsr_dot(matrix_p(ispin)%matrix, dipmat(idir)%matrix, tmp) + ener_field = ener_field + fieldpol(idir)*tmp + END DO + END DO + ec_env%efield_elec = ener_field + CASE DEFAULT + CPABORT("unknown energy correction") + END SELECT END IF ! forces from the efield contribution diff --git a/src/ec_env_types.F b/src/ec_env_types.F index e72990b917..4172979f91 100644 --- a/src/ec_env_types.F +++ b/src/ec_env_types.F @@ -70,8 +70,9 @@ MODULE ec_env_types INTEGER :: mao_iolevel ! energy components REAL(KIND=dp) :: etotal - REAL(KIND=dp) :: eband, exc, ehartree, vhxc - REAL(KIND=dp) :: edispersion, efield_nuclear + REAL(KIND=dp) :: eband, ecore, exc, ehartree, vhxc + REAL(KIND=dp) :: edispersion, efield_elec, & + efield_nuclear, exc_aux_fit ! forces TYPE(qs_force_type), DIMENSION(:), POINTER :: force => Null() ! full neighbor lists and corresponding task list diff --git a/src/ec_environment.F b/src/ec_environment.F index d076998043..2603ddf962 100644 --- a/src/ec_environment.F +++ b/src/ec_environment.F @@ -30,11 +30,12 @@ MODULE ec_environment USE dm_ls_scf_types, ONLY: ls_scf_env_type USE ec_env_types, ONLY: energy_correction_type USE input_constants, ONLY: & - ec_diagonalization, ec_functional_harris, ec_matrix_sign, ec_matrix_tc2, ec_matrix_trs4, & - ec_ot_atomic, ec_ot_diag, ec_ot_gs, kg_cholesky, ls_cluster_atomic, ls_cluster_molecular, & - ls_s_inversion_hotelling, ls_s_inversion_none, ls_s_inversion_sign_sqrt, & - ls_s_preconditioner_atomic, ls_s_preconditioner_molecular, ls_s_preconditioner_none, & - ls_s_sqrt_ns, ls_s_sqrt_proot, xc_vdw_fun_nonloc, xc_vdw_fun_pairpot + ec_diagonalization, ec_functional_dc, ec_functional_harris, ec_matrix_sign, ec_matrix_tc2, & + ec_matrix_trs4, ec_ot_atomic, ec_ot_diag, ec_ot_gs, kg_cholesky, ls_cluster_atomic, & + ls_cluster_molecular, ls_s_inversion_hotelling, ls_s_inversion_none, & + ls_s_inversion_sign_sqrt, ls_s_preconditioner_atomic, ls_s_preconditioner_molecular, & + ls_s_preconditioner_none, ls_s_sqrt_ns, ls_s_sqrt_proot, xc_vdw_fun_nonloc, & + xc_vdw_fun_pairpot USE input_cp2k_check, ONLY: xc_functionals_expand USE input_section_types, ONLY: section_get_ival,& section_vals_get,& @@ -76,10 +77,10 @@ CONTAINS ! ************************************************************************************************** !> \brief Allocates and intitializes ec_env -!> \param qs_env ... -!> \param ec_env the object to create -!> \param dft_section ... -!> \param ec_section Parameters for the energy correction and method to solve it +!> \param qs_env The QS environment +!> \param ec_env The energy correction environment (the object to create) +!> \param dft_section The DFT section +!> \param ec_section The energy correction input section !> \par History !> 2019.09 created !> \author JGH @@ -97,11 +98,11 @@ CONTAINS END SUBROUTINE ec_env_create ! ************************************************************************************************** -!> \brief Initializes ec_env -!> \param qs_env ... -!> \param ec_env ... -!> \param dft_section ... -!> \param ec_section Parameters for the energy correction and method to solve it +!> \brief Initializes energy correction environment +!> \param qs_env The QS environment +!> \param ec_env The energy correction environment +!> \param dft_section The DFT section +!> \param ec_section The energy correction input section !> \par History !> 2019.09 created !> \author JGH @@ -246,10 +247,30 @@ CONTAINS ! CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto, basis_type="HARRIS") CALL init_orbital_pointers(maxlgto + 1) + ! + ! Density-corrected DFT must be performed with the same basis as ground-state + CALL uppercase(ec_env%basis) + ! Basis may only differ from ground-state if explicitly added + IF (ec_env%energy_functional == ec_functional_dc .AND. ec_env%basis == "HARRIS") THEN + DO ikind = 1, nkind + qs_kind => qs_kind_set(ikind) + ! Basis sets of ground-state + CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type="ORB") + ! Basis sets of energy correction + CALL get_qs_kind(qs_kind=qs_kind, basis_set=harris_basis, basis_type="HARRIS") + + IF (basis_set%name .NE. harris_basis%name) THEN + CPABORT("DC-DFT: Correction and ground state need to use the same basis") + END IF + END DO + END IF + ! ! set functional SELECT CASE (ec_env%energy_functional) CASE (ec_functional_harris) ec_env%ec_name = "Harris" + CASE (ec_functional_dc) + ec_env%ec_name = "DC-DFT" CASE DEFAULT CPABORT("unknown energy correction") END SELECT @@ -429,8 +450,8 @@ CONTAINS END SUBROUTINE ec_ls_create ! ************************************************************************************************** -!> \brief Initializes linear scaling environment for LS based solver of -!> Harris energy functional and parses input section +!> \brief Print out the energy correction input section +!> !> \param ec_env ... !> \param unit_nr ... !> \par History @@ -451,119 +472,134 @@ CONTAINS WRITE (unit_nr, '(T2,A)') & "!"//REPEAT("-", 29)//" Energy Correction "//REPEAT("-", 29)//"!" - ! Algorithm - SELECT CASE (ec_env%ks_solver) - CASE (ec_diagonalization) - WRITE (unit_nr, '(T2,A,T61,A20)') "Algorithm: ", "DIAGONALIZATION" - CASE (ec_ot_diag) - WRITE (unit_nr, '(T2,A,T61,A20)') "Algorithm: ", "OT DIAGONALIZATION" - CASE (ec_matrix_sign) - WRITE (unit_nr, '(T2,A,T61,A20)') "Algorithm: ", "MATRIX_SIGN" - CASE (ec_matrix_trs4) - WRITE (unit_nr, '(T2,A,T61,A20)') "Algorithm: ", "TRS4" - CALL cite_reference(Niklasson2003) - CASE (ec_matrix_tc2) - WRITE (unit_nr, '(T2,A,T61,A20)') "Algorithm: ", "TC2" - CALL cite_reference(Niklasson2014) + ! Type of energy correction + SELECT CASE (ec_env%energy_functional) + CASE (ec_functional_harris) + WRITE (unit_nr, '(T2,A,T61,A20)') "Energy Correction: ", "HARRIS FUNCTIONAL" + CASE (ec_functional_dc) + WRITE (unit_nr, '(T2,A,T61,A20)') "Energy Correction: ", "DC-DFT" END SELECT + WRITE (unit_nr, '()') - ! Harris functional parameters + ! Energy correction parameters WRITE (unit_nr, '(T2,A,T61,E20.3)') "eps_default:", ec_env%eps_default CALL uppercase(ec_env%basis) SELECT CASE (ec_env%basis) CASE ("ORBITAL") - WRITE (unit_nr, '(T2,A,T61,A20)') "Harris basis: ", "ORBITAL" + WRITE (unit_nr, '(T2,A,T61,A20)') "EC basis: ", "ORBITAL" CASE ("PRIMITIVE") - WRITE (unit_nr, '(T2,A,T61,A20)') "Harris basis: ", "PRIMITIVE" + WRITE (unit_nr, '(T2,A,T61,A20)') "EC basis: ", "PRIMITIVE" CASE ("HARRIS") - WRITE (unit_nr, '(T2,A,T61,A20)') "Harris Basis: ", "HARRIS" + WRITE (unit_nr, '(T2,A,T61,A20)') "EC Basis: ", "HARRIS" END SELECT - WRITE (unit_nr, '()') - ! MAO - IF (ec_env%mao) THEN - WRITE (unit_nr, '(T2,A,T61,L20)') "MAO:", ec_env%mao - WRITE (unit_nr, '(T2,A,T61,I20)') "MAO_MAX_ITER:", ec_env%mao_max_iter - WRITE (unit_nr, '(T2,A,T61,E20.3)') "MAO_EPS_GRAD:", ec_env%mao_eps_grad - WRITE (unit_nr, '()') - END IF - - ! Parameters depending on chosen algorithm to solve Harris functional - IF (.NOT. ec_env%use_ls_solver) THEN - - WRITE (unit_nr, '(T2,A)') "MO Solver" - WRITE (unit_nr, '()') + ! Parameters for Harris functional solver + IF (ec_env%energy_functional == ec_functional_harris) THEN + ! Algorithm SELECT CASE (ec_env%ks_solver) CASE (ec_diagonalization) - - SELECT CASE (ec_env%factorization) - CASE (kg_cholesky) - WRITE (unit_nr, '(T2,A,T61,A20)') "Factorization: ", "CHOLESKY" - END SELECT - + WRITE (unit_nr, '(T2,A,T61,A20)') "Algorithm: ", "DIAGONALIZATION" CASE (ec_ot_diag) - - ! OT Diagonalization - ! Initial guess : 1) block diagonal initial guess - ! 2) GS-density matrix (might require trafo if basis diff) - - SELECT CASE (ec_env%ec_initial_guess) - CASE (ec_ot_atomic) - WRITE (unit_nr, '(T2,A,T61,A20)') "OT Diag initial guess: ", "ATOMIC" - CASE (ec_ot_gs) - WRITE (unit_nr, '(T2,A,T61,A20)') "OT Diag initial guess: ", "GROUND STATE DM" - END SELECT - - CASE DEFAULT - CPABORT("Unknown Diagonalization algorithm for Harris functional") + WRITE (unit_nr, '(T2,A,T61,A20)') "Algorithm: ", "OT DIAGONALIZATION" + CASE (ec_matrix_sign) + WRITE (unit_nr, '(T2,A,T61,A20)') "Algorithm: ", "MATRIX_SIGN" + CASE (ec_matrix_trs4) + WRITE (unit_nr, '(T2,A,T61,A20)') "Algorithm: ", "TRS4" + CALL cite_reference(Niklasson2003) + CASE (ec_matrix_tc2) + WRITE (unit_nr, '(T2,A,T61,A20)') "Algorithm: ", "TC2" + CALL cite_reference(Niklasson2014) END SELECT - - ELSE - - WRITE (unit_nr, '(T2,A)') "AO Solver" WRITE (unit_nr, '()') - ls_env => ec_env%ls_env - WRITE (unit_nr, '(T2,A,T61,E20.3)') "eps_filter:", ls_env%eps_filter - WRITE (unit_nr, '(T2,A,T61,L20)') "fixed chemical potential (mu)", ls_env%fixed_mu - WRITE (unit_nr, '(T2,A,T61,L20)') "Computing inv(S):", ls_env%needs_s_inv - WRITE (unit_nr, '(T2,A,T61,L20)') "Computing sqrt(S):", ls_env%use_s_sqrt - WRITE (unit_nr, '(T2,A,T61,L20)') "Computing S preconditioner ", ls_env%has_s_preconditioner - WRITE (unit_nr, '(T2,A,T61,L20)') "has unit metric:", ls_env%has_unit_metric - WRITE (unit_nr, '(T2,A,T61,L20)') "Use single precision matrices", & - ls_env%ls_mstruct%single_precision - - IF (ls_env%use_s_sqrt) THEN - SELECT CASE (ls_env%s_sqrt_method) - CASE (ls_s_sqrt_ns) - WRITE (unit_nr, '(T2,A,T61,A20)') "S sqrt method:", "NEWTONSCHULZ" - CASE (ls_s_sqrt_proot) - WRITE (unit_nr, '(T2,A,T61,A20)') "S sqrt method:", "PROOT" - CASE DEFAULT - CPABORT("Unknown sqrt method.") - END SELECT - WRITE (unit_nr, '(T2,A,T61,I20)') "S sqrt order:", ls_env%s_sqrt_order + ! MAO + IF (ec_env%mao) THEN + WRITE (unit_nr, '(T2,A,T61,L20)') "MAO:", ec_env%mao + WRITE (unit_nr, '(T2,A,T61,L20)') "MAO_IOLEVEL:", ec_env%mao_iolevel + WRITE (unit_nr, '(T2,A,T61,I20)') "MAO_MAX_ITER:", ec_env%mao_max_iter + WRITE (unit_nr, '(T2,A,T61,E20.3)') "MAO_EPS_GRAD:", ec_env%mao_eps_grad + WRITE (unit_nr, '(T2,A,T61,E20.3)') "MAO_EPS1:", ec_env%mao_eps1 + WRITE (unit_nr, '()') END IF - SELECT CASE (ls_env%s_preconditioner_type) - CASE (ls_s_preconditioner_none) - WRITE (unit_nr, '(T2,A,T61,A20)') "S preconditioner type ", "NONE" - CASE (ls_s_preconditioner_atomic) - WRITE (unit_nr, '(T2,A,T61,A20)') "S preconditioner type ", "ATOMIC" - CASE (ls_s_preconditioner_molecular) - WRITE (unit_nr, '(T2,A,T61,A20)') "S preconditioner type ", "MOLECULAR" - END SELECT + ! Parameters for linear response solver + IF (.NOT. ec_env%use_ls_solver) THEN - SELECT CASE (ls_env%ls_mstruct%cluster_type) - CASE (ls_cluster_atomic) - WRITE (unit_nr, '(T2,A,T61,A20)') "Cluster type", ADJUSTR("ATOMIC") - CASE (ls_cluster_molecular) - WRITE (unit_nr, '(T2,A,T61,A20)') "Cluster type", ADJUSTR("MOLECULAR") - CASE DEFAULT - CPABORT("Unknown cluster type") - END SELECT + WRITE (unit_nr, '(T2,A)') "MO Solver" + WRITE (unit_nr, '()') + + SELECT CASE (ec_env%ks_solver) + CASE (ec_diagonalization) + + SELECT CASE (ec_env%factorization) + CASE (kg_cholesky) + WRITE (unit_nr, '(T2,A,T61,A20)') "Factorization: ", "CHOLESKY" + END SELECT + + CASE (ec_ot_diag) + + ! OT Diagonalization + ! Initial guess : 1) block diagonal initial guess + ! 2) GS-density matrix (might require trafo if basis diff) + + SELECT CASE (ec_env%ec_initial_guess) + CASE (ec_ot_atomic) + WRITE (unit_nr, '(T2,A,T61,A20)') "OT Diag initial guess: ", "ATOMIC" + CASE (ec_ot_gs) + WRITE (unit_nr, '(T2,A,T61,A20)') "OT Diag initial guess: ", "GROUND STATE DM" + END SELECT + + CASE DEFAULT + CPABORT("Unknown Diagonalization algorithm for Harris functional") + END SELECT + + ELSE + + WRITE (unit_nr, '(T2,A)') "AO Solver" + WRITE (unit_nr, '()') + + ls_env => ec_env%ls_env + WRITE (unit_nr, '(T2,A,T61,E20.3)') "eps_filter:", ls_env%eps_filter + WRITE (unit_nr, '(T2,A,T61,L20)') "fixed chemical potential (mu)", ls_env%fixed_mu + WRITE (unit_nr, '(T2,A,T61,L20)') "Computing inv(S):", ls_env%needs_s_inv + WRITE (unit_nr, '(T2,A,T61,L20)') "Computing sqrt(S):", ls_env%use_s_sqrt + WRITE (unit_nr, '(T2,A,T61,L20)') "Computing S preconditioner ", ls_env%has_s_preconditioner + WRITE (unit_nr, '(T2,A,T61,L20)') "Use single precision matrices", & + ls_env%ls_mstruct%single_precision + + IF (ls_env%use_s_sqrt) THEN + SELECT CASE (ls_env%s_sqrt_method) + CASE (ls_s_sqrt_ns) + WRITE (unit_nr, '(T2,A,T61,A20)') "S sqrt method:", "NEWTONSCHULZ" + CASE (ls_s_sqrt_proot) + WRITE (unit_nr, '(T2,A,T61,A20)') "S sqrt method:", "PROOT" + CASE DEFAULT + CPABORT("Unknown sqrt method.") + END SELECT + WRITE (unit_nr, '(T2,A,T61,I20)') "S sqrt order:", ls_env%s_sqrt_order + END IF + + SELECT CASE (ls_env%s_preconditioner_type) + CASE (ls_s_preconditioner_none) + WRITE (unit_nr, '(T2,A,T61,A20)') "S preconditioner type ", "NONE" + CASE (ls_s_preconditioner_atomic) + WRITE (unit_nr, '(T2,A,T61,A20)') "S preconditioner type ", "ATOMIC" + CASE (ls_s_preconditioner_molecular) + WRITE (unit_nr, '(T2,A,T61,A20)') "S preconditioner type ", "MOLECULAR" + END SELECT + + SELECT CASE (ls_env%ls_mstruct%cluster_type) + CASE (ls_cluster_atomic) + WRITE (unit_nr, '(T2,A,T61,A20)') "Cluster type", ADJUSTR("ATOMIC") + CASE (ls_cluster_molecular) + WRITE (unit_nr, '(T2,A,T61,A20)') "Cluster type", ADJUSTR("MOLECULAR") + CASE DEFAULT + CPABORT("Unknown cluster type") + END SELECT + + END IF END IF diff --git a/src/ec_methods.F b/src/ec_methods.F index 9ab318c0f5..80e3287a7a 100644 --- a/src/ec_methods.F +++ b/src/ec_methods.F @@ -120,7 +120,7 @@ CONTAINS ! folding of second deriv with density in rho1_set CALL xc_calc_2nd_deriv(v_xc=vxc, & ! XC-potential, rho-dependent - v_xc_tau=vxc_tau, & ! XC-potential, tau-dependent + v_xc_tau=vxc_tau, & ! XC-potential, tau-dependent deriv_set=deriv_set, & ! deriv of xc-potential rho_set=rho_set, & ! density at which deriv are calculated rho1_r=rho1_r, & ! density with which to fold diff --git a/src/ec_orth_solver.F b/src/ec_orth_solver.F index c914b6925a..75fc2e2f06 100644 --- a/src/ec_orth_solver.F +++ b/src/ec_orth_solver.F @@ -6,34 +6,33 @@ !--------------------------------------------------------------------------------------------------! ! ************************************************************************************************** -!> \brief AO-based conjugate-graadient response solver routines +!> \brief AO-based conjugate-gradient response solver routines !> !> !> \date 09.2019 !> \author Fabian Belleflamme ! ************************************************************************************************** MODULE ec_orth_solver + USE admm_types, ONLY: admm_type,& + get_admm_env USE cp_control_types, ONLY: dft_control_type USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set,& dbcsr_deallocate_matrix_set USE cp_external_control, ONLY: external_control - USE cp_log_handling, ONLY: cp_get_default_logger,& - cp_logger_get_default_unit_nr,& - cp_logger_type USE cp_para_types, ONLY: cp_para_env_type USE dbcsr_api, ONLY: & dbcsr_add, dbcsr_add_on_diag, dbcsr_checksum, dbcsr_copy, dbcsr_create, & dbcsr_desymmetrize, dbcsr_dot, dbcsr_filter, dbcsr_finalize, dbcsr_get_info, & - dbcsr_multiply, dbcsr_p_type, dbcsr_release, dbcsr_scale, dbcsr_transposed, dbcsr_type, & - dbcsr_type_no_symmetry + dbcsr_multiply, dbcsr_p_type, dbcsr_release, dbcsr_scale, dbcsr_set, dbcsr_transposed, & + dbcsr_type, dbcsr_type_no_symmetry USE ec_methods, ONLY: create_kernel - USE input_constants, ONLY: kg_tnadd_embed,& + USE input_constants, ONLY: do_admm_aux_exch_func_none,& + kg_tnadd_embed,& kg_tnadd_embed_ri,& ls_s_sqrt_ns,& ls_s_sqrt_proot,& precond_mlp - USE input_section_types, ONLY: section_vals_get,& - section_vals_get_subs_vals,& + USE input_section_types, ONLY: section_vals_get_subs_vals,& section_vals_type,& section_vals_val_get USE iterate_matrix, ONLY: matrix_sqrt_Newton_Schulz,& @@ -63,13 +62,17 @@ MODULE ec_orth_solver USE qs_environment_types, ONLY: get_qs_env,& qs_environment_type USE qs_integrate_potential, ONLY: integrate_v_rspace + USE qs_kpp1_env_types, ONLY: qs_kpp1_env_type + USE qs_linres_kernel, ONLY: apply_hfx,& + apply_xc_admm USE qs_linres_types, ONLY: linres_control_type - USE qs_rho_methods, ONLY: qs_rho_rebuild,& - qs_rho_update_rho - USE qs_rho_types, ONLY: qs_rho_create,& - qs_rho_get,& - qs_rho_release,& + USE qs_p_env_methods, ONLY: p_env_check_i_alloc,& + p_env_finish_kpp1,& + p_env_update_rho + USE qs_p_env_types, ONLY: qs_p_env_type + USE qs_rho_types, ONLY: qs_rho_get,& qs_rho_type + USE xc, ONLY: xc_prep_2nd_deriv #include "./base/base_uses.f90" IMPLICIT NONE @@ -339,12 +342,13 @@ CONTAINS END SUBROUTINE preconditioner ! ************************************************************************************************** -!> \brief AO-based conjugate radient linear response solver. +!> \brief AO-based conjugate gradient linear response solver. !> In goes the right hand side B of the equation AZ=B, and the linear transformation of the -!> Hessian matrix A on trial matries is iteratively solved. Result are +!> Hessian matrix A on trial matrices is iteratively solved. Result are !> the response density matrix_pz, and the energy-weighted response density matrix_wz. !> !> \param qs_env ... +!> \param p_env ... !> \param matrix_hz Right hand-side of linear response equation !> \param matrix_pz Response density !> \param matrix_wz Energy-weighted response density matrix @@ -354,9 +358,10 @@ CONTAINS !> \date 01.2020 !> \author Fabian Belleflamme ! ************************************************************************************************** - SUBROUTINE ec_response_ao(qs_env, matrix_hz, matrix_pz, matrix_wz, iounit, should_stop) + SUBROUTINE ec_response_ao(qs_env, p_env, matrix_hz, matrix_pz, matrix_wz, iounit, should_stop) TYPE(qs_environment_type), POINTER :: qs_env + TYPE(qs_p_env_type), POINTER :: p_env TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN), & POINTER :: matrix_hz TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT), & @@ -535,6 +540,7 @@ CONTAINS ! Ax0 CALL build_hessian_op(qs_env=qs_env, & + p_env=p_env, & matrix_ks=matrix_ks, & matrix_p=matrix_p, & ! p matrix_s_sqrt_inv=matrix_s_sqrt_inv, & @@ -646,6 +652,7 @@ CONTAINS ! Hessian Ax = [F,B] + [G(B),P] CALL build_hessian_op(qs_env=qs_env, & + p_env=p_env, & matrix_ks=matrix_ks, & matrix_p=matrix_p, & ! p matrix_s_sqrt_inv=matrix_s_sqrt_inv, & @@ -680,6 +687,7 @@ CONTAINS END IF CALL build_hessian_op(qs_env=qs_env, & + p_env=p_env, & matrix_ks=matrix_ks, & matrix_p=matrix_p, & ! p matrix_s_sqrt_inv=matrix_s_sqrt_inv, & @@ -719,6 +727,7 @@ CONTAINS ! ! r_j+1 = b - A * x_j+1 CALL build_hessian_op(qs_env=qs_env, & + p_env=p_env, & matrix_ks=matrix_ks, & matrix_p=matrix_p, & matrix_s_sqrt_inv=matrix_s_sqrt_inv, & @@ -799,7 +808,7 @@ CONTAINS END DO iteration ! Matrix projector - CALL projector(qs_env, matrix_p, matrix_Ax, eps_filter) + CALL projector(qs_env, matrix_p, matrix_cg_z, eps_filter) ! Z = [cg_z,P] CALL commutator(matrix_cg_z, matrix_p, matrix_z, eps_filter, .TRUE., alpha=0.5_dp) @@ -861,9 +870,8 @@ CONTAINS CHARACTER(len=*), PARAMETER :: routineN = 'ec_wz_matrix', routineP = moduleN//':'//routineN - INTEGER :: handle, iounit, ispin, nspins + INTEGER :: handle, ispin, nspins REAL(KIND=dp) :: scaling - TYPE(cp_logger_type), POINTER :: logger TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_p, matrix_s TYPE(dbcsr_type) :: matrix_tmp, matrix_tmp2 TYPE(dft_control_type), POINTER :: dft_control @@ -881,8 +889,6 @@ CONTAINS matrix_s=matrix_s, & rho=rho) nspins = dft_control%nspins - logger => cp_get_default_logger() - iounit = cp_logger_get_default_unit_nr(logger) CALL qs_rho_get(rho, rho_ao=matrix_p) @@ -978,11 +984,12 @@ CONTAINS END SUBROUTINE hessian_op1 ! ************************************************************************************************** -!> \brief calculate lin transformation of Hessian matrix on a trial vector (matrix) matrix_cg +!> \brief calculate linear transformation of Hessian matrix on a trial matrix matrix_cg !> which is stored in response density B = [cg,P] = cg*P - P*cg = cg*P + (cg*P)^T !> Ax = [F, B] + [G(B), Pin] in orthonormal basis !> !> \param qs_env ... +!> \param p_env ... !> \param matrix_ks Ground-state Kohn-Sham matrix !> \param matrix_p Ground-state Density matrix !> \param matrix_s_sqrt_inv S^(-1/2) needed for transformation to/from orthonormal basis @@ -993,10 +1000,11 @@ CONTAINS !> \date 12.2019 !> \author Fabian Belleflamme ! ************************************************************************************************** - SUBROUTINE build_hessian_op(qs_env, matrix_ks, matrix_p, matrix_s_sqrt_inv, & + SUBROUTINE build_hessian_op(qs_env, p_env, matrix_ks, matrix_p, matrix_s_sqrt_inv, & matrix_cg, matrix_Ax, eps_filter) TYPE(qs_environment_type), POINTER :: qs_env + TYPE(qs_p_env_type), POINTER :: p_env TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN), & POINTER :: matrix_ks, matrix_p TYPE(dbcsr_type), INTENT(IN) :: matrix_s_sqrt_inv @@ -1014,7 +1022,7 @@ CONTAINS TYPE(cp_para_env_type), POINTER :: para_env TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_b, rho1_ao TYPE(dft_control_type), POINTER :: dft_control - TYPE(qs_rho_type), POINTER :: rho, rho1 + TYPE(qs_rho_type), POINTER :: rho CALL timeset(routineN, handle) @@ -1055,13 +1063,13 @@ CONTAINS ! skip the kernel if the DM is very small IF (chksum .GT. 1.0E-14_dp) THEN - CALL get_qs_env(qs_env=qs_env, rho=rho) - ! Bring B as density on grid - NULLIFY (rho1) - CALL qs_rho_create(rho1) - CALL qs_rho_rebuild(rho1, qs_env, rebuild_ao=.TRUE., rebuild_grids=.TRUE.) + ! Bring matrix B as density on grid + + ! prepare perturbation environment + CALL p_env_check_i_alloc(p_env, qs_env) + ! Get response density matrix - CALL qs_rho_get(rho1, rho_ao=rho1_ao) + CALL qs_rho_get(p_env%rho1, rho_ao=rho1_ao) DO ispin = 1, nspins ! Transform B into NON-ortho basis for collocation @@ -1070,22 +1078,21 @@ CONTAINS CALL dbcsr_filter(matrix_b(ispin)%matrix, eps_filter) ! Keep symmetry of density matrix CALL dbcsr_copy(rho1_ao(ispin)%matrix, matrix_b(ispin)%matrix, keep_sparsity=.TRUE.) + CALL dbcsr_copy(p_env%p1(ispin)%matrix, matrix_b(ispin)%matrix, keep_sparsity=.TRUE.) END DO ! Updates densities on grid wrt density matrix - CALL qs_rho_update_rho(rho1, qs_env) + CALL p_env_update_rho(p_env, qs_env) + + DO ispin = 1, nspins + CALL dbcsr_set(p_env%kpp1(ispin)%matrix, 0.0_dp) + IF (ASSOCIATED(p_env%kpp1_admm)) CALL dbcsr_set(p_env%kpp1_admm(ispin)%matrix, 0.0_dp) + END DO ! Calculate kernel ! Ax = F*B - B*F + G(B)*P - P*G(B) - ! IN/OUT IN IN IN - CALL hessian_op2(qs_env, matrix_Ax, matrix_p, matrix_s_sqrt_inv, rho1, eps_filter) - - ! Add contribution from ADMM to kernel -! IF (dft_control%do_admm) THEN -! CALL apply_xc_admm(qs_env, p_env, matrix_Ax, matrix_p, p_env%rho1, eps_filter) -! END IF - - CALL qs_rho_release(rho1) + ! IN/OUT IN IN IN + CALL hessian_op2(qs_env, p_env, matrix_Ax, matrix_p, matrix_s_sqrt_inv, eps_filter) END IF @@ -1096,55 +1103,58 @@ CONTAINS END SUBROUTINE build_hessian_op ! ************************************************************************************************** -!> \brief Calculate lin transformation of Hessian matrix on a trial vector (matrix) matrix_cg +!> \brief Calculate lin transformation of Hessian matrix on a trial matrix matrix_cg !> which is stored in response density B = [cg,P] = cg*P - P*cg = cg*P + (cg*P)^T !> Ax = [F, B] + [G(B), Pin] in orthonormal basis !> !> \param qs_env ... +!> \param p_env p-environment with trial density environment !> \param matrix_Ax contains first part of Hessian linear transformation, kernel contribution !> is calculated and added in this routine !> \param matrix_p Density matrix in orthogonal basis !> \param matrix_s_sqrt_inv contains matrix S^(-1/2) for switching to orthonormal Lowdin basis -!> \param rho1 trial density environment !> \param eps_filter ... !> !> \date 12.2019 !> \author Fabian Belleflamme ! ************************************************************************************************** - SUBROUTINE hessian_op2(qs_env, matrix_Ax, matrix_p, matrix_s_sqrt_inv, rho1, eps_filter) + SUBROUTINE hessian_op2(qs_env, p_env, matrix_Ax, matrix_p, matrix_s_sqrt_inv, eps_filter) TYPE(qs_environment_type), POINTER :: qs_env + TYPE(qs_p_env_type), POINTER :: p_env TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT), & POINTER :: matrix_Ax TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN), & POINTER :: matrix_p TYPE(dbcsr_type), INTENT(IN) :: matrix_s_sqrt_inv - TYPE(qs_rho_type), INTENT(IN), POINTER :: rho1 REAL(KIND=dp), INTENT(IN) :: eps_filter CHARACTER(len=*), PARAMETER :: routineN = 'hessian_op2', routineP = moduleN//':'//routineN INTEGER :: handle, ispin, nspins - LOGICAL :: do_hfx REAL(KIND=dp) :: ekin_mol + TYPE(admm_type), POINTER :: admm_env TYPE(cp_para_env_type), POINTER :: para_env TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_G, matrix_s, rho1_ao, rho_ao TYPE(dft_control_type), POINTER :: dft_control TYPE(pw_env_type), POINTER :: pw_env TYPE(pw_p_type) :: rho_tot_gspace, v_hartree_gspace, & v_hartree_rspace - TYPE(pw_p_type), DIMENSION(:), POINTER :: rho1_g, rho1_r, tau1_r, v_xc, v_xc_tau + TYPE(pw_p_type), DIMENSION(:), POINTER :: rho1_g, rho1_r, rho_r, tau1_r, v_xc, & + v_xc_tau TYPE(pw_poisson_type), POINTER :: poisson_env TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools TYPE(pw_pool_type), POINTER :: auxbas_pw_pool - TYPE(qs_rho_type), POINTER :: rho - TYPE(section_vals_type), POINTER :: hfx_section, input, xc_section + TYPE(qs_kpp1_env_type), POINTER :: kpp1_env + TYPE(qs_rho_type), POINTER :: rho, rho_aux + TYPE(section_vals_type), POINTER :: input, xc_section, xc_section_aux CALL timeset(routineN, handle) - NULLIFY (dft_control, input, matrix_s, para_env, rho, rho1_g, rho1_r) + NULLIFY (admm_env, dft_control, input, matrix_s, para_env, rho, rho_r, rho1_g, rho1_r) CALL get_qs_env(qs_env=qs_env, & + admm_env=admm_env, & dft_control=dft_control, & input=input, & matrix_s=matrix_s, & @@ -1152,10 +1162,14 @@ CONTAINS rho=rho) nspins = dft_control%nspins + CPASSERT(ASSOCIATED(p_env%kpp1)) + CPASSERT(ASSOCIATED(p_env%kpp1_env)) + kpp1_env => p_env%kpp1_env + ! Get non-ortho input density matrix on grid CALL qs_rho_get(rho, rho_ao=rho_ao) - ! Get non-ortho trial density - CALL qs_rho_get(rho1, rho_g=rho1_g, rho_r=rho1_r, tau_r=tau1_r) + ! Get non-ortho trial density stored in p_env + CALL qs_rho_get(p_env%rho1, rho_g=rho1_g, rho_r=rho1_r, tau_r=tau1_r) NULLIFY (pw_env) CALL get_qs_env(qs_env, pw_env=pw_env) @@ -1185,15 +1199,11 @@ CONTAINS ! XC-Kernel NULLIFY (v_xc, v_xc_tau, xc_section) - xc_section => section_vals_get_subs_vals(input, "DFT%XC") - ! No HFX allowed - hfx_section => section_vals_get_subs_vals(xc_section, "HF") - CALL section_vals_get(hfx_section, explicit=do_hfx) - IF (do_hfx) THEN - CALL cp_warn(__LOCATION__, "HFX not possible with AO based response solver. "// & - "Use the MO solver: RESPONSE_SOLVER/METOD MO_SOLVER") - CPABORT("hessian_op2@ec_orth_solver") + IF (dft_control%do_admm) THEN + xc_section => admm_env%xc_section_primary + ELSE + xc_section => section_vals_get_subs_vals(input, "DFT%XC") END IF ! add xc-kernel @@ -1208,8 +1218,26 @@ CONTAINS DO ispin = 1, nspins CALL pw_scale(v_xc(ispin)%pw, v_xc(ispin)%pw%pw_grid%dvol) + IF (ASSOCIATED(v_xc_tau)) THEN + CALL pw_scale(v_xc_tau(ispin)%pw, v_xc_tau(ispin)%pw%pw_grid%dvol) + END IF END DO + ! ADMM Correction + IF (dft_control%do_admm) THEN + IF (admm_env%aux_exch_func /= do_admm_aux_exch_func_none) THEN + IF (.NOT. ASSOCIATED(kpp1_env%deriv_set_admm)) THEN + xc_section_aux => admm_env%xc_section_aux + CALL get_admm_env(qs_env%admm_env, rho_aux_fit=rho_aux) + CALL qs_rho_get(rho_aux, rho_r=rho_r) + ALLOCATE (kpp1_env%deriv_set_admm, kpp1_env%rho_set_admm) + CALL xc_prep_2nd_deriv(kpp1_env%deriv_set_admm, kpp1_env%rho_set_admm, & + rho_r, auxbas_pw_pool, & + xc_section=xc_section_aux) + END IF + END IF + END IF + ! take trial density to build G^{H}[B] CALL pw_zero(rho_tot_gspace%pw) DO ispin = 1, nspins @@ -1231,27 +1259,18 @@ CONTAINS IF (ASSOCIATED(v_xc_tau)) CALL pw_scale(v_xc_tau(1)%pw, 2.0_dp) END IF - ! Init response kernel matrix - ! matrix G(B) - NULLIFY (matrix_G) - CALL dbcsr_allocate_matrix_set(matrix_G, nspins) DO ispin = 1, nspins - ALLOCATE (matrix_G(ispin)%matrix) - CALL dbcsr_copy(matrix_G(ispin)%matrix, matrix_s(1)%matrix, & - name="MATRIX Kernel") - ! Integrate with ground-state density matrix, in non-orthogonal basis CALL integrate_v_rspace(v_rspace=v_xc(ispin)%pw, & pmat=rho_ao(ispin), & - hmat=matrix_G(ispin), & + hmat=p_env%kpp1(ispin), & qs_env=qs_env, & calculate_forces=.FALSE., & basis_type="ORB") - IF (ASSOCIATED(v_xc_tau)) THEN CALL integrate_v_rspace(v_rspace=v_xc_tau(ispin)%pw, & pmat=rho_ao(ispin), & - hmat=matrix_G(ispin), & + hmat=p_env%kpp1(ispin), & qs_env=qs_env, & compute_tau=.TRUE., & calculate_forces=.FALSE., & @@ -1259,6 +1278,13 @@ CONTAINS END IF END DO + ! Hartree-Fock contribution + CALL apply_hfx(qs_env, p_env) + ! Calculate ADMM exchange correction to kernel + CALL apply_xc_admm(qs_env, p_env) + ! Add contribution from ADMM exchange correction to kernel + CALL p_env_finish_kpp1(qs_env, p_env) + ! Calculate KG correction to kernel IF (dft_control%qs_control%do_kg) THEN IF (qs_env%kg_env%tnadd_method == kg_tnadd_embed .OR. & @@ -1266,9 +1292,9 @@ CONTAINS CPASSERT(dft_control%nimages == 1) ekin_mol = 0.0_dp - CALL qs_rho_get(rho1, rho_ao=rho1_ao) + CALL qs_rho_get(p_env%rho1, rho_ao=rho1_ao) CALL kg_ekin_subset(qs_env=qs_env, & - ks_matrix=matrix_G, & + ks_matrix=p_env%kpp1, & ekin_mol=ekin_mol, & calc_force=.FALSE., & do_kernel=.TRUE., & @@ -1276,7 +1302,18 @@ CONTAINS END IF END IF + ! Init response kernel matrix + ! matrix G(B) + NULLIFY (matrix_G) + CALL dbcsr_allocate_matrix_set(matrix_G, nspins) + DO ispin = 1, nspins + ALLOCATE (matrix_G(ispin)%matrix) + CALL dbcsr_copy(matrix_G(ispin)%matrix, p_env%kpp1(ispin)%matrix, & + name="MATRIX Kernel") + END DO + ! Transforming G(B) into orthonormal basis + ! Careful, this de-symmetrizes matrix_G DO ispin = 1, nspins CALL transform_m_orth(matrix_G(ispin)%matrix, matrix_s_sqrt_inv, eps_filter) CALL dbcsr_filter(matrix_G(ispin)%matrix, eps_filter) diff --git a/src/energy_corrections.F b/src/energy_corrections.F index 9f3650f0b1..69c613cdb6 100644 --- a/src/energy_corrections.F +++ b/src/energy_corrections.F @@ -6,14 +6,15 @@ !--------------------------------------------------------------------------------------------------! ! ************************************************************************************************** -!> \brief Routines for a Harris type energy correction on top of a -!> Kohn-Sham calculation +!> \brief Routines for an energy correction on top of a Kohn-Sham calculation !> \par History !> 03.2014 created !> 09.2019 Moved from KG to Kohn-Sham +!> 08.2022 Add Density-Corrected DFT methods !> \author JGH ! ************************************************************************************************** MODULE energy_corrections + USE admm_methods, ONLY: admm_mo_merge_ks_matrix USE atomic_kind_types, ONLY: atomic_kind_type,& get_atomic_kind USE basis_set_types, ONLY: get_gto_basis_set,& @@ -81,11 +82,13 @@ MODULE energy_corrections gth_potential_type,& sgp_potential_type USE input_constants, ONLY: & - ec_diagonalization, ec_functional_harris, ec_matrix_sign, ec_matrix_tc2, ec_matrix_trs4, & - ec_ot_atomic, ec_ot_diag, ec_ot_gs, ot_precond_full_single_inverse, & + do_admm_basis_projection, do_admm_exch_scaling_none, do_admm_purify_none, & + ec_diagonalization, ec_functional_dc, ec_functional_harris, ec_matrix_sign, ec_matrix_tc2, & + ec_matrix_trs4, ec_ot_atomic, ec_ot_diag, ec_ot_gs, ot_precond_full_single_inverse, & ot_precond_solver_default, vdw_pairpot_dftd3, vdw_pairpot_dftd3bj, xc_vdw_fun_pairpot USE input_section_types, ONLY: section_get_ival,& section_get_lval,& + section_vals_get,& section_vals_get_subs_vals,& section_vals_type,& section_vals_val_get @@ -108,6 +111,7 @@ MODULE energy_corrections preconditioner_type USE pw_env_types, ONLY: pw_env_get,& pw_env_type + USE pw_grid_types, ONLY: pw_grid_type USE pw_methods, ONLY: pw_axpy,& pw_copy,& pw_integral_ab,& @@ -124,6 +128,7 @@ MODULE energy_corrections REALDATA3D,& REALSPACE,& RECIPROCALSPACE,& + pw_create,& pw_p_type,& pw_type USE qs_atomic_block, ONLY: calculate_atomic_block_dm @@ -148,6 +153,7 @@ MODULE energy_corrections USE qs_kinetic, ONLY: build_kinetic_matrix USE qs_ks_methods, ONLY: calc_rho_tot_gspace USE qs_ks_types, ONLY: qs_ks_env_type + USE qs_linres_kernel, ONLY: hfx_matrix USE qs_mo_methods, ONLY: calculate_subspace_eigenvalues,& make_basis_sm USE qs_mo_types, ONLY: deallocate_mo_set,& @@ -199,7 +205,9 @@ MODULE energy_corrections CONTAINS ! ************************************************************************************************** -!> \brief Energy correction to a KG simulation +!> \brief Energy Correction to a Kohn-Sham simulation +!> Available energy corrections: (1) Harris energy functional +!> (2) Density-corrected DFT !> !> \param qs_env ... !> \param ec_init ... @@ -233,7 +241,7 @@ CONTAINS NULLIFY (ec_env) CALL get_qs_env(qs_env, ec_env=ec_env) - ! Skip Harris functional calculation if ground-state is NOT converged + ! Skip energy correction if ground-state is NOT converged IF (.NOT. ec_env%do_skip) THEN ec_env%should_update = .TRUE. @@ -249,6 +257,7 @@ CONTAINS ec_env%exc = 0.0_dp ec_env%vhxc = 0.0_dp ec_env%edispersion = 0.0_dp + ec_env%exc_aux_fit = 0.0_dp END IF IF (my_calc_forces) THEN @@ -263,23 +272,16 @@ CONTAINS WRITE (unit_nr, '(T2,A,A,A,A,A)') "!", REPEAT("-", 29), & " Energy Correction ", REPEAT("-", 29), "!" END IF - END IF + ! Perform the energy correction + CALL energy_correction_low(qs_env, ec_env, my_calc_forces, unit_nr) + CALL get_qs_env(qs_env, energy=energy) - SELECT CASE (ec_env%energy_functional) - CASE (ec_functional_harris) - ! - CALL harris_energy(qs_env, ec_env, my_calc_forces, unit_nr) - ! - IF (ec_env%should_update) THEN - energy%nonscf_correction = ec_env%etotal - energy%total - CALL evaluate_ec_core_matrix_traces(qs_env, ec_env) - energy%total = ec_env%etotal - END IF - CASE DEFAULT - CPABORT("unknown energy correction") - END SELECT + IF (ec_env%should_update) THEN + energy%nonscf_correction = ec_env%etotal - energy%total + energy%total = ec_env%etotal + END IF IF (.NOT. my_calc_forces .AND. unit_nr > 0) THEN WRITE (unit_nr, '(T3,A,T56,F25.15)') "Energy Correction ", energy%nonscf_correction @@ -305,6 +307,165 @@ CONTAINS END SUBROUTINE energy_correction +! ************************************************************************************************** +!> \brief Energy Correction to a Kohn-Sham simulation +!> +!> \param qs_env ... +!> \param ec_env ... +!> \param calculate_forces ... +!> \param unit_nr ... +!> \par History +!> 03.2014 created +!> \author JGH +! ************************************************************************************************** + SUBROUTINE energy_correction_low(qs_env, ec_env, calculate_forces, unit_nr) + TYPE(qs_environment_type), POINTER :: qs_env + TYPE(energy_correction_type), POINTER :: ec_env + LOGICAL, INTENT(IN) :: calculate_forces + INTEGER, INTENT(IN) :: unit_nr + + INTEGER :: ispin, nspins + REAL(KIND=dp) :: exc + TYPE(dft_control_type), POINTER :: dft_control + TYPE(pw_env_type), POINTER :: pw_env + TYPE(pw_pool_type), POINTER :: auxbas_pw_pool + + IF (ec_env%should_update) THEN + CALL ec_build_neighborlist(qs_env, ec_env) + IF (.NOT. ASSOCIATED(ec_env%vh_rspace)) ALLOCATE (ec_env%vh_rspace) + IF (.NOT. ASSOCIATED(ec_env%vh_rspace%pw)) ALLOCATE (ec_env%vh_rspace%pw) + CALL ks_ref_potential(qs_env, & + ec_env%vh_rspace%pw, & + ec_env%vxc_rspace, & + ec_env%vtau_rspace, & + ec_env%vadmm_rspace, & + ec_env%ehartree, exc) + + SELECT CASE (ec_env%energy_functional) + CASE (ec_functional_harris) + + CALL ec_build_core_hamiltonian(qs_env, ec_env) + CALL ec_build_ks_matrix(qs_env, ec_env) + + IF (ec_env%mao) THEN + ! MAO basis + IF (ASSOCIATED(ec_env%mao_coef)) CALL dbcsr_deallocate_matrix_set(ec_env%mao_coef) + NULLIFY (ec_env%mao_coef) + CALL mao_generate_basis(qs_env, ec_env%mao_coef, ref_basis_set="HARRIS", molecular=.TRUE., & + max_iter=ec_env%mao_max_iter, eps_grad=ec_env%mao_eps_grad, & + eps1_mao=ec_env%mao_eps1, iolevel=ec_env%mao_iolevel, unit_nr=unit_nr) + END IF + + CALL ec_ks_solver(qs_env, ec_env) + + CALL evaluate_ec_core_matrix_traces(qs_env, ec_env) + + CASE (ec_functional_dc) + + ! Prepare Density-corrected DFT (DC-DFT) calculation + CALL ec_dc_energy(qs_env, ec_env, calculate_forces=.FALSE.) + + ! Rebuild KS matrix with DC-DFT XC functional evaluated in ground-state density. + ! KS matrix might contain unwanted contributions + ! Calculate Hartree and XC related energies here + CALL ec_build_ks_matrix(qs_env, ec_env) + + CASE DEFAULT + CPABORT("unknown energy correction") + END SELECT + + ! dispersion through pairpotentials + CALL ec_disp(qs_env, ec_env, calculate_forces=.FALSE.) + + ! Calculate total energy + CALL ec_energy(ec_env, unit_nr) + + END IF + + IF (calculate_forces) THEN + + CALL ec_disp(qs_env, ec_env, calculate_forces=.TRUE.) + + SELECT CASE (ec_env%energy_functional) + CASE (ec_functional_harris) + + CALL ec_build_core_hamiltonian_force(qs_env, ec_env, & + ec_env%matrix_p, & + ec_env%matrix_s, & + ec_env%matrix_w) + CALL ec_build_ks_matrix_force(qs_env, ec_env) + CASE (ec_functional_dc) + + ! Prepare Density-corrected DFT (DC-DFT) calculation + ! by getting ground-state matrices + CALL ec_dc_energy(qs_env, ec_env, calculate_forces=.TRUE.) + + CALL ec_build_core_hamiltonian_force(qs_env, ec_env, & + ec_env%matrix_p, & + ec_env%matrix_s, & + ec_env%matrix_w) + CALL ec_dc_build_ks_matrix_force(qs_env, ec_env) + + CASE DEFAULT + CPABORT("unknown energy correction") + END SELECT + + CALL response_calculation(qs_env, ec_env) + + ! Allocate response density on real space grid for use in properties + ! Calculated in response_force + CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, pw_env=pw_env) + nspins = dft_control%nspins + + CPASSERT(ASSOCIATED(pw_env)) + CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool) + ALLOCATE (ec_env%rhoz_r(nspins)) + DO ispin = 1, nspins + ALLOCATE (ec_env%rhoz_r(nspins)%pw) + CALL pw_pool_create_pw(auxbas_pw_pool, ec_env%rhoz_r(ispin)%pw, & + use_data=REALDATA3D, in_space=REALSPACE) + END DO + + CALL response_force(qs_env, & + vh_rspace=ec_env%vh_rspace%pw, & + vxc_rspace=ec_env%vxc_rspace, & + vtau_rspace=ec_env%vtau_rspace, & + vadmm_rspace=ec_env%vadmm_rspace, & + matrix_hz=ec_env%matrix_hz, & + matrix_pz=ec_env%matrix_z, & + matrix_pz_admm=ec_env%z_admm, & + matrix_wz=ec_env%matrix_wz, & + rhopz_r=ec_env%rhoz_r, & + zehartree=ec_env%ehartree, & + zexc=ec_env%exc, & + zexc_aux_fit=ec_env%exc_aux_fit, & + p_env=ec_env%p_env) + + CALL ec_properties(qs_env, ec_env) + + ! Deallocate Harris density and response density on grid + DO ispin = 1, nspins + CALL pw_pool_give_back_pw(auxbas_pw_pool, ec_env%rhoout_r(ispin)%pw) + CALL pw_pool_give_back_pw(auxbas_pw_pool, ec_env%rhoz_r(ispin)%pw) + DEALLOCATE (ec_env%rhoout_r(ispin)%pw, ec_env%rhoz_r(ispin)%pw) + END DO + DEALLOCATE (ec_env%rhoout_r, ec_env%rhoz_r) + + ! Deallocate matrices + IF (ASSOCIATED(ec_env%matrix_ks)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_ks) + IF (ASSOCIATED(ec_env%matrix_h)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_h) + IF (ASSOCIATED(ec_env%matrix_s)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_s) + IF (ASSOCIATED(ec_env%matrix_t)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_t) + IF (ASSOCIATED(ec_env%matrix_p)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_p) + IF (ASSOCIATED(ec_env%matrix_w)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_w) + IF (ASSOCIATED(ec_env%matrix_hz)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_hz) + IF (ASSOCIATED(ec_env%matrix_wz)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_wz) + IF (ASSOCIATED(ec_env%matrix_z)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_z) + + END IF + + END SUBROUTINE energy_correction_low + ! ************************************************************************************************** !> \brief Calculates the traces of the core matrices and the density matrix. !> \param qs_env ... @@ -338,112 +499,621 @@ CONTAINS END SUBROUTINE evaluate_ec_core_matrix_traces ! ************************************************************************************************** -!> \brief Harris Energy Correction to a Kohn-Sham simulation +!> \brief Prepare DC-DFT calculation by copying unaffected ground-state matrices (core Hamiltonian, +!> density matrix) into energy correction environment and rebuild the overlap matrix !> !> \param qs_env ... !> \param ec_env ... !> \param calculate_forces ... -!> \param unit_nr ... !> \par History -!> 03.2014 created -!> \author JGH +!> 07.2022 created +!> \author fbelle ! ************************************************************************************************** - SUBROUTINE harris_energy(qs_env, ec_env, calculate_forces, unit_nr) + SUBROUTINE ec_dc_energy(qs_env, ec_env, calculate_forces) TYPE(qs_environment_type), POINTER :: qs_env TYPE(energy_correction_type), POINTER :: ec_env LOGICAL, INTENT(IN) :: calculate_forces - INTEGER, INTENT(IN) :: unit_nr - INTEGER :: ispin, nspins - REAL(KIND=dp) :: exc + CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_dc_energy' + + CHARACTER(LEN=default_string_length) :: headline + INTEGER :: handle, ispin, nspins + TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p, matrix_s, matrix_w TYPE(dft_control_type), POINTER :: dft_control - TYPE(pw_env_type), POINTER :: pw_env - TYPE(pw_pool_type), POINTER :: auxbas_pw_pool + TYPE(qs_ks_env_type), POINTER :: ks_env + TYPE(qs_rho_type), POINTER :: rho - IF (ec_env%should_update) THEN - CALL ec_build_neighborlist(qs_env, ec_env) - ! - CALL ec_build_core_hamiltonian(qs_env, ec_env) - IF (.NOT. ASSOCIATED(ec_env%vh_rspace)) ALLOCATE (ec_env%vh_rspace) - IF (.NOT. ASSOCIATED(ec_env%vh_rspace%pw)) ALLOCATE (ec_env%vh_rspace%pw) - CALL ks_ref_potential(qs_env, ec_env%vh_rspace%pw, ec_env%vxc_rspace, ec_env%vtau_rspace, ec_env%vadmm_rspace, & - ec_env%ehartree, exc) - CALL ec_build_ks_matrix(qs_env, ec_env) - ! - IF (ec_env%mao) THEN - ! MAO basis - IF (ASSOCIATED(ec_env%mao_coef)) CALL dbcsr_deallocate_matrix_set(ec_env%mao_coef) - NULLIFY (ec_env%mao_coef) - CALL mao_generate_basis(qs_env, ec_env%mao_coef, ref_basis_set="HARRIS", molecular=.TRUE., & - max_iter=ec_env%mao_max_iter, eps_grad=ec_env%mao_eps_grad, & - eps1_mao=ec_env%mao_eps1, iolevel=ec_env%mao_iolevel, unit_nr=unit_nr) - END IF - ! - CALL ec_ks_solver(qs_env, ec_env) - ! - ! dispersion through pairpotentials - CALL ec_disp(qs_env, ec_env, calculate_forces=.FALSE.) - CALL ec_energy(ec_env, unit_nr) + CALL timeset(routineN, handle) + + NULLIFY (dft_control, ks_env, matrix_h, matrix_p, matrix_s, matrix_w, rho) + CALL get_qs_env(qs_env=qs_env, & + dft_control=dft_control, & + ks_env=ks_env, & + matrix_h_kp=matrix_h, & + matrix_s_kp=matrix_s, & + matrix_w_kp=matrix_w, & + rho=rho) + CALL qs_rho_get(rho, rho_ao_kp=matrix_p) + nspins = dft_control%nspins + + ! For density-corrected DFT only the ground-state matrices are required + ! Comply with ec_env environment for property calculations later + CALL build_overlap_matrix(ks_env, matrixkp_s=ec_env%matrix_s, & + matrix_name="OVERLAP MATRIX", & + basis_type_a="HARRIS", & + basis_type_b="HARRIS", & + sab_nl=ec_env%sab_orb) + + ! Core Hamiltonian matrix + IF (ASSOCIATED(ec_env%matrix_h)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_h) + CALL dbcsr_allocate_matrix_set(ec_env%matrix_h, 1, 1) + headline = "CORE HAMILTONIAN MATRIX" + ALLOCATE (ec_env%matrix_h(1, 1)%matrix) + CALL dbcsr_create(ec_env%matrix_h(1, 1)%matrix, name=TRIM(headline), & + template=matrix_s(1, 1)%matrix, matrix_type=dbcsr_type_symmetric) + CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_h(1, 1)%matrix, ec_env%sab_orb) + CALL dbcsr_copy(ec_env%matrix_h(1, 1)%matrix, matrix_h(1, 1)%matrix) + + ! Density matrix + IF (ASSOCIATED(ec_env%matrix_p)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_p) + CALL dbcsr_allocate_matrix_set(ec_env%matrix_p, nspins, 1) + headline = "DENSITY MATRIX" + DO ispin = 1, nspins + ALLOCATE (ec_env%matrix_p(ispin, 1)%matrix) + CALL dbcsr_create(ec_env%matrix_p(ispin, 1)%matrix, name=TRIM(headline), & + template=matrix_s(1, 1)%matrix, matrix_type=dbcsr_type_symmetric) + CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_p(ispin, 1)%matrix, ec_env%sab_orb) + CALL dbcsr_copy(ec_env%matrix_p(ispin, 1)%matrix, matrix_p(ispin, 1)%matrix) + END DO - END IF IF (calculate_forces) THEN - CALL ec_disp(qs_env, ec_env, calculate_forces=.TRUE.) - CALL ec_build_core_hamiltonian_force(qs_env, ec_env) - CALL ec_build_ks_matrix_force(qs_env, ec_env) - CALL response_calculation(qs_env, ec_env) - - ! Allocate response density on real space grid for use in properties - ! Calculated in response_force - CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, pw_env=pw_env) - nspins = dft_control%nspins - CPASSERT(ASSOCIATED(pw_env)) - CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool) - ALLOCATE (ec_env%rhoz_r(nspins)) - DO ispin = 1, nspins - ALLOCATE (ec_env%rhoz_r(nspins)%pw) - CALL pw_pool_create_pw(auxbas_pw_pool, ec_env%rhoz_r(ispin)%pw, & - use_data=REALDATA3D, in_space=REALSPACE) - END DO - - CALL response_force(qs_env, & - vh_rspace=ec_env%vh_rspace%pw, & - vxc_rspace=ec_env%vxc_rspace, & - vtau_rspace=ec_env%vtau_rspace, & - vadmm_rspace=ec_env%vadmm_rspace, & - matrix_hz=ec_env%matrix_hz, & - matrix_pz=ec_env%matrix_z, & - matrix_pz_admm=ec_env%z_admm, & - matrix_wz=ec_env%matrix_wz, & - rhopz_r=ec_env%rhoz_r, & - zehartree=ec_env%ehartree, & - p_env=ec_env%p_env, & - zexc=ec_env%exc) - - CALL ec_properties(qs_env, ec_env) - - ! Deallocate Harris density and response density on grid - DO ispin = 1, nspins - CALL pw_pool_give_back_pw(auxbas_pw_pool, ec_env%rhoout_r(ispin)%pw) - CALL pw_pool_give_back_pw(auxbas_pw_pool, ec_env%rhoz_r(ispin)%pw) - DEALLOCATE (ec_env%rhoout_r(ispin)%pw, ec_env%rhoz_r(ispin)%pw) - END DO - DEALLOCATE (ec_env%rhoout_r, ec_env%rhoz_r) - - ! Deallocate matrices - IF (ASSOCIATED(ec_env%matrix_ks)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_ks) - IF (ASSOCIATED(ec_env%matrix_h)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_h) - IF (ASSOCIATED(ec_env%matrix_s)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_s) - IF (ASSOCIATED(ec_env%matrix_t)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_t) - IF (ASSOCIATED(ec_env%matrix_p)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_p) + ! Energy-weighted density matrix IF (ASSOCIATED(ec_env%matrix_w)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_w) - IF (ASSOCIATED(ec_env%matrix_hz)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_hz) - IF (ASSOCIATED(ec_env%matrix_wz)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_wz) - IF (ASSOCIATED(ec_env%matrix_z)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_z) + CALL dbcsr_allocate_matrix_set(ec_env%matrix_w, nspins, 1) + headline = "ENERGY-WEIGHTED DENSITY MATRIX" + DO ispin = 1, nspins + ALLOCATE (ec_env%matrix_w(ispin, 1)%matrix) + CALL dbcsr_create(ec_env%matrix_w(ispin, 1)%matrix, name=TRIM(headline), & + template=matrix_s(1, 1)%matrix, matrix_type=dbcsr_type_symmetric) + CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_w(ispin, 1)%matrix, ec_env%sab_orb) + CALL dbcsr_copy(ec_env%matrix_w(ispin, 1)%matrix, matrix_w(ispin, 1)%matrix) + END DO END IF - END SUBROUTINE harris_energy + ! External field (nonperiodic case) + ec_env%efield_nuclear = 0.0_dp + ec_env%efield_elec = 0.0_dp + CALL ec_efield_local_operator(qs_env, ec_env, calculate_forces=.FALSE.) + + CALL timestop(handle) + + END SUBROUTINE ec_dc_energy + +! ************************************************************************************************** +!> \brief Kohn-Sham matrix contributions to force in DC-DFT +!> also calculate right-hand-side matrix B for response equations AX=B +!> \param qs_env ... +!> \param ec_env ... +!> \par History +!> 08.2022 adapted from qs_ks_build_kohn_sham_matrix +!> \author fbelle +! ************************************************************************************************** + SUBROUTINE ec_dc_build_ks_matrix_force(qs_env, ec_env) + TYPE(qs_environment_type), POINTER :: qs_env + TYPE(energy_correction_type), POINTER :: ec_env + + CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_dc_build_ks_matrix_force' + + INTEGER :: handle, i, iounit, ispin, natom, nspins + LOGICAL :: do_adiabatic_rescaling, do_hfx, & + use_virial + REAL(dp) :: ehartree, eovrl, exc, fconv + REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: ftot + REAL(dp), DIMENSION(3) :: fodeb + REAL(KIND=dp), DIMENSION(3, 3) :: h_stress, pv_loc, stdeb, sttot + TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set + TYPE(cell_type), POINTER :: cell + TYPE(cp_logger_type), POINTER :: logger + TYPE(cp_para_env_type), POINTER :: para_env + TYPE(dbcsr_p_type) :: scrm + TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: dbcsr_work, matrix_ks, matrix_p, matrix_s + TYPE(dft_control_type), POINTER :: dft_control + TYPE(neighbor_list_set_p_type), DIMENSION(:), & + POINTER :: sab_orb + TYPE(pw_env_type), POINTER :: pw_env + TYPE(pw_grid_type), POINTER :: pw_grid + TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_r, v_rspace, v_rspace_in, & + v_tau_rspace + TYPE(pw_poisson_type), POINTER :: poisson_env + TYPE(pw_pool_type), POINTER :: auxbas_pw_pool + TYPE(pw_type) :: rho_tot_gspace, v_hartree_gspace, & + v_hartree_rspace + TYPE(qs_force_type), DIMENSION(:), POINTER :: force + TYPE(qs_ks_env_type), POINTER :: ks_env + TYPE(qs_rho_type), POINTER :: rho + TYPE(section_vals_type), POINTER :: adiabatic_rescaling_section, & + hfx_sections, input + TYPE(virial_type), POINTER :: virdeb, virial + + CALL timeset(routineN, handle) + + logger => cp_get_default_logger() + IF (logger%para_env%ionode) THEN + iounit = cp_logger_get_default_unit_nr(logger, local=.TRUE.) + ELSE + iounit = -1 + END IF + + NULLIFY (atomic_kind_set, cell, dft_control, force, ks_env, matrix_ks, & + matrix_p, matrix_s, para_env, pw_env, rho, sab_orb, virial) + CALL get_qs_env(qs_env=qs_env, & + cell=cell, & + dft_control=dft_control, & + force=force, & + input=input, & + ks_env=ks_env, & + matrix_ks=matrix_ks, & + matrix_s=matrix_s, & + para_env=para_env, & + pw_env=pw_env, & + rho=rho, & + sab_orb=sab_orb, & + virial=virial) + CPASSERT(ASSOCIATED(pw_env)) + + nspins = dft_control%nspins + use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer) + + fconv = 1.0E-9_dp*pascal/cell%deth + IF (debug_stress .AND. use_virial) THEN + sttot = virial%pv_virial + END IF + + NULLIFY (auxbas_pw_pool, poisson_env) + ! gets the tmp grids + CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, & + poisson_env=poisson_env) + + ! Calculate the Hartree potential + CALL pw_pool_create_pw(auxbas_pw_pool, v_hartree_gspace, & + use_data=COMPLEXDATA1D, in_space=RECIPROCALSPACE) + CALL pw_pool_create_pw(auxbas_pw_pool, rho_tot_gspace, & + use_data=COMPLEXDATA1D, in_space=RECIPROCALSPACE) + CALL pw_pool_create_pw(auxbas_pw_pool, v_hartree_rspace, & + use_data=REALDATA3D, in_space=REALSPACE) + + ! Get the total input density in g-space [ions + electrons] + CALL calc_rho_tot_gspace(rho_tot_gspace, qs_env, rho) + + ! v_H[n_in] + IF (use_virial) THEN + + ! Stress tensor - Volume and Green function contribution + h_stress(:, :) = 0.0_dp + CALL pw_poisson_solve(poisson_env, & + density=rho_tot_gspace, & + ehartree=ehartree, & + vhartree=v_hartree_gspace, & + h_stress=h_stress) + + virial%pv_ehartree = virial%pv_ehartree + h_stress/REAL(para_env%num_pe, dp) + virial%pv_virial = virial%pv_virial + h_stress/REAL(para_env%num_pe, dp) + + IF (debug_stress) THEN + stdeb = fconv*(h_stress/REAL(para_env%num_pe, dp)) + CALL mp_sum(stdeb, para_env%group) + IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") & + 'STRESS| GREEN 1st V_H[n_in]*n_in ', one_third_sum_diag(stdeb), det_3x3(stdeb) + END IF + + ELSE + CALL pw_poisson_solve(poisson_env, rho_tot_gspace, ehartree, & + v_hartree_gspace) + END IF + + CALL pw_transfer(v_hartree_gspace, v_hartree_rspace) + CALL pw_scale(v_hartree_rspace, v_hartree_rspace%pw_grid%dvol) + + ! Save Harris on real space grid for use in properties + CALL qs_rho_get(rho, rho_r=rho_r) + ALLOCATE (ec_env%rhoout_r(nspins)) + DO ispin = 1, nspins + ALLOCATE (ec_env%rhoout_r(ispin)%pw) + CALL pw_pool_create_pw(auxbas_pw_pool, ec_env%rhoout_r(ispin)%pw, & + use_data=REALDATA3D, in_space=REALSPACE) + CALL pw_copy(rho_r(ispin)%pw, ec_env%rhoout_r(ispin)%pw) + END DO + + ! Getting nuclear force contribution from the core charge density + ! Vh(rho_c + rho_in) + IF (debug_forces) fodeb(1:3) = force(1)%rho_core(1:3, 1) + IF (debug_stress .AND. use_virial) stdeb = virial%pv_ehartree + CALL integrate_v_core_rspace(v_hartree_rspace, qs_env) + IF (debug_forces) THEN + fodeb(1:3) = force(1)%rho_core(1:3, 1) - fodeb(1:3) + CALL mp_sum(fodeb, para_env%group) + IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Vtot*dncore", fodeb + END IF + IF (debug_stress .AND. use_virial) THEN + stdeb = fconv*(virial%pv_ehartree - stdeb) + CALL mp_sum(stdeb, para_env%group) + IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") & + 'STRESS| Vtot*dncore', one_third_sum_diag(stdeb), det_3x3(stdeb) + END IF + + ! v_XC[n_in]_DC + ! v_rspace and v_tau_rspace are generated from the auxbas pool + NULLIFY (v_rspace, v_tau_rspace) + + ! only activate stress calculation if + IF (use_virial) virial%pv_calculate = .TRUE. + + CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho, xc_section=ec_env%xc_section, & + vxc_rho=v_rspace, vxc_tau=v_tau_rspace, exc=exc, just_energy=.FALSE.) + + IF (.NOT. ASSOCIATED(v_rspace)) THEN + ALLOCATE (v_rspace(nspins)) + DO ispin = 1, nspins + ALLOCATE (v_rspace(ispin)%pw) + CALL pw_pool_create_pw(auxbas_pw_pool, v_rspace(ispin)%pw, & + use_data=REALDATA3D, in_space=REALSPACE) + CALL pw_zero(v_rspace(ispin)%pw) + END DO + END IF + + IF (use_virial) THEN + virial%pv_exc = virial%pv_exc - virial%pv_xc + virial%pv_virial = virial%pv_virial - virial%pv_xc + ! virial%pv_xc will be zeroed in the xc routines + END IF + + pw_grid => v_hartree_rspace%pw_grid + ALLOCATE (v_rspace_in(nspins)) + DO ispin = 1, nspins + ALLOCATE (v_rspace_in(ispin)%pw) + CALL pw_create(v_rspace_in(ispin)%pw, pw_grid, & + use_data=REALDATA3D, in_space=REALSPACE) + END DO + + ! v_rspace_in = v_H[n_in] + v_xc[n_in] calculated in ks_ref_potential + DO ispin = 1, nspins + ! v_xc[n_in]_GS + CALL pw_transfer(ec_env%vxc_rspace(ispin)%pw, v_rspace_in(ispin)%pw) + ! add v_H[n_in] + CALL pw_axpy(ec_env%vh_rspace%pw, v_rspace_in(ispin)%pw) + END DO + + ! initialize src matrix + NULLIFY (scrm%matrix) + ALLOCATE (scrm%matrix) + CALL dbcsr_create(scrm%matrix, template=matrix_s(1)%matrix) + CALL cp_dbcsr_alloc_block_from_nbl(scrm%matrix, ec_env%sab_orb) + + ! Stress-tensor contribution derivative of integrand + ! int v_Hxc[n^în]*n^out + IF (use_virial) THEN + pv_loc = virial%pv_virial + END IF + + CALL qs_rho_get(rho, rho_ao=matrix_p) + IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1) + IF (debug_stress .AND. use_virial) stdeb = virial%pv_virial + + DO ispin = 1, nspins + ! Add v_H[n_in] + v_xc[n_in] = v_rspace + CALL pw_scale(v_rspace(ispin)%pw, v_rspace(ispin)%pw%pw_grid%dvol) + CALL pw_axpy(v_hartree_rspace, v_rspace(ispin)%pw) + ! integrate over potential + CALL integrate_v_rspace(v_rspace=v_rspace(ispin)%pw, & + hmat=scrm, & + pmat=matrix_p(ispin), & + qs_env=qs_env, & + calculate_forces=.TRUE., & + basis_type="HARRIS", & + task_list_external=ec_env%task_list) + END DO + + IF (debug_forces) THEN + fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3) + CALL mp_sum(fodeb, para_env%group) + IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pout*dVhxc ", fodeb + END IF + IF (debug_stress .AND. use_virial) THEN + stdeb = fconv*(virial%pv_virial - stdeb) + CALL mp_sum(stdeb, para_env%group) + IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") & + 'STRESS| INT Pout*dVhxc ', one_third_sum_diag(stdeb), det_3x3(stdeb) + END IF + + IF (ASSOCIATED(v_tau_rspace)) THEN + IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1) + IF (debug_stress .AND. use_virial) stdeb = virial%pv_virial + DO ispin = 1, nspins + CALL pw_scale(v_tau_rspace(ispin)%pw, v_tau_rspace(ispin)%pw%pw_grid%dvol) + ! integrate over Tau-potential + CALL integrate_v_rspace(v_rspace=v_tau_rspace(ispin)%pw, & + hmat=scrm, & + pmat=matrix_p(ispin), & + qs_env=qs_env, & + calculate_forces=.TRUE., & + compute_tau=.TRUE., & + basis_type="HARRIS", & + task_list_external=ec_env%task_list) + END DO + + IF (debug_forces) THEN + fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3) + CALL mp_sum(fodeb, para_env%group) + IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pout*dVhxc_tau ", fodeb + END IF + IF (debug_stress .AND. use_virial) THEN + stdeb = fconv*(virial%pv_virial - stdeb) + CALL mp_sum(stdeb, para_env%group) + IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") & + 'STRESS| INT Pout*dVhxc_tau ', one_third_sum_diag(stdeb), det_3x3(stdeb) + END IF + END IF + + ! Stress-tensor + IF (use_virial) THEN + virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc) + END IF + + ! delete scr matrix + CALL dbcsr_release(scrm%matrix) + DEALLOCATE (scrm%matrix) + + !---------------------------------------------------- + ! Right-hand-side matrix B for linear response equations AX = B + !---------------------------------------------------- + + NULLIFY (ec_env%matrix_hz) + CALL dbcsr_allocate_matrix_set(ec_env%matrix_hz, nspins) + DO ispin = 1, nspins + ALLOCATE (ec_env%matrix_hz(ispin)%matrix) + CALL dbcsr_create(ec_env%matrix_hz(ispin)%matrix, template=matrix_s(1)%matrix) + CALL dbcsr_copy(ec_env%matrix_hz(ispin)%matrix, matrix_s(1)%matrix) + CALL dbcsr_set(ec_env%matrix_hz(ispin)%matrix, 0.0_dp) + END DO + + DO ispin = 1, nspins + ! v_rspace = v_Hxc[n_in]_DC - v_Hxc[n_in]_GS + CALL pw_axpy(v_rspace_in(ispin)%pw, v_rspace(ispin)%pw, -1.0_dp) + END DO + + DO ispin = 1, nspins + CALL integrate_v_rspace(v_rspace=v_rspace(ispin)%pw, & + hmat=ec_env%matrix_hz(ispin), & + pmat=matrix_p(ispin), & + qs_env=qs_env, & + calculate_forces=.FALSE., & + basis_type="HARRIS", & + task_list_external=ec_env%task_list) + END DO + + ! Check if mGGA functionals are used + IF (dft_control%use_kinetic_energy_density) THEN + + ! If DC-DFT without mGGA functional, this needs to be allocated now. + IF (.NOT. ASSOCIATED(v_tau_rspace)) THEN + ALLOCATE (v_tau_rspace(nspins)) + DO ispin = 1, nspins + ALLOCATE (v_tau_rspace(ispin)%pw) + CALL pw_pool_create_pw(auxbas_pw_pool, v_tau_rspace(ispin)%pw, & + use_data=REALDATA3D, in_space=REALSPACE) + CALL pw_zero(v_tau_rspace(ispin)%pw) + END DO + END IF + + DO ispin = 1, nspins + ! v_tau_rspace = v_Hxc_tau[n_in]_DC - v_Hxc_tau[n_in]_GS + IF (ASSOCIATED(ec_env%vtau_rspace)) THEN + CALL pw_axpy(ec_env%vtau_rspace(ispin)%pw, v_tau_rspace(ispin)%pw, -1.0_dp) + END IF + ! integrate over Tau-potential + CALL integrate_v_rspace(v_rspace=v_tau_rspace(ispin)%pw, & + hmat=ec_env%matrix_hz(ispin), & + pmat=matrix_p(ispin), & + qs_env=qs_env, & + calculate_forces=.FALSE., compute_tau=.TRUE., & + basis_type="HARRIS", & + task_list_external=ec_env%task_list) + END DO + END IF + + ! Need to also subtract HFX contribution from ec_env%matrix_hz + hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%HF") + CALL section_vals_get(hfx_sections, explicit=do_hfx) + adiabatic_rescaling_section => section_vals_get_subs_vals(input, "DFT%XC%ADIABATIC_RESCALING") + CALL section_vals_get(adiabatic_rescaling_section, explicit=do_adiabatic_rescaling) + + IF (do_hfx .AND. .NOT. do_adiabatic_rescaling) THEN + + ! HFX-ADMM + IF (dft_control%do_admm) THEN + + IF (dft_control%admm_control%purification_method /= do_admm_purify_none) THEN + CPABORT("ADMM: Linear Response needs purification_method=none") + END IF + IF (dft_control%admm_control%scaling_model /= do_admm_exch_scaling_none) THEN + CPABORT("ADMM: Linear Response needs scaling_model=none") + END IF + IF (dft_control%admm_control%method /= do_admm_basis_projection) THEN + CPABORT("ADMM: Linear Response needs admm_method=basis_projection") + END IF + + ! HFX matrix in ADMM environment already calculated during ground-state + ! Here, need to subtract it from RHS matrix for the response equations + + ! make a work matrix to isolate the HFX matrix + ! dbcsr_work = old KS matrix + ALLOCATE (dbcsr_work(nspins)) + DO ispin = 1, nspins + ALLOCATE (dbcsr_work(ispin)%matrix) + CALL dbcsr_copy(dbcsr_work(ispin)%matrix, matrix_ks(ispin)%matrix) + END DO + + ! Can actually only call subroutine merge_ks_matrix_none + ! Adds to matrix_ks of qs_env + CALL admm_mo_merge_ks_matrix(qs_env) + + ! dbcsr_work = HFX matrix = (KS+HFX) - KS + DO ispin = 1, nspins + CALL dbcsr_add(dbcsr_work(ispin)%matrix, matrix_ks(ispin)%matrix, -1.0_dp, 1.0_dp) + END DO + + ! Subtract from RHS matrix + DO ispin = 1, nspins + CALL dbcsr_add(ec_env%matrix_hz(ispin)%matrix, dbcsr_work(ispin)%matrix, & + 1.0_dp, -1.0_dp) + END DO + + ! Restore Kohn-Sham matrix, to make it consistent for response calculation + DO ispin = 1, nspins + CALL dbcsr_copy(matrix_ks(ispin)%matrix, ec_env%matrix_ks(ispin, 1)%matrix) + END DO + + CALL dbcsr_deallocate_matrix_set(dbcsr_work) + + ! conventional HFX + ELSE + ALLOCATE (dbcsr_work(nspins)) + DO ispin = 1, nspins + ALLOCATE (dbcsr_work(ispin)%matrix) + CALL dbcsr_copy(dbcsr_work(ispin)%matrix, matrix_s(ispin)%matrix) + CALL dbcsr_set(dbcsr_work(ispin)%matrix, 0.0_dp) + END DO + + ! Get the HFX contribution to the Kohn-Sham matrix + CALL hfx_matrix(dbcsr_work, matrix_p, qs_env, hfx_sections) + + ! Subtract from RHS matrix + DO ispin = 1, nspins + CALL dbcsr_add(ec_env%matrix_hz(ispin)%matrix, dbcsr_work(ispin)%matrix, & + 1.0_dp, -1.0_dp) + END DO + CALL dbcsr_deallocate_matrix_set(dbcsr_work) + + END IF ! do_admm + + END IF ! do_fhx + + ! Core overlap + IF (debug_forces) fodeb(1:3) = force(1)%core_overlap(1:3, 1) + IF (debug_stress .AND. use_virial) stdeb = virial%pv_ecore_overlap + CALL calculate_ecore_overlap(qs_env, para_env, .TRUE., E_overlap_core=eovrl) + IF (debug_forces) THEN + fodeb(1:3) = force(1)%core_overlap(1:3, 1) - fodeb(1:3) + CALL mp_sum(fodeb, para_env%group) + IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: CoreOverlap", fodeb + END IF + IF (debug_stress .AND. use_virial) THEN + stdeb = fconv*(stdeb - virial%pv_ecore_overlap) + CALL mp_sum(stdeb, para_env%group) + IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") & + 'STRESS| CoreOverlap ', one_third_sum_diag(stdeb), det_3x3(stdeb) + END IF + + IF (debug_forces) THEN + CALL get_qs_env(qs_env, natom=natom, atomic_kind_set=atomic_kind_set) + ALLOCATE (ftot(3, natom)) + CALL total_qs_force(ftot, force, atomic_kind_set) + fodeb(1:3) = ftot(1:3, 1) + DEALLOCATE (ftot) + CALL mp_sum(fodeb, para_env%group) + IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Force Explicit", fodeb + END IF + + ! return pw grids + DO ispin = 1, nspins + CALL pw_pool_give_back_pw(auxbas_pw_pool, v_rspace(ispin)%pw) + CALL pw_pool_give_back_pw(auxbas_pw_pool, v_rspace_in(ispin)%pw) + DEALLOCATE (v_rspace(ispin)%pw, v_rspace_in(ispin)%pw) + IF (ASSOCIATED(v_tau_rspace)) THEN + CALL pw_pool_give_back_pw(auxbas_pw_pool, v_tau_rspace(ispin)%pw) + DEALLOCATE (v_tau_rspace(ispin)%pw) + END IF + END DO + + DEALLOCATE (v_rspace, v_rspace_in) + IF (ASSOCIATED(v_tau_rspace)) DEALLOCATE (v_tau_rspace) + ! + CALL pw_pool_give_back_pw(auxbas_pw_pool, rho_tot_gspace) + CALL pw_pool_give_back_pw(auxbas_pw_pool, v_hartree_gspace) + CALL pw_pool_give_back_pw(auxbas_pw_pool, v_hartree_rspace) + + ! Stress tensor - volume terms need to be stored, + ! for a sign correction in QS at the end of qs_force + IF (use_virial) THEN + IF (qs_env%energy_correction) THEN + ec_env%ehartree = ehartree + ec_env%exc = exc + END IF + END IF + + IF (debug_stress .AND. use_virial) THEN + ! In total: -1.0*E_H + stdeb = -1.0_dp*fconv*ehartree + IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") & + 'STRESS| VOL 1st v_H[n_in]*n_in', one_third_sum_diag(stdeb), det_3x3(stdeb) + + stdeb = -1.0_dp*fconv*exc + IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") & + 'STRESS| VOL 1st E_XC_DC[n_in]', one_third_sum_diag(stdeb), det_3x3(stdeb) + + ! For debugging, create a second virial environment, + ! apply volume terms immediately + CALL virial_create(virdeb) + CALL cp_virial(virial, virdeb) + + CALL mp_sum(virdeb%pv_overlap, para_env%group) + CALL mp_sum(virdeb%pv_ekinetic, para_env%group) + CALL mp_sum(virdeb%pv_ppl, para_env%group) + CALL mp_sum(virdeb%pv_ppnl, para_env%group) + CALL mp_sum(virdeb%pv_ecore_overlap, para_env%group) + CALL mp_sum(virdeb%pv_ehartree, para_env%group) + CALL mp_sum(virdeb%pv_exc, para_env%group) + CALL mp_sum(virdeb%pv_exx, para_env%group) + CALL mp_sum(virdeb%pv_vdw, para_env%group) + CALL mp_sum(virdeb%pv_mp2, para_env%group) + CALL mp_sum(virdeb%pv_nlcc, para_env%group) + CALL mp_sum(virdeb%pv_gapw, para_env%group) + CALL mp_sum(virdeb%pv_lrigpw, para_env%group) + CALL mp_sum(virdeb%pv_virial, para_env%group) + CALL symmetrize_virial(virdeb) + + ! apply stress-tensor 1st terms + DO i = 1, 3 + virdeb%pv_ehartree(i, i) = virdeb%pv_ehartree(i, i) - 2.0_dp*ehartree + virdeb%pv_virial(i, i) = virdeb%pv_virial(i, i) - exc & + - 2.0_dp*ehartree + virdeb%pv_exc(i, i) = virdeb%pv_exc(i, i) - exc + ! The factor 2 is a hack. It compensates the plus sign in h_stress/pw_poisson_solve. + ! The sign in pw_poisson_solve is correct for FIST, but not for QS. + ! There should be a more elegant solution to that ... + END DO + + CALL mp_sum(sttot, para_env%group) + stdeb = fconv*(virdeb%pv_virial - sttot) + IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") & + 'STRESS| Explicit electronic stress ', one_third_sum_diag(stdeb), det_3x3(stdeb) + + stdeb = fconv*(virdeb%pv_virial) + IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") & + 'STRESS| Explicit total stress ', one_third_sum_diag(stdeb), det_3x3(stdeb) + + CALL write_stress_tensor_components(virdeb, iounit, cell) + CALL write_stress_tensor(virdeb%pv_virial, iounit, cell, .FALSE.) + + CALL virial_release(virdeb) + + END IF + + CALL timestop(handle) + + END SUBROUTINE ec_dc_build_ks_matrix_force ! ************************************************************************************************** !> \brief ... @@ -476,13 +1146,11 @@ CONTAINS CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_build_core_hamiltonian' - INTEGER :: handle, iounit, nder, nimages + INTEGER :: handle, nder, nimages INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index LOGICAL :: calculate_forces, use_virial REAL(KIND=dp) :: eps_ppnl TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set - TYPE(cp_logger_type), POINTER :: logger - TYPE(cp_para_env_type), POINTER :: para_env TYPE(dft_control_type), POINTER :: dft_control TYPE(neighbor_list_set_p_type), DIMENSION(:), & POINTER :: sab_orb, sac_ppl, sap_ppnl @@ -494,18 +1162,16 @@ CONTAINS CALL timeset(routineN, handle) - logger => cp_get_default_logger() - IF (logger%para_env%ionode) THEN - iounit = cp_logger_get_default_unit_nr(logger, local=.TRUE.) - ELSE - iounit = -1 - END IF + NULLIFY (atomic_kind_set, cell_to_index, dft_control, ks_env, particle_set, qs_kind_set, virial) + + CALL get_qs_env(qs_env=qs_env, & + atomic_kind_set=atomic_kind_set, & + dft_control=dft_control, & + particle_set=particle_set, & + qs_kind_set=qs_kind_set, & + ks_env=ks_env) ! no k-points possible - CALL get_qs_env(qs_env=qs_env, & - dft_control=dft_control, & - ks_env=ks_env, & - para_env=para_env) nimages = dft_control%nimages IF (nimages /= 1) THEN CPABORT("K-points for Harris functional not implemented") @@ -516,6 +1182,10 @@ CONTAINS CPABORT("Harris functional for GAPW not implemented") END IF + ! Do not calculate forces or stress tensor here + use_virial = .FALSE. + calculate_forces = .FALSE. + ! get neighbor lists, we need the full sab_orb list from the ec_env NULLIFY (sab_orb, sac_ppl, sap_ppnl) sab_orb => ec_env%sab_orb @@ -545,11 +1215,6 @@ CONTAINS keep_sparsity=.TRUE., name="CORE HAMILTONIAN MATRIX") ! compute the ppl contribution to the core hamiltonian - CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, particle_set=particle_set, & - atomic_kind_set=atomic_kind_set) - NULLIFY (cell_to_index, virial) - use_virial = .FALSE. - calculate_forces = .FALSE. IF (ASSOCIATED(sac_ppl)) THEN CALL build_core_ppl(ec_env%matrix_h, ec_env%matrix_p, force, & virial, calculate_forces, use_virial, nder, & @@ -590,39 +1255,26 @@ CONTAINS CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_build_ks_matrix' CHARACTER(LEN=default_string_length) :: headline - INTEGER :: handle, iounit, ispin, nspins + INTEGER :: handle, ispin, nspins LOGICAL :: calculate_forces, use_virial REAL(dp) :: eexc, evhxc - TYPE(cp_blacs_env_type), POINTER :: blacs_env - TYPE(cp_logger_type), POINTER :: logger - TYPE(cp_para_env_type), POINTER :: para_env TYPE(dft_control_type), POINTER :: dft_control TYPE(pw_env_type), POINTER :: pw_env TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_r, tau_r, v_rspace, v_tau_rspace TYPE(pw_pool_type), POINTER :: auxbas_pw_pool TYPE(qs_ks_env_type), POINTER :: ks_env TYPE(qs_rho_type), POINTER :: rho - TYPE(virial_type), POINTER :: virial CALL timeset(routineN, handle) - logger => cp_get_default_logger() - IF (logger%para_env%ionode) THEN - iounit = cp_logger_get_default_unit_nr(logger, local=.TRUE.) - ELSE - iounit = -1 - END IF - - calculate_forces = .FALSE. - ! get all information on the electronic density - NULLIFY (rho, ks_env) - CALL get_qs_env(qs_env=qs_env, rho=rho, virial=virial, dft_control=dft_control, & - para_env=para_env, blacs_env=blacs_env, ks_env=ks_env) - + NULLIFY (auxbas_pw_pool, dft_control, ks_env, rho, rho_r, tau_r) + CALL get_qs_env(qs_env=qs_env, & + dft_control=dft_control, & + ks_env=ks_env, & + rho=rho) nspins = dft_control%nspins - use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer) - + calculate_forces = .FALSE. use_virial = .FALSE. ! Kohn-Sham matrix @@ -646,6 +1298,7 @@ CONTAINS NULLIFY (v_rspace, v_tau_rspace) CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho, xc_section=ec_env%xc_section, & vxc_rho=v_rspace, vxc_tau=v_tau_rspace, exc=eexc, just_energy=.FALSE.) + IF (.NOT. ASSOCIATED(v_rspace)) THEN ALLOCATE (v_rspace(nspins)) DO ispin = 1, nspins @@ -676,8 +1329,11 @@ CONTAINS IF (ASSOCIATED(v_tau_rspace)) THEN ! integrate over Tau-potential CALL pw_scale(v_tau_rspace(ispin)%pw, v_tau_rspace(ispin)%pw%pw_grid%dvol) - CALL integrate_v_rspace(v_rspace=v_tau_rspace(ispin)%pw, hmat=ec_env%matrix_ks(ispin, 1), & - qs_env=qs_env, calculate_forces=.FALSE., compute_tau=.TRUE., & + CALL integrate_v_rspace(v_rspace=v_tau_rspace(ispin)%pw, & + hmat=ec_env%matrix_ks(ispin, 1), & + qs_env=qs_env, & + calculate_forces=.FALSE., & + compute_tau=.TRUE., & basis_type="HARRIS", & task_list_external=ec_env%task_list) END IF @@ -725,11 +1381,15 @@ CONTAINS !> Short version of qs_core_hamiltonian !> \param qs_env ... !> \param ec_env ... +!> \param matrix_p ... +!> \param matrix_s ... +!> \param matrix_w ... !> \author Creation (03.2014,JGH) ! ************************************************************************************************** - SUBROUTINE ec_build_core_hamiltonian_force(qs_env, ec_env) + SUBROUTINE ec_build_core_hamiltonian_force(qs_env, ec_env, matrix_p, matrix_s, matrix_w) TYPE(qs_environment_type), POINTER :: qs_env TYPE(energy_correction_type), POINTER :: ec_env + TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s, matrix_w CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_build_core_hamiltonian_force' @@ -762,19 +1422,23 @@ CONTAINS iounit = -1 END IF - IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "ec_build_core_hamiltonian_force - START" - calculate_forces = .TRUE. ! no k-points possible - CALL get_qs_env(qs_env=qs_env, para_env=para_env, dft_control=dft_control) + NULLIFY (cell, dft_control, force, ks_env, para_env, virial) + CALL get_qs_env(qs_env=qs_env, & + cell=cell, & + dft_control=dft_control, & + force=force, & + ks_env=ks_env, & + para_env=para_env, & + virial=virial) nimages = dft_control%nimages IF (nimages /= 1) THEN CPABORT("K-points for Harris functional not implemented") END IF - ! check for virial (currently no stress tensor available) - CALL get_qs_env(qs_env=qs_env, cell=cell, virial=virial) + ! check for virial use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer) fconv = 1.0E-9_dp*pascal/cell%deth @@ -797,15 +1461,14 @@ CONTAINS NULLIFY (scrm) CALL dbcsr_allocate_matrix_set(scrm, 1, 1) ALLOCATE (scrm(1, 1)%matrix) - CALL dbcsr_create(scrm(1, 1)%matrix, template=ec_env%matrix_s(1, 1)%matrix) + CALL dbcsr_create(scrm(1, 1)%matrix, template=matrix_s(1, 1)%matrix) CALL cp_dbcsr_alloc_block_from_nbl(scrm(1, 1)%matrix, sab_orb) - CALL get_qs_env(qs_env=qs_env, ks_env=ks_env, force=force) nder = 1 - IF (SIZE(ec_env%matrix_p, 1) == 2) THEN - CALL dbcsr_add(ec_env%matrix_p(1, 1)%matrix, ec_env%matrix_p(2, 1)%matrix, & + IF (SIZE(matrix_p, 1) == 2) THEN + CALL dbcsr_add(matrix_p(1, 1)%matrix, matrix_p(2, 1)%matrix, & alpha_scalar=1.0_dp, beta_scalar=1.0_dp) - CALL dbcsr_add(ec_env%matrix_w(1, 1)%matrix, ec_env%matrix_w(2, 1)%matrix, & + CALL dbcsr_add(matrix_w(1, 1)%matrix, matrix_w(2, 1)%matrix, & alpha_scalar=1.0_dp, beta_scalar=1.0_dp) END IF @@ -813,12 +1476,12 @@ CONTAINS IF (debug_forces) fodeb(1:3) = force(1)%overlap(1:3, 1) IF (debug_stress .AND. use_virial) stdeb = virial%pv_overlap CALL build_overlap_matrix(ks_env, matrixkp_s=scrm, & - !CALL build_overlap_matrix(ks_env, matrixkp_s=ec_env%matrix_s, & matrix_name="OVERLAP MATRIX", & basis_type_a="HARRIS", & basis_type_b="HARRIS", & sab_nl=sab_orb, calculate_forces=.TRUE., & - matrixkp_p=ec_env%matrix_w) + matrixkp_p=matrix_w) + IF (debug_forces) THEN fodeb(1:3) = force(1)%overlap(1:3, 1) - fodeb(1:3) CALL mp_sum(fodeb, para_env%group) @@ -836,7 +1499,7 @@ CONTAINS matrix_name="KINETIC ENERGY MATRIX", & basis_type="HARRIS", & sab_nl=sab_orb, calculate_forces=.TRUE., & - matrixkp_p=ec_env%matrix_p) + matrixkp_p=matrix_p) IF (debug_forces) THEN fodeb(1:3) = force(1)%kinetic(1:3, 1) - fodeb(1:3) CALL mp_sum(fodeb, para_env%group) @@ -848,19 +1511,20 @@ CONTAINS IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") & 'STRESS| Pout*dT', one_third_sum_diag(stdeb), det_3x3(stdeb) END IF - IF (SIZE(ec_env%matrix_p, 1) == 2) THEN - CALL dbcsr_add(ec_env%matrix_p(1, 1)%matrix, ec_env%matrix_p(2, 1)%matrix, & + IF (SIZE(matrix_p, 1) == 2) THEN + CALL dbcsr_add(matrix_p(1, 1)%matrix, matrix_p(2, 1)%matrix, & alpha_scalar=1.0_dp, beta_scalar=-1.0_dp) END IF ! compute the ppl contribution to the core hamiltonian + NULLIFY (atomic_kind_set, particle_set, qs_kind_set) CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, particle_set=particle_set, & atomic_kind_set=atomic_kind_set) IF (ASSOCIATED(sac_ppl)) THEN IF (calculate_forces .AND. debug_forces) fodeb(1:3) = force(1)%gth_ppl(1:3, 1) IF (debug_stress .AND. use_virial) stdeb = virial%pv_ppl - CALL build_core_ppl(scrm, ec_env%matrix_p, force, & + CALL build_core_ppl(scrm, matrix_p, force, & virial, calculate_forces, use_virial, nder, & qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ppl, & nimages, cell_to_index, "HARRIS") @@ -882,7 +1546,7 @@ CONTAINS IF (ASSOCIATED(sap_ppnl)) THEN IF (calculate_forces .AND. debug_forces) fodeb(1:3) = force(1)%gth_ppnl(1:3, 1) IF (debug_stress .AND. use_virial) stdeb = virial%pv_ppnl - CALL build_core_ppnl(scrm, ec_env%matrix_p, force, & + CALL build_core_ppnl(scrm, matrix_p, force, & virial, calculate_forces, use_virial, nder, & qs_kind_set, atomic_kind_set, particle_set, & sab_orb, sap_ppnl, eps_ppnl, & @@ -948,7 +1612,6 @@ CONTAINS REAL(KIND=dp), DIMENSION(3, 3) :: h_stress, pv_loc, stdeb, sttot TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set TYPE(cell_type), POINTER :: cell - TYPE(cp_blacs_env_type), POINTER :: blacs_env TYPE(cp_logger_type), POINTER :: logger TYPE(cp_para_env_type), POINTER :: para_env TYPE(dbcsr_p_type) :: scrm @@ -982,12 +1645,13 @@ CONTAINS END IF ! get all information on the electronic density - NULLIFY (atomic_kind_set, blacs_env, cell, dft_control, force, ks_env, & - matrix_p, matrix_s, para_env, rho, sab_orb, virial) + NULLIFY (atomic_kind_set, cell, dft_control, force, ks_env, & + matrix_p, matrix_s, para_env, rho, rho_core, rho_g, rho_r, sab_orb, & + tau_r, virial) CALL get_qs_env(qs_env=qs_env, & - blacs_env=blacs_env, & cell=cell, & dft_control=dft_control, & + force=force, & ks_env=ks_env, & para_env=para_env, & rho=rho, & @@ -1021,7 +1685,6 @@ CONTAINS CALL pw_transfer(ec_env%vh_rspace%pw, v_hartree_rspace) - CALL get_qs_env(qs_env=qs_env, force=force) ! calculate output density on grid ! rho_in(R): CALL qs_rho_get(rho, rho_r=rho_r) ! rho_in(G): CALL qs_rho_get(rho, rho_g=rho_g) @@ -1060,9 +1723,6 @@ CONTAINS NULLIFY (tauout_r, tauout_g) IF (dft_control%use_kinetic_energy_density) THEN - IF (use_virial) & - CALL cp_abort(__LOCATION__, "Stress tensor of the Harris functional with "// & - "tau-dependent functionals not implemented!") ALLOCATE (tauout_r(nspins), tauout_g(nspins)) DO ispin = 1, nspins ALLOCATE (tauout_r(ispin)%pw, tauout_g(ispin)%pw) @@ -1160,6 +1820,7 @@ CONTAINS CALL pw_axpy(rhoout_g(ispin)%pw, rhodn_tot_gspace) IF (dft_control%use_kinetic_energy_density) CALL pw_axpy(tau_r(ispin)%pw, tauout_r(ispin)%pw, -1.0_dp) END DO + ! calculate associated hartree potential IF (use_virial) THEN @@ -1214,7 +1875,6 @@ CONTAINS ! Pulay force from Tr P_in (V_H(drho)+ Fxc(rho_in)*drho) ! RHS of CPKS equations: (V_H(drho)+ Fxc(rho_in)*drho)*C0 ! Fxc*drho term - NULLIFY (v_xc) xc_section => ec_env%xc_section IF (use_virial) virial%pv_xc = 0.0_dp @@ -1264,9 +1924,10 @@ CONTAINS DO ispin = 1, nspins CALL pw_scale(v_xc(ispin)%pw, v_xc(ispin)%pw%pw_grid%dvol) CALL pw_axpy(dv_hartree_rspace, v_xc(ispin)%pw) - CALL integrate_v_rspace(qs_env=qs_env, v_rspace=v_xc(ispin)%pw, & + CALL integrate_v_rspace(v_rspace=v_xc(ispin)%pw, & hmat=ec_env%matrix_hz(ispin), & pmat=matrix_p(ispin, 1), & + qs_env=qs_env, & calculate_forces=.TRUE.) END DO @@ -1288,9 +1949,10 @@ CONTAINS DO ispin = 1, nspins CALL pw_scale(v_xc_tau(ispin)%pw, v_xc_tau(ispin)%pw%pw_grid%dvol) - CALL integrate_v_rspace(qs_env=qs_env, v_rspace=v_xc_tau(ispin)%pw, & + CALL integrate_v_rspace(v_rspace=v_xc_tau(ispin)%pw, & hmat=ec_env%matrix_hz(ispin), & pmat=matrix_p(ispin, 1), & + qs_env=qs_env, & compute_tau=.TRUE., & calculate_forces=.TRUE.) END DO @@ -1370,18 +2032,8 @@ CONTAINS calculate_forces=.TRUE., & basis_type="HARRIS", & task_list_external=ec_env%task_list) - - IF (ASSOCIATED(v_tau_rspace)) THEN - ! integrate over Tau-potential - CALL pw_scale(v_tau_rspace(ispin)%pw, v_tau_rspace(ispin)%pw%pw_grid%dvol) - CALL integrate_v_rspace(v_rspace=v_tau_rspace(ispin)%pw, hmat=scrm, & - pmat=ec_env%matrix_p(ispin, 1), & - qs_env=qs_env, calculate_forces=.TRUE., compute_tau=.TRUE., & - basis_type="HARRIS", & - task_list_external=ec_env%task_list) - END IF - END DO + IF (debug_forces) THEN fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3) CALL mp_sum(fodeb, para_env%group) @@ -1399,6 +2051,27 @@ CONTAINS virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc) END IF + IF (ASSOCIATED(v_tau_rspace)) THEN + IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1) + DO ispin = 1, nspins + ! integrate over Tau-potential + CALL pw_scale(v_tau_rspace(ispin)%pw, v_tau_rspace(ispin)%pw%pw_grid%dvol) + CALL integrate_v_rspace(v_rspace=v_tau_rspace(ispin)%pw, & + hmat=scrm, & + pmat=ec_env%matrix_p(ispin, 1), & + qs_env=qs_env, & + calculate_forces=.TRUE., & + compute_tau=.TRUE., & + basis_type="HARRIS", & + task_list_external=ec_env%task_list) + END DO + IF (debug_forces) THEN + fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3) + CALL mp_sum(fodeb, para_env%group) + IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pout*dVhxc_tau ", fodeb + END IF + END IF + ! delete scr matrix CALL dbcsr_release(scrm%matrix) DEALLOCATE (scrm%matrix) @@ -1538,7 +2211,6 @@ CONTAINS IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") & 'STRESS| Explicit total stress ', one_third_sum_diag(stdeb), det_3x3(stdeb) - IF (iounit > 0) WRITE (iounit, *) "Harris_ks_force - END" CALL write_stress_tensor_components(virdeb, iounit, cell) CALL write_stress_tensor(virdeb%pv_virial, iounit, cell, .FALSE.) @@ -1893,19 +2565,25 @@ CONTAINS CALL timeset(routineN, handle) + nspins = SIZE(ec_env%matrix_ks, 1) + DO ispin = 1, nspins + CALL dbcsr_dot(ec_env%matrix_p(ispin, 1)%matrix, ec_env%matrix_s(1, 1)%matrix, trace) + IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T65,F16.10)') 'Tr[PS] ', trace + END DO + + ! Total energy depends on energy correction method SELECT CASE (ec_env%energy_functional) CASE (ec_functional_harris) - nspins = SIZE(ec_env%matrix_ks, 1) + + ! Get energy of "band structure" term eband = 0.0_dp DO ispin = 1, nspins - - CALL dbcsr_dot(ec_env%matrix_p(ispin, 1)%matrix, ec_env%matrix_s(1, 1)%matrix, trace) - IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T65,F16.10)') 'Tr[PS] ', trace - CALL dbcsr_dot(ec_env%matrix_ks(ispin, 1)%matrix, ec_env%matrix_p(ispin, 1)%matrix, trace) eband = eband + trace END DO ec_env%eband = eband + ec_env%efield_nuclear + + ! Add Harris functional "correction" terms ec_env%etotal = ec_env%eband + ec_env%ehartree + ec_env%exc - ec_env%vhxc + ec_env%edispersion IF (unit_nr > 0) THEN WRITE (unit_nr, '(T3,A,T56,F25.15)') "Eband ", ec_env%eband @@ -1916,6 +2594,21 @@ CONTAINS WRITE (unit_nr, '(T3,A,T56,F25.15)') "Etotal Harris Functional ", ec_env%etotal END IF + CASE (ec_functional_dc) + + ! Core hamiltonian energy + CALL calculate_ptrace(ec_env%matrix_h, ec_env%matrix_p, ec_env%ecore, SIZE(ec_env%matrix_p, 1)) + + ec_env%etotal = ec_env%ecore + ec_env%ehartree + ec_env%exc + ec_env%edispersion & + + ec_env%efield_elec + ec_env%efield_nuclear + IF (unit_nr > 0) THEN + WRITE (unit_nr, '(T3,A,T56,F25.15)') "Ecore ", ec_env%ecore + WRITE (unit_nr, '(T3,A,T56,F25.15)') "Ehartree ", ec_env%ehartree + WRITE (unit_nr, '(T3,A,T56,F25.15)') "Exc ", ec_env%exc + WRITE (unit_nr, '(T3,A,T56,F25.15)') "Edisp ", ec_env%edispersion + WRITE (unit_nr, '(T3,A,T56,F25.15)') "Etotal Energy Functional ", ec_env%etotal + END IF + CASE DEFAULT CPASSERT(.FALSE.) @@ -1933,7 +2626,7 @@ CONTAINS !> \param ec_env ... !> \par History !> 2012.07 created [Martin Haeufel] -!> 2016.07 Adapted for Harris functional {JGH] +!> 2016.07 Adapted for Harris functional [JGH] !> \author Martin Haeufel ! ************************************************************************************************** SUBROUTINE ec_build_neighborlist(qs_env, ec_env) @@ -2148,6 +2841,7 @@ CONTAINS NULLIFY (dft_control) CALL get_qs_env(qs_env, dft_control=dft_control) + nspins = dft_control%nspins ec_section => section_vals_get_subs_vals(qs_env%input, "DFT%ENERGY_CORRECTION") print_key => section_vals_get_subs_vals(section_vals=ec_section, & @@ -2210,7 +2904,6 @@ CONTAINS ! direct density contribution CALL ec_efield_integrals(qs_env, ec_env, rcc) ! - nspins = SIZE(ec_env%matrix_p, 1) pdip = 0.0_dp DO ispin = 1, nspins DO idir = 1, 3 @@ -2237,6 +2930,7 @@ CONTAINS DO ispin = 1, nspins DO idir = 1, 3 CALL dbcsr_dot(ec_env%matrix_z(ispin)%matrix, moments(idir)%matrix, tmp) + IF (logger%para_env%mepos == 1) WRITE (*, *) "ec_properties rdip(idir)", idir, tmp rdip(idir) = rdip(idir) + tmp END DO END DO @@ -2658,7 +3352,7 @@ CONTAINS ! Harris functional to use the same basis set CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, nkind=nkind) CALL uppercase(ec_env%basis) - ! Harris basis may differ from ground-state basis only if explicitly added + ! Harris basis only differs from ground-state basis if explicitly added ! thus only two cases that need to be tested ! 1) explicit Harris basis present? IF (ec_env%basis == "HARRIS") THEN diff --git a/src/input_constants.F b/src/input_constants.F index f6a190196e..24acf4b1c0 100644 --- a/src/input_constants.F +++ b/src/input_constants.F @@ -1097,13 +1097,16 @@ MODULE input_constants INTEGER, PARAMETER, PUBLIC :: kg_cholesky = 3001 ! non-scf energy corrections + INTEGER, PARAMETER, PUBLIC :: ec_functional_harris = 2001, & + ec_functional_dc = 2002 + + ! Energy correction solver INTEGER, PARAMETER, PUBLIC :: ec_diagonalization = 1001, & ec_curvy_steps = 1002, & ec_matrix_sign = 1003, & ec_matrix_trs4 = 1004, & ec_matrix_tc2 = 1005, & ec_ot_diag = 1006 - INTEGER, PARAMETER, PUBLIC :: ec_functional_harris = 2001 ! response solver for energy correction INTEGER, PARAMETER, PUBLIC :: ec_ot_atomic = 1, & diff --git a/src/input_cp2k_ec.F b/src/input_cp2k_ec.F index bddebeed95..7e9b65bd3a 100644 --- a/src/input_cp2k_ec.F +++ b/src/input_cp2k_ec.F @@ -20,9 +20,9 @@ MODULE input_cp2k_ec high_print_level USE input_constants, ONLY: & bqb_opt_exhaustive, bqb_opt_normal, bqb_opt_off, bqb_opt_patient, bqb_opt_quick, & - ec_diagonalization, ec_functional_harris, ec_ls_solver, ec_matrix_sign, ec_matrix_tc2, & - ec_matrix_trs4, ec_mo_solver, ec_ot_atomic, ec_ot_diag, ec_ot_gs, kg_cholesky, & - ls_cluster_atomic, ls_cluster_molecular, ls_s_inversion_hotelling, & + ec_diagonalization, ec_functional_dc, ec_functional_harris, ec_ls_solver, ec_matrix_sign, & + ec_matrix_tc2, ec_matrix_trs4, ec_mo_solver, ec_ot_atomic, ec_ot_diag, ec_ot_gs, & + kg_cholesky, ls_cluster_atomic, ls_cluster_molecular, ls_s_inversion_hotelling, & ls_s_inversion_sign_sqrt, ls_s_preconditioner_atomic, ls_s_preconditioner_molecular, & ls_s_preconditioner_none, ls_s_sqrt_ns, ls_s_sqrt_proot, ls_scf_sign_ns, & ls_scf_sign_proot, ot_precond_full_all, ot_precond_full_kinetic, ot_precond_full_single, & @@ -101,9 +101,10 @@ CONTAINS description="Functional used in energy correction", & usage="ENERGY_FUNCTIONAL HARRIS", & default_i_val=ec_functional_harris, & - enum_c_vals=s2a("HARRIS"), & - enum_desc=s2a("Harris functional"), & - enum_i_vals=(/ec_functional_harris/)) + enum_c_vals=s2a("HARRIS", "DCDFT"), & + enum_desc=s2a("Harris functional", & + "Density-corrected DFT"), & + enum_i_vals=(/ec_functional_harris, ec_functional_dc/)) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) diff --git a/src/qs_force.F b/src/qs_force.F index 81b7f83b80..91fa63e500 100644 --- a/src/qs_force.F +++ b/src/qs_force.F @@ -290,7 +290,8 @@ CONTAINS CALL build_xtb_matrices(qs_env=qs_env, para_env=para_env, & calculate_forces=.TRUE.) ELSEIF (perform_ec) THEN - ! + ! Calculates core and grid based forces + CALL energy_correction(qs_env, ec_init=.FALSE., calculate_forces=.TRUE.) ELSE ! Dispersion energy and forces are calculated in qs_energy? CALL build_core_hamiltonian_matrix(qs_env=qs_env, calculate_forces=.TRUE.) @@ -342,10 +343,6 @@ CONTAINS CALL qs_ks_update_qs_env(qs_env, calculate_forces=.TRUE.) END IF - IF (perform_ec) THEN - CALL energy_correction(qs_env, ec_init=.FALSE., calculate_forces=.TRUE.) - END IF - ! Excited state forces CALL excited_state_energy(qs_env, calculate_forces=.TRUE.) @@ -398,6 +395,9 @@ CONTAINS CALL get_qs_env(qs_env, ec_env=ec_env) energy%hartree = ec_env%ehartree energy%exc = ec_env%exc + IF (dft_control%do_admm) THEN + energy%exc_aux_fit = ec_env%exc_aux_fit + END IF END IF DO i = 1, 3 virial%pv_ehartree(i, i) = virial%pv_ehartree(i, i) & diff --git a/src/qs_linres_kernel.F b/src/qs_linres_kernel.F index 3660b629eb..141381a3b8 100644 --- a/src/qs_linres_kernel.F +++ b/src/qs_linres_kernel.F @@ -117,7 +117,10 @@ MODULE qs_linres_kernel PRIVATE ! *** Public subroutines *** + PUBLIC :: apply_xc_admm + PUBLIC :: apply_hfx PUBLIC :: apply_op_2 + PUBLIC :: hfx_matrix CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_linres_kernel' diff --git a/src/qs_linres_methods.F b/src/qs_linres_methods.F index b5ae2635b3..2c9ae9537f 100644 --- a/src/qs_linres_methods.F +++ b/src/qs_linres_methods.F @@ -97,7 +97,6 @@ MODULE qs_linres_methods PUBLIC :: linres_localize, linres_solver PUBLIC :: linres_write_restart, linres_read_restart PUBLIC :: build_dm_response - PUBLIC :: p_env_check_i_alloc CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_linres_methods' diff --git a/src/response_solver.F b/src/response_solver.F index 253222c33b..f696e0f4f8 100644 --- a/src/response_solver.F +++ b/src/response_solver.F @@ -61,8 +61,7 @@ MODULE response_solver kg_tnadd_embed_ri, ls_s_sqrt_ns, ls_s_sqrt_proot, ot_precond_full_all, & ot_precond_full_kinetic, ot_precond_full_single, ot_precond_full_single_inverse, & ot_precond_none, ot_precond_s_inverse, precond_mlp - USE input_section_types, ONLY: section_get_lval,& - section_vals_get,& + USE input_section_types, ONLY: section_vals_get,& section_vals_get_subs_vals,& section_vals_type,& section_vals_val_get @@ -111,8 +110,6 @@ MODULE response_solver set_qs_env USE qs_force_types, ONLY: qs_force_type,& total_qs_force - USE qs_fxc, ONLY: qs_fxc_analytic,& - qs_fxc_fdiff USE qs_integrate_potential, ONLY: integrate_v_core_rspace,& integrate_v_rspace USE qs_kind_types, ONLY: get_qs_kind,& @@ -136,9 +133,7 @@ MODULE response_solver p_env_psi0_changed USE qs_p_env_types, ONLY: p_env_release,& qs_p_env_type - USE qs_rho_types, ONLY: qs_rho_create,& - qs_rho_get,& - qs_rho_set,& + USE qs_rho_types, ONLY: qs_rho_get,& qs_rho_type USE qs_vxc, ONLY: qs_vxc_create USE task_list_types, ONLY: task_list_type @@ -170,10 +165,11 @@ MODULE response_solver CONTAINS ! ************************************************************************************************** -!> \brief Initializes solver of linear response equation for Harris functional Lagrangian +!> \brief Initializes solver of linear response equation for energy correction +!> \brief Call AO or MO based linear response solver for energy correction !> -!> \param qs_env ... -!> \param ec_env Environment of Harris energy correction +!> \param qs_env The quickstep environment +!> \param ec_env The energy correction environment !> !> \date 01.2020 !> \author Fabian Belleflamme @@ -184,8 +180,8 @@ CONTAINS CHARACTER(LEN=*), PARAMETER :: routineN = 'response_calculation' - INTEGER :: handle, homo, ispin, nao, nmo, nspins, & - solver_method, unit_nr + INTEGER :: handle, homo, ispin, nao, nao_aux, nmo, & + nspins, solver_method, unit_nr LOGICAL :: should_stop REAL(KIND=dp) :: focc TYPE(admm_type), POINTER :: admm_env @@ -209,7 +205,8 @@ CONTAINS CALL timeset(routineN, handle) - NULLIFY (dft_control, logger, para_env, sab_orb, solver_section) + NULLIFY (admm_env, dft_control, energy, logger, matrix_s, matrix_s_aux, mo_coeff, mos, para_env, & + rho_ao, sab_orb, solver_section) ! Get useful output unit logger => cp_get_default_logger() @@ -249,8 +246,9 @@ CONTAINS ! Write input section of response solver CALL response_solver_write_input(solver_section, linres_control, unit_nr) - ! Allocate and initialize Z matrix, and energy weighted z-matrix - ! Template is ground-state overlap matrix + ! Allocate and initialize response density matrix Z, + ! and the energy weighted response density matrix + ! Template is the ground-state overlap matrix CALL dbcsr_allocate_matrix_set(ec_env%matrix_wz, nspins) CALL dbcsr_allocate_matrix_set(ec_env%matrix_z, nspins) DO ispin = 1, nspins @@ -266,86 +264,88 @@ CONTAINS CALL dbcsr_set(ec_env%matrix_z(ispin)%matrix, 0.0_dp) END DO + ! MO solver requires MO's of the ground-state calculation, + ! The MOs environment is not allocated if LS-DFT has been used. + ! Introduce MOs here + ! Remark: MOS environment also required for creation of p_env + IF (dft_control%qs_control%do_ls_scf) THEN + + ! Allocate and initialize MO environment + CALL ec_mos_init(qs_env, matrix_s(1)%matrix) + CALL get_qs_env(qs_env, mos=mos, rho=rho) + + ! Get ground-state density matrix + CALL qs_rho_get(rho, rho_ao=rho_ao) + + DO ispin = 1, nspins + CALL get_mo_set(mo_set=mos(ispin)%mo_set, & + mo_coeff=mo_coeff, & + nmo=nmo, nao=nao, homo=homo) + + CALL cp_fm_set_all(mo_coeff, 0.0_dp) + CALL cp_fm_init_random(mo_coeff, nmo) + + CALL cp_fm_create(sv, mo_coeff%matrix_struct, "SV") + ! multiply times PS + ! PS*C(:,1:nomo)+C(:,nomo+1:nmo) (nomo=NINT(nelectron/maxocc)) + CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, mo_coeff, sv, nmo) + CALL cp_dbcsr_sm_fm_multiply(rho_ao(ispin)%matrix, sv, mo_coeff, homo) + CALL cp_fm_release(sv) + ! and ortho the result + CALL make_basis_sm(mo_coeff, nmo, matrix_s(1)%matrix) + + ! rebuilds fm_pools + ! originally done in qs_env_setup, only when mos associated + NULLIFY (blacs_env) + CALL get_qs_env(qs_env, blacs_env=blacs_env) + CALL mpools_rebuild_fm_pools(qs_env%mpools, mos=mos, & + blacs_env=blacs_env, para_env=para_env) + END DO + END IF + + ! initialize p_env + ! Remark: mos environment is needed for this + IF (ASSOCIATED(ec_env%p_env)) THEN + CALL p_env_release(ec_env%p_env) + DEALLOCATE (ec_env%p_env) + NULLIFY (ec_env%p_env) + END IF + ALLOCATE (ec_env%p_env) + CALL p_env_create(ec_env%p_env, qs_env, orthogonal_orbitals=.TRUE., & + linres_control=linres_control) + CALL set_qs_env(qs_env, linres_control=linres_control) + CALL p_env_psi0_changed(ec_env%p_env, qs_env) + ! Total energy overwritten, replace with Etot from energy correction + CALL get_qs_env(qs_env, energy=energy) + energy%total = ec_env%etotal + ! + p_env => ec_env%p_env + ! + CALL dbcsr_allocate_matrix_set(p_env%p1, nspins) + CALL dbcsr_allocate_matrix_set(p_env%w1, nspins) + DO ispin = 1, nspins + ALLOCATE (p_env%p1(ispin)%matrix, p_env%w1(ispin)%matrix) + CALL dbcsr_create(matrix=p_env%p1(ispin)%matrix, template=matrix_s(1)%matrix) + CALL dbcsr_create(matrix=p_env%w1(ispin)%matrix, template=matrix_s(1)%matrix) + CALL cp_dbcsr_alloc_block_from_nbl(p_env%p1(ispin)%matrix, sab_orb) + CALL cp_dbcsr_alloc_block_from_nbl(p_env%w1(ispin)%matrix, sab_orb) + END DO + IF (dft_control%do_admm) THEN + CALL get_admm_env(qs_env%admm_env, matrix_s_aux_fit=matrix_s_aux) + CALL dbcsr_allocate_matrix_set(p_env%p1_admm, nspins) + DO ispin = 1, nspins + ALLOCATE (p_env%p1_admm(ispin)%matrix) + CALL dbcsr_create(p_env%p1_admm(ispin)%matrix, & + template=matrix_s_aux(1)%matrix) + CALL dbcsr_copy(p_env%p1_admm(ispin)%matrix, matrix_s_aux(1)%matrix) + CALL dbcsr_set(p_env%p1_admm(ispin)%matrix, 0.0_dp) + END DO + END IF + ! Choose between MO-solver and AO-solver SELECT CASE (solver_method) CASE (ec_mo_solver) - ! MO solver requires MO's calculated during ground-state, - ! which are not available if linear scaling schemes are used. - ! If ground-state calculation was performed with linear scaling methods, - ! introduce MOs here - IF (dft_control%qs_control%do_ls_scf) THEN - - ! Allocate and initialize MO environment - CALL ec_mos_init(qs_env, matrix_s(1)%matrix) - CALL get_qs_env(qs_env, mos=mos, rho=rho) - - ! Get ground-state density matrix - CALL qs_rho_get(rho, rho_ao=rho_ao) - - DO ispin = 1, nspins - CALL get_mo_set(mo_set=mos(ispin)%mo_set, & - mo_coeff=mo_coeff, & - nmo=nmo, nao=nao, homo=homo) - - CALL cp_fm_set_all(mo_coeff, 0.0_dp) - CALL cp_fm_init_random(mo_coeff, nmo) - - CALL cp_fm_create(sv, mo_coeff%matrix_struct, "SV") - ! multiply times PS - ! PS*C(:,1:nomo)+C(:,nomo+1:nmo) (nomo=NINT(nelectron/maxocc)) - CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, mo_coeff, sv, nmo) - CALL cp_dbcsr_sm_fm_multiply(rho_ao(ispin)%matrix, sv, mo_coeff, homo) - CALL cp_fm_release(sv) - ! and ortho the result - CALL make_basis_sm(mo_coeff, nmo, matrix_s(1)%matrix) - - ! rebuilds fm_pools - ! originally done in qs_env_setup, only when mos associated - CALL get_qs_env(qs_env, blacs_env=blacs_env) - CALL mpools_rebuild_fm_pools(qs_env%mpools, mos=mos, & - blacs_env=blacs_env, para_env=para_env) - END DO - END IF - - ! initialized p_env - IF (ASSOCIATED(ec_env%p_env)) THEN - CALL p_env_release(ec_env%p_env) - DEALLOCATE (ec_env%p_env) - NULLIFY (ec_env%p_env) - END IF - ALLOCATE (ec_env%p_env) - CALL p_env_create(ec_env%p_env, qs_env, orthogonal_orbitals=.TRUE., & - linres_control=linres_control) - CALL set_qs_env(qs_env, linres_control=linres_control) - CALL p_env_psi0_changed(ec_env%p_env, qs_env) - CALL get_qs_env(qs_env, energy=energy) - energy%total = ec_env%etotal - ! - p_env => ec_env%p_env - ! - CALL get_qs_env(qs_env, matrix_s=matrix_s, sab_orb=sab_orb) - CALL dbcsr_allocate_matrix_set(p_env%p1, nspins) - CALL dbcsr_allocate_matrix_set(p_env%w1, nspins) - DO ispin = 1, nspins - ALLOCATE (p_env%p1(ispin)%matrix, p_env%w1(ispin)%matrix) - CALL dbcsr_create(matrix=p_env%p1(ispin)%matrix, template=matrix_s(1)%matrix) - CALL dbcsr_create(matrix=p_env%w1(ispin)%matrix, template=matrix_s(1)%matrix) - CALL cp_dbcsr_alloc_block_from_nbl(p_env%p1(ispin)%matrix, sab_orb) - CALL cp_dbcsr_alloc_block_from_nbl(p_env%w1(ispin)%matrix, sab_orb) - END DO - IF (dft_control%do_admm) THEN - CALL get_admm_env(qs_env%admm_env, matrix_s_aux_fit=matrix_s_aux) - CALL dbcsr_allocate_matrix_set(p_env%p1_admm, nspins) - DO ispin = 1, nspins - ALLOCATE (p_env%p1_admm(ispin)%matrix) - CALL dbcsr_create(p_env%p1_admm(ispin)%matrix, & - template=matrix_s_aux(1)%matrix) - CALL dbcsr_copy(p_env%p1_admm(ispin)%matrix, matrix_s_aux(1)%matrix) - CALL dbcsr_set(p_env%p1_admm(ispin)%matrix, 0.0_dp) - END DO - END IF - ! CPKS vector cpmos - RHS of response equation as Ax + b = 0 (sign of b) ! Sign is changed in linres_solver! ! Projector Q applied in linres_solver! @@ -377,15 +377,6 @@ CONTAINS CALL dbcsr_copy(ec_env%matrix_z(ispin)%matrix, p_env%p1(ispin)%matrix) CALL dbcsr_copy(ec_env%matrix_wz(ispin)%matrix, p_env%w1(ispin)%matrix) END DO - IF (dft_control%do_admm) THEN - CALL dbcsr_allocate_matrix_set(ec_env%z_admm, nspins) - DO ispin = 1, nspins - ALLOCATE (ec_env%z_admm(ispin)%matrix) - CALL dbcsr_create(matrix=ec_env%z_admm(ispin)%matrix, template=matrix_s_aux(1)%matrix) - CALL get_qs_env(qs_env, admm_env=admm_env) - CALL dbcsr_copy(ec_env%z_admm(ispin)%matrix, p_env%p1_admm(ispin)%matrix) - END DO - END IF DO ispin = 1, nspins IF (ASSOCIATED(cpmos(ispin)%matrix)) THEN @@ -395,36 +386,64 @@ CONTAINS END DO DEALLOCATE (cpmos) - ! Get rid of MO environment again - IF (dft_control%qs_control%do_ls_scf) THEN + CASE (ec_ls_solver) + + ! AO ortho solver + CALL ec_response_ao(qs_env=qs_env, & + p_env=p_env, & + matrix_hz=ec_env%matrix_hz, & + matrix_pz=ec_env%matrix_z, & + matrix_wz=ec_env%matrix_wz, & + iounit=unit_nr, & + should_stop=should_stop) + + IF (dft_control%do_admm) THEN + CALL get_qs_env(qs_env, admm_env=admm_env) + CPASSERT(ASSOCIATED(admm_env%work_orb_orb)) + CPASSERT(ASSOCIATED(admm_env%work_aux_orb)) + CPASSERT(ASSOCIATED(admm_env%work_aux_aux)) + nao = admm_env%nao_orb + nao_aux = admm_env%nao_aux_fit DO ispin = 1, nspins - CALL deallocate_mo_set(mos(ispin)%mo_set) + CALL copy_dbcsr_to_fm(ec_env%matrix_z(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, & + admm_env%work_aux_orb) + CALL parallel_gemm('N', 'T', nao_aux, nao_aux, nao, & + 1.0_dp, admm_env%work_aux_orb, admm_env%A, 0.0_dp, & + admm_env%work_aux_aux) + CALL copy_fm_to_dbcsr(admm_env%work_aux_aux, p_env%p1_admm(ispin)%matrix, & + keep_sparsity=.TRUE.) END DO - IF (ASSOCIATED(qs_env%mos)) THEN - DO ispin = 1, SIZE(qs_env%mos) - CALL deallocate_mo_set(qs_env%mos(ispin)%mo_set) - END DO - DEALLOCATE (qs_env%mos) - END IF END IF - CASE (ec_ls_solver) - IF (dft_control%do_admm) THEN - CALL cp_warn(__LOCATION__, "ADMM not possible with AO based response solver. "// & - "Use the MO solver: RESPONSE_SOLVER/METOD MO_SOLVER") - CPABORT("response_calculation") - END IF - ! AO ortho solver - CALL ec_response_ao(qs_env, & - ec_env%matrix_hz, & - ec_env%matrix_z, & - ec_env%matrix_wz, & - unit_nr, & - should_stop) CASE DEFAULT CPABORT("Unknown solver for response equation requested") END SELECT + IF (dft_control%do_admm) THEN + CALL dbcsr_allocate_matrix_set(ec_env%z_admm, nspins) + DO ispin = 1, nspins + ALLOCATE (ec_env%z_admm(ispin)%matrix) + CALL dbcsr_create(matrix=ec_env%z_admm(ispin)%matrix, template=matrix_s_aux(1)%matrix) + CALL get_qs_env(qs_env, admm_env=admm_env) + CALL dbcsr_copy(ec_env%z_admm(ispin)%matrix, p_env%p1_admm(ispin)%matrix) + END DO + END IF + + ! Get rid of MO environment again + IF (dft_control%qs_control%do_ls_scf) THEN + DO ispin = 1, nspins + CALL deallocate_mo_set(mos(ispin)%mo_set) + END DO + IF (ASSOCIATED(qs_env%mos)) THEN + DO ispin = 1, SIZE(qs_env%mos) + CALL deallocate_mo_set(qs_env%mos(ispin)%mo_set) + END DO + DEALLOCATE (qs_env%mos) + END IF + END IF + CALL linres_control_release(linres_control) CALL timestop(handle) @@ -583,6 +602,7 @@ CONTAINS DO ispin = 1, nspins CALL dbcsr_scale(p_env%p1(ispin)%matrix, 0.5_dp) END DO + IF (dft_control%do_admm) THEN CALL get_qs_env(qs_env, admm_env=admm_env) CPASSERT(ASSOCIATED(admm_env%work_orb_orb)) @@ -794,19 +814,20 @@ CONTAINS !> \param matrix_wz Energy-weighted linear response density !> \param zehartree Hartree volume response contribution to stress tensor !> \param zexc XC volume response contribution to stress tensor +!> \param zexc_aux_fit ADMM XC volume response contribution to stress tensor !> \param rhopz_r Response density on real space grid !> \param p_env ... !> \param ex_env ... ! ************************************************************************************************** SUBROUTINE response_force(qs_env, vh_rspace, vxc_rspace, vtau_rspace, vadmm_rspace, & matrix_hz, matrix_pz, matrix_pz_admm, matrix_wz, & - zehartree, zexc, rhopz_r, p_env, ex_env) + zehartree, zexc, zexc_aux_fit, rhopz_r, p_env, ex_env) TYPE(qs_environment_type), POINTER :: qs_env TYPE(pw_type), INTENT(IN) :: vh_rspace TYPE(pw_p_type), DIMENSION(:), POINTER :: vxc_rspace, vtau_rspace, vadmm_rspace TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_hz, matrix_pz, matrix_pz_admm, & matrix_wz - REAL(KIND=dp), OPTIONAL :: zehartree, zexc + REAL(KIND=dp), OPTIONAL :: zehartree, zexc, zexc_aux_fit TYPE(pw_p_type), DIMENSION(:), INTENT(IN), & OPTIONAL :: rhopz_r TYPE(qs_p_env_type), OPTIONAL :: p_env @@ -818,10 +839,11 @@ CONTAINS nao, nao_aux, natom, nder, nimages, & nspins INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index - LOGICAL :: deriv2_analytic, distribute_fock_matrix, do_ex, do_hfx, hfx_treat_lsd_in_core, & - resp_only, s_mstruct_changed, use_virial + LOGICAL :: distribute_fock_matrix, do_ex, do_hfx, & + hfx_treat_lsd_in_core, resp_only, & + s_mstruct_changed, use_virial REAL(KIND=dp) :: eh1, ehartree, ekin_mol, eps_filter, & - eps_ppnl, exc, fconv, focc + eps_ppnl, exc, exc_aux_fit, fconv, focc REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: ftot1, ftot2, ftot3 REAL(KIND=dp), DIMENSION(3) :: fodeb REAL(KIND=dp), DIMENSION(3, 3) :: h_stress, pv_loc, stdeb, sttot, sttot2 @@ -856,7 +878,7 @@ CONTAINS TYPE(qs_force_type), DIMENSION(:), POINTER :: force TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set TYPE(qs_ks_env_type), POINTER :: ks_env - TYPE(qs_rho_type), POINTER :: rho, rho_aux_fit, rhoz_aux + TYPE(qs_rho_type), POINTER :: rho, rho_aux_fit TYPE(section_vals_type), POINTER :: hfx_section, xc_section TYPE(task_list_type), POINTER :: task_list_aux_fit TYPE(virial_type), POINTER :: virial @@ -877,7 +899,7 @@ CONTAINS CPASSERT(PRESENT(p_env)) END IF - NULLIFY (ks_env, sab_orb, sac_ppl, sap_ppnl, virial) + NULLIFY (ks_env, sab_orb, sac_ppl, sap_ppnl, tauz_r, virial) CALL get_qs_env(qs_env=qs_env, & cell=cell, & force=force, & @@ -1073,7 +1095,7 @@ CONTAINS CALL total_qs_force(ftot2, force, atomic_kind_set) fodeb(1:3) = ftot2(1:3, 1) - ftot1(1:3, 1) CALL mp_sum(fodeb, para_env%group) - IF (iounit > 0) WRITE (iounit, "(T3,A,T30,3F16.8)") "DEBUG:: Force Pz*dHcore", fodeb + IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pz*dHcore", fodeb END IF IF (debug_stress .AND. use_virial) THEN stdeb = fconv*(virial%pv_virial - sttot) @@ -1123,9 +1145,6 @@ CONTAINS END IF IF (ASSOCIATED(vtau_rspace)) THEN - IF (use_virial) & - CALL cp_abort(__LOCATION__, "Stress tensor of the Harris functional with "// & - "tau-dependent functionals not implemented!") IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1) IF (debug_stress .AND. use_virial) stdeb = virial%pv_virial DO ispin = 1, nspins @@ -1264,7 +1283,7 @@ CONTAINS h_stress(:, :) = 0.0_dp ! calculate associated hartree potential ! This term appears twice in the derivation of the equations - ! v_H{n_in]*n_z and v_H[n_z]*n_in + ! v_H[n_in]*n_z and v_H[n_z]*n_in ! due to symmetry we only need to call this routine once, ! and count the Volume and Green function contribution ! which is stored in h_stress twice @@ -1520,19 +1539,18 @@ CONTAINS CALL total_qs_force(ftot3, force, atomic_kind_set) fodeb(1:3) = ftot3(1:3, 1) - ftot2(1:3, 1) CALL mp_sum(fodeb, para_env%group) - IF (iounit > 0) WRITE (iounit, "(T3,A,T30,3F16.8)") "DEBUG:: Force Pin*V(rhoz)", fodeb + IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pin*V(rhoz)", fodeb END IF CALL dbcsr_deallocate_matrix_set(scrm) - ! Stress-tensor volume contributions - ! These need to be applied at the end of qs_force - IF (use_virial) THEN - ! Adding mixed Hartree energy twice, due to symmetry - zehartree = zehartree + 2.0_dp*ehartree - zexc = zexc + exc - END IF + ! ----------------------------------------- + ! Apply ADMM exchange correction + ! ----------------------------------------- IF (dft_control%do_admm) THEN + ! volume term + exc_aux_fit = 0.0_dp + IF (qs_env%admm_env%aux_exch_func == do_admm_aux_exch_func_none) THEN ! nothing to do NULLIFY (mpz, mhz, mhx, mhy) @@ -1568,7 +1586,13 @@ CONTAINS END DO ! xc_section => admm_env%xc_section_aux + ! Stress-tensor: integration contribution direct term + ! int Pz*v_xc[rho_admm] + IF (use_virial) THEN + pv_loc = virial%pv_virial + END IF IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1) + IF (debug_stress .AND. use_virial) stdeb = virial%pv_virial DO ispin = 1, nspins CALL integrate_v_rspace(v_rspace=vadmm_rspace(ispin)%pw, & hmat=mhx(ispin, 1), pmat=mpz(ispin, 1), & @@ -1580,6 +1604,16 @@ CONTAINS CALL mp_sum(fodeb, para_env%group) IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pz*Vxc(rho_admm)", fodeb END IF + IF (debug_stress .AND. use_virial) THEN + stdeb = fconv*(virial%pv_virial - pv_loc) + CALL mp_sum(stdeb, para_env%group) + IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") & + 'STRESS| INT 1st Pz*dVxc(rho_admm) ', one_third_sum_diag(stdeb), det_3x3(stdeb) + END IF + ! Stress-tensor Pz_admm*v_xc[rho_admm] + IF (use_virial) THEN + virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc) + END IF ! NULLIFY (rho_g_aux, rho_r_aux, tau_r_aux, rhoz_g_aux, rhoz_r_aux) CALL qs_rho_get(rho_aux_fit, rho_r=rho_r_aux, rho_g=rho_g_aux, tau_r=tau_r_aux) @@ -1600,21 +1634,59 @@ CONTAINS task_list_external=task_list_aux_fit) END DO ! + ! Add ADMM volume contribution to stress tensor + IF (use_virial) THEN + + ! Stress tensor volume term: \int v_xc[n_in_admm]*n_z_admm + ! vadmm_rspace already scaled, we need to unscale it! + DO ispin = 1, nspins + exc_aux_fit = exc_aux_fit + pw_integral_ab(rhoz_r_aux(ispin)%pw, vadmm_rspace(ispin)%pw)/ & + vadmm_rspace(ispin)%pw%pw_grid%dvol + END DO + + IF (debug_stress) THEN + stdeb = -1.0_dp*fconv*exc_aux_fit + IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T43,2(1X,ES19.11))") & + 'STRESS| VOL 1st eps_XC[n_in_admm]*n_z_admm', one_third_sum_diag(stdeb), det_3x3(stdeb) + END IF + + END IF + ! NULLIFY (v_xc) - deriv2_analytic = section_get_lval(xc_section, "2ND_DERIV_ANALYTICAL") - IF (deriv2_analytic) THEN - CALL qs_fxc_analytic(rho_aux_fit, rhoz_r_aux, Null(), xc_section, auxbas_pw_pool, .FALSE., v_xc, v_xc_tau) - ELSE - CPABORT("NYA 00005") - NULLIFY (rhoz_aux) - CALL qs_rho_create(rhoz_aux) - CALL qs_rho_set(rhoz_aux, rho_r=rhoz_r_aux, rho_g=rhoz_g_aux) - CALL qs_fxc_fdiff(ks_env, rho_aux_fit, rhoz_aux, xc_section, 6, .FALSE., v_xc, tau_r) - DEALLOCATE (rhoz_aux) + + IF (use_virial) virial%pv_xc = 0.0_dp + + CALL create_kernel(qs_env=qs_env, & + vxc=v_xc, & + vxc_tau=v_xc_tau, & + rho=rho_aux_fit, & + rho1_r=rhoz_r_aux, & + rho1_g=rhoz_g_aux, & + tau1_r=tau_r_aux, & + xc_section=xc_section, & + compute_virial=use_virial, & + virial_xc=virial%pv_xc) + + ! Stress-tensor ADMM-kernel GGA contribution + IF (use_virial) THEN + virial%pv_exc = virial%pv_exc + virial%pv_xc + virial%pv_virial = virial%pv_virial + virial%pv_xc + END IF + + IF (debug_stress .AND. use_virial) THEN + stdeb = 1.0_dp*fconv*virial%pv_xc + CALL mp_sum(stdeb, para_env%group) + IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") & + 'STRESS| GGA 2nd Pin_admm*dK*rhoz_admm', one_third_sum_diag(stdeb), det_3x3(stdeb) END IF ! CALL qs_rho_get(rho_aux_fit, rho_ao_kp=matrix_p) + ! Stress-tensor Pin*dK*rhoz_admm + IF (use_virial) THEN + virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc) + END IF IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1) + IF (debug_stress .AND. use_virial) stdeb = virial%pv_virial DO ispin = 1, nspins CALL dbcsr_set(mhy(ispin, 1)%matrix, 0.0_dp) CALL pw_scale(v_xc(ispin)%pw, v_xc(ispin)%pw%pw_grid%dvol) @@ -1628,6 +1700,16 @@ CONTAINS CALL mp_sum(fodeb, para_env%group) IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pin*dK*rhoz_admm ", fodeb END IF + IF (debug_stress .AND. use_virial) THEN + stdeb = fconv*(virial%pv_virial - pv_loc) + CALL mp_sum(stdeb, para_env%group) + IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") & + 'STRESS| INT 2nd Pin*dK*rhoz_admm ', one_third_sum_diag(stdeb), det_3x3(stdeb) + END IF + ! Stress-tensor Pin*dK*rhoz_admm + IF (use_virial) THEN + virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc) + END IF DO ispin = 1, nspins CALL pw_pool_give_back_pw(auxbas_pw_pool, v_xc(ispin)%pw) CALL pw_pool_give_back_pw(auxbas_pw_pool, rhoz_r_aux(ispin)%pw) @@ -1654,8 +1736,12 @@ CONTAINS CALL dbcsr_release(dbwork) DEALLOCATE (dbwork) CALL dbcsr_deallocate_matrix_set(mpz) - END IF - END IF + END IF ! qs_env%admm_env%aux_exch_func == do_admm_aux_exch_func_none + END IF ! do_admm + + ! ----------------------------------------- + ! HFX + ! ----------------------------------------- ! HFX hfx_section => section_vals_get_subs_vals(xc_section, "HF") @@ -1667,10 +1753,15 @@ CONTAINS i_rep_section=1) mspin = 1 IF (hfx_treat_lsd_in_core) mspin = nspins + IF (use_virial) virial%pv_fock_4c = 0.0_dp ! CALL get_qs_env(qs_env=qs_env, rho=rho, x_data=x_data, & s_mstruct_changed=s_mstruct_changed) distribute_fock_matrix = .TRUE. + + ! ----------------------------------------- + ! HFX-ADMM + ! ----------------------------------------- IF (dft_control%do_admm) THEN CALL get_qs_env(qs_env=qs_env, admm_env=admm_env) CALL get_admm_env(admm_env, matrix_s_aux_fit=scrm, rho_aux_fit=rho_aux_fit) @@ -1773,6 +1864,9 @@ CONTAINS END IF DEALLOCATE (mpd) ELSE + ! ----------------------------------------- + ! conventional HFX + ! ----------------------------------------- ALLOCATE (mpz(nspins, 1), mhz(nspins, 1)) DO ispin = 1, nspins mhz(ispin, 1)%matrix => matrix_hz(ispin)%matrix @@ -1796,9 +1890,16 @@ CONTAINS DEALLOCATE (mhz, mpz) END IF + ! ----------------------------------------- + ! HFX FORCES + ! ----------------------------------------- + resp_only = .TRUE. IF (debug_forces) fodeb(1:3) = force(1)%fock_4c(1:3, 1) IF (dft_control%do_admm) THEN + ! ----------------------------------------- + ! HFX-ADMM FORCES + ! ----------------------------------------- CALL qs_rho_get(rho_aux_fit, rho_ao_kp=matrix_p) NULLIFY (matrix_pza) CALL dbcsr_allocate_matrix_set(matrix_pza, nspins) @@ -1826,6 +1927,9 @@ CONTAINS END IF CALL dbcsr_deallocate_matrix_set(matrix_pza) ELSE + ! ----------------------------------------- + ! conventional HFX FORCES + ! ----------------------------------------- CALL qs_rho_get(rho, rho_ao_kp=matrix_p) IF (x_data(1, 1)%do_hfx_ri) THEN @@ -1837,12 +1941,37 @@ CONTAINS CALL derivatives_four_center(qs_env, matrix_p, mpa, hfx_section, para_env, & 1, use_virial, resp_only=resp_only) END IF + END IF ! do_admm + + IF (use_virial) THEN + virial%pv_exx = virial%pv_exx - virial%pv_fock_4c + virial%pv_virial = virial%pv_virial - virial%pv_fock_4c + virial%pv_calculate = .FALSE. END IF + IF (debug_forces) THEN fodeb(1:3) = force(1)%fock_4c(1:3, 1) - fodeb(1:3) CALL mp_sum(fodeb, para_env%group) IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pz*hfx ", fodeb END IF + IF (debug_stress .AND. use_virial) THEN + stdeb = -1.0_dp*fconv*virial%pv_fock_4c + CALL mp_sum(stdeb, para_env%group) + IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") & + 'STRESS| Pz*hfx ', one_third_sum_diag(stdeb), det_3x3(stdeb) + END IF + END IF ! do_hfx + + ! Stress-tensor volume contributions + ! These need to be applied at the end of qs_force + IF (use_virial) THEN + ! Adding mixed Hartree energy twice, due to symmetry + zehartree = zehartree + 2.0_dp*ehartree + zexc = zexc + exc + ! ADMM contribution handled differently in qs_force + IF (dft_control%do_admm) THEN + zexc_aux_fit = zexc_aux_fit + exc_aux_fit + END IF END IF ! Overlap matrix @@ -1900,10 +2029,10 @@ CONTAINS CALL total_qs_force(ftot2, force, atomic_kind_set) fodeb(1:3) = ftot2(1:3, 1) - ftot1(1:3, 1) CALL mp_sum(fodeb, para_env%group) - IF (iounit > 0) WRITE (iounit, "(T3,A,T30,3F16.8)") "DEBUG:: Response Force", fodeb + IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Response Force", fodeb fodeb(1:3) = ftot2(1:3, 1) CALL mp_sum(fodeb, para_env%group) - IF (iounit > 0) WRITE (iounit, "(T3,A,T30,3F16.8)") "DEBUG:: Total Force ", fodeb + IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Total Force ", fodeb DEALLOCATE (ftot1, ftot2, ftot3) END IF @@ -2325,6 +2454,12 @@ CONTAINS END DO DEALLOCATE (v_admm_rspace) END IF + IF (ASSOCIATED(v_admm_tau_rspace)) THEN + DO ispin = 1, nspins + CALL pw_pool_give_back_pw(auxbas_pw_pool, v_admm_tau_rspace(ispin)%pw) + END DO + DEALLOCATE (v_admm_tau_rspace) + END IF CALL timestop(handle) diff --git a/tests/QS/regtest-dcdft-force/N2_t01.inp b/tests/QS/regtest-dcdft-force/N2_t01.inp new file mode 100644 index 0000000000..65b04b21e0 --- /dev/null +++ b/tests/QS/regtest-dcdft-force/N2_t01.inp @@ -0,0 +1,60 @@ +&FORCE_EVAL + METHOD Quickstep + &DFT + &PRINT + &DERIVATIVES + &END + &END + BASIS_SET_FILE_NAME BASIS_SET + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 300 + &END MGRID + &QS + EPS_DEFAULT 1.E-12 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL DCDFT + HARRIS_BASIS ORBITAL + &RESPONSE_SOLVER + METHOD MO_SOLVER + PRECONDITIONER FULL_SINGLE_INVERSE + &END + &XC + &XC_FUNCTIONAL PADE + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-7 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL PADE + &END + &END XC + &END DFT + &SUBSYS + &CELL + ABC 6.0 6.0 6.0 + &END CELL + &COORD + N 0.400000 0.000000 0.500000 + N -0.400000 0.000000 -0.500000 + &END COORD + &KIND N + BASIS_SET ORB DZV-GTH-BLYP + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE GEO_OPT + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &GEO_OPT + MAX_ITER 2 + &END +&END diff --git a/tests/QS/regtest-dcdft-force/N2_t02.inp b/tests/QS/regtest-dcdft-force/N2_t02.inp new file mode 100644 index 0000000000..e748cb2a8f --- /dev/null +++ b/tests/QS/regtest-dcdft-force/N2_t02.inp @@ -0,0 +1,60 @@ +&FORCE_EVAL + METHOD Quickstep + &DFT + &PRINT + &DERIVATIVES + &END + &END + BASIS_SET_FILE_NAME BASIS_SET + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 300 + &END MGRID + &QS + EPS_DEFAULT 1.E-12 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL DCDFT + HARRIS_BASIS ORBITAL + &RESPONSE_SOLVER + METHOD MO_SOLVER + PRECONDITIONER FULL_SINGLE_INVERSE + &END + &XC + &XC_FUNCTIONAL PBE + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-6 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL PADE + &END + &END XC + &END DFT + &SUBSYS + &CELL + ABC 6.0 6.0 6.0 + &END CELL + &COORD + N 0.400000 0.000000 0.500000 + N -0.400000 0.000000 -0.500000 + &END COORD + &KIND N + BASIS_SET ORB DZV-GTH-BLYP + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE GEO_OPT + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &GEO_OPT + MAX_ITER 2 + &END +&END diff --git a/tests/QS/regtest-dcdft-force/N2_t03.inp b/tests/QS/regtest-dcdft-force/N2_t03.inp new file mode 100644 index 0000000000..ae1b1911a3 --- /dev/null +++ b/tests/QS/regtest-dcdft-force/N2_t03.inp @@ -0,0 +1,60 @@ +&FORCE_EVAL + METHOD Quickstep + &DFT + &PRINT + &DERIVATIVES + &END + &END + BASIS_SET_FILE_NAME BASIS_SET + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 200 + &END MGRID + &QS + EPS_DEFAULT 1.E-10 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL DCDFT + HARRIS_BASIS ORBITAL + &RESPONSE_SOLVER + METHOD MO_SOLVER + PRECONDITIONER FULL_SINGLE_INVERSE + &END + &XC + &XC_FUNCTIONAL PADE + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-6 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL PBE + &END + &END XC + &END DFT + &SUBSYS + &CELL + ABC 5.0 5.0 5.0 + &END CELL + &COORD + N 0.400000 0.000000 0.500000 + N -0.400000 0.000000 -0.500000 + &END COORD + &KIND N + BASIS_SET ORB DZV-GTH-BLYP + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE GEO_OPT + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &GEO_OPT + MAX_ITER 1 + &END +&END diff --git a/tests/QS/regtest-dcdft-force/N2_t04.inp b/tests/QS/regtest-dcdft-force/N2_t04.inp new file mode 100644 index 0000000000..c8cd303e30 --- /dev/null +++ b/tests/QS/regtest-dcdft-force/N2_t04.inp @@ -0,0 +1,61 @@ +&FORCE_EVAL + METHOD Quickstep + &DFT + &PRINT + &DERIVATIVES + &END + &END + BASIS_SET_FILE_NAME BASIS_SET + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 200 + REL_CUTOFF 40 + &END MGRID + &QS + EPS_DEFAULT 1.E-10 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL DCDFT + HARRIS_BASIS ORBITAL + &RESPONSE_SOLVER + METHOD MO_SOLVER + PRECONDITIONER FULL_SINGLE_INVERSE + &END + &XC + &XC_FUNCTIONAL PBE + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-6 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL PBE + &END + &END XC + &END DFT + &SUBSYS + &CELL + ABC 5.0 5.0 5.0 + &END CELL + &COORD + N 0.400000 0.000000 0.500000 + N -0.400000 0.000000 -0.500000 + &END COORD + &KIND N + BASIS_SET ORB DZV-GTH-BLYP + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE GEO_OPT + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &GEO_OPT + MAX_ITER 2 + &END +&END diff --git a/tests/QS/regtest-dcdft-force/N2_t05.inp b/tests/QS/regtest-dcdft-force/N2_t05.inp new file mode 100644 index 0000000000..ed43005d96 --- /dev/null +++ b/tests/QS/regtest-dcdft-force/N2_t05.inp @@ -0,0 +1,58 @@ +&FORCE_EVAL + METHOD Quickstep + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 200 + &END MGRID + &QS + EPS_DEFAULT 1.E-10 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL DCDFT + HARRIS_BASIS ORBITAL + &RESPONSE_SOLVER + METHOD MO_SOLVER + PRECONDITIONER FULL_SINGLE_INVERSE + &END + &XC + &XC_FUNCTIONAL + &MGGA_X_TPSS + &END + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-6 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL PADE + &END + &END XC + &END DFT + &SUBSYS + &CELL + ABC 5.0 5.0 5.0 + &END CELL + &COORD + N 0.400000 0.000000 0.500000 + N -0.400000 0.000000 -0.500000 + &END COORD + &KIND N + BASIS_SET ORB DZVP-GTH + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE GEO_OPT + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &GEO_OPT + MAX_ITER 1 + &END +&END diff --git a/tests/QS/regtest-dcdft-force/N2_t06.inp b/tests/QS/regtest-dcdft-force/N2_t06.inp new file mode 100644 index 0000000000..d76f586c69 --- /dev/null +++ b/tests/QS/regtest-dcdft-force/N2_t06.inp @@ -0,0 +1,63 @@ +&FORCE_EVAL + METHOD Quickstep + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 200 + &END MGRID + &QS + EPS_DEFAULT 1.E-10 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL DCDFT + HARRIS_BASIS ORBITAL + &RESPONSE_SOLVER + METHOD MO_SOLVER + PRECONDITIONER FULL_SINGLE_INVERSE + &END + &XC + &XC_FUNCTIONAL PADE + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-6 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL + &MGGA_X_TPSS + &END + &END + &END XC + &POISSON + PERIODIC NONE + POISSON_SOLVER MT + &END + &END DFT + &SUBSYS + &CELL + ABC 4.0 4.0 4.0 + PERIODIC NONE + &END CELL + &COORD + N 0.400000 0.000000 0.500000 + N -0.400000 0.000000 -0.500000 + &END COORD + &KIND N + BASIS_SET ORB DZVP-GTH + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE GEO_OPT + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &GEO_OPT + MAX_ITER 1 + &END +&END diff --git a/tests/QS/regtest-dcdft-force/N2_t07.inp b/tests/QS/regtest-dcdft-force/N2_t07.inp new file mode 100644 index 0000000000..e5425aeaef --- /dev/null +++ b/tests/QS/regtest-dcdft-force/N2_t07.inp @@ -0,0 +1,63 @@ +&FORCE_EVAL + METHOD Quickstep + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 200 + &END MGRID + &QS + EPS_DEFAULT 1.E-10 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL DCDFT + HARRIS_BASIS ORBITAL + &RESPONSE_SOLVER + METHOD AO_ORTHO + PRECONDITIONER MULTI_LEVEL + &END + &XC + &XC_FUNCTIONAL + &MGGA_X_BR89 + &END + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-6 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL PADE + &END + &END XC + &POISSON + PERIODIC NONE + POISSON_SOLVER MT + &END + &END DFT + &SUBSYS + &CELL + ABC 4.0 4.0 4.0 + PERIODIC NONE + &END CELL + &COORD + N 0.400000 0.000000 0.500000 + N -0.400000 0.000000 -0.500000 + &END COORD + &KIND N + BASIS_SET ORB DZVP-GTH + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE GEO_OPT + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &GEO_OPT + MAX_ITER 1 + &END +&END diff --git a/tests/QS/regtest-dcdft-force/N2_t08.inp b/tests/QS/regtest-dcdft-force/N2_t08.inp new file mode 100644 index 0000000000..cc05144da2 --- /dev/null +++ b/tests/QS/regtest-dcdft-force/N2_t08.inp @@ -0,0 +1,64 @@ +&FORCE_EVAL + METHOD Quickstep + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 200 + &END MGRID + &QS + EPS_DEFAULT 1.E-10 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL DCDFT + HARRIS_BASIS ORBITAL + &RESPONSE_SOLVER + METHOD AO_ORTHO + PRECONDITIONER MULTI_LEVEL + MAX_ITER 2 + &END + &XC + &XC_FUNCTIONAL PADE + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-6 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL + &MGGA_X_BR89 + &END + &END + &END XC + &POISSON + PERIODIC NONE + POISSON_SOLVER MT + &END + &END DFT + &SUBSYS + &CELL + ABC 4.0 4.0 4.0 + PERIODIC NONE + &END CELL + &COORD + N 0.400000 0.000000 0.500000 + N -0.400000 0.000000 -0.500000 + &END COORD + &KIND N + BASIS_SET ORB DZVP-GTH + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE GEO_OPT + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &GEO_OPT + MAX_ITER 1 + &END +&END diff --git a/tests/QS/regtest-dcdft-force/N2_t09.inp b/tests/QS/regtest-dcdft-force/N2_t09.inp new file mode 100644 index 0000000000..56f723ac7a --- /dev/null +++ b/tests/QS/regtest-dcdft-force/N2_t09.inp @@ -0,0 +1,66 @@ +&FORCE_EVAL + METHOD Quickstep + &DFT + &PRINT + &DERIVATIVES + &END + &END + BASIS_SET_FILE_NAME BASIS_SET + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 200 + REL_CUTOFF 40 + &END MGRID + &QS + EPS_DEFAULT 1.E-10 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL DCDFT + HARRIS_BASIS ORBITAL + &RESPONSE_SOLVER + METHOD MO_SOLVER + PRECONDITIONER FULL_SINGLE_INVERSE + &END + &XC + &XC_FUNCTIONAL PBE + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-6 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL NONE + &END + &HF + &SCREENING + EPS_SCHWARZ 1.0E-14 + &END + &END + &END XC + &END DFT + &SUBSYS + &CELL + ABC 5.0 5.0 5.0 + &END CELL + &COORD + N 0.400000 0.000000 0.500000 + N -0.400000 0.000000 -0.500000 + &END COORD + &KIND N + BASIS_SET ORB DZV-GTH-BLYP + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE GEO_OPT + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &GEO_OPT + MAX_ITER 1 + &END +&END diff --git a/tests/QS/regtest-dcdft-force/N2_t10.inp b/tests/QS/regtest-dcdft-force/N2_t10.inp new file mode 100644 index 0000000000..ba343d7158 --- /dev/null +++ b/tests/QS/regtest-dcdft-force/N2_t10.inp @@ -0,0 +1,67 @@ +&FORCE_EVAL + METHOD Quickstep + &DFT + &PRINT + &DERIVATIVES + &END + &END + BASIS_SET_FILE_NAME BASIS_SET + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 200 + REL_CUTOFF 20 + &END MGRID + &QS + EPS_DEFAULT 1.E-10 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL DCDFT + HARRIS_BASIS ORBITAL + &RESPONSE_SOLVER + METHOD AO_ORTHO + PRECONDITIONER MULTI_LEVEL + MAX_ITER 2 + &END + &XC + &XC_FUNCTIONAL PBE + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-6 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL NONE + &END + &HF + &SCREENING + EPS_SCHWARZ 1.0E-14 + &END + &END + &END XC + &END DFT + &SUBSYS + &CELL + ABC 5.0 5.0 5.0 + &END CELL + &COORD + N 0.400000 0.000000 0.500000 + N -0.400000 0.000000 -0.500000 + &END COORD + &KIND N + BASIS_SET ORB DZV-GTH-BLYP + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE GEO_OPT + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &GEO_OPT + MAX_ITER 1 + &END +&END diff --git a/tests/QS/regtest-dcdft-force/N2_t11.inp b/tests/QS/regtest-dcdft-force/N2_t11.inp new file mode 100644 index 0000000000..5b1dcc779f --- /dev/null +++ b/tests/QS/regtest-dcdft-force/N2_t11.inp @@ -0,0 +1,60 @@ +&FORCE_EVAL + METHOD Quickstep + &DFT + &PRINT + &DERIVATIVES + &END + &END + BASIS_SET_FILE_NAME BASIS_SET + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 200 + &END MGRID + &QS + EPS_DEFAULT 1.E-10 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL DCDFT + HARRIS_BASIS ORBITAL + &RESPONSE_SOLVER + METHOD MO_SOLVER + PRECONDITIONER FULL_SINGLE_INVERSE + &END + &XC + &XC_FUNCTIONAL PADE + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-6 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL PBE0 + &END + &END XC + &END DFT + &SUBSYS + &CELL + ABC 4.0 4.0 4.0 + &END CELL + &COORD + N 0.400000 0.000000 0.500000 + N -0.400000 0.000000 -0.500000 + &END COORD + &KIND N + BASIS_SET ORB DZV-GTH-BLYP + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE GEO_OPT + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &GEO_OPT + MAX_ITER 1 + &END +&END diff --git a/tests/QS/regtest-dcdft-force/N2_t12.inp b/tests/QS/regtest-dcdft-force/N2_t12.inp new file mode 100644 index 0000000000..f13e153d82 --- /dev/null +++ b/tests/QS/regtest-dcdft-force/N2_t12.inp @@ -0,0 +1,61 @@ +&FORCE_EVAL + METHOD Quickstep + &DFT + &PRINT + &DERIVATIVES + &END + &END + BASIS_SET_FILE_NAME BASIS_SET + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 200 + &END MGRID + &QS + EPS_DEFAULT 1.E-10 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL DCDFT + HARRIS_BASIS ORBITAL + &RESPONSE_SOLVER + METHOD AO_ORTHO + PRECONDITIONER MULTI_LEVEL + MAX_ITER 2 + &END + &XC + &XC_FUNCTIONAL PADE + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-6 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL PBE0 + &END + &END XC + &END DFT + &SUBSYS + &CELL + ABC 4.0 4.0 4.0 + &END CELL + &COORD + N 0.400000 0.000000 0.500000 + N -0.400000 0.000000 -0.500000 + &END COORD + &KIND N + BASIS_SET ORB DZV-GTH-BLYP + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE GEO_OPT + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &GEO_OPT + MAX_ITER 1 + &END +&END diff --git a/tests/QS/regtest-dcdft-force/N2_t13.inp b/tests/QS/regtest-dcdft-force/N2_t13.inp new file mode 100644 index 0000000000..c5d0e86538 --- /dev/null +++ b/tests/QS/regtest-dcdft-force/N2_t13.inp @@ -0,0 +1,79 @@ +&FORCE_EVAL + METHOD Quickstep + &DFT +# &PRINT +# &DERIVATIVES +# &END +# &END + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME GTH_POTENTIALS + BASIS_SET_FILE_NAME BASIS_ADMM + &AUXILIARY_DENSITY_MATRIX_METHOD + ADMM_PURIFICATION_METHOD NONE + EXCH_CORRECTION_FUNC NONE + EXCH_SCALING_MODEL NONE + METHOD BASIS_PROJECTION + &END + &MGRID + CUTOFF 300 + &END MGRID + &QS + EPS_DEFAULT 1.E-10 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL DCDFT + HARRIS_BASIS ORBITAL + &RESPONSE_SOLVER + METHOD MO_SOLVER + PRECONDITIONER FULL_SINGLE_INVERSE + &END + &XC + &XC_FUNCTIONAL PBE + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-6 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL NONE + &END + &HF + &SCREENING + EPS_SCHWARZ 1.0E-14 + &END + &END + &END XC + &POISSON + PERIODIC NONE + POISSON_SOLVER MT + &END + &END DFT + &SUBSYS + &CELL + ABC 4.0 4.0 4.0 + PERIODIC NONE + &END CELL + &COORD + N 0.400000 0.000000 0.500000 + N -0.400000 0.000000 -0.500000 + &END COORD + &KIND N + BASIS_SET ORB DZVP-GTH + BASIS_SET AUX_FIT FIT3 + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE GEO_OPT + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &GEO_OPT + MAX_ITER 1 + &END +&END + diff --git a/tests/QS/regtest-dcdft-force/N2_t14.inp b/tests/QS/regtest-dcdft-force/N2_t14.inp new file mode 100644 index 0000000000..3ea6d1b276 --- /dev/null +++ b/tests/QS/regtest-dcdft-force/N2_t14.inp @@ -0,0 +1,80 @@ +&FORCE_EVAL + METHOD Quickstep + &DFT +# &PRINT +# &DERIVATIVES +# &END +# &END + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME GTH_POTENTIALS + BASIS_SET_FILE_NAME BASIS_ADMM + &AUXILIARY_DENSITY_MATRIX_METHOD + ADMM_PURIFICATION_METHOD NONE + EXCH_CORRECTION_FUNC NONE + EXCH_SCALING_MODEL NONE + METHOD BASIS_PROJECTION + &END + &MGRID + CUTOFF 300 + &END MGRID + &QS + EPS_DEFAULT 1.E-10 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL DCDFT + HARRIS_BASIS ORBITAL + &RESPONSE_SOLVER + METHOD AO_ORTHO + PRECONDITIONER MULTI_LEVEL + MAX_ITER 2 + &END + &XC + &XC_FUNCTIONAL PBE + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-6 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL NONE + &END + &HF + &SCREENING + EPS_SCHWARZ 1.0E-14 + &END + &END + &END XC + &POISSON + PERIODIC NONE + POISSON_SOLVER MT + &END + &END DFT + &SUBSYS + &CELL + ABC 4.0 4.0 4.0 + PERIODIC NONE + &END CELL + &COORD + N 0.400000 0.000000 0.500000 + N -0.400000 0.000000 -0.500000 + &END COORD + &KIND N + BASIS_SET ORB DZVP-GTH + BASIS_SET AUX_FIT FIT3 + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE GEO_OPT + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &GEO_OPT + MAX_ITER 1 + &END +&END + diff --git a/tests/QS/regtest-dcdft-force/N2_t15.inp b/tests/QS/regtest-dcdft-force/N2_t15.inp new file mode 100644 index 0000000000..dc04a9fa2c --- /dev/null +++ b/tests/QS/regtest-dcdft-force/N2_t15.inp @@ -0,0 +1,80 @@ +&FORCE_EVAL + METHOD Quickstep + &DFT +# &PRINT +# &DERIVATIVES +# &END +# &END + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME GTH_POTENTIALS + BASIS_SET_FILE_NAME BASIS_ADMM + &AUXILIARY_DENSITY_MATRIX_METHOD + ADMM_PURIFICATION_METHOD NONE + EXCH_CORRECTION_FUNC DEFAULT + EXCH_SCALING_MODEL NONE + METHOD BASIS_PROJECTION + &END + &MGRID + CUTOFF 200 + &END MGRID + &QS + EPS_DEFAULT 1.E-10 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL DCDFT + HARRIS_BASIS ORBITAL + &RESPONSE_SOLVER + METHOD MO_SOLVER + PRECONDITIONER FULL_SINGLE_INVERSE + MAX_ITER 5 + &END + &XC + &XC_FUNCTIONAL PBE + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-5 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL NONE + &END + &HF + &SCREENING + EPS_SCHWARZ 1.0E-8 + &END + &END + &END XC + &POISSON + PERIODIC NONE + POISSON_SOLVER MT + &END + &END DFT + &SUBSYS + &CELL + ABC 4.0 4.0 4.0 + PERIODIC NONE + &END CELL + &COORD + N 0.400000 0.000000 0.500000 + N -0.400000 0.000000 -0.500000 + &END COORD + &KIND N + BASIS_SET ORB DZVP-GTH + BASIS_SET AUX_FIT FIT3 + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE GEO_OPT + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &GEO_OPT + MAX_ITER 1 + &END +&END + diff --git a/tests/QS/regtest-dcdft-force/N2_t16.inp b/tests/QS/regtest-dcdft-force/N2_t16.inp new file mode 100644 index 0000000000..6957865bb3 --- /dev/null +++ b/tests/QS/regtest-dcdft-force/N2_t16.inp @@ -0,0 +1,80 @@ +&FORCE_EVAL + METHOD Quickstep + &DFT +# &PRINT +# &DERIVATIVES +# &END +# &END + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME GTH_POTENTIALS + BASIS_SET_FILE_NAME BASIS_ADMM + &AUXILIARY_DENSITY_MATRIX_METHOD + ADMM_PURIFICATION_METHOD NONE + EXCH_CORRECTION_FUNC DEFAULT + EXCH_SCALING_MODEL NONE + METHOD BASIS_PROJECTION + &END + &MGRID + CUTOFF 200 + &END MGRID + &QS + EPS_DEFAULT 1.E-10 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL DCDFT + HARRIS_BASIS ORBITAL + &RESPONSE_SOLVER + METHOD AO_ORTHO + PRECONDITIONER MULTI_LEVEL + MAX_ITER 2 + &END + &XC + &XC_FUNCTIONAL PBE + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-5 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL NONE + &END + &HF + &SCREENING + EPS_SCHWARZ 1.0E-14 + &END + &END + &END XC + &POISSON + PERIODIC NONE + POISSON_SOLVER MT + &END + &END DFT + &SUBSYS + &CELL + ABC 4.0 4.0 4.0 + PERIODIC NONE + &END CELL + &COORD + N 0.400000 0.000000 0.500000 + N -0.400000 0.000000 -0.500000 + &END COORD + &KIND N + BASIS_SET ORB DZVP-GTH + BASIS_SET AUX_FIT FIT3 + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE GEO_OPT + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &GEO_OPT + MAX_ITER 1 + &END +&END + diff --git a/tests/QS/regtest-dcdft-force/TEST_FILES b/tests/QS/regtest-dcdft-force/TEST_FILES new file mode 100644 index 0000000000..5ba775e0fc --- /dev/null +++ b/tests/QS/regtest-dcdft-force/TEST_FILES @@ -0,0 +1,37 @@ +# 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 +# see regtest/TEST_FILES +# Density-corrected DFT Force tests +# LDA - LDA +N2_t01.inp 11 1e-12 -19.811978822953012 +# LDA - PBE +N2_t02.inp 11 1e-11 -19.918999952069047 +# PBE - LDA +N2_t03.inp 11 1e-11 -19.810879217896087 +# PBE - PBE +N2_t04.inp 11 1e-12 -19.919507647923325 +# LDA - TPSS MO solver +N2_t05.inp 11 5e-12 -19.597072615859805 +# TPSS - LDA MO solver +N2_t06.inp 11 5e-12 -19.794531050533937 +# LDA - BR89 AO solver +N2_t07.inp 11 1e-12 -19.518445253058768 +# BR89 - LDA AO solver +N2_t08.inp 11 1e-12 -19.794774602971181 +# HFX - PBE +N2_t09.inp 11 1e-12 -19.910974562359936 +# HFX - PBE AO solver +N2_t10.inp 11 1e-12 -19.910894021392888 +# PBE0 - LDA +N2_t11.inp 11 1e-12 -19.807038516382899 +# PBE0 - LDA AO solver +N2_t12.inp 11 1e-12 -19.807038021121674 +# HFX-ADMM-NONE - PBE +N2_t13.inp 11 1e-12 -19.890454327406964 +# HFX-ADMM-NONE - PBE AO solver +N2_t14.inp 11 1e-12 -19.890468780431696 +# HFX-ADMM-DEFAULT - PBE +N2_t15.inp 11 1e-12 -19.894462675351377 +# HFX-ADMM-DEFAULT - PBE AO solver +N2_t16.inp 11 1e-12 -19.894473078176457 +#EOF diff --git a/tests/QS/regtest-dcdft-stress/N2_t01.inp b/tests/QS/regtest-dcdft-stress/N2_t01.inp new file mode 100644 index 0000000000..0395de30c1 --- /dev/null +++ b/tests/QS/regtest-dcdft-stress/N2_t01.inp @@ -0,0 +1,64 @@ +&FORCE_EVAL + METHOD Quickstep + STRESS_TENSOR ANALYTICAL + &PRINT + &FORCES + &END FORCES + &STRESS_TENSOR + COMPONENTS + &END STRESS_TENSOR + &END + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 200 + &END MGRID + &QS + EPS_DEFAULT 1.E-10 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL DCDFT + HARRIS_BASIS ORBITAL + &RESPONSE_SOLVER + METHOD MO_SOLVER + PRECONDITIONER FULL_SINGLE_INVERSE + &END + &XC + &XC_FUNCTIONAL PBE + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-6 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL PADE + &END + &END XC + &END DFT + &SUBSYS + &CELL + ABC 4.0 4.0 4.0 + &END CELL + &COORD + N 0.400000 0.000000 0.500000 + N -0.400000 0.000000 -0.500000 + &END COORD + &KIND N + BASIS_SET ORB DZVP-GTH + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE CELL_OPT + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &CELL_OPT + MAX_ITER 1 + &END +&END diff --git a/tests/QS/regtest-dcdft-stress/N2_t02.inp b/tests/QS/regtest-dcdft-stress/N2_t02.inp new file mode 100644 index 0000000000..2c57e77dba --- /dev/null +++ b/tests/QS/regtest-dcdft-stress/N2_t02.inp @@ -0,0 +1,65 @@ +&FORCE_EVAL + METHOD Quickstep + STRESS_TENSOR ANALYTICAL + &PRINT + &FORCES + &END FORCES + &STRESS_TENSOR + COMPONENTS + &END STRESS_TENSOR + &END + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 200 + &END MGRID + &QS + EPS_DEFAULT 1.E-10 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL DCDFT + HARRIS_BASIS ORBITAL + &RESPONSE_SOLVER + METHOD AO_ORTHO + PRECONDITIONER MULTI_LEVEL + MAX_ITER 2 + &END + &XC + &XC_FUNCTIONAL PBE + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-6 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL PADE + &END + &END XC + &END DFT + &SUBSYS + &CELL + ABC 4.0 4.0 4.0 + &END CELL + &COORD + N 0.400000 0.000000 0.500000 + N -0.400000 0.000000 -0.500000 + &END COORD + &KIND N + BASIS_SET ORB DZVP-GTH + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE CELL_OPT + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &CELL_OPT + MAX_ITER 1 + &END +&END diff --git a/tests/QS/regtest-dcdft-stress/N2_t03.inp b/tests/QS/regtest-dcdft-stress/N2_t03.inp new file mode 100644 index 0000000000..31c13e1c6c --- /dev/null +++ b/tests/QS/regtest-dcdft-stress/N2_t03.inp @@ -0,0 +1,64 @@ +&FORCE_EVAL + METHOD Quickstep + STRESS_TENSOR ANALYTICAL + &PRINT + &FORCES + &END FORCES + &STRESS_TENSOR + COMPONENTS + &END STRESS_TENSOR + &END + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 200 + &END MGRID + &QS + EPS_DEFAULT 1.E-10 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL DCDFT + HARRIS_BASIS ORBITAL + &RESPONSE_SOLVER + METHOD MO_SOLVER + PRECONDITIONER FULL_SINGLE_INVERSE + &END + &XC + &XC_FUNCTIONAL PADE + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-6 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL PBE + &END + &END XC + &END DFT + &SUBSYS + &CELL + ABC 4.0 4.0 4.0 + &END CELL + &COORD + N 0.400000 0.000000 0.500000 + N -0.400000 0.000000 -0.500000 + &END COORD + &KIND N + BASIS_SET ORB DZVP-GTH + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE CELL_OPT + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &CELL_OPT + MAX_ITER 1 + &END +&END diff --git a/tests/QS/regtest-dcdft-stress/N2_t04.inp b/tests/QS/regtest-dcdft-stress/N2_t04.inp new file mode 100644 index 0000000000..2c57e77dba --- /dev/null +++ b/tests/QS/regtest-dcdft-stress/N2_t04.inp @@ -0,0 +1,65 @@ +&FORCE_EVAL + METHOD Quickstep + STRESS_TENSOR ANALYTICAL + &PRINT + &FORCES + &END FORCES + &STRESS_TENSOR + COMPONENTS + &END STRESS_TENSOR + &END + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 200 + &END MGRID + &QS + EPS_DEFAULT 1.E-10 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL DCDFT + HARRIS_BASIS ORBITAL + &RESPONSE_SOLVER + METHOD AO_ORTHO + PRECONDITIONER MULTI_LEVEL + MAX_ITER 2 + &END + &XC + &XC_FUNCTIONAL PBE + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-6 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL PADE + &END + &END XC + &END DFT + &SUBSYS + &CELL + ABC 4.0 4.0 4.0 + &END CELL + &COORD + N 0.400000 0.000000 0.500000 + N -0.400000 0.000000 -0.500000 + &END COORD + &KIND N + BASIS_SET ORB DZVP-GTH + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE CELL_OPT + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &CELL_OPT + MAX_ITER 1 + &END +&END diff --git a/tests/QS/regtest-dcdft-stress/N2_t05.inp b/tests/QS/regtest-dcdft-stress/N2_t05.inp new file mode 100644 index 0000000000..0e5058b362 --- /dev/null +++ b/tests/QS/regtest-dcdft-stress/N2_t05.inp @@ -0,0 +1,66 @@ +&FORCE_EVAL + METHOD Quickstep + STRESS_TENSOR ANALYTICAL + &PRINT + &FORCES + &END FORCES + &STRESS_TENSOR + COMPONENTS + &END STRESS_TENSOR + &END + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 200 + &END MGRID + &QS + EPS_DEFAULT 1.E-10 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL DCDFT + HARRIS_BASIS ORBITAL + &RESPONSE_SOLVER + METHOD MO_SOLVER + PRECONDITIONER FULL_SINGLE_INVERSE + &END + &XC + &XC_FUNCTIONAL + &MGGA_X_TPSS + &END + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-6 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL PADE + &END + &END XC + &END DFT + &SUBSYS + &CELL + ABC 4.0 4.0 4.0 + &END CELL + &COORD + N 0.400000 0.000000 0.500000 + N -0.400000 0.000000 -0.500000 + &END COORD + &KIND N + BASIS_SET ORB DZVP-GTH + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE CELL_OPT + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &CELL_OPT + MAX_ITER 1 + &END +&END diff --git a/tests/QS/regtest-dcdft-stress/N2_t06.inp b/tests/QS/regtest-dcdft-stress/N2_t06.inp new file mode 100644 index 0000000000..4547556ab1 --- /dev/null +++ b/tests/QS/regtest-dcdft-stress/N2_t06.inp @@ -0,0 +1,71 @@ +&FORCE_EVAL + METHOD Quickstep + STRESS_TENSOR ANALYTICAL + &PRINT + &FORCES + &END FORCES + &STRESS_TENSOR + COMPONENTS + &END STRESS_TENSOR + &END + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 200 + &END MGRID + &QS + EPS_DEFAULT 1.E-10 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL DCDFT + HARRIS_BASIS ORBITAL + &RESPONSE_SOLVER + METHOD MO_SOLVER + PRECONDITIONER FULL_SINGLE_INVERSE + &END + &XC + &XC_FUNCTIONAL PADE + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-6 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL + &MGGA_X_TPSS + &END + &END + &END XC + &POISSON + PERIODIC NONE + POISSON_SOLVER MT + &END + &END DFT + &SUBSYS + &CELL + ABC 4.0 4.0 4.0 + PERIODIC NONE + &END CELL + &COORD + N 0.400000 0.000000 0.500000 + N -0.400000 0.000000 -0.500000 + &END COORD + &KIND N + BASIS_SET ORB DZVP-GTH + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE CELL_OPT + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &CELL_OPT + MAX_ITER 1 + &END +&END diff --git a/tests/QS/regtest-dcdft-stress/N2_t07.inp b/tests/QS/regtest-dcdft-stress/N2_t07.inp new file mode 100644 index 0000000000..dd7cd9506e --- /dev/null +++ b/tests/QS/regtest-dcdft-stress/N2_t07.inp @@ -0,0 +1,71 @@ +&FORCE_EVAL + METHOD Quickstep + STRESS_TENSOR ANALYTICAL + &PRINT + &FORCES + &END FORCES + &STRESS_TENSOR + COMPONENTS + &END STRESS_TENSOR + &END + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 200 + &END MGRID + &QS + EPS_DEFAULT 1.E-10 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL DCDFT + HARRIS_BASIS ORBITAL + &RESPONSE_SOLVER + METHOD AO_ORTHO + PRECONDITIONER MULTI_LEVEL + &END + &XC + &XC_FUNCTIONAL + &MGGA_X_BR89 + &END + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-6 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL PADE + &END + &END XC + &POISSON + PERIODIC NONE + POISSON_SOLVER MT + &END + &END DFT + &SUBSYS + &CELL + ABC 4.0 4.0 4.0 + PERIODIC NONE + &END CELL + &COORD + N 0.400000 0.000000 0.500000 + N -0.400000 0.000000 -0.500000 + &END COORD + &KIND N + BASIS_SET ORB DZVP-GTH + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE CELL_OPT + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &CELL_OPT + MAX_ITER 1 + &END +&END diff --git a/tests/QS/regtest-dcdft-stress/N2_t08.inp b/tests/QS/regtest-dcdft-stress/N2_t08.inp new file mode 100644 index 0000000000..10cfe19ede --- /dev/null +++ b/tests/QS/regtest-dcdft-stress/N2_t08.inp @@ -0,0 +1,72 @@ +&FORCE_EVAL + METHOD Quickstep + STRESS_TENSOR ANALYTICAL + &PRINT + &FORCES + &END FORCES + &STRESS_TENSOR + COMPONENTS + &END STRESS_TENSOR + &END + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 200 + &END MGRID + &QS + EPS_DEFAULT 1.E-10 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL DCDFT + HARRIS_BASIS ORBITAL + &RESPONSE_SOLVER + METHOD AO_ORTHO + PRECONDITIONER MULTI_LEVEL + MAX_ITER 2 + &END + &XC + &XC_FUNCTIONAL PADE + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-6 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL + &MGGA_X_BR89 + &END + &END + &END XC + &POISSON + PERIODIC NONE + POISSON_SOLVER MT + &END + &END DFT + &SUBSYS + &CELL + ABC 4.0 4.0 4.0 + PERIODIC NONE + &END CELL + &COORD + N 0.400000 0.000000 0.500000 + N -0.400000 0.000000 -0.500000 + &END COORD + &KIND N + BASIS_SET ORB DZVP-GTH + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE CELL_OPT + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &CELL_OPT + MAX_ITER 1 + &END +&END diff --git a/tests/QS/regtest-dcdft-stress/N2_t09.inp b/tests/QS/regtest-dcdft-stress/N2_t09.inp new file mode 100644 index 0000000000..113696d9ce --- /dev/null +++ b/tests/QS/regtest-dcdft-stress/N2_t09.inp @@ -0,0 +1,71 @@ +&FORCE_EVAL + METHOD Quickstep + STRESS_TENSOR ANALYTICAL + &PRINT + &FORCES + &END FORCES + &STRESS_TENSOR + COMPONENTS + &END STRESS_TENSOR + &END + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 200 + REL_CUTOFF 20 + &END MGRID + &QS + EPS_DEFAULT 1.E-10 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL DCDFT + HARRIS_BASIS ORBITAL + &RESPONSE_SOLVER + METHOD MO_SOLVER + PRECONDITIONER FULL_SINGLE_INVERSE + MAX_ITER 3 + &END + &XC + &XC_FUNCTIONAL PBE + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-5 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL NONE + &END + &HF + &SCREENING + EPS_SCHWARZ 1.0E-14 + &END + &END + &END XC + &END DFT + &SUBSYS + &CELL + ABC 4.0 4.0 4.0 + &END CELL + &COORD + N 0.400000 0.000000 0.500000 + N -0.400000 0.000000 -0.500000 + &END COORD + &KIND N + BASIS_SET ORB DZVP-GTH + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE CELL_OPT + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &CELL_OPT + MAX_ITER 1 + &END +&END diff --git a/tests/QS/regtest-dcdft-stress/N2_t10.inp b/tests/QS/regtest-dcdft-stress/N2_t10.inp new file mode 100644 index 0000000000..4a41b21375 --- /dev/null +++ b/tests/QS/regtest-dcdft-stress/N2_t10.inp @@ -0,0 +1,71 @@ +&FORCE_EVAL + METHOD Quickstep + STRESS_TENSOR ANALYTICAL + &PRINT + &FORCES + &END FORCES + &STRESS_TENSOR + COMPONENTS + &END STRESS_TENSOR + &END + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 200 + REL_CUTOFF 20 + &END MGRID + &QS + EPS_DEFAULT 1.E-10 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL DCDFT + HARRIS_BASIS ORBITAL + &RESPONSE_SOLVER + METHOD AO_ORTHO + PRECONDITIONER MULTI_LEVEL + MAX_ITER 2 + &END + &XC + &XC_FUNCTIONAL PBE + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-5 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL NONE + &END + &HF + &SCREENING + EPS_SCHWARZ 1.0E-14 + &END + &END + &END XC + &END DFT + &SUBSYS + &CELL + ABC 4.0 4.0 4.0 + &END CELL + &COORD + N 0.400000 0.000000 0.500000 + N -0.400000 0.000000 -0.500000 + &END COORD + &KIND N + BASIS_SET ORB DZVP-GTH + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE CELL_OPT + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &CELL_OPT + MAX_ITER 1 + &END +&END diff --git a/tests/QS/regtest-dcdft-stress/N2_t11.inp b/tests/QS/regtest-dcdft-stress/N2_t11.inp new file mode 100644 index 0000000000..53cc33bd03 --- /dev/null +++ b/tests/QS/regtest-dcdft-stress/N2_t11.inp @@ -0,0 +1,64 @@ +&FORCE_EVAL + METHOD Quickstep + STRESS_TENSOR ANALYTICAL + &PRINT + &FORCES + &END FORCES + &STRESS_TENSOR + COMPONENTS + &END STRESS_TENSOR + &END + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 200 + &END MGRID + &QS + EPS_DEFAULT 1.E-10 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL DCDFT + HARRIS_BASIS ORBITAL + &RESPONSE_SOLVER + METHOD MO_SOLVER + PRECONDITIONER FULL_SINGLE_INVERSE + &END + &XC + &XC_FUNCTIONAL PADE + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-6 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL PBE0 + &END + &END XC + &POISSON + PERIODIC NONE + POISSON_SOLVER MT + &END + &END DFT + &SUBSYS + &CELL + ABC 4.0 4.0 4.0 + PERIODIC NONE + &END CELL + &COORD + N 0.400000 0.000000 0.500000 + N -0.400000 0.000000 -0.500000 + &END COORD + &KIND N + BASIS_SET ORB DZVP-GTH + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE ENERGY_FORCE + PRINT_LEVEL LOW +&END GLOBAL diff --git a/tests/QS/regtest-dcdft-stress/N2_t12.inp b/tests/QS/regtest-dcdft-stress/N2_t12.inp new file mode 100644 index 0000000000..27ce705c34 --- /dev/null +++ b/tests/QS/regtest-dcdft-stress/N2_t12.inp @@ -0,0 +1,65 @@ +&FORCE_EVAL + METHOD Quickstep + STRESS_TENSOR ANALYTICAL + &PRINT + &FORCES + &END FORCES + &STRESS_TENSOR + COMPONENTS + &END STRESS_TENSOR + &END + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 200 + &END MGRID + &QS + EPS_DEFAULT 1.E-10 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL DCDFT + HARRIS_BASIS ORBITAL + &RESPONSE_SOLVER + METHOD AO_ORTHO + PRECONDITIONER MULTI_LEVEL + MAX_ITER 2 + &END + &XC + &XC_FUNCTIONAL PADE + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-6 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL PBE0 + &END + &END XC + &POISSON + PERIODIC NONE + POISSON_SOLVER MT + &END + &END DFT + &SUBSYS + &CELL + ABC 4.0 4.0 4.0 + PERIODIC NONE + &END CELL + &COORD + N 0.400000 0.000000 0.500000 + N -0.400000 0.000000 -0.500000 + &END COORD + &KIND N + BASIS_SET ORB DZVP-GTH + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE ENERGY_FORCE + PRINT_LEVEL LOW +&END GLOBAL diff --git a/tests/QS/regtest-dcdft-stress/N2_t13.inp b/tests/QS/regtest-dcdft-stress/N2_t13.inp new file mode 100644 index 0000000000..da409d6e40 --- /dev/null +++ b/tests/QS/regtest-dcdft-stress/N2_t13.inp @@ -0,0 +1,83 @@ +&FORCE_EVAL + METHOD Quickstep + STRESS_TENSOR ANALYTICAL + &PRINT + &FORCES + &END FORCES + &STRESS_TENSOR + COMPONENTS + &END STRESS_TENSOR + &END + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME GTH_POTENTIALS + BASIS_SET_FILE_NAME BASIS_ADMM + &AUXILIARY_DENSITY_MATRIX_METHOD + ADMM_PURIFICATION_METHOD NONE + EXCH_CORRECTION_FUNC NONE + EXCH_SCALING_MODEL NONE + METHOD BASIS_PROJECTION + &END + &MGRID + CUTOFF 200 + REL_CUTOFF 20 + &END MGRID + &QS + EPS_DEFAULT 1.E-10 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL DCDFT + HARRIS_BASIS ORBITAL + &RESPONSE_SOLVER + METHOD MO_SOLVER + PRECONDITIONER FULL_SINGLE_INVERSE + &END + &XC + &XC_FUNCTIONAL PBE + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-6 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL NONE + &END + &HF + &SCREENING + EPS_SCHWARZ 1.0E-14 + &END + &END + &END XC + &POISSON + PERIODIC NONE + POISSON_SOLVER MT + &END + &END DFT + &SUBSYS + &CELL + ABC 4.0 4.0 4.0 + PERIODIC NONE + &END CELL + &COORD + N 0.400000 0.000000 0.500000 + N -0.400000 0.000000 -0.500000 + &END COORD + &KIND N + BASIS_SET ORB DZVP-GTH + BASIS_SET AUX_FIT FIT3 + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE CELL_OPT + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &CELL_OPT + MAX_ITER 1 + &END +&END diff --git a/tests/QS/regtest-dcdft-stress/N2_t14.inp b/tests/QS/regtest-dcdft-stress/N2_t14.inp new file mode 100644 index 0000000000..aa398047a5 --- /dev/null +++ b/tests/QS/regtest-dcdft-stress/N2_t14.inp @@ -0,0 +1,83 @@ +&FORCE_EVAL + METHOD Quickstep + STRESS_TENSOR ANALYTICAL + &PRINT + &FORCES + &END FORCES + &STRESS_TENSOR + COMPONENTS + &END STRESS_TENSOR + &END + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME GTH_POTENTIALS + BASIS_SET_FILE_NAME BASIS_ADMM + &AUXILIARY_DENSITY_MATRIX_METHOD + ADMM_PURIFICATION_METHOD NONE + EXCH_CORRECTION_FUNC NONE + EXCH_SCALING_MODEL NONE + METHOD BASIS_PROJECTION + &END + &MGRID + CUTOFF 200 + REL_CUTOFF 20 + &END MGRID + &QS + EPS_DEFAULT 1.E-10 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL DCDFT + HARRIS_BASIS ORBITAL + &RESPONSE_SOLVER + METHOD AO_ORTHO + PRECONDITIONER MULTI_LEVEL + &END + &XC + &XC_FUNCTIONAL PBE + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-6 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL NONE + &END + &HF + &SCREENING + EPS_SCHWARZ 1.0E-14 + &END + &END + &END XC + &POISSON + PERIODIC NONE + POISSON_SOLVER MT + &END + &END DFT + &SUBSYS + &CELL + ABC 4.0 4.0 4.0 + PERIODIC NONE + &END CELL + &COORD + N 0.400000 0.000000 0.500000 + N -0.400000 0.000000 -0.500000 + &END COORD + &KIND N + BASIS_SET ORB DZVP-GTH + BASIS_SET AUX_FIT FIT3 + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE CELL_OPT + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &CELL_OPT + MAX_ITER 1 + &END +&END diff --git a/tests/QS/regtest-dcdft-stress/N2_t15.inp b/tests/QS/regtest-dcdft-stress/N2_t15.inp new file mode 100644 index 0000000000..f821009cc3 --- /dev/null +++ b/tests/QS/regtest-dcdft-stress/N2_t15.inp @@ -0,0 +1,83 @@ +&FORCE_EVAL + METHOD Quickstep + STRESS_TENSOR ANALYTICAL + &PRINT + &FORCES + &END FORCES + &STRESS_TENSOR + COMPONENTS + &END STRESS_TENSOR + &END + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME GTH_POTENTIALS + BASIS_SET_FILE_NAME BASIS_ADMM + &AUXILIARY_DENSITY_MATRIX_METHOD + ADMM_PURIFICATION_METHOD NONE + EXCH_CORRECTION_FUNC DEFAULT + EXCH_SCALING_MODEL NONE + METHOD BASIS_PROJECTION + &END + &MGRID + CUTOFF 200 + REL_CUTOFF 20 + &END MGRID + &QS + EPS_DEFAULT 1.E-10 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL DCDFT + HARRIS_BASIS ORBITAL + &RESPONSE_SOLVER + METHOD MO_SOLVER + PRECONDITIONER FULL_SINGLE_INVERSE + &END + &XC + &XC_FUNCTIONAL PBE + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-6 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL NONE + &END + &HF + &SCREENING + EPS_SCHWARZ 1.0E-14 + &END + &END + &END XC + &POISSON + PERIODIC NONE + POISSON_SOLVER MT + &END + &END DFT + &SUBSYS + &CELL + ABC 4.0 4.0 4.0 + PERIODIC NONE + &END CELL + &COORD + N 0.400000 0.000000 0.500000 + N -0.400000 0.000000 -0.500000 + &END COORD + &KIND N + BASIS_SET ORB DZVP-GTH + BASIS_SET AUX_FIT FIT3 + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE CELL_OPT + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &CELL_OPT + MAX_ITER 1 + &END +&END diff --git a/tests/QS/regtest-dcdft-stress/N2_t16.inp b/tests/QS/regtest-dcdft-stress/N2_t16.inp new file mode 100644 index 0000000000..39cf011d78 --- /dev/null +++ b/tests/QS/regtest-dcdft-stress/N2_t16.inp @@ -0,0 +1,83 @@ +&FORCE_EVAL + METHOD Quickstep + STRESS_TENSOR ANALYTICAL + &PRINT + &FORCES + &END FORCES + &STRESS_TENSOR + COMPONENTS + &END STRESS_TENSOR + &END + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME GTH_POTENTIALS + BASIS_SET_FILE_NAME BASIS_ADMM + &AUXILIARY_DENSITY_MATRIX_METHOD + ADMM_PURIFICATION_METHOD NONE + EXCH_CORRECTION_FUNC DEFAULT + EXCH_SCALING_MODEL NONE + METHOD BASIS_PROJECTION + &END + &MGRID + CUTOFF 200 + REL_CUTOFF 20 + &END MGRID + &QS + EPS_DEFAULT 1.E-10 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL DCDFT + HARRIS_BASIS ORBITAL + &RESPONSE_SOLVER + METHOD AO_ORTHO + PRECONDITIONER MULTI_LEVEL + &END + &XC + &XC_FUNCTIONAL PBE + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-6 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL NONE + &END + &HF + &SCREENING + EPS_SCHWARZ 1.0E-14 + &END + &END + &END XC + &POISSON + PERIODIC NONE + POISSON_SOLVER MT + &END + &END DFT + &SUBSYS + &CELL + ABC 4.0 4.0 4.0 + PERIODIC NONE + &END CELL + &COORD + N 0.400000 0.000000 0.500000 + N -0.400000 0.000000 -0.500000 + &END COORD + &KIND N + BASIS_SET ORB DZVP-GTH + BASIS_SET AUX_FIT FIT3 + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE CELL_OPT + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &CELL_OPT + MAX_ITER 1 + &END +&END diff --git a/tests/QS/regtest-dcdft-stress/TEST_FILES b/tests/QS/regtest-dcdft-stress/TEST_FILES new file mode 100644 index 0000000000..8e1aa2d219 --- /dev/null +++ b/tests/QS/regtest-dcdft-stress/TEST_FILES @@ -0,0 +1,37 @@ +# 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 +# see regtest/TEST_FILES +# Density-corrected DFT Stress tensor tests +# LDA - PBE MO solver +N2_t01.inp 31 1e-05 3.27841013742E+01 +# LDA - PBE AO solver +N2_t02.inp 31 1e-05 3.27915201801E+01 +# PBE - LDA MO solver +N2_t03.inp 31 1e-05 3.27261434593E+01 +# PBE - LDA AO solver +N2_t04.inp 31 1e-05 3.27915201801E+01 +# LDA - TPSS MO solver +N2_t05.inp 31 1e-05 2.76973145058E+01 +# TPSS - LDA MO solver +N2_t06.inp 31 1e-05 -5.02530331310E+02 +# LDA - BR89 AO solver +N2_t07.inp 31 1e-05 -5.03370759578E+02 +# BR89 - LDA AO solver +N2_t08.inp 31 1e-05 -5.02524932385E+02 +# HFX - PBE MO solver +N2_t09.inp 31 1e-05 6.87213341923E+01 +# HFX - PBE AO solver +N2_t10.inp 31 1e-05 6.90331816280E+01 +# PBE0 - PADE MO solver +N2_t11.inp 31 1e-05 -4.60746606987E+02 +# PBE0 - PADE AO solver +N2_t12.inp 31 1e-05 -4.60764171148E+02 +# HFX-ADMM-NONE - PBE MO solver +N2_t13.inp 31 1e-05 -4.97852690104E+02 +# HFX-ADMM-NONE - PBE AO solver +N2_t14.inp 31 1e-05 -4.97852717316E+02 +# HFX-ADMM-DEFAULT - PBE MO solver +N2_t15.inp 31 1e-05 -4.98121970425E+02 +# HFX-ADMM-DEFAULT - PBE AO solver +N2_t16.inp 31 1e-05 -4.98121969330E+02 +#EOF diff --git a/tests/QS/regtest-ec-force-meta/TEST_FILES b/tests/QS/regtest-ec-force-meta/TEST_FILES deleted file mode 100644 index dd31b07ce8..0000000000 --- a/tests/QS/regtest-ec-force-meta/TEST_FILES +++ /dev/null @@ -1,8 +0,0 @@ -# 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 -# see regtest/TEST_FILES -N2_tpss_pade.inp 1 1e-12 -19.52888605962090 -N2_pade_tpss.inp 1 1e-12 -19.79534215334482 -N2_br89_pade.inp 1 1e-12 -19.51922270784055 -N2_pade_br89.inp 1 1e-12 -19.79514891375846 -#EOF diff --git a/tests/QS/regtest-ec-force/N2_t03.inp b/tests/QS/regtest-ec-force/N2_t03.inp index 0e2c8208e1..e20d34ec95 100644 --- a/tests/QS/regtest-ec-force/N2_t03.inp +++ b/tests/QS/regtest-ec-force/N2_t03.inp @@ -8,7 +8,7 @@ BASIS_SET_FILE_NAME BASIS_SET POTENTIAL_FILE_NAME GTH_POTENTIALS &MGRID - CUTOFF 300 + CUTOFF 200 &END MGRID &QS EPS_DEFAULT 1.E-12 @@ -21,7 +21,7 @@ PRECONDITIONER FULL_SINGLE_INVERSE &END &XC - &XC_FUNCTIONAL PBE + &XC_FUNCTIONAL PADE &END &END XC &END ENERGY_CORRECTION @@ -30,13 +30,13 @@ SCF_GUESS ATOMIC &END &XC - &XC_FUNCTIONAL PADE + &XC_FUNCTIONAL PBE &END &END XC &END DFT &SUBSYS &CELL - ABC 6.0 6.0 6.0 + ABC 5.0 5.0 5.0 &END CELL &COORD N 0.400000 0.000000 0.500000 diff --git a/tests/QS/regtest-ec-force/TEST_FILES b/tests/QS/regtest-ec-force/TEST_FILES index a63784aa1d..5809afc8b5 100644 --- a/tests/QS/regtest-ec-force/TEST_FILES +++ b/tests/QS/regtest-ec-force/TEST_FILES @@ -1,19 +1,35 @@ # 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 # see regtest/TEST_FILES -N2_t01.inp 1 1e-12 -19.81197882293305 -N2_t02.inp 1 5e-12 -19.81193992209264 -N2_t03.inp 1 5e-12 -19.81193992209264 -N2_t04.inp 1 1e-12 -19.91950764923067 -N2_t05.inp 1 5e-11 -19.45601327358312 -N2_t06.inp 1 5e-11 -19.89228099201857 -N2_t07.inp 1 1e-11 -19.44508292698253 -N2_t08.inp 1 5e-11 -19.54242231073764 -N2_t08b.inp 1 5e-11 -19.54242225427242 -N2_t09.inp 1 5e-12 -19.52266035514095 -N2_t10.inp 1 1e-12 -19.10767135870040 -N2_t11.inp 1 1e-12 -19.87201569403329 -N2_t12.inp 1 1e-12 -19.87215956274964 -N2_t13.inp 1 1e-12 -19.87400375270014 -N2_t14.inp 1 5e-10 -19.87377227302795 +# Harris functional energy correction force tests +# PADE - PADE MO solver +N2_t01.inp 11 1e-12 -19.811978822934346 +# PADE - PBE +N2_t02.inp 11 5e-12 -19.920157485812176 +# PBE - PADE +N2_t03.inp 11 5e-11 -19.812276059625223 +# PBE - PBE +N2_t04.inp 11 1e-12 -19.919507649230656 +# NONE - PBE +N2_t05.inp 11 5e-11 -19.921963323203908 +# PBE0 - PADE +N2_t06.inp 11 5e-11 -19.812023408523121 +# PADE - PADE w/ Harris basis +N2_t07.inp 11 1e-11 -19.971002008270300 +# PBE - PBE w/ Harris basis +N2_t08.inp 11 5e-11 -20.085007089026508 +# PBE0 - PBE w/ Harris basis AO solver +N2_t08b.inp 11 5e-11 -20.085007124872416 +# PBE0 - PBE w/ Harris basis +N2_t09.inp 11 5e-12 -20.088188602650991 +# HFX - PBE w/ Harris basis +N2_t10.inp 11 1e-12 -20.097436479043111 +# HFX-ADMM-NONE - PBE +N2_t11.inp 11 1e-12 -19.903215469949565 +# PBE0-ADMM-NONE - PBE w/ Harris basis +N2_t12.inp 11 1e-12 -19.927659891608339 +# PBE0-ADMM-DEFAULT - PBE +N2_t13.inp 11 1e-12 -19.903137808157719 +# PBE0-ADMM-DEFAULT - PBE w/ Harris basis +N2_t14.inp 11 5e-10 -19.927232210625430 #EOF diff --git a/tests/QS/regtest-ec-force-meta/N2_br89_pade.inp b/tests/QS/regtest-ec-meta/N2_br89_pade.inp similarity index 92% rename from tests/QS/regtest-ec-force-meta/N2_br89_pade.inp rename to tests/QS/regtest-ec-meta/N2_br89_pade.inp index 45b0bd2eaa..d192b508b9 100644 --- a/tests/QS/regtest-ec-force-meta/N2_br89_pade.inp +++ b/tests/QS/regtest-ec-meta/N2_br89_pade.inp @@ -26,10 +26,10 @@ SCF_GUESS ATOMIC &END &XC - &XC_FUNCTIONAL - &MGGA_X_BR89 - &END - &END + &XC_FUNCTIONAL + &MGGA_X_BR89 + &END + &END &END XC &POISSON PERIODIC NONE @@ -58,7 +58,7 @@ &END GLOBAL &MOTION &GEO_OPT - MAX_ITER 1 + MAX_ITER 1 &END &END diff --git a/tests/QS/regtest-ec-meta/N2_br89_pade_stress.inp b/tests/QS/regtest-ec-meta/N2_br89_pade_stress.inp new file mode 100644 index 0000000000..873a6b9778 --- /dev/null +++ b/tests/QS/regtest-ec-meta/N2_br89_pade_stress.inp @@ -0,0 +1,71 @@ +&FORCE_EVAL + METHOD Quickstep + STRESS_TENSOR ANALYTICAL + &PRINT + &FORCES + &END FORCES + &STRESS_TENSOR + &END STRESS_TENSOR + &END + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 200 + &END MGRID + &QS + EPS_DEFAULT 1.E-10 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL HARRIS + HARRIS_BASIS ORBITAL + &RESPONSE_SOLVER + METHOD MO_SOLVER + PRECONDITIONER FULL_SINGLE_INVERSE + &END + &XC + &XC_FUNCTIONAL PADE + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-7 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL + &MGGA_X_BR89 + &END + &END + &END XC + &POISSON + PERIODIC NONE + POISSON_SOLVER MT + &END + &END DFT + &SUBSYS + &CELL + ABC 4.0 4.0 4.0 + PERIODIC NONE + &END CELL + &COORD + N 0.400000 0.000000 0.500000 + N -0.400000 0.000000 -0.500000 + &END COORD + &KIND N + BASIS_SET ORB DZVP-GTH + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2_br89_pade + RUN_TYPE CELL_OPT + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &CELL_OPT + MAX_ITER 1 + &END +&END + diff --git a/tests/QS/regtest-ec-force-meta/N2_pade_br89.inp b/tests/QS/regtest-ec-meta/N2_pade_br89.inp similarity index 91% rename from tests/QS/regtest-ec-force-meta/N2_pade_br89.inp rename to tests/QS/regtest-ec-meta/N2_pade_br89.inp index 5c07c195c3..a85d33ac04 100644 --- a/tests/QS/regtest-ec-force-meta/N2_pade_br89.inp +++ b/tests/QS/regtest-ec-meta/N2_pade_br89.inp @@ -13,8 +13,9 @@ ENERGY_FUNCTIONAL HARRIS HARRIS_BASIS ORBITAL &RESPONSE_SOLVER - METHOD MO_SOLVER - PRECONDITIONER FULL_SINGLE_INVERSE + METHOD AO_ORTHO + PRECONDITIONER MULTI_LEVEL + MAX_ITER 2 &END &XC &XC_FUNCTIONAL @@ -58,7 +59,7 @@ &END GLOBAL &MOTION &GEO_OPT - MAX_ITER 1 + MAX_ITER 1 &END &END diff --git a/tests/QS/regtest-ec-meta/N2_pade_br89_stress.inp b/tests/QS/regtest-ec-meta/N2_pade_br89_stress.inp new file mode 100644 index 0000000000..27c9013ab6 --- /dev/null +++ b/tests/QS/regtest-ec-meta/N2_pade_br89_stress.inp @@ -0,0 +1,72 @@ +&FORCE_EVAL + METHOD Quickstep + STRESS_TENSOR ANALYTICAL + &PRINT + &FORCES + &END FORCES + &STRESS_TENSOR + &END STRESS_TENSOR + &END + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 200 + &END MGRID + &QS + EPS_DEFAULT 1.E-10 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL HARRIS + HARRIS_BASIS ORBITAL + &RESPONSE_SOLVER + METHOD AO_ORTHO + PRECONDITIONER MULTI_LEVEL + MAX_ITER 2 + &END + &XC + &XC_FUNCTIONAL + &MGGA_X_BR89 + &END + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-7 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL PADE + &END + &END XC + &POISSON + PERIODIC NONE + POISSON_SOLVER MT + &END + &END DFT + &SUBSYS + &CELL + ABC 4.0 4.0 4.0 + PERIODIC NONE + &END CELL + &COORD + N 0.400000 0.000000 0.500000 + N -0.400000 0.000000 -0.500000 + &END COORD + &KIND N + BASIS_SET ORB DZVP-GTH + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2_pade_br89 + RUN_TYPE CELL_OPT + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &CELL_OPT + MAX_ITER 1 + &END +&END + diff --git a/tests/QS/regtest-ec-force-meta/N2_pade_tpss.inp b/tests/QS/regtest-ec-meta/N2_pade_tpss.inp similarity index 98% rename from tests/QS/regtest-ec-force-meta/N2_pade_tpss.inp rename to tests/QS/regtest-ec-meta/N2_pade_tpss.inp index cc2e9c21c7..acf5bb133b 100644 --- a/tests/QS/regtest-ec-force-meta/N2_pade_tpss.inp +++ b/tests/QS/regtest-ec-meta/N2_pade_tpss.inp @@ -58,7 +58,7 @@ &END GLOBAL &MOTION &GEO_OPT - MAX_ITER 1 + MAX_ITER 1 &END &END diff --git a/tests/QS/regtest-ec-meta/N2_pade_tpss_stress.inp b/tests/QS/regtest-ec-meta/N2_pade_tpss_stress.inp new file mode 100644 index 0000000000..d27b19420f --- /dev/null +++ b/tests/QS/regtest-ec-meta/N2_pade_tpss_stress.inp @@ -0,0 +1,71 @@ +&FORCE_EVAL + METHOD Quickstep + STRESS_TENSOR ANALYTICAL + &PRINT + &FORCES + &END FORCES + &STRESS_TENSOR + &END STRESS_TENSOR + &END + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 200 + &END MGRID + &QS + EPS_DEFAULT 1.E-10 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL HARRIS + HARRIS_BASIS ORBITAL + &RESPONSE_SOLVER + METHOD MO_SOLVER + PRECONDITIONER FULL_SINGLE_INVERSE + &END + &XC + &XC_FUNCTIONAL + &MGGA_X_TPSS + &END + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-6 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL PADE + &END + &END XC + &POISSON + PERIODIC NONE + POISSON_SOLVER MT + &END + &END DFT + &SUBSYS + &CELL + ABC 4.0 4.0 4.0 + PERIODIC NONE + &END CELL + &COORD + N 0.400000 0.000000 0.500000 + N -0.400000 0.000000 -0.500000 + &END COORD + &KIND N + BASIS_SET ORB DZVP-GTH + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2_pade_tpss + RUN_TYPE CELL_OPT + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &CELL_OPT + MAX_ITER 1 + &END +&END + diff --git a/tests/QS/regtest-ec-force-meta/N2_tpss_pade.inp b/tests/QS/regtest-ec-meta/N2_tpss_pade.inp similarity index 92% rename from tests/QS/regtest-ec-force-meta/N2_tpss_pade.inp rename to tests/QS/regtest-ec-meta/N2_tpss_pade.inp index 77f08372d6..5d2da8b0f7 100644 --- a/tests/QS/regtest-ec-force-meta/N2_tpss_pade.inp +++ b/tests/QS/regtest-ec-meta/N2_tpss_pade.inp @@ -26,10 +26,10 @@ SCF_GUESS ATOMIC &END &XC - &XC_FUNCTIONAL - &MGGA_X_TPSS - &END - &END + &XC_FUNCTIONAL + &MGGA_X_TPSS + &END + &END &END XC &POISSON PERIODIC NONE @@ -58,7 +58,7 @@ &END GLOBAL &MOTION &GEO_OPT - MAX_ITER 1 + MAX_ITER 1 &END &END diff --git a/tests/QS/regtest-ec-meta/N2_tpss_pade_stress.inp b/tests/QS/regtest-ec-meta/N2_tpss_pade_stress.inp new file mode 100644 index 0000000000..b1b0ad1ad5 --- /dev/null +++ b/tests/QS/regtest-ec-meta/N2_tpss_pade_stress.inp @@ -0,0 +1,71 @@ +&FORCE_EVAL + METHOD Quickstep + STRESS_TENSOR ANALYTICAL + &PRINT + &FORCES + &END FORCES + &STRESS_TENSOR + &END STRESS_TENSOR + &END + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 200 + &END MGRID + &QS + EPS_DEFAULT 1.E-10 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL HARRIS + HARRIS_BASIS ORBITAL + &RESPONSE_SOLVER + METHOD MO_SOLVER + PRECONDITIONER FULL_SINGLE_INVERSE + &END + &XC + &XC_FUNCTIONAL PADE + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-7 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL + &MGGA_X_TPSS + &END + &END + &END XC + &POISSON + PERIODIC NONE + POISSON_SOLVER MT + &END + &END DFT + &SUBSYS + &CELL + ABC 4.0 4.0 4.0 + PERIODIC NONE + &END CELL + &COORD + N 0.400000 0.000000 0.500000 + N -0.400000 0.000000 -0.500000 + &END COORD + &KIND N + BASIS_SET ORB DZVP-GTH + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2_tpss_pade + RUN_TYPE CELL_OPT + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &CELL_OPT + MAX_ITER 1 + &END +&END + diff --git a/tests/QS/regtest-ec-meta/TEST_FILES b/tests/QS/regtest-ec-meta/TEST_FILES new file mode 100644 index 0000000000..d0501ccde7 --- /dev/null +++ b/tests/QS/regtest-ec-meta/TEST_FILES @@ -0,0 +1,20 @@ +# 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 +# see regtest/TEST_FILES +# FORCE: PADE - TPSS MO solver +N2_pade_tpss.inp 11 1e-12 -19.529546812850718 +# FORCE: TPSS - PADE MO solver +N2_tpss_pade.inp 11 1e-12 -19.795492393131141 +# FORCE: PADE - BR89 AO solver +N2_pade_br89.inp 11 1e-12 -19.519044036224891 +# FORCE: BR89 - PADE AO solver +N2_br89_pade.inp 11 1e-12 -19.795370307536491 +# STRESS: PADE - TPSS MO solver +N2_pade_tpss_stress.inp 31 1e-05 -5.01969865527E+02 +# STRESS: TPSS - PADE MO solver +N2_tpss_pade_stress.inp 31 1e-05 -5.02652401254E+02 +# STRESS: PADE - BR89 AO solver +N2_pade_br89_stress.inp 31 1e-05 -5.03273252086E+02 +# STRESS: BR89 - PADE AO solver +N2_br89_pade_stress.inp 31 1e-05 -5.02543292201E+02 +#EOF diff --git a/tests/QS/regtest-ec-stress/N2_t01.inp b/tests/QS/regtest-ec-stress/N2_t01.inp new file mode 100644 index 0000000000..c8a8c78994 --- /dev/null +++ b/tests/QS/regtest-ec-stress/N2_t01.inp @@ -0,0 +1,63 @@ +&FORCE_EVAL + METHOD Quickstep + STRESS_TENSOR ANALYTICAL + &PRINT + &FORCES + &END FORCES + &STRESS_TENSOR + &END STRESS_TENSOR + &END + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 200 + &END MGRID + &QS + EPS_DEFAULT 1.E-10 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL HARRIS + HARRIS_BASIS ORBITAL + &RESPONSE_SOLVER + METHOD MO_SOLVER + PRECONDITIONER FULL_SINGLE_INVERSE + &END + &XC + &XC_FUNCTIONAL PADE + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-6 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL PADE + &END + &END XC + &END DFT + &SUBSYS + &CELL + ABC 4.0 4.0 4.0 + &END CELL + &COORD + N 0.400000 0.000000 0.500000 + N -0.400000 0.000000 -0.500000 + &END COORD + &KIND N + BASIS_SET ORB DZVP-GTH + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE CELL_OPT + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &CELL_OPT + MAX_ITER 1 + &END +&END diff --git a/tests/QS/regtest-ec-stress/N2_t02.inp b/tests/QS/regtest-ec-stress/N2_t02.inp new file mode 100644 index 0000000000..4858436541 --- /dev/null +++ b/tests/QS/regtest-ec-stress/N2_t02.inp @@ -0,0 +1,63 @@ +&FORCE_EVAL + METHOD Quickstep + STRESS_TENSOR ANALYTICAL + &PRINT + &FORCES + &END FORCES + &STRESS_TENSOR + &END STRESS_TENSOR + &END + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 200 + &END MGRID + &QS + EPS_DEFAULT 1.E-10 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL HARRIS + HARRIS_BASIS ORBITAL + &RESPONSE_SOLVER + METHOD MO_SOLVER + PRECONDITIONER FULL_SINGLE_INVERSE + &END + &XC + &XC_FUNCTIONAL PBE + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-6 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL PADE + &END + &END XC + &END DFT + &SUBSYS + &CELL + ABC 4.0 4.0 4.0 + &END CELL + &COORD + N 0.400000 0.000000 0.500000 + N -0.400000 0.000000 -0.500000 + &END COORD + &KIND N + BASIS_SET ORB DZVP-GTH + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE CELL_OPT + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &CELL_OPT + MAX_ITER 1 + &END +&END diff --git a/tests/QS/regtest-ec-stress/N2_t03.inp b/tests/QS/regtest-ec-stress/N2_t03.inp new file mode 100644 index 0000000000..11feecdb64 --- /dev/null +++ b/tests/QS/regtest-ec-stress/N2_t03.inp @@ -0,0 +1,63 @@ +&FORCE_EVAL + METHOD Quickstep + STRESS_TENSOR ANALYTICAL + &PRINT + &FORCES + &END FORCES + &STRESS_TENSOR + &END STRESS_TENSOR + &END + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 200 + &END MGRID + &QS + EPS_DEFAULT 1.E-10 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL HARRIS + HARRIS_BASIS ORBITAL + &RESPONSE_SOLVER + METHOD MO_SOLVER + PRECONDITIONER FULL_SINGLE_INVERSE + &END + &XC + &XC_FUNCTIONAL PADE + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-6 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL PBE + &END + &END XC + &END DFT + &SUBSYS + &CELL + ABC 4.0 4.0 4.0 + &END CELL + &COORD + N 0.400000 0.000000 0.500000 + N -0.400000 0.000000 -0.500000 + &END COORD + &KIND N + BASIS_SET ORB DZVP-GTH + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE CELL_OPT + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &CELL_OPT + MAX_ITER 1 + &END +&END diff --git a/tests/QS/regtest-ec-stress/N2_t04.inp b/tests/QS/regtest-ec-stress/N2_t04.inp new file mode 100644 index 0000000000..790ad56927 --- /dev/null +++ b/tests/QS/regtest-ec-stress/N2_t04.inp @@ -0,0 +1,63 @@ +&FORCE_EVAL + METHOD Quickstep + STRESS_TENSOR ANALYTICAL + &PRINT + &FORCES + &END FORCES + &STRESS_TENSOR + &END STRESS_TENSOR + &END + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 200 + &END MGRID + &QS + EPS_DEFAULT 1.E-10 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL HARRIS + HARRIS_BASIS ORBITAL + &RESPONSE_SOLVER + METHOD MO_SOLVER + PRECONDITIONER FULL_SINGLE_INVERSE + &END + &XC + &XC_FUNCTIONAL PBE + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-6 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL PBE + &END + &END XC + &END DFT + &SUBSYS + &CELL + ABC 4.0 4.0 4.0 + &END CELL + &COORD + N 0.400000 0.000000 0.500000 + N -0.400000 0.000000 -0.500000 + &END COORD + &KIND N + BASIS_SET ORB DZVP-GTH + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE CELL_OPT + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &CELL_OPT + MAX_ITER 1 + &END +&END diff --git a/tests/QS/regtest-ec-stress/N2_t05.inp b/tests/QS/regtest-ec-stress/N2_t05.inp new file mode 100644 index 0000000000..114576a205 --- /dev/null +++ b/tests/QS/regtest-ec-stress/N2_t05.inp @@ -0,0 +1,74 @@ +&FORCE_EVAL + METHOD Quickstep + STRESS_TENSOR ANALYTICAL + &PRINT + &FORCES + &END FORCES + &STRESS_TENSOR + COMPONENTS + &END STRESS_TENSOR + &END + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 200 + &END MGRID + &QS + EPS_DEFAULT 1.E-10 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL HARRIS + HARRIS_BASIS ORBITAL + &RESPONSE_SOLVER + METHOD MO_SOLVER + PRECONDITIONER FULL_SINGLE_INVERSE + &END + &XC + &XC_FUNCTIONAL PBE + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-6 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL NONE + &END + &HF + &SCREENING + EPS_SCHWARZ 1.0E-14 + &END + &END + &END XC + &POISSON + PERIODIC NONE + POISSON_SOLVER MT + &END + &END DFT + &SUBSYS + &CELL + ABC 4.0 4.0 4.0 + PERIODIC NONE + &END CELL + &COORD + N 0.400000 0.000000 0.500000 + N -0.400000 0.000000 -0.500000 + &END COORD + &KIND N + BASIS_SET ORB DZVP-GTH + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE CELL_OPT + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &CELL_OPT + MAX_ITER 1 + &END +&END diff --git a/tests/QS/regtest-ec-stress/N2_t06.inp b/tests/QS/regtest-ec-stress/N2_t06.inp new file mode 100644 index 0000000000..94f1789aa4 --- /dev/null +++ b/tests/QS/regtest-ec-stress/N2_t06.inp @@ -0,0 +1,68 @@ +&FORCE_EVAL + METHOD Quickstep + STRESS_TENSOR ANALYTICAL + &PRINT + &FORCES + &END FORCES + &STRESS_TENSOR + &END STRESS_TENSOR + &END + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 200 + &END MGRID + &QS + EPS_DEFAULT 1.E-10 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL HARRIS + HARRIS_BASIS ORBITAL + &RESPONSE_SOLVER + METHOD MO_SOLVER + PRECONDITIONER FULL_SINGLE_INVERSE + &END + &XC + &XC_FUNCTIONAL PADE + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-6 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL PBE0 + &END + &END XC + &POISSON + PERIODIC NONE + POISSON_SOLVER MT + &END + &END DFT + &SUBSYS + &CELL + ABC 4.0 4.0 4.0 + PERIODIC NONE + &END CELL + &COORD + N 0.400000 0.000000 0.500000 + N -0.400000 0.000000 -0.500000 + &END COORD + &KIND N + BASIS_SET ORB DZVP-GTH + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE CELL_OPT + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &CELL_OPT + MAX_ITER 1 + &END +&END diff --git a/tests/QS/regtest-ec-stress/N2_t07.inp b/tests/QS/regtest-ec-stress/N2_t07.inp new file mode 100644 index 0000000000..ef849947b9 --- /dev/null +++ b/tests/QS/regtest-ec-stress/N2_t07.inp @@ -0,0 +1,71 @@ +&FORCE_EVAL + METHOD Quickstep + STRESS_TENSOR ANALYTICAL + &PRINT + &FORCES + &END FORCES + &STRESS_TENSOR + COMPONENTS + &END STRESS_TENSOR + &END + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 200 + &END MGRID + &QS + EPS_DEFAULT 1.E-10 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL HARRIS + HARRIS_BASIS HARRIS + &RESPONSE_SOLVER + METHOD MO_SOLVER + PRECONDITIONER FULL_SINGLE_INVERSE + &END + &XC + &XC_FUNCTIONAL PADE + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-8 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL PADE + &END + &END XC + &POISSON + PERIODIC NONE + POISSON_SOLVER MT + &END + &END DFT + &SUBSYS + &CELL + ABC 4.0 4.0 4.0 + PERIODIC NONE + &END CELL + &COORD + N 0.400000 0.000000 0.500000 + N -0.400000 0.000000 -0.500000 + &END COORD + &KIND N + BASIS_SET ORB SZV-GTH + BASIS_SET HARRIS DZVP-GTH + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE CELL_OPT + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &CELL_OPT + MAX_ITER 2 + &END +&END + diff --git a/tests/QS/regtest-ec-stress/N2_t08.inp b/tests/QS/regtest-ec-stress/N2_t08.inp new file mode 100644 index 0000000000..835e30b8a1 --- /dev/null +++ b/tests/QS/regtest-ec-stress/N2_t08.inp @@ -0,0 +1,69 @@ +&FORCE_EVAL + METHOD Quickstep + STRESS_TENSOR ANALYTICAL + &PRINT + &FORCES + &END FORCES + &STRESS_TENSOR + &END STRESS_TENSOR + &END + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 200 + &END MGRID + &QS + EPS_DEFAULT 1.E-10 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL HARRIS + HARRIS_BASIS HARRIS + &RESPONSE_SOLVER + METHOD MO_SOLVER + PRECONDITIONER FULL_SINGLE_INVERSE + &END + &XC + &XC_FUNCTIONAL PBE + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-6 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL PBE + &END + &END XC + &POISSON + PERIODIC NONE + POISSON_SOLVER MT + &END + &END DFT + &SUBSYS + &CELL + ABC 4.0 4.0 4.0 + PERIODIC NONE + &END CELL + &COORD + N 0.400000 0.000000 0.500000 + N -0.400000 0.000000 -0.500000 + &END COORD + &KIND N + BASIS_SET ORB SZV-GTH + BASIS_SET HARRIS DZVP-GTH + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE CELL_OPT + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &CELL_OPT + MAX_ITER 1 + &END +&END diff --git a/tests/QS/regtest-ec-stress/N2_t08b.inp b/tests/QS/regtest-ec-stress/N2_t08b.inp new file mode 100644 index 0000000000..ef158fbb9e --- /dev/null +++ b/tests/QS/regtest-ec-stress/N2_t08b.inp @@ -0,0 +1,65 @@ +&FORCE_EVAL + METHOD Quickstep + STRESS_TENSOR ANALYTICAL + &PRINT + &FORCES + &END FORCES + &STRESS_TENSOR + &END STRESS_TENSOR + &END + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 200 + &END MGRID + &QS + EPS_DEFAULT 1.E-10 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL HARRIS + HARRIS_BASIS HARRIS + &XC + &XC_FUNCTIONAL PBE + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-7 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL PBE + &END + &END XC + &POISSON + PERIODIC NONE + POISSON_SOLVER MT + &END + &END DFT + &SUBSYS + &CELL + ABC 4.0 4.0 4.0 + PERIODIC NONE + &END CELL + &COORD + N 0.400000 0.000000 0.500000 + N -0.400000 0.000000 -0.500000 + &END COORD + &KIND N + BASIS_SET ORB SZV-GTH + BASIS_SET HARRIS DZVP-GTH + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE CELL_OPT + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &CELL_OPT + MAX_ITER 2 + &END +&END diff --git a/tests/QS/regtest-ec-stress/N2_t09.inp b/tests/QS/regtest-ec-stress/N2_t09.inp new file mode 100644 index 0000000000..53866edbf2 --- /dev/null +++ b/tests/QS/regtest-ec-stress/N2_t09.inp @@ -0,0 +1,70 @@ +&FORCE_EVAL + METHOD Quickstep + STRESS_TENSOR ANALYTICAL + &PRINT + &FORCES + &END FORCES + &STRESS_TENSOR + COMPONENTS + &END STRESS_TENSOR + &END + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 200 + &END MGRID + &QS + EPS_DEFAULT 1.E-10 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL HARRIS + HARRIS_BASIS HARRIS + &RESPONSE_SOLVER + METHOD MO_SOLVER + PRECONDITIONER FULL_SINGLE_INVERSE + &END + &XC + &XC_FUNCTIONAL PBE + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-7 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL PBE0 + &END + &END XC + &POISSON + PERIODIC NONE + POISSON_SOLVER MT + &END + &END DFT + &SUBSYS + &CELL + ABC 4.0 4.0 4.0 + PERIODIC NONE + &END CELL + &COORD + N 0.400000 0.000000 0.511111 + N -0.400000 0.000000 -0.500000 + &END COORD + &KIND N + BASIS_SET ORB SZV-GTH + BASIS_SET HARRIS DZVP-GTH + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE CELL_OPT + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &CELL_OPT + MAX_ITER 2 + &END +&END diff --git a/tests/QS/regtest-ec-stress/N2_t10.inp b/tests/QS/regtest-ec-stress/N2_t10.inp new file mode 100644 index 0000000000..af9d13e1a9 --- /dev/null +++ b/tests/QS/regtest-ec-stress/N2_t10.inp @@ -0,0 +1,76 @@ +&FORCE_EVAL + METHOD Quickstep + STRESS_TENSOR ANALYTICAL + &PRINT + &FORCES + &END FORCES + &STRESS_TENSOR + COMPONENTS + &END STRESS_TENSOR + &END + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 400 + &END MGRID + &QS + EPS_DEFAULT 1.E-14 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL HARRIS + HARRIS_BASIS HARRIS + &RESPONSE_SOLVER + METHOD MO_SOLVER + PRECONDITIONER FULL_SINGLE_INVERSE + &END + &XC + &XC_FUNCTIONAL PBE + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-8 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL NONE + &END + &HF + &SCREENING + EPS_SCHWARZ 1.0E-14 + &END + &END + &END XC + &POISSON + PERIODIC NONE + POISSON_SOLVER MT + &END + &END DFT + &SUBSYS + &CELL + ABC 4.0 4.0 4.0 + PERIODIC NONE + &END CELL + &COORD + N 0.400000 0.000000 0.500000 + N -0.400000 0.000000 -0.500000 + &END COORD + &KIND N + BASIS_SET ORB SZV-GTH + BASIS_SET HARRIS DZVP-GTH + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE CELL_OPT + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &CELL_OPT + MAX_ITER 1 + &END +&END + diff --git a/tests/QS/regtest-ec-stress/N2_t11.inp b/tests/QS/regtest-ec-stress/N2_t11.inp new file mode 100644 index 0000000000..ee3279bf0d --- /dev/null +++ b/tests/QS/regtest-ec-stress/N2_t11.inp @@ -0,0 +1,78 @@ +&FORCE_EVAL + METHOD Quickstep + STRESS_TENSOR ANALYTICAL + &PRINT + &FORCES + &END FORCES + &STRESS_TENSOR + COMPONENTS + &END STRESS_TENSOR + &END + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME GTH_POTENTIALS + BASIS_SET_FILE_NAME BASIS_ADMM + &AUXILIARY_DENSITY_MATRIX_METHOD + ADMM_PURIFICATION_METHOD NONE + EXCH_CORRECTION_FUNC NONE + EXCH_SCALING_MODEL NONE + METHOD BASIS_PROJECTION + &END + &MGRID + CUTOFF 200 + &END MGRID + &QS + EPS_DEFAULT 1.E-10 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL HARRIS + HARRIS_BASIS ORBITAL + &RESPONSE_SOLVER + METHOD MO_SOLVER + PRECONDITIONER FULL_SINGLE_INVERSE + &END + &XC + &XC_FUNCTIONAL PBE + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-7 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL PBE0 + &END + &END XC + &POISSON + PERIODIC NONE + POISSON_SOLVER MT + &END + &END DFT + &SUBSYS + &CELL + ABC 4.0 4.0 4.0 + PERIODIC NONE + &END CELL + &COORD + N 0.400000 0.000000 0.500000 + N -0.400000 0.000000 -0.500000 + &END COORD + &KIND N + BASIS_SET ORB DZVP-GTH + BASIS_SET AUX_FIT FIT3 + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE CELL_OPT + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &CELL_OPT + MAX_ITER 1 + &END +&END + diff --git a/tests/QS/regtest-ec-stress/N2_t12.inp b/tests/QS/regtest-ec-stress/N2_t12.inp new file mode 100644 index 0000000000..86c5b8b4c0 --- /dev/null +++ b/tests/QS/regtest-ec-stress/N2_t12.inp @@ -0,0 +1,78 @@ +&FORCE_EVAL + METHOD Quickstep + STRESS_TENSOR ANALYTICAL + &PRINT + &FORCES + &END FORCES + &STRESS_TENSOR + &END STRESS_TENSOR + &END + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME GTH_POTENTIALS + BASIS_SET_FILE_NAME BASIS_ADMM + &AUXILIARY_DENSITY_MATRIX_METHOD + ADMM_PURIFICATION_METHOD NONE + EXCH_CORRECTION_FUNC NONE + EXCH_SCALING_MODEL NONE + METHOD BASIS_PROJECTION + &END + &MGRID + CUTOFF 300 + &END MGRID + &QS + EPS_DEFAULT 1.E-12 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL HARRIS + HARRIS_BASIS HARRIS + &RESPONSE_SOLVER + METHOD MO_SOLVER + PRECONDITIONER FULL_SINGLE_INVERSE + &END + &XC + &XC_FUNCTIONAL PBE + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-7 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL PBE0 + &END + &END XC + &POISSON + PERIODIC NONE + POISSON_SOLVER MT + &END + &END DFT + &SUBSYS + &CELL + ABC 4.0 4.0 4.0 + PERIODIC NONE + &END CELL + &COORD + N 0.400000 0.000000 0.500000 + N -0.400000 0.000000 -0.500000 + &END COORD + &KIND N + BASIS_SET ORB DZVP-GTH + BASIS_SET AUX_FIT FIT3 + BASIS_SET HARRIS TZV2P-GTH + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE CELL_OPT + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &CELL_OPT + MAX_ITER 2 + &END +&END + diff --git a/tests/QS/regtest-ec-stress/N2_t13.inp b/tests/QS/regtest-ec-stress/N2_t13.inp new file mode 100644 index 0000000000..4e9b4b10c1 --- /dev/null +++ b/tests/QS/regtest-ec-stress/N2_t13.inp @@ -0,0 +1,76 @@ +&FORCE_EVAL + METHOD Quickstep + STRESS_TENSOR ANALYTICAL + &PRINT + &FORCES + &END FORCES + &STRESS_TENSOR + &END STRESS_TENSOR + &END + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME GTH_POTENTIALS + BASIS_SET_FILE_NAME BASIS_ADMM + &AUXILIARY_DENSITY_MATRIX_METHOD + ADMM_PURIFICATION_METHOD NONE + EXCH_CORRECTION_FUNC DEFAULT + EXCH_SCALING_MODEL NONE + METHOD BASIS_PROJECTION + &END + &MGRID + CUTOFF 200 + &END MGRID + &QS + EPS_DEFAULT 1.E-10 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL HARRIS + HARRIS_BASIS ORBITAL + &RESPONSE_SOLVER + METHOD MO_SOLVER + PRECONDITIONER FULL_SINGLE_INVERSE + &END + &XC + &XC_FUNCTIONAL PBE + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-6 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL PBE0 + &END + &END XC + &POISSON + PERIODIC NONE + POISSON_SOLVER MT + &END + &END DFT + &SUBSYS + &CELL + ABC 4.0 4.0 4.0 + PERIODIC NONE + &END CELL + &COORD + N 0.400000 0.000000 0.500000 + N -0.400000 0.000000 -0.500000 + &END COORD + &KIND N + BASIS_SET ORB DZVP-GTH + BASIS_SET AUX_FIT FIT3 + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE CELL_OPT + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &CELL_OPT + MAX_ITER 2 + &END +&END diff --git a/tests/QS/regtest-ec-stress/N2_t14.inp b/tests/QS/regtest-ec-stress/N2_t14.inp new file mode 100644 index 0000000000..3bf98ca241 --- /dev/null +++ b/tests/QS/regtest-ec-stress/N2_t14.inp @@ -0,0 +1,78 @@ +&FORCE_EVAL + METHOD Quickstep + STRESS_TENSOR ANALYTICAL + &PRINT + &FORCES + &END FORCES + &STRESS_TENSOR + &END STRESS_TENSOR + &END + &DFT + BASIS_SET_FILE_NAME GTH_BASIS_SETS + POTENTIAL_FILE_NAME GTH_POTENTIALS + BASIS_SET_FILE_NAME BASIS_ADMM + &AUXILIARY_DENSITY_MATRIX_METHOD + ADMM_PURIFICATION_METHOD NONE + EXCH_CORRECTION_FUNC DEFAULT + EXCH_SCALING_MODEL NONE + METHOD BASIS_PROJECTION + &END + &MGRID + CUTOFF 200 + &END MGRID + &QS + EPS_DEFAULT 1.E-10 + &END QS + &ENERGY_CORRECTION + ENERGY_FUNCTIONAL HARRIS + HARRIS_BASIS HARRIS + &RESPONSE_SOLVER + METHOD MO_SOLVER + PRECONDITIONER FULL_SINGLE_INVERSE + &END + &XC + &XC_FUNCTIONAL PBE + &END + &END XC + &END ENERGY_CORRECTION + &SCF + EPS_SCF 1.0E-7 + SCF_GUESS ATOMIC + &END + &XC + &XC_FUNCTIONAL PBE0 + &END + &END XC + &POISSON + PERIODIC NONE + POISSON_SOLVER MT + &END + &END DFT + &SUBSYS + &CELL + ABC 4.0 4.0 4.0 + PERIODIC NONE + &END CELL + &COORD + N 0.400000 0.000000 0.500000 + N -0.400000 0.000000 -0.500000 + &END COORD + &KIND N + BASIS_SET ORB DZVP-GTH + BASIS_SET AUX_FIT FIT3 + BASIS_SET HARRIS TZV2P-GTH + POTENTIAL GTH-PADE-q5 + &END KIND + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT N2 + RUN_TYPE CELL_OPT + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &CELL_OPT + MAX_ITER 2 + &END +&END + diff --git a/tests/QS/regtest-ec-stress/TEST_FILES b/tests/QS/regtest-ec-stress/TEST_FILES new file mode 100644 index 0000000000..0af5f558cb --- /dev/null +++ b/tests/QS/regtest-ec-stress/TEST_FILES @@ -0,0 +1,35 @@ +# 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 +# see regtest/TEST_FILES +# Harris functional energy correction stress tensor tests +# PADE - PADE MO solver +N2_t01.inp 31 1e-05 3.26443566797E+01 +# PADE - PBE +N2_t02.inp 31 5e-05 3.26450578945E+01 +# PBE - PADE +N2_t03.inp 31 5e-05 3.25833150889E+01 +# PBE - PBE +N2_t04.inp 31 1e-05 3.27026970587E+01 +# HFX - PBE +N2_t05.inp 31 5e-05 -4.99918272376E+02 +# PBE0 - PADE +N2_t06.inp 31 5e-05 -5.02611808106E+02 +# PADE - PADE w/ Harris basis +N2_t07.inp 31 1e-05 -5.09706561950E+02 +# PBE - PBE w/ Harris basis +N2_t08.inp 31 5e-05 -4.40020265431E+02 +# PBE0 - PBE w/ Harris basis AO solver +N2_t08b.inp 31 5e-05 -5.04243430029E+02 +# PBE0 - PBE w/ Harris basis +N2_t09.inp 31 5e-05 -5.06200216764E+02 +# HFX - PBE w/ Harris basis +N2_t10.inp 31 1e-05 -4.39582087759E+02 +# HFX-ADMM-NONE - PBE +N2_t11.inp 31 1e-05 -4.99281462627E+02 +# PBE0-ADMM-NONE - PBE w/ Harris basis +N2_t12.inp 31 1e-05 -5.69367982783E+02 +# PBE0-ADMM-DEFAULT - PBE +N2_t13.inp 31 1e-05 -5.68694042522E+02 +# PBE0-ADMM-DEFAULT - PBE w/ Harris basis +N2_t14.inp 31 1e-05 -5.69209232515E+02 +#EOF diff --git a/tests/TEST_DIRS b/tests/TEST_DIRS index c8a4675fc8..bc2b58bd5d 100644 --- a/tests/TEST_DIRS +++ b/tests/TEST_DIRS @@ -3,6 +3,9 @@ # 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-dcdft-force libxc +QS/regtest-dcdft-stress libxc +QS/regtest-ec-stress libxc QS/regtest-ri-rpa-grad libint QS/regtest-sos-mp2-grad libint QS/regtest-ri-rpa libint @@ -11,7 +14,7 @@ QS/regtest-mp2-grad-solvers libint QS/regtest-rpa-cubic-scaling libint QS/regtest-rpa-cubic-scaling-2 libint QS/regtest-double-hybrid-stress-numer-laplace libxc -QS/regtest-ec-force-meta libxc +QS/regtest-ec-meta libxc QS/regtest-double-hybrid-grad-laplace libxc QS/regtest-double-hybrid-stress-numer-meta libxc QS/regtest-double-hybrid-grad-numer-meta libxc