diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 7aaed0a540..3cbc932962 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -246,6 +246,7 @@ list( input_cp2k_colvar.F input_cp2k_constraints.F input_cp2k_dft.F + input_cp2k_harris.F input_cp2k_print_dft.F input_cp2k_scf.F input_cp2k_as.F @@ -575,6 +576,8 @@ list( qs_gcp_utils.F qs_grid_atom.F qs_gspace_mixing.F + qs_harris_types.F + qs_harris_utils.F qs_harmonics_atom.F qs_hash_table_functions.F qs_initial_guess.F diff --git a/src/aobasis/aux_basis_set.F b/src/aobasis/aux_basis_set.F index 3cdf96fce7..8b466d75a4 100644 --- a/src/aobasis/aux_basis_set.F +++ b/src/aobasis/aux_basis_set.F @@ -220,7 +220,6 @@ SUBROUTINE create_aux_basis(aux_basis, bsname, nsets, lmin, lmax, nl, npgf, zet) ALLOCATE (aux_basis%scon(maxco, nsgf)) ALLOCATE (aux_basis%norm_cgf(ncgf)) aux_basis%norm_type = 2 -! CALL init_orb_basis_set(aux_basis) END SUBROUTINE create_aux_basis diff --git a/src/aobasis/basis_set_container_types.F b/src/aobasis/basis_set_container_types.F index 79f2868244..5c7f1047cb 100644 --- a/src/aobasis/basis_set_container_types.F +++ b/src/aobasis/basis_set_container_types.F @@ -45,7 +45,8 @@ MODULE basis_set_container_types p_lri_aux_basis = 116, & aux_opt_basis = 117, & min_basis = 118, & - tda_k_basis = 119 + tda_k_basis = 119, & + rhoin_basis = 120 ! ************************************************************************************************** TYPE basis_set_container_type PRIVATE @@ -133,6 +134,8 @@ FUNCTION get_basis_type(basis_set_type) RESULT(basis_type_nr) basis_type_nr = ri_xas_basis CASE ("AUX_OPT") basis_type_nr = aux_opt_basis + CASE ("RHOIN") + basis_type_nr = rhoin_basis CASE DEFAULT basis_type_nr = unknown_basis END SELECT diff --git a/src/aobasis/basis_set_types.F b/src/aobasis/basis_set_types.F index fca9410a0b..c5e4071eaf 100644 --- a/src/aobasis/basis_set_types.F +++ b/src/aobasis/basis_set_types.F @@ -279,13 +279,15 @@ END SUBROUTINE copy_gto_basis_set !> \brief ... !> \param basis_set ... !> \param pbasis ... +!> \param lmax ... ! ************************************************************************************************** - SUBROUTINE create_primitive_basis_set(basis_set, pbasis) + SUBROUTINE create_primitive_basis_set(basis_set, pbasis, lmax) ! Create a primitives only basis set TYPE(gto_basis_set_type), INTENT(IN) :: basis_set TYPE(gto_basis_set_type), POINTER :: pbasis + INTEGER, INTENT(IN), OPTIONAL :: lmax INTEGER :: i, ico, ip, ipgf, iset, ishell, l, lm, & lshell, m, maxco, mpgf, nc, ncgf, ns, & @@ -332,6 +334,11 @@ SUBROUTINE create_primitive_basis_set(basis_set, pbasis) CALL allocate_gto_basis_set(pbasis) + IF (PRESENT(lmax)) THEN + ! if requested, reduce max l val + lm = MIN(lm, lmax) + END IF + pbasis%name = basis_set%name//"_primitive" pbasis%kind_radius = basis_set%kind_radius pbasis%short_kind_radius = basis_set%short_kind_radius diff --git a/src/atom_fit.F b/src/atom_fit.F index 23737e9053..9c1f489003 100644 --- a/src/atom_fit.F +++ b/src/atom_fit.F @@ -47,6 +47,7 @@ MODULE atom_fit evolt USE powell, ONLY: opt_state_type,& powell_optimize + USE qs_grid_atom, ONLY: grid_atom_type #include "./base/base_uses.f90" IMPLICIT NONE @@ -69,24 +70,27 @@ MODULE atom_fit !> \param num_gto number of Gaussian basis functions !> \param norder ... !> \param iunit output file unit +!> \param agto Gaussian exponents !> \param powell_section POWELL input section !> \param results (output) array of num_gto+2 real numbers in the following format: !> starting exponent, factor of geometrical series, expansion coefficients(1:num_gto) ! ************************************************************************************************** - SUBROUTINE atom_fit_density(atom, num_gto, norder, iunit, powell_section, results) + SUBROUTINE atom_fit_density(atom, num_gto, norder, iunit, agto, powell_section, results) TYPE(atom_type), POINTER :: atom INTEGER, INTENT(IN) :: num_gto, norder, iunit + REAL(KIND=dp), DIMENSION(:), OPTIONAL :: agto TYPE(section_vals_type), OPTIONAL, POINTER :: powell_section REAL(KIND=dp), DIMENSION(:), OPTIONAL :: results - INTEGER :: n10 - REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: co + INTEGER :: i, n10 + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: co, pe REAL(KIND=dp), DIMENSION(2) :: x TYPE(opgrid_type), POINTER :: density TYPE(opt_state_type) :: ostate - ALLOCATE (co(num_gto)) + ALLOCATE (co(num_gto), pe(num_gto)) co = 0._dp + pe = 0._dp NULLIFY (density) CALL create_opgrid(density, atom%basis%grid) CALL atom_denmat(atom%orbitals%pmat, atom%orbitals%wfn, atom%basis%nbas, atom%state%occupation, & @@ -98,52 +102,60 @@ SUBROUTINE atom_fit_density(atom, num_gto, norder, iunit, powell_section, result density%op = density%op*atom%basis%grid%rad**norder END IF - ostate%nf = 0 - ostate%nvar = 2 - x(1) = 0.10_dp !starting point of geometric series - x(2) = 2.00_dp !factor of series - IF (PRESENT(powell_section)) THEN - CALL section_vals_val_get(powell_section, "ACCURACY", r_val=ostate%rhoend) - CALL section_vals_val_get(powell_section, "STEP_SIZE", r_val=ostate%rhobeg) - CALL section_vals_val_get(powell_section, "MAX_FUN", i_val=ostate%maxfun) + IF (PRESENT(agto)) THEN + CPASSERT(num_gto <= SIZE(agto)) + pe(1:num_gto) = agto(1:num_gto) ELSE - ostate%rhoend = 1.e-8_dp - ostate%rhobeg = 5.e-2_dp - ostate%maxfun = 1000 - END IF - ostate%iprint = 1 - ostate%unit = iunit + ostate%nf = 0 + ostate%nvar = 2 + x(1) = 0.10_dp !starting point of geometric series + x(2) = 2.00_dp !factor of series + IF (PRESENT(powell_section)) THEN + CALL section_vals_val_get(powell_section, "ACCURACY", r_val=ostate%rhoend) + CALL section_vals_val_get(powell_section, "STEP_SIZE", r_val=ostate%rhobeg) + CALL section_vals_val_get(powell_section, "MAX_FUN", i_val=ostate%maxfun) + ELSE + ostate%rhoend = 1.e-8_dp + ostate%rhobeg = 5.e-2_dp + ostate%maxfun = 1000 + END IF + ostate%iprint = 1 + ostate%unit = iunit - ostate%state = 0 - IF (iunit > 0) THEN - WRITE (iunit, '(/," POWELL| Start optimization procedure")') - END IF - n10 = ostate%maxfun/10 + ostate%state = 0 + IF (iunit > 0) THEN + WRITE (iunit, '(/," POWELL| Start optimization procedure")') + END IF + n10 = ostate%maxfun/10 - DO + DO - IF (ostate%state == 2) THEN - CALL density_fit(density, atom, num_gto, x(1), x(2), co, ostate%f) - END IF + IF (ostate%state == 2) THEN + pe(1:num_gto) = [(x(1)*x(2)**i, i=1, num_gto)] + CALL density_fit(density%op, density%grid, num_gto, pe, co, ostate%f) + END IF - IF (ostate%state == -1) EXIT + IF (ostate%state == -1) EXIT - CALL powell_optimize(ostate%nvar, x, ostate) + CALL powell_optimize(ostate%nvar, x, ostate) - IF (MOD(ostate%nf, n10) == 0 .AND. iunit > 0) THEN - WRITE (iunit, '(" POWELL| Reached",i4,"% of maximal function calls")') & - INT(REAL(ostate%nf, dp)/REAL(ostate%maxfun, dp)*100._dp) - END IF + IF (MOD(ostate%nf, n10) == 0 .AND. iunit > 0) THEN + WRITE (iunit, '(" POWELL| Reached",i4,"% of maximal function calls")') & + INT(REAL(ostate%nf, dp)/REAL(ostate%maxfun, dp)*100._dp) + END IF - END DO + END DO - ostate%state = 8 - CALL powell_optimize(ostate%nvar, x, ostate) - CALL density_fit(density, atom, num_gto, x(1), x(2), co, ostate%f) + ostate%state = 8 + CALL powell_optimize(ostate%nvar, x, ostate) + pe(1:num_gto) = [(x(1)*x(2)**i, i=1, num_gto)] + END IF + + CALL density_fit(density%op, density%grid, num_gto, pe, co, ostate%f) CALL release_opgrid(density) - IF (iunit > 0) THEN + IF (iunit > 0 .AND. .NOT. PRESENT(agto)) THEN WRITE (iunit, '(" POWELL| Number of function evaluations",T71,I10)') ostate%nf WRITE (iunit, '(" POWELL| Final value of function",T61,G20.10)') ostate%fopt WRITE (iunit, '(" Optimized representation of density using a Geometrical Gaussian basis")') @@ -154,13 +166,18 @@ SUBROUTINE atom_fit_density(atom, num_gto, norder, iunit, powell_section, result END IF IF (PRESENT(results)) THEN - CPASSERT(SIZE(results) >= num_gto + 2) - results(1) = x(1) - results(2) = x(2) - results(3:2 + num_gto) = co(1:num_gto) + IF (PRESENT(agto)) THEN + CPASSERT(SIZE(results) >= num_gto) + results(1:num_gto) = co(1:num_gto) + ELSE + CPASSERT(SIZE(results) >= num_gto + 2) + results(1) = x(1) + results(2) = x(2) + results(3:2 + num_gto) = co(1:num_gto) + END IF END IF - DEALLOCATE (co) + DEALLOCATE (co, pe) END SUBROUTINE atom_fit_density ! ************************************************************************************************** @@ -1286,50 +1303,46 @@ END SUBROUTINE opt_nlcc_param ! ************************************************************************************************** !> \brief Low level routine to fit an atomic electron density. !> \param density electron density on an atomic radial grid 'atom%basis%grid' -!> \param atom information about the atomic kind -!> (only 'atom%basis%grid' is accessed, why not to pass it instead?) +!> \param grid information about the atomic grid !> \param n number of Gaussian basis functions -!> \param aval exponent of the first Gaussian basis function in the series -!> \param cval factor of geometrical series +!> \param pe exponents of Gaussian basis functions !> \param co computed expansion coefficients !> \param aerr error in fitted density \f[\sqrt{\sum_{i}^{nr} (density_fitted(i)-density(i))^2 }\f] ! ************************************************************************************************** - SUBROUTINE density_fit(density, atom, n, aval, cval, co, aerr) - TYPE(opgrid_type), POINTER :: density - TYPE(atom_type), POINTER :: atom + SUBROUTINE density_fit(density, grid, n, pe, co, aerr) + REAL(dp), DIMENSION(:), INTENT(IN) :: density + TYPE(grid_atom_type) :: grid INTEGER, INTENT(IN) :: n - REAL(dp), INTENT(IN) :: aval, cval + REAL(dp), DIMENSION(:), INTENT(IN) :: pe REAL(dp), DIMENSION(:), INTENT(INOUT) :: co REAL(dp), INTENT(OUT) :: aerr INTEGER :: i, info, ip, iq, k, nr INTEGER, ALLOCATABLE, DIMENSION(:) :: ipiv REAL(dp) :: a, rk, zval - REAL(dp), ALLOCATABLE, DIMENSION(:) :: den, pe, uval + REAL(dp), ALLOCATABLE, DIMENSION(:) :: den, uval REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: bf, smat, tval - ALLOCATE (pe(n)) - nr = atom%basis%grid%nr + nr = grid%nr ALLOCATE (bf(nr, n), den(nr)) bf = 0._dp DO i = 1, n - pe(i) = aval*cval**i a = pe(i) DO k = 1, nr - rk = atom%basis%grid%rad(k) + rk = grid%rad(k) bf(k, i) = EXP(-a*rk**2) END DO END DO ! total charge - zval = integrate_grid(density%op, atom%basis%grid) + zval = integrate_grid(density, grid) ! allocate vectors and matrices for overlaps ALLOCATE (tval(n + 1, 1), uval(n), smat(n + 1, n + 1)) DO i = 1, n uval(i) = (pi/pe(i))**1.5_dp - tval(i, 1) = integrate_grid(density%op, bf(:, i), atom%basis%grid) + tval(i, 1) = integrate_grid(density, bf(:, i), grid) END DO tval(n + 1, 1) = zval @@ -1354,10 +1367,10 @@ SUBROUTINE density_fit(density, atom, n, aval, cval, co, aerr) den(:) = den(:) + co(i)*bf(:, i) END DO den(:) = den(:)*fourpi - den(:) = (den(:) - density%op(:))**2 - aerr = SQRT(integrate_grid(den, atom%basis%grid)) + den(:) = (den(:) - density(:))**2 + aerr = SQRT(integrate_grid(den, grid)) - DEALLOCATE (pe, bf, den) + DEALLOCATE (bf, den) DEALLOCATE (tval, uval, smat) diff --git a/src/atom_kind_orbitals.F b/src/atom_kind_orbitals.F index 5e135fb7b1..8ec7022762 100644 --- a/src/atom_kind_orbitals.F +++ b/src/atom_kind_orbitals.F @@ -472,25 +472,27 @@ END SUBROUTINE calculate_atomic_orbitals !> \param qs_kind ... !> \param ngto ... !> \param iunit ... +!> \param optbasis ... Default=T, if basis should be optimized, if not basis is given in input (density) !> \param allelectron ... !> \param confine ... ! ************************************************************************************************** - SUBROUTINE calculate_atomic_density(density, atomic_kind, qs_kind, ngto, iunit, allelectron, confine) + SUBROUTINE calculate_atomic_density(density, atomic_kind, qs_kind, ngto, iunit, & + optbasis, allelectron, confine) REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: density TYPE(atomic_kind_type), POINTER :: atomic_kind TYPE(qs_kind_type), POINTER :: qs_kind INTEGER, INTENT(IN) :: ngto INTEGER, INTENT(IN), OPTIONAL :: iunit - LOGICAL, INTENT(IN), OPTIONAL :: allelectron, confine + LOGICAL, INTENT(IN), OPTIONAL :: optbasis, allelectron, confine INTEGER, PARAMETER :: num_gto = 40 - INTEGER :: i, ii, k, l, ll, m, mb, mo, ngp, nn, nr, & - quadtype, relativistic, z + INTEGER :: i, ii, iw, k, l, ll, m, mb, mo, ngp, nn, & + nr, quadtype, relativistic, z INTEGER, DIMENSION(0:lmat) :: starti INTEGER, DIMENSION(0:lmat, 10) :: ncalc, ncore, nelem INTEGER, DIMENSION(:), POINTER :: econf - LOGICAL :: ecp_semi_local + LOGICAL :: do_basopt, ecp_semi_local REAL(KIND=dp) :: al, aval, cc, cval, ear, rk, xx, zeff REAL(KIND=dp), DIMENSION(num_gto+2) :: results TYPE(all_potential_type), POINTER :: all_potential @@ -513,6 +515,12 @@ SUBROUTINE calculate_atomic_density(density, atomic_kind, qs_kind, ngto, iunit, gth_potential=gth_potential, & sgp_potential=sgp_potential) + IF (PRESENT(iunit)) THEN + iw = iunit + ELSE + iw = -1 + END IF + IF (PRESENT(allelectron)) THEN IF (allelectron) THEN NULLIFY (gth_potential) @@ -520,6 +528,11 @@ SUBROUTINE calculate_atomic_density(density, atomic_kind, qs_kind, ngto, iunit, END IF END IF + do_basopt = .TRUE. + IF (PRESENT(optbasis)) THEN + do_basopt = optbasis + END IF + CPASSERT(ngto <= num_gto) IF (ASSOCIATED(gth_potential) .OR. ASSOCIATED(sgp_potential)) THEN @@ -715,21 +728,21 @@ SUBROUTINE calculate_atomic_density(density, atomic_kind, qs_kind, ngto, iunit, CALL create_atom_orbs(orbitals, mb, mo) CALL set_atom(atom, orbitals=orbitals) - IF (PRESENT(iunit)) THEN - CALL calculate_atom(atom, iunit) - CALL atom_fit_density(atom, ngto, 0, iunit, results=results) + CALL calculate_atom(atom, iw) + + IF (do_basopt) THEN + CALL atom_fit_density(atom, ngto, 0, iw, results=results) + xx = results(1) + cc = results(2) + DO i = 1, ngto + density(i, 1) = xx*cc**i + density(i, 2) = results(2 + i) + END DO ELSE - CALL calculate_atom(atom, -1) - CALL atom_fit_density(atom, ngto, 0, -1, results=results) + CALL atom_fit_density(atom, ngto, 0, iw, agto=density(:, 1), results=results) + density(1:ngto, 2) = results(1:ngto) END IF - xx = results(1) - cc = results(2) - DO i = 1, ngto - density(i, 1) = xx*cc**i - density(i, 2) = results(2 + i) - END DO - ! clean up CALL atom_int_release(integrals) CALL atom_ppint_release(integrals) diff --git a/src/input_constants.F b/src/input_constants.F index 2a15824b27..300b88093a 100644 --- a/src/input_constants.F +++ b/src/input_constants.F @@ -1145,6 +1145,11 @@ MODULE input_constants INTEGER, PARAMETER, PUBLIC :: kg_cholesky = 3001 + ! Harris method + INTEGER, PARAMETER, PUBLIC :: hfun_harris = 1 + INTEGER, PARAMETER, PUBLIC :: hden_atomic = 10 + INTEGER, PARAMETER, PUBLIC :: horb_default = 100 + ! non-scf energy corrections INTEGER, PARAMETER, PUBLIC :: ec_functional_harris = 2001, & ec_functional_dc = 2002 diff --git a/src/input_cp2k_dft.F b/src/input_cp2k_dft.F index c3f4dd875a..1ce2aa7c8f 100644 --- a/src/input_cp2k_dft.F +++ b/src/input_cp2k_dft.F @@ -60,6 +60,7 @@ MODULE input_cp2k_dft create_ext_vxc_section USE input_cp2k_field, ONLY: create_efield_section,& create_per_efield_section + USE input_cp2k_harris, ONLY: create_harris_section USE input_cp2k_kpoints, ONLY: create_kpoints_section USE input_cp2k_loc, ONLY: create_localize_section USE input_cp2k_ls, ONLY: create_ls_scf_section @@ -338,6 +339,10 @@ SUBROUTINE create_dft_section(section) CALL section_add_subsection(section, subsection) CALL section_release(subsection) + CALL create_harris_section(subsection) + CALL section_add_subsection(section, subsection) + CALL section_release(subsection) + CALL create_ec_section(subsection) CALL section_add_subsection(section, subsection) CALL section_release(subsection) diff --git a/src/input_cp2k_harris.F b/src/input_cp2k_harris.F new file mode 100644 index 0000000000..0b68cd873d --- /dev/null +++ b/src/input_cp2k_harris.F @@ -0,0 +1,100 @@ +!--------------------------------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations ! +! Copyright 2000-2024 CP2K developers group ! +! ! +! SPDX-License-Identifier: GPL-2.0-or-later ! +!--------------------------------------------------------------------------------------------------! + +! ************************************************************************************************** +!> \brief Harris input section +! ************************************************************************************************** +MODULE input_cp2k_harris + USE input_constants, ONLY: hden_atomic,& + hfun_harris,& + horb_default + USE input_keyword_types, ONLY: keyword_create,& + keyword_release,& + keyword_type + USE input_section_types, ONLY: section_add_keyword,& + section_create,& + section_type + USE string_utilities, ONLY: s2a +#include "./base/base_uses.f90" + + IMPLICIT NONE + PRIVATE + + CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'input_cp2k_harris' + + PUBLIC :: create_harris_section + +CONTAINS + +! ************************************************************************************************** +!> \brief creates the HARRIS_METHOD section +!> \param section ... +!> \author JGH +! ************************************************************************************************** + SUBROUTINE create_harris_section(section) + TYPE(section_type), POINTER :: section + + TYPE(keyword_type), POINTER :: keyword + + CPASSERT(.NOT. ASSOCIATED(section)) + + NULLIFY (keyword) + CALL section_create(section, __LOCATION__, name="HARRIS_METHOD", & + description="Sets the various options for the Harris method", & + n_keywords=5, n_subsections=0, repeats=.FALSE.) + + CALL keyword_create(keyword, __LOCATION__, name="_SECTION_PARAMETERS_", & + description="Controls the activation of the Harris method", & + usage="&HARRIS_METHOD T", & + default_l_val=.FALSE., & + lone_keyword_l_val=.TRUE.) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="ENERGY_FUNCTIONAL", & + description="Functional used in energy correction", & + usage="ENERGY_FUNCTIONAL HARRIS", & + default_i_val=hfun_harris, & + enum_c_vals=s2a("HARRIS"), & + enum_desc=s2a("Harris functional"), & + enum_i_vals=(/hfun_harris/)) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="DENSITY_SOURCE", & + description="Method to create the input density", & + usage="DENSITY_SOURCE ATOMIC", & + default_i_val=hden_atomic, & + enum_c_vals=s2a("ATOMIC"), & + enum_desc=s2a("Atomic densities"), & + enum_i_vals=(/hden_atomic/)) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="ORBITAL_BASIS", & + description="Specifies the type of basis to be used for the energy functional. ", & + default_i_val=horb_default, & + enum_c_vals=s2a("ATOMIC_KIND_BASIS"), & + enum_desc=s2a("Atomic kind orbital basis"), & + enum_i_vals=(/horb_default/)) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="DEBUG_FORCES", & + description="Additional output to debug Harris method forces.", & + usage="DEBUG_FORCES T", default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + CALL keyword_create(keyword, __LOCATION__, name="DEBUG_STRESS", & + description="Additional output to debug Harris method stress.", & + usage="DEBUG_STRESS T", default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + END SUBROUTINE create_harris_section + +END MODULE input_cp2k_harris diff --git a/src/input_cp2k_scf.F b/src/input_cp2k_scf.F index 9778454f0b..0ca6139a07 100644 --- a/src/input_cp2k_scf.F +++ b/src/input_cp2k_scf.F @@ -215,7 +215,7 @@ SUBROUTINE create_scf_section(section) "Extrapolated from previous RESTART files.", & "Use same guess as MOPAC for semi-empirical methods or a simple diagonal density matrix for other methods", & "Generate a sparse wavefunction using the atomic code (for OT based methods)", & - "Skip initial guess (only for NON-SCC DFTB)."), & + "Skip initial guess (only for non-self consistent methods)."), & enum_i_vals=(/atomic_guess, restart_guess, random_guess, core_guess, & history_guess, mopac_guess, sparse_guess, no_guess/)) CALL section_add_keyword(section, keyword) diff --git a/src/pw_env_methods.F b/src/pw_env_methods.F index fbafa5036e..b544212f68 100644 --- a/src/pw_env_methods.F +++ b/src/pw_env_methods.F @@ -115,7 +115,9 @@ MODULE pw_env_methods CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pw_env_methods' PUBLIC :: pw_env_create, pw_env_rebuild -!*** + +! ************************************************************************************************** + CONTAINS ! ************************************************************************************************** @@ -778,8 +780,8 @@ SUBROUTINE compute_max_radius(radius, pw_env, qs_env) CHARACTER(len=*), PARAMETER :: routineN = 'compute_max_radius' CHARACTER(LEN=8), DIMENSION(4), PARAMETER :: & pbas = (/"ORB ", "AUX_FIT ", "MAO ", "HARRIS "/) - CHARACTER(LEN=8), DIMENSION(8), PARAMETER :: sbas = (/"ORB ", "AUX ", "RI_AUX ", & - "MAO ", "HARRIS ", "RI_HXC ", "RI_K ", "LRI_AUX "/) + CHARACTER(LEN=8), DIMENSION(9), PARAMETER :: sbas = (/"ORB ", "AUX ", "RI_AUX ", & + "MAO ", "HARRIS ", "RI_HXC ", "RI_K ", "LRI_AUX ", "RHOIN "/) REAL(KIND=dp), PARAMETER :: safety_factor = 1.015_dp INTEGER :: handle, ibasis_set_type, igrid_level, igrid_zet0_s, ikind, ipgf, iset, ishell, & @@ -1061,7 +1063,6 @@ SUBROUTINE setup_diel_rs_grid(diel_rs_grid, method, input, pw_grid) border_points=border_points) ALLOCATE (diel_rs_grid) CALL rs_grid_create(diel_rs_grid, rs_desc) -! CALL rs_grid_print(diel_rs_grid,6) CALL rs_grid_release_descriptor(rs_desc) CALL timestop(handle) diff --git a/src/qs_energy.F b/src/qs_energy.F index 3d9c2b43b3..7fd1278fd1 100644 --- a/src/qs_energy.F +++ b/src/qs_energy.F @@ -31,6 +31,7 @@ MODULE qs_energy USE qs_environment_methods, ONLY: qs_env_rebuild_pw_env USE qs_environment_types, ONLY: get_qs_env,& qs_environment_type + USE qs_harris_utils, ONLY: harris_energy_correction USE qs_ks_methods, ONLY: qs_ks_update_qs_env USE qs_matrix_w, ONLY: qs_energies_compute_matrix_w USE qs_nonscf, ONLY: nonscf @@ -133,6 +134,11 @@ SUBROUTINE qs_energies(qs_env, consistent_energies, calc_forces) END IF + ! Check for energy correction + IF (qs_env%harris_method) THEN + CALL harris_energy_correction(qs_env, my_calc_forces) + END IF + ! Do active space calculation CALL active_space_main(qs_env) diff --git a/src/qs_energy_types.F b/src/qs_energy_types.F index aae52bb139..8952633347 100644 --- a/src/qs_energy_types.F +++ b/src/qs_energy_types.F @@ -58,6 +58,7 @@ MODULE qs_energy_types sccs_sol = 0.0_dp, & ! SCCS solvation free energy ktS = 0.0_dp, & ! electronic entropic contribution efermi = 0.0_dp, & ! Fermi energy + band = 0.0_dp, & ! Band energy (Tr PH) dftb3 = 0.0_dp, & ! DFTB 3rd order correction nonscf_correction = 0.0_dp, & ! e.g. Harris correction mp2 = 0.0_dp, & diff --git a/src/qs_environment.F b/src/qs_environment.F index 5d1e99519d..19f987aa80 100644 --- a/src/qs_environment.F +++ b/src/qs_environment.F @@ -24,6 +24,8 @@ MODULE qs_environment create_ri_aux_basis_set USE basis_set_container_types, ONLY: add_basis_set_to_container USE basis_set_types, ONLY: basis_sort_zet,& + create_primitive_basis_set,& + deallocate_gto_basis_set,& get_gto_basis_set,& gto_basis_set_type USE bibliography, ONLY: Iannuzzi2006,& @@ -92,8 +94,8 @@ MODULE qs_environment do_method_gapw, do_method_gapw_xc, do_method_gpw, do_method_lrigpw, do_method_mndo, & do_method_mndod, do_method_ofgpw, do_method_pdg, do_method_pm3, do_method_pm6, & do_method_pm6fm, do_method_pnnl, do_method_rigpw, do_method_rm1, do_method_xtb, & - do_qmmm_gauss, do_qmmm_swave, general_roks, kg_tnadd_embed_ri, rel_none, rel_trans_atom, & - vdw_pairpot_dftd2, vdw_pairpot_dftd3, vdw_pairpot_dftd3bj, wfi_aspc_nr, & + do_qmmm_gauss, do_qmmm_swave, general_roks, hden_atomic, kg_tnadd_embed_ri, rel_none, & + rel_trans_atom, vdw_pairpot_dftd2, vdw_pairpot_dftd3, vdw_pairpot_dftd3bj, wfi_aspc_nr, & wfi_linear_ps_method_nr, wfi_linear_wf_method_nr, wfi_ps_method_nr, & wfi_use_guess_method_nr, xc_vdw_fun_none, xc_vdw_fun_nonloc, xc_vdw_fun_pairpot USE input_section_types, ONLY: section_vals_get,& @@ -150,6 +152,10 @@ MODULE qs_environment USE qs_gcp_types, ONLY: qs_gcp_type USE qs_gcp_utils, ONLY: qs_gcp_env_set,& qs_gcp_init + USE qs_harris_types, ONLY: harris_rhoin_init,& + harris_type + USE qs_harris_utils, ONLY: harris_env_create,& + harris_write_input USE qs_interactions, ONLY: init_interaction_radii,& init_se_nlradius,& write_core_charge_radii,& @@ -258,8 +264,10 @@ SUBROUTINE qs_init(qs_env, para_env, root_section, globenv, cp_subsys, kpoint_en TYPE(cell_type), POINTER :: my_cell, my_cell_ref TYPE(cp_blacs_env_type), POINTER :: blacs_env TYPE(dft_control_type), POINTER :: dft_control + TYPE(distribution_1d_type), POINTER :: local_particles TYPE(energy_correction_type), POINTER :: ec_env TYPE(excited_energy_type), POINTER :: exstate_env + TYPE(harris_type), POINTER :: harris_env TYPE(kpoint_type), POINTER :: kpoints TYPE(lri_environment_type), POINTER :: lri_env TYPE(particle_type), DIMENSION(:), POINTER :: particle_set @@ -479,6 +487,16 @@ SUBROUTINE qs_init(qs_env, para_env, root_section, globenv, cp_subsys, kpoint_en CALL ls_scf_create(qs_env) END IF + CALL get_qs_env(qs_env, harris_env=harris_env) + IF (qs_env%harris_method) THEN + ! initialize the Harris input density and potential integrals + CALL get_qs_env(qs_env, local_particles=local_particles) + CALL harris_rhoin_init(harris_env%rhoin, "RHOIN", qs_kind_set, atomic_kind_set, & + local_particles, dft_control%nspins) + ! Print information of the HARRIS section + CALL harris_write_input(harris_env) + END IF + NULLIFY (ec_env) dft_section => section_vals_get_subs_vals(qs_env%input, "DFT") CALL section_vals_val_get(dft_section, "ENERGY_CORRECTION%_SECTION_PARAMETERS_", & @@ -597,8 +615,10 @@ SUBROUTINE qs_init_subsys(qs_env, para_env, subsys, cell, cell_ref, use_ref_cell TYPE(fist_nonbond_env_type), POINTER :: se_nonbond_env TYPE(gapw_control_type), POINTER :: gapw_control TYPE(gto_basis_set_type), POINTER :: aux_fit_basis, lri_aux_basis, & - ri_aux_basis_set, ri_hfx_basis, & - ri_xas_basis, tmp_basis_set + rhoin_basis, ri_aux_basis_set, & + ri_hfx_basis, ri_xas_basis, & + tmp_basis_set + TYPE(harris_type), POINTER :: harris_env TYPE(local_rho_type), POINTER :: local_rho_set TYPE(lri_environment_type), POINTER :: lri_env TYPE(mo_set_type), DIMENSION(:), POINTER :: mos, mos_last_converged @@ -624,8 +644,8 @@ SUBROUTINE qs_init_subsys(qs_env, para_env, subsys, cell, cell_ref, use_ref_cell TYPE(scf_control_type), POINTER :: scf_control TYPE(se_taper_type), POINTER :: se_taper TYPE(section_vals_type), POINTER :: dft_section, et_coupling_section, et_ddapc_section, & - ewald_section, lri_section, mp2_section, nl_section, poisson_section, pp_section, & - print_section, qs_section, se_section, tddfpt_section, xc_section, xtb_section + ewald_section, harris_section, lri_section, mp2_section, nl_section, poisson_section, & + pp_section, print_section, qs_section, se_section, tddfpt_section, xc_section, xtb_section TYPE(semi_empirical_control_type), POINTER :: se_control TYPE(semi_empirical_si_type), POINTER :: se_store_int_env TYPE(xtb_control_type), POINTER :: xtb_control @@ -919,9 +939,6 @@ SUBROUTINE qs_init_subsys(qs_env, para_env, subsys, cell, cell_ref, use_ref_cell qs_kind => qs_kind_set(ikind) CALL get_qs_kind(qs_kind, basis_set=ri_hfx_basis, basis_type="RI_HXC") IF (.NOT. (ASSOCIATED(ri_hfx_basis))) THEN - ! RI_HXC basis set is not yet loaded - CALL cp_warn(__LOCATION__, "Automatic Generation of RI_HXC basis. "// & - "This is experimental code.") ! Generate a default basis CALL create_ri_aux_basis_set(ri_hfx_basis, qs_kind, dft_control%auto_basis_ri_hxc) CALL add_basis_set_to_container(qs_kind%basis_sets, ri_hfx_basis, "RI_HXC") @@ -929,6 +946,35 @@ SUBROUTINE qs_init_subsys(qs_env, para_env, subsys, cell, cell_ref, use_ref_cell END DO END IF + ! Harris method + NULLIFY (harris_env) + CALL section_vals_val_get(dft_section, "HARRIS_METHOD%_SECTION_PARAMETERS_", & + l_val=qs_env%harris_method) + harris_section => section_vals_get_subs_vals(dft_section, "HARRIS_METHOD") + CALL harris_env_create(qs_env, harris_env, harris_section) + CALL set_qs_env(qs_env, harris_env=harris_env) + ! + IF (qs_env%harris_method) THEN + CALL get_qs_env(qs_env, nkind=nkind) + ! Check if RI_HXC basis is available, auto-generate if needed + DO ikind = 1, nkind + NULLIFY (tmp_basis_set) + qs_kind => qs_kind_set(ikind) + CALL get_qs_kind(qs_kind, basis_set=rhoin_basis, basis_type="RHOIN") + IF (.NOT. (ASSOCIATED(rhoin_basis))) THEN + ! Generate a default basis + CALL create_ri_aux_basis_set(tmp_basis_set, qs_kind, dft_control%auto_basis_ri_hxc) + IF (qs_env%harris_env%density_source == hden_atomic) THEN + CALL create_primitive_basis_set(tmp_basis_set, rhoin_basis, lmax=0) + CALL deallocate_gto_basis_set(tmp_basis_set) + ELSE + rhoin_basis => tmp_basis_set + END IF + CALL add_basis_set_to_container(qs_kind%basis_sets, rhoin_basis, "RHOIN") + END IF + END DO + END IF + mp2_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%WF_CORRELATION") CALL section_vals_get(mp2_section, explicit=mp2_present) IF (mp2_present) THEN @@ -1075,9 +1121,11 @@ SUBROUTINE qs_init_subsys(qs_env, para_env, subsys, cell, cell_ref, use_ref_cell IF (dft_control%qs_control%dftb) THEN scf_control%non_selfconsistent = .NOT. dft_control%qs_control%dftb_control%self_consistent END IF + IF (qs_env%harris_method) THEN + scf_control%non_selfconsistent = .TRUE. + END IF CALL scf_c_create(scf_control) CALL scf_c_read_parameters(scf_control, dft_section) - ! *** Allocate the data structure for Quickstep energies *** CALL allocate_qs_energy(energy) diff --git a/src/qs_environment_types.F b/src/qs_environment_types.F index 98d8f3d674..cc43487430 100644 --- a/src/qs_environment_types.F +++ b/src/qs_environment_types.F @@ -98,6 +98,8 @@ MODULE qs_environment_types USE qs_force_types, ONLY: qs_force_type USE qs_gcp_types, ONLY: qs_gcp_release,& qs_gcp_type + USE qs_harris_types, ONLY: harris_env_release,& + harris_type USE qs_kind_types, ONLY: qs_kind_type USE qs_ks_qmmm_types, ONLY: qs_ks_qmmm_env_type,& qs_ks_qmmm_release @@ -224,6 +226,7 @@ MODULE qs_environment_types LOGICAL :: single_point_run = .FALSE. LOGICAL :: given_embed_pot = .FALSE. LOGICAL :: energy_correction = .FALSE. + LOGICAL :: harris_method = .FALSE. REAL(KIND=dp) :: sim_time = -1.0_dp REAL(KIND=dp) :: start_time = -1.0_dp, target_time = -1.0_dp REAL(KIND=dp), DIMENSION(:, :), POINTER :: image_matrix => NULL() @@ -234,15 +237,15 @@ MODULE qs_environment_types TYPE(almo_scf_env_type), POINTER :: almo_scf_env => NULL() TYPE(transport_env_type), POINTER :: transport_env => NULL() TYPE(cell_type), POINTER :: super_cell => NULL() - TYPE(mo_set_type), DIMENSION(:), POINTER :: mos => NULL() - TYPE(cp_fm_type), DIMENSION(:), POINTER :: mo_loc_history => NULL() + TYPE(mo_set_type), DIMENSION(:), POINTER :: mos => NULL() + TYPE(cp_fm_type), DIMENSION(:), POINTER :: mo_loc_history => NULL() TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mo_derivs => NULL() TYPE(scf_control_type), POINTER :: scf_control => NULL() TYPE(rel_control_type), POINTER :: rel_control => NULL() ! ZMP adding variables TYPE(qs_rho_type), POINTER :: rho_external => NULL() - TYPE(pw_r3d_rs_type), POINTER :: external_vxc => NULL() - TYPE(pw_r3d_rs_type), POINTER :: mask => NULL() + TYPE(pw_r3d_rs_type), POINTER :: external_vxc => NULL() + TYPE(pw_r3d_rs_type), POINTER :: mask => NULL() TYPE(qs_charges_type), POINTER :: qs_charges => NULL() TYPE(qs_ks_env_type), POINTER :: ks_env => NULL() TYPE(qs_ks_qmmm_env_type), POINTER :: ks_qmmm_env => NULL() @@ -270,6 +273,8 @@ MODULE qs_environment_types ! LRI TYPE(lri_environment_type), POINTER :: lri_env => NULL() TYPE(lri_density_type), POINTER :: lri_density => NULL() + ! Harris model + TYPE(harris_type), POINTER :: harris_env => NULL() ! Energy correction TYPE(energy_correction_type), POINTER :: ec_env => NULL() ! Excited States @@ -299,15 +304,15 @@ MODULE qs_environment_types ! Subsystem densities TYPE(qs_rho_p_type), DIMENSION(:), POINTER :: subsys_dens => NULL() ! Embedding potential - TYPE(pw_r3d_rs_type), POINTER :: embed_pot => NULL() - TYPE(pw_r3d_rs_type), POINTER :: spin_embed_pot => NULL() + TYPE(pw_r3d_rs_type), POINTER :: embed_pot => NULL() + TYPE(pw_r3d_rs_type), POINTER :: spin_embed_pot => NULL() ! Polarizability tensor TYPE(polar_env_type), POINTER :: polar_env => NULL() ! Resp charges REAL(KIND=dp), DIMENSION(:), POINTER :: rhs => NULL() REAL(KIND=dp) :: total_zeff_corr = -1.0_dp, surface_dipole_moment = -1.0_dp LOGICAL :: surface_dipole_switch_off = .FALSE. - TYPE(mo_set_type), DIMENSION(:), POINTER :: mos_last_converged => NULL() + TYPE(mo_set_type), DIMENSION(:), POINTER :: mos_last_converged => NULL() END TYPE qs_environment_type CONTAINS @@ -441,6 +446,7 @@ MODULE qs_environment_types !> \param lri_density ... !> \param exstate_env ... !> \param ec_env ... +!> \param harris_env ... !> \param dispersion_env ... !> \param gcp_env ... !> \param vee ... @@ -494,7 +500,7 @@ SUBROUTINE get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, ce neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, & outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, & se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, & - lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, & + lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, & rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, & WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, & s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, & @@ -605,6 +611,7 @@ SUBROUTINE get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, ce TYPE(lri_density_type), OPTIONAL, POINTER :: lri_density TYPE(excited_energy_type), OPTIONAL, POINTER :: exstate_env TYPE(energy_correction_type), OPTIONAL, POINTER :: ec_env + TYPE(harris_type), OPTIONAL, POINTER :: harris_env TYPE(qs_dispersion_type), OPTIONAL, POINTER :: dispersion_env TYPE(qs_gcp_type), OPTIONAL, POINTER :: gcp_env TYPE(pw_r3d_rs_type), OPTIONAL, POINTER :: vee @@ -675,6 +682,7 @@ SUBROUTINE get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, ce IF (PRESENT(se_nonbond_env)) se_nonbond_env => qs_env%se_nonbond_env IF (PRESENT(lri_env)) lri_env => qs_env%lri_env IF (PRESENT(lri_density)) lri_density => qs_env%lri_density + IF (PRESENT(harris_env)) harris_env => qs_env%harris_env IF (PRESENT(ec_env)) ec_env => qs_env%ec_env IF (PRESENT(exstate_env)) exstate_env => qs_env%exstate_env IF (PRESENT(dispersion_env)) dispersion_env => qs_env%dispersion_env @@ -887,6 +895,7 @@ SUBROUTINE init_qs_env(qs_env, globenv) NULLIFY (qs_env%admm_env) NULLIFY (qs_env%efield) NULLIFY (qs_env%lri_env) + NULLIFY (qs_env%harris_env) NULLIFY (qs_env%ec_env) NULLIFY (qs_env%exstate_env) NULLIFY (qs_env%lri_density) @@ -895,7 +904,6 @@ SUBROUTINE init_qs_env(qs_env, globenv) NULLIFY (qs_env%mp2_env) NULLIFY (qs_env%bs_env) NULLIFY (qs_env%kg_env) - NULLIFY (qs_env%ec_env) NULLIFY (qs_env%WannierCentres) qs_env%outer_scf_ihistory = 0 @@ -997,6 +1005,7 @@ END SUBROUTINE init_qs_env !> \param exstate_env ... !> \param ec_env ... !> \param dispersion_env ... +!> \param harris_env ... !> \param gcp_env ... !> \param mp2_env ... !> \param bs_env ... @@ -1029,7 +1038,7 @@ SUBROUTINE set_qs_env(qs_env, super_cell, & outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, & se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, ls_scf_env, & do_transport, transport_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, & - gcp_env, mp2_env, bs_env, kg_env, force, & + harris_env, gcp_env, mp2_env, bs_env, kg_env, force, & kpoints, WannierCentres, almo_scf_env, gradient_history, variable_history, embed_pot, & spin_embed_pot, polar_env, mos_last_converged, rhs) @@ -1091,6 +1100,7 @@ SUBROUTINE set_qs_env(qs_env, super_cell, & TYPE(excited_energy_type), OPTIONAL, POINTER :: exstate_env TYPE(energy_correction_type), OPTIONAL, POINTER :: ec_env TYPE(qs_dispersion_type), OPTIONAL, POINTER :: dispersion_env + TYPE(harris_type), OPTIONAL, POINTER :: harris_env TYPE(qs_gcp_type), OPTIONAL, POINTER :: gcp_env TYPE(mp2_type), OPTIONAL, POINTER :: mp2_env TYPE(post_scf_bandstructure_type), OPTIONAL, & @@ -1322,6 +1332,7 @@ SUBROUTINE set_qs_env(qs_env, super_cell, & IF (PRESENT(admm_env)) qs_env%admm_env => admm_env IF (PRESENT(lri_env)) qs_env%lri_env => lri_env IF (PRESENT(lri_density)) qs_env%lri_density => lri_density + IF (PRESENT(harris_env)) qs_env%harris_env => harris_env IF (PRESENT(ec_env)) qs_env%ec_env => ec_env IF (PRESENT(exstate_env)) qs_env%exstate_env => exstate_env IF (PRESENT(dispersion_env)) qs_env%dispersion_env => dispersion_env @@ -1548,6 +1559,9 @@ SUBROUTINE qs_env_release(qs_env) CALL lri_density_release(qs_env%lri_density) DEALLOCATE (qs_env%lri_density) END IF + IF (ASSOCIATED(qs_env%harris_env)) THEN + CALL harris_env_release(qs_env%harris_env) + END IF IF (ASSOCIATED(qs_env%ec_env)) THEN CALL ec_env_release(qs_env%ec_env) END IF @@ -1765,6 +1779,9 @@ SUBROUTINE qs_env_part_release(qs_env) CALL lri_density_release(qs_env%lri_density) DEALLOCATE (qs_env%lri_density) END IF + IF (ASSOCIATED(qs_env%harris_env)) THEN + CALL harris_env_release(qs_env%harris_env) + END IF IF (ASSOCIATED(qs_env%ec_env)) THEN CALL ec_env_release(qs_env%ec_env) END IF diff --git a/src/qs_force.F b/src/qs_force.F index e3929ec575..eddaa674d6 100644 --- a/src/qs_force.F +++ b/src/qs_force.F @@ -360,7 +360,11 @@ SUBROUTINE qs_forces(qs_env) t3=dummy_real) END IF END IF - ELSEIF (.NOT. perform_ec) THEN + ELSEIF (perform_ec) THEN + ! energy correction forces postponed + ELSEIF (qs_env%harris_method) THEN + ! Harris method forces already done in harris_energy_correction + ELSE ! Compute grid-based forces CALL qs_ks_update_qs_env(qs_env, calculate_forces=.TRUE.) END IF diff --git a/src/qs_harris_types.F b/src/qs_harris_types.F new file mode 100644 index 0000000000..23fbf14e2b --- /dev/null +++ b/src/qs_harris_types.F @@ -0,0 +1,255 @@ +!--------------------------------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations ! +! Copyright 2000-2024 CP2K developers group ! +! ! +! SPDX-License-Identifier: GPL-2.0-or-later ! +!--------------------------------------------------------------------------------------------------! + +! ************************************************************************************************** +!> \brief Types needed for a for a Harris model calculation +!> \par History +!> 2024.07 created +!> \author JGH +! ************************************************************************************************** +MODULE qs_harris_types + USE atomic_kind_types, ONLY: atomic_kind_type,& + get_atomic_kind_set + USE basis_set_types, ONLY: get_gto_basis_set,& + gto_basis_set_type + USE distribution_1d_types, ONLY: distribution_1d_type + USE kinds, ONLY: default_string_length,& + dp + USE pw_types, ONLY: pw_r3d_rs_type + USE qs_kind_types, ONLY: get_qs_kind,& + qs_kind_type +#include "./base/base_uses.f90" + + IMPLICIT NONE + + PRIVATE + + CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_harris_types' + +! ***************************************************************************** + TYPE rho_vec_type + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: rvecs + END TYPE rho_vec_type + + TYPE harris_rhoin_type + CHARACTER(LEN=default_string_length) :: basis_type = "NDef" + TYPE(rho_vec_type), ALLOCATABLE, DIMENSION(:, :) :: rhovec + TYPE(rho_vec_type), ALLOCATABLE, DIMENSION(:, :) :: intvec + INTEGER :: nspin = 0 + INTEGER :: nbas = 0 + INTEGER, ALLOCATABLE, DIMENSION(:, :) :: basptr + LOGICAL :: frozen = .FALSE. + END TYPE harris_rhoin_type + + TYPE harris_energy_type + REAL(KIND=dp) :: eharris = 0.0_dp + REAL(KIND=dp) :: eband = 0.0_dp + REAL(KIND=dp) :: exc_correction = 0.0_dp + REAL(KIND=dp) :: eh_correction = 0.0_dp + REAL(KIND=dp) :: ewald_correction = 0.0_dp + REAL(KIND=dp) :: dispersion = 0.0_dp + END TYPE harris_energy_type + +! ***************************************************************************** +!> \brief Contains information on the Harris method +!> \par History +!> 07.2024 created +!> \author JGH +! ***************************************************************************** + TYPE harris_type + INTEGER :: energy_functional = 0 + INTEGER :: density_source = 0 + INTEGER :: orbital_basis = 0 + ! + TYPE(harris_energy_type) :: energy + ! + TYPE(harris_rhoin_type) :: rhoin + ! + TYPE(pw_r3d_rs_type) :: vh_rspace = pw_r3d_rs_type() + TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: vxc_rspace => Null() + + ! + LOGICAL :: debug_forces = .FALSE. + LOGICAL :: debug_stress = .FALSE. + END TYPE harris_type +! ************************************************************************************************** + + PUBLIC :: harris_type, harris_energy_type, harris_env_release, & + harris_print_energy, harris_rhoin_type, harris_rhoin_init + +! ************************************************************************************************** + +CONTAINS + +! ************************************************************************************************** + +! ************************************************************************************************** +!> \brief ... +!> \param iounit ... +!> \param energy ... +! ************************************************************************************************** + SUBROUTINE harris_print_energy(iounit, energy) + INTEGER, INTENT(IN) :: iounit + TYPE(harris_energy_type) :: energy + + IF (iounit > 0) THEN + WRITE (UNIT=iounit, FMT="(/,(T2,A))") "HARRIS MODEL ENERGY INFORMATION" + WRITE (UNIT=iounit, FMT="((T3,A,T56,F25.14))") & + "Harris model energy: ", energy%eharris, & + "Band energy: ", energy%eband, & + "Hartree correction energy: ", energy%eh_correction, & + "XC correction energy: ", energy%exc_correction, & + "Ewald sum correction energy: ", energy%ewald_correction, & + "Dispersion energy (pair potential): ", energy%dispersion + END IF + + END SUBROUTINE harris_print_energy + +! ************************************************************************************************** +!> \brief ... +!> \param rhoin ... +!> \param basis_type ... +!> \param qs_kind_set ... +!> \param atomic_kind_set ... +!> \param local_particles ... +!> \param nspin ... +! ************************************************************************************************** + SUBROUTINE harris_rhoin_init(rhoin, basis_type, qs_kind_set, atomic_kind_set, & + local_particles, nspin) + TYPE(harris_rhoin_type) :: rhoin + CHARACTER(LEN=*) :: basis_type + TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set + TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set + TYPE(distribution_1d_type), POINTER :: local_particles + INTEGER, INTENT(IN) :: nspin + + INTEGER :: iatom, ikind, iptr, ispin, natom, nkind, & + nparticle_local, nsgf + INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of, nbasf + TYPE(gto_basis_set_type), POINTER :: basis_set + TYPE(qs_kind_type), POINTER :: qs_kind + + CALL harris_rhoin_release(rhoin) + + rhoin%basis_type = basis_type + rhoin%nspin = nspin + + CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, & + atom_of_kind=atom_of_kind, kind_of=kind_of) + natom = SIZE(atom_of_kind) + nkind = SIZE(qs_kind_set) + + ALLOCATE (nbasf(nkind)) + DO ikind = 1, nkind + qs_kind => qs_kind_set(ikind) + CALL get_qs_kind(qs_kind, basis_set=basis_set, basis_type=basis_type) + CALL get_gto_basis_set(basis_set, nsgf=nsgf) + nbasf(ikind) = nsgf + END DO + + ALLOCATE (rhoin%basptr(natom, 2)) + iptr = 1 + DO iatom = 1, natom + ikind = kind_of(iatom) + rhoin%basptr(iatom, 1) = iptr + iptr = iptr + nbasf(ikind) + rhoin%basptr(iatom, 2) = iptr - 1 + END DO + rhoin%nbas = iptr - 1 + + ALLOCATE (rhoin%rhovec(nkind, nspin)) + DO ispin = 1, nspin + DO ikind = 1, nkind + nsgf = nbasf(ikind) + nparticle_local = local_particles%n_el(ikind) + ALLOCATE (rhoin%rhovec(ikind, ispin)%rvecs(nsgf, nparticle_local)) + END DO + END DO + + ALLOCATE (rhoin%intvec(nkind, nspin)) + DO ispin = 1, nspin + DO ikind = 1, nkind + nsgf = nbasf(ikind) + nparticle_local = local_particles%n_el(ikind) + ALLOCATE (rhoin%intvec(ikind, ispin)%rvecs(nsgf, nparticle_local)) + END DO + END DO + + DEALLOCATE (nbasf) + + END SUBROUTINE harris_rhoin_init + +! ************************************************************************************************** +!> \brief ... +!> \param harris_env ... +! ************************************************************************************************** + SUBROUTINE harris_env_release(harris_env) + TYPE(harris_type), POINTER :: harris_env + + INTEGER :: iab + + IF (ASSOCIATED(harris_env)) THEN + ! + CALL harris_rhoin_release(harris_env%rhoin) + ! + IF (ASSOCIATED(harris_env%vh_rspace%pw_grid)) THEN + CALL harris_env%vh_rspace%release() + END IF + IF (ASSOCIATED(harris_env%vxc_rspace)) THEN + DO iab = 1, SIZE(harris_env%vxc_rspace) + CALL harris_env%vxc_rspace(iab)%release() + END DO + DEALLOCATE (harris_env%vxc_rspace) + END IF + ! + DEALLOCATE (harris_env) + END IF + + NULLIFY (harris_env) + + END SUBROUTINE harris_env_release + +! ************************************************************************************************** +!> \brief ... +!> \param rhoin ... +! ************************************************************************************************** + SUBROUTINE harris_rhoin_release(rhoin) + TYPE(harris_rhoin_type) :: rhoin + + INTEGER :: i, j + + IF (ALLOCATED(rhoin%rhovec)) THEN + DO i = 1, SIZE(rhoin%rhovec, 2) + DO j = 1, SIZE(rhoin%rhovec, 1) + IF (ALLOCATED(rhoin%rhovec(j, i)%rvecs)) THEN + DEALLOCATE (rhoin%rhovec(j, i)%rvecs) + END IF + END DO + END DO + DEALLOCATE (rhoin%rhovec) + END IF + IF (ALLOCATED(rhoin%intvec)) THEN + DO i = 1, SIZE(rhoin%intvec, 2) + DO j = 1, SIZE(rhoin%intvec, 1) + IF (ALLOCATED(rhoin%intvec(j, i)%rvecs)) THEN + DEALLOCATE (rhoin%intvec(j, i)%rvecs) + END IF + END DO + END DO + DEALLOCATE (rhoin%intvec) + END IF + IF (ALLOCATED(rhoin%basptr)) THEN + DEALLOCATE (rhoin%basptr) + END IF + rhoin%basis_type = "NDef" + rhoin%nspin = 0 + rhoin%nbas = 0 + rhoin%frozen = .FALSE. + + END SUBROUTINE harris_rhoin_release + +END MODULE qs_harris_types diff --git a/src/qs_harris_utils.F b/src/qs_harris_utils.F new file mode 100644 index 0000000000..88184bdc62 --- /dev/null +++ b/src/qs_harris_utils.F @@ -0,0 +1,736 @@ +!--------------------------------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations ! +! Copyright 2000-2024 CP2K developers group ! +! ! +! SPDX-License-Identifier: GPL-2.0-or-later ! +!--------------------------------------------------------------------------------------------------! + +! ************************************************************************************************** +!> \brief Harris method environment setup and handling +!> \par History +!> 2024.07 created +!> \author JGH +! ************************************************************************************************** +MODULE qs_harris_utils + USE atom_kind_orbitals, ONLY: calculate_atomic_density + USE atomic_kind_types, ONLY: atomic_kind_type + USE basis_set_types, ONLY: get_gto_basis_set,& + gto_basis_set_type + USE cell_types, ONLY: cell_type + USE cp_control_types, ONLY: dft_control_type + USE cp_dbcsr_api, ONLY: dbcsr_copy,& + dbcsr_create,& + dbcsr_p_type,& + dbcsr_release,& + dbcsr_set + USE cp_log_handling, ONLY: cp_get_default_logger,& + cp_logger_get_default_io_unit,& + cp_logger_get_default_unit_nr,& + cp_logger_type + USE distribution_1d_types, ONLY: distribution_1d_type + USE ec_methods, ONLY: create_kernel + USE input_constants, ONLY: hden_atomic,& + hfun_harris,& + horb_default + USE input_section_types, ONLY: section_vals_get_subs_vals,& + section_vals_type,& + section_vals_val_get + USE kinds, ONLY: dp + USE message_passing, ONLY: mp_para_env_type + USE particle_types, ONLY: particle_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,& + pw_integrate_function,& + pw_scale,& + pw_transfer,& + pw_zero + USE pw_poisson_methods, ONLY: pw_poisson_solve + USE pw_poisson_types, ONLY: pw_poisson_type + USE pw_pool_types, ONLY: pw_pool_type + USE pw_types, ONLY: pw_c1d_gs_type,& + pw_r3d_rs_type + USE qs_collocate_density, ONLY: calculate_rho_elec,& + collocate_function + USE qs_energy_types, ONLY: qs_energy_type + USE qs_environment_types, ONLY: get_qs_env,& + qs_environment_type + USE qs_force_types, ONLY: qs_force_type + USE qs_harris_types, ONLY: harris_energy_type,& + harris_print_energy,& + harris_rhoin_type,& + harris_type + USE qs_integrate_potential, ONLY: integrate_function,& + integrate_v_core_rspace,& + integrate_v_rspace + USE qs_kind_types, ONLY: get_qs_kind,& + qs_kind_type + USE qs_ks_types, ONLY: qs_ks_env_type + USE qs_rho_types, ONLY: qs_rho_get,& + qs_rho_type +#include "./base/base_uses.f90" + + IMPLICIT NONE + + PRIVATE + + CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_harris_utils' + + PUBLIC :: harris_env_create, harris_write_input, harris_density_update, calculate_harris_density, & + harris_energy_correction, harris_set_potentials + +CONTAINS + +! ************************************************************************************************** +!> \brief Allocates and intitializes harris_env +!> \param qs_env The QS environment +!> \param harris_env The Harris method environment (the object to create) +!> \param harris_section The Harris method input section +!> \par History +!> 2024.07 created +!> \author JGH +! ************************************************************************************************** + SUBROUTINE harris_env_create(qs_env, harris_env, harris_section) + TYPE(qs_environment_type), POINTER :: qs_env + TYPE(harris_type), POINTER :: harris_env + TYPE(section_vals_type), OPTIONAL, POINTER :: harris_section + + CPASSERT(.NOT. ASSOCIATED(harris_env)) + ALLOCATE (harris_env) + CALL init_harris_env(qs_env, harris_env, harris_section) + + END SUBROUTINE harris_env_create + +! ************************************************************************************************** +!> \brief Initializes Harris method environment +!> \param qs_env The QS environment +!> \param harris_env The Harris method environment +!> \param harris_section The Harris method input section +!> \par History +!> 2024.07 created +!> \author JGH +! ************************************************************************************************** + SUBROUTINE init_harris_env(qs_env, harris_env, harris_section) + TYPE(qs_environment_type), POINTER :: qs_env + TYPE(harris_type), POINTER :: harris_env + TYPE(section_vals_type), OPTIONAL, POINTER :: harris_section + + CHARACTER(LEN=*), PARAMETER :: routineN = 'init_harris_env' + + INTEGER :: handle, unit_nr + TYPE(cp_logger_type), POINTER :: logger + + CALL timeset(routineN, handle) + + IF (qs_env%harris_method) THEN + + CPASSERT(PRESENT(harris_section)) + ! get a useful output_unit + logger => cp_get_default_logger() + IF (logger%para_env%is_source()) THEN + unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.) + ELSE + unit_nr = -1 + END IF + + CALL section_vals_val_get(harris_section, "ENERGY_FUNCTIONAL", & + i_val=harris_env%energy_functional) + CALL section_vals_val_get(harris_section, "DENSITY_SOURCE", & + i_val=harris_env%density_source) + CALL section_vals_val_get(harris_section, "ORBITAL_BASIS", & + i_val=harris_env%orbital_basis) + ! + CALL section_vals_val_get(harris_section, "DEBUG_FORCES", & + l_val=harris_env%debug_forces) + CALL section_vals_val_get(harris_section, "DEBUG_STRESS", & + l_val=harris_env%debug_stress) + + END IF + + CALL timestop(handle) + + END SUBROUTINE init_harris_env + +! ************************************************************************************************** +!> \brief Print out the Harris method input section +!> +!> \param harris_env ... +!> \par History +!> 2024.07 created [JGH] +!> \author JGH +! ************************************************************************************************** + SUBROUTINE harris_write_input(harris_env) + TYPE(harris_type), POINTER :: harris_env + + CHARACTER(LEN=*), PARAMETER :: routineN = 'harris_write_input' + + INTEGER :: handle, unit_nr + TYPE(cp_logger_type), POINTER :: logger + + CALL timeset(routineN, handle) + + logger => cp_get_default_logger() + IF (logger%para_env%is_source()) THEN + unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.) + ELSE + unit_nr = -1 + END IF + + IF (unit_nr > 0) THEN + + WRITE (unit_nr, '(/,T2,A)') & + "!"//REPEAT("-", 29)//" Harris Model "//REPEAT("-", 29)//"!" + + ! Type of energy functional + SELECT CASE (harris_env%energy_functional) + CASE (hfun_harris) + WRITE (unit_nr, '(T2,A,T61,A20)') "Energy Functional: ", "Harris" + END SELECT + ! density source + SELECT CASE (harris_env%density_source) + CASE (hden_atomic) + WRITE (unit_nr, '(T2,A,T61,A20)') "Harris model density: Type", " Atomic kind density" + END SELECT + WRITE (unit_nr, '(T2,A,T71,A10)') "Harris model density: Basis type", & + ADJUSTR(TRIM(harris_env%rhoin%basis_type)) + WRITE (unit_nr, '(T2,A,T71,I10)') "Harris model density: Number of basis functions", & + harris_env%rhoin%nbas + ! orbital basis + SELECT CASE (harris_env%orbital_basis) + CASE (horb_default) + WRITE (unit_nr, '(T2,A,T61,A20)') "Harris model basis: ", "Atomic kind orbitals" + END SELECT + + WRITE (unit_nr, '(T2,A)') REPEAT("-", 79) + WRITE (unit_nr, '()') + + END IF ! unit_nr + + CALL timestop(handle) + + END SUBROUTINE harris_write_input + +! ************************************************************************************************** +!> \brief ... +!> \param qs_env ... +!> \param harris_env ... +! ************************************************************************************************** + SUBROUTINE harris_density_update(qs_env, harris_env) + TYPE(qs_environment_type), POINTER :: qs_env + TYPE(harris_type), POINTER :: harris_env + + CHARACTER(LEN=*), PARAMETER :: routineN = 'harris_density_update' + + INTEGER :: handle, i, ikind, ngto, nkind, nset, nsgf + INTEGER, DIMENSION(:), POINTER :: lmax, npgf + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: coef + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: density + REAL(KIND=dp), DIMENSION(:), POINTER :: norm + REAL(KIND=dp), DIMENSION(:, :), POINTER :: zet + REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: gcc + TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set + TYPE(atomic_kind_type), POINTER :: atomic_kind + TYPE(gto_basis_set_type), POINTER :: basis_set + TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set + TYPE(qs_kind_type), POINTER :: qs_kind + + CALL timeset(routineN, handle) + + SELECT CASE (harris_env%density_source) + CASE (hden_atomic) + IF (.NOT. harris_env%rhoin%frozen) THEN + CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, & + nkind=nkind) + DO ikind = 1, nkind + atomic_kind => atomic_kind_set(ikind) + qs_kind => qs_kind_set(ikind) + CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, & + basis_type=harris_env%rhoin%basis_type) + CALL get_gto_basis_set(gto_basis_set=basis_set, nset=nset, lmax=lmax, nsgf=nsgf, & + npgf=npgf, norm_cgf=norm, zet=zet, gcc=gcc) + IF (nset /= 1 .OR. lmax(1) /= 0 .OR. npgf(1) /= nsgf) THEN + CPABORT("RHOIN illegal basis type") + END IF + DO i = 1, npgf(1) + IF (SUM(ABS(gcc(1:npgf(1), i, 1))) /= MAXVAL(ABS(gcc(1:npgf(1), i, 1)))) THEN + CPABORT("RHOIN illegal basis type") + END IF + END DO + ! + ngto = npgf(1) + ALLOCATE (density(ngto, 2)) + density(1:ngto, 1) = zet(1:ngto, 1) + density(1:ngto, 2) = 0.0_dp + CALL calculate_atomic_density(density, atomic_kind, qs_kind, ngto, & + optbasis=.FALSE., confine=.TRUE.) + ALLOCATE (coef(ngto)) + DO i = 1, ngto + coef(i) = density(i, 2)/gcc(i, i, 1)/norm(i) + END DO + IF (harris_env%rhoin%nspin == 2) THEN + DO i = 1, SIZE(harris_env%rhoin%rhovec(ikind, 1)%rvecs, 2) + harris_env%rhoin%rhovec(ikind, 1)%rvecs(1:ngto, i) = coef(1:ngto)*0.5_dp + harris_env%rhoin%rhovec(ikind, 2)%rvecs(1:ngto, i) = coef(1:ngto)*0.5_dp + END DO + ELSE + DO i = 1, SIZE(harris_env%rhoin%rhovec(ikind, 1)%rvecs, 2) + harris_env%rhoin%rhovec(ikind, 1)%rvecs(1:ngto, i) = coef(1:ngto) + END DO + END IF + DEALLOCATE (density, coef) + END DO + END IF + CASE DEFAULT + CPABORT("Illeagal value of harris_env%density_source") + END SELECT + + CALL timestop(handle) + + END SUBROUTINE harris_density_update + +! ************************************************************************************************** +!> \brief ... +!> \param qs_env ... +!> \param rhoin ... +!> \param rho_struct ... +! ************************************************************************************************** + SUBROUTINE calculate_harris_density(qs_env, rhoin, rho_struct) + TYPE(qs_environment_type), POINTER :: qs_env + TYPE(harris_rhoin_type), INTENT(IN) :: rhoin + TYPE(qs_rho_type), INTENT(INOUT) :: rho_struct + + CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_harris_density' + + INTEGER :: handle, i1, i2, iatom, ikind, ilocal, & + ispin, n, nkind, nlocal, nspin + REAL(KIND=dp) :: eps_rho_rspace + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: vector + REAL(KIND=dp), DIMENSION(:), POINTER :: total_rho + TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set + TYPE(cell_type), POINTER :: cell + TYPE(dft_control_type), POINTER :: dft_control + TYPE(distribution_1d_type), POINTER :: local_particles + TYPE(mp_para_env_type), POINTER :: para_env + TYPE(particle_type), DIMENSION(:), POINTER :: particle_set + TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_gspace + TYPE(pw_env_type), POINTER :: pw_env + TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_rspace + TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set + + CALL timeset(routineN, handle) + + CALL get_qs_env(qs_env, dft_control=dft_control, para_env=para_env) + eps_rho_rspace = dft_control%qs_control%eps_rho_rspace + CALL get_qs_env(qs_env, & + atomic_kind_set=atomic_kind_set, particle_set=particle_set, & + local_particles=local_particles, & + qs_kind_set=qs_kind_set, cell=cell, pw_env=pw_env) + + CALL qs_rho_get(rho_struct, rho_r=rho_rspace, rho_g=rho_gspace, & + tot_rho_r=total_rho) + + ALLOCATE (vector(rhoin%nbas)) + + nkind = SIZE(rhoin%rhovec, 1) + nspin = SIZE(rhoin%rhovec, 2) + + DO ispin = 1, nspin + vector = 0.0_dp + DO ikind = 1, nkind + nlocal = local_particles%n_el(ikind) + DO ilocal = 1, nlocal + iatom = local_particles%list(ikind)%array(ilocal) + i1 = rhoin%basptr(iatom, 1) + i2 = rhoin%basptr(iatom, 2) + n = i2 - i1 + 1 + vector(i1:i2) = rhoin%rhovec(ikind, ispin)%rvecs(1:n, ilocal) + END DO + END DO + CALL para_env%sum(vector) + ! + CALL collocate_function(vector, rho_rspace(ispin), rho_gspace(ispin), & + atomic_kind_set, qs_kind_set, cell, particle_set, pw_env, & + eps_rho_rspace, rhoin%basis_type) + total_rho(ispin) = pw_integrate_function(rho_rspace(ispin), isign=-1) + END DO + + DEALLOCATE (vector) + + CALL timestop(handle) + + END SUBROUTINE calculate_harris_density + +! ************************************************************************************************** +!> \brief ... +!> \param qs_env ... +!> \param rhoin ... +!> \param v_rspace ... +!> \param calculate_forces ... +! ************************************************************************************************** + SUBROUTINE calculate_harris_integrals(qs_env, rhoin, v_rspace, calculate_forces) + TYPE(qs_environment_type), POINTER :: qs_env + TYPE(harris_rhoin_type), INTENT(INOUT) :: rhoin + TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(IN) :: v_rspace + LOGICAL, INTENT(IN) :: calculate_forces + + CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_harris_integrals' + + INTEGER :: handle, i1, i2, iatom, ikind, ilocal, & + ispin, n, nkind, nlocal, nspin + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: integral, vector + TYPE(distribution_1d_type), POINTER :: local_particles + TYPE(mp_para_env_type), POINTER :: para_env + + CALL timeset(routineN, handle) + + CALL get_qs_env(qs_env, para_env=para_env, local_particles=local_particles) + + ALLOCATE (vector(rhoin%nbas)) + ALLOCATE (integral(rhoin%nbas)) + + nkind = SIZE(rhoin%rhovec, 1) + nspin = SIZE(rhoin%rhovec, 2) + + DO ispin = 1, nspin + vector = 0.0_dp + integral = 0.0_dp + DO ikind = 1, nkind + nlocal = local_particles%n_el(ikind) + DO ilocal = 1, nlocal + iatom = local_particles%list(ikind)%array(ilocal) + i1 = rhoin%basptr(iatom, 1) + i2 = rhoin%basptr(iatom, 2) + n = i2 - i1 + 1 + vector(i1:i2) = rhoin%rhovec(ikind, ispin)%rvecs(1:n, ilocal) + END DO + END DO + CALL para_env%sum(vector) + ! + CALL integrate_function(qs_env, v_rspace(ispin), vector, integral, & + calculate_forces, rhoin%basis_type) + DO ikind = 1, nkind + nlocal = local_particles%n_el(ikind) + DO ilocal = 1, nlocal + iatom = local_particles%list(ikind)%array(ilocal) + i1 = rhoin%basptr(iatom, 1) + i2 = rhoin%basptr(iatom, 2) + n = i2 - i1 + 1 + rhoin%intvec(ikind, ispin)%rvecs(1:n, ilocal) = integral(i1:i2) + END DO + END DO + END DO + + DEALLOCATE (vector, integral) + + CALL timestop(handle) + + END SUBROUTINE calculate_harris_integrals + +! ************************************************************************************************** +!> \brief ... +!> \param harris_env ... +!> \param vh_rspace ... +!> \param vxc_rspace ... +! ************************************************************************************************** + SUBROUTINE harris_set_potentials(harris_env, vh_rspace, vxc_rspace) + TYPE(harris_type), POINTER :: harris_env + TYPE(pw_r3d_rs_type), INTENT(IN) :: vh_rspace + TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: vxc_rspace + + INTEGER :: iab, ispin, nspins + TYPE(pw_grid_type), POINTER :: pw_grid + + ! release possible old potentials + IF (ASSOCIATED(harris_env%vh_rspace%pw_grid)) THEN + CALL harris_env%vh_rspace%release() + END IF + IF (ASSOCIATED(harris_env%vxc_rspace)) THEN + DO iab = 1, SIZE(harris_env%vxc_rspace) + CALL harris_env%vxc_rspace(iab)%release() + END DO + DEALLOCATE (harris_env%vxc_rspace) + END IF + + ! generate new potential data structures + nspins = harris_env%rhoin%nspin + ALLOCATE (harris_env%vxc_rspace(nspins)) + + pw_grid => vh_rspace%pw_grid + CALL harris_env%vh_rspace%create(pw_grid) + DO ispin = 1, nspins + CALL harris_env%vxc_rspace(ispin)%create(pw_grid) + END DO + + ! copy potentials + CALL pw_transfer(vh_rspace, harris_env%vh_rspace) + IF (ASSOCIATED(vxc_rspace)) THEN + DO ispin = 1, nspins + CALL pw_transfer(vxc_rspace(ispin), harris_env%vxc_rspace(ispin)) + CALL pw_scale(harris_env%vxc_rspace(ispin), vxc_rspace(ispin)%pw_grid%dvol) + END DO + ELSE + DO ispin = 1, nspins + CALL pw_zero(harris_env%vxc_rspace(ispin)) + END DO + END IF + + END SUBROUTINE harris_set_potentials + +! ************************************************************************************************** +!> \brief ... +!> \param qs_env ... +!> \param calculate_forces ... +! ************************************************************************************************** + SUBROUTINE harris_energy_correction(qs_env, calculate_forces) + TYPE(qs_environment_type), POINTER :: qs_env + LOGICAL, INTENT(IN) :: calculate_forces + + CHARACTER(LEN=*), PARAMETER :: routineN = 'harris_energy_correction' + + INTEGER :: handle, iounit, ispin, nspins + REAL(KIND=dp) :: dvol, ec, eh, exc, vxc + TYPE(cp_logger_type), POINTER :: logger + TYPE(harris_energy_type), POINTER :: energy + TYPE(harris_type), POINTER :: harris_env + TYPE(pw_c1d_gs_type), POINTER :: rho_core + TYPE(pw_env_type), POINTER :: pw_env + TYPE(pw_pool_type), POINTER :: auxbas_pw_pool + TYPE(pw_r3d_rs_type) :: core_rspace + TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r + TYPE(qs_energy_type), POINTER :: ks_energy + TYPE(qs_rho_type), POINTER :: rho + + MARK_USED(calculate_forces) + + CALL timeset(routineN, handle) + + CALL get_qs_env(qs_env, harris_env=harris_env, energy=ks_energy) + energy => harris_env%energy + energy%eband = ks_energy%band + energy%ewald_correction = ks_energy%core_overlap + ks_energy%core_self + energy%dispersion = ks_energy%dispersion + + nspins = harris_env%rhoin%nspin + + CALL get_qs_env(qs_env, rho=rho, rho_core=rho_core) + CALL qs_rho_get(rho, rho_r=rho_r) + + CALL get_qs_env(qs_env=qs_env, pw_env=pw_env) + CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool) + CALL auxbas_pw_pool%create_pw(core_rspace) + CALL pw_transfer(rho_core, core_rspace) + + dvol = harris_env%vh_rspace%pw_grid%dvol + eh = 0.0_dp + DO ispin = 1, nspins + eh = eh + pw_integral_ab(rho_r(ispin), harris_env%vh_rspace)/dvol + END DO + ec = pw_integral_ab(core_rspace, harris_env%vh_rspace)/dvol + eh = 0.5_dp*(eh + ec) + energy%eh_correction = ec - eh + + exc = ks_energy%exc + vxc = 0.0_dp + IF (ASSOCIATED(harris_env%vxc_rspace)) THEN + DO ispin = 1, nspins + vxc = vxc + pw_integral_ab(rho_r(ispin), harris_env%vxc_rspace(ispin))/ & + harris_env%vxc_rspace(ispin)%pw_grid%dvol + END DO + END IF + energy%exc_correction = exc - vxc + + ! Total Harris model energy + energy%eharris = energy%eband + energy%eh_correction + energy%exc_correction + & + energy%ewald_correction + energy%dispersion + + CALL auxbas_pw_pool%give_back_pw(core_rspace) + + ks_energy%total = ks_energy%total + ks_energy%core + ks_energy%nonscf_correction = energy%eharris - ks_energy%total + ks_energy%total = energy%eharris + + logger => cp_get_default_logger() + iounit = cp_logger_get_default_io_unit(logger) + + CALL harris_print_energy(iounit, energy) + + IF (calculate_forces) THEN + CALL harris_forces(qs_env, iounit) + END IF + + CALL timestop(handle) + + END SUBROUTINE harris_energy_correction + +! ************************************************************************************************** +!> \brief ... +!> \param qs_env ... +!> \param iounit ... +! ************************************************************************************************** + SUBROUTINE harris_forces(qs_env, iounit) + TYPE(qs_environment_type), POINTER :: qs_env + INTEGER, INTENT(IN) :: iounit + + CHARACTER(LEN=*), PARAMETER :: routineN = 'harris_forces' + LOGICAL, PARAMETER :: debug_forces = .TRUE. + + INTEGER :: handle, ispin, nspins + REAL(KIND=dp) :: ehartree + REAL(KIND=dp), DIMENSION(3) :: fodeb + TYPE(dbcsr_p_type) :: scrm + TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rhoh_ao, smat + TYPE(harris_type), POINTER :: harris_env + TYPE(mp_para_env_type), POINTER :: para_env + TYPE(pw_c1d_gs_type) :: rhoh_tot_gspace, vhout_gspace + TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g, rhoh_g + TYPE(pw_c1d_gs_type), POINTER :: rho_core + TYPE(pw_env_type), POINTER :: pw_env + TYPE(pw_poisson_type), POINTER :: poisson_env + TYPE(pw_pool_type), POINTER :: auxbas_pw_pool + TYPE(pw_r3d_rs_type) :: vhout_rspace, vhxc_rspace + TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: fhxc_rspace, ftau, fxc, rho_r, rhoh_r, & + tauh_r + 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 :: xc_section + + CALL timeset(routineN, handle) + + IF (debug_forces) THEN + IF (iounit > 0) WRITE (iounit, "(/,T3,A)") & + "DEBUG:: Harris Method Forces (density dependent)" + END IF + + CALL get_qs_env(qs_env, harris_env=harris_env, force=force, para_env=para_env) + nspins = harris_env%rhoin%nspin + + CALL get_qs_env(qs_env, rho=rho, rho_core=rho_core, matrix_s=smat) + ! Warning: rho_ao = output DM; rho_r = rhoin + CALL qs_rho_get(rho, rho_ao=rhoh_ao, rho_r=rho_r, rho_g=rho_g) + ALLOCATE (scrm%matrix) + CALL dbcsr_create(scrm%matrix, template=rhoh_ao(1)%matrix) + CALL dbcsr_copy(scrm%matrix, smat(1)%matrix) + CALL dbcsr_set(scrm%matrix, 0.0_dp) + + CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, ks_env=ks_env) + CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool) + CALL auxbas_pw_pool%create_pw(vhxc_rspace) + + IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1) + DO ispin = 1, nspins + CALL pw_copy(harris_env%vh_rspace, vhxc_rspace) + CALL pw_axpy(harris_env%vxc_rspace(ispin), vhxc_rspace) + CALL integrate_v_rspace(v_rspace=vhxc_rspace, & + hmat=scrm, pmat=rhoh_ao(ispin), & + qs_env=qs_env, calculate_forces=.TRUE.) + END DO + IF (debug_forces) THEN + fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3) + CALL para_env%sum(fodeb) + IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*(Vh[in]+Vxc)", fodeb + END IF + + CALL dbcsr_release(scrm%matrix) + DEALLOCATE (scrm%matrix) + CALL auxbas_pw_pool%give_back_pw(vhxc_rspace) + + ALLOCATE (rhoh_r(nspins), rhoh_g(nspins)) + DO ispin = 1, nspins + CALL auxbas_pw_pool%create_pw(rhoh_r(ispin)) + CALL auxbas_pw_pool%create_pw(rhoh_g(ispin)) + END DO + CALL auxbas_pw_pool%create_pw(rhoh_tot_gspace) + CALL pw_copy(rho_core, rhoh_tot_gspace) + DO ispin = 1, nspins + CALL calculate_rho_elec(ks_env=ks_env, matrix_p=rhoh_ao(ispin)%matrix, & + rho=rhoh_r(ispin), rho_gspace=rhoh_g(ispin)) + CALL pw_axpy(rhoh_g(ispin), rhoh_tot_gspace) + END DO + ! no meta functionals here + NULLIFY (tauh_r) + + CALL auxbas_pw_pool%create_pw(vhout_rspace) + CALL auxbas_pw_pool%create_pw(vhout_gspace) + CALL pw_env_get(pw_env, poisson_env=poisson_env) + ! + CALL pw_poisson_solve(poisson_env, rhoh_tot_gspace, ehartree, vhout_gspace) + ! + CALL pw_transfer(vhout_gspace, vhout_rspace) + CALL pw_scale(vhout_rspace, vhout_rspace%pw_grid%dvol) + + IF (debug_forces) fodeb(1:3) = force(1)%rho_core(1:3, 1) + CALL integrate_v_core_rspace(vhout_rspace, qs_env) + IF (debug_forces) THEN + fodeb(1:3) = force(1)%rho_core(1:3, 1) - fodeb(1:3) + CALL para_env%sum(fodeb) + IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Vh[out]*dncore ", fodeb + END IF + + ALLOCATE (fhxc_rspace(nspins)) + DO ispin = 1, nspins + CALL auxbas_pw_pool%create_pw(fhxc_rspace(ispin)) + END DO + ! vh = vh[out] - vh[in] + CALL pw_axpy(harris_env%vh_rspace, vhout_rspace, alpha=-1._dp, beta=1.0_dp) + ! kernel fxc + ! drho = rho[out] - rho[in] + DO ispin = 1, nspins + CALL pw_axpy(rho_r(ispin), rhoh_r(ispin), alpha=-1._dp, beta=1.0_dp) + CALL pw_axpy(rho_g(ispin), rhoh_g(ispin), alpha=-1._dp, beta=1.0_dp) + END DO + xc_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC") + NULLIFY (fxc, ftau) + CALL create_kernel(qs_env, vxc=fxc, vxc_tau=ftau, & + rho=rho, rho1_r=rhoh_r, rho1_g=rhoh_g, tau1_r=tauh_r, & + xc_section=xc_section) + CPASSERT(.NOT. ASSOCIATED(ftau)) + + DO ispin = 1, nspins + CALL pw_copy(vhout_rspace, fhxc_rspace(ispin)) + IF (ASSOCIATED(fxc)) THEN + CALL pw_scale(fxc(ispin), fxc(ispin)%pw_grid%dvol) + CALL pw_axpy(fxc(ispin), fhxc_rspace(ispin)) + END IF + END DO + + IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1) + CALL calculate_harris_integrals(qs_env, harris_env%rhoin, fhxc_rspace, .TRUE.) + IF (debug_forces) THEN + fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3) + CALL para_env%sum(fodeb) + IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: (dVh+fxc)*dn[in] ", fodeb + END IF + + IF (ASSOCIATED(fxc)) THEN + DO ispin = 1, nspins + CALL auxbas_pw_pool%give_back_pw(fxc(ispin)) + END DO + DEALLOCATE (fxc) + END IF + IF (ASSOCIATED(ftau)) THEN + DO ispin = 1, nspins + CALL auxbas_pw_pool%give_back_pw(ftau(ispin)) + END DO + DEALLOCATE (ftau) + END IF + + CALL auxbas_pw_pool%give_back_pw(rhoh_tot_gspace) + CALL auxbas_pw_pool%give_back_pw(vhout_rspace) + CALL auxbas_pw_pool%give_back_pw(vhout_gspace) + + DO ispin = 1, nspins + CALL auxbas_pw_pool%give_back_pw(rhoh_r(ispin)) + CALL auxbas_pw_pool%give_back_pw(rhoh_g(ispin)) + CALL auxbas_pw_pool%give_back_pw(fhxc_rspace(ispin)) + END DO + DEALLOCATE (rhoh_r, rhoh_g, fhxc_rspace) + + CALL timestop(handle) + + END SUBROUTINE harris_forces + +END MODULE qs_harris_utils diff --git a/src/qs_integrate_potential.F b/src/qs_integrate_potential.F index 90d46fb4df..4cd0c1e984 100644 --- a/src/qs_integrate_potential.F +++ b/src/qs_integrate_potential.F @@ -23,7 +23,8 @@ MODULE qs_integrate_potential USE grid_api, ONLY: integrate_pgf_product USE qs_integrate_potential_product, ONLY: integrate_v_dbasis,& integrate_v_rspace - USE qs_integrate_potential_single, ONLY: integrate_ppl_rspace,& + USE qs_integrate_potential_single, ONLY: integrate_function,& + integrate_ppl_rspace,& integrate_rho_nlcc,& integrate_v_core_rspace,& integrate_v_gaussian_rspace,& @@ -48,6 +49,7 @@ MODULE qs_integrate_potential integrate_v_rspace_diagonal, & integrate_v_core_rspace, & integrate_v_gaussian_rspace, & + integrate_function, & integrate_ppl_rspace, & integrate_rho_nlcc diff --git a/src/qs_integrate_potential_single.F b/src/qs_integrate_potential_single.F index f56793886a..119105337f 100644 --- a/src/qs_integrate_potential_single.F +++ b/src/qs_integrate_potential_single.F @@ -30,7 +30,8 @@ MODULE qs_integrate_potential_single USE ao_util, ONLY: exp_radius_very_extended USE atomic_kind_types, ONLY: atomic_kind_type,& - get_atomic_kind + get_atomic_kind,& + get_atomic_kind_set USE atprop_types, ONLY: atprop_array_init,& atprop_type USE basis_set_types, ONLY: get_gto_basis_set,& @@ -59,6 +60,7 @@ MODULE qs_integrate_potential_single qs_environment_type USE qs_force_types, ONLY: qs_force_type USE qs_kind_types, ONLY: get_qs_kind,& + get_qs_kind_set,& qs_kind_type USE realspace_grid_types, ONLY: map_gaussian_here,& realspace_grid_type,& @@ -84,6 +86,7 @@ MODULE qs_integrate_potential_single integrate_v_rspace_diagonal, & integrate_v_core_rspace, & integrate_v_gaussian_rspace, & + integrate_function, & integrate_ppl_rspace, & integrate_rho_nlcc @@ -1175,4 +1178,184 @@ SUBROUTINE integrate_v_rspace_diagonal(v_rspace, ksmat, pmat, qs_env, calculate_ END SUBROUTINE integrate_v_rspace_diagonal +! ************************************************************************************************** +!> \brief computes integrals of product of v_rspace times a basis function (vector function) +!> and possible forces +!> \param qs_env ... +!> \param v_rspace ... +!> \param f_coef ... +!> \param f_integral ... +!> \param calculate_forces ... +!> \param basis_type ... +!> \author JGH [8.2024] +! ************************************************************************************************** + SUBROUTINE integrate_function(qs_env, v_rspace, f_coef, f_integral, calculate_forces, basis_type) + TYPE(qs_environment_type), POINTER :: qs_env + TYPE(pw_r3d_rs_type), INTENT(IN) :: v_rspace + REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: f_coef + REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: f_integral + LOGICAL, INTENT(IN) :: calculate_forces + CHARACTER(len=*), INTENT(IN) :: basis_type + + CHARACTER(len=*), PARAMETER :: routineN = 'integrate_function' + + INTEGER :: atom_a, group_size, handle, i, iatom, igrid_level, ikind, ipgf, iset, maxco, & + maxsgf_set, my_pos, na1, natom, ncoa, nkind, nseta, offset, sgfa + INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind + INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgfa, nsgfa + INTEGER, DIMENSION(:, :), POINTER :: first_sgfa + LOGICAL :: use_virial + REAL(KIND=dp) :: eps_rho_rspace, radius + REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, ra + REAL(KIND=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b + REAL(KIND=dp), DIMENSION(:, :), POINTER :: hab, pab, sphi_a, work, zeta + TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set + TYPE(cell_type), POINTER :: cell + TYPE(dft_control_type), POINTER :: dft_control + TYPE(gridlevel_info_type), POINTER :: gridlevel_info + TYPE(gto_basis_set_type), POINTER :: orb_basis_set + TYPE(particle_type), DIMENSION(:), POINTER :: particle_set + TYPE(pw_env_type), POINTER :: pw_env + TYPE(qs_force_type), DIMENSION(:), POINTER :: force + TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set + TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_v + TYPE(realspace_grid_type), POINTER :: rs_grid + TYPE(virial_type), POINTER :: virial + + CALL timeset(routineN, handle) + + CALL get_qs_env(qs_env=qs_env, pw_env=pw_env) + gridlevel_info => pw_env%gridlevel_info + + CALL pw_env_get(pw_env, rs_grids=rs_v) + CALL potential_pw2rs(rs_v, v_rspace, pw_env) + + CALL get_qs_env(qs_env=qs_env, & + atomic_kind_set=atomic_kind_set, & + qs_kind_set=qs_kind_set, & + force=force, & + cell=cell, & + dft_control=dft_control, & + nkind=nkind, & + natom=natom, & + particle_set=particle_set, & + pw_env=pw_env, & + virial=virial) + CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind) + + use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer) + IF (use_virial) THEN + CPABORT("Virial NYA") + END IF + + eps_rho_rspace = dft_control%qs_control%eps_rho_rspace + + CALL get_qs_kind_set(qs_kind_set, & + maxco=maxco, maxsgf_set=maxsgf_set, basis_type=basis_type) + ALLOCATE (hab(maxco, 1), pab(maxco, 1), work(maxco, 1)) + + offset = 0 + my_pos = v_rspace%pw_grid%para%group%mepos + group_size = v_rspace%pw_grid%para%group%num_pe + + DO iatom = 1, natom + ikind = particle_set(iatom)%atomic_kind%kind_number + atom_a = atom_of_kind(iatom) + CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, basis_type=basis_type) + CALL get_gto_basis_set(gto_basis_set=orb_basis_set, & + first_sgf=first_sgfa, & + lmax=la_max, & + lmin=la_min, & + npgf=npgfa, & + nset=nseta, & + nsgf_set=nsgfa, & + sphi=sphi_a, & + zet=zeta) + ra(:) = pbc(particle_set(iatom)%r, cell) + + force_a(:) = 0._dp + force_b(:) = 0._dp + my_virial_a(:, :) = 0._dp + my_virial_b(:, :) = 0._dp + + DO iset = 1, nseta + + ncoa = npgfa(iset)*ncoset(la_max(iset)) + sgfa = first_sgfa(1, iset) + + hab = 0._dp + pab = 0._dp + + DO i = 1, nsgfa(iset) + work(i, 1) = f_coef(offset + i) + END DO + + CALL dgemm("N", "N", ncoa, 1, nsgfa(iset), & + 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), & + work(1, 1), SIZE(work, 1), & + 0.0_dp, pab(1, 1), SIZE(pab, 1)) + + DO ipgf = 1, npgfa(iset) + + na1 = (ipgf - 1)*ncoset(la_max(iset)) + + igrid_level = gaussian_gridlevel(gridlevel_info, zeta(ipgf, iset)) + rs_grid => rs_v(igrid_level) + + IF (map_gaussian_here(rs_grid, cell%h_inv, ra, offset, group_size, my_pos)) THEN + radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), & + lb_min=0, lb_max=0, ra=ra, rb=ra, rp=ra, & + zetp=zeta(ipgf, iset), eps=eps_rho_rspace, & + prefactor=1.0_dp, cutoff=1.0_dp) + + IF (calculate_forces) THEN + CALL integrate_pgf_product(la_max=la_max(iset), & + zeta=zeta(ipgf, iset), la_min=la_min(iset), & + lb_max=0, zetb=0.0_dp, lb_min=0, & + ra=ra, rab=(/0.0_dp, 0.0_dp, 0.0_dp/), & + rsgrid=rs_grid, & + hab=hab, pab=pab, o1=na1, o2=0, radius=radius, & + calculate_forces=.TRUE., & + force_a=force_a, force_b=force_b, & + use_virial=use_virial, & + my_virial_a=my_virial_a, my_virial_b=my_virial_b) + ELSE + CALL integrate_pgf_product(la_max=la_max(iset), & + zeta=zeta(ipgf, iset), la_min=la_min(iset), & + lb_max=0, zetb=0.0_dp, lb_min=0, & + ra=ra, rab=(/0.0_dp, 0.0_dp, 0.0_dp/), & + rsgrid=rs_grid, & + hab=hab, o1=na1, o2=0, radius=radius, & + calculate_forces=.FALSE.) + END IF + + END IF + + END DO + ! + CALL dgemm("T", "N", nsgfa(iset), 1, ncoa, 1.0_dp, sphi_a(1, sgfa), & + SIZE(sphi_a, 1), hab(1, 1), SIZE(hab, 1), 0.0_dp, work(1, 1), SIZE(work, 1)) + DO i = 1, nsgfa(iset) + f_integral(offset + i) = work(i, 1) + END DO + + offset = offset + nsgfa(iset) + + END DO + + IF (calculate_forces) THEN + force(ikind)%rho_elec(:, atom_a) = force(ikind)%rho_elec(:, atom_a) + force_a(:) + IF (use_virial) THEN + virial%pv_virial = virial%pv_virial + my_virial_a + END IF + END IF + + END DO + + DEALLOCATE (hab, pab, work) + + CALL timestop(handle) + + END SUBROUTINE integrate_function + END MODULE qs_integrate_potential_single diff --git a/src/qs_interactions.F b/src/qs_interactions.F index 33e45c69f4..d43ac6c207 100644 --- a/src/qs_interactions.F +++ b/src/qs_interactions.F @@ -101,8 +101,8 @@ SUBROUTINE init_interaction_radii(qs_control, qs_kind_set) TYPE(gth_potential_type), POINTER :: gth_potential TYPE(gto_basis_set_type), POINTER :: aux_basis_set, aux_fit_basis_set, aux_gw_basis, & aux_opt_basis_set, gapw_1c_basis, harris_basis, lri_basis, mao_basis, min_basis_set, & - orb_basis_set, p_lri_basis, ri_aux_basis_set, ri_basis, ri_xas_basis, soft_basis, & - tda_k_basis + orb_basis_set, p_lri_basis, rhoin_basis, ri_aux_basis_set, ri_basis, ri_xas_basis, & + soft_basis, tda_k_basis TYPE(paw_proj_set_type), POINTER :: paw_proj_set TYPE(sgp_potential_type), POINTER :: sgp_potential @@ -136,6 +136,7 @@ SUBROUTINE init_interaction_radii(qs_control, qs_kind_set) CALL get_qs_kind(qs_kind_set(ikind), basis_set=soft_basis, basis_type="ORB_SOFT") CALL get_qs_kind(qs_kind_set(ikind), basis_set=gapw_1c_basis, basis_type="GAPW_1C") CALL get_qs_kind(qs_kind_set(ikind), basis_set=tda_k_basis, basis_type="TDA_HFX") + CALL get_qs_kind(qs_kind_set(ikind), basis_set=rhoin_basis, basis_type="RHOIN") CALL get_qs_kind(qs_kind_set(ikind), & paw_proj_set=paw_proj_set, & paw_atom=paw_atom, & @@ -381,6 +382,11 @@ SUBROUTINE init_interaction_radii(qs_control, qs_kind_set) CALL init_interaction_radii_orb_basis(p_lri_basis, qs_control%eps_pgf_orb) END IF + ! Calculate the density input basis function radii + IF (ASSOCIATED(rhoin_basis)) THEN + CALL init_interaction_radii_orb_basis(rhoin_basis, qs_control%eps_pgf_orb) + END IF + ! Calculate the paw projector basis function radii IF (ASSOCIATED(paw_proj_set)) THEN CALL get_paw_proj_set(paw_proj_set=paw_proj_set, & diff --git a/src/qs_ks_methods.F b/src/qs_ks_methods.F index 9a23526fe9..1cd051f063 100644 --- a/src/qs_ks_methods.F +++ b/src/qs_ks_methods.F @@ -99,6 +99,8 @@ MODULE qs_ks_methods USE qs_environment_types, ONLY: get_qs_env,& qs_environment_type USE qs_gapw_densities, ONLY: prepare_gapw_den + USE qs_harris_types, ONLY: harris_type + USE qs_harris_utils, ONLY: harris_set_potentials USE qs_integrate_potential, ONLY: integrate_ppl_rspace,& integrate_rho_nlcc,& integrate_v_core_rspace @@ -195,6 +197,7 @@ SUBROUTINE qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces, just_energy, & matrix_h_im, matrix_s, my_rho, rho_ao TYPE(dft_control_type), POINTER :: dft_control TYPE(ecoul_1center_type), DIMENSION(:), POINTER :: ecoul_1c + TYPE(harris_type), POINTER :: harris_env TYPE(local_rho_type), POINTER :: local_rho_set TYPE(lri_density_type), POINTER :: lri_density TYPE(lri_environment_type), POINTER :: lri_env @@ -586,6 +589,12 @@ SUBROUTINE qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces, just_energy, & END IF + ! set hartree and xc potentials for use in Harris method + IF (qs_env%harris_method) THEN + CALL get_qs_env(qs_env, harris_env=harris_env) + CALL harris_set_potentials(harris_env, v_hartree_rspace, v_rspace_new) + END IF + NULLIFY (rho_struct) IF (use_virial .AND. calculate_forces) THEN virial%pv_exc = virial%pv_exc - virial%pv_xc diff --git a/src/qs_nonscf.F b/src/qs_nonscf.F index 2072389e41..b688da168a 100644 --- a/src/qs_nonscf.F +++ b/src/qs_nonscf.F @@ -23,8 +23,10 @@ MODULE qs_nonscf section_vals_type USE kinds, ONLY: dp USE kpoint_types, ONLY: kpoint_type - USE machine, ONLY: m_walltime + USE machine, ONLY: m_flush,& + m_walltime USE message_passing, ONLY: mp_para_env_type + USE qs_charges_types, ONLY: qs_charges_type USE qs_core_energies, ONLY: calculate_ptrace USE qs_energy_types, ONLY: qs_energy_type USE qs_environment_types, ONLY: get_qs_env,& @@ -41,7 +43,6 @@ MODULE qs_nonscf USE qs_scf_loop_utils, ONLY: qs_scf_new_mos,& qs_scf_new_mos_kp USE qs_scf_output, ONLY: qs_scf_loop_print,& - qs_scf_print_summary,& qs_scf_write_mos USE qs_scf_post_scf, ONLY: qs_scf_compute_properties USE qs_scf_types, ONLY: qs_scf_env_type @@ -125,11 +126,11 @@ SUBROUTINE do_nonscf(qs_env, scf_env, scf_control) CHARACTER(LEN=*), PARAMETER :: routineN = 'do_nonscf' - INTEGER :: handle, img, ispin, output_unit + INTEGER :: handle, img, iounit, ispin LOGICAL :: diis_step, do_kpoints - REAL(KIND=dp) :: pc_ener, qmmm_el, t1, t2 + REAL(KIND=dp) :: pc_ener, qmmm_el, t1, t2, tdiag TYPE(cp_logger_type), POINTER :: logger - TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrixkp_ks, rho_ao_kp + TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrixkp_ks, rho_ao_kp TYPE(dft_control_type), POINTER :: dft_control TYPE(kpoint_type), POINTER :: kpoints TYPE(mo_set_type), DIMENSION(:), POINTER :: mos @@ -144,7 +145,7 @@ SUBROUTINE do_nonscf(qs_env, scf_env, scf_control) t1 = m_walltime() logger => cp_get_default_logger() - output_unit = cp_logger_get_default_io_unit(logger) + iounit = cp_logger_get_default_io_unit(logger) CALL get_qs_env(qs_env=qs_env, & energy=energy, & @@ -190,10 +191,17 @@ SUBROUTINE do_nonscf(qs_env, scf_env, scf_control) CALL dbcsr_copy(rho_ao_kp(ispin, img)%matrix, scf_env%p_mix_new(ispin, img)%matrix) END DO END DO + CALL qs_ks_did_change(ks_env, rho_changed=.TRUE., potential_changed=.TRUE.) - ! core energy : Tr(PH) + ! band energy : Tr(PH) CALL get_qs_env(qs_env, matrix_ks_kp=matrixkp_ks) + CALL calculate_ptrace(matrixkp_ks, rho_ao_kp, energy%band, dft_control%nspins) + ! core energy : Tr(Ph) + energy%total = energy%total - energy%core + CALL get_qs_env(qs_env, matrix_h_kp=matrix_h) + CALL calculate_ptrace(matrix_h, rho_ao_kp, energy%core, dft_control%nspins) + IF (qs_env%qmmm) THEN ! Compute QM/MM Energy CPASSERT(SIZE(matrixkp_ks, 2) == 1) @@ -207,22 +215,137 @@ SUBROUTINE do_nonscf(qs_env, scf_env, scf_control) ELSE energy%qmmm_el = 0.0_dp END IF - energy%total = energy%total - energy%core - CALL calculate_ptrace(matrixkp_ks, rho_ao_kp, energy%core, dft_control%nspins) - energy%total = energy%total + energy%core + energy%qmmm_el t2 = m_walltime() + tdiag = t2 - t1 - IF (output_unit > 0) THEN - WRITE (UNIT=output_unit, & - FMT="(T2,A,T40,A,F10.2,T61,F20.10)") & - "Diagonalization", "Time:", t2 - t1, energy%total - END IF - - CALL qs_scf_print_summary(output_unit, qs_env) + CALL qs_nonscf_print_summary(qs_env, tdiag, iounit) CALL timestop(handle) END SUBROUTINE do_nonscf +! ************************************************************************************************** +!> \brief writes a summary of information after diagonalization +!> \param qs_env ... +!> \param tdiag ... +!> \param iounit ... +! ************************************************************************************************** + SUBROUTINE qs_nonscf_print_summary(qs_env, tdiag, iounit) + TYPE(qs_environment_type), POINTER :: qs_env + REAL(KIND=dp), INTENT(IN) :: tdiag + INTEGER, INTENT(IN) :: iounit + + INTEGER :: nelectron_total + REAL(KIND=dp), DIMENSION(:), POINTER :: tot_rho_r + TYPE(dft_control_type), POINTER :: dft_control + TYPE(qs_charges_type), POINTER :: qs_charges + TYPE(qs_energy_type), POINTER :: energy + TYPE(qs_rho_type), POINTER :: rho + TYPE(qs_scf_env_type), POINTER :: scf_env + + IF (iounit > 0) THEN + CALL get_qs_env(qs_env=qs_env, energy=energy, dft_control=dft_control) + IF (qs_env%harris_method) THEN + CPASSERT(.NOT. dft_control%qs_control%gapw) + CPASSERT(.NOT. dft_control%qs_control%gapw_xc) + + CALL get_qs_env(qs_env=qs_env, rho=rho, scf_env=scf_env, qs_charges=qs_charges) + nelectron_total = scf_env%nelectron + CALL qs_rho_get(rho, tot_rho_r=tot_rho_r) + WRITE (UNIT=iounit, FMT="(/,(T3,A,T41,2F20.10))") & + "Electronic density on regular grids: ", & + SUM(tot_rho_r), & + SUM(tot_rho_r) + nelectron_total, & + "Core density on regular grids:", & + qs_charges%total_rho_core_rspace, & + qs_charges%total_rho_core_rspace - REAL(nelectron_total + dft_control%charge, dp) + WRITE (UNIT=iounit, FMT="(T3,A,T41,F20.10)") & + "Total charge density on r-space grids: ", & + SUM(tot_rho_r) + & + qs_charges%total_rho_core_rspace, & + "Total charge density g-space grids: ", & + qs_charges%total_rho_gspace + + WRITE (UNIT=iounit, FMT="(/,T2,A,T40,A,F10.2,T61,F20.10)") & + "Diagonalization", "Time:", tdiag, energy%band + + WRITE (UNIT=iounit, FMT="((T3,A,T56,F25.14))") & + "Core Hamiltonian energy: ", energy%core, & + "Overlap energy of the core charge distribution:", energy%core_overlap, & + "Self energy of the core charge distribution: ", energy%core_self, & + "Hartree energy: ", energy%hartree, & + "Exchange-correlation energy: ", energy%exc + IF (energy%dispersion /= 0.0_dp) & + WRITE (UNIT=iounit, FMT="(T3,A,T56,F25.14)") & + "Dispersion energy: ", energy%dispersion + IF (energy%gcp /= 0.0_dp) & + WRITE (UNIT=iounit, FMT="(T3,A,T56,F25.14)") & + "gCP energy: ", energy%gcp + IF (energy%efield /= 0.0_dp) & + WRITE (UNIT=iounit, FMT="(T3,A,T56,F25.14)") & + "Electric field interaction energy: ", energy%efield + + ELSEIF (dft_control%qs_control%semi_empirical) THEN + CPABORT("NONSCF not available") + ELSEIF (dft_control%qs_control%dftb) THEN + CPASSERT(energy%dftb3 == 0.0_dp) + energy%total = energy%total + energy%band + energy%qmmm_el + WRITE (UNIT=iounit, FMT="(/,T2,A,T40,A,F10.2,T61,F20.10)") & + "Diagonalization", "Time:", tdiag, energy%total + WRITE (UNIT=iounit, FMT="((T3,A,T56,F25.14))") & + "Core Hamiltonian energy: ", energy%core, & + "Repulsive potential energy: ", energy%repulsive, & + "Dispersion energy: ", energy%dispersion + IF (energy%efield /= 0.0_dp) & + WRITE (UNIT=iounit, FMT="(T3,A,T56,F25.14)") & + "Electric field interaction energy: ", energy%efield + ELSEIF (dft_control%qs_control%xtb) THEN + energy%total = energy%total + energy%band + energy%qmmm_el + WRITE (UNIT=iounit, FMT="(/,T2,A,T40,A,F10.2,T61,F20.10)") & + "Diagonalization", "Time:", tdiag, energy%total + WRITE (UNIT=iounit, FMT="((T3,A,T56,F25.14))") & + "Core Hamiltonian energy: ", energy%core, & + "Repulsive potential energy: ", energy%repulsive, & + "Electronic energy: ", energy%hartree, & + "Dispersion energy: ", energy%dispersion + IF (dft_control%qs_control%xtb_control%xb_interaction) & + WRITE (UNIT=iounit, FMT="(T3,A,T56,F25.14)") & + "Correction for halogen bonding: ", energy%xtb_xb_inter + IF (dft_control%qs_control%xtb_control%do_nonbonded) & + WRITE (UNIT=iounit, FMT="(T3,A,T56,F25.14)") & + "Correction for nonbonded interactions: ", energy%xtb_nonbonded + IF (energy%efield /= 0.0_dp) & + WRITE (UNIT=iounit, FMT="(T3,A,T56,F25.14)") & + "Electric field interaction energy: ", energy%efield + ELSE + CPABORT("NONSCF not available") + END IF + IF (dft_control%smear) THEN + WRITE (UNIT=iounit, FMT="((T3,A,T56,F25.14))") & + "Electronic entropic energy:", energy%kTS + WRITE (UNIT=iounit, FMT="((T3,A,T56,F25.14))") & + "Fermi energy:", energy%efermi + END IF + IF (energy%qmmm_el /= 0.0_dp) THEN + WRITE (UNIT=iounit, FMT="(T3,A,T56,F25.14)") & + "QM/MM Electrostatic energy: ", energy%qmmm_el + IF (qs_env%qmmm_env_qm%image_charge) THEN + WRITE (UNIT=iounit, FMT="(T3,A,T56,F25.14)") & + "QM/MM image charge energy: ", energy%image_charge + END IF + END IF + + IF (dft_control%qs_control%semi_empirical) THEN + CPABORT("NONSCF not available") + ELSEIF (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb) THEN + WRITE (UNIT=iounit, FMT="(/,(T3,A,T56,F25.14))") & + "Total energy: ", energy%total + END IF + + CALL m_flush(iounit) + END IF + + END SUBROUTINE qs_nonscf_print_summary + END MODULE qs_nonscf diff --git a/src/qs_rho_methods.F b/src/qs_rho_methods.F index 6d9d209557..9a74182dec 100644 --- a/src/qs_rho_methods.F +++ b/src/qs_rho_methods.F @@ -45,6 +45,8 @@ MODULE qs_rho_methods USE qs_environment_types, ONLY: get_qs_env,& qs_environment_type,& set_qs_env + USE qs_harris_types, ONLY: harris_type + USE qs_harris_utils, ONLY: calculate_harris_density USE qs_kind_types, ONLY: qs_kind_type USE qs_ks_types, ONLY: get_ks_env,& qs_ks_env_type @@ -386,6 +388,7 @@ SUBROUTINE qs_rho_update_rho(rho_struct, qs_env, rho_xc_external, local_rho_set, TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp TYPE(dft_control_type), POINTER :: dft_control + TYPE(harris_type), POINTER :: harris_env TYPE(kpoint_type), POINTER :: kpoints TYPE(lri_density_type), POINTER :: lri_density TYPE(lri_environment_type), POINTER :: lri_env @@ -397,8 +400,13 @@ SUBROUTINE qs_rho_update_rho(rho_struct, qs_env, rho_xc_external, local_rho_set, para_env=para_env) IF (PRESENT(para_env_external)) para_env => para_env_external - IF (dft_control%qs_control%semi_empirical .OR. & - dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb) THEN + IF (qs_env%harris_method) THEN + CALL get_qs_env(qs_env, harris_env=harris_env) + CALL calculate_harris_density(qs_env, harris_env%rhoin, rho_struct) + CALL qs_rho_set(rho_struct, rho_r_valid=.TRUE., rho_g_valid=.TRUE.) + + ELSEIF (dft_control%qs_control%semi_empirical .OR. & + dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb) THEN CALL qs_rho_set(rho_struct, rho_r_valid=.FALSE., rho_g_valid=.FALSE.) diff --git a/src/qs_scf_initialization.F b/src/qs_scf_initialization.F index a394c3982e..6e8c706546 100644 --- a/src/qs_scf_initialization.F +++ b/src/qs_scf_initialization.F @@ -73,6 +73,8 @@ MODULE qs_scf_initialization fb_env_write_info USE qs_fb_env_types, ONLY: fb_env_create,& fb_env_has_data + USE qs_harris_types, ONLY: harris_type + USE qs_harris_utils, ONLY: harris_density_update USE qs_initial_guess, ONLY: calculate_first_density_matrix USE qs_kind_types, ONLY: get_qs_kind,& qs_kind_type,& @@ -1045,10 +1047,11 @@ SUBROUTINE scf_env_initial_rho_setup(scf_env, qs_env, scf_section, scf_control) INTEGER :: extrapolation_method_nr, handle, ispin, & nmo, output_unit - LOGICAL :: orthogonal_wf + LOGICAL :: do_harris, orthogonal_wf TYPE(cp_fm_type), POINTER :: mo_coeff TYPE(cp_logger_type), POINTER :: logger TYPE(dft_control_type), POINTER :: dft_control + TYPE(harris_type), POINTER :: harris_env TYPE(mo_set_type), DIMENSION(:), POINTER :: mos TYPE(mp_para_env_type), POINTER :: para_env TYPE(qs_rho_type), POINTER :: rho @@ -1066,6 +1069,8 @@ SUBROUTINE scf_env_initial_rho_setup(scf_env, qs_env, scf_section, scf_control) dft_control=dft_control, & para_env=para_env) + do_harris = qs_env%harris_method + extrapolation_method_nr = wfi_use_guess_method_nr IF (ASSOCIATED(qs_env%wf_history)) THEN CALL wfi_extrapolate(qs_env%wf_history, & @@ -1084,22 +1089,30 @@ SUBROUTINE scf_env_initial_rho_setup(scf_env, qs_env, scf_section, scf_control) END DO END IF END IF - output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%PROGRAM_RUN_INFO", & - extension=".scfLog") - IF (output_unit > 0) THEN - WRITE (UNIT=output_unit, FMT="(/,T2,A,I0)") & - "Extrapolation method: "// & - TRIM(wfi_get_method_label(extrapolation_method_nr)) - IF (extrapolation_method_nr == wfi_ps_method_nr) THEN - WRITE (UNIT=output_unit, FMT="(T2,A,I0,A)") & - "Extrapolation order: ", & - MAX((MIN(qs_env%wf_history%memory_depth, qs_env%wf_history%snapshot_count) - 1), 0) + + IF (.NOT. do_harris) THEN + output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%PROGRAM_RUN_INFO", & + extension=".scfLog") + IF (output_unit > 0) THEN + WRITE (UNIT=output_unit, FMT="(/,T2,A,I0)") & + "Extrapolation method: "// & + TRIM(wfi_get_method_label(extrapolation_method_nr)) + IF (extrapolation_method_nr == wfi_ps_method_nr) THEN + WRITE (UNIT=output_unit, FMT="(T2,A,I0,A)") & + "Extrapolation order: ", & + MAX((MIN(qs_env%wf_history%memory_depth, qs_env%wf_history%snapshot_count) - 1), 0) + END IF END IF + CALL cp_print_key_finished_output(output_unit, logger, scf_section, & + "PRINT%PROGRAM_RUN_INFO") END IF - CALL cp_print_key_finished_output(output_unit, logger, scf_section, & - "PRINT%PROGRAM_RUN_INFO") - IF (extrapolation_method_nr == wfi_use_guess_method_nr) THEN + IF (do_harris) THEN + CALL get_qs_env(qs_env, harris_env=harris_env) + CALL harris_density_update(qs_env, harris_env) + CALL qs_rho_update_rho(rho, qs_env=qs_env) + CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.) + ELSE IF (extrapolation_method_nr == wfi_use_guess_method_nr) THEN CALL calculate_first_density_matrix(scf_env=scf_env, qs_env=qs_env) CALL qs_rho_update_rho(rho, qs_env=qs_env) CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.) diff --git a/tests/QS/regtest-harris/2H2O_t01.inp b/tests/QS/regtest-harris/2H2O_t01.inp new file mode 100644 index 0000000000..0d268ce7c6 --- /dev/null +++ b/tests/QS/regtest-harris/2H2O_t01.inp @@ -0,0 +1,58 @@ +&GLOBAL + PRINT_LEVEL LOW + PROJECT TEST1 + RUN_TYPE ENERGY +&END GLOBAL + +&FORCE_EVAL + METHOD Quickstep + &DFT + BASIS_SET_FILE_NAME BASIS_MOLOPT_UZH + POTENTIAL_FILE_NAME POTENTIAL_UZH + &HARRIS_METHOD ON + &END HARRIS_METHOD + &MGRID + CUTOFF 300 + REL_CUTOFF 60 + &END MGRID + &QS + EPS_DEFAULT 1.E-12 + &END QS + &SCF + SCF_GUESS NONE + &END SCF + &XC + &VDW_POTENTIAL + DISPERSION_FUNCTIONAL PAIR_POTENTIAL + &PAIR_POTENTIAL + PARAMETER_FILE_NAME dftd3.dat + REFERENCE_FUNCTIONAL BLYP + TYPE DFTD3(BJ) + &END PAIR_POTENTIAL + &END VDW_POTENTIAL + &XC_FUNCTIONAL BLYP + &END XC_FUNCTIONAL + &END XC + &END DFT + &SUBSYS + &CELL + ABC 8.0 8.0 8.0 + &END CELL + &COORD + O 0.000000 0.000000 0.000000 H2O1 + H 0.000000 0.000000 1.000000 H2O1 + H 0.942809 0.000000 -0.333333 H2O1 + O -1.617979 -0.948062 -2.341650 H2O2 + H -2.529195 -1.296822 -2.122437 H2O2 + H -1.534288 -0.833088 -3.331486 H2O2 + &END COORD + &KIND H + BASIS_SET ORB DZVP-MOLOPT-GGA-GTH-q1 + POTENTIAL GTH-GGA-q1 + &END KIND + &KIND O + BASIS_SET ORB DZVP-MOLOPT-GGA-GTH-q6 + POTENTIAL GTH-GGA-q6 + &END KIND + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/QS/regtest-harris/TEST_FILES b/tests/QS/regtest-harris/TEST_FILES new file mode 100644 index 0000000000..8c4582bb83 --- /dev/null +++ b/tests/QS/regtest-harris/TEST_FILES @@ -0,0 +1,6 @@ +# 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 +2H2O_t01.inp 11 1e-011 -34.570515851298339 +ft1.inp 0 +#EOF diff --git a/tests/QS/regtest-harris/ft1.inp b/tests/QS/regtest-harris/ft1.inp new file mode 100644 index 0000000000..85417d577e --- /dev/null +++ b/tests/QS/regtest-harris/ft1.inp @@ -0,0 +1,72 @@ +&GLOBAL + PRINT_LEVEL LOW + PROJECT cc02 + RUN_TYPE DEBUG +&END GLOBAL + +&DEBUG + CHECK_ATOM_FORCE 1 xyz + DE 0.001 + DEBUG_DIPOLE .FALSE. + DEBUG_FORCES .TRUE. + DEBUG_POLARIZABILITY .FALSE. + DEBUG_STRESS_TENSOR .FALSE. + EPS_NO_ERROR_CHECK 1.E-4 + MAX_RELATIVE_ERROR 0.05 + STOP_ON_MISMATCH T +&END DEBUG + +&FORCE_EVAL + METHOD Quickstep + &DFT + BASIS_SET_FILE_NAME BASIS_MOLOPT_UZH + POTENTIAL_FILE_NAME POTENTIAL_UZH + &HARRIS_METHOD ON + &END HARRIS_METHOD + &MGRID + CUTOFF 300 + REL_CUTOFF 60 + &END MGRID + &QS + EPS_DEFAULT 1.E-14 + &END QS + &SCF + SCF_GUESS NONE + &END SCF + &XC + &VDW_POTENTIAL + DISPERSION_FUNCTIONAL PAIR_POTENTIAL + &PAIR_POTENTIAL + PARAMETER_FILE_NAME dftd3.dat + REFERENCE_FUNCTIONAL PBE + TYPE DFTD3(BJ) + &END PAIR_POTENTIAL + &END VDW_POTENTIAL + &XC_FUNCTIONAL PBE + &END XC_FUNCTIONAL + &END XC + &END DFT + &SUBSYS + &CELL + ABC 6.0 6.0 6.0 + &END CELL + &COORD + O 0.094933 -0.000368 0.895642 + C -0.031077 -0.000121 -0.307326 + H -0.090437 0.947608 -0.895642 + H -0.094933 -0.947608 -0.895562 + &END COORD + &KIND H + BASIS_SET ORB DZVP-MOLOPT-GGA-GTH-q1 + POTENTIAL GTH-GGA-q1 + &END KIND + &KIND C + BASIS_SET ORB DZVP-MOLOPT-GGA-GTH-q4 + POTENTIAL GTH-GGA-q4 + &END KIND + &KIND O + BASIS_SET ORB DZVP-MOLOPT-GGA-GTH-q6 + POTENTIAL GTH-GGA-q6 + &END KIND + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/TEST_DIRS b/tests/TEST_DIRS index 532534fcf7..597f85a04b 100644 --- a/tests/TEST_DIRS +++ b/tests/TEST_DIRS @@ -346,6 +346,7 @@ QS/regtest-linearscaling QS/regtest-spgr spglib Fist/regtest-field QS/regtest-iao +QS/regtest-harris TAMC/regtest QS/regtest-fftw-wisdom fftw3 FE/regtest-2