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

Harris Method (#3665)

Browse source
This commit is contained in:
Juerg Hutter 2024-09-09 09:57:35 +02:00 • committed by GitHub
parent 17fb60b92f
commit f2bc97ddae
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
29 changed files with 1845 additions and 148 deletions
Download patch file Download diff file Expand all files Collapse all files

3
src/CMakeLists.txt
View file

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

1
src/aobasis/aux_basis_set.F
View file

@ -220,7 +220,6 @@ CONTAINS
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

5
src/aobasis/basis_set_container_types.F
View file

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

9
src/aobasis/basis_set_types.F
View file

@ -279,13 +279,15 @@ CONTAINS
!> \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 @@ CONTAINS
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

141
src/atom_fit.F
View file

@ -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 @@ CONTAINS
!> \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 @@ CONTAINS
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%state = 0
IF (iunit > 0) THEN
WRITE (iunit, '(/," POWELL| Start optimization procedure")')
END IF
n10 = ostate%maxfun/10
DO
IF (ostate%state == 2) THEN
CALL density_fit(density, atom, num_gto, x(1), x(2), co, ostate%f)
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
IF (ostate%state == -1) EXIT
ostate%state = 0
IF (iunit > 0) THEN
WRITE (iunit, '(/," POWELL| Start optimization procedure")')
END IF
n10 = ostate%maxfun/10
DO
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
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
END DO
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
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
ostate%state = 8
CALL powell_optimize(ostate%nvar, x, ostate)
CALL density_fit(density, atom, num_gto, x(1), x(2), co, ostate%f)
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 @@ CONTAINS
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 @@ CONTAINS
! **************************************************************************************************
!> \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 @@ CONTAINS
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)

49
src/atom_kind_orbitals.F
View file

@ -472,25 +472,27 @@ CONTAINS
!> \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 @@ CONTAINS
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 @@ CONTAINS
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,20 +728,20 @@ CONTAINS
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)
ELSE
CALL calculate_atom(atom, -1)
CALL atom_fit_density(atom, ngto, 0, -1, results=results)
END IF
CALL calculate_atom(atom, iw)
xx = results(1)
cc = results(2)
DO i = 1, ngto
density(i, 1) = xx*cc**i
density(i, 2) = results(2 + i)
END DO
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 atom_fit_density(atom, ngto, 0, iw, agto=density(:, 1), results=results)
density(1:ngto, 2) = results(1:ngto)
END IF
! clean up
CALL atom_int_release(integrals)

5
src/input_constants.F
View file

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

5
src/input_cp2k_dft.F
View file

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

100
src/input_cp2k_harris.F Normal file
View file

@ -0,0 +1,100 @@
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
! !
! 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

2
src/input_cp2k_scf.F
View file

@ -215,7 +215,7 @@ CONTAINS
"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)

9
src/pw_env_methods.F
View file

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

6
src/qs_energy.F
View file

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

1
src/qs_energy_types.F
View file

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

68
src/qs_environment.F
View file

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

37
src/qs_environment_types.F
View file

@ -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 @@ CONTAINS
!> \param lri_density ...
!> \param exstate_env ...
!> \param ec_env ...
!> \param harris_env ...
!> \param dispersion_env ...
!> \param gcp_env ...
!> \param vee ...
@ -494,7 +500,7 @@ CONTAINS
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 @@ CONTAINS
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 @@ CONTAINS
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 @@ CONTAINS
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 @@ CONTAINS
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 @@ CONTAINS
!> \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 @@ CONTAINS
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 @@ CONTAINS
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 @@ CONTAINS
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 @@ CONTAINS
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 @@ CONTAINS
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

6
src/qs_force.F
View file

@ -360,7 +360,11 @@ CONTAINS
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

255
src/qs_harris_types.F Normal file
View file

@ -0,0 +1,255 @@
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
! !
! 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

736
src/qs_harris_utils.F Normal file
View file

@ -0,0 +1,736 @@
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
! !
! 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

4
src/qs_integrate_potential.F
View file

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

185
src/qs_integrate_potential_single.F
View file

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

10
src/qs_interactions.F
View file

@ -101,8 +101,8 @@ CONTAINS
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 @@ CONTAINS
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 @@ CONTAINS
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, &

9
src/qs_ks_methods.F
View file

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

157
src/qs_nonscf.F
View file

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

12
src/qs_rho_methods.F
View file

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

41
src/qs_scf_initialization.F
View file

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

58
tests/QS/regtest-harris/2H2O_t01.inp Normal file
View file

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

6
tests/QS/regtest-harris/TEST_FILES Normal file
View file

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

72
tests/QS/regtest-harris/ft1.inp Normal file
View file

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

1
tests/TEST_DIRS
View file

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