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

Response properties and Harris functionals (#706)

Browse source 
* Energy Correction moved to DFT section, new regtests

Towards Harris Functional Forces

Harris force: CPKS equations

Harris Forces Debug Code

Harris functional + external field and dipole

Harris forces: refactoring, CPKS solver

* Response: HFX and ADMM-HFX

Linear response: HFX and ADMM

* Debug Linres regtest adjustments

* Refactoring of hfx_derivatives: from mp2 specific to response

* Harris: HFX+ADMM ground state (1. part)

Harris forces for ADMM (not debugged)

* Harris functional using linear scaling methods

* Harris: solver and MAO
This commit is contained in:
Juerg Hutter 2020-01-06 21:34:36 +01:00 • committed by GitHub
parent 63f2e15396
commit 8a724e3728
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
63 changed files with 6199 additions and 1534 deletions
Download patch file Download diff file Expand all files Collapse all files

8
data/BASIS_SET
View file

@ -3288,6 +3288,14 @@ C DZVP-GTH-BLYP
3 2 2 1 1
0.6000000000 1.0000000000
#
N DZV-GTH-BLYP
1
2 0 1 4 2 2
6.1514293421 0.1529382304 0.0000000000 -0.0952593918 0.0000000000
1.8080207176 -0.0415737876 0.0000000000 -0.2946953365 0.0000000000
0.5656159036 -0.7060780663 0.0000000000 -0.4740427661 0.0000000000
0.1593207770 -0.3734107042 1.0000000000 -0.3884930037 1.0000000000
#
N DZVP-GTH-BLYP
2
2 0 1 4 2 2

377
src/ec_efield_local.F Normal file
View file

@ -0,0 +1,377 @@
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 - 2020 CP2K developers group !
!--------------------------------------------------------------------------------------------------!
! **************************************************************************************************
!> \brief Calculates the energy contribution and the mo_derivative of
!> a static electric field (nonperiodic)
!> \par History
!> Adjusted from qs_efield_local
!> \author JGH (10.2019)
! **************************************************************************************************
MODULE ec_efield_local
USE ai_moments, ONLY: dipole_force
USE atomic_kind_types, ONLY: atomic_kind_type,&
get_atomic_kind,&
get_atomic_kind_set
USE basis_set_types, ONLY: gto_basis_set_p_type,&
gto_basis_set_type
USE cell_types, ONLY: cell_type,&
pbc
USE cp_control_types, ONLY: dft_control_type
USE cp_para_types, ONLY: cp_para_env_type
USE dbcsr_api, ONLY: dbcsr_add,&
dbcsr_copy,&
dbcsr_get_block_p,&
dbcsr_p_type,&
dbcsr_set
USE ec_env_types, ONLY: energy_correction_type
USE kinds, ONLY: dp
USE orbital_pointers, ONLY: ncoset
USE particle_types, ONLY: particle_type
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_kind_types, ONLY: get_qs_kind,&
qs_kind_type
USE qs_moments, ONLY: build_local_moment_matrix
USE qs_neighbor_list_types, ONLY: get_iterator_info,&
neighbor_list_iterate,&
neighbor_list_iterator_create,&
neighbor_list_iterator_p_type,&
neighbor_list_iterator_release,&
neighbor_list_set_p_type
USE qs_period_efield_types, ONLY: efield_berry_type,&
init_efield_matrices,&
set_efield_matrices
#include "./base/base_uses.f90"
IMPLICIT NONE
PRIVATE
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ec_efield_local'
! *** Public subroutines ***
PUBLIC :: ec_efield_local_operator, ec_efield_integrals
! **************************************************************************************************
CONTAINS
! **************************************************************************************************
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param ec_env ...
!> \param calculate_forces ...
! **************************************************************************************************
SUBROUTINE ec_efield_local_operator(qs_env, ec_env, calculate_forces)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(energy_correction_type), POINTER :: ec_env
LOGICAL, INTENT(IN) :: calculate_forces
CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_efield_local_operator', &
routineP = moduleN//':'//routineN
INTEGER :: handle
REAL(dp), DIMENSION(3) :: rpoint
TYPE(dft_control_type), POINTER :: dft_control
CALL timeset(routineN, handle)
NULLIFY (dft_control)
CALL get_qs_env(qs_env, dft_control=dft_control)
IF (dft_control%apply_efield) THEN
rpoint = 0.0_dp
CALL ec_efield_integrals(qs_env, ec_env, rpoint)
CALL ec_efield_mo_derivatives(qs_env, ec_env, rpoint, calculate_forces)
END IF
CALL timestop(handle)
END SUBROUTINE ec_efield_local_operator
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param ec_env ...
!> \param rpoint ...
! **************************************************************************************************
SUBROUTINE ec_efield_integrals(qs_env, ec_env, rpoint)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(energy_correction_type), POINTER :: ec_env
REAL(dp), DIMENSION(3), INTENT(IN) :: rpoint
CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_efield_integrals', &
routineP = moduleN//':'//routineN
INTEGER :: handle, i
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: dipmat, matrix_s
TYPE(dft_control_type), POINTER :: dft_control
TYPE(efield_berry_type), POINTER :: efield, efieldref
CALL timeset(routineN, handle)
CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, efield=efieldref)
efield => ec_env%efield
CALL init_efield_matrices(efield)
matrix_s => ec_env%matrix_s(:, 1)
ALLOCATE (dipmat(3))
DO i = 1, 3
ALLOCATE (dipmat(i)%matrix)
CALL dbcsr_copy(dipmat(i)%matrix, matrix_s(1)%matrix, 'DIP MAT')
CALL dbcsr_set(dipmat(i)%matrix, 0.0_dp)
END DO
CALL build_local_moment_matrix(qs_env, dipmat, 1, rpoint, basis_type="HARRIS")
CALL set_efield_matrices(efield=efield, dipmat=dipmat)
ec_env%efield => efield
CALL timestop(handle)
END SUBROUTINE ec_efield_integrals
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param ec_env ...
!> \param rpoint ...
!> \param calculate_forces ...
! **************************************************************************************************
SUBROUTINE ec_efield_mo_derivatives(qs_env, ec_env, rpoint, calculate_forces)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(energy_correction_type), POINTER :: ec_env
REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rpoint
LOGICAL :: calculate_forces
CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_efield_mo_derivatives', &
routineP = moduleN//':'//routineN
INTEGER :: atom_a, atom_b, handle, i, ia, iatom, icol, idir, ikind, irow, iset, ispin, &
jatom, jkind, jset, ldab, natom, ncoa, ncob, nkind, nseta, nsetb, sgfa, sgfb
INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind
INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
npgfb, nsgfa, nsgfb
INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
LOGICAL :: found, trans
REAL(dp) :: charge, dab, fdir
REAL(dp), DIMENSION(3) :: ci, fieldpol, ra, rab, rac, rbc, ria
REAL(dp), DIMENSION(3, 3) :: forcea, forceb
REAL(dp), DIMENSION(:, :), POINTER :: p_block_a, p_block_b, pblock, pmat, work
REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(cell_type), POINTER :: cell
TYPE(cp_para_env_type), POINTER :: para_env
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: dipmat, matrix_ks, matrix_p
TYPE(dft_control_type), POINTER :: dft_control
TYPE(efield_berry_type), POINTER :: efield
TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
TYPE(neighbor_list_iterator_p_type), &
DIMENSION(:), POINTER :: nl_iterator
TYPE(neighbor_list_set_p_type), DIMENSION(:), &
POINTER :: sab_orb
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(qs_energy_type), POINTER :: energy
TYPE(qs_force_type), DIMENSION(:), POINTER :: force
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(qs_kind_type), POINTER :: qs_kind
CALL timeset(routineN, handle)
CALL get_qs_env(qs_env, dft_control=dft_control, cell=cell, particle_set=particle_set)
CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, &
energy=energy, para_env=para_env, sab_orb=sab_orb)
efield => ec_env%efield
fieldpol = dft_control%efield_fields(1)%efield%polarisation* &
dft_control%efield_fields(1)%efield%strength
! nuclear contribution
natom = SIZE(particle_set)
IF (calculate_forces) THEN
CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, force=force)
ALLOCATE (atom_of_kind(natom))
CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind)
END IF
ci = 0.0_dp
DO ia = 1, natom
CALL get_atomic_kind(particle_set(ia)%atomic_kind, kind_number=ikind)
CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge)
ria = particle_set(ia)%r - rpoint
ria = pbc(ria, cell)
ci(:) = ci(:) + charge*ria(:)
IF (calculate_forces) THEN
IF (para_env%mepos == 0) THEN
iatom = atom_of_kind(ia)
DO idir = 1, 3
force(ikind)%efield(idir, iatom) = force(ikind)%efield(idir, iatom) - fieldpol(idir)*charge
END DO
END IF
END IF
END DO
IF (ec_env%should_update) THEN
ec_env%efield_nuclear = -SUM(ci(:)*fieldpol(:))
! Update KS matrix
matrix_ks => ec_env%matrix_h(:, 1)
dipmat => efield%dipmat
DO ispin = 1, SIZE(matrix_ks)
DO idir = 1, 3
CALL dbcsr_add(matrix_ks(ispin)%matrix, dipmat(idir)%matrix, &
alpha_scalar=1.0_dp, beta_scalar=fieldpol(idir))
END DO
END DO
END IF
! forces from the efield contribution
IF (calculate_forces) THEN
nkind = SIZE(qs_kind_set)
natom = SIZE(particle_set)
ALLOCATE (basis_set_list(nkind))
DO ikind = 1, nkind
qs_kind => qs_kind_set(ikind)
CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type="HARRIS")
IF (ASSOCIATED(basis_set_a)) THEN
basis_set_list(ikind)%gto_basis_set => basis_set_a
ELSE
NULLIFY (basis_set_list(ikind)%gto_basis_set)
END IF
END DO
!
CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
iatom=iatom, jatom=jatom, r=rab)
basis_set_a => basis_set_list(ikind)%gto_basis_set
IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
basis_set_b => basis_set_list(jkind)%gto_basis_set
IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
! basis ikind
first_sgfa => basis_set_a%first_sgf
la_max => basis_set_a%lmax
la_min => basis_set_a%lmin
npgfa => basis_set_a%npgf
nseta = basis_set_a%nset
nsgfa => basis_set_a%nsgf_set
rpgfa => basis_set_a%pgf_radius
set_radius_a => basis_set_a%set_radius
sphi_a => basis_set_a%sphi
zeta => basis_set_a%zet
! basis jkind
first_sgfb => basis_set_b%first_sgf
lb_max => basis_set_b%lmax
lb_min => basis_set_b%lmin
npgfb => basis_set_b%npgf
nsetb = basis_set_b%nset
nsgfb => basis_set_b%nsgf_set
rpgfb => basis_set_b%pgf_radius
set_radius_b => basis_set_b%set_radius
sphi_b => basis_set_b%sphi
zetb => basis_set_b%zet
atom_a = atom_of_kind(iatom)
atom_b = atom_of_kind(jatom)
ra(:) = particle_set(iatom)%r(:) - rpoint(:)
rac(:) = pbc(ra(:), cell)
rbc(:) = rac(:) + rab(:)
dab = SQRT(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3))
IF (iatom <= jatom) THEN
irow = iatom
icol = jatom
trans = .FALSE.
ELSE
irow = jatom
icol = iatom
trans = .TRUE.
END IF
fdir = 2.0_dp
IF (iatom == jatom .AND. dab < 1.e-10_dp) fdir = 1.0_dp
! density matrix
NULLIFY (p_block_a)
CALL dbcsr_get_block_p(ec_env%matrix_p(1, 1)%matrix, irow, icol, p_block_a, found)
IF (.NOT. found) CYCLE
IF (SIZE(matrix_p, 1) > 1) THEN
NULLIFY (p_block_b)
CALL dbcsr_get_block_p(ec_env%matrix_p(2, 1)%matrix, irow, icol, p_block_b, found)
CPASSERT(found)
END IF
forcea = 0.0_dp
forceb = 0.0_dp
DO iset = 1, nseta
ncoa = npgfa(iset)*ncoset(la_max(iset))
sgfa = first_sgfa(1, iset)
DO jset = 1, nsetb
IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
ncob = npgfb(jset)*ncoset(lb_max(jset))
sgfb = first_sgfb(1, jset)
! Calculate the primitive integrals (da|O|b) and (a|O|db)
ldab = MAX(ncoa, ncob)
ALLOCATE (work(ldab, ldab), pmat(ncoa, ncob))
! Decontract P matrix block
pmat = 0.0_dp
DO i = 1, SIZE(matrix_p, 1)
IF (i == 1) THEN
pblock => p_block_a
ELSE
pblock => p_block_b
END IF
IF (trans) THEN
CALL dgemm("N", "T", ncoa, nsgfb(jset), nsgfa(iset), &
1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
pblock(sgfb, sgfa), SIZE(pblock, 1), &
0.0_dp, work(1, 1), ldab)
ELSE
CALL dgemm("N", "N", ncoa, nsgfb(jset), nsgfa(iset), &
1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
pblock(sgfa, sgfb), SIZE(pblock, 1), &
0.0_dp, work(1, 1), ldab)
END IF
CALL dgemm("N", "T", ncoa, ncob, nsgfb(jset), &
1.0_dp, work(1, 1), ldab, &
sphi_b(1, sgfb), SIZE(sphi_b, 1), &
1.0_dp, pmat(1, 1), ncoa)
END DO
CALL dipole_force(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), &
1, rac, rbc, pmat, forcea, forceb)
DEALLOCATE (work, pmat)
END DO
END DO
DO idir = 1, 3
force(ikind)%efield(1:3, atom_a) = force(ikind)%efield(1:3, atom_a) &
+ fdir*fieldpol(idir)*forcea(idir, 1:3)
force(jkind)%efield(1:3, atom_b) = force(jkind)%efield(1:3, atom_b) &
+ fdir*fieldpol(idir)*forceb(idir, 1:3)
END DO
END DO
CALL neighbor_list_iterator_release(nl_iterator)
DEALLOCATE (basis_set_list)
DEALLOCATE (atom_of_kind)
END IF
CALL timestop(handle)
END SUBROUTINE ec_efield_mo_derivatives
END MODULE ec_efield_local

172
src/ec_env_types.F Normal file
View file

@ -0,0 +1,172 @@
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 - 2020 CP2K developers group !
!--------------------------------------------------------------------------------------------------!
! **************************************************************************************************
!> \brief Types needed for a for a Energy Correction
!> \par History
!> 2019.09 created
!> \author JGH
! **************************************************************************************************
MODULE ec_env_types
USE cp_dbcsr_operations, ONLY: dbcsr_deallocate_matrix_set
USE cp_fm_types, ONLY: cp_fm_p_type,&
cp_fm_release
USE dbcsr_api, ONLY: dbcsr_p_type
USE input_section_types, ONLY: section_vals_type
USE kinds, ONLY: dp
USE pw_types, ONLY: pw_p_type,&
pw_release
USE qs_dispersion_types, ONLY: qs_dispersion_release,&
qs_dispersion_type
USE qs_force_types, ONLY: deallocate_qs_force,&
qs_force_type
USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type,&
release_neighbor_list_sets
USE qs_p_env_types, ONLY: p_env_release,&
qs_p_env_type
USE qs_period_efield_types, ONLY: efield_berry_release,&
efield_berry_type
USE task_list_types, ONLY: deallocate_task_list,&
task_list_type
#include "./base/base_uses.f90"
IMPLICIT NONE
PRIVATE
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ec_env_types'
PUBLIC :: energy_correction_type, ec_env_release
! *****************************************************************************
!> \brief Contains information on the energy correction functional for KG
!> \par History
!> 03.2014 created
!> \author JGH
! *****************************************************************************
TYPE energy_correction_type
CHARACTER(len=20) :: ec_name
INTEGER :: energy_functional
INTEGER :: ks_solver
INTEGER :: factorization
REAL(KIND=dp) :: eps_default
LOGICAL :: should_update
! basis set
CHARACTER(len=20) :: basis
LOGICAL :: mao
INTEGER :: mao_max_iter
REAL(KIND=dp) :: mao_eps_grad
! energy components
REAL(KIND=dp) :: etotal
REAL(KIND=dp) :: eband, exc, ehartree, vhxc
REAL(KIND=dp) :: edispersion, efield_nuclear
! forces
TYPE(qs_force_type), DIMENSION(:), POINTER :: force => Null()
! full neighbor lists and corresponding task list
TYPE(neighbor_list_set_p_type), &
DIMENSION(:), POINTER :: sab_orb, sac_ppl, sap_ppnl
TYPE(task_list_type), POINTER :: task_list
! the XC function to be used for the correction, dispersion info
TYPE(section_vals_type), POINTER :: xc_section
TYPE(qs_dispersion_type), POINTER :: dispersion_env
! matrices in complete basis
! KS: Kohn-Sham; H: Core; S: overlap; T: kinetic energy;
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_t
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_w
! reduce basis
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mao_coef
! CP equations
TYPE(qs_p_env_type), POINTER :: p_env
TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: cpmos
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_hz
! potentials from input density
TYPE(pw_p_type), POINTER :: vh_rspace
TYPE(pw_p_type), DIMENSION(:), POINTER :: vxc_rspace, vtau_rspace
! efield
TYPE(efield_berry_type), POINTER :: efield => NULL()
END TYPE energy_correction_type
CONTAINS
! **************************************************************************************************
!> \brief ...
!> \param ec_env ...
! **************************************************************************************************
SUBROUTINE ec_env_release(ec_env)
TYPE(energy_correction_type), POINTER :: ec_env
CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_env_release', routineP = moduleN//':'//routineN
INTEGER :: handle, iab
CALL timeset(routineN, handle)
IF (ASSOCIATED(ec_env)) THEN
! neighbor lists
CALL release_neighbor_list_sets(ec_env%sab_orb)
CALL release_neighbor_list_sets(ec_env%sac_ppl)
CALL release_neighbor_list_sets(ec_env%sap_ppnl)
! forces
IF (ASSOCIATED(ec_env%force)) CALL deallocate_qs_force(ec_env%force)
! operator matrices
IF (ASSOCIATED(ec_env%matrix_ks)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_ks)
IF (ASSOCIATED(ec_env%matrix_h)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_h)
IF (ASSOCIATED(ec_env%matrix_s)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_s)
IF (ASSOCIATED(ec_env%matrix_t)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_t)
IF (ASSOCIATED(ec_env%matrix_p)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_p)
IF (ASSOCIATED(ec_env%matrix_w)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_w)
IF (ASSOCIATED(ec_env%task_list)) THEN
CALL deallocate_task_list(ec_env%task_list)
END IF
! reduced basis
IF (ASSOCIATED(ec_env%mao_coef)) CALL dbcsr_deallocate_matrix_set(ec_env%mao_coef)
! dispersion environment
IF (ASSOCIATED(ec_env%dispersion_env)) THEN
CALL qs_dispersion_release(ec_env%dispersion_env)
END IF
! CP env
IF (ASSOCIATED(ec_env%cpmos)) THEN
DO iab = 1, SIZE(ec_env%cpmos)
CALL cp_fm_release(ec_env%cpmos(iab)%matrix)
END DO
DEALLOCATE (ec_env%cpmos)
NULLIFY (ec_env%cpmos)
END IF
IF (ASSOCIATED(ec_env%matrix_hz)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_hz)
IF (ASSOCIATED(ec_env%p_env)) THEN
CALL p_env_release(ec_env%p_env)
END IF
! potential
IF (ASSOCIATED(ec_env%vh_rspace)) THEN
CALL pw_release(ec_env%vh_rspace%pw)
DEALLOCATE (ec_env%vh_rspace)
END IF
IF (ASSOCIATED(ec_env%vxc_rspace)) THEN
DO iab = 1, SIZE(ec_env%vxc_rspace)
CALL pw_release(ec_env%vxc_rspace(iab)%pw)
END DO
DEALLOCATE (ec_env%vxc_rspace)
END IF
IF (ASSOCIATED(ec_env%vtau_rspace)) THEN
DO iab = 1, SIZE(ec_env%vtau_rspace)
CALL pw_release(ec_env%vtau_rspace(iab)%pw)
END DO
DEALLOCATE (ec_env%vtau_rspace)
END IF
CALL efield_berry_release(ec_env%efield)
DEALLOCATE (ec_env)
END IF
CALL timestop(handle)
END SUBROUTINE ec_env_release
END MODULE ec_env_types

234
src/ec_environment.F Normal file
View file

@ -0,0 +1,234 @@
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 - 2020 CP2K developers group !
!--------------------------------------------------------------------------------------------------!
! **************************************************************************************************
!> \brief Energy correction environment setup and handling
!> \par History
!> 2019.09 created
!> \author JGH
! **************************************************************************************************
MODULE ec_environment
USE atomic_kind_types, ONLY: atomic_kind_type
USE basis_set_container_types, ONLY: add_basis_set_to_container,&
remove_basis_from_container
USE basis_set_types, ONLY: copy_gto_basis_set,&
create_primitive_basis_set,&
gto_basis_set_type
USE cp_control_types, ONLY: dft_control_type
USE cp_para_types, ONLY: cp_para_env_type
USE ec_env_types, ONLY: energy_correction_type
USE input_constants, ONLY: ec_functional_harris,&
xc_vdw_fun_nonloc,&
xc_vdw_fun_pairpot
USE input_cp2k_check, ONLY: xc_functionals_expand
USE input_section_types, ONLY: section_vals_get,&
section_vals_get_subs_vals,&
section_vals_type,&
section_vals_val_get
USE kinds, ONLY: dp
USE orbital_pointers, ONLY: init_orbital_pointers
USE qs_dispersion_nonloc, ONLY: qs_dispersion_nonloc_init
USE qs_dispersion_pairpot, ONLY: qs_dispersion_pairpot_init
USE qs_dispersion_types, ONLY: qs_dispersion_type
USE qs_dispersion_utils, ONLY: qs_dispersion_env_set
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
USE qs_interactions, ONLY: init_interaction_radii_orb_basis
USE qs_kind_types, ONLY: get_qs_kind,&
get_qs_kind_set,&
qs_kind_type
USE string_utilities, ONLY: uppercase
#include "./base/base_uses.f90"
IMPLICIT NONE
PRIVATE
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ec_environment'
PUBLIC :: ec_env_create
CONTAINS
! **************************************************************************************************
!> \brief Allocates and intitializes ec_env
!> \param qs_env ...
!> \param ec_env the object to create
!> \param dft_section ...
!> \par History
!> 2019.09 created
!> \author JGH
! **************************************************************************************************
SUBROUTINE ec_env_create(qs_env, ec_env, dft_section)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(energy_correction_type), POINTER :: ec_env
TYPE(section_vals_type), POINTER :: dft_section
CHARACTER(len=*), PARAMETER :: routineN = 'ec_env_create', routineP = moduleN//':'//routineN
CPASSERT(.NOT. ASSOCIATED(ec_env))
ALLOCATE (ec_env)
CALL init_ec_env(qs_env, ec_env, dft_section)
END SUBROUTINE ec_env_create
! **************************************************************************************************
!> \brief Initializes ec_env
!> \param qs_env ...
!> \param ec_env ...
!> \param dft_section ...
!> \par History
!> 2019.09 created
!> \author JGH
! **************************************************************************************************
SUBROUTINE init_ec_env(qs_env, ec_env, dft_section)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(energy_correction_type), POINTER :: ec_env
TYPE(section_vals_type), OPTIONAL, POINTER :: dft_section
CHARACTER(LEN=*), PARAMETER :: routineN = 'init_ec_env', routineP = moduleN//':'//routineN
INTEGER :: ikind, maxlgto, nkind
LOGICAL :: explicit
REAL(KIND=dp) :: eps_pgf_orb
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(cp_para_env_type), POINTER :: para_env
TYPE(dft_control_type), POINTER :: dft_control
TYPE(gto_basis_set_type), POINTER :: basis_set, harris_basis
TYPE(qs_dispersion_type), POINTER :: dispersion_env
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(qs_kind_type), POINTER :: qs_kind
TYPE(section_vals_type), POINTER :: nl_section, pp_section, section1, &
section2, xc_section
NULLIFY (atomic_kind_set, dispersion_env, para_env)
NULLIFY (ec_env%sab_orb, ec_env%sac_ppl, ec_env%sap_ppnl)
NULLIFY (ec_env%matrix_ks, ec_env%matrix_h, ec_env%matrix_s)
NULLIFY (ec_env%matrix_t, ec_env%matrix_p, ec_env%matrix_w)
NULLIFY (ec_env%task_list)
NULLIFY (ec_env%force)
NULLIFY (ec_env%mao_coef)
NULLIFY (ec_env%dispersion_env)
NULLIFY (ec_env%xc_section)
NULLIFY (ec_env%cpmos)
NULLIFY (ec_env%matrix_hz)
NULLIFY (ec_env%p_env)
NULLIFY (ec_env%vh_rspace)
NULLIFY (ec_env%vxc_rspace)
NULLIFY (ec_env%vtau_rspace)
ec_env%should_update = .TRUE.
ec_env%mao = .FALSE.
IF (qs_env%energy_correction) THEN
CALL section_vals_val_get(dft_section, "ENERGY_CORRECTION%ALGORITHM", &
i_val=ec_env%ks_solver)
CALL section_vals_val_get(dft_section, "ENERGY_CORRECTION%ENERGY_FUNCTIONAL", &
i_val=ec_env%energy_functional)
CALL section_vals_val_get(dft_section, "ENERGY_CORRECTION%FACTORIZATION", &
i_val=ec_env%factorization)
CALL section_vals_val_get(dft_section, "ENERGY_CORRECTION%EPS_DEFAULT", &
r_val=ec_env%eps_default)
CALL section_vals_val_get(dft_section, "ENERGY_CORRECTION%HARRIS_BASIS", &
c_val=ec_env%basis)
CALL section_vals_val_get(dft_section, "ENERGY_CORRECTION%MAO", &
l_val=ec_env%mao)
CALL section_vals_val_get(dft_section, "ENERGY_CORRECTION%MAO_MAX_ITER", &
i_val=ec_env%mao_max_iter)
CALL section_vals_val_get(dft_section, "ENERGY_CORRECTION%MAO_EPS_GRAD", &
r_val=ec_env%mao_eps_grad)
! set basis
CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, nkind=nkind)
CALL uppercase(ec_env%basis)
SELECT CASE (ec_env%basis)
CASE ("ORBITAL")
DO ikind = 1, nkind
qs_kind => qs_kind_set(ikind)
CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type="ORB")
IF (ASSOCIATED(basis_set)) THEN
NULLIFY (harris_basis)
CALL get_qs_kind(qs_kind=qs_kind, basis_set=harris_basis, basis_type="HARRIS")
IF (ASSOCIATED(harris_basis)) THEN
CALL remove_basis_from_container(qs_kind%basis_sets, basis_type="HARRIS")
END IF
NULLIFY (harris_basis)
CALL copy_gto_basis_set(basis_set, harris_basis)
CALL add_basis_set_to_container(qs_kind%basis_sets, harris_basis, "HARRIS")
END IF
END DO
CASE ("PRIMITIVE")
DO ikind = 1, nkind
qs_kind => qs_kind_set(ikind)
CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type="ORB")
IF (ASSOCIATED(basis_set)) THEN
NULLIFY (harris_basis)
CALL get_qs_kind(qs_kind=qs_kind, basis_set=harris_basis, basis_type="HARRIS")
IF (ASSOCIATED(harris_basis)) THEN
CALL remove_basis_from_container(qs_kind%basis_sets, basis_type="HARRIS")
END IF
NULLIFY (harris_basis)
CALL create_primitive_basis_set(basis_set, harris_basis)
CALL get_qs_env(qs_env, dft_control=dft_control)
eps_pgf_orb = dft_control%qs_control%eps_pgf_orb
CALL init_interaction_radii_orb_basis(harris_basis, eps_pgf_orb)
harris_basis%kind_radius = basis_set%kind_radius
CALL add_basis_set_to_container(qs_kind%basis_sets, harris_basis, "HARRIS")
END IF
END DO
CASE ("HARRIS")
DO ikind = 1, nkind
qs_kind => qs_kind_set(ikind)
NULLIFY (harris_basis)
CALL get_qs_kind(qs_kind=qs_kind, basis_set=harris_basis, basis_type="HARRIS")
IF (.NOT. ASSOCIATED(harris_basis)) THEN
CPWARN("Harris Basis not defined for all types of atoms.")
END IF
END DO
CASE DEFAULT
CPABORT("Unknown basis set for energy correction (Harris functional)")
END SELECT
!
CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto, basis_type="HARRIS")
CALL init_orbital_pointers(maxlgto + 1)
! set functional
SELECT CASE (ec_env%energy_functional)
CASE (ec_functional_harris)
ec_env%ec_name = "Harris"
CASE DEFAULT
CPABORT("unknown energy correction")
END SELECT
! select the XC section
NULLIFY (xc_section)
xc_section => section_vals_get_subs_vals(dft_section, "XC")
section1 => section_vals_get_subs_vals(dft_section, "ENERGY_CORRECTION%XC")
section2 => section_vals_get_subs_vals(dft_section, "ENERGY_CORRECTION%XC%XC_FUNCTIONAL")
CALL section_vals_get(section2, explicit=explicit)
IF (explicit) THEN
CALL xc_functionals_expand(section2, section1)
ec_env%xc_section => section1
ELSE
ec_env%xc_section => xc_section
END IF
! dispersion
ALLOCATE (dispersion_env)
NULLIFY (xc_section)
xc_section => ec_env%xc_section
CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, para_env=para_env)
CALL qs_dispersion_env_set(dispersion_env, xc_section)
IF (dispersion_env%type == xc_vdw_fun_pairpot) THEN
NULLIFY (pp_section)
pp_section => section_vals_get_subs_vals(xc_section, "VDW_POTENTIAL%PAIR_POTENTIAL")
CALL qs_dispersion_pairpot_init(atomic_kind_set, qs_kind_set, dispersion_env, pp_section, para_env)
ELSE IF (dispersion_env%type == xc_vdw_fun_nonloc) THEN
NULLIFY (nl_section)
nl_section => section_vals_get_subs_vals(xc_section, "VDW_POTENTIAL%NON_LOCAL")
CALL qs_dispersion_nonloc_init(dispersion_env, para_env)
END IF
ec_env%dispersion_env => dispersion_env
END IF
END SUBROUTINE init_ec_env
END MODULE ec_environment

1920
src/energy_corrections.F Normal file
View file

File diff suppressed because it is too large Load diff

63
src/hfx_admm_utils.F
View file

@ -24,6 +24,7 @@ MODULE hfx_admm_utils
USE cp_fm_types, ONLY: cp_fm_get_info,&
cp_fm_type
USE cp_log_handling, ONLY: cp_get_default_logger,&
cp_logger_get_default_io_unit,&
cp_logger_type
USE cp_output_handling, ONLY: cp_print_key_finished_output,&
cp_print_key_unit_nr
@ -223,10 +224,8 @@ CONTAINS
REAL(dp) :: eh1, ehfx, ehfxrt
REAL(dp), ALLOCATABLE, DIMENSION(:) :: hf_energy
TYPE(cp_para_env_type), POINTER :: para_env
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks_1d, matrix_ks_aux_fit, &
matrix_ks_aux_fit_hfx, &
matrix_ks_aux_fit_im, matrix_ks_im, &
rho_ao_1d
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks_1d, matrix_ks_aux_fit, &
matrix_ks_aux_fit_hfx, matrix_ks_aux_fit_im, matrix_ks_im, rho_ao_1d, rho_ao_resp
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_ks_orb, rho_ao_orb
TYPE(dft_control_type), POINTER :: dft_control
TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
@ -331,7 +330,12 @@ CONTAINS
IF (dft_control%do_admm) THEN
CALL scale_dm(qs_env, rho_ao_orb, scale_back=.FALSE.)
END IF
CALL derivatives_four_center(qs_env, rho_ao_orb, hfx_sections, &
IF (ASSOCIATED(qs_env%mp2_env)) THEN
CALL get_qs_env(qs_env, matrix_p_mp2=rho_ao_resp)
ELSE
NULLIFY (rho_ao_resp)
END IF
CALL derivatives_four_center(qs_env, rho_ao_orb, rho_ao_resp, hfx_sections, &
para_env, irep, use_virial)
!Scale auxiliary density matrix for ADMMP back with 1/gsi(ispin)
IF (dft_control%do_admm) THEN
@ -376,7 +380,8 @@ CONTAINS
END DO
IF (calculate_forces .AND. .NOT. do_adiabatic_rescaling) THEN
CALL derivatives_four_center(qs_env, rho_ao_orb, hfx_sections, &
NULLIFY (rho_ao_resp)
CALL derivatives_four_center(qs_env, rho_ao_orb, rho_ao_resp, hfx_sections, &
para_env, irep, use_virial)
END IF
@ -618,14 +623,17 @@ CONTAINS
CHARACTER(LEN=*), PARAMETER :: routineN = 'create_admm_xc_section', &
routineP = moduleN//':'//routineN
LOGICAL, PARAMETER :: debug_functional = .FALSE.
CHARACTER(LEN=20) :: name_x_func
INTEGER :: hfx_potential_type, ifun, nfun
INTEGER :: hfx_potential_type, ifun, iounit, nfun
LOGICAL :: funct_found
REAL(dp) :: cutoff_radius, hfx_fraction, omega, &
scale_x
TYPE(cp_logger_type), POINTER :: logger
TYPE(section_vals_type), POINTER :: xc_fun, xc_fun_section
logger => cp_get_default_logger()
NULLIFY (admm_env%xc_section_aux, admm_env%xc_section_primary)
CALL get_qs_env(qs_env)
@ -905,8 +913,11 @@ CONTAINS
END IF
IF (1 == 0) THEN
WRITE (*, *) "primary"
IF (debug_functional) THEN
iounit = cp_logger_get_default_io_unit(logger)
IF (iounit > 0) THEN
WRITE (iounit, "(A)") " ADMM Primary Basis Set Functional"
END IF
xc_fun_section => section_vals_get_subs_vals(admm_env%xc_section_primary, "XC_FUNCTIONAL")
ifun = 0
funct_found = .FALSE.
@ -917,20 +928,23 @@ CONTAINS
scale_x = -1000.0_dp
IF (xc_fun%section%name /= "LYP" .AND. xc_fun%section%name /= "VWN") THEN
CALL section_vals_val_get(xc_fun, "SCALE_X", &
r_val=scale_x)
CALL section_vals_val_get(xc_fun, "SCALE_X", r_val=scale_x)
END IF
IF (xc_fun%section%name == "XWPBE") THEN
CALL section_vals_val_get(xc_fun, "SCALE_X0", &
r_val=hfx_fraction)
WRITE (*, *) xc_fun%section%name, scale_x, hfx_fraction
CALL section_vals_val_get(xc_fun, "SCALE_X0", r_val=hfx_fraction)
IF (iounit > 0) THEN
WRITE (iounit, "(T5,A,T25,2F10.3)") TRIM(xc_fun%section%name), scale_x, hfx_fraction
END IF
ELSE
WRITE (*, *) xc_fun%section%name, scale_x
IF (iounit > 0) THEN
WRITE (iounit, "(T5,A,T25,F10.3)") TRIM(xc_fun%section%name), scale_x
END IF
END IF
END DO
WRITE (*, *) "auxiliary"
IF (iounit > 0) THEN
WRITE (iounit, "(A)") " Auxiliary Basis Set Functional"
END IF
xc_fun_section => section_vals_get_subs_vals(admm_env%xc_section_aux, "XC_FUNCTIONAL")
ifun = 0
funct_found = .FALSE.
@ -940,16 +954,17 @@ CONTAINS
IF (.NOT. ASSOCIATED(xc_fun)) EXIT
scale_x = -1000.0_dp
IF (xc_fun%section%name /= "LYP" .AND. xc_fun%section%name /= "VWN") THEN
CALL section_vals_val_get(xc_fun, "SCALE_X", &
r_val=scale_x)
CALL section_vals_val_get(xc_fun, "SCALE_X", r_val=scale_x)
END IF
IF (xc_fun%section%name == "XWPBE") THEN
CALL section_vals_val_get(xc_fun, "SCALE_X0", &
r_val=hfx_fraction)
WRITE (*, *) xc_fun%section%name, scale_x, hfx_fraction
CALL section_vals_val_get(xc_fun, "SCALE_X0", r_val=hfx_fraction)
IF (iounit > 0) THEN
WRITE (iounit, "(T5,A,T25,2F10.3)") TRIM(xc_fun%section%name), scale_x, hfx_fraction
END IF
ELSE
WRITE (*, *) xc_fun%section%name, scale_x
IF (iounit > 0) THEN
WRITE (iounit, "(T5,A,T25,F10.3)") TRIM(xc_fun%section%name), scale_x
END IF
END IF
END DO
END IF

306
src/hfx_derivatives.F
View file

@ -94,11 +94,13 @@ CONTAINS
!> forces%fock_4c arrays. Uses all 8 eri symmetries
!> \param qs_env ...
!> \param rho_ao density matrix
!> \param rho_ao_resp relaxed density matrix from response
!> \param hfx_section HFX input section
!> \param para_env para_env
!> \param irep ID of HFX replica
!> \param use_virial ...
!> \param adiabatic_rescale_factor parameter used for MCY3 hybrid
!> \param resp_only Calculate only forces from response density
!> \par History
!> 06.2007 created [Manuel Guidon]
!> 08.2007 optimized load balance [Manuel Guidon]
@ -108,16 +110,18 @@ CONTAINS
!> [Tobias Binninger + Valery Weber]
!> \author Manuel Guidon
! **************************************************************************************************
SUBROUTINE derivatives_four_center(qs_env, rho_ao, hfx_section, para_env, &
irep, use_virial, adiabatic_rescale_factor)
SUBROUTINE derivatives_four_center(qs_env, rho_ao, rho_ao_resp, hfx_section, para_env, &
irep, use_virial, adiabatic_rescale_factor, resp_only)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao_resp
TYPE(section_vals_type), POINTER :: hfx_section
TYPE(cp_para_env_type), POINTER :: para_env
INTEGER :: irep
LOGICAL :: use_virial
INTEGER, INTENT(IN) :: irep
LOGICAL, INTENT(IN) :: use_virial
REAL(dp), INTENT(IN), OPTIONAL :: adiabatic_rescale_factor
LOGICAL, INTENT(IN), OPTIONAL :: resp_only
CHARACTER(LEN=*), PARAMETER :: routineN = 'derivatives_four_center', &
routineP = moduleN//':'//routineN
@ -152,8 +156,8 @@ CONTAINS
INTEGER, DIMENSION(:, :, :, :), POINTER :: shm_set_offset
INTEGER, SAVE :: shm_number_of_p_entries
LOGICAL :: bins_left, buffer_overflow, do_dynamic_load_balancing, do_it, do_periodic, &
do_print_load_balance_info, is_anti_symmetric, screen_pmat_forces, treat_forces_in_core, &
use_disk_storage, with_mp2_density
do_print_load_balance_info, is_anti_symmetric, my_resp_only, screen_pmat_forces, &
treat_forces_in_core, use_disk_storage, with_resp_density
LOGICAL, DIMENSION(:, :), POINTER :: shm_atomic_pair_list
REAL(dp) :: bintime_start, bintime_stop, cartesian_estimate, cartesian_estimate_virial, &
compression_factor, eps_schwarz, eps_storage, fac, hf_fraction, ln_10, log10_eps_schwarz, &
@ -163,12 +167,12 @@ CONTAINS
tmp_virial(3, 3)
REAL(dp), ALLOCATABLE, DIMENSION(:) :: ede_buffer1, ede_buffer2, ede_primitives_tmp, &
ede_primitives_tmp_virial, ede_work, ede_work2, ede_work2_virial, ede_work_forces, &
ede_work_virial, pac_buf, pac_buf_beta, pac_buf_mp2, pac_buf_mp2_beta, pad_buf, &
pad_buf_beta, pad_buf_mp2, pad_buf_mp2_beta, pbc_buf, pbc_buf_beta, pbc_buf_mp2, &
pbc_buf_mp2_beta, pbd_buf, pbd_buf_beta, pbd_buf_mp2, pbd_buf_mp2_beta
ede_work_virial, pac_buf, pac_buf_beta, pac_buf_resp, pac_buf_resp_beta, pad_buf, &
pad_buf_beta, pad_buf_resp, pad_buf_resp_beta, pbc_buf, pbc_buf_beta, pbc_buf_resp, &
pbc_buf_resp_beta, pbd_buf, pbd_buf_beta, pbd_buf_resp, pbd_buf_resp_beta
REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: primitive_forces, primitive_forces_virial
REAL(dp), DIMENSION(:), POINTER :: full_density_mp2, full_density_mp2_beta, &
T2
REAL(dp), DIMENSION(:), POINTER :: full_density_resp, &
full_density_resp_beta, T2
REAL(dp), DIMENSION(:, :), POINTER :: full_density_alpha, full_density_beta, &
max_contraction, ptr_p_1, ptr_p_2, ptr_p_3, ptr_p_4, shm_pmax_atom, shm_pmax_block, &
sphi_b, zeta, zetb, zetc, zetd
@ -182,7 +186,6 @@ CONTAINS
TYPE(cell_type), POINTER :: cell
TYPE(cp_libint_t) :: private_deriv
TYPE(cp_logger_type), POINTER :: logger
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao_mp2
TYPE(dft_control_type), POINTER :: dft_control
TYPE(hfx_basis_info_type), POINTER :: basis_info
TYPE(hfx_basis_type), DIMENSION(:), POINTER :: basis_parameter
@ -218,10 +221,12 @@ CONTAINS
NULLIFY (dft_control, admm_env)
with_mp2_density = .FALSE.
IF (ASSOCIATED(qs_env%mp2_env)) with_mp2_density = .TRUE.
with_resp_density = .FALSE.
IF (ASSOCIATED(rho_ao_resp)) with_resp_density = .TRUE.
my_resp_only = .FALSE.
IF (PRESENT(resp_only)) my_resp_only = resp_only
IF (with_mp2_density) CALL get_qs_env(qs_env, matrix_p_mp2=rho_ao_mp2)
is_anti_symmetric = dbcsr_get_matrix_type(rho_ao(1, 1)%matrix) .EQ. dbcsr_type_antisymmetric
CALL get_qs_env(qs_env, &
@ -280,18 +285,19 @@ CONTAINS
!$OMP PARALLEL DEFAULT(NONE) SHARED(qs_env,&
!$OMP rho_ao,&
!$OMP rho_ao_mp2,&
!$OMP with_mp2_density,&
!$OMP rho_ao_resp,&
!$OMP with_resp_density,&
!$OMP hfx_section,&
!$OMP para_env,&
!$OMP irep,&
!$OMP my_resp_only,&
!$OMP ncoset,&
!$OMP nco,&
!$OMP nso,&
!$OMP n_threads,&
!$OMP full_density_alpha,&
!$OMP full_density_mp2,&
!$OMP full_density_mp2_beta,&
!$OMP full_density_resp,&
!$OMP full_density_resp_beta,&
!$OMP full_density_beta,&
!$OMP shm_initial_p,&
!$OMP shm_is_assoc_atomic_block,&
@ -347,9 +353,9 @@ CONTAINS
!$OMP neris_total,nimages,nints,nneighbors,npgfa,npgfb,npgfc,npgfd,&
!$OMP nprim_ints,n_processes,nseta,nsetb,nsgfa,nsgfb,nsgfc,nsgfd,&
!$OMP nsgfl_a,nsgfl_b,nsgfl_c,nsgfl_d,nsgf_max,offset_ac_set,offset_ad_set,&
!$OMP offset_bc_set,offset_bd_set,pac_buf,pac_buf_beta,pac_buf_mp2,pac_buf_mp2_beta,pad_buf,pad_buf_beta,pad_buf_mp2,&
!$OMP pad_buf_mp2_beta,particle_set,pbc_buf,pbc_buf_beta,pbc_buf_mp2,pbc_buf_mp2_beta,pbd_buf,pbd_buf_beta, &
!$OMP pbd_buf_mp2, pbd_buf_mp2_beta, pgf_list_ij,&
!$OMP offset_bc_set,offset_bd_set,pac_buf,pac_buf_beta,pac_buf_resp,pac_buf_resp_beta,pad_buf,pad_buf_beta,pad_buf_resp,&
!$OMP pad_buf_resp_beta,particle_set,pbc_buf,pbc_buf_beta,pbc_buf_resp,pbc_buf_resp_beta,pbd_buf,pbd_buf_beta, &
!$OMP pbd_buf_resp, pbd_buf_resp_beta, pgf_list_ij,&
!$OMP pgf_list_kl,pgf_product_list,pmax_atom,pmax_blocks,pmax_entry,pmax_tmp,potential_parameter,primitive_forces,&
!$OMP primitive_forces_virial,private_deriv,ptr_p_1,ptr_p_2,ptr_p_3,ptr_p_4,ra,rab2,&
!$OMP rb,rc,rcd2,rd,screening_parameter,screen_pmat_forces,set_list_ij,set_list_kl,&
@ -376,7 +382,7 @@ CONTAINS
load_balance_parameter => actual_x_data%load_balance_parameter
basis_parameter => actual_x_data%basis_parameter
! IF(with_mp2_density) THEN
! IF(with_resp_density) THEN
! ! here we also do a copy of the original load_balance_parameter
! ! since if MP2 forces are required after the calculation of the HFX
! ! derivatives we need to calculate again the SCF energy.
@ -491,10 +497,10 @@ CONTAINS
kind_of, basis_parameter, get_max_vals_spin=.FALSE., &
antisymmetric=is_anti_symmetric)
! for now only closed shell case
IF (with_mp2_density) THEN
NULLIFY (full_density_mp2)
ALLOCATE (full_density_mp2(shm_block_offset(ncpu + 1)))
CALL get_full_density(para_env, full_density_mp2, rho_ao_mp2(1)%matrix, shm_number_of_p_entries, &
IF (with_resp_density) THEN
NULLIFY (full_density_resp)
ALLOCATE (full_density_resp(shm_block_offset(ncpu + 1)))
CALL get_full_density(para_env, full_density_resp, rho_ao_resp(1)%matrix, shm_number_of_p_entries, &
shm_block_offset, &
kind_of, basis_parameter, get_max_vals_spin=.FALSE., &
antisymmetric=is_anti_symmetric)
@ -506,11 +512,11 @@ CONTAINS
shm_block_offset, &
kind_of, basis_parameter, get_max_vals_spin=.FALSE., &
antisymmetric=is_anti_symmetric)
! With mp2 density
IF (with_mp2_density) THEN
NULLIFY (full_density_mp2_beta)
ALLOCATE (full_density_mp2_beta(shm_block_offset(ncpu + 1)))
CALL get_full_density(para_env, full_density_mp2_beta, rho_ao_mp2(2)%matrix, shm_number_of_p_entries, &
! With resp density
IF (with_resp_density) THEN
NULLIFY (full_density_resp_beta)
ALLOCATE (full_density_resp_beta(shm_block_offset(ncpu + 1)))
CALL get_full_density(para_env, full_density_resp_beta, rho_ao_resp(2)%matrix, shm_number_of_p_entries, &
shm_block_offset, &
kind_of, basis_parameter, get_max_vals_spin=.FALSE., &
antisymmetric=is_anti_symmetric)
@ -544,13 +550,13 @@ CONTAINS
! restore as full density the HF density
! maybe in the future
IF (with_mp2_density) THEN
full_density_alpha(:, 1) = full_density_alpha(:, 1) - full_density_mp2
IF (with_resp_density) THEN
full_density_alpha(:, 1) = full_density_alpha(:, 1) - full_density_resp
IF (nspins == 2) THEN
full_density_beta(:, 1) = &
full_density_beta(:, 1) - full_density_mp2_beta
full_density_beta(:, 1) - full_density_resp_beta
ENDIF
! full_density_mp2=full_density+full_density_mp2
! full_density_resp=full_density+full_density_resp
END IF
screen_coeffs_set => actual_x_data%screen_funct_coeffs_set
@ -623,22 +629,22 @@ CONTAINS
ALLOCATE (pbc_buf(nsgf_max**2))
ALLOCATE (pad_buf(nsgf_max**2))
ALLOCATE (pac_buf(nsgf_max**2))
IF (with_mp2_density) THEN
ALLOCATE (pbd_buf_mp2(nsgf_max**2))
ALLOCATE (pbc_buf_mp2(nsgf_max**2))
ALLOCATE (pad_buf_mp2(nsgf_max**2))
ALLOCATE (pac_buf_mp2(nsgf_max**2))
IF (with_resp_density) THEN
ALLOCATE (pbd_buf_resp(nsgf_max**2))
ALLOCATE (pbc_buf_resp(nsgf_max**2))
ALLOCATE (pad_buf_resp(nsgf_max**2))
ALLOCATE (pac_buf_resp(nsgf_max**2))
END IF
IF (nspins == 2) THEN
ALLOCATE (pbd_buf_beta(nsgf_max**2))
ALLOCATE (pbc_buf_beta(nsgf_max**2))
ALLOCATE (pad_buf_beta(nsgf_max**2))
ALLOCATE (pac_buf_beta(nsgf_max**2))
IF (with_mp2_density) THEN
ALLOCATE (pbd_buf_mp2_beta(nsgf_max**2))
ALLOCATE (pbc_buf_mp2_beta(nsgf_max**2))
ALLOCATE (pad_buf_mp2_beta(nsgf_max**2))
ALLOCATE (pac_buf_mp2_beta(nsgf_max**2))
IF (with_resp_density) THEN
ALLOCATE (pbd_buf_resp_beta(nsgf_max**2))
ALLOCATE (pbc_buf_resp_beta(nsgf_max**2))
ALLOCATE (pad_buf_resp_beta(nsgf_max**2))
ALLOCATE (pac_buf_resp_beta(nsgf_max**2))
END IF
END IF
@ -1035,8 +1041,10 @@ CONTAINS
symm_fac = 0.25_dp
IF (iatom == jatom) symm_fac = symm_fac*2.0_dp
IF (katom == latom) symm_fac = symm_fac*2.0_dp
IF (iatom == katom .AND. jatom == latom .AND. iatom /= jatom .AND. katom /= latom) symm_fac = symm_fac*2.0_dp
IF (iatom == katom .AND. iatom == jatom .AND. katom == latom) symm_fac = symm_fac*2.0_dp
IF (iatom == katom .AND. jatom == latom .AND. &
iatom /= jatom .AND. katom /= latom) symm_fac = symm_fac*2.0_dp
IF (iatom == katom .AND. iatom == jatom .AND. &
katom == latom) symm_fac = symm_fac*2.0_dp
symm_fac = 1.0_dp/symm_fac
fac = fac*symm_fac
@ -1273,7 +1281,8 @@ CONTAINS
tmp_screen_pgf1 => screen_coeffs_pgf(:, :, jset, iset, jkind, ikind)
tmp_screen_pgf2 => screen_coeffs_pgf(:, :, lset, kset, lkind, kkind)
CALL forces4(private_deriv, ra, rb, rc, rd, npgfa(iset), npgfb(jset), npgfc(kset), npgfd(lset), &
CALL forces4(private_deriv, ra, rb, rc, rd, npgfa(iset), npgfb(jset), &
npgfc(kset), npgfd(lset), &
la_min(iset), la_max(iset), lb_min(jset), lb_max(jset), &
lc_min(kset), lc_max(kset), ld_min(lset), ld_max(lset), &
nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), &
@ -1285,7 +1294,8 @@ CONTAINS
zetc(1:npgfc(kset), kset), zetd(1:npgfd(lset), lset), &
primitive_forces, &
potential_parameter, &
actual_x_data%neighbor_cells, screen_coeffs_set(jset, iset, jkind, ikind)%x, &
actual_x_data%neighbor_cells, &
screen_coeffs_set(jset, iset, jkind, ikind)%x, &
screen_coeffs_set(lset, kset, lkind, kkind)%x, eps_schwarz, &
max_contraction_val, cartesian_estimate, cell, neris_tmp, &
log10_pmax, log10_eps_schwarz, &
@ -1297,7 +1307,8 @@ CONTAINS
ede_work, ede_work2, ede_work_forces, &
ede_buffer1, ede_buffer2, ede_primitives_tmp, &
nimages, do_periodic, use_virial, ede_work_virial, ede_work2_virial, &
ede_primitives_tmp_virial, primitive_forces_virial, cartesian_estimate_virial)
ede_primitives_tmp_virial, primitive_forces_virial, &
cartesian_estimate_virial)
nints = nsgfa(iset)*nsgfb(jset)*nsgfc(kset)*nsgfd(lset)*12
neris_total = neris_total + nints
@ -1379,7 +1390,7 @@ CONTAINS
END IF
END IF
END IF
IF (with_mp2_density) THEN
IF (with_resp_density) THEN
CALL prefetch_density_matrix(nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), &
full_density_alpha(:, 1), pbd_buf, pbc_buf, pad_buf, pac_buf, &
iatom, jatom, katom, latom, &
@ -1387,7 +1398,7 @@ CONTAINS
offset_ad_set, offset_ac_set, atomic_offset_bd, atomic_offset_bc, &
atomic_offset_ad, atomic_offset_ac, is_anti_symmetric)
CALL prefetch_density_matrix(nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), &
full_density_mp2, pbd_buf_mp2, pbc_buf_mp2, pad_buf_mp2, pac_buf_mp2, &
full_density_resp, pbd_buf_resp, pbc_buf_resp, pad_buf_resp, pac_buf_resp, &
iatom, jatom, katom, latom, &
iset, jset, kset, lset, offset_bd_set, offset_bc_set, &
offset_ad_set, offset_ac_set, atomic_offset_bd, atomic_offset_bc, &
@ -1401,7 +1412,7 @@ CONTAINS
atomic_offset_ad, atomic_offset_ac, is_anti_symmetric)
END IF
IF (nspins == 2) THEN
IF (with_mp2_density) THEN
IF (with_resp_density) THEN
CALL prefetch_density_matrix(nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), &
full_density_beta(:, 1), pbd_buf_beta, pbc_buf_beta, &
pad_buf_beta, pac_buf_beta, iatom, jatom, katom, latom, &
@ -1409,8 +1420,8 @@ CONTAINS
offset_ad_set, offset_ac_set, atomic_offset_bd, atomic_offset_bc, &
atomic_offset_ad, atomic_offset_ac, is_anti_symmetric)
CALL prefetch_density_matrix(nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), &
full_density_mp2_beta, pbd_buf_mp2_beta, pbc_buf_mp2_beta, &
pad_buf_mp2_beta, pac_buf_mp2_beta, iatom, jatom, katom, latom, &
full_density_resp_beta, pbd_buf_resp_beta, pbc_buf_resp_beta, &
pad_buf_resp_beta, pac_buf_resp_beta, iatom, jatom, katom, latom, &
iset, jset, kset, lset, offset_bd_set, offset_bc_set, &
offset_ad_set, offset_ac_set, atomic_offset_bd, atomic_offset_bc, &
atomic_offset_ad, atomic_offset_ac, is_anti_symmetric)
@ -1428,23 +1439,25 @@ CONTAINS
T2 => primitive_forces((coord - 1)*nsgfa(iset)*nsgfb(jset)*nsgfc(kset)*nsgfd(lset) + 1: &
coord*nsgfa(iset)*nsgfb(jset)*nsgfc(kset)*nsgfd(lset))
IF (with_mp2_density) THEN
IF (with_resp_density) THEN
CALL update_forces(nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), &
pbd_buf, pbc_buf, pad_buf, pac_buf, fac, &
T2, force, forces_map, coord, &
pbd_buf_mp2, pbc_buf_mp2, pad_buf_mp2, pac_buf_mp2)
pbd_buf_resp, pbc_buf_resp, pad_buf_resp, pac_buf_resp, &
my_resp_only)
ELSE
CALL update_forces(nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), &
pbd_buf, pbc_buf, pad_buf, pac_buf, fac, &
T2, force, forces_map, coord)
END IF
IF (nspins == 2) THEN
IF (with_mp2_density) THEN
IF (with_resp_density) THEN
CALL update_forces(nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), &
pbd_buf_beta, pbc_buf_beta, pad_buf_beta, pac_buf_beta, &
fac, T2, force, forces_map, coord, &
pbd_buf_mp2_beta, pbc_buf_mp2_beta, &
pad_buf_mp2_beta, pac_buf_mp2_beta)
pbd_buf_resp_beta, pbc_buf_resp_beta, &
pad_buf_resp_beta, pac_buf_resp_beta, &
my_resp_only)
ELSE
CALL update_forces(nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), &
pbd_buf_beta, pbc_buf_beta, pad_buf_beta, pac_buf_beta, fac, &
@ -1459,26 +1472,27 @@ CONTAINS
T2 => primitive_forces_virial( &
((i - 1)*12 + coord - 1)*nsgfa(iset)*nsgfb(jset)*nsgfc(kset)*nsgfd(lset) + 1: &
((i - 1)*12 + coord)*nsgfa(iset)*nsgfb(jset)*nsgfc(kset)*nsgfd(lset))
IF (with_mp2_density) THEN
IF (with_resp_density) THEN
CALL update_virial(nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), &
pbd_buf, pbc_buf, pad_buf, pac_buf, fac, &
T2, tmp_virial, coord, i, &
pbd_buf_mp2, pbc_buf_mp2, pad_buf_mp2, pac_buf_mp2)
pbd_buf_resp, pbc_buf_resp, pad_buf_resp, pac_buf_resp)
ELSE
CALL update_virial(nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), &
pbd_buf, pbc_buf, pad_buf, pac_buf, fac, &
T2, tmp_virial, coord, i)
END IF
IF (nspins == 2) THEN
IF (with_mp2_density) THEN
IF (with_resp_density) THEN
CALL update_virial(nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), &
pbd_buf_beta, pbc_buf_beta, pad_buf_beta, pac_buf_beta, fac, &
T2, tmp_virial, coord, i, &
pbd_buf_mp2_beta, pbc_buf_mp2_beta, pad_buf_mp2_beta, pac_buf_mp2_beta)
pbd_buf_beta, pbc_buf_beta, pad_buf_beta, pac_buf_beta, &
fac, T2, tmp_virial, coord, i, &
pbd_buf_resp_beta, pbc_buf_resp_beta, &
pad_buf_resp_beta, pac_buf_resp_beta)
ELSE
CALL update_virial(nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), &
pbd_buf_beta, pbc_buf_beta, pad_buf_beta, pac_buf_beta, fac, &
T2, tmp_virial, coord, i)
pbd_buf_beta, pbc_buf_beta, pad_buf_beta, pac_buf_beta, &
fac, T2, tmp_virial, coord, i)
ENDIF
END IF
END DO
@ -1551,7 +1565,7 @@ CONTAINS
shm_stor_count_max_val = shm_stor_count_max_val + stor_count_max_val
!$OMP BARRIER
! IF(with_mp2_density) THEN
! IF(with_resp_density) THEN
! ! restore original load_balance_parameter
! actual_x_data%load_balance_parameter = load_balance_parameter_energy
! DEALLOCATE(load_balance_parameter_energy)
@ -1662,13 +1676,13 @@ CONTAINS
DEALLOCATE (last_sgf_global)
!$OMP MASTER
DEALLOCATE (full_density_alpha)
IF (with_mp2_density) THEN
DEALLOCATE (full_density_mp2)
IF (with_resp_density) THEN
DEALLOCATE (full_density_resp)
END IF
IF (nspins == 2) THEN
DEALLOCATE (full_density_beta)
IF (with_mp2_density) THEN
DEALLOCATE (full_density_mp2_beta)
IF (with_resp_density) THEN
DEALLOCATE (full_density_resp_beta)
END IF
END IF
IF (do_dynamic_load_balancing) THEN
@ -1683,11 +1697,11 @@ CONTAINS
DEALLOCATE (pad_buf)
DEALLOCATE (pac_buf)
IF (with_mp2_density) THEN
DEALLOCATE (pbd_buf_mp2)
DEALLOCATE (pbc_buf_mp2)
DEALLOCATE (pad_buf_mp2)
DEALLOCATE (pac_buf_mp2)
IF (with_resp_density) THEN
DEALLOCATE (pbd_buf_resp)
DEALLOCATE (pbc_buf_resp)
DEALLOCATE (pad_buf_resp)
DEALLOCATE (pac_buf_resp)
END IF
DO i = 1, max_pgf**2
@ -1711,11 +1725,11 @@ CONTAINS
IF (nspins == 2) THEN
DEALLOCATE (pbd_buf_beta, pbc_buf_beta, pad_buf_beta, pac_buf_beta)
IF (with_mp2_density) THEN
DEALLOCATE (pbd_buf_mp2_beta)
DEALLOCATE (pbc_buf_mp2_beta)
DEALLOCATE (pad_buf_mp2_beta)
DEALLOCATE (pac_buf_mp2_beta)
IF (with_resp_density) THEN
DEALLOCATE (pbd_buf_resp_beta)
DEALLOCATE (pbc_buf_resp_beta)
DEALLOCATE (pad_buf_resp_beta)
DEALLOCATE (pac_buf_resp_beta)
END IF
END IF
@ -2165,10 +2179,11 @@ CONTAINS
!> \param force storage loacation for forces
!> \param forces_map index table
!> \param coord which of the 12 coords to be updated
!> \param pbd_mp2 ...
!> \param pbc_mp2 ...
!> \param pad_mp2 ...
!> \param pac_mp2 ...
!> \param pbd_resp ...
!> \param pbc_resp ...
!> \param pad_resp ...
!> \param pac_resp ...
!> \param resp_only ...
!> \par History
!> 03.2009 created [Manuel Guidon]
!> \author Manuel Guidon
@ -2176,7 +2191,7 @@ CONTAINS
SUBROUTINE update_forces(ma_max, mb_max, mc_max, md_max, &
pbd, pbc, pad, pac, fac, &
prim, force, forces_map, coord, &
pbd_mp2, pbc_mp2, pad_mp2, pac_mp2)
pbd_resp, pbc_resp, pad_resp, pac_resp, resp_only)
INTEGER, INTENT(IN) :: ma_max, mb_max, mc_max, md_max
REAL(dp), DIMENSION(*), INTENT(IN) :: pbd, pbc, pad, pac
@ -2185,20 +2200,23 @@ CONTAINS
INTENT(IN) :: prim
TYPE(qs_force_type), DIMENSION(:), POINTER :: force
INTEGER, INTENT(IN) :: forces_map(4, 2), coord
REAL(dp), DIMENSION(*), INTENT(IN), OPTIONAL :: pbd_mp2, pbc_mp2, pad_mp2, pac_mp2
REAL(dp), DIMENSION(*), INTENT(IN), OPTIONAL :: pbd_resp, pbc_resp, pad_resp, pac_resp
LOGICAL, INTENT(IN), OPTIONAL :: resp_only
INTEGER :: ma, mb, mc, md, p_index
LOGICAL :: with_mp2_density
REAL(dp) :: temp1, temp1_mp2, temp2, temp3, &
temp3_mp2, temp4
LOGICAL :: full_force, with_resp_density
REAL(dp) :: temp1, temp1_resp, temp3, temp3_resp, &
temp4, teresp
with_mp2_density = .FALSE.
IF (PRESENT(pbd_mp2) .AND. &
PRESENT(pbc_mp2) .AND. &
PRESENT(pad_mp2) .AND. &
PRESENT(pac_mp2)) with_mp2_density = .TRUE.
with_resp_density = .FALSE.
IF (PRESENT(pbd_resp) .AND. &
PRESENT(pbc_resp) .AND. &
PRESENT(pad_resp) .AND. &
PRESENT(pac_resp)) with_resp_density = .TRUE.
IF (with_mp2_density) THEN
IF (with_resp_density) THEN
full_force = .TRUE.
IF (PRESENT(resp_only)) full_force = .NOT. resp_only
p_index = 0
temp4 = 0.0_dp
DO md = 1, md_max
@ -2206,20 +2224,24 @@ CONTAINS
DO mb = 1, mb_max
temp1 = pbc((mc - 1)*mb_max + mb)*fac
temp3 = pbd((md - 1)*mb_max + mb)*fac
temp1_mp2 = pbc_mp2((mc - 1)*mb_max + mb)*fac
temp3_mp2 = pbd_mp2((md - 1)*mb_max + mb)*fac
temp1_resp = pbc_resp((mc - 1)*mb_max + mb)*fac
temp3_resp = pbd_resp((md - 1)*mb_max + mb)*fac
DO ma = 1, ma_max
p_index = p_index + 1
! HF-SCF
temp2 = temp1*pad((md - 1)*ma_max + ma) + &
temp3*pac((mc - 1)*ma_max + ma)
! MP2+HF
temp2 = temp2 + &
pac((mc - 1)*ma_max + ma)*temp3_mp2 + &
pac_mp2((mc - 1)*ma_max + ma)*temp3 + &
pad((md - 1)*ma_max + ma)*temp1_mp2 + &
pad_mp2((md - 1)*ma_max + ma)*temp1
temp4 = temp4 + temp2*prim(p_index)
IF (full_force) THEN
teresp = temp1*pad((md - 1)*ma_max + ma) + &
temp3*pac((mc - 1)*ma_max + ma)
ELSE
teresp = 0.0_dp
END IF
! RESP+HF
teresp = teresp + &
pac((mc - 1)*ma_max + ma)*temp3_resp + &
pac_resp((mc - 1)*ma_max + ma)*temp3 + &
pad((md - 1)*ma_max + ma)*temp1_resp + &
pad_resp((md - 1)*ma_max + ma)*temp1
temp4 = temp4 + teresp*prim(p_index)
END DO !ma
END DO !mb
END DO !mc
@ -2234,9 +2256,9 @@ CONTAINS
temp3 = pbd((md - 1)*mb_max + mb)*fac
DO ma = 1, ma_max
p_index = p_index + 1
temp2 = temp1*pad((md - 1)*ma_max + ma) + &
temp3*pac((mc - 1)*ma_max + ma)
temp4 = temp4 + temp2*prim(p_index)
teresp = temp1*pad((md - 1)*ma_max + ma) + &
temp3*pac((mc - 1)*ma_max + ma)
temp4 = temp4 + teresp*prim(p_index)
END DO !ma
END DO !mb
END DO !mc
@ -2267,10 +2289,10 @@ CONTAINS
!> \param tmp_virial ...
!> \param coord which of the 12 coords to be updated
!> \param l ...
!> \param pbd_mp2 ...
!> \param pbc_mp2 ...
!> \param pad_mp2 ...
!> \param pac_mp2 ...
!> \param pbd_resp ...
!> \param pbc_resp ...
!> \param pad_resp ...
!> \param pac_resp ...
!> \par History
!> 03.2009 created [Manuel Guidon]
!> \author Manuel Guidon
@ -2278,7 +2300,7 @@ CONTAINS
SUBROUTINE update_virial(ma_max, mb_max, mc_max, md_max, &
pbd, pbc, pad, pac, fac, &
prim, tmp_virial, coord, l, &
pbd_mp2, pbc_mp2, pad_mp2, pac_mp2)
pbd_resp, pbc_resp, pad_resp, pac_resp)
INTEGER, INTENT(IN) :: ma_max, mb_max, mc_max, md_max
REAL(dp), DIMENSION(*), INTENT(IN) :: pbd, pbc, pad, pac
@ -2287,20 +2309,20 @@ CONTAINS
INTENT(IN) :: prim
REAL(dp) :: tmp_virial(3, 3)
INTEGER, INTENT(IN) :: coord, l
REAL(dp), DIMENSION(*), INTENT(IN), OPTIONAL :: pbd_mp2, pbc_mp2, pad_mp2, pac_mp2
REAL(dp), DIMENSION(*), INTENT(IN), OPTIONAL :: pbd_resp, pbc_resp, pad_resp, pac_resp
INTEGER :: i, j, ma, mb, mc, md, p_index
LOGICAL :: with_mp2_density
REAL(dp) :: temp1, temp1_mp2, temp2, temp3, &
temp3_mp2, temp4
LOGICAL :: with_resp_density
REAL(dp) :: temp1, temp1_resp, temp3, temp3_resp, &
temp4, teresp
with_mp2_density = .FALSE.
IF (PRESENT(pbd_mp2) .AND. &
PRESENT(pbc_mp2) .AND. &
PRESENT(pad_mp2) .AND. &
PRESENT(pac_mp2)) with_mp2_density = .TRUE.
with_resp_density = .FALSE.
IF (PRESENT(pbd_resp) .AND. &
PRESENT(pbc_resp) .AND. &
PRESENT(pad_resp) .AND. &
PRESENT(pac_resp)) with_resp_density = .TRUE.
IF (with_mp2_density) THEN
IF (with_resp_density) THEN
p_index = 0
temp4 = 0.0_dp
DO md = 1, md_max
@ -2308,20 +2330,20 @@ CONTAINS
DO mb = 1, mb_max
temp1 = pbc((mc - 1)*mb_max + mb)*fac
temp3 = pbd((md - 1)*mb_max + mb)*fac
temp1_mp2 = pbc_mp2((mc - 1)*mb_max + mb)*fac
temp3_mp2 = pbd_mp2((md - 1)*mb_max + mb)*fac
temp1_resp = pbc_resp((mc - 1)*mb_max + mb)*fac
temp3_resp = pbd_resp((md - 1)*mb_max + mb)*fac
DO ma = 1, ma_max
p_index = p_index + 1
! HF-SCF
temp2 = temp1*pad((md - 1)*ma_max + ma) + &
temp3*pac((mc - 1)*ma_max + ma)
! MP2+HF
temp2 = temp2 + &
pac((mc - 1)*ma_max + ma)*temp3_mp2 + &
pac_mp2((mc - 1)*ma_max + ma)*temp3 + &
pad((md - 1)*ma_max + ma)*temp1_mp2 + &
pad_mp2((md - 1)*ma_max + ma)*temp1
temp4 = temp4 + temp2*prim(p_index)
teresp = temp1*pad((md - 1)*ma_max + ma) + &
temp3*pac((mc - 1)*ma_max + ma)
! RESP+HF
teresp = teresp + &
pac((mc - 1)*ma_max + ma)*temp3_resp + &
pac_resp((mc - 1)*ma_max + ma)*temp3 + &
pad((md - 1)*ma_max + ma)*temp1_resp + &
pad_resp((md - 1)*ma_max + ma)*temp1
temp4 = temp4 + teresp*prim(p_index)
END DO !ma
END DO !mb
END DO !mc
@ -2336,9 +2358,9 @@ CONTAINS
temp3 = pbd((md - 1)*mb_max + mb)*fac
DO ma = 1, ma_max
p_index = p_index + 1
temp2 = temp1*pad((md - 1)*ma_max + ma) + &
temp3*pac((mc - 1)*ma_max + ma)
temp4 = temp4 + temp2*prim(p_index)
teresp = temp1*pad((md - 1)*ma_max + ma) + &
temp3*pac((mc - 1)*ma_max + ma)
temp4 = temp4 + teresp*prim(p_index)
END DO !ma
END DO !mb
END DO !mc

13
src/input_constants.F
View file

@ -1020,13 +1020,16 @@ MODULE input_constants
kg_tnadd_atomic = 200, &
kg_tnadd_none = 300
! kg energy corrections
INTEGER, PARAMETER, PUBLIC :: kg_ec_diagonalization = 1001
INTEGER, PARAMETER, PUBLIC :: kg_ec_functional_harris = 2001
INTEGER, PARAMETER, PUBLIC :: kg_cholesky = 3001
! non-scf energy corrections
INTEGER, PARAMETER, PUBLIC :: ec_diagonalization = 1001, &
ec_curvy_steps = 1002, &
ec_matrix_sign = 1003, &
ec_matrix_trs4 = 1004, &
ec_matrix_tc2 = 1005
INTEGER, PARAMETER, PUBLIC :: ec_functional_harris = 2001
! swarm parameters
INTEGER, PARAMETER, PUBLIC :: swarm_do_glbopt = 1

103
src/input_cp2k_dft.F
View file

@ -62,11 +62,10 @@ MODULE input_cp2k_dft
eri_operator_coulomb, eri_operator_erf, eri_operator_erfc, eri_operator_gaussian, &
eri_operator_yukawa, gaussian, general_roks, gto_cartesian, gto_spherical, hf_model, &
high_spin_roks, history_guess, jacobian_fd1, jacobian_fd1_backward, jacobian_fd1_central, &
jacobian_fd2, jacobian_fd2_backward, kg_cholesky, kg_color_dsatur, kg_color_greedy, &
kg_ec_diagonalization, kg_ec_functional_harris, kg_tnadd_atomic, kg_tnadd_embed, &
kg_tnadd_embed_ri, kg_tnadd_none, ls_2pnt, ls_3pnt, ls_gold, ls_none, mao_basis_ext, &
mao_basis_orb, mao_basis_prim, mao_projection, mopac_guess, no_excitations, no_guess, &
numerical, oe_gllb, oe_lb, oe_none, oe_saop, oe_sic, op_loc_berry, op_loc_boys, &
jacobian_fd2, jacobian_fd2_backward, kg_color_dsatur, kg_color_greedy, kg_tnadd_atomic, &
kg_tnadd_embed, kg_tnadd_embed_ri, kg_tnadd_none, ls_2pnt, ls_3pnt, ls_gold, ls_none, &
mao_basis_ext, mao_basis_orb, mao_basis_prim, mao_projection, mopac_guess, no_excitations, &
no_guess, numerical, oe_gllb, oe_lb, oe_none, oe_saop, oe_sic, op_loc_berry, op_loc_boys, &
op_loc_pipek, orb_dx2, orb_dxy, orb_dy2, orb_dyz, orb_dz2, orb_dzx, orb_px, orb_py, &
orb_pz, orb_s, ot_algo_irac, ot_algo_taylor_or_diag, ot_chol_irac, ot_lwdn_irac, &
ot_mini_broyden, ot_mini_cg, ot_mini_diis, ot_mini_sd, ot_poly_irac, ot_precond_full_all, &
@ -99,6 +98,7 @@ MODULE input_cp2k_dft
xas_tp_xfh, xas_tp_xhh, xes_tp_val
USE input_cp2k_almo, ONLY: create_almo_scf_section
USE input_cp2k_distribution, ONLY: create_distribution_section
USE input_cp2k_ec, ONLY: create_ec_section
USE input_cp2k_field, ONLY: create_efield_section,&
create_per_efield_section
USE input_cp2k_kpoints, ONLY: create_kpoints_section
@ -343,6 +343,10 @@ CONTAINS
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)
CALL create_admm_section(subsection)
CALL section_add_subsection(section, subsection)
CALL section_release(subsection)
@ -5163,7 +5167,7 @@ CONTAINS
routineP = moduleN//':'//routineN
TYPE(keyword_type), POINTER :: keyword
TYPE(section_type), POINTER :: print_key, subsection, subsubsection
TYPE(section_type), POINTER :: print_key, subsection
CPASSERT(.NOT. ASSOCIATED(section))
CALL section_create(section, __LOCATION__, name="KG_METHOD", &
@ -5256,93 +5260,6 @@ CONTAINS
CALL section_add_subsection(section, subsection)
CALL section_release(subsection)
!!!!!!!!!!!!!!! ENERGY_CORRECTION SECTION !!!!!!!!!!!!!!!!!!!
NULLIFY (keyword, subsection, print_key)
CALL section_create(subsection, __LOCATION__, name="ENERGY_CORRECTION", &
description="Sets the various options for the Energy Correction", &
n_keywords=0, n_subsections=1, repeats=.FALSE.)
CALL keyword_create(keyword, __LOCATION__, name="_SECTION_PARAMETERS_", &
description="Controls the activation of the energy_correction", &
usage="&ENERGY_CORRECTION T", &
default_l_val=.FALSE., &
lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(subsection, keyword)
CALL keyword_release(keyword)
! add a special XC section
NULLIFY (subsubsection)
CALL create_xc_section(subsubsection)
CALL section_add_subsection(subsection, subsubsection)
CALL section_release(subsubsection)
CALL keyword_create(keyword, __LOCATION__, name="ENERGY_FUNCTIONAL", &
description="Functional used in energy correction", &
usage="ENERGY_FUNCTIONAL HARRIS", &
default_i_val=kg_ec_functional_harris, &
enum_c_vals=s2a("HARRIS"), &
enum_desc=s2a("Harris functional"), &
enum_i_vals=(/kg_ec_functional_harris/))
CALL section_add_keyword(subsection, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="HARRIS_BASIS", &
description="Specifies the type of basis to be used for the KG energy correction."// &
"Options are: (1) the default orbital basis (ORBITAL);"// &
"(2) the primitive functions of the default orbital basis (PRIMITIVE);"// &
"(3) the basis set labeled in Kind section (HARRIS)", &
usage="HARRIS_BASIS ORBITAL", &
type_of_var=char_t, default_c_val="ORBITAL", n_var=-1)
CALL section_add_keyword(subsection, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="MAO", &
description="Use modified atomic orbitals (MAO) to solve Harris equation", &
usage="MAO T", default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(subsection, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="MAO_MAX_ITER", &
description="Maximum iterations in MAO optimization. ", &
usage="MAO_MAX_ITER 100 ", default_i_val=0)
CALL section_add_keyword(subsection, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="MAO_EPS_GRAD", &
description="Threshold used for MAO iterations. ", &
usage="MAO_EPS_GRAD 1.0E-4 ", default_r_val=1.0E-5_dp)
CALL section_add_keyword(subsection, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="ALGORITHM", &
description="Algorithm used to solve KS equation", &
usage="ALGORITHM DIAGONALIZATION", &
default_i_val=kg_ec_diagonalization, &
enum_c_vals=s2a("DIAGONALIZATION"), &
enum_desc=s2a("Diagonalization of KS matrix."), &
enum_i_vals=(/kg_ec_diagonalization/))
CALL section_add_keyword(subsection, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="FACTORIZATION", &
description="Algorithm used to calculate factorization of overlap matrix", &
usage="FACTORIZATION CHOLESKY", &
default_i_val=kg_cholesky, &
enum_c_vals=s2a("CHOLESKY"), &
enum_desc=s2a("Cholesky factorization of overlap matrix"), &
enum_i_vals=(/kg_cholesky/))
CALL section_add_keyword(subsection, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="EPS_DEFAULT", &
description="Threshold used for accuracy estimates within energy correction algorithms. ", &
usage="EPS_DEFAULT 1.0E-6 ", default_r_val=1.0E-12_dp)
CALL section_add_keyword(subsection, keyword)
CALL keyword_release(keyword)
CALL section_add_subsection(section, subsection)
CALL section_release(subsection)
END SUBROUTINE create_kg_section
! **************************************************************************************************

333
src/input_cp2k_ec.F Normal file
View file

@ -0,0 +1,333 @@
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 - 2020 CP2K developers group !
!--------------------------------------------------------------------------------------------------!
! **************************************************************************************************
!> \brief function that build the dft section of the input
!> \par History
!> 10.2005 moved out of input_cp2k [fawzi]
!> \author fawzi
! **************************************************************************************************
MODULE input_cp2k_ec
USE bibliography, ONLY: Niklasson2003,&
VandeVondele2012
USE input_constants, ONLY: &
ec_curvy_steps, ec_diagonalization, ec_functional_harris, ec_matrix_sign, ec_matrix_tc2, &
ec_matrix_trs4, kg_cholesky, ls_cluster_atomic, ls_cluster_molecular, &
ls_s_inversion_hotelling, ls_s_inversion_sign_sqrt, ls_s_preconditioner_atomic, &
ls_s_preconditioner_molecular, ls_s_preconditioner_none, ls_s_sqrt_ns, ls_s_sqrt_proot, &
ls_scf_sign_ns, ls_scf_sign_proot
USE input_cp2k_xc, ONLY: create_xc_section
USE input_keyword_types, ONLY: keyword_create,&
keyword_release,&
keyword_type
USE input_section_types, ONLY: section_add_keyword,&
section_add_subsection,&
section_create,&
section_release,&
section_type
USE input_val_types, ONLY: char_t
USE kinds, ONLY: dp
USE string_utilities, ONLY: s2a
#include "./base/base_uses.f90"
IMPLICIT NONE
PRIVATE
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'input_cp2k_ec'
PUBLIC :: create_ec_section
CONTAINS
! **************************************************************************************************
!> \brief creates the ENERGY CORRECTION section
!> \param section ...
!> \author JGH
! **************************************************************************************************
SUBROUTINE create_ec_section(section)
TYPE(section_type), POINTER :: section
CHARACTER(len=*), PARAMETER :: routineN = 'create_ec_section', &
routineP = moduleN//':'//routineN
TYPE(keyword_type), POINTER :: keyword
TYPE(section_type), POINTER :: subsection
CPASSERT(.NOT. ASSOCIATED(section))
NULLIFY (keyword)
CALL section_create(section, __LOCATION__, name="ENERGY_CORRECTION", &
description="Sets the various options for the Energy Correction", &
n_keywords=0, n_subsections=1, repeats=.FALSE.)
CALL keyword_create(keyword, __LOCATION__, name="_SECTION_PARAMETERS_", &
description="Controls the activation of the energy_correction", &
usage="&ENERGY_CORRECTION T", &
default_l_val=.FALSE., &
lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
! add a special XC section
NULLIFY (subsection)
CALL create_xc_section(subsection)
CALL section_add_subsection(section, subsection)
CALL section_release(subsection)
! add a section for solver keywords
NULLIFY (subsection)
CALL create_ec_solver_section(subsection)
CALL section_add_subsection(section, subsection)
CALL section_release(subsection)
CALL keyword_create(keyword, __LOCATION__, name="ENERGY_FUNCTIONAL", &
description="Functional used in energy correction", &
usage="ENERGY_FUNCTIONAL HARRIS", &
default_i_val=ec_functional_harris, &
enum_c_vals=s2a("HARRIS"), &
enum_desc=s2a("Harris functional"), &
enum_i_vals=(/ec_functional_harris/))
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="HARRIS_BASIS", &
description="Specifies the type of basis to be used for the KG energy correction."// &
"Options are: (1) the default orbital basis (ORBITAL);"// &
"(2) the primitive functions of the default orbital basis (PRIMITIVE);"// &
"(3) the basis set labeled in Kind section (HARRIS)", &
usage="HARRIS_BASIS ORBITAL", &
type_of_var=char_t, default_c_val="ORBITAL", n_var=-1)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="MAO", &
description="Use modified atomic orbitals (MAO) to solve Harris equation", &
usage="MAO 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="MAO_MAX_ITER", &
description="Maximum iterations in MAO optimization. ", &
usage="MAO_MAX_ITER 100 ", default_i_val=0)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="MAO_EPS_GRAD", &
description="Threshold used for MAO iterations. ", &
usage="MAO_EPS_GRAD 1.0E-4 ", default_r_val=1.0E-5_dp)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="ALGORITHM", &
description="Algorithm used to solve KS equation", &
usage="ALGORITHM DIAGONALIZATION", &
default_i_val=ec_diagonalization, &
enum_c_vals=s2a("DIAGONALIZATION", "CURVY_STEPS", &
"MATRIX_SIGN", "TRS4", "TC2"), &
enum_desc=s2a("Diagonalization of KS matrix.", &
"Head-Grodon curvy step algorithm", &
"Matrix Sign algorithm", &
"Trace resetting trs4 algorithm", &
"Trace resetting tc2 algorithm"), &
enum_i_vals=(/ec_diagonalization, ec_curvy_steps, ec_matrix_sign, &
ec_matrix_trs4, ec_matrix_tc2/))
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="FACTORIZATION", &
description="Algorithm used to calculate factorization of overlap matrix", &
usage="FACTORIZATION CHOLESKY", &
default_i_val=kg_cholesky, &
enum_c_vals=s2a("CHOLESKY"), &
enum_desc=s2a("Cholesky factorization of overlap matrix"), &
enum_i_vals=(/kg_cholesky/))
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="EPS_DEFAULT", &
description="Threshold used for accuracy estimates within energy correction. ", &
usage="EPS_DEFAULT 1.0E-6 ", default_r_val=1.0E-12_dp)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
END SUBROUTINE create_ec_section
! **************************************************************************************************
!> \brief creates the linear scaling solver section
!> \param section ...
!> \author Joost VandeVondele [2010-10], JGH [2019-12]
! **************************************************************************************************
SUBROUTINE create_ec_solver_section(section)
TYPE(section_type), POINTER :: section
CHARACTER(len=*), PARAMETER :: routineN = 'create_ec_solver_section', &
routineP = moduleN//':'//routineN
TYPE(keyword_type), POINTER :: keyword
CPASSERT(.NOT. ASSOCIATED(section))
CALL section_create(section, __LOCATION__, name="LS_SOLVER", &
description="Specifies the parameters of the linear scaling solver routines", &
n_keywords=24, n_subsections=3, repeats=.FALSE., &
citations=(/VandeVondele2012/))
NULLIFY (keyword)
CALL keyword_create(keyword, __LOCATION__, name="EPS_FILTER", &
description="Threshold used for filtering matrix operations.", &
usage="EPS_FILTER 1.0E-7", default_r_val=1.0E-6_dp)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="EPS_LANCZOS", &
description="Threshold used for lanczos estimates.", &
usage="EPS_LANCZOS 1.0E-4", default_r_val=1.0E-3_dp)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="MAX_ITER_LANCZOS", &
description="Maximum number of lanczos iterations.", &
usage="MAX_ITER_LANCZOS ", default_i_val=128)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="MU", &
description="Value (or initial guess) for the chemical potential,"// &
" i.e. some suitable energy between HOMO and LUMO energy.", &
usage="MU 0.0", default_r_val=-0.1_dp)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="FIXED_MU", &
description="Should the calculation be performed at fixed chemical potential,"// &
" or should it be found fixing the number of electrons", &
usage="FIXED_MU .TRUE.", 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="EXTRAPOLATION_ORDER", &
description="Number of previous matrices used for the ASPC extrapolation of the initial guess. "// &
"0 implies that an atomic guess is used at each step. "// &
"low (1-2) will result in a drift of the constant of motion during MD. "// &
"high (>5) might be somewhat unstable, leading to more SCF iterations.", &
usage="EXTRAPOLATION_ORDER 3", default_i_val=4)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="S_PRECONDITIONER", &
description="Preconditions S with some appropriate form.", &
usage="S_PRECONDITIONER MOLECULAR", &
default_i_val=ls_s_preconditioner_atomic, &
enum_c_vals=s2a("NONE", "ATOMIC", "MOLECULAR"), &
enum_desc=s2a("No preconditioner", &
"Using atomic blocks", &
"Using molecular sub-blocks. Recommended if molecules are defined and not too large."), &
enum_i_vals=(/ls_s_preconditioner_none, ls_s_preconditioner_atomic, ls_s_preconditioner_molecular/))
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="S_SQRT_METHOD", &
description="Method for the caclulation of the sqrt of S.", &
usage="S_SQRT_METHOD NEWTONSCHULZ", &
default_i_val=ls_s_sqrt_ns, &
enum_c_vals=s2a("NEWTONSCHULZ", "PROOT"), &
enum_desc=s2a("Using a Newton-Schulz-like iteration", &
"Using the p-th root method."), &
enum_i_vals=(/ls_s_sqrt_ns, ls_s_sqrt_proot/))
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="S_SQRT_ORDER", &
variants=s2a("SIGN_SQRT_ORDER"), &
description="Order of the iteration method for the calculation of the sqrt of S.", &
usage="S_SQRT_ORDER 3", default_i_val=3)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="SIGN_METHOD", &
description="Method used for the computation of the sign matrix.", &
usage="SIGN_METHOD NEWTONSCHULZ", &
default_i_val=ls_scf_sign_ns, &
citations=(/VandeVondele2012, Niklasson2003/), &
enum_c_vals=s2a("NEWTONSCHULZ", "PROOT"), &
enum_desc=s2a("Newton-Schulz iteration.", &
"p-th order root iteration"), &
enum_i_vals=(/ls_scf_sign_ns, ls_scf_sign_proot/))
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="SIGN_ORDER", &
description="Order of the method used for the computation of the sign matrix.", &
usage="SIGN_ORDER 2", &
default_i_val=2)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="DYNAMIC_THRESHOLD", &
description="Should the threshold for the purification be chosen dynamically", &
usage="DYNAMIC_THRESHOLD .TRUE.", 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="NON_MONOTONIC", &
description="Should the purification be performed non-monotonically. Relevant for TC2 only.", &
usage="NON_MONOTONIC .TRUE.", default_l_val=.TRUE., lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create( &
keyword, __LOCATION__, name="MATRIX_CLUSTER_TYPE", &
description="Specify how atomic blocks should be clustered in the used matrices, in order to improve flop rate, "// &
"and possibly speedup the matrix multiply. Note that the atomic s_preconditioner can not be used."// &
"Furthermore, since screening is on matrix blocks, "// &
"slightly more accurate results can be expected with molecular.", &
usage="MATRIX_CLUSTER_TYPE MOLECULAR", &
default_i_val=ls_cluster_atomic, &
enum_c_vals=s2a("ATOMIC", "MOLECULAR"), &
enum_desc=s2a("Using atomic blocks", &
"Using molecular blocks."), &
enum_i_vals=(/ls_cluster_atomic, ls_cluster_molecular/))
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="SINGLE_PRECISION_MATRICES", &
description="Matrices used within the LS code can be either double or single precision.", &
usage="SINGLE_PRECISION_MATRICES", 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="S_INVERSION", &
description="Method used to compute the inverse of S.", &
usage="S_PRECONDITIONER MOLECULAR", &
default_i_val=ls_s_inversion_sign_sqrt, &
enum_c_vals=s2a("SIGN_SQRT", "HOTELLING"), &
enum_desc=s2a("Using the inverse sqrt as obtained from sign function iterations.", &
"Using the Hotellign iteration."), &
enum_i_vals=(/ls_s_inversion_sign_sqrt, ls_s_inversion_hotelling/))
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="REPORT_ALL_SPARSITIES", &
description="Run the sparsity report at the end of the SCF", &
usage="REPORT_ALL_SPARSITIES", default_l_val=.TRUE., lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="CHECK_S_INV", &
description="Perform an accuracy check on the inverse/sqrt of the s matrix.", &
usage="CHECK_S_INV", default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
! CALL create_ls_curvy_section(subsection)
! CALL section_add_subsection(section, subsection)
! CALL section_release(subsection)
! CALL create_chebyshev_section(subsection)
! CALL section_add_subsection(section, subsection)
! CALL section_release(subsection)
END SUBROUTINE create_ec_solver_section
END MODULE input_cp2k_ec

958
src/kg_energy_corrections.F
View file

@ -1,958 +0,0 @@
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 - 2020 CP2K developers group !
!--------------------------------------------------------------------------------------------------!
! **************************************************************************************************
!> \brief Routines for a Harris type energy correction on top of a
!> Kim-Gordon calculation
!> \par History
!> 03.2014 created
!> \author JGH
! **************************************************************************************************
MODULE kg_energy_corrections
USE atomic_kind_types, ONLY: atomic_kind_type,&
get_atomic_kind
USE basis_set_types, ONLY: get_gto_basis_set,&
gto_basis_set_type
USE cell_types, ONLY: cell_type
USE core_ppl, ONLY: build_core_ppl
USE core_ppnl, ONLY: build_core_ppnl
USE cp_blacs_env, ONLY: cp_blacs_env_type
USE cp_control_types, ONLY: dft_control_type
USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
copy_fm_to_dbcsr,&
dbcsr_allocate_matrix_set,&
dbcsr_deallocate_matrix_set
USE cp_fm_basic_linalg, ONLY: cp_fm_triangular_invert
USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose
USE cp_fm_diag, ONLY: choose_eigv_solver
USE cp_fm_struct, ONLY: cp_fm_struct_create,&
cp_fm_struct_release,&
cp_fm_struct_type
USE cp_fm_types, ONLY: cp_fm_create,&
cp_fm_release,&
cp_fm_set_all,&
cp_fm_to_fm_triangular,&
cp_fm_type
USE cp_log_handling, ONLY: cp_get_default_logger,&
cp_logger_get_default_unit_nr,&
cp_logger_type
USE cp_para_types, ONLY: cp_para_env_type
USE dbcsr_api, ONLY: &
dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_desymmetrize, dbcsr_distribution_type, &
dbcsr_dot, dbcsr_filter, dbcsr_get_info, dbcsr_init_p, dbcsr_multiply, dbcsr_p_type, &
dbcsr_release, dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry, dbcsr_type_symmetric
USE distribution_1d_types, ONLY: distribution_1d_type
USE distribution_2d_types, ONLY: distribution_2d_type
USE external_potential_types, ONLY: get_potential,&
gth_potential_type,&
sgp_potential_type
USE input_constants, ONLY: kg_ec_diagonalization,&
kg_ec_functional_harris
USE input_section_types, ONLY: section_vals_val_get
USE kg_environment_types, ONLY: energy_correction_type,&
kg_environment_type
USE kinds, ONLY: default_string_length,&
dp
USE mao_basis, ONLY: mao_generate_basis
USE molecule_types, ONLY: molecule_type
USE particle_types, ONLY: particle_type
USE pw_env_types, ONLY: pw_env_get,&
pw_env_type
USE pw_methods, ONLY: pw_axpy,&
pw_integral_ab,&
pw_scale,&
pw_transfer
USE pw_poisson_methods, ONLY: pw_poisson_solve
USE pw_poisson_types, ONLY: pw_poisson_type
USE pw_pool_types, ONLY: pw_pool_create_pw,&
pw_pool_give_back_pw,&
pw_pool_p_type,&
pw_pool_type
USE pw_types, ONLY: COMPLEXDATA1D,&
REALDATA3D,&
REALSPACE,&
RECIPROCALSPACE,&
pw_p_type
USE qs_core_energies, ONLY: calculate_ecore_overlap,&
calculate_ecore_self
USE qs_dispersion_pairpot, ONLY: calculate_dispersion_pairpot
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
USE qs_force_types, ONLY: qs_force_type
USE qs_integrate_potential, ONLY: integrate_v_rspace
USE qs_kind_types, ONLY: get_qs_kind,&
get_qs_kind_set,&
qs_kind_type
USE qs_kinetic, ONLY: build_kinetic_matrix
USE qs_ks_methods, ONLY: calc_rho_tot_gspace
USE qs_ks_types, ONLY: qs_ks_env_type
USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
USE qs_neighbor_lists, ONLY: atom2d_build,&
atom2d_cleanup,&
build_neighbor_lists,&
local_atoms_type,&
pair_radius_setup
USE qs_overlap, ONLY: build_overlap_matrix
USE qs_rho_types, ONLY: qs_rho_get,&
qs_rho_type
USE qs_vxc, ONLY: qs_vxc_create
USE task_list_methods, ONLY: generate_qs_task_list
USE task_list_types, ONLY: allocate_task_list,&
deallocate_task_list
USE virial_types, ONLY: virial_type
#include "./base/base_uses.f90"
IMPLICIT NONE
PRIVATE
! *** Global parameters ***
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'kg_energy_corrections'
PUBLIC :: kg_energy_correction
CONTAINS
! **************************************************************************************************
!> \brief Energy correction to a KG simulation
!>
!> \param qs_env ...
!> \param calculate_forces ...
!> \par History
!> 03.2014 created
!> \author JGH
! **************************************************************************************************
SUBROUTINE kg_energy_correction(qs_env, calculate_forces)
TYPE(qs_environment_type), POINTER :: qs_env
LOGICAL, INTENT(IN), OPTIONAL :: calculate_forces
CHARACTER(len=*), PARAMETER :: routineN = 'kg_energy_correction', &
routineP = moduleN//':'//routineN
INTEGER :: handle, unit_nr
LOGICAL :: my_calc_forces
TYPE(cp_logger_type), POINTER :: logger
TYPE(energy_correction_type), POINTER :: ec_env
TYPE(kg_environment_type), POINTER :: kg_env
CALL timeset(routineN, handle)
my_calc_forces = .FALSE.
IF (PRESENT(calculate_forces)) my_calc_forces = calculate_forces
NULLIFY (ec_env, kg_env)
CALL get_qs_env(qs_env=qs_env, kg_env=kg_env)
! Check for energy correction
IF (kg_env%energy_correction) THEN
ec_env => kg_env%ec_env
ec_env%etotal = 0.0_dp
ec_env%eband = 0.0_dp
ec_env%ehartree = 0.0_dp
ec_env%exc = 0.0_dp
ec_env%vhxc = 0.0_dp
ec_env%edispersion = 0.0_dp
logger => cp_get_default_logger()
IF (logger%para_env%ionode) THEN
unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
ELSE
unit_nr = -1
ENDIF
IF (unit_nr > 0) THEN
WRITE (unit_nr, '(/,T2,A)') '!-----------------------------------------------------------------------------!'
WRITE (unit_nr, '(T2,A,A,A,A,A)') "!", REPEAT("-", 27), " KG energy correction ", REPEAT("-", 28), "!"
END IF
! build neighbor and task lists
CALL ec_build_neighborlist(qs_env, ec_env)
CALL ec_build_core_hamiltonian(qs_env, ec_env, my_calc_forces)
CALL ec_build_ks_matrix(qs_env, ec_env)
! MAO basis
IF (ec_env%mao) THEN
IF (ASSOCIATED(ec_env%mao_coef)) CALL dbcsr_deallocate_matrix_set(ec_env%mao_coef)
NULLIFY (ec_env%mao_coef)
CALL mao_generate_basis(qs_env, ec_env%mao_coef, ref_basis_set="HARRIS", molecular=.TRUE., &
max_iter=ec_env%mao_max_iter, eps_grad=ec_env%mao_eps_grad, unit_nr=unit_nr)
END IF
CALL ec_ks_solver(qs_env, ec_env)
CALL ec_energy(qs_env, ec_env, unit_nr)
IF (calculate_forces) THEN
! CALL kg_response_solver(qs_env,ec_env)
END IF
IF (unit_nr > 0) THEN
WRITE (unit_nr, '(/,T2,A)') '!-----------------------------------------------------------------------------!'
END IF
END IF
CALL timestop(handle)
END SUBROUTINE kg_energy_correction
! **************************************************************************************************
!> \brief Construction of the Core Hamiltonian Matrix
!> Short version of qs_core_hamiltonian
!> \param qs_env ...
!> \param ec_env ...
!> \param calculate_forces ...
!> \author Creation (03.2014,JGH)
! **************************************************************************************************
SUBROUTINE ec_build_core_hamiltonian(qs_env, ec_env, calculate_forces)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(energy_correction_type), POINTER :: ec_env
LOGICAL, INTENT(IN) :: calculate_forces
CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_build_core_hamiltonian', &
routineP = moduleN//':'//routineN
INTEGER :: handle, nder, nimages
INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
LOGICAL :: use_virial
REAL(KIND=dp) :: eps_filter, eps_ppnl
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p
TYPE(dft_control_type), POINTER :: dft_control
TYPE(neighbor_list_set_p_type), DIMENSION(:), &
POINTER :: sab_orb, sac_ppl, sap_ppnl
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(qs_force_type), DIMENSION(:), POINTER :: force
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(qs_ks_env_type), POINTER :: ks_env
TYPE(virial_type), POINTER :: virial
IF (calculate_forces) THEN
CALL timeset(routineN//"_forces", handle)
ELSE
CALL timeset(routineN, handle)
ENDIF
! no forces yet
CPASSERT(.NOT. calculate_forces)
! no k-points possible
CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
nimages = dft_control%nimages
CPASSERT(nimages == 1)
! check for virial (currently no stress tensor available)
CALL get_qs_env(qs_env=qs_env, virial=virial)
use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
CPASSERT(.NOT. use_virial)
! get neighbor lists, we need the full sab_orb list from the ec_env
NULLIFY (sab_orb, sac_ppl, sap_ppnl)
sab_orb => ec_env%sab_orb
sac_ppl => ec_env%sac_ppl
sap_ppnl => ec_env%sap_ppnl
! Overlap and kinetic energy matrices
CALL get_qs_env(qs_env=qs_env, ks_env=ks_env)
eps_filter = dft_control%qs_control%eps_filter_matrix
nder = 0
CALL build_overlap_matrix(ks_env, nderivative=nder, matrixkp_s=ec_env%matrix_s, &
matrix_name="OVERLAP MATRIX", &
basis_type_a="HARRIS", &
basis_type_b="HARRIS", &
sab_nl=sab_orb)
CALL build_kinetic_matrix(ks_env, matrixkp_t=ec_env%matrix_t, &
matrix_name="KINETIC ENERGY MATRIX", &
basis_type="HARRIS", &
sab_nl=sab_orb, &
eps_filter=eps_filter)
! initialize H matrix
CALL dbcsr_allocate_matrix_set(ec_env%matrix_h, 1, 1)
ALLOCATE (ec_env%matrix_h(1, 1)%matrix)
CALL dbcsr_create(ec_env%matrix_h(1, 1)%matrix, template=ec_env%matrix_s(1, 1)%matrix)
CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_h(1, 1)%matrix, sab_orb)
! add kinetic energy
CALL dbcsr_copy(ec_env%matrix_h(1, 1)%matrix, ec_env%matrix_t(1, 1)%matrix, &
keep_sparsity=.TRUE., name="CORE HAMILTONIAN MATRIX")
! compute the ppl contribution to the core hamiltonian
CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, particle_set=particle_set, &
atomic_kind_set=atomic_kind_set)
NULLIFY (cell_to_index, matrix_p, force, virial)
use_virial = .FALSE.
IF (ASSOCIATED(sac_ppl)) THEN
CALL build_core_ppl(ec_env%matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, &
qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ppl, &
nimages, cell_to_index, "HARRIS")
END IF
! compute the ppnl contribution to the core hamiltonian ***
eps_ppnl = dft_control%qs_control%eps_ppnl
IF (ASSOCIATED(sap_ppnl)) THEN
CALL build_core_ppnl(ec_env%matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, &
qs_kind_set, atomic_kind_set, particle_set, sab_orb, sap_ppnl, eps_ppnl, &
nimages, cell_to_index, "HARRIS")
END IF
CALL timestop(handle)
END SUBROUTINE ec_build_core_hamiltonian
! **************************************************************************************************
!> \brief calculate the complete KS matrix
!> \param qs_env ...
!> \param ec_env ...
!> \par History
!> 03.2014 adapted from qs_ks_build_kohn_sham_matrix [JGH]
!> \author JGH
! **************************************************************************************************
SUBROUTINE ec_build_ks_matrix(qs_env, ec_env)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(energy_correction_type), POINTER :: ec_env
CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_build_ks_matrix', &
routineP = moduleN//':'//routineN
CHARACTER(LEN=default_string_length) :: headline
INTEGER :: handle, ispin, nspins
LOGICAL :: use_virial
REAL(dp) :: eexc, ehartree, eovrl, eself, evhxc
TYPE(cp_para_env_type), POINTER :: para_env
TYPE(dft_control_type), POINTER :: dft_control
TYPE(pw_env_type), POINTER :: pw_env
TYPE(pw_p_type) :: rho_tot_gspace, v_hartree_gspace, &
v_hartree_rspace
TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_r, tau_r, v_rspace, v_tau_rspace
TYPE(pw_poisson_type), POINTER :: poisson_env
TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
TYPE(qs_ks_env_type), POINTER :: ks_env
TYPE(qs_rho_type), POINTER :: rho
TYPE(virial_type), POINTER :: virial
CALL timeset(routineN, handle)
! get all information on the electronic density
NULLIFY (rho, ks_env)
CALL get_qs_env(qs_env=qs_env, rho=rho, virial=virial, dft_control=dft_control, &
para_env=para_env, ks_env=ks_env)
nspins = dft_control%nspins
use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
CPASSERT(.NOT. use_virial)
! Kohn-Sham matrix
IF (ASSOCIATED(ec_env%matrix_ks)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_ks)
CALL dbcsr_allocate_matrix_set(ec_env%matrix_ks, nspins, 1)
DO ispin = 1, nspins
headline = "KOHN-SHAM MATRIX"
ALLOCATE (ec_env%matrix_ks(ispin, 1)%matrix)
CALL dbcsr_create(ec_env%matrix_ks(ispin, 1)%matrix, name=TRIM(headline), &
template=ec_env%matrix_s(1, 1)%matrix, matrix_type=dbcsr_type_symmetric)
CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_ks(ispin, 1)%matrix, ec_env%sab_orb)
CALL dbcsr_set(ec_env%matrix_ks(ispin, 1)%matrix, 0.0_dp)
ENDDO
NULLIFY (pw_env)
CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
CPASSERT(ASSOCIATED(pw_env))
NULLIFY (auxbas_pw_pool, poisson_env, pw_pools)
! gets the tmp grids
CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
pw_pools=pw_pools, poisson_env=poisson_env)
! Calculate the Hartree potential
CALL pw_pool_create_pw(auxbas_pw_pool, &
v_hartree_gspace%pw, &
use_data=COMPLEXDATA1D, &
in_space=RECIPROCALSPACE)
CALL pw_pool_create_pw(auxbas_pw_pool, &
rho_tot_gspace%pw, &
use_data=COMPLEXDATA1D, &
in_space=RECIPROCALSPACE)
CALL pw_pool_create_pw(auxbas_pw_pool, &
v_hartree_rspace%pw, &
use_data=REALDATA3D, &
in_space=REALSPACE)
! Get the total density in g-space [ions + electrons]
CALL calc_rho_tot_gspace(rho_tot_gspace, qs_env, rho)
CALL pw_poisson_solve(poisson_env, rho_tot_gspace%pw, ehartree, &
v_hartree_gspace%pw)
CALL pw_transfer(v_hartree_gspace%pw, v_hartree_rspace%pw)
CALL pw_scale(v_hartree_rspace%pw, v_hartree_rspace%pw%pw_grid%dvol)
CALL pw_pool_give_back_pw(auxbas_pw_pool, v_hartree_gspace%pw)
CALL pw_pool_give_back_pw(auxbas_pw_pool, rho_tot_gspace%pw)
! v_rspace and v_tau_rspace are generated from the auxbas pool
NULLIFY (v_rspace, v_tau_rspace)
CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho, xc_section=ec_env%xc_section, &
vxc_rho=v_rspace, vxc_tau=v_tau_rspace, exc=eexc, just_energy=.FALSE.)
evhxc = 0.0_dp
CALL qs_rho_get(rho, rho_r=rho_r)
IF (ASSOCIATED(v_tau_rspace)) THEN
CALL qs_rho_get(rho, tau_r=tau_r)
END IF
DO ispin = 1, nspins
! Add v_hartree + v_xc = v_rspace
CALL pw_scale(v_rspace(ispin)%pw, v_rspace(ispin)%pw%pw_grid%dvol)
CALL pw_axpy(v_hartree_rspace%pw, v_rspace(ispin)%pw)
! integrate over potential <a|V|b>
CALL integrate_v_rspace(v_rspace=v_rspace(ispin), &
hmat=ec_env%matrix_ks(ispin, 1), &
qs_env=qs_env, &
calculate_forces=.FALSE., &
basis_type="HARRIS", &
task_list_external=ec_env%task_list)
IF (ASSOCIATED(v_tau_rspace)) THEN
! integrate over Tau-potential <nabla.a|V|nabla.b>
CALL pw_scale(v_tau_rspace(ispin)%pw, v_tau_rspace(ispin)%pw%pw_grid%dvol)
CALL integrate_v_rspace(v_rspace=v_tau_rspace(ispin), hmat=ec_env%matrix_ks(ispin, 1), &
qs_env=qs_env, calculate_forces=.FALSE., compute_tau=.TRUE., &
basis_type="HARRIS", &
task_list_external=ec_env%task_list)
END IF
! calclulate Int(vhxc*rho)dr and Int(vtau*tau)dr
evhxc = evhxc + pw_integral_ab(rho_r(ispin)%pw, v_rspace(ispin)%pw)/v_rspace(1)%pw%pw_grid%dvol
IF (ASSOCIATED(v_tau_rspace)) THEN
evhxc = evhxc + pw_integral_ab(tau_r(ispin)%pw, v_tau_rspace(ispin)%pw)/v_tau_rspace(ispin)%pw%pw_grid%dvol
END IF
END DO
! return pw grids
CALL pw_pool_give_back_pw(auxbas_pw_pool, v_hartree_rspace%pw)
DO ispin = 1, nspins
CALL pw_pool_give_back_pw(auxbas_pw_pool, v_rspace(ispin)%pw)
IF (ASSOCIATED(v_tau_rspace)) THEN
CALL pw_pool_give_back_pw(auxbas_pw_pool, v_tau_rspace(ispin)%pw)
END IF
ENDDO
! energies
CALL calculate_ecore_self(qs_env, E_self_core=eself)
CALL calculate_ecore_overlap(qs_env, para_env, calculate_forces=.FALSE., E_overlap_core=eovrl)
ec_env%exc = eexc
ec_env%ehartree = ehartree + eovrl + eself
ec_env%vhxc = evhxc
! add the core matrix
DO ispin = 1, nspins
CALL dbcsr_add(ec_env%matrix_ks(ispin, 1)%matrix, ec_env%matrix_h(1, 1)%matrix, &
alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
END DO
! At this point the ks matrix is up to date, filter it if requested
DO ispin = 1, nspins
CALL dbcsr_filter(ec_env%matrix_ks(ispin, 1)%matrix, &
dft_control%qs_control%eps_filter_matrix)
ENDDO
DEALLOCATE (v_rspace)
CALL timestop(handle)
END SUBROUTINE ec_build_ks_matrix
! **************************************************************************************************
!> \brief Solve KS equation for a given matrix
!> \param qs_env ...
!> \param ec_env ...
!> \par History
!> 03.2014 created [JGH]
!> \author JGH
! **************************************************************************************************
SUBROUTINE ec_ks_solver(qs_env, ec_env)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(energy_correction_type), POINTER :: ec_env
CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_ks_solver', routineP = moduleN//':'//routineN
CHARACTER(LEN=default_string_length) :: headline
INTEGER :: handle, ispin, nspins
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ksmat, pmat, smat
TYPE(dft_control_type), POINTER :: dft_control
CALL timeset(routineN, handle)
CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
nspins = dft_control%nspins
! create density matrix
IF (.NOT. ASSOCIATED(ec_env%matrix_p)) THEN
headline = "DENSITY MATRIX"
CALL dbcsr_allocate_matrix_set(ec_env%matrix_p, nspins, 1)
DO ispin = 1, nspins
ALLOCATE (ec_env%matrix_p(ispin, 1)%matrix)
CALL dbcsr_create(ec_env%matrix_p(ispin, 1)%matrix, name=TRIM(headline), &
template=ec_env%matrix_s(1, 1)%matrix)
CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_p(ispin, 1)%matrix, ec_env%sab_orb)
END DO
END IF
IF (ec_env%mao) THEN
CALL mao_create_matrices(ec_env, ksmat, smat, pmat)
ELSE
ksmat => ec_env%matrix_ks
smat => ec_env%matrix_s
pmat => ec_env%matrix_p
ENDIF
SELECT CASE (ec_env%ks_solver)
CASE (kg_ec_diagonalization)
CALL ec_diag_solver(qs_env, ksmat, smat, pmat)
CASE DEFAULT
CPASSERT(.FALSE.)
END SELECT
IF (ec_env%mao) THEN
CALL mao_release_matrices(ec_env, ksmat, smat, pmat)
ENDIF
CALL timestop(handle)
END SUBROUTINE ec_ks_solver
! **************************************************************************************************
!> \brief Create matrices with MAO sizes
!> \param ec_env ...
!> \param ksmat ...
!> \param smat ...
!> \param pmat ...
!> \par History
!> 08.2016 created [JGH]
!> \author JGH
! **************************************************************************************************
SUBROUTINE mao_create_matrices(ec_env, ksmat, smat, pmat)
TYPE(energy_correction_type), POINTER :: ec_env
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ksmat, smat, pmat
CHARACTER(LEN=*), PARAMETER :: routineN = 'mao_create_matrices', &
routineP = moduleN//':'//routineN
INTEGER :: handle, ispin, nspins
INTEGER, DIMENSION(:), POINTER :: col_blk_sizes
TYPE(dbcsr_distribution_type) :: dbcsr_dist
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mao_coef
TYPE(dbcsr_type) :: cgmat
CALL timeset(routineN, handle)
mao_coef => ec_env%mao_coef
NULLIFY (ksmat, smat, pmat)
nspins = SIZE(ec_env%matrix_ks, 1)
CALL dbcsr_get_info(mao_coef(1)%matrix, col_blk_size=col_blk_sizes, distribution=dbcsr_dist)
CALL dbcsr_allocate_matrix_set(ksmat, nspins, 1)
CALL dbcsr_allocate_matrix_set(smat, nspins, 1)
DO ispin = 1, nspins
ALLOCATE (ksmat(ispin, 1)%matrix)
CALL dbcsr_create(ksmat(ispin, 1)%matrix, dist=dbcsr_dist, name="MAO KS mat", &
matrix_type=dbcsr_type_symmetric, row_blk_size=col_blk_sizes, &
col_blk_size=col_blk_sizes, nze=0)
ALLOCATE (smat(ispin, 1)%matrix)
CALL dbcsr_create(smat(ispin, 1)%matrix, dist=dbcsr_dist, name="MAO S mat", &
matrix_type=dbcsr_type_symmetric, row_blk_size=col_blk_sizes, &
col_blk_size=col_blk_sizes, nze=0)
END DO
!
CALL dbcsr_create(cgmat, name="TEMP matrix", template=mao_coef(1)%matrix)
DO ispin = 1, nspins
CALL dbcsr_multiply("N", "N", 1.0_dp, ec_env%matrix_s(1, 1)%matrix, mao_coef(ispin)%matrix, &
0.0_dp, cgmat)
CALL dbcsr_multiply("T", "N", 1.0_dp, mao_coef(ispin)%matrix, cgmat, 0.0_dp, smat(ispin, 1)%matrix)
CALL dbcsr_multiply("N", "N", 1.0_dp, ec_env%matrix_ks(1, 1)%matrix, mao_coef(ispin)%matrix, &
0.0_dp, cgmat)
CALL dbcsr_multiply("T", "N", 1.0_dp, mao_coef(ispin)%matrix, cgmat, 0.0_dp, ksmat(ispin, 1)%matrix)
END DO
CALL dbcsr_release(cgmat)
CALL dbcsr_allocate_matrix_set(pmat, nspins, 1)
DO ispin = 1, nspins
ALLOCATE (pmat(ispin, 1)%matrix)
CALL dbcsr_create(pmat(ispin, 1)%matrix, template=smat(1, 1)%matrix)
CALL cp_dbcsr_alloc_block_from_nbl(pmat(ispin, 1)%matrix, ec_env%sab_orb)
END DO
CALL timestop(handle)
END SUBROUTINE mao_create_matrices
! **************************************************************************************************
!> \brief Release matrices with MAO sizes
!> \param ec_env ...
!> \param ksmat ...
!> \param smat ...
!> \param pmat ...
!> \par History
!> 08.2016 created [JGH]
!> \author JGH
! **************************************************************************************************
SUBROUTINE mao_release_matrices(ec_env, ksmat, smat, pmat)
TYPE(energy_correction_type), POINTER :: ec_env
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ksmat, smat, pmat
CHARACTER(LEN=*), PARAMETER :: routineN = 'mao_release_matrices', &
routineP = moduleN//':'//routineN
INTEGER :: handle, ispin, nspins
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mao_coef
TYPE(dbcsr_type) :: cgmat
CALL timeset(routineN, handle)
mao_coef => ec_env%mao_coef
nspins = SIZE(mao_coef, 1)
! save pmat in full basis format
CALL dbcsr_create(cgmat, name="TEMP matrix", template=mao_coef(1)%matrix)
DO ispin = 1, nspins
CALL dbcsr_multiply("N", "N", 1.0_dp, mao_coef(ispin)%matrix, pmat(ispin, 1)%matrix, 0.0_dp, cgmat)
CALL dbcsr_multiply("N", "T", 1.0_dp, mao_coef(ispin)%matrix, cgmat, 0.0_dp, &
ec_env%matrix_p(ispin, 1)%matrix, retain_sparsity=.TRUE.)
END DO
CALL dbcsr_release(cgmat)
CALL dbcsr_deallocate_matrix_set(ksmat)
CALL dbcsr_deallocate_matrix_set(smat)
CALL dbcsr_deallocate_matrix_set(pmat)
CALL timestop(handle)
END SUBROUTINE mao_release_matrices
! **************************************************************************************************
!> \brief Solve KS equation using diagonalization
!> \param qs_env ...
!> \param matrix_ks ...
!> \param matrix_s ...
!> \param matrix_p ...
!> \par History
!> 03.2014 created [JGH]
!> \author JGH
! **************************************************************************************************
SUBROUTINE ec_diag_solver(qs_env, matrix_ks, matrix_s, matrix_p)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_s, matrix_p
CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_diag_solver', routineP = moduleN//':'//routineN
INTEGER :: handle, info, ispin, nmo(2), nsize, &
nspins
REAL(KIND=dp) :: eps_filter, focc(2)
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues
TYPE(cp_blacs_env_type), POINTER :: blacs_env
TYPE(cp_fm_struct_type), POINTER :: fm_struct
TYPE(cp_fm_type), POINTER :: fm_ks, fm_mo, fm_ortho
TYPE(cp_para_env_type), POINTER :: para_env
TYPE(dbcsr_type), POINTER :: buf1_dbcsr, buf2_dbcsr, ortho_dbcsr, &
ref_matrix
TYPE(dft_control_type), POINTER :: dft_control
CALL timeset(routineN, handle)
NULLIFY (blacs_env, para_env)
CALL get_qs_env(qs_env=qs_env, blacs_env=blacs_env, para_env=para_env)
CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
eps_filter = dft_control%qs_control%eps_filter_matrix
nspins = dft_control%nspins
nmo = 0
CALL get_qs_env(qs_env=qs_env, nelectron_spin=nmo)
focc = 1._dp
IF (nspins == 1) THEN
focc = 2._dp
nmo(1) = nmo(1)/2
END IF
CALL dbcsr_get_info(matrix_ks(1, 1)%matrix, nfullrows_total=nsize)
ALLOCATE (eigenvalues(nsize))
NULLIFY (fm_ortho, fm_ks, fm_mo, fm_struct, ref_matrix)
CALL cp_fm_struct_create(fm_struct, context=blacs_env, nrow_global=nsize, &
ncol_global=nsize, para_env=para_env)
CALL cp_fm_create(fm_ortho, fm_struct)
CALL cp_fm_create(fm_ks, fm_struct)
CALL cp_fm_create(fm_mo, fm_struct)
CALL cp_fm_struct_release(fm_struct)
! factorization
ref_matrix => matrix_s(1, 1)%matrix
NULLIFY (ortho_dbcsr, buf1_dbcsr, buf2_dbcsr)
CALL dbcsr_init_p(ortho_dbcsr)
CALL dbcsr_create(ortho_dbcsr, template=ref_matrix, &
matrix_type=dbcsr_type_no_symmetry)
CALL dbcsr_init_p(buf1_dbcsr)
CALL dbcsr_create(buf1_dbcsr, template=ref_matrix, &
matrix_type=dbcsr_type_no_symmetry)
CALL dbcsr_init_p(buf2_dbcsr)
CALL dbcsr_create(buf2_dbcsr, template=ref_matrix, &
matrix_type=dbcsr_type_no_symmetry)
DO ispin = 1, nspins
ref_matrix => matrix_s(ispin, 1)%matrix
CALL copy_dbcsr_to_fm(ref_matrix, fm_ortho)
CALL cp_fm_cholesky_decompose(fm_ortho)
CALL cp_fm_triangular_invert(fm_ortho)
CALL cp_fm_set_all(fm_ks, 0.0_dp)
CALL cp_fm_to_fm_triangular(fm_ortho, fm_ks, "U")
CALL copy_fm_to_dbcsr(fm_ks, ortho_dbcsr)
CALL cp_fm_set_all(fm_ks, 0.0_dp)
! calculate ZHZ(T)
! calculate Z(T)HZ
CALL dbcsr_desymmetrize(matrix_ks(ispin, 1)%matrix, buf1_dbcsr)
CALL dbcsr_multiply("N", "N", 1.0_dp, buf1_dbcsr, ortho_dbcsr, &
0.0_dp, buf2_dbcsr, filter_eps=eps_filter)
CALL dbcsr_multiply("T", "N", 1.0_dp, ortho_dbcsr, buf2_dbcsr, &
0.0_dp, buf1_dbcsr, filter_eps=eps_filter)
! copy to fm format
CALL copy_dbcsr_to_fm(buf1_dbcsr, fm_ks)
CALL choose_eigv_solver(fm_ks, fm_mo, eigenvalues, info)
CPASSERT(info == 0)
! back transform of mos c = Z(T)*c
CALL copy_fm_to_dbcsr(fm_mo, buf1_dbcsr)
CALL dbcsr_multiply("N", "N", 1.0_dp, ortho_dbcsr, buf1_dbcsr, &
0.0_dp, buf2_dbcsr, filter_eps=eps_filter)
! density matrix
CALL dbcsr_set(matrix_p(ispin, 1)%matrix, 0.0_dp)
CALL dbcsr_multiply("N", "T", focc(ispin), buf2_dbcsr, buf2_dbcsr, &
1.0_dp, matrix_p(ispin, 1)%matrix, retain_sparsity=.TRUE., last_k=nmo(ispin))
END DO
CALL cp_fm_release(fm_ks)
CALL cp_fm_release(fm_mo)
CALL cp_fm_release(fm_ortho)
CALL dbcsr_release(ortho_dbcsr)
CALL dbcsr_release(buf1_dbcsr)
CALL dbcsr_release(buf2_dbcsr)
DEALLOCATE (ortho_dbcsr, buf1_dbcsr, buf2_dbcsr)
DEALLOCATE (eigenvalues)
CALL timestop(handle)
END SUBROUTINE ec_diag_solver
! **************************************************************************************************
!> \brief Calculate the energy correction
!> \param qs_env ...
!> \param ec_env ...
!> \param unit_nr ...
!> \author Creation (03.2014,JGH)
! **************************************************************************************************
SUBROUTINE ec_energy(qs_env, ec_env, unit_nr)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(energy_correction_type) :: ec_env
INTEGER, INTENT(IN) :: unit_nr
CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_energy', routineP = moduleN//':'//routineN
INTEGER :: handle, ispin, nspins
REAL(KIND=dp) :: eband, energy, trace
CALL timeset(routineN, handle)
! dispersion through pairpotentials
CALL calculate_dispersion_pairpot(qs_env, ec_env%dispersion_env, energy, .FALSE.)
ec_env%edispersion = ec_env%edispersion + energy
!deb ! gCP pairpotentials
!deb CALL calculate_gcp_pairpot(qs_env, ec_env%gcp_env, energy, .FALSE.)
!deb ec_env%edispersion = ec_env%edispersion+energy
SELECT CASE (ec_env%energy_functional)
CASE (kg_ec_functional_harris)
nspins = SIZE(ec_env%matrix_ks, 1)
eband = 0.0_dp
DO ispin = 1, nspins
!dbg
CALL dbcsr_dot(ec_env%matrix_p(ispin, 1)%matrix, ec_env%matrix_s(1, 1)%matrix, trace)
IF (unit_nr > 0) WRITE (unit_nr, '(T2,A,T16,F16.10)') 'Tr[PS] ', trace
!dbg
CALL dbcsr_dot(ec_env%matrix_ks(ispin, 1)%matrix, ec_env%matrix_p(ispin, 1)%matrix, trace)
eband = eband + trace
END DO
ec_env%eband = eband
ec_env%etotal = ec_env%eband + ec_env%ehartree + ec_env%exc - ec_env%vhxc + ec_env%edispersion
IF (unit_nr > 0) THEN
WRITE (unit_nr, '(T2,A,T16,F16.10)') "HF Etotal ", ec_env%etotal
WRITE (unit_nr, '(T2,A,T16,F16.10)') "Eband ", ec_env%eband
WRITE (unit_nr, '(T2,A,T16,F16.10)') "Ehartree ", ec_env%ehartree
WRITE (unit_nr, '(T2,A,T16,F16.10)') "Exc ", ec_env%exc
WRITE (unit_nr, '(T2,A,T16,F16.10)') "Evhxc ", ec_env%vhxc
WRITE (unit_nr, '(T2,A,T16,F16.10)') "Edisp ", ec_env%edispersion
END IF
CASE DEFAULT
CPASSERT(.FALSE.)
END SELECT
CALL timestop(handle)
END SUBROUTINE ec_energy
! **************************************************************************************************
!> \brief builds either the full neighborlist or neighborlists of molecular
!> \brief subsets, depending on parameter values
!> \param qs_env ...
!> \param ec_env ...
!> \par History
!> 2012.07 created [Martin Haeufel]
!> 2016.07 Adapted for Harris functional {JGH]
!> \author Martin Haeufel
! **************************************************************************************************
SUBROUTINE ec_build_neighborlist(qs_env, ec_env)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(energy_correction_type), POINTER :: ec_env
CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_build_neighborlist', &
routineP = moduleN//':'//routineN
INTEGER :: handle, ikind, nkind
LOGICAL :: gth_potential_present, &
sgp_potential_present, &
skip_load_balance_distributed
LOGICAL, ALLOCATABLE, DIMENSION(:) :: orb_present, ppl_present, ppnl_present
REAL(dp) :: subcells
REAL(dp), ALLOCATABLE, DIMENSION(:) :: orb_radius, ppl_radius, ppnl_radius
REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: pair_radius
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 :: distribution_1d
TYPE(distribution_2d_type), POINTER :: distribution_2d
TYPE(gth_potential_type), POINTER :: gth_potential
TYPE(gto_basis_set_type), POINTER :: basis_set
TYPE(local_atoms_type), ALLOCATABLE, DIMENSION(:) :: atom2d
TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(qs_kind_type), POINTER :: qs_kind
TYPE(qs_ks_env_type), POINTER :: ks_env
TYPE(sgp_potential_type), POINTER :: sgp_potential
CALL timeset(routineN, handle)
CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set)
CALL get_qs_kind_set(qs_kind_set, gth_potential_present=gth_potential_present, &
sgp_potential_present=sgp_potential_present)
nkind = SIZE(qs_kind_set)
ALLOCATE (orb_radius(nkind), ppl_radius(nkind), ppnl_radius(nkind))
ALLOCATE (orb_present(nkind), ppl_present(nkind), ppnl_present(nkind))
ALLOCATE (pair_radius(nkind, nkind))
ALLOCATE (atom2d(nkind))
CALL get_qs_env(qs_env, &
atomic_kind_set=atomic_kind_set, &
cell=cell, &
distribution_2d=distribution_2d, &
local_particles=distribution_1d, &
particle_set=particle_set, &
molecule_set=molecule_set)
CALL atom2d_build(atom2d, distribution_1d, distribution_2d, atomic_kind_set, &
molecule_set, .FALSE., particle_set)
DO ikind = 1, nkind
CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom2d(ikind)%list)
qs_kind => qs_kind_set(ikind)
CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type="HARRIS")
IF (ASSOCIATED(basis_set)) THEN
orb_present(ikind) = .TRUE.
CALL get_gto_basis_set(gto_basis_set=basis_set, kind_radius=orb_radius(ikind))
ELSE
orb_present(ikind) = .FALSE.
orb_radius(ikind) = 0.0_dp
END IF
CALL get_qs_kind(qs_kind, gth_potential=gth_potential, sgp_potential=sgp_potential)
IF (gth_potential_present .OR. sgp_potential_present) THEN
IF (ASSOCIATED(gth_potential)) THEN
CALL get_potential(potential=gth_potential, &
ppl_present=ppl_present(ikind), &
ppl_radius=ppl_radius(ikind), &
ppnl_present=ppnl_present(ikind), &
ppnl_radius=ppnl_radius(ikind))
ELSE IF (ASSOCIATED(sgp_potential)) THEN
CALL get_potential(potential=sgp_potential, &
ppl_present=ppl_present(ikind), &
ppl_radius=ppl_radius(ikind), &
ppnl_present=ppnl_present(ikind), &
ppnl_radius=ppnl_radius(ikind))
ELSE
ppl_present(ikind) = .FALSE.
ppl_radius(ikind) = 0.0_dp
ppnl_present(ikind) = .FALSE.
ppnl_radius(ikind) = 0.0_dp
END IF
END IF
END DO
CALL section_vals_val_get(qs_env%input, "DFT%SUBCELLS", r_val=subcells)
! overlap
CALL pair_radius_setup(orb_present, orb_present, orb_radius, orb_radius, pair_radius)
CALL build_neighbor_lists(ec_env%sab_orb, particle_set, atom2d, cell, pair_radius, &
subcells=subcells, nlname="sab_orb")
! pseudopotential
IF (gth_potential_present .OR. sgp_potential_present) THEN
IF (ANY(ppl_present)) THEN
CALL pair_radius_setup(orb_present, ppl_present, orb_radius, ppl_radius, pair_radius)
CALL build_neighbor_lists(ec_env%sac_ppl, particle_set, atom2d, cell, pair_radius, &
subcells=subcells, operator_type="ABC", nlname="sac_ppl")
END IF
IF (ANY(ppnl_present)) THEN
CALL pair_radius_setup(orb_present, ppnl_present, orb_radius, ppnl_radius, pair_radius)
CALL build_neighbor_lists(ec_env%sap_ppnl, particle_set, atom2d, cell, pair_radius, &
subcells=subcells, operator_type="ABBA", nlname="sap_ppnl")
END IF
END IF
! Release work storage
CALL atom2d_cleanup(atom2d)
DEALLOCATE (atom2d)
DEALLOCATE (orb_present, ppl_present, ppnl_present)
DEALLOCATE (orb_radius, ppl_radius, ppnl_radius)
DEALLOCATE (pair_radius)
! Task list
CALL get_qs_env(qs_env, ks_env=ks_env, dft_control=dft_control)
skip_load_balance_distributed = dft_control%qs_control%skip_load_balance_distributed
IF (ASSOCIATED(ec_env%task_list)) CALL deallocate_task_list(ec_env%task_list)
CALL allocate_task_list(ec_env%task_list)
CALL generate_qs_task_list(ks_env, ec_env%task_list, &
reorder_rs_grid_ranks=.FALSE., soft_valid=.FALSE., &
skip_load_balance_distributed=skip_load_balance_distributed, &
basis_type="HARRIS", sab_orb_external=ec_env%sab_orb)
CALL timestop(handle)
END SUBROUTINE ec_build_neighborlist
! **************************************************************************************************
END MODULE kg_energy_corrections

144
src/kg_environment.F
View file

@ -12,16 +12,11 @@
MODULE kg_environment
USE atomic_kind_types, ONLY: atomic_kind_type,&
get_atomic_kind
USE basis_set_container_types, ONLY: add_basis_set_to_container,&
remove_basis_from_container
USE basis_set_types, ONLY: copy_gto_basis_set,&
create_primitive_basis_set,&
get_gto_basis_set,&
USE basis_set_types, ONLY: get_gto_basis_set,&
gto_basis_set_type
USE bibliography, ONLY: Andermatt2016,&
cite_reference
USE cell_types, ONLY: cell_type
USE cp_control_types, ONLY: dft_control_type
USE cp_files, ONLY: close_file,&
open_file
USE cp_log_handling, ONLY: cp_get_default_logger,&
@ -32,13 +27,10 @@ MODULE kg_environment
USE distribution_2d_types, ONLY: distribution_2d_type
USE external_potential_types, ONLY: get_potential,&
local_potential_type
USE input_constants, ONLY: kg_ec_functional_harris,&
kg_tnadd_atomic,&
USE input_constants, ONLY: kg_tnadd_atomic,&
kg_tnadd_embed,&
kg_tnadd_embed_ri,&
kg_tnadd_none,&
xc_vdw_fun_nonloc,&
xc_vdw_fun_pairpot
kg_tnadd_none
USE input_section_types, ONLY: section_vals_get_subs_vals,&
section_vals_type,&
section_vals_val_get
@ -57,14 +49,9 @@ MODULE kg_environment
mp_max
USE molecule_types, ONLY: molecule_type
USE particle_types, ONLY: particle_type
USE qs_dispersion_nonloc, ONLY: qs_dispersion_nonloc_init
USE qs_dispersion_pairpot, ONLY: qs_dispersion_pairpot_init
USE qs_dispersion_types, ONLY: qs_dispersion_type
USE qs_dispersion_utils, ONLY: qs_dispersion_env_set
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
USE qs_grid_atom, ONLY: initialize_atomic_grid
USE qs_interactions, ONLY: init_interaction_radii_orb_basis
USE qs_kind_types, ONLY: get_qs_kind,&
qs_kind_type
USE qs_neighbor_list_types, ONLY: get_iterator_info,&
@ -140,23 +127,19 @@ CONTAINS
INTEGER :: handle, i, iatom, ib, ikind, iunit, n, &
na, natom, nbatch, nkind, np, nr
INTEGER, ALLOCATABLE, DIMENSION(:, :) :: bid
REAL(KIND=dp) :: eps_pgf_orb, load, radb, rmax
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
REAL(KIND=dp) :: load, radb, rmax
TYPE(cp_logger_type), POINTER :: logger
TYPE(cp_para_env_type), POINTER :: para_env
TYPE(dft_control_type), POINTER :: dft_control
TYPE(gto_basis_set_type), POINTER :: basis_set, harris_basis, lri_aux_basis
TYPE(gto_basis_set_type), POINTER :: lri_aux_basis
TYPE(integration_grid_type), POINTER :: ig_full, ig_mol
TYPE(qs_dispersion_type), POINTER :: dispersion_env
TYPE(qs_kind_type), POINTER :: qs_kind
TYPE(section_vals_type), POINTER :: lri_section, nl_section, pp_section, &
xc_section, xc_section_kg
TYPE(section_vals_type), POINTER :: lri_section
CALL timeset(routineN, handle)
CALL cite_reference(Andermatt2016)
NULLIFY (atomic_kind_set, dispersion_env, para_env)
NULLIFY (para_env)
NULLIFY (kg_env%sab_orb_full)
NULLIFY (kg_env%sac_kin)
NULLIFY (kg_env%subset_of_mol)
@ -167,15 +150,6 @@ CONTAINS
NULLIFY (kg_env%int_grid_molecules)
NULLIFY (kg_env%int_grid_full)
NULLIFY (kg_env%lri_density)
NULLIFY (kg_env%ec_env%sab_orb, kg_env%ec_env%sac_ppl, kg_env%ec_env%sap_ppnl)
NULLIFY (kg_env%ec_env%matrix_ks, kg_env%ec_env%matrix_h, kg_env%ec_env%matrix_s)
NULLIFY (kg_env%ec_env%matrix_t, kg_env%ec_env%matrix_p)
NULLIFY (kg_env%ec_env%task_list)
NULLIFY (kg_env%ec_env%mao_coef)
NULLIFY (kg_env%ec_env%dispersion_env)
NULLIFY (kg_env%ec_env%xc_section)
kg_env%ec_env%mao = .FALSE.
kg_env%nsubsets = 0
@ -184,110 +158,6 @@ CONTAINS
! get method for nonadditive kinetic energy embedding potential
CALL section_vals_val_get(input, "DFT%KG_METHOD%TNADD_METHOD", i_val=kg_env%tnadd_method)
!
CALL section_vals_val_get(input, "DFT%KG_METHOD%ENERGY_CORRECTION%_SECTION_PARAMETERS_", &
l_val=kg_env%energy_correction)
IF (kg_env%energy_correction) THEN
CALL section_vals_val_get(input, "DFT%KG_METHOD%ENERGY_CORRECTION%ALGORITHM", &
i_val=kg_env%ec_env%ks_solver)
CALL section_vals_val_get(input, "DFT%KG_METHOD%ENERGY_CORRECTION%ENERGY_FUNCTIONAL", &
i_val=kg_env%ec_env%energy_functional)
CALL section_vals_val_get(input, "DFT%KG_METHOD%ENERGY_CORRECTION%FACTORIZATION", &
i_val=kg_env%ec_env%factorization)
CALL section_vals_val_get(input, "DFT%KG_METHOD%ENERGY_CORRECTION%EPS_DEFAULT", &
r_val=kg_env%ec_env%eps_default)
CALL section_vals_val_get(input, "DFT%KG_METHOD%ENERGY_CORRECTION%HARRIS_BASIS", &
c_val=kg_env%ec_env%basis)
CALL section_vals_val_get(input, "DFT%KG_METHOD%ENERGY_CORRECTION%MAO", &
l_val=kg_env%ec_env%mao)
CALL section_vals_val_get(input, "DFT%KG_METHOD%ENERGY_CORRECTION%MAO_MAX_ITER", &
i_val=kg_env%ec_env%mao_max_iter)
CALL section_vals_val_get(input, "DFT%KG_METHOD%ENERGY_CORRECTION%MAO_EPS_GRAD", &
r_val=kg_env%ec_env%mao_eps_grad)
! set basis
nkind = SIZE(qs_kind_set)
CALL uppercase(kg_env%ec_env%basis)
SELECT CASE (kg_env%ec_env%basis)
CASE ("ORBITAL")
DO ikind = 1, nkind
qs_kind => qs_kind_set(ikind)
CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type="ORB")
IF (ASSOCIATED(basis_set)) THEN
NULLIFY (harris_basis)
CALL get_qs_kind(qs_kind=qs_kind, basis_set=harris_basis, basis_type="HARRIS")
IF (ASSOCIATED(harris_basis)) THEN
CALL remove_basis_from_container(qs_kind%basis_sets, basis_type="HARRIS")
END IF
NULLIFY (harris_basis)
CALL copy_gto_basis_set(basis_set, harris_basis)
CALL add_basis_set_to_container(qs_kind%basis_sets, harris_basis, "HARRIS")
END IF
END DO
CASE ("PRIMITIVE")
DO ikind = 1, nkind
qs_kind => qs_kind_set(ikind)
CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type="ORB")
IF (ASSOCIATED(basis_set)) THEN
NULLIFY (harris_basis)
CALL get_qs_kind(qs_kind=qs_kind, basis_set=harris_basis, basis_type="HARRIS")
IF (ASSOCIATED(harris_basis)) THEN
CALL remove_basis_from_container(qs_kind%basis_sets, basis_type="HARRIS")
END IF
NULLIFY (harris_basis)
CALL create_primitive_basis_set(basis_set, harris_basis)
CALL get_qs_env(qs_env, dft_control=dft_control)
eps_pgf_orb = dft_control%qs_control%eps_pgf_orb
CALL init_interaction_radii_orb_basis(harris_basis, eps_pgf_orb)
harris_basis%kind_radius = basis_set%kind_radius
CALL add_basis_set_to_container(qs_kind%basis_sets, harris_basis, "HARRIS")
END IF
END DO
CASE ("HARRIS")
DO ikind = 1, nkind
qs_kind => qs_kind_set(ikind)
NULLIFY (harris_basis)
CALL get_qs_kind(qs_kind=qs_kind, basis_set=harris_basis, basis_type="HARRIS")
IF (.NOT. ASSOCIATED(harris_basis)) THEN
CPWARN("Harris Basis not defined for all types of atoms.")
END IF
END DO
CASE DEFAULT
CPABORT("Unknown KG energy correction basis")
END SELECT
! set functional
SELECT CASE (kg_env%ec_env%energy_functional)
CASE (kg_ec_functional_harris)
kg_env%ec_env%ec_name = "Harris"
CASE DEFAULT
CPABORT("unknown kg energy correction")
END SELECT
! select the XC section
NULLIFY (xc_section, xc_section_kg)
xc_section => section_vals_get_subs_vals(input, "DFT%XC")
xc_section_kg => section_vals_get_subs_vals(input, "DFT%KG_METHOD%ENERGY_CORRECTION%XC")
IF (ASSOCIATED(xc_section_kg)) THEN
kg_env%ec_env%xc_section => xc_section_kg
ELSE
kg_env%ec_env%xc_section => xc_section
END IF
! dispersion
ALLOCATE (dispersion_env)
NULLIFY (xc_section)
xc_section => kg_env%ec_env%xc_section
CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, para_env=para_env)
CALL qs_dispersion_env_set(dispersion_env, xc_section)
IF (dispersion_env%type == xc_vdw_fun_pairpot) THEN
NULLIFY (pp_section)
pp_section => section_vals_get_subs_vals(xc_section, "VDW_POTENTIAL%PAIR_POTENTIAL")
CALL qs_dispersion_pairpot_init(atomic_kind_set, qs_kind_set, dispersion_env, pp_section, para_env)
ELSE IF (dispersion_env%type == xc_vdw_fun_nonloc) THEN
NULLIFY (nl_section)
nl_section => section_vals_get_subs_vals(xc_section, "VDW_POTENTIAL%NON_LOCAL")
CALL qs_dispersion_nonloc_init(dispersion_env, para_env)
END IF
kg_env%ec_env%dispersion_env => dispersion_env
END IF
SELECT CASE (kg_env%tnadd_method)
CASE (kg_tnadd_embed, kg_tnadd_embed_ri)
! kinetic energy functional

29
src/kg_environment_types.F
View file

@ -22,8 +22,7 @@ MODULE kg_environment_types
lri_env_release,&
lri_environment_type
USE molecule_types, ONLY: molecule_type
USE qs_dispersion_types, ONLY: qs_dispersion_release,&
qs_dispersion_type
USE qs_dispersion_types, ONLY: qs_dispersion_type
USE qs_grid_atom, ONLY: atom_integration_grid_type,&
deallocate_atom_int_grid
USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type,&
@ -115,9 +114,6 @@ MODULE kg_environment_types
INTEGER :: coloring_method
!
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: tnadd_mat
!
LOGICAL :: energy_correction
TYPE(energy_correction_type) :: ec_env
! LRI
TYPE(lri_environment_type), POINTER :: lri_env
TYPE(lri_density_type), POINTER :: lri_density
@ -179,29 +175,6 @@ CONTAINS
CALL deallocate_intgrid(kg_env%int_grid_full)
END IF
! energy correction
IF (kg_env%energy_correction) THEN
! neighbor lists
CALL release_neighbor_list_sets(kg_env%ec_env%sab_orb)
CALL release_neighbor_list_sets(kg_env%ec_env%sac_ppl)
CALL release_neighbor_list_sets(kg_env%ec_env%sap_ppnl)
! operator matrices
IF (ASSOCIATED(kg_env%ec_env%matrix_ks)) CALL dbcsr_deallocate_matrix_set(kg_env%ec_env%matrix_ks)
IF (ASSOCIATED(kg_env%ec_env%matrix_h)) CALL dbcsr_deallocate_matrix_set(kg_env%ec_env%matrix_h)
IF (ASSOCIATED(kg_env%ec_env%matrix_s)) CALL dbcsr_deallocate_matrix_set(kg_env%ec_env%matrix_s)
IF (ASSOCIATED(kg_env%ec_env%matrix_t)) CALL dbcsr_deallocate_matrix_set(kg_env%ec_env%matrix_t)
IF (ASSOCIATED(kg_env%ec_env%matrix_p)) CALL dbcsr_deallocate_matrix_set(kg_env%ec_env%matrix_p)
IF (ASSOCIATED(kg_env%ec_env%task_list)) THEN
CALL deallocate_task_list(kg_env%ec_env%task_list)
END IF
! reduced basis
IF (ASSOCIATED(kg_env%ec_env%mao_coef)) CALL dbcsr_deallocate_matrix_set(kg_env%ec_env%mao_coef)
! dispersion environment
IF (ASSOCIATED(kg_env%ec_env%dispersion_env)) THEN
CALL qs_dispersion_release(kg_env%ec_env%dispersion_env)
END IF
END IF
DEALLOCATE (kg_env)
CALL timestop(handle)

4
src/pw/pw_poisson_methods.F
View file

@ -297,11 +297,7 @@ CONTAINS
in_space=RECIPROCALSPACE)
CALL pw_copy(rhog, tmpg)
END IF
!$OMP PARALLEL DEFAULT(NONE) SHARED(rhog, poisson_env)
!$OMP WORKSHARE
rhog%cc(:) = rhog%cc(:)*poisson_env%green_fft%influence_fn%cc(:)
!$OMP END WORKSHARE
!$OMP END PARALLEL
IF (PRESENT(vhartree)) THEN
CALL pw_transfer(rhog, vhartree)
IF (PRESENT(ehartree)) THEN

24
src/qs_core_hamiltonian.F
View file

@ -287,19 +287,17 @@ CONTAINS
sab_nl=sab_orb, calculate_forces=.TRUE., &
matrixkp_p=matrix_p, &
eps_filter=eps_filter)
IF (calculate_forces) THEN
! *** If LSD, then recover alpha density and beta density ***
! *** from the total density (1) and the spin density (2) ***
! *** The W matrix is neglected, since it will be destroyed ***
! *** in the calling force routine after leaving this routine ***
IF (SIZE(matrix_p, 1) == 2) THEN
DO img = 1, nimages
CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
alpha_scalar=0.5_dp, beta_scalar=0.5_dp)
CALL dbcsr_add(matrix_p(2, img)%matrix, matrix_p(1, img)%matrix, &
alpha_scalar=-1.0_dp, beta_scalar=1.0_dp)
END DO
END IF
! *** If LSD, then recover alpha density and beta density ***
! *** from the total density (1) and the spin density (2) ***
! *** The W matrix is neglected, since it will be destroyed ***
! *** in the calling force routine after leaving this routine ***
IF (SIZE(matrix_p, 1) == 2) THEN
DO img = 1, nimages
CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
alpha_scalar=0.5_dp, beta_scalar=0.5_dp)
CALL dbcsr_add(matrix_p(2, img)%matrix, matrix_p(1, img)%matrix, &
alpha_scalar=-1.0_dp, beta_scalar=1.0_dp)
END DO
END IF
ELSE
! S matrix

10
src/qs_energy.F
View file

@ -13,7 +13,7 @@ MODULE qs_energy
USE almo_scf, ONLY: almo_entry_scf
USE cp_control_types, ONLY: dft_control_type
USE dm_ls_scf, ONLY: ls_scf
USE kg_energy_corrections, ONLY: kg_energy_correction
USE energy_corrections, ONLY: energy_correction
USE lri_environment_methods, ONLY: lri_print_stat
USE qs_energy_init, ONLY: qs_energies_init
USE qs_energy_types, ONLY: qs_energy_type
@ -93,11 +93,6 @@ CONTAINS
END IF
IF (dft_control%qs_control%do_kg) THEN
! Check for energy correction
CALL kg_energy_correction(qs_env, calculate_forces=my_calc_forces)
END IF
IF (PRESENT(consistent_energies)) THEN
IF (consistent_energies) THEN
CALL qs_ks_update_qs_env(qs_env, calculate_forces=.FALSE., just_energy=.TRUE.)
@ -108,6 +103,9 @@ CONTAINS
END IF
END IF
! Check for energy correction
CALL energy_correction(qs_env, ec_init=.TRUE., calculate_forces=.FALSE.)
CALL qs_energies_properties(qs_env)
IF (dft_control%qs_control%lrigpw) THEN

2
src/qs_energy_types.F
View file

@ -56,6 +56,7 @@ MODULE qs_energy_types
ktS, & ! electronic entropic contribution
efermi, & ! Fermi energy
dftb3, & ! DFTB 3rd order correction
nonscf_correction, & ! e.g. Harris correction
mp2, &
! single excitations correction for all
! non-scf orbital(density) corrections
@ -188,6 +189,7 @@ CONTAINS
qs_energy%surf_dipole = 0.0_dp
qs_energy%total = 0.0_dp
qs_energy%singles_corr = 0.0_dp
qs_energy%nonscf_correction = 0.0_dp
IF (.NOT. ASSOCIATED(qs_energy%ddapc_restraint)) THEN
ALLOCATE (qs_energy%ddapc_restraint(1))
END IF

10
src/qs_environment.F
View file

@ -63,6 +63,8 @@ MODULE qs_environment
distribution_1d_type
USE distribution_methods, ONLY: distribute_molecules_1d
USE dm_ls_scf_create, ONLY: ls_scf_create
USE ec_env_types, ONLY: energy_correction_type
USE ec_environment, ONLY: ec_env_create
USE et_coupling_types, ONLY: et_coupling_create
USE ewald_environment_types, ONLY: ewald_env_create,&
ewald_env_get,&
@ -257,6 +259,7 @@ 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(energy_correction_type), POINTER :: ec_env
TYPE(kpoint_type), POINTER :: kpoints
TYPE(lri_environment_type), POINTER :: lri_env
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
@ -424,6 +427,13 @@ CONTAINS
CALL kg_env_create(qs_env, qs_env%kg_env, qs_kind_set, qs_env%input)
END IF
dft_section => section_vals_get_subs_vals(qs_env%input, "DFT")
CALL section_vals_val_get(dft_section, "ENERGY_CORRECTION%_SECTION_PARAMETERS_", &
l_val=qs_env%energy_correction)
NULLIFY (ec_env)
CALL ec_env_create(qs_env, ec_env, dft_section)
CALL set_qs_env(qs_env, ec_env=ec_env)
et_coupling_section => section_vals_get_subs_vals(qs_env%input, &
"PROPERTIES%ET_COUPLING")
CALL section_vals_get(et_coupling_section, explicit=do_et)

33
src/qs_environment_types.F
View file

@ -39,6 +39,8 @@ MODULE qs_environment_types
USE distribution_2d_types, ONLY: distribution_2d_type
USE dm_ls_scf_types, ONLY: ls_scf_env_type,&
ls_scf_release
USE ec_env_types, ONLY: ec_env_release,&
energy_correction_type
USE et_coupling_types, ONLY: et_coupling_release,&
et_coupling_type
USE ewald_environment_types, ONLY: ewald_env_release,&
@ -136,7 +138,8 @@ MODULE qs_environment_types
USE qs_scf_types, ONLY: qs_scf_env_type,&
scf_env_release,&
scf_env_retain
USE qs_subsys_types, ONLY: qs_subsys_type
USE qs_subsys_types, ONLY: qs_subsys_set,&
qs_subsys_type
USE qs_wf_history_types, ONLY: qs_wf_history_type,&
wfi_release,&
wfi_retain
@ -277,6 +280,9 @@ MODULE qs_environment_types
! LRI
TYPE(lri_environment_type), POINTER :: lri_env
TYPE(lri_density_type), POINTER :: lri_density
! Energy correction
LOGICAL :: energy_correction
TYPE(energy_correction_type), POINTER :: ec_env
! Empirical dispersion
TYPE(qs_dispersion_type), POINTER :: dispersion_env
! Empirical geometrical BSSE correction
@ -459,6 +465,7 @@ CONTAINS
!> \param admm_dm ...
!> \param lri_env ...
!> \param lri_density ...
!> \param ec_env ...
!> \param dispersion_env ...
!> \param gcp_env ...
!> \param vee ...
@ -512,7 +519,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, admm_dm, &
lri_env, lri_density, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, &
lri_env, lri_density, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, &
mp2_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, &
gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, rhs)
@ -625,6 +632,7 @@ CONTAINS
TYPE(admm_dm_type), OPTIONAL, POINTER :: admm_dm
TYPE(lri_environment_type), OPTIONAL, POINTER :: lri_env
TYPE(lri_density_type), OPTIONAL, POINTER :: lri_density
TYPE(energy_correction_type), OPTIONAL, POINTER :: ec_env
TYPE(qs_dispersion_type), OPTIONAL, POINTER :: dispersion_env
TYPE(qs_gcp_type), OPTIONAL, POINTER :: gcp_env
TYPE(pw_p_type), OPTIONAL, POINTER :: vee
@ -699,6 +707,7 @@ CONTAINS
IF (PRESENT(admm_env)) admm_env => qs_env%admm_env
IF (PRESENT(lri_env)) lri_env => qs_env%lri_env
IF (PRESENT(lri_density)) lri_density => qs_env%lri_density
IF (PRESENT(ec_env)) ec_env => qs_env%ec_env
IF (PRESENT(dispersion_env)) dispersion_env => qs_env%dispersion_env
IF (PRESENT(gcp_env)) gcp_env => qs_env%gcp_env
IF (PRESENT(run_rtp)) run_rtp = qs_env%run_rtp
@ -915,6 +924,7 @@ CONTAINS
NULLIFY (qs_env%admm_env)
NULLIFY (qs_env%efield)
NULLIFY (qs_env%lri_env)
NULLIFY (qs_env%ec_env)
NULLIFY (qs_env%lri_density)
NULLIFY (qs_env%gcp_env)
NULLIFY (qs_env%rtp)
@ -1018,10 +1028,12 @@ CONTAINS
!> \param transport_env ...
!> \param lri_env ...
!> \param lri_density ...
!> \param ec_env ...
!> \param dispersion_env ...
!> \param gcp_env ...
!> \param mp2_env ...
!> \param kg_env ...
!> \param force ...
!> \param kpoints ...
!> \param WannierCentres ...
!> \param almo_scf_env ...
@ -1047,7 +1059,8 @@ CONTAINS
linres_control, xas_env, cp_ddapc_env, cp_ddapc_ewald, &
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, dispersion_env, gcp_env, mp2_env, kg_env, &
do_transport, transport_env, lri_env, lri_density, ec_env, dispersion_env, &
gcp_env, mp2_env, kg_env, force, &
kpoints, WannierCentres, almo_scf_env, gradient_history, variable_history, embed_pot, &
spin_embed_pot, polar_env, rhs)
@ -1108,10 +1121,13 @@ CONTAINS
TYPE(transport_env_type), OPTIONAL, POINTER :: transport_env
TYPE(lri_environment_type), OPTIONAL, POINTER :: lri_env
TYPE(lri_density_type), OPTIONAL, POINTER :: lri_density
TYPE(energy_correction_type), OPTIONAL, POINTER :: ec_env
TYPE(qs_dispersion_type), OPTIONAL, POINTER :: dispersion_env
TYPE(qs_gcp_type), OPTIONAL, POINTER :: gcp_env
TYPE(mp2_type), OPTIONAL, POINTER :: mp2_env
TYPE(kg_environment_type), OPTIONAL, POINTER :: kg_env
TYPE(qs_force_type), DIMENSION(:), OPTIONAL, &
POINTER :: force
TYPE(kpoint_type), OPTIONAL, POINTER :: kpoints
TYPE(wannier_centres_type), DIMENSION(:), &
OPTIONAL, POINTER :: WannierCentres
@ -1123,6 +1139,8 @@ CONTAINS
CHARACTER(len=*), PARAMETER :: routineN = 'set_qs_env', routineP = moduleN//':'//routineN
TYPE(qs_subsys_type), POINTER :: subsys
! CPPrecondition(ASSOCIATED(qs_env),cp_failure_level,routineP,failure)
CPASSERT(qs_env%ref_count > 0)
@ -1290,6 +1308,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(ec_env)) qs_env%ec_env => ec_env
IF (PRESENT(dispersion_env)) qs_env%dispersion_env => dispersion_env
IF (PRESENT(gcp_env)) qs_env%gcp_env => gcp_env
IF (PRESENT(WannierCentres)) qs_env%WannierCentres => WannierCentres
@ -1298,6 +1317,11 @@ CONTAINS
! Resp charges
IF (PRESENT(rhs)) qs_env%rhs => rhs
IF (PRESENT(force)) THEN
CALL get_qs_env(qs_env, subsys=subsys)
CALL qs_subsys_set(subsys, force=force)
END IF
END SUBROUTINE set_qs_env
! **************************************************************************************************
@ -1511,6 +1535,9 @@ CONTAINS
IF (ASSOCIATED(qs_env%lri_density)) THEN
CALL lri_density_release(qs_env%lri_density)
END IF
IF (ASSOCIATED(qs_env%ec_env)) THEN
CALL ec_env_release(qs_env%ec_env)
END IF
IF (ASSOCIATED(qs_env%mp2_env)) THEN
CALL mp2_env_release(qs_env%mp2_env)
END IF

8
src/qs_force.F
View file

@ -31,6 +31,7 @@ MODULE qs_force
dbcsr_set
USE dft_plus_u, ONLY: plus_u
USE efield_utils, ONLY: calculate_ecore_efield
USE energy_corrections, ONLY: energy_correction
USE input_constants, ONLY: do_admm_purify_none
USE input_section_types, ONLY: section_vals_get_subs_vals,&
section_vals_type,&
@ -380,7 +381,12 @@ CONTAINS
CALL calc_mixed_overlap_force(qs_env)
END IF
! *** replicate forces ***
! Energy_correction
CALL energy_correction(qs_env, ec_init=.FALSE., calculate_forces=.TRUE.)
! replicate forces (get current pointer)
NULLIFY (force)
CALL get_qs_env(qs_env=qs_env, force=force)
CALL replicate_qs_force(force, para_env)
DO iatom = 1, natom

209
src/qs_force_types.F
View file

@ -56,6 +56,10 @@ MODULE qs_force_types
add_qs_force, &
deallocate_qs_force, &
replicate_qs_force, &
sum_qs_force, &
get_qs_force, &
put_qs_force, &
total_qs_force, &
zero_qs_force
CONTAINS
@ -289,6 +293,76 @@ CONTAINS
END SUBROUTINE zero_qs_force
! **************************************************************************************************
!> \brief Sum up two qs_force entities qs_force_out = qs_force_out + qs_force_in
!> \param qs_force_out ...
!> \param qs_force_in ...
!> \author JGH
! **************************************************************************************************
SUBROUTINE sum_qs_force(qs_force_out, qs_force_in)
TYPE(qs_force_type), DIMENSION(:), POINTER :: qs_force_out, qs_force_in
CHARACTER(len=*), PARAMETER :: routineN = 'sum_qs_force', routineP = moduleN//':'//routineN
INTEGER :: ikind
CPASSERT(ASSOCIATED(qs_force_out))
CPASSERT(ASSOCIATED(qs_force_in))
DO ikind = 1, SIZE(qs_force_out)
qs_force_out(ikind)%all_potential(:, :) = qs_force_out(ikind)%all_potential(:, :) + &
qs_force_in(ikind)%all_potential(:, :)
qs_force_out(ikind)%core_overlap(:, :) = qs_force_out(ikind)%core_overlap(:, :) + &
qs_force_in(ikind)%core_overlap(:, :)
qs_force_out(ikind)%gth_ppl(:, :) = qs_force_out(ikind)%gth_ppl(:, :) + &
qs_force_in(ikind)%gth_ppl(:, :)
qs_force_out(ikind)%gth_nlcc(:, :) = qs_force_out(ikind)%gth_nlcc(:, :) + &
qs_force_in(ikind)%gth_nlcc(:, :)
qs_force_out(ikind)%gth_ppnl(:, :) = qs_force_out(ikind)%gth_ppnl(:, :) + &
qs_force_in(ikind)%gth_ppnl(:, :)
qs_force_out(ikind)%kinetic(:, :) = qs_force_out(ikind)%kinetic(:, :) + &
qs_force_in(ikind)%kinetic(:, :)
qs_force_out(ikind)%overlap(:, :) = qs_force_out(ikind)%overlap(:, :) + &
qs_force_in(ikind)%overlap(:, :)
qs_force_out(ikind)%overlap_admm(:, :) = qs_force_out(ikind)%overlap_admm(:, :) + &
qs_force_in(ikind)%overlap_admm(:, :)
qs_force_out(ikind)%rho_core(:, :) = qs_force_out(ikind)%rho_core(:, :) + &
qs_force_in(ikind)%rho_core(:, :)
qs_force_out(ikind)%rho_elec(:, :) = qs_force_out(ikind)%rho_elec(:, :) + &
qs_force_in(ikind)%rho_elec(:, :)
qs_force_out(ikind)%rho_lri_elec(:, :) = qs_force_out(ikind)%rho_lri_elec(:, :) + &
qs_force_in(ikind)%rho_lri_elec(:, :)
qs_force_out(ikind)%vhxc_atom(:, :) = qs_force_out(ikind)%vhxc_atom(:, :) + &
qs_force_in(ikind)%vhxc_atom(:, :)
qs_force_out(ikind)%g0s_Vh_elec(:, :) = qs_force_out(ikind)%g0s_Vh_elec(:, :) + &
qs_force_in(ikind)%g0s_Vh_elec(:, :)
qs_force_out(ikind)%repulsive(:, :) = qs_force_out(ikind)%repulsive(:, :) + &
qs_force_in(ikind)%repulsive(:, :)
qs_force_out(ikind)%dispersion(:, :) = qs_force_out(ikind)%dispersion(:, :) + &
qs_force_in(ikind)%dispersion(:, :)
qs_force_out(ikind)%gcp(:, :) = qs_force_out(ikind)%gcp(:, :) + &
qs_force_in(ikind)%gcp(:, :)
qs_force_out(ikind)%other(:, :) = qs_force_out(ikind)%other(:, :) + &
qs_force_in(ikind)%other(:, :)
qs_force_out(ikind)%fock_4c(:, :) = qs_force_out(ikind)%fock_4c(:, :) + &
qs_force_in(ikind)%fock_4c(:, :)
qs_force_out(ikind)%ehrenfest(:, :) = qs_force_out(ikind)%ehrenfest(:, :) + &
qs_force_in(ikind)%ehrenfest(:, :)
qs_force_out(ikind)%efield(:, :) = qs_force_out(ikind)%efield(:, :) + &
qs_force_in(ikind)%efield(:, :)
qs_force_out(ikind)%eev(:, :) = qs_force_out(ikind)%eev(:, :) + &
qs_force_in(ikind)%eev(:, :)
qs_force_out(ikind)%mp2_non_sep(:, :) = qs_force_out(ikind)%mp2_non_sep(:, :) + &
qs_force_in(ikind)%mp2_non_sep(:, :)
qs_force_out(ikind)%mp2_sep(:, :) = qs_force_out(ikind)%mp2_sep(:, :) + &
qs_force_in(ikind)%mp2_sep(:, :)
qs_force_out(ikind)%total(:, :) = qs_force_out(ikind)%total(:, :) + &
qs_force_in(ikind)%total(:, :)
END DO
END SUBROUTINE sum_qs_force
! **************************************************************************************************
!> \brief Replicate and sum up the force
!> \param qs_force ...
@ -399,4 +473,139 @@ CONTAINS
END SUBROUTINE add_qs_force
! **************************************************************************************************
!> \brief Put force to a force_type variable.
!> \param force Input force, dimension (3,natom)
!> \param qs_force The force type variable to be used
!> \param forcetype ...
!> \param atomic_kind_set ...
!> \par History
!> 09.2019 JGH
!> \author JGH
! **************************************************************************************************
SUBROUTINE put_qs_force(force, qs_force, forcetype, atomic_kind_set)
REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: force
TYPE(qs_force_type), DIMENSION(:), POINTER :: qs_force
CHARACTER(LEN=*), INTENT(IN) :: forcetype
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
CHARACTER(len=*), PARAMETER :: routineN = 'put_qs_force', routineP = moduleN//':'//routineN
INTEGER :: ia, iatom, ikind, natom_kind
TYPE(atomic_kind_type), POINTER :: atomic_kind
! ------------------------------------------------------------------------
SELECT CASE (forcetype)
CASE ("dispersion")
DO ikind = 1, SIZE(atomic_kind_set, 1)
atomic_kind => atomic_kind_set(ikind)
CALL get_atomic_kind(atomic_kind=atomic_kind, natom=natom_kind)
DO ia = 1, natom_kind
iatom = atomic_kind%atom_list(ia)
qs_force(ikind)%dispersion(:, ia) = force(:, iatom)
END DO
END DO
CASE DEFAULT
CPABORT("")
END SELECT
END SUBROUTINE put_qs_force
! **************************************************************************************************
!> \brief Get force from a force_type variable.
!> \param force Input force, dimension (3,natom)
!> \param qs_force The force type variable to be used
!> \param forcetype ...
!> \param atomic_kind_set ...
!> \par History
!> 09.2019 JGH
!> \author JGH
! **************************************************************************************************
SUBROUTINE get_qs_force(force, qs_force, forcetype, atomic_kind_set)
REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: force
TYPE(qs_force_type), DIMENSION(:), POINTER :: qs_force
CHARACTER(LEN=*), INTENT(IN) :: forcetype
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
CHARACTER(len=*), PARAMETER :: routineN = 'get_qs_force', routineP = moduleN//':'//routineN
INTEGER :: ia, iatom, ikind, natom_kind
TYPE(atomic_kind_type), POINTER :: atomic_kind
! ------------------------------------------------------------------------
SELECT CASE (forcetype)
CASE ("dispersion")
DO ikind = 1, SIZE(atomic_kind_set, 1)
atomic_kind => atomic_kind_set(ikind)
CALL get_atomic_kind(atomic_kind=atomic_kind, natom=natom_kind)
DO ia = 1, natom_kind
iatom = atomic_kind%atom_list(ia)
force(:, iatom) = qs_force(ikind)%dispersion(:, ia)
END DO
END DO
CASE DEFAULT
CPABORT("")
END SELECT
END SUBROUTINE get_qs_force
! **************************************************************************************************
!> \brief Get current total force
!> \param force Input force, dimension (3,natom)
!> \param qs_force The force type variable to be used
!> \param atomic_kind_set ...
!> \par History
!> 09.2019 JGH
!> \author JGH
! **************************************************************************************************
SUBROUTINE total_qs_force(force, qs_force, atomic_kind_set)
REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: force
TYPE(qs_force_type), DIMENSION(:), POINTER :: qs_force
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
CHARACTER(len=*), PARAMETER :: routineN = 'total_qs_force', routineP = moduleN//':'//routineN
INTEGER :: ia, iatom, ikind, natom_kind
TYPE(atomic_kind_type), POINTER :: atomic_kind
! ------------------------------------------------------------------------
force(:, :) = 0.0_dp
DO ikind = 1, SIZE(atomic_kind_set, 1)
atomic_kind => atomic_kind_set(ikind)
CALL get_atomic_kind(atomic_kind=atomic_kind, natom=natom_kind)
DO ia = 1, natom_kind
iatom = atomic_kind%atom_list(ia)
force(:, iatom) = qs_force(ikind)%core_overlap(:, ia) + &
qs_force(ikind)%gth_ppl(:, ia) + &
qs_force(ikind)%gth_nlcc(:, ia) + &
qs_force(ikind)%gth_ppnl(:, ia) + &
qs_force(ikind)%all_potential(:, ia) + &
qs_force(ikind)%kinetic(:, ia) + &
qs_force(ikind)%overlap(:, ia) + &
qs_force(ikind)%overlap_admm(:, ia) + &
qs_force(ikind)%rho_core(:, ia) + &
qs_force(ikind)%rho_elec(:, ia) + &
qs_force(ikind)%rho_lri_elec(:, ia) + &
qs_force(ikind)%vhxc_atom(:, ia) + &
qs_force(ikind)%g0s_Vh_elec(:, ia) + &
qs_force(ikind)%fock_4c(:, ia) + &
qs_force(ikind)%mp2_non_sep(:, ia) + &
qs_force(ikind)%mp2_sep(:, ia) + &
qs_force(ikind)%repulsive(:, ia) + &
qs_force(ikind)%dispersion(:, ia) + &
qs_force(ikind)%gcp(:, ia) + &
qs_force(ikind)%ehrenfest(:, ia) + &
qs_force(ikind)%efield(:, ia) + &
qs_force(ikind)%eev(:, ia)
END DO
END DO
END SUBROUTINE total_qs_force
END MODULE qs_force_types

19
src/qs_integrate_potential_single.F
View file

@ -503,7 +503,7 @@ CONTAINS
TYPE(virial_type), POINTER :: virial
CALL timeset(routineN, handle)
NULLIFY (virial, atprop, dft_control)
NULLIFY (virial, force, atprop, dft_control)
CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
@ -553,16 +553,9 @@ CONTAINS
alpha_core_charge=alpha_core_charge, &
ccore_charge=ccore_charge)
IF (paw_atom) THEN
force(ikind)%rho_core(:, :) = 0.0_dp
CYCLE
END IF
IF (paw_atom) CYCLE
pab(1, 1) = -ccore_charge
IF (ASSOCIATED(force)) THEN
force(ikind)%rho_core = 0.0_dp
ENDIF
IF (alpha_core_charge == 0.0_dp .OR. pab(1, 1) == 0.0_dp) CYCLE
CALL reallocate(cores, 1, natom_of_kind)
@ -606,10 +599,8 @@ CONTAINS
my_virial_b=my_virial_b, use_subpatch=.TRUE., subpatch_pattern=0_int_8)
IF (ASSOCIATED(force)) THEN
force(ikind)%rho_core(:, iatom) = &
force(ikind)%rho_core(:, iatom) + force_a(:)
ENDIF
force(ikind)%rho_core(:, iatom) = force(ikind)%rho_core(:, iatom) + force_a(:)
END IF
IF (use_virial) THEN
virial%pv_virial = virial%pv_virial + my_virial_a
virial%pv_hartree = virial%pv_hartree + my_virial_a

8
src/qs_kpp1_env_methods.F
View file

@ -114,8 +114,8 @@ CONTAINS
ALLOCATE (kpp1_env)
NULLIFY (kpp1_env%v_rspace, kpp1_env%v_ao, kpp1_env%drho_r, &
kpp1_env%rho_set, &
kpp1_env%deriv_set, kpp1_env%spin_pot, kpp1_env%grad_pot, &
kpp1_env%rho_set, kpp1_env%deriv_set, kpp1_env%spin_pot, kpp1_env%grad_pot, &
kpp1_env%rho_set_admm, kpp1_env%deriv_set_admm, &
kpp1_env%ndiag_term)
kpp1_env%ref_count = 1
last_kpp1_id_nr = last_kpp1_id_nr + 1
@ -660,11 +660,9 @@ CONTAINS
END IF
xc_section => section_vals_get_subs_vals(input, "DFT%XC")
CALL section_vals_val_get(input, "DFT%EXCITATIONS", &
i_val=excitations)
CALL section_vals_val_get(input, "DFT%EXCITATIONS", i_val=excitations)
IF (excitations == tddfpt_excitations) THEN
xc_section => section_vals_get_subs_vals(input, "DFT%TDDFPT%XC")
!FM this check should already had happened and section made explicit, give an error?
CALL section_vals_get(xc_section, explicit=explicit)
IF (.NOT. explicit) THEN
xc_section => section_vals_get_subs_vals(input, "DFT%XC")

10
src/qs_kpp1_env_types.F
View file

@ -60,6 +60,8 @@ MODULE qs_kpp1_env_types
TYPE(pw_p_type), DIMENSION(:, :), POINTER :: drho_r
TYPE(xc_derivative_set_type), POINTER :: deriv_set
TYPE(xc_rho_set_type), POINTER :: rho_set
TYPE(xc_derivative_set_type), POINTER :: deriv_set_admm
TYPE(xc_rho_set_type), POINTER :: rho_set_admm
INTEGER, DIMENSION(:, :), POINTER :: spin_pot
LOGICAL, DIMENSION(:, :), POINTER :: grad_pot
LOGICAL, DIMENSION(:), POINTER :: ndiag_term
@ -121,6 +123,14 @@ CONTAINS
CALL xc_rho_set_release(kpp1_env%rho_set)
NULLIFY (kpp1_env%rho_set)
END IF
IF (ASSOCIATED(kpp1_env%deriv_set_admm)) THEN
CALL xc_dset_release(kpp1_env%deriv_set_admm)
NULLIFY (kpp1_env%deriv_set_admm)
END IF
IF (ASSOCIATED(kpp1_env%rho_set_admm)) THEN
CALL xc_rho_set_release(kpp1_env%rho_set_admm)
NULLIFY (kpp1_env%rho_set_admm)
END IF
IF (ASSOCIATED(kpp1_env%spin_pot)) THEN
DEALLOCATE (kpp1_env%spin_pot)
END IF

514
src/qs_linres_methods.F
View file

@ -11,12 +11,14 @@
!> \author MI
! **************************************************************************************************
MODULE qs_linres_methods
USE admm_types, ONLY: admm_type
USE atomic_kind_types, ONLY: atomic_kind_type,&
get_atomic_kind
USE cp_control_types, ONLY: dft_control_type
USE cp_dbcsr_operations, ONLY: cp_dbcsr_plus_fm_fm_t,&
cp_dbcsr_sm_fm_multiply,&
dbcsr_allocate_matrix_set
dbcsr_allocate_matrix_set,&
dbcsr_deallocate_matrix_set
USE cp_external_control, ONLY: external_control
USE cp_files, ONLY: close_file,&
open_file
@ -46,19 +48,23 @@ MODULE qs_linres_methods
USE cp_para_types, ONLY: cp_para_env_type
USE dbcsr_api, ONLY: dbcsr_checksum,&
dbcsr_copy,&
dbcsr_create,&
dbcsr_deallocate_matrix,&
dbcsr_filter,&
dbcsr_p_type,&
dbcsr_set,&
dbcsr_type
USE hartree_local_methods, ONLY: Vh_1c_gg_integrals
USE input_constants, ONLY: do_loc_none,&
op_loc_berry,&
ot_precond_none,&
ot_precond_solver_default,&
state_loc_all
USE hfx_energy_potential, ONLY: integrate_four_center
USE hfx_types, ONLY: hfx_type
USE input_constants, ONLY: &
do_admm_aux_exch_func_none, do_admm_basis_projection, do_admm_exch_scaling_none, &
do_admm_purify_none, do_loc_none, op_loc_berry, ot_precond_none, &
ot_precond_solver_default, state_loc_all
USE input_section_types, ONLY: section_get_ival,&
section_get_lval,&
section_get_rval,&
section_vals_get,&
section_vals_get_subs_vals,&
section_vals_type,&
section_vals_val_get
@ -87,7 +93,6 @@ MODULE qs_linres_methods
USE pw_poisson_types, ONLY: pw_poisson_type
USE pw_pool_types, ONLY: pw_pool_create_pw,&
pw_pool_give_back_pw,&
pw_pool_p_type,&
pw_pool_type
USE pw_types, ONLY: COMPLEXDATA1D,&
REALDATA3D,&
@ -97,6 +102,7 @@ MODULE qs_linres_methods
pw_p_type,&
pw_release,&
pw_retain
USE qs_collocate_density, ONLY: calculate_rho_elec
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
USE qs_gapw_densities, ONLY: prepare_gapw_den
@ -108,6 +114,7 @@ MODULE qs_linres_methods
qs_kind_type
USE qs_kpp1_env_types, ONLY: qs_kpp1_env_type
USE qs_ks_atom, ONLY: update_ks_atom
USE qs_ks_types, ONLY: qs_ks_env_type
USE qs_linres_types, ONLY: linres_control_type
USE qs_loc_methods, ONLY: qs_loc_driver
USE qs_loc_types, ONLY: get_qs_loc_env,&
@ -129,6 +136,7 @@ MODULE qs_linres_methods
qs_rho_type
USE qs_vxc_atom, ONLY: calculate_xc_2nd_deriv_atom
USE string_utilities, ONLY: xstring
USE task_list_types, ONLY: task_list_type
USE xc, ONLY: xc_calc_2nd_deriv,&
xc_prep_2nd_deriv,&
xc_vxc_pw_create
@ -151,6 +159,7 @@ MODULE qs_linres_methods
! *** Public subroutines ***
PUBLIC :: linres_localize, linres_solver
PUBLIC :: linres_write_restart, linres_read_restart
PUBLIC :: build_dm_response
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_linres_methods'
@ -379,7 +388,7 @@ CONTAINS
IF (iounit > 0) THEN
WRITE (iounit, "(/,T3,A,T16,A,T25,A,T38,A,T52,A,T72,A,/,T3,A)") &
"Iteration", "Method", "Restart", "Stepsize", "Convergence", "Time", &
REPEAT("-", 80)
REPEAT("-", 78)
ENDIF
!
! orthogonalize x with respect to the psi0
@ -707,7 +716,7 @@ CONTAINS
CALL get_qs_env(qs_env, rho=rho) ! that could be called before
CALL qs_rho_update_rho(rho, qs_env=qs_env) ! that could be called before
CALL apply_op_2(qs_env, p_env, c0, v, Av, chc)
CALL apply_op_2(qs_env, p_env, c0, v, Av)
ENDIF
@ -806,13 +815,12 @@ CONTAINS
!> \param c0 ...
!> \param v ...
!> \param Av ...
!> \param chc ...
! **************************************************************************************************
SUBROUTINE apply_op_2(qs_env, p_env, c0, v, Av, chc)
SUBROUTINE apply_op_2(qs_env, p_env, c0, v, Av)
!
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(qs_p_env_type), POINTER :: p_env
TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: c0, v, Av, chc
TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: c0, v, Av
CHARACTER(len=*), PARAMETER :: routineN = 'apply_op_2', routineP = moduleN//':'//routineN
@ -824,9 +832,11 @@ CONTAINS
ELSEIF (dft_control%qs_control%dftb) THEN
CPABORT("Linear response not available with DFTB")
ELSEIF (dft_control%qs_control%xtb) THEN
CALL apply_op_2_xtb(qs_env, p_env, c0, v, Av, chc)
CALL apply_op_2_xtb(qs_env, p_env, c0, v, Av)
ELSE
CALL apply_op_2_dft(qs_env, p_env, c0, v, Av, chc)
CALL apply_op_2_dft(qs_env, p_env, c0, v, Av)
CALL apply_hfx(qs_env, p_env, c0, v, Av)
CALL apply_xc_admm(qs_env, p_env, c0, v, Av)
END IF
END SUBROUTINE apply_op_2
@ -838,12 +848,11 @@ CONTAINS
!> \param c0 ...
!> \param v ...
!> \param Av ...
!> \param chc ...
! **************************************************************************************************
SUBROUTINE apply_op_2_dft(qs_env, p_env, c0, v, Av, chc)
SUBROUTINE apply_op_2_dft(qs_env, p_env, c0, v, Av)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(qs_p_env_type), POINTER :: p_env
TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: c0, v, Av, chc
TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: c0, v, Av
CHARACTER(len=*), PARAMETER :: routineN = 'apply_op_2_dft', routineP = moduleN//':'//routineN
REAL(KIND=dp), PARAMETER :: h = 0.001_dp
@ -857,6 +866,7 @@ CONTAINS
fac
REAL(KIND=dp), DIMENSION(3, 3) :: virial_xc
REAL(kind=dp), DIMENSION(:, :, :), POINTER :: rho3, rhoa, rhob
TYPE(admm_type), POINTER :: admm_env
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(cp_logger_type), POINTER :: logger
TYPE(cp_para_env_type), POINTER :: para_env
@ -873,19 +883,17 @@ CONTAINS
TYPE(pw_p_type), DIMENSION(:), POINTER :: rho1_g, rho1_g_pw, rho1_r, rho1_r_pw, rho_g, &
rho_r, tau, tau_pw, v_rspace_new, v_xc, vxc_rho_1, vxc_rho_2, vxc_rho_3, vxc_rho_4
TYPE(pw_poisson_type), POINTER :: poisson_env
TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
TYPE(qs_rho_type), POINTER :: rho, rho1, rho1_xc
TYPE(section_vals_type), POINTER :: input, scf_section, xc_fun_section, &
xc_section
TYPE(section_vals_type), POINTER :: input, xc_fun_section, xc_section
TYPE(xc_rho_cflags_type) :: needs
TYPE(xc_rho_set_type), POINTER :: rho1_set
CALL timeset(routineN, handle)
NULLIFY (auxbas_pw_pool, pw_pools, pw_env, v_rspace_new, &
NULLIFY (auxbas_pw_pool, pw_env, v_rspace_new, &
rho1_r, rho1_g_pw, tau_pw, v_xc, rho1_set, rho1_ao, rho_ao, &
poisson_env, input, scf_section, rho, dft_control, logger, rho1_g)
poisson_env, input, rho, dft_control, logger, rho1_g)
logger => cp_get_default_logger()
energy_hartree = 0.0_dp
@ -894,7 +902,6 @@ CONTAINS
CPASSERT(ASSOCIATED(c0))
CPASSERT(ASSOCIATED(v))
CPASSERT(ASSOCIATED(Av))
CPASSERT(ASSOCIATED(chc))
CPASSERT(ASSOCIATED(p_env%kpp1_env))
CPASSERT(ASSOCIATED(p_env%kpp1))
@ -932,15 +939,19 @@ CONTAINS
nspins = SIZE(p_env%kpp1)
lsd = (nspins == 2)
xc_section => section_vals_get_subs_vals(input, "DFT%XC")
scf_section => section_vals_get_subs_vals(input, "DFT%SCF")
IF (dft_control%do_admm) THEN
CALL get_qs_env(qs_env, admm_env=admm_env)
xc_section => admm_env%xc_section_primary
ELSE
xc_section => section_vals_get_subs_vals(input, "DFT%XC")
END IF
p_env%kpp1_env%iter = p_env%kpp1_env%iter + 1
! gets the tmp grids
CPASSERT(ASSOCIATED(pw_env))
CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
pw_pools=pw_pools, poisson_env=poisson_env)
poisson_env=poisson_env)
ALLOCATE (v_rspace_new(nspins))
CALL pw_pool_create_pw(auxbas_pw_pool, v_hartree_gspace%pw, &
use_data=COMPLEXDATA1D, &
@ -1323,9 +1334,9 @@ CONTAINS
Av(ispin)%matrix, &
ncol=ncol, alpha=1.0_dp, beta=1.0_dp)
ENDDO
!
CALL timestop(handle)
!
END SUBROUTINE apply_op_2_dft
! **************************************************************************************************
@ -1335,12 +1346,11 @@ CONTAINS
!> \param c0 ...
!> \param v ...
!> \param Av ...
!> \param chc ...
! **************************************************************************************************
SUBROUTINE apply_op_2_xtb(qs_env, p_env, c0, v, Av, chc)
SUBROUTINE apply_op_2_xtb(qs_env, p_env, c0, v, Av)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(qs_p_env_type), POINTER :: p_env
TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: c0, v, Av, chc
TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: c0, v, Av
CHARACTER(len=*), PARAMETER :: routineN = 'apply_op_2_xtb', routineP = moduleN//':'//routineN
@ -1370,7 +1380,6 @@ CONTAINS
CPASSERT(ASSOCIATED(c0))
CPASSERT(ASSOCIATED(v))
CPASSERT(ASSOCIATED(Av))
CPASSERT(ASSOCIATED(chc))
CPASSERT(ASSOCIATED(p_env%kpp1_env))
CPASSERT(ASSOCIATED(p_env%kpp1))
@ -1457,6 +1466,421 @@ CONTAINS
END SUBROUTINE apply_op_2_xtb
! **************************************************************************************************
!> \brief Update action of TDDFPT operator on trial vectors by adding exact-exchange term.
!> \param qs_env ...
!> \param p_env ...
!> \param c0 ...
!> \param v ...
!> \param Av ...
!> \par History
!> * 11.2019 adapted from tddfpt_apply_hfx
! **************************************************************************************************
SUBROUTINE apply_hfx(qs_env, p_env, c0, v, Av)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(qs_p_env_type), POINTER :: p_env
TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: c0, v, Av
CHARACTER(LEN=*), PARAMETER :: routineN = 'apply_hfx', routineP = moduleN//':'//routineN
INTEGER :: handle, ispin, nao, nao_aux, ncol, nspins
LOGICAL :: do_hfx
REAL(KIND=dp) :: alpha
TYPE(admm_type), POINTER :: admm_env
TYPE(cp_fm_type), POINTER :: tc0, tv
TYPE(cp_logger_type), POINTER :: logger
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, rho1_ao, work
TYPE(dft_control_type), POINTER :: dft_control
TYPE(section_vals_type), POINTER :: hfx_section, input
CALL timeset(routineN, handle)
logger => cp_get_default_logger()
CPASSERT(ASSOCIATED(c0))
CPASSERT(ASSOCIATED(v))
CPASSERT(ASSOCIATED(Av))
CALL get_qs_env(qs_env=qs_env, &
input=input, &
matrix_s=matrix_s, &
dft_control=dft_control)
nspins = dft_control%nspins
hfx_section => section_vals_get_subs_vals(input, "DFT%XC%HF")
CALL section_vals_get(hfx_section, explicit=do_hfx)
IF (do_hfx) THEN
IF (dft_control%do_admm) THEN
IF (dft_control%admm_control%purification_method /= do_admm_purify_none) THEN
CPABORT("ADMM: Linear Response needs purification_method=none")
END IF
IF (dft_control%admm_control%scaling_model /= do_admm_exch_scaling_none) THEN
CPABORT("ADMM: Linear Response needs scaling_model=none")
END IF
IF (dft_control%admm_control%method /= do_admm_basis_projection) THEN
CPABORT("ADMM: Linear Response needs admm_method=basis_projection")
END IF
!
CALL get_qs_env(qs_env, admm_env=admm_env)
CPASSERT(ASSOCIATED(admm_env%A))
CPASSERT(ASSOCIATED(admm_env%work_aux_orb))
CPASSERT(ASSOCIATED(admm_env%work_aux_orb2))
CALL cp_fm_get_info(admm_env%A, nrow_global=nao_aux)
CALL get_qs_env(qs_env=qs_env, matrix_s_aux_fit=matrix_s)
NULLIFY (work, rho1_ao)
CALL dbcsr_allocate_matrix_set(work, nspins)
CALL dbcsr_allocate_matrix_set(rho1_ao, nspins)
DO ispin = 1, nspins
ALLOCATE (work(ispin)%matrix, rho1_ao(ispin)%matrix)
CALL dbcsr_create(work(ispin)%matrix, template=matrix_s(1)%matrix)
CALL dbcsr_copy(work(ispin)%matrix, matrix_s(1)%matrix)
CALL dbcsr_set(work(ispin)%matrix, 0.0_dp)
CALL dbcsr_create(rho1_ao(ispin)%matrix, template=matrix_s(1)%matrix)
CALL dbcsr_copy(rho1_ao(ispin)%matrix, matrix_s(1)%matrix)
CALL dbcsr_set(rho1_ao(ispin)%matrix, 0.0_dp)
END DO
! P1 -> AUX BASIS
DO ispin = 1, nspins
CALL cp_fm_get_info(c0(ispin)%matrix, nrow_global=nao, ncol_global=ncol)
tv => admm_env%work_aux_orb
tc0 => admm_env%work_aux_orb2
CALL cp_gemm('N', 'N', nao_aux, ncol, nao, 1.0_dp, admm_env%A, &
v(ispin)%matrix, 0.0_dp, tv)
CALL cp_gemm('N', 'N', nao_aux, ncol, nao, 1.0_dp, admm_env%A, &
c0(ispin)%matrix, 0.0_dp, tc0)
CALL cp_dbcsr_plus_fm_fm_t(rho1_ao(ispin)%matrix, matrix_v=tv, matrix_g=tc0, &
ncol=ncol, alpha=1.0_dp)
CALL cp_dbcsr_plus_fm_fm_t(rho1_ao(ispin)%matrix, matrix_v=tc0, matrix_g=tv, &
ncol=ncol, alpha=1.0_dp)
ENDDO
ELSE
NULLIFY (work, rho1_ao)
CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s)
CALL dbcsr_allocate_matrix_set(work, nspins)
DO ispin = 1, nspins
ALLOCATE (work(ispin)%matrix)
CALL dbcsr_create(work(ispin)%matrix, template=matrix_s(1)%matrix)
CALL dbcsr_copy(work(ispin)%matrix, matrix_s(1)%matrix)
CALL dbcsr_set(work(ispin)%matrix, 0.0_dp)
END DO
rho1_ao => p_env%p1
END IF
CALL hfx_matrix(work, rho1_ao, qs_env, hfx_section)
alpha = 2.0_dp
IF (nspins == 2) alpha = 1.0_dp
IF (dft_control%do_admm) THEN
DO ispin = 1, nspins
CALL cp_fm_get_info(c0(ispin)%matrix, nrow_global=nao, ncol_global=ncol)
CALL cp_gemm('N', 'N', nao_aux, ncol, nao, 1.0_dp, admm_env%A, &
c0(ispin)%matrix, 0.0_dp, tc0)
CALL cp_dbcsr_sm_fm_multiply(work(ispin)%matrix, tc0, tv, &
ncol=ncol, alpha=alpha, beta=0.0_dp)
CALL cp_gemm('T', 'N', nao, ncol, nao_aux, 1.0_dp, admm_env%A, &
tv, 1.0_dp, Av(ispin)%matrix)
END DO
CALL dbcsr_deallocate_matrix_set(rho1_ao)
CALL dbcsr_deallocate_matrix_set(work)
ELSE
DO ispin = 1, nspins
CALL cp_fm_get_info(c0(ispin)%matrix, ncol_global=ncol)
CALL cp_dbcsr_sm_fm_multiply(work(ispin)%matrix, c0(ispin)%matrix, Av(ispin)%matrix, &
ncol=ncol, alpha=alpha, beta=1.0_dp)
END DO
CALL dbcsr_deallocate_matrix_set(work)
END IF
END IF
CALL timestop(handle)
END SUBROUTINE apply_hfx
! **************************************************************************************************
!> \brief Add the hfx contributions to the Hamiltonian
!>
!> \param matrix_ks ...
!> \param rho_ao ...
!> \param qs_env ...
!> \param hfx_sections ...
!> \note
!> Simplified version of subroutine hfx_ks_matrix()
! **************************************************************************************************
SUBROUTINE hfx_matrix(matrix_ks, rho_ao, qs_env, hfx_sections)
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, rho_ao
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(section_vals_type), POINTER :: hfx_sections
CHARACTER(LEN=*), PARAMETER :: routineN = 'hfx_matrix', routineP = moduleN//':'//routineN
INTEGER :: handle, irep, ispin, mspin, n_rep_hf, &
nspins
LOGICAL :: distribute_fock_matrix, &
hfx_treat_lsd_in_core, &
s_mstruct_changed
REAL(KIND=dp) :: eh1
TYPE(cp_para_env_type), POINTER :: para_env
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp, rho_ao_kp
TYPE(dft_control_type), POINTER :: dft_control
TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
CALL timeset(routineN, handle)
NULLIFY (dft_control, para_env, matrix_ks_kp, rho_ao_kp)
CALL get_qs_env(qs_env=qs_env, &
dft_control=dft_control, &
para_env=para_env, &
s_mstruct_changed=s_mstruct_changed, &
x_data=x_data)
CPASSERT(dft_control%nimages == 1)
nspins = dft_control%nspins
CALL section_vals_get(hfx_sections, n_repetition=n_rep_hf)
CALL section_vals_val_get(hfx_sections, "TREAT_LSD_IN_CORE", l_val=hfx_treat_lsd_in_core, &
i_rep_section=1)
CALL section_vals_get(hfx_sections, n_repetition=n_rep_hf)
distribute_fock_matrix = .TRUE.
mspin = 1
IF (hfx_treat_lsd_in_core) mspin = nspins
matrix_ks_kp(1:nspins, 1:1) => matrix_ks(1:nspins)
rho_ao_kp(1:nspins, 1:1) => rho_ao(1:nspins)
DO irep = 1, n_rep_hf
DO ispin = 1, mspin
CALL integrate_four_center(qs_env, x_data, matrix_ks_kp, eh1, rho_ao_kp, hfx_sections, para_env, &
s_mstruct_changed, irep, distribute_fock_matrix, ispin=ispin)
END DO
END DO
CALL timestop(handle)
END SUBROUTINE hfx_matrix
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param p_env ...
!> \param c0 ...
!> \param v ...
!> \param Av ...
! **************************************************************************************************
SUBROUTINE apply_xc_admm(qs_env, p_env, c0, v, Av)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(qs_p_env_type), POINTER :: p_env
TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: c0, v, Av
CHARACTER(len=*), PARAMETER :: routineN = 'apply_xc_admm', routineP = moduleN//':'//routineN
INTEGER :: handle, ispin, nao, nao_aux, ncol, nspins
INTEGER, DIMENSION(2, 3) :: bo
LOGICAL :: lsd
REAL(KIND=dp) :: fac
TYPE(admm_type), POINTER :: admm_env
TYPE(cp_fm_type), POINTER :: tc0, tc1
TYPE(dbcsr_p_type) :: xcmat
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
TYPE(dft_control_type), POINTER :: dft_control
TYPE(linres_control_type), POINTER :: linres_control
TYPE(pw_env_type), POINTER :: pw_env
TYPE(pw_p_type), DIMENSION(:), POINTER :: rho1_aux_g, rho1_aux_r, tau_pw, v_xc
TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
TYPE(section_vals_type), POINTER :: xc_fun_section, xc_section
TYPE(xc_rho_cflags_type) :: needs
TYPE(xc_rho_set_type), POINTER :: rho1_set
CALL timeset(routineN, handle)
CPASSERT(ASSOCIATED(c0))
CPASSERT(ASSOCIATED(v))
CPASSERT(ASSOCIATED(Av))
CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
IF (dft_control%do_admm) THEN
IF (dft_control%admm_control%aux_exch_func == do_admm_aux_exch_func_none) THEN
! nothing to do
ELSE
CALL get_qs_env(qs_env=qs_env, linres_control=linres_control)
CPASSERT(.NOT. dft_control%qs_control%gapw)
CPASSERT(.NOT. dft_control%qs_control%gapw_xc)
CPASSERT(.NOT. dft_control%qs_control%lrigpw)
CPASSERT(.NOT. linres_control%lr_triplet)
nspins = dft_control%nspins
! AUX basis contribution
CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
CPASSERT(ASSOCIATED(pw_env))
CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
ALLOCATE (v_xc(nspins))
DO ispin = 1, nspins
NULLIFY (v_xc(ispin)%pw)
CALL pw_pool_create_pw(auxbas_pw_pool, v_xc(ispin)%pw, &
use_data=REALDATA3D, in_space=REALSPACE)
CALL pw_zero(v_xc(ispin)%pw)
END DO
NULLIFY (tau_pw)
! calculate the xc potential
lsd = (nspins == 2)
CALL get_qs_env(qs_env=qs_env, matrix_s_aux_fit=matrix_s)
ALLOCATE (xcmat%matrix)
CALL dbcsr_create(xcmat%matrix, template=matrix_s(1)%matrix)
ALLOCATE (rho1_aux_r(nspins), rho1_aux_g(nspins))
DO ispin = 1, nspins
NULLIFY (rho1_aux_r(ispin)%pw, rho1_aux_g(ispin)%pw)
CALL pw_pool_create_pw(auxbas_pw_pool, rho1_aux_r(ispin)%pw, &
use_data=REALDATA3D, in_space=REALSPACE)
CALL pw_pool_create_pw(auxbas_pw_pool, rho1_aux_g(ispin)%pw, &
in_space=RECIPROCALSPACE, use_data=COMPLEXDATA1D)
END DO
CALL admm_aux_reponse_density(qs_env, c0, v, rho1_aux_r, rho1_aux_g)
CALL get_qs_env(qs_env, admm_env=admm_env)
xc_section => admm_env%xc_section_aux
bo = rho1_aux_r(1)%pw%pw_grid%bounds_local
! create the place where to store the argument for the functionals
NULLIFY (rho1_set)
CALL xc_rho_set_create(rho1_set, bo, &
rho_cutoff=section_get_rval(xc_section, "DENSITY_CUTOFF"), &
drho_cutoff=section_get_rval(xc_section, "GRADIENT_CUTOFF"), &
tau_cutoff=section_get_rval(xc_section, "TAU_CUTOFF"))
xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
needs = xc_functionals_get_needs(xc_fun_section, lsd, .TRUE.)
! calculate the arguments needed by the functionals
CALL xc_rho_set_update(rho1_set, rho1_aux_r, rho1_aux_g, tau_pw, needs, &
section_get_ival(xc_section, "XC_GRID%XC_DERIV"), &
section_get_ival(xc_section, "XC_GRID%XC_SMOOTH_RHO"), &
auxbas_pw_pool)
fac = 0._dp
CALL xc_calc_2nd_deriv(v_xc, p_env%kpp1_env%deriv_set_admm, &
p_env%kpp1_env%rho_set_admm, &
rho1_set, auxbas_pw_pool, xc_section=xc_section, &
tddfpt_fac=fac)
CALL xc_rho_set_release(rho1_set)
tc0 => admm_env%work_aux_orb
tc1 => admm_env%work_aux_orb2
CALL cp_fm_get_info(admm_env%A, nrow_global=nao_aux)
DO ispin = 1, nspins
v_xc(ispin)%pw%cr3d = v_xc(ispin)%pw%cr3d*v_xc(ispin)%pw%pw_grid%dvol
IF (nspins == 1) THEN
v_xc(ispin)%pw%cr3d = 2.0_dp*v_xc(ispin)%pw%cr3d
END IF
CALL dbcsr_copy(xcmat%matrix, matrix_s(1)%matrix)
CALL dbcsr_set(xcmat%matrix, 0.0_dp)
CALL integrate_v_rspace(v_rspace=v_xc(ispin), hmat=xcmat, qs_env=qs_env, &
calculate_forces=.FALSE., basis_type="AUX_FIT")
CALL cp_fm_get_info(c0(ispin)%matrix, nrow_global=nao, ncol_global=ncol)
CALL cp_gemm('N', 'N', nao_aux, ncol, nao, 1.0_dp, admm_env%A, &
c0(ispin)%matrix, 0.0_dp, tc0)
CALL cp_dbcsr_sm_fm_multiply(xcmat%matrix, tc0, tc1, &
ncol=ncol, alpha=1.0_dp, beta=0.0_dp)
CALL cp_gemm('T', 'N', nao, ncol, nao_aux, 1.0_dp, admm_env%A, &
tc1, 1.0_dp, Av(ispin)%matrix)
END DO
DO ispin = 1, nspins
CALL pw_pool_give_back_pw(auxbas_pw_pool, v_xc(ispin)%pw)
CALL pw_pool_give_back_pw(auxbas_pw_pool, rho1_aux_r(ispin)%pw)
CALL pw_pool_give_back_pw(auxbas_pw_pool, rho1_aux_g(ispin)%pw)
END DO
DEALLOCATE (v_xc, rho1_aux_r, rho1_aux_g)
CALL dbcsr_deallocate_matrix(xcmat%matrix)
END IF
END IF
CALL timestop(handle)
END SUBROUTINE apply_xc_admm
! **************************************************************************************************
!> \brief Calculate ADMM auxiliary response density
!> \param qs_env ...
!> \param c0 ...
!> \param c1 ...
!> \param rho1_aux_r ...
!> \param rho1_aux_g ...
! **************************************************************************************************
SUBROUTINE admm_aux_reponse_density(qs_env, c0, c1, rho1_aux_r, rho1_aux_g)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: c0, c1
TYPE(pw_p_type), DIMENSION(:), POINTER :: rho1_aux_r, rho1_aux_g
CHARACTER(LEN=*), PARAMETER :: routineN = 'admm_aux_reponse_density', &
routineP = moduleN//':'//routineN
INTEGER :: handle, ispin, nao, nao_aux, ncol, nspins
REAL(KIND=dp) :: tot_rho_aux
TYPE(admm_type), POINTER :: admm_env
TYPE(cp_fm_type), POINTER :: tc0, tc1
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, rho1_ao
TYPE(dft_control_type), POINTER :: dft_control
TYPE(qs_ks_env_type), POINTER :: ks_env
TYPE(task_list_type), POINTER :: task_list_aux_fit
CALL timeset(routineN, handle)
CALL get_qs_env(qs_env=qs_env, &
ks_env=ks_env, &
matrix_s=matrix_s, &
dft_control=dft_control)
nspins = dft_control%nspins
CALL get_qs_env(qs_env, admm_env=admm_env, &
task_list_aux_fit=task_list_aux_fit)
CPASSERT(ASSOCIATED(admm_env%A))
CPASSERT(ASSOCIATED(admm_env%work_aux_orb))
CPASSERT(ASSOCIATED(admm_env%work_aux_orb2))
CALL cp_fm_get_info(admm_env%A, nrow_global=nao_aux)
CALL get_qs_env(qs_env=qs_env, matrix_s_aux_fit=matrix_s)
NULLIFY (rho1_ao)
CALL dbcsr_allocate_matrix_set(rho1_ao, nspins)
DO ispin = 1, nspins
ALLOCATE (rho1_ao(ispin)%matrix)
CALL dbcsr_create(rho1_ao(ispin)%matrix, template=matrix_s(1)%matrix)
CALL dbcsr_copy(rho1_ao(ispin)%matrix, matrix_s(1)%matrix)
CALL dbcsr_set(rho1_ao(ispin)%matrix, 0.0_dp)
END DO
! P1 -> AUX BASIS
DO ispin = 1, nspins
CALL cp_fm_get_info(c0(ispin)%matrix, nrow_global=nao, ncol_global=ncol)
tc0 => admm_env%work_aux_orb
tc1 => admm_env%work_aux_orb2
CALL cp_gemm('N', 'N', nao_aux, ncol, nao, 1.0_dp, admm_env%A, &
c1(ispin)%matrix, 0.0_dp, tc1)
CALL cp_gemm('N', 'N', nao_aux, ncol, nao, 1.0_dp, admm_env%A, &
c0(ispin)%matrix, 0.0_dp, tc0)
CALL cp_dbcsr_plus_fm_fm_t(rho1_ao(ispin)%matrix, matrix_v=tc1, matrix_g=tc0, &
ncol=ncol, alpha=1.0_dp)
CALL cp_dbcsr_plus_fm_fm_t(rho1_ao(ispin)%matrix, matrix_v=tc0, matrix_g=tc1, &
ncol=ncol, alpha=1.0_dp)
ENDDO
DO ispin = 1, nspins
CALL calculate_rho_elec(matrix_p=rho1_ao(ispin)%matrix, &
rho=rho1_aux_r(ispin), rho_gspace=rho1_aux_g(ispin), &
total_rho=tot_rho_aux, ks_env=ks_env, &
basis_type="AUX_FIT", &
task_list_external=task_list_aux_fit)
END DO
CALL dbcsr_deallocate_matrix_set(rho1_ao)
CALL timestop(handle)
END SUBROUTINE admm_aux_reponse_density
! **************************************************************************************************
!> \brief ...
!> \param p_env ...
@ -1530,7 +1954,9 @@ CONTAINS
routineP = moduleN//':'//routineN
INTEGER :: ispin, nspins
TYPE(admm_type), POINTER :: admm_env
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
TYPE(dft_control_type), POINTER :: dft_control
TYPE(pw_env_type), POINTER :: pw_env
TYPE(pw_p_type), DIMENSION(:), POINTER :: my_rho_r, rho_r
TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
@ -1568,8 +1994,9 @@ CONTAINS
END DO
END IF
IF (.NOT. ASSOCIATED(kpp1_env%deriv_set)) THEN
CALL get_qs_env(qs_env, admm_env=admm_env, dft_control=dft_control)
IF (.NOT. ASSOCIATED(kpp1_env%deriv_set)) THEN
IF (nspins == 1 .AND. lr_triplet) THEN
ALLOCATE (my_rho_r(2))
DO ispin = 1, 2
@ -1585,7 +2012,11 @@ CONTAINS
END DO
END IF
xc_section => section_vals_get_subs_vals(input, "DFT%XC")
IF (dft_control%do_admm) THEN
xc_section => admm_env%xc_section_primary
ELSE
xc_section => section_vals_get_subs_vals(input, "DFT%XC")
END IF
CALL xc_prep_2nd_deriv(kpp1_env%deriv_set, kpp1_env%rho_set, &
my_rho_r, auxbas_pw_pool, &
@ -1597,6 +2028,21 @@ CONTAINS
DEALLOCATE (my_rho_r)
ENDIF
! ADMM Correction
IF (dft_control%do_admm) THEN
IF (dft_control%admm_control%aux_exch_func /= do_admm_aux_exch_func_none) THEN
IF (.NOT. ASSOCIATED(kpp1_env%deriv_set_admm)) THEN
CPASSERT(.NOT. lr_triplet)
xc_section => admm_env%xc_section_aux
CALL get_qs_env(qs_env=qs_env, rho_aux_fit=rho)
CALL qs_rho_get(rho, rho_r=rho_r)
CALL xc_prep_2nd_deriv(kpp1_env%deriv_set_admm, kpp1_env%rho_set_admm, &
rho_r, auxbas_pw_pool, &
xc_section=xc_section)
END IF
END IF
END IF
END SUBROUTINE kpp1_check_i_alloc
! **************************************************************************************************

21
src/qs_moments.F
View file

@ -107,8 +107,9 @@ CONTAINS
!> \param nmoments ...
!> \param ref_point ...
!> \param ref_points ...
!> \param basis_type ...
! **************************************************************************************************
SUBROUTINE build_local_moment_matrix(qs_env, moments, nmoments, ref_point, ref_points)
SUBROUTINE build_local_moment_matrix(qs_env, moments, nmoments, ref_point, ref_points, basis_type)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: moments
@ -116,6 +117,7 @@ CONTAINS
REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: ref_point
REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
OPTIONAL :: ref_points
CHARACTER(len=*), OPTIONAL :: basis_type
CHARACTER(LEN=*), PARAMETER :: routineN = 'build_local_moment_matrix', &
routineP = moduleN//':'//routineN
@ -164,7 +166,8 @@ CONTAINS
! Allocate work storage
CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
maxco=maxco, maxsgf=maxsgf)
maxco=maxco, maxsgf=maxsgf, &
basis_type=basis_type)
ALLOCATE (mab(maxco, maxco, nm))
mab(:, :, :) = 0.0_dp
@ -180,7 +183,7 @@ CONTAINS
ALLOCATE (basis_set_list(nkind))
DO ikind = 1, nkind
qs_kind => qs_kind_set(ikind)
CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a)
CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type=basis_type)
IF (ASSOCIATED(basis_set_a)) THEN
basis_set_list(ikind)%gto_basis_set => basis_set_a
ELSE
@ -332,8 +335,9 @@ CONTAINS
!> \param nmoments ...
!> \param ref_point ...
!> \param ref_points ...
!> \param basis_type ...
! **************************************************************************************************
SUBROUTINE build_local_magmom_matrix(qs_env, magmom, nmoments, ref_point, ref_points)
SUBROUTINE build_local_magmom_matrix(qs_env, magmom, nmoments, ref_point, ref_points, basis_type)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: magmom
@ -341,6 +345,7 @@ CONTAINS
REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: ref_point
REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
OPTIONAL :: ref_points
CHARACTER(len=*), OPTIONAL :: basis_type
CHARACTER(LEN=*), PARAMETER :: routineN = 'build_local_magmom_matrix', &
routineP = moduleN//':'//routineN
@ -416,7 +421,7 @@ CONTAINS
ALLOCATE (basis_set_list(nkind))
DO ikind = 1, nkind
qs_kind => qs_kind_set(ikind)
CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a)
CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type=basis_type)
IF (ASSOCIATED(basis_set_a)) THEN
basis_set_list(ikind)%gto_basis_set => basis_set_a
ELSE
@ -723,12 +728,14 @@ CONTAINS
!> \param cosmat ...
!> \param sinmat ...
!> \param kvec ...
!> \param basis_type ...
! **************************************************************************************************
SUBROUTINE build_berry_kpoint_matrix(qs_env, cosmat, sinmat, kvec)
SUBROUTINE build_berry_kpoint_matrix(qs_env, cosmat, sinmat, kvec, basis_type)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: cosmat, sinmat
REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: kvec
CHARACTER(len=*), OPTIONAL :: basis_type
CHARACTER(LEN=*), PARAMETER :: routineN = 'build_berry_kpoint_matrix', &
routineP = moduleN//':'//routineN
@ -785,7 +792,7 @@ CONTAINS
ALLOCATE (basis_set_list(nkind))
DO ikind = 1, nkind
qs_kind => qs_kind_set(ikind)
CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set)
CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type=basis_type)
IF (ASSOCIATED(basis_set)) THEN
basis_set_list(ikind)%gto_basis_set => basis_set
ELSE

1
src/qs_p_env_methods.F
View file

@ -152,6 +152,7 @@ CONTAINS
ALLOCATE (p_env)
NULLIFY (p_env%kpp1, &
p_env%p1, &
p_env%w1, &
p_env%m_epsilon, &
p_env%psi0d, &
p_env%S_psi0, &

5
src/qs_p_env_types.F
View file

@ -60,7 +60,8 @@ MODULE qs_p_env_types
LOGICAL :: orthogonal_orbitals
INTEGER :: id_nr, ref_count, iter
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: kpp1, p1
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: kpp1, p1, w1
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p1_admm => NULL()
TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: m_epsilon, &
psi0d, S_psi0, Smo_inv
TYPE(qs_kpp1_env_type), POINTER :: kpp1_env
@ -149,6 +150,8 @@ CONTAINS
CALL qs_rho_release(p_env%rho1)
IF (ASSOCIATED(p_env%kpp1)) CALL dbcsr_deallocate_matrix_set(p_env%kpp1)
IF (ASSOCIATED(p_env%p1)) CALL dbcsr_deallocate_matrix_set(p_env%p1)
IF (ASSOCIATED(p_env%w1)) CALL dbcsr_deallocate_matrix_set(p_env%w1)
IF (ASSOCIATED(p_env%p1_admm)) CALL dbcsr_deallocate_matrix_set(p_env%p1_admm)
IF (ASSOCIATED(p_env%local_rho_set)) THEN
CALL local_rho_set_release(p_env%local_rho_set)
END IF

1196
src/response_solver.F Normal file
View file

File diff suppressed because it is too large Load diff

9
src/xc_adiabatic_utils.F
View file

@ -84,6 +84,7 @@ CONTAINS
scale_dDFA, scale_ddW0, scale_dEx1, &
scale_dEx2, total_energy_xc
TYPE(cp_para_env_type), POINTER :: para_env
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao_resp
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao
TYPE(dft_control_type), POINTER :: dft_control
TYPE(qs_ks_env_type), POINTER :: ks_env
@ -93,7 +94,7 @@ CONTAINS
CALL timeset(routineN, handle)
NULLIFY (para_env, dft_control, adiabatic_rescaling_section, hfx_sections, &
input, xc_section, rho_xc, ks_env, rho_ao)
input, xc_section, rho_xc, ks_env, rho_ao, rho_ao_resp)
CALL get_qs_env(qs_env, &
dft_control=dft_control, &
@ -146,9 +147,11 @@ CONTAINS
IF (calculate_forces) THEN
CPASSERT(.NOT. use_virial)
!! we also have to scale the forces!!!!
CALL derivatives_four_center(qs_env, rho_ao, hfx_sections, para_env, 1, use_virial, &
CALL derivatives_four_center(qs_env, rho_ao, rho_ao_resp, hfx_sections, &
para_env, 1, use_virial, &
adiabatic_rescale_factor=scale_dEx1)
CALL derivatives_four_center(qs_env, rho_ao, hfx_sections, para_env, 2, use_virial, &
CALL derivatives_four_center(qs_env, rho_ao, rho_ao_resp, hfx_sections, &
para_env, 2, use_virial, &
adiabatic_rescale_factor=scale_dEx2)
END IF

2
tests/QS/regtest-debug/TEST_FILES → tests/QS/regtest-debug-1/TEST_FILES
View file

@ -9,6 +9,4 @@ h2o_pdip.inp 86 1e-05
h2o_periodic.inp 87 1e-05 0.139741440657E+02
h2o_gga.inp 87 1e-05 0.163373158129E+02
h2o_gapw.inp 87 1e-05 0.163123128994E+02
h2o_lri.inp 87 1e-05 0.216331245542E+02
h2o_pade_fd.inp 87 1e-05 0.164539007391E+02
#EOF

0
tests/QS/regtest-debug/TEST_FILES_RESET → tests/QS/regtest-debug-1/TEST_FILES_RESET
View file

0
tests/QS/regtest-debug/h2o_dipole.inp → tests/QS/regtest-debug-1/h2o_dipole.inp
View file

0
tests/QS/regtest-debug/h2o_gapw.inp → tests/QS/regtest-debug-1/h2o_gapw.inp
View file

0
tests/QS/regtest-debug/h2o_gga.inp → tests/QS/regtest-debug-1/h2o_gga.inp
View file

0
tests/QS/regtest-debug/h2o_pdip.inp → tests/QS/regtest-debug-1/h2o_pdip.inp
View file

0
tests/QS/regtest-debug/h2o_periodic.inp → tests/QS/regtest-debug-1/h2o_periodic.inp
View file

0
tests/QS/regtest-debug/h2o_polar.inp → tests/QS/regtest-debug-1/h2o_polar.inp
View file

10
tests/QS/regtest-debug-2/TEST_FILES Normal file
View file

@ -0,0 +1,10 @@
# runs are executed in the same order as in this file
# the second field tells which test should be run in order to compare with the last available output
# e.g. 0 means do not compare anything, running is enough
# 1 compares the last total energy in the file
# for details see cp2k/tools/do_regtest
h2o_lri.inp 87 1e-05 0.216331245542E+02
h2o_pade_fd.inp 87 1e-05 0.164539007391E+02
h2o_hfx.inp 87 1e-05 0.155040875037E+02
h2o_hfx_admm.inp 87 1e-05 0.151709407272E+02
#EOF

1
tests/QS/regtest-debug-2/TEST_FILES_RESET Normal file
View file

@ -0,0 +1 @@
#

88
tests/QS/regtest-debug-2/h2o_hfx.inp Normal file
View file

@ -0,0 +1,88 @@
&FORCE_EVAL
METHOD Quickstep
&PROPERTIES
&LINRES
PRECONDITIONER FULL_ALL
EPS 1.e-10
&POLAR
DO_RAMAN T
PERIODIC_DIPOLE_OPERATOR F
&END
&END
&END
&DFT
&QS
METHOD GPW
EPS_DEFAULT 1.e-10
&END QS
&EFIELD
&END
&SCF
SCF_GUESS ATOMIC
&OT
PRECONDITIONER FULL_SINGLE_INVERSE
MINIMIZER DIIS
&END
&OUTER_SCF
MAX_SCF 10
EPS_SCF 1.0E-7
&END
MAX_SCF 10
EPS_SCF 1.0E-7
&END SCF
&XC
&XC_FUNCTIONAL NONE
&END XC_FUNCTIONAL
&HF
&SCREENING
EPS_SCHWARZ 1.0E-10
&END
&END
&END XC
&PRINT
&MOMENTS ON
PERIODIC .FALSE.
REFERENCE COM
&END
&END
&END DFT
&SUBSYS
&CELL
ABC [angstrom] 6.0 6.0 6.0
PERIODIC NONE
&END
&COORD
O 0.000000 0.000000 -0.065587
H 0.000000 -0.757136 0.520545
H 0.000000 0.757136 0.520545
&END COORD
&TOPOLOGY
&CENTER_COORDINATES
&END
&END
&KIND H
BASIS_SET DZV-GTH-PADE
POTENTIAL GTH-PADE-q1
&END KIND
&KIND O
BASIS_SET DZVP-GTH-PADE
POTENTIAL GTH-PADE-q6
&END KIND
&END SUBSYS
&END FORCE_EVAL
&GLOBAL
PRINT_LEVEL LOW
PROJECT dipole
RUN_TYPE DEBUG
&END GLOBAL
&DEBUG
DEBUG_FORCES .FALSE.
DEBUG_STRESS_TENSOR .FALSE.
DEBUG_DIPOLE .FALSE.
DEBUG_POLARIZABILITY .TRUE.
DE 0.0002
EPS_NO_ERROR_CHECK 5.e-5
&END

100
tests/QS/regtest-debug-2/h2o_hfx_admm.inp Normal file
View file

@ -0,0 +1,100 @@
&FORCE_EVAL
METHOD Quickstep
&PROPERTIES
&LINRES
PRECONDITIONER FULL_ALL
EPS 1.e-10
&POLAR
DO_RAMAN T
PERIODIC_DIPOLE_OPERATOR F
&END
&END
&END
&DFT
BASIS_SET_FILE_NAME BASIS_SET
BASIS_SET_FILE_NAME BASIS_ADMM
&QS
METHOD GPW
EPS_DEFAULT 1.e-10
&END QS
&AUXILIARY_DENSITY_MATRIX_METHOD
ADMM_PURIFICATION_METHOD NONE
EXCH_CORRECTION_FUNC NONE
EXCH_SCALING_MODEL NONE
METHOD BASIS_PROJECTION
&END
&EFIELD
&END
&SCF
SCF_GUESS ATOMIC
&OT OFF
PRECONDITIONER FULL_SINGLE_INVERSE
MINIMIZER DIIS
&END
&OUTER_SCF
MAX_SCF 10
EPS_SCF 1.0E-7
&END
MAX_SCF 100
EPS_SCF 1.0E-7
&END SCF
&XC
&XC_FUNCTIONAL NONE
&END XC_FUNCTIONAL
&HF
&SCREENING
EPS_SCHWARZ 1.0E-10
&END
&END
&END XC
&PRINT
&MOMENTS ON
PERIODIC .FALSE.
REFERENCE COM
&END
&END
&END DFT
&SUBSYS
&CELL
ABC [angstrom] 6.0 6.0 6.0
PERIODIC NONE
&END
&COORD
O 0.000000 0.000000 -0.065587
H 0.000000 -0.757136 0.520545
H 0.000000 0.757136 0.520545
&END COORD
&TOPOLOGY
&CENTER_COORDINATES
&END
&END
&KIND H
BASIS_SET DZV-GTH-PADE
# BASIS_SET AUX_FIT DZV-GTH-PADE
BASIS_SET AUX_FIT FIT3
POTENTIAL GTH-PADE-q1
&END KIND
&KIND O
BASIS_SET DZVP-GTH-PADE
# BASIS_SET AUX_FIT DZVP-GTH-PADE
BASIS_SET AUX_FIT FIT3
POTENTIAL GTH-PADE-q6
&END KIND
&END SUBSYS
&END FORCE_EVAL
&GLOBAL
PRINT_LEVEL LOW
PROJECT dipole
RUN_TYPE DEBUG
&END GLOBAL
&DEBUG
DEBUG_FORCES .FALSE.
DEBUG_STRESS_TENSOR .FALSE.
DEBUG_DIPOLE .FALSE.
DEBUG_POLARIZABILITY .TRUE.
DE 0.0002
EPS_NO_ERROR_CHECK 5.e-5
&END

0
tests/QS/regtest-debug/h2o_lri.inp → tests/QS/regtest-debug-2/h2o_lri.inp
View file

0
tests/QS/regtest-debug/h2o_pade_fd.inp → tests/QS/regtest-debug-2/h2o_pade_fd.inp
View file

8
tests/QS/regtest-debug-3/TEST_FILES Normal file
View file

@ -0,0 +1,8 @@
# runs are executed in the same order as in this file
# the second field tells which test should be run in order to compare with the last available output
# e.g. 0 means do not compare anything, running is enough
# 1 compares the last total energy in the file
# for details see cp2k/tools/do_regtest
h2o_pbe0.inp 87 1e-05 0.157066941183E+02
h2o_pbe0_admm.inp 87 1e-05 0.159060491447E+02
#EOF

1
tests/QS/regtest-debug-3/TEST_FILES_RESET Normal file
View file

@ -0,0 +1 @@
#

87
tests/QS/regtest-debug-3/h2o_pbe0.inp Normal file
View file

@ -0,0 +1,87 @@
&FORCE_EVAL
METHOD Quickstep
&PROPERTIES
&LINRES
PRECONDITIONER FULL_ALL
EPS 1.e-10
&POLAR
DO_RAMAN T
PERIODIC_DIPOLE_OPERATOR F
&END
&END
&END
&DFT
&QS
METHOD GPW
EPS_DEFAULT 1.e-10
&END QS
&EFIELD
&END
&SCF
SCF_GUESS ATOMIC
&OT
PRECONDITIONER FULL_SINGLE_INVERSE
MINIMIZER DIIS
&END
&OUTER_SCF
MAX_SCF 10
EPS_SCF 1.0E-7
&END
MAX_SCF 10
EPS_SCF 1.0E-7
&END SCF
&XC
&XC_FUNCTIONAL PBE0
&END XC_FUNCTIONAL
&END XC
&PRINT
&MOMENTS ON
PERIODIC .FALSE.
REFERENCE COM
&END
&END
&POISSON
POISSON_SOLVER MT
PERIODIC NONE
&END
&END DFT
&SUBSYS
&CELL
ABC [angstrom] 6.0 6.0 6.0
PERIODIC NONE
&END
&COORD
O 0.000000 0.000000 -0.065587
H 0.000000 -0.757136 0.520545
H 0.000000 0.757136 0.520545
&END COORD
&TOPOLOGY
&CENTER_COORDINATES
&END
&END
&KIND H
BASIS_SET DZV-GTH-PADE
POTENTIAL GTH-PADE-q1
&END KIND
&KIND O
BASIS_SET DZVP-GTH-PADE
POTENTIAL GTH-PADE-q6
&END KIND
&END SUBSYS
&END FORCE_EVAL
&GLOBAL
PRINT_LEVEL LOW
PROJECT dipole
RUN_TYPE DEBUG
&END GLOBAL
&DEBUG
DEBUG_FORCES .FALSE.
DEBUG_STRESS_TENSOR .FALSE.
DEBUG_DIPOLE .FALSE.
DEBUG_POLARIZABILITY .TRUE.
DE 0.002
EPS_NO_ERROR_CHECK 5.e-5
&END

101
tests/QS/regtest-debug-3/h2o_pbe0_admm.inp Normal file
View file

@ -0,0 +1,101 @@
&FORCE_EVAL
METHOD Quickstep
&PROPERTIES
&LINRES
PRECONDITIONER FULL_ALL
EPS 1.e-10
&POLAR
DO_RAMAN T
PERIODIC_DIPOLE_OPERATOR F
&END
&END
&END
&DFT
BASIS_SET_FILE_NAME BASIS_SET
BASIS_SET_FILE_NAME BASIS_ADMM
&QS
METHOD GPW
EPS_DEFAULT 1.e-10
&END QS
&AUXILIARY_DENSITY_MATRIX_METHOD
ADMM_PURIFICATION_METHOD NONE
EXCH_CORRECTION_FUNC BECKE88X
EXCH_SCALING_MODEL NONE
METHOD BASIS_PROJECTION
&END
&EFIELD
&END
&SCF
SCF_GUESS RESTART
&OT OFF
PRECONDITIONER FULL_SINGLE_INVERSE
MINIMIZER DIIS
&END
&OUTER_SCF
MAX_SCF 10
EPS_SCF 1.0E-7
&END
MAX_SCF 100
EPS_SCF 1.0E-7
&END SCF
&XC
&XC_FUNCTIONAL PBE0
&END XC_FUNCTIONAL
&HF
&SCREENING
EPS_SCHWARZ 1.0E-10
&END
&END
&END XC
&PRINT
&MOMENTS ON
PERIODIC .FALSE.
REFERENCE COM
&END
&END
&END DFT
&SUBSYS
&CELL
ABC [angstrom] 6.0 6.0 6.0
PERIODIC NONE
&END
&COORD
O 0.000000 0.000000 -0.065587
H 0.000000 -0.757136 0.520545
H 0.000000 0.757136 0.520545
&END COORD
&TOPOLOGY
&CENTER_COORDINATES
&END
&END
&KIND H
BASIS_SET DZV-GTH-PADE
# BASIS_SET AUX_FIT DZV-GTH-PADE
BASIS_SET AUX_FIT FIT3
POTENTIAL GTH-PADE-q1
&END KIND
&KIND O
BASIS_SET DZVP-GTH-PADE
# BASIS_SET AUX_FIT DZVP-GTH-PADE
BASIS_SET AUX_FIT FIT3
POTENTIAL GTH-PADE-q6
&END KIND
&END SUBSYS
&END FORCE_EVAL
&GLOBAL
PRINT_LEVEL LOW
PROJECT dipole
RUN_TYPE DEBUG
# RUN_TYPE ENERGY
&END GLOBAL
&DEBUG
DEBUG_FORCES .FALSE.
DEBUG_STRESS_TENSOR .FALSE.
DEBUG_DIPOLE .FALSE.
DEBUG_POLARIZABILITY .TRUE.
DE 0.0002
EPS_NO_ERROR_CHECK 5.e-5
&END

20
tests/QS/regtest-kg/2H2O_ecmao.inp → tests/QS/regtest-ec/2H2O_ecmao.inp
View file

@ -26,17 +26,17 @@
&END
&END
&END XC
&ENERGY_CORRECTION
HARRIS_BASIS HARRIS
MAO
&XC
&XC_FUNCTIONAL
&PBE
&END
&END
&END XC
&END ENERGY_CORRECTION
&END KG_METHOD
&ENERGY_CORRECTION
HARRIS_BASIS HARRIS
MAO
&XC
&XC_FUNCTIONAL
&PBE
&END
&END
&END XC
&END ENERGY_CORRECTION
&SCF
EPS_SCF 1.0E-5
SCF_GUESS ATOMIC

24
tests/QS/regtest-kg/2H2O_ecmao2.inp → tests/QS/regtest-ec/2H2O_ecmao2.inp
View file

@ -26,19 +26,19 @@
&END
&END
&END XC
&ENERGY_CORRECTION
HARRIS_BASIS HARRIS
MAO
MAO_MAX_ITER 2000
MAO_EPS_GRAD 1.e-7
&XC
&XC_FUNCTIONAL
&PBE
&END
&END
&END XC
&END ENERGY_CORRECTION
&END KG_METHOD
&ENERGY_CORRECTION
HARRIS_BASIS HARRIS
MAO
MAO_MAX_ITER 2000
MAO_EPS_GRAD 1.e-7
&XC
&XC_FUNCTIONAL
&PBE
&END
&END
&END XC
&END ENERGY_CORRECTION
&SCF
EPS_SCF 1.0E-5
SCF_GUESS ATOMIC

61
tests/QS/regtest-ec/H2O_ec.inp Normal file
View file

@ -0,0 +1,61 @@
&FORCE_EVAL
METHOD Quickstep
&DFT
BASIS_SET_FILE_NAME BASIS_SET
POTENTIAL_FILE_NAME GTH_POTENTIALS
&MGRID
CUTOFF 200
&END MGRID
&QS
EPS_DEFAULT 1.E-10
EPS_KG_ORB 1.0E-5
&END QS
&ENERGY_CORRECTION
ENERGY_FUNCTIONAL HARRIS
HARRIS_BASIS ORBITAL
&XC
&XC_FUNCTIONAL
&PBE
&END
&END
&END XC
&END ENERGY_CORRECTION
&SCF
EPS_SCF 1.0E-5
SCF_GUESS ATOMIC
&END
&XC
&XC_FUNCTIONAL
&PADE
&END
&END
&END XC
&END DFT
&SUBSYS
&CELL
ABC 6.0 6.0 6.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
&END COORD
&KIND H
BASIS_SET ORB DZVP-GTH-BLYP
BASIS_SET HARRIS TZVDD3DF3PD-GTH-BLYP
POTENTIAL GTH-PADE-q1
MAO 1
&END KIND
&KIND O
BASIS_SET ORB DZVP-GTH-BLYP
BASIS_SET HARRIS TZVDD3DF3PD-GTH-BLYP
POTENTIAL GTH-PADE-q6
MAO 4
&END KIND
&END SUBSYS
&END FORCE_EVAL
&GLOBAL
PROJECT H2O
RUN_TYPE ENERGY_FORCE
PRINT_LEVEL LOW
&END GLOBAL

16
tests/QS/regtest-kg/H2_H2O-xcLLP_ec.inp → tests/QS/regtest-ec/H2_H2O-xcLLP_ec.inp
View file

@ -34,15 +34,15 @@
XC_DERIV SPLINE2
&END
&END XC
&ENERGY_CORRECTION
&XC
&XC_FUNCTIONAL
&PADE
&END
&END
&END XC
&END ENERGY_CORRECTION
&END KG_METHOD
&ENERGY_CORRECTION
&XC
&XC_FUNCTIONAL
&PADE
&END
&END
&END XC
&END ENERGY_CORRECTION
&SCF
EPS_SCF 1.0E-8
SCF_GUESS ATOMIC

18
tests/QS/regtest-kg/H2_H2O_ec.inp → tests/QS/regtest-ec/H2_H2O_ec.inp
View file

@ -26,16 +26,16 @@
&END
&END
&END XC
&ENERGY_CORRECTION
HARRIS_BASIS HARRIS
&XC
&XC_FUNCTIONAL
&PBE
&END
&END
&END XC
&END ENERGY_CORRECTION
&END KG_METHOD
&ENERGY_CORRECTION
HARRIS_BASIS HARRIS
&XC
&XC_FUNCTIONAL
&PBE
&END
&END
&END XC
&END ENERGY_CORRECTION
&SCF
EPS_SCF 1.0E-5
SCF_GUESS ATOMIC

18
tests/QS/regtest-kg/H2_H2O_ecprim.inp → tests/QS/regtest-ec/H2_H2O_ecprim.inp
View file

@ -26,16 +26,16 @@
&END
&END
&END XC
&ENERGY_CORRECTION
HARRIS_BASIS PRIMITIVE
&XC
&XC_FUNCTIONAL
&PBE
&END
&END
&END XC
&END ENERGY_CORRECTION
&END KG_METHOD
&ENERGY_CORRECTION
HARRIS_BASIS PRIMITIVE
&XC
&XC_FUNCTIONAL
&PBE
&END
&END
&END XC
&END ENERGY_CORRECTION
&SCF
EPS_SCF 1.0E-5
SCF_GUESS ATOMIC

61
tests/QS/regtest-ec/HF_ec_dipole.inp Normal file
View file

@ -0,0 +1,61 @@
&FORCE_EVAL
METHOD Quickstep
&DFT
&PRINT
&MOMENTS ON
PERIODIC .FALSE.
REFERENCE COM
&END
&END
BASIS_SET_FILE_NAME BASIS_SET
POTENTIAL_FILE_NAME GTH_POTENTIALS
&MGRID
CUTOFF 400
&END MGRID
&QS
EPS_DEFAULT 1.E-14
&END QS
&ENERGY_CORRECTION
ENERGY_FUNCTIONAL HARRIS
HARRIS_BASIS ORBITAL
&XC
&XC_FUNCTIONAL
&PBE
&END
&END
&END XC
&END ENERGY_CORRECTION
&SCF
EPS_SCF 1.0E-7
SCF_GUESS ATOMIC
&END
&XC
&XC_FUNCTIONAL
&PADE
&END
&END
&END XC
&END DFT
&SUBSYS
&CELL
ABC 6.0 6.0 6.0
&END CELL
&COORD
H 0.000000 0.000000 0.000000
F 0.000000 0.000000 1.050000
&END COORD
&KIND H
BASIS_SET ORB DZVP-GTH-BLYP
POTENTIAL GTH-PADE-q1
&END KIND
&KIND F
BASIS_SET ORB DZVP-GTH-BLYP
POTENTIAL GTH-PADE-q7
&END KIND
&END SUBSYS
&END FORCE_EVAL
&GLOBAL
PROJECT HF
RUN_TYPE ENERGY_FORCE
PRINT_LEVEL LOW
&END GLOBAL

65
tests/QS/regtest-ec/HF_ec_field.inp Normal file
View file

@ -0,0 +1,65 @@
&FORCE_EVAL
METHOD Quickstep
&DFT
&EFIELD
INTENSITY 0.001
POLARISATION 0.0 0.0 1.0
&END
&PRINT
&MOMENTS ON
PERIODIC .FALSE.
REFERENCE COM
&END
&END
BASIS_SET_FILE_NAME BASIS_SET
POTENTIAL_FILE_NAME GTH_POTENTIALS
&MGRID
CUTOFF 400
&END MGRID
&QS
EPS_DEFAULT 1.E-14
&END QS
&ENERGY_CORRECTION
ENERGY_FUNCTIONAL HARRIS
HARRIS_BASIS ORBITAL
&XC
&XC_FUNCTIONAL
&PBE
&END
&END
&END XC
&END ENERGY_CORRECTION
&SCF
EPS_SCF 1.0E-7
SCF_GUESS ATOMIC
&END
&XC
&XC_FUNCTIONAL
&PADE
&END
&END
&END XC
&END DFT
&SUBSYS
&CELL
ABC 6.0 6.0 6.0
&END CELL
&COORD
H 0.000000 0.000000 0.000000
F 0.000000 0.000000 1.050000
&END COORD
&KIND H
BASIS_SET ORB DZVP-GTH-BLYP
POTENTIAL GTH-PADE-q1
&END KIND
&KIND F
BASIS_SET ORB DZVP-GTH-BLYP
POTENTIAL GTH-PADE-q7
&END KIND
&END SUBSYS
&END FORCE_EVAL
&GLOBAL
PROJECT HF
RUN_TYPE ENERGY_FORCE
PRINT_LEVEL LOW
&END GLOBAL

75
tests/QS/regtest-ec/N2_ec-hfx-admm.inp Normal file
View file

@ -0,0 +1,75 @@
&FORCE_EVAL
METHOD Quickstep
&DFT
&PRINT
&DERIVATIVES
&END
&END
BASIS_SET_FILE_NAME BASIS_SET
BASIS_SET_FILE_NAME BASIS_ADMM
POTENTIAL_FILE_NAME GTH_POTENTIALS
&MGRID
CUTOFF 400
&END MGRID
&QS
EPS_DEFAULT 1.E-14
&END QS
&AUXILIARY_DENSITY_MATRIX_METHOD
ADMM_PURIFICATION_METHOD NONE
EXCH_CORRECTION_FUNC NONE
EXCH_SCALING_MODEL NONE
METHOD BASIS_PROJECTION
&END
&ENERGY_CORRECTION
ENERGY_FUNCTIONAL HARRIS
HARRIS_BASIS ORBITAL
&XC
&XC_FUNCTIONAL
&PADE
&END
&END
&END XC
&END ENERGY_CORRECTION
&SCF
EPS_SCF 1.0E-7
SCF_GUESS RESTART
&END
&XC
&XC_FUNCTIONAL NONE
&END
&HF
&SCREENING
EPS_SCHWARZ 1.0E-14
&END
&END
&END XC
&END DFT
&SUBSYS
&CELL
ABC 6.0 6.0 6.0
&END CELL
&COORD
N 0.000000 0.000000 0.650000
N 0.000000 0.000000 -0.650000
&END COORD
&KIND N
BASIS_SET ORB DZVP-GTH-BLYP
BASIS_SET HARRIS DZVP-GTH-BLYP
### BASIS_SET AUX_FIT FIT3
BASIS_SET AUX_FIT DZVP-GTH-BLYP
POTENTIAL GTH-PADE-q5
&END KIND
&END SUBSYS
&END FORCE_EVAL
&GLOBAL
PROJECT N2
RUN_TYPE DEBUG
## RUN_TYPE ENERGY_FORCE
PRINT_LEVEL LOW
&END GLOBAL
&DEBUG
DEBUG_FORCES T
DEBUG_STRESS_TENSOR F
STOP_ON_MISMATCH F
DX 0.001
&END

66
tests/QS/regtest-ec/N2_ec-hfx.inp Normal file
View file

@ -0,0 +1,66 @@
&FORCE_EVAL
METHOD Quickstep
&DFT
&PRINT
&DERIVATIVES
&END
&END
BASIS_SET_FILE_NAME BASIS_SET
POTENTIAL_FILE_NAME GTH_POTENTIALS
&MGRID
CUTOFF 400
&END MGRID
&QS
EPS_DEFAULT 1.E-14
&END QS
&ENERGY_CORRECTION
ENERGY_FUNCTIONAL HARRIS
HARRIS_BASIS ORBITAL
&XC
&XC_FUNCTIONAL
&PADE
&END
&END
&END XC
&END ENERGY_CORRECTION
&SCF
EPS_SCF 1.0E-8
SCF_GUESS ATOMIC
&END
&XC
&XC_FUNCTIONAL NONE
&END
&HF
&SCREENING
EPS_SCHWARZ 1.0E-14
&END
&END
&END XC
&END DFT
&SUBSYS
&CELL
ABC 6.0 6.0 6.0
&END CELL
&COORD
N 0.000000 0.000000 0.650000
N 0.000000 0.000000 -0.650000
&END COORD
&KIND N
BASIS_SET ORB DZVP-GTH-BLYP
BASIS_SET HARRIS DZVP-GTH-BLYP
POTENTIAL GTH-PADE-q5
&END KIND
&END SUBSYS
&END FORCE_EVAL
&GLOBAL
PROJECT N2
## RUN_TYPE DEBUG
RUN_TYPE ENERGY_FORCE
PRINT_LEVEL LOW
&END GLOBAL
&DEBUG
DEBUG_FORCES T
DEBUG_STRESS_TENSOR F
STOP_ON_MISMATCH F
DX 0.001
&END

63
tests/QS/regtest-ec/N2_ec.inp Normal file
View file

@ -0,0 +1,63 @@
&FORCE_EVAL
METHOD Quickstep
&DFT
&PRINT
&DERIVATIVES
&END
&END
BASIS_SET_FILE_NAME BASIS_SET
POTENTIAL_FILE_NAME GTH_POTENTIALS
&MGRID
CUTOFF 400
&END MGRID
&QS
EPS_DEFAULT 1.E-14
&END QS
&ENERGY_CORRECTION
ENERGY_FUNCTIONAL HARRIS
HARRIS_BASIS ORBITAL
&XC
&XC_FUNCTIONAL
&PBE
&END
&END
&END XC
&END ENERGY_CORRECTION
&SCF
EPS_SCF 1.0E-8
SCF_GUESS ATOMIC
&END
&XC
&XC_FUNCTIONAL
&PADE
&END
&END
&END XC
&END DFT
&SUBSYS
&CELL
ABC 6.0 6.0 6.0
&END CELL
&COORD
N 0.000000 0.000000 0.650000
N 0.000000 0.000000 -0.650000
&END COORD
&KIND N
BASIS_SET ORB DZVP-GTH-BLYP
BASIS_SET HARRIS DZVP-GTH-BLYP
POTENTIAL GTH-PADE-q5
&END KIND
&END SUBSYS
&END FORCE_EVAL
&GLOBAL
PROJECT N2
## RUN_TYPE DEBUG
RUN_TYPE ENERGY_FORCE
PRINT_LEVEL LOW
&END GLOBAL
&DEBUG
DEBUG_FORCES T
DEBUG_STRESS_TENSOR F
STOP_ON_MISMATCH F
DX 0.001
&END

12
tests/QS/regtest-ec/TEST_FILES Normal file
View file

@ -0,0 +1,12 @@
# 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
H2_H2O-xcLLP_ec.inp 11 1e-08 -18.1323810074
H2_H2O_ec.inp 11 1e-10 -18.4641086282
H2_H2O_ecprim.inp 11 1e-10 -18.4071025433
2H2O_ecmao.inp 11 1e-10 -34.0838901711
2H2O_ecmao2.inp 11 1e-08 -34.5024586832
H2O_ec.inp 11 1e-08 -17.2629537942
HF_ec_dipole.inp 11 1e-08 -24.8917060368
HF_ec_field.inp 11 1e-08 -24.8908706031
#EOF

83
tests/QS/regtest-ec/TEST_FILES_RESET Normal file
View file

@ -0,0 +1,83 @@
#
# add files to be reset here
#
H2_MD.inp
# redefine libxc interface
H2-libxc.inp
# faster testcase
H2-libxc.inp
# reset numerics
H2_MD-2.inp
H2_MD.inp
#
H2_MD.inp
#
H2_MD-2.inp
#
H2_MD-3.inp
#
H2_MD.inp
#
H2_MD-2.inp
#
H2_MD-3.inp
# new coloring algo
H2_MD-2.inp
# fix extrapolation for ls scf
H2_MD.inp
H2_MD-2.inp
H2_MD-3.inp
# fix extrapolation with S_PRECONDITIONER
H2_MD.inp
H2_MD-2.inp
H2_MD-3.inp
#
H2-libxc.inp
#
H2_H2O-xcLC.inp
#
H2_H2O-xcLLP.inp
#
H2_H2O-xcPW86.inp
#
H2_H2O-xcPW91.inp
#
H2_H2O-xcT92.inp
#
H2-libxc.inp
#
H2_H2O-xcLC.inp
#
H2_H2O-xcLLP.inp
#
H2_H2O-xcPW86.inp
#
H2_H2O-xcPW91.inp
#
H2_H2O-xcT92.inp
#
H2_MD.inp
#
H2_MD-2.inp
#
H2_MD-3.inp
#
H2-libxc.inp
#
H2-libxc-ot.inp
#
H2-libxc-diag.inp
#
H2_KG-1.inp
#
H2_KG-2.inp
#
H2_H2O-xcLC.inp
#
H2_H2O-xcLLP.inp
#
H2_H2O-xcPW86.inp
#
H2_H2O-xcPW91.inp
#
H2_H2O-xcT92.inp

5
tests/QS/regtest-kg/TEST_FILES
View file

@ -16,13 +16,8 @@ H2_H2O-xcLLP.inp 11 1e-10 -
H2_H2O-xcPW86.inp 11 1e-10 -18.143043340247019
H2_H2O-xcPW91.inp 11 1e-10 -18.142668216945395
H2_H2O-xcT92.inp 11 1e-10 -18.143724634599963
H2_H2O-xcLLP_ec.inp 66 1e-08 -18.1323810074
H2_H2O_ks.inp 11 1e-10 -18.130229622911152
H2_H2O_lsks.inp 11 1e-10 -18.130225854904928
H2_H2O_ec.inp 66 1e-10 -18.4641086282
H2_H2O_ecprim.inp 66 1e-10 -18.4071025433
2H2O_ecmao.inp 66 1e-10 -34.0838901711
2H2O_ecmao2.inp 66 1e-08 -34.5024586832
H2-none.inp 11 7e-07 -3.359680469888914
H2_H2O-vdW.inp 11 1e-10 -18.149003390914348
H2_H2O-lri.inp 72 3e-02 0.00011458

5
tests/TEST_DIRS
View file

@ -17,7 +17,9 @@ QS/regtest-cdft-hirshfeld
SIRIUS/regtest-1 sirius
QS/regtest-embed libint
QS/regtest-pod
QS/regtest-debug
QS/regtest-debug-1
QS/regtest-debug-2 libint
QS/regtest-debug-3 libint
QS/regtest-cdft-diag
xTB/regtest-1
xTB/regtest-2
@ -37,6 +39,7 @@ QS/regtest-p-efield
QS/regtest-mp2-stress libint
QS/regtest-ri-rpa libint
QS/regtest-kg libxc
QS/regtest-ec
QS/regtest-gpw-4
QS/regtest-gpw-8
QS/regtest-gpw-2-3