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

Add analytical derivatives of MO coefficients wrt nuclear coordinates

Browse source
This commit is contained in:
Edward Ditler 2021-11-17 15:21:58 +01:00 • committed by GitHub
parent 33488dafac
commit 67208cef2f
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
33 changed files with 5023 additions and 115 deletions
Download patch file Download diff file Expand all files Collapse all files

105
src/aobasis/ai_moments.F
View file

@ -46,9 +46,114 @@ MODULE ai_moments
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_moments'
PUBLIC :: cossin, moment, diffop, diff_momop, contract_cossin, dipole_force
PUBLIC :: diff_momop2
CONTAINS
! *****************************************************************************
!> \brief This returns the derivative of the moment integrals [a|\mu|b], with respect
!> to the position of the primitive on the left and right, i.e.
!> [da/dR_ai|\mu|b] + [a|\mu|d/dR_bi]
!> [da/dR_ai|\mu|b] = 2*zeta*[a+1i|\mu|b] - Ni(a)[a-1i|\mu|b]
!> [a|\mu|d/dR_bi] = 2*zetb*[a|\mu|b+1i] - Ni(b)[a|\mu|b-1i]
!> order indicates the max order of the moment operator to be calculated
!> 1: dipole
!> 2: quadrupole
!> ...
!> \param la_max ...
!> \param npgfa ...
!> \param zeta ...
!> \param rpgfa ...
!> \param la_min ...
!> \param lb_max ...
!> \param npgfb ...
!> \param zetb ...
!> \param rpgfb ...
!> \param lb_min ...
!> \param order ...
!> \param rac ...
!> \param rbc ...
!> \param difmab ...
!> \param mab_ext ...
!> \param deltaR needed for weighted derivative
!> \param iatom ...
!> \param jatom ...
!> SL August 2015, ED 2021
! **************************************************************************************************
SUBROUTINE diff_momop2(la_max, npgfa, zeta, rpgfa, la_min, &
lb_max, npgfb, zetb, rpgfb, lb_min, &
order, rac, rbc, difmab, mab_ext, deltaR, iatom, jatom)
INTEGER, INTENT(IN) :: la_max, npgfa
REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
INTEGER, INTENT(IN) :: la_min, lb_max, npgfb
REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
INTEGER, INTENT(IN) :: lb_min, order
REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rac, rbc
REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(OUT) :: difmab
REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, &
POINTER :: mab_ext
REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
OPTIONAL, POINTER :: deltaR
INTEGER, INTENT(IN), OPTIONAL :: iatom, jatom
INTEGER :: imom, istat, lda, lda_min, ldb, ldb_min
REAL(KIND=dp) :: dab, rab(3)
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: difmab_tmp
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: mab
rab = rbc - rac
dab = SQRT(SUM(rab**2))
lda_min = MAX(0, la_min - 1)
ldb_min = MAX(0, lb_min - 1)
lda = ncoset(la_max)*npgfa
ldb = ncoset(lb_max)*npgfb
ALLOCATE (difmab_tmp(lda, ldb, 3))
IF (PRESENT(mab_ext)) THEN
mab => mab_ext
ELSE
ALLOCATE (mab(npgfa*ncoset(la_max + 1), npgfb*ncoset(lb_max + 1), &
ncoset(order) - 1), STAT=istat)
mab = 0.0_dp
! *** Calculate the primitive moment integrals ***
CALL moment(la_max + 1, npgfa, zeta, rpgfa, lda_min, &
lb_max + 1, npgfb, zetb, rpgfb, &
order, rac, rbc, mab)
END IF
DO imom = 1, ncoset(order) - 1
difmab(:, :, imom, :) = 0.0_dp
difmab_tmp = 0.0_dp
CALL adbdr(la_max, npgfa, rpgfa, la_min, &
lb_max, npgfb, zetb, rpgfb, lb_min, &
dab, mab(:, :, imom), difmab_tmp(:, :, 1), &
difmab_tmp(:, :, 2), difmab_tmp(:, :, 3))
difmab(:, :, imom, 1) = difmab_tmp(:, :, 1)*deltaR(1, jatom)
difmab(:, :, imom, 2) = difmab_tmp(:, :, 2)*deltaR(2, jatom)
difmab(:, :, imom, 3) = difmab_tmp(:, :, 3)*deltaR(3, jatom)
difmab_tmp = 0.0_dp
CALL dabdr(la_max, npgfa, zeta, rpgfa, la_min, &
lb_max, npgfb, rpgfb, lb_min, &
dab, mab(:, :, imom), difmab_tmp(:, :, 1), &
difmab_tmp(:, :, 2), difmab_tmp(:, :, 3))
difmab(:, :, imom, 1) = difmab(:, :, imom, 1) + difmab_tmp(:, :, 1)*deltaR(1, iatom)
difmab(:, :, imom, 2) = difmab(:, :, imom, 2) + difmab_tmp(:, :, 2)*deltaR(2, iatom)
difmab(:, :, imom, 3) = difmab(:, :, imom, 3) + difmab_tmp(:, :, 3)*deltaR(3, iatom)
END DO
IF (PRESENT(mab_ext)) THEN
NULLIFY (mab)
ELSE
DEALLOCATE (mab)
END IF
DEALLOCATE (difmab_tmp)
END SUBROUTINE diff_momop2
! **************************************************************************************************
!> \brief ...
!> \param cos_block ...

60
src/aobasis/ai_oneelectron.F
View file

@ -78,13 +78,21 @@ CONTAINS
!> \param force_a ...
!> \param force_b ...
!> \param fs ...
!> \param vab2 The derivative of the 3-center integrals according to the weighting factors.
!> \param vab2_work ...
!> \param deltaR DIMENSION(3, natoms), weighting factors of the derivatives for each atom and direction
!> \param iatom ...
!> \param jatom ...
!> \param katom ...
!> \date May 2011
!> \author Juerg Hutter
!> \version 1.0
!> \note Extended by the derivatives for DFPT [Sandra Luber, Edward Ditler, 2021]
! **************************************************************************************************
SUBROUTINE os_3center(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
lb_max_set, lb_min_set, npgfb, rpgfb, zetb, auxint, rpgfc, &
rab, dab, rac, dac, rbc, dbc, vab, s, pab, force_a, force_b, fs)
rab, dab, rac, dac, rbc, dbc, vab, s, pab, force_a, force_b, fs, &
vab2, vab2_work, deltaR, iatom, jatom, katom)
INTEGER, INTENT(IN) :: la_max_set, la_min_set, npgfa
REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
INTEGER, INTENT(IN) :: lb_max_set, lb_min_set, npgfb
@ -103,14 +111,17 @@ CONTAINS
OPTIONAL :: pab
REAL(KIND=dp), DIMENSION(3), INTENT(OUT), OPTIONAL :: force_a, force_b
REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT), &
OPTIONAL :: fs
OPTIONAL :: fs, vab2, vab2_work
REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
OPTIONAL :: deltaR
INTEGER, INTENT(IN), OPTIONAL :: iatom, jatom, katom
INTEGER :: ax, ay, az, bx, by, bz, cda, cdax, cday, cdaz, cdb, cdbx, cdby, cdbz, coa, coamx, &
coamy, coamz, coapx, coapy, coapz, cob, cobmx, cobmy, cobmz, cobpx, cobpy, cobpz, da, &
da_max, dax, day, daz, db, db_max, dbx, dby, dbz, i, ia, iap, iax, iay, iaz, ib, ibm, &
ibx, iby, ibz, ii(3), iim(3), ij, ipgf, ir, ir1, ir2, irm(3), irr(3), irx, iry, irz, ix, &
ixx(1), j, jj(3), jjp(3), jpgf, la, la_max, la_min, lb, lb_max, lb_min, llr, m, ma, mb, &
mmax, na, nb
ibx, iby, ibz, idir, ii(3), iim(3), ij, ipgf, ir, ir1, ir2, irm(3), irr(3), irx, iry, &
irz, ix, ixx(1), j, jj(3), jjp(3), jpgf, la, la_max, la_min, lb, lb_max, lb_min, llr, m, &
ma, mb, mmax, na, nb
INTEGER, ALLOCATABLE, DIMENSION(:, :) :: iiap
LOGICAL :: calculate_force_a, calculate_force_b
REAL(KIND=dp) :: aai, abx, fax, fay, faz, fbx, fby, fbz, &
@ -148,6 +159,11 @@ CONTAINS
db_max = 0
END IF
IF (PRESENT(vab2)) THEN
da_max = 1
db_max = 1
END IF
la_max = la_max_set + da_max
la_min = MAX(0, la_min_set - da_max)
@ -456,6 +472,17 @@ CONTAINS
END DO
END DO
! DFPT for APTs
IF (PRESENT(vab2_work)) THEN
DO j = ncoset(lb_min_set - 1) + 1, ncoset(lb_max_set)
DO i = ncoset(la_min_set - 1) + 1, ncoset(la_max_set)
vab2_work(na + i, nb + j, 1) = vab2_work(na + i, nb + j, 1) + fs(i, j, 2)
vab2_work(na + i, nb + j, 2) = vab2_work(na + i, nb + j, 2) + fs(i, j, 3)
vab2_work(na + i, nb + j, 3) = vab2_work(na + i, nb + j, 3) + fs(i, j, 4)
END DO
END DO
END IF
! *** Calculate the force contribution for the atomic center a ***
IF (calculate_force_a) THEN
@ -510,6 +537,29 @@ CONTAINS
END DO
END DO
! DFPT for APTs
IF (PRESENT(vab2_work)) THEN
DO j = ncoset(lb_min_set - 1) + 1, ncoset(lb_max_set)
DO i = ncoset(la_min_set - 1) + 1, ncoset(la_max_set)
vab2_work(na + i, nb + j, 4) = vab2_work(na + i, nb + j, 4) + fs(i, j, 2)
vab2_work(na + i, nb + j, 5) = vab2_work(na + i, nb + j, 5) + fs(i, j, 3)
vab2_work(na + i, nb + j, 6) = vab2_work(na + i, nb + j, 6) + fs(i, j, 4)
END DO
END DO
DO idir = 1, 3
DO j = ncoset(lb_min_set - 1) + 1, ncoset(lb_max_set)
DO i = ncoset(la_min_set - 1) + 1, ncoset(la_max_set)
vab2(na + i, nb + j, idir) = vab2(na + i, nb + j, idir) &
+ vab2_work(na + i, nb + j, idir)*deltaR(idir, iatom) &
- vab2_work(na + i, nb + j, idir)*deltaR(idir, katom) &
+ vab2_work(na + i, nb + j, idir + 3)*deltaR(idir, jatom) &
- vab2_work(na + i, nb + j, idir + 3)*deltaR(idir, katom)
END DO
END DO
END DO
END IF
! *** Calculate the force contribution for the atomic center b ***
IF (calculate_force_b) THEN

21
src/aobasis/ai_overlap_ppl.F
View file

@ -85,13 +85,21 @@ CONTAINS
!> \param force_a ...
!> \param force_b ...
!> \param fs ...
!> \param hab2 The derivative of the ppl integrals according to the weighting factors deltaR
!> \param hab2_work ...
!> \param deltaR Weighting factors for the derivatives wrt. nuclear positions
!> \param iatom ...
!> \param jatom ...
!> \param katom ...
!> \date May 2011
!> \author Juerg Hutter
!> \version 1.0
!> \note Extended by the derivatives for DFPT [Sandra Luber, Edward Ditler, 2021]
! **************************************************************************************************
SUBROUTINE ppl_integral(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
lb_max_set, lb_min_set, npgfb, rpgfb, zetb, nexp_ppl, alpha_ppl, nct_ppl, cexp_ppl, rpgfc, &
rab, dab, rac, dac, rbc, dbc, vab, s, pab, force_a, force_b, fs)
rab, dab, rac, dac, rbc, dbc, vab, s, pab, force_a, force_b, fs, &
hab2, hab2_work, deltaR, iatom, jatom, katom)
INTEGER, INTENT(IN) :: la_max_set, la_min_set, npgfa
REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
INTEGER, INTENT(IN) :: lb_max_set, lb_min_set, npgfb
@ -113,7 +121,10 @@ CONTAINS
OPTIONAL :: pab
REAL(KIND=dp), DIMENSION(3), INTENT(OUT), OPTIONAL :: force_a, force_b
REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT), &
OPTIONAL :: fs
OPTIONAL :: fs, hab2, hab2_work
REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
OPTIONAL :: deltaR
INTEGER, INTENT(IN), OPTIONAL :: iatom, jatom, katom
INTEGER :: iexp, ij, ipgf, jpgf, mmax, nexp
REAL(KIND=dp) :: rho, sab, t, zetc
@ -127,6 +138,8 @@ CONTAINS
mmax = la_max_set + lb_max_set + 2
force_a(:) = 0.0_dp
force_b(:) = 0.0_dp
ELSE IF (PRESENT(hab2)) THEN
mmax = la_max_set + lb_max_set + 2
ELSE
mmax = la_max_set + lb_max_set
END IF
@ -162,7 +175,9 @@ CONTAINS
CALL os_3center(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
lb_max_set, lb_min_set, npgfb, rpgfb, zetb, auxint, rpgfc, &
rab, dab, rac, dac, rbc, dbc, vab, s, pab, force_a, force_b, fs)
rab, dab, rac, dac, rbc, dbc, vab, s, pab, force_a, force_b, fs, &
vab2=hab2, vab2_work=hab2_work, &
deltaR=deltaR, iatom=iatom, jatom=jatom, katom=katom)
DEALLOCATE (auxint)

18
src/common/bibliography.F
View file

@ -89,7 +89,7 @@ MODULE bibliography
Wilhelm2016a, Wilhelm2016b, Wilhelm2017, Wilhelm2018, Lass2018, cp2kqs2020, &
Behler2007, Behler2011, Schran2020a, Schran2020b, &
Rycroft2009, Thomas2015, Brehm2018, Brehm2020, Shigeta2001, Heinecke2016, &
Brehm2021, Bussy2021a, Bussy2021b
Brehm2021, Bussy2021a, Bussy2021b, Ditler2021
CONTAINS
@ -4703,6 +4703,22 @@ CONTAINS
"ER"), &
DOI="10.3390/molecules26071875")
CALL add_reference(key=Ditler2021, ISI_record=s2a( &
"TY JOUR", &
"PT J", &
"AU Ditler, Edward", &
" Kumar, Chandan", &
" Luber, Sandra", &
"TI Analytic calculation and analysis of atomic polar tensors", &
" for molecules and materials using the Gaussian and plane waves approach", &
"SO The Journal of Chemical Physics", &
"PY 2021", &
"VL 154", &
"AR 104121", &
"DI 10.1063/5.0041056", &
"ER"), &
DOI="10.1063/5.0041056")
END SUBROUTINE add_all_references
END MODULE bibliography

115
src/core_ppl.F
View file

@ -13,6 +13,7 @@
!> - OpenMP added [Iain Bethune, Fiona Reid, 2013-11-13]
!> - Bug fix: correct orbital pointer range [07.2014,JGH]
!> - k-point aware [07.2015,JGH]
!> - Extended by the derivatives for DFPT [Sandra Luber, Edward Ditler, 2021]
! **************************************************************************************************
MODULE core_ppl
@ -82,10 +83,11 @@ CONTAINS
!> \param nimages ...
!> \param cell_to_index ...
!> \param basis_type ...
!> \param deltaR Weighting factors of the derivatives wrt. nuclear positions
! **************************************************************************************************
SUBROUTINE build_core_ppl(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, basis_type)
nimages, cell_to_index, basis_type, deltaR)
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p
TYPE(qs_force_type), DIMENSION(:), POINTER :: force
@ -101,6 +103,8 @@ CONTAINS
INTEGER, INTENT(IN) :: nimages
INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
CHARACTER(LEN=*), INTENT(IN) :: basis_type
REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
OPTIONAL :: deltaR
CHARACTER(LEN=*), PARAMETER :: routineN = 'build_core_ppl'
INTEGER, PARAMETER :: nexp_max = 30
@ -116,11 +120,14 @@ CONTAINS
nct_lpot, npgfa, npgfb, nsgfa, nsgfb
INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
INTEGER, DIMENSION(nexp_max) :: nct_ppl
LOGICAL :: dokp, ecp_local, found, lpotextended
LOGICAL :: do_dR, dokp, ecp_local, found, &
lpotextended
REAL(KIND=dp) :: alpha, dab, dac, dbc, f0, ppl_radius
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: work
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: ppl_fwork, ppl_work
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: qab, work
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: hab2_w, ppl_fwork, ppl_work
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: hab, pab
REAL(KIND=dp), ALLOCATABLE, &
DIMENSION(:, :, :, :, :) :: hab2
REAL(KIND=dp), DIMENSION(1:10) :: aloc, bloc
REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, rab, rac, rbc
REAL(KIND=dp), DIMENSION(3, 3) :: pv_thread
@ -130,7 +137,8 @@ CONTAINS
TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
TYPE(gth_potential_type), POINTER :: gth_potential
REAL(KIND=dp), DIMENSION(nexp_max) :: alpha_ppl
REAL(KIND=dp), DIMENSION(:, :), POINTER :: cval_lpot, h_block, p_block, rpgfa, &
REAL(KIND=dp), DIMENSION(:, :), POINTER :: cval_lpot, h1_1block, h1_2block, &
h1_3block, h_block, p_block, rpgfa, &
rpgfb, sphi_a, sphi_b, zeta, zetb
REAL(KIND=dp), DIMENSION(:), POINTER :: a_local, alpha_lpot, c_local, cexp_ppl, &
set_radius_a, set_radius_b
@ -144,6 +152,7 @@ CONTAINS
!$ INTEGER(KIND=int_8) :: iatom8
!$ INTEGER, PARAMETER :: nlock = 501
do_dR = PRESENT(deltaR)
MARK_USED(int_8)
IF (calculate_forces) THEN
@ -204,13 +213,14 @@ CONTAINS
!$OMP SHARED (ap_iterator, basis_set_list, calculate_forces, use_virial, &
!$OMP matrix_h, matrix_p, atomic_kind_set, qs_kind_set, particle_set, &
!$OMP sab_orb, sac_ppl, nthread, ncoset, nkind, cell_to_index, &
!$OMP ldsab, maxnset, maxder, &
!$OMP ldsab, maxnset, maxder, do_dR, deltaR, &
!$OMP maxlgto, nder, maxco, dokp, locks, natom) &
!$OMP PRIVATE (ikind, jkind, iatom, jatom, rab, basis_set_a, basis_set_b, &
!$OMP first_sgfa, la_max, la_min, npgfa, nsgfa, sphi_a, &
!$OMP zeta, first_sgfb, lb_max, lb_min, npgfb, nsetb, rpgfb, set_radius_b, sphi_b, &
!$OMP zetb, dab, irow, icol, h_block, found, iset, ncoa, &
!$OMP sgfa, jset, ncob, sgfb, nsgfb, p_block, work, pab, hab, kkind, nseta, &
!$OMP sgfa, jset, ncob, sgfb, nsgfb, p_block, work, pab, hab, hab2, hab2_w, qab, &
!$OMP h1_1block, h1_2block, h1_3block, kkind, nseta, &
!$OMP gth_potential, sgp_potential, alpha, cexp_ppl, lpotextended, &
!$OMP ppl_radius, nexp_lpot, nexp_ppl, alpha_ppl, alpha_lpot, nct_ppl, &
!$OMP nct_lpot, cval_ppl, cval_lpot, rac, dac, rbc, dbc, &
@ -243,6 +253,13 @@ CONTAINS
!$OMP DO SCHEDULE(GUIDED)
DO slot = 1, sab_orb(1)%nl_size
!SL
IF (do_dR) THEN
ALLOCATE (hab2(ldsab, ldsab, 4, maxnset, maxnset))
ALLOCATE (hab2_w(ldsab, ldsab, 6))
ALLOCATE (qab(ldsab, ldsab))
ALLOCATE (ppl_fwork(ldai, ldai, maxder))
END IF
ikind = sab_orb(1)%nlist_task(slot)%ikind
jkind = sab_orb(1)%nlist_task(slot)%jkind
@ -306,6 +323,18 @@ CONTAINS
icol = iatom
END IF
NULLIFY (h_block)
IF (do_dR) THEN
NULLIFY (h1_1block, h1_2block, h1_3block)
CALL dbcsr_get_block_p(matrix=matrix_h(1, img)%matrix, &
row=irow, col=icol, BLOCK=h1_1block, found=found)
CALL dbcsr_get_block_p(matrix=matrix_h(2, img)%matrix, &
row=irow, col=icol, BLOCK=h1_2block, found=found)
CALL dbcsr_get_block_p(matrix=matrix_h(3, img)%matrix, &
row=irow, col=icol, BLOCK=h1_3block, found=found)
END IF
CALL dbcsr_get_block_p(matrix_h(1, img)%matrix, irow, icol, h_block, found)
CPASSERT(found)
IF (calculate_forces) THEN
@ -336,6 +365,7 @@ CONTAINS
END IF
hab = 0._dp
IF (do_dr) hab2 = 0._dp
! loop over all kinds for pseudopotential atoms
DO kkind = 1, nkind
@ -441,7 +471,18 @@ CONTAINS
CALL virial_pair_force(pv_thread, f0, force_a, rac)
CALL virial_pair_force(pv_thread, f0, force_b, rbc)
END IF
ELSEIF (do_dR) THEN
hab2_w = 0._dp
CALL ppl_integral( &
la_max(iset), la_min(iset), npgfa(iset), &
rpgfa(:, iset), zeta(:, iset), &
lb_max(jset), lb_min(jset), npgfb(jset), &
rpgfb(:, jset), zetb(:, jset), &
nexp_ppl, alpha_ppl, nct_ppl, cval_ppl, ppl_radius, &
rab, dab, rac, dac, rbc, dbc, &
vab=hab(:, :, iset, jset), s=ppl_work, &
hab2=hab2(:, :, :, iset, jset), hab2_work=hab2_w, fs=ppl_fwork, &
deltaR=deltaR, iatom=iatom, jatom=jatom, katom=katom)
ELSE
CALL ppl_integral( &
la_max(iset), la_min(iset), npgfa(iset), &
@ -457,6 +498,7 @@ CONTAINS
END DO
! *** Contract PPL integrals
IF (.NOT. do_dR) THEN
DO iset = 1, nseta
ncoa = npgfa(iset)*ncoset(la_max(iset))
sgfa = first_sgfa(1, iset)
@ -483,7 +525,62 @@ CONTAINS
END DO
END DO
END DO
ELSE ! do_dr == .true.
DO iset = 1, nseta
ncoa = npgfa(iset)*ncoset(la_max(iset))
sgfa = first_sgfa(1, iset)
DO jset = 1, nsetb
ncob = npgfb(jset)*ncoset(lb_max(jset))
sgfb = first_sgfb(1, jset)
work(1:ncoa, 1:nsgfb(jset)) = MATMUL(hab2(1:ncoa, 1:ncob, 1, iset, jset), &
sphi_b(1:ncob, sgfb:sgfb + nsgfb(jset) - 1))
!$OMP CRITICAL(h1_1block_critical)
IF (iatom <= jatom) THEN
h1_1block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1) = &
h1_1block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1) + &
MATMUL(TRANSPOSE(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1)), work(1:ncoa, 1:nsgfb(jset)))
ELSE
h1_1block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) = &
h1_1block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) + &
MATMUL(TRANSPOSE(work(1:ncoa, 1:nsgfb(jset))), sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1))
END IF
!$OMP END CRITICAL(h1_1block_critical)
work(1:ncoa, 1:nsgfb(jset)) = MATMUL(hab2(1:ncoa, 1:ncob, 2, iset, jset), &
sphi_b(1:ncob, sgfb:sgfb + nsgfb(jset) - 1))
!$OMP CRITICAL(h1_2block_critical)
IF (iatom <= jatom) THEN
h1_2block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1) = &
h1_2block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1) + &
MATMUL(TRANSPOSE(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1)), work(1:ncoa, 1:nsgfb(jset)))
ELSE
h1_2block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) = &
h1_2block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) + &
MATMUL(TRANSPOSE(work(1:ncoa, 1:nsgfb(jset))), sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1))
END IF
!$OMP END CRITICAL(h1_2block_critical)
work(1:ncoa, 1:nsgfb(jset)) = MATMUL(hab2(1:ncoa, 1:ncob, 3, iset, jset), &
sphi_b(1:ncob, sgfb:sgfb + nsgfb(jset) - 1))
!$OMP CRITICAL(h1_3block_critical)
IF (iatom <= jatom) THEN
h1_3block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1) = &
h1_3block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1) + &
MATMUL(TRANSPOSE(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1)), work(1:ncoa, 1:nsgfb(jset)))
ELSE
h1_3block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) = &
h1_3block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) + &
MATMUL(TRANSPOSE(work(1:ncoa, 1:nsgfb(jset))), sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1))
END IF
!$OMP END CRITICAL(h1_3block_critical)
END DO
END DO
END IF
IF (do_dR) DEALLOCATE (qab, hab2, ppl_fwork, hab2_w)
END DO ! slot
DEALLOCATE (hab, work, ppl_work)
IF (calculate_forces) THEN

103
src/core_ppnl.F
View file

@ -10,6 +10,7 @@
!> \par History
!> - refactered from qs_core_hamiltian [Joost VandeVondele, 2008-11-01]
!> - full rewrite [jhu, 2009-01-23]
!> - Extended by the derivatives for DFPT [Sandra Luber, Edward Ditler, 2021]
! **************************************************************************************************
MODULE core_ppnl
USE ai_overlap, ONLY: overlap
@ -78,10 +79,11 @@ CONTAINS
!> \param nimages ...
!> \param cell_to_index ...
!> \param basis_type ...
!> \param deltaR Weighting factors of the derivatives wrt. nuclear positions
! **************************************************************************************************
SUBROUTINE build_core_ppnl(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, basis_type)
nimages, cell_to_index, basis_type, deltaR)
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p
TYPE(qs_force_type), DIMENSION(:), POINTER :: force
@ -98,6 +100,8 @@ CONTAINS
INTEGER, INTENT(IN) :: nimages
INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
CHARACTER(LEN=*), INTENT(IN) :: basis_type
REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
OPTIONAL :: deltaR
CHARACTER(LEN=*), PARAMETER :: routineN = 'build_core_ppnl'
@ -110,7 +114,7 @@ CONTAINS
INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgfa, nprj_ppnl, &
nsgf_seta
INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
LOGICAL :: dogth, dokp, found, ppnl_present
LOGICAL :: do_dR, dogth, dokp, found, ppnl_present
LOGICAL, DIMENSION(0:9) :: is_nonlocal
REAL(KIND=dp) :: dac, f0, ppnl_radius
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: radp
@ -126,8 +130,9 @@ CONTAINS
TYPE(clist_type), POINTER :: clist
TYPE(alist_type), POINTER :: alist_ac, alist_bc
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: achint, acint, bchint, bcint, c_nl
REAL(KIND=dp), DIMENSION(:, :), POINTER :: cprj, h_block, h_nl, p_block, rpgfa, &
sphi_a, vprj_ppnl, zeta
REAL(KIND=dp), DIMENSION(:, :), POINTER :: cprj, h_block, h_nl, p_block, r_2block, &
r_3block, rpgfa, sphi_a, vprj_ppnl, &
zeta
REAL(KIND=dp), DIMENSION(:), POINTER :: a_nl, alpha_ppnl, hprj, set_radius_a
REAL(KIND=dp), DIMENSION(3, SIZE(particle_set)) :: force_thread
TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
@ -142,6 +147,9 @@ CONTAINS
MARK_USED(int_8)
do_dR = .FALSE.
IF (PRESENT(deltaR)) do_dR = .TRUE.
IF (calculate_forces) THEN
CALL timeset(routineN//"_forces", handle)
ELSE
@ -432,12 +440,14 @@ CONTAINS
!$OMP DEFAULT (NONE) &
!$OMP SHARED (dokp, basis_set, matrix_h, cell_to_index,&
!$OMP sab_orb, matrix_p, sap_int, nkind, eps_ppnl, force, &
!$OMP locks, natom, virial, use_virial, calculate_forces) &
!$OMP do_dR, deltaR, maxder, nder, &
!$OMP locks, virial, use_virial, calculate_forces) &
!$OMP PRIVATE (ikind, jkind, iatom, jatom, cell_b, rab, &
!$OMP slot, iab, atom_a, f0, irow, icol, h_block, &
!$OMP r_2block, r_3block, &
!$OMP found,p_block, iac, ibc, alist_ac, alist_bc, acint, bcint, &
!$OMP achint, bchint, na, np, nb, katom, j, fa, fb, rbc, rac, &
!$OMP kkind, kac, kbc, i, img, hash, iatom8) &
!$OMP kkind, kac, kbc, i, img, hash, iatom8, natom) &
!$OMP REDUCTION (+ : pv_thread, force_thread )
!$OMP SINGLE
@ -488,6 +498,13 @@ CONTAINS
END IF
NULLIFY (h_block)
CALL dbcsr_get_block_p(matrix_h(1, img)%matrix, irow, icol, h_block, found)
IF (do_dR) THEN
NULLIFY (r_2block, r_3block)
CALL dbcsr_get_block_p(matrix_h(2, img)%matrix, irow, icol, r_2block, found)
CALL dbcsr_get_block_p(matrix_h(3, img)%matrix, irow, icol, r_3block, found)
END IF
IF (calculate_forces) THEN
NULLIFY (p_block)
CALL dbcsr_get_block_p(matrix_p(1, img)%matrix, irow, icol, p_block, found)
@ -519,12 +536,14 @@ CONTAINS
!$ iatom8 = INT(iatom - 1, int_8)*INT(natom, int_8) + INT(jatom, int_8)
!$ hash = INT(MOD(iatom8, INT(nlock, int_8)) + 1)
!$ CALL omp_set_lock(locks(hash))
IF (iatom <= jatom) THEN
h_block(1:na, 1:nb) = h_block(1:na, 1:nb) + &
MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 1)))
ELSE
h_block(1:nb, 1:na) = h_block(1:nb, 1:na) + &
MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 1)))
IF (.NOT. do_dR) THEN
IF (iatom <= jatom) THEN
h_block(1:na, 1:nb) = h_block(1:na, 1:nb) + &
MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 1)))
ELSE
h_block(1:nb, 1:na) = h_block(1:nb, 1:na) + &
MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 1)))
END IF
END IF
!$ CALL omp_unset_lock(locks(hash))
IF (calculate_forces) THEN
@ -557,6 +576,66 @@ CONTAINS
END IF
END IF
END IF
IF (do_dR) THEN
i = 1; j = 2;
katom = alist_ac%clist(kac)%catom
IF (iatom <= jatom) THEN
h_block(1:na, 1:nb) = h_block(1:na, 1:nb) + &
(deltaR(i, iatom) - deltaR(i, katom))* &
MATMUL(acint(1:na, 1:np, j), TRANSPOSE(bchint(1:nb, 1:np, 1)))
h_block(1:na, 1:nb) = h_block(1:na, 1:nb) + &
(deltaR(i, jatom) - deltaR(i, katom))* &
MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, j)))
ELSE
h_block(1:nb, 1:na) = h_block(1:nb, 1:na) + &
(deltaR(i, iatom) - deltaR(i, katom))* &
MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, j)))
h_block(1:nb, 1:na) = h_block(1:nb, 1:na) + &
(deltaR(i, jatom) - deltaR(i, katom))* &
MATMUL(bcint(1:nb, 1:np, j), TRANSPOSE(achint(1:na, 1:np, 1)))
END IF
i = 2; j = 3;
katom = alist_ac%clist(kac)%catom
IF (iatom <= jatom) THEN
r_2block(1:na, 1:nb) = r_2block(1:na, 1:nb) + &
(deltaR(i, iatom) - deltaR(i, katom))* &
MATMUL(acint(1:na, 1:np, j), TRANSPOSE(bchint(1:nb, 1:np, 1)))
r_2block(1:na, 1:nb) = r_2block(1:na, 1:nb) + &
(deltaR(i, jatom) - deltaR(i, katom))* &
MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, j)))
ELSE
r_2block(1:nb, 1:na) = r_2block(1:nb, 1:na) + &
(deltaR(i, iatom) - deltaR(i, katom))* &
MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, j)))
r_2block(1:nb, 1:na) = r_2block(1:nb, 1:na) + &
(deltaR(i, jatom) - deltaR(i, katom))* &
MATMUL(bcint(1:nb, 1:np, j), TRANSPOSE(achint(1:na, 1:np, 1)))
END IF
i = 3; j = 4;
katom = alist_ac%clist(kac)%catom
IF (iatom <= jatom) THEN
r_3block(1:na, 1:nb) = r_3block(1:na, 1:nb) + &
(deltaR(i, iatom) - deltaR(i, katom))* &
MATMUL(acint(1:na, 1:np, j), TRANSPOSE(bchint(1:nb, 1:np, 1)))
r_3block(1:na, 1:nb) = r_3block(1:na, 1:nb) + &
(deltaR(i, jatom) - deltaR(i, katom))* &
MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, j)))
ELSE
r_3block(1:nb, 1:na) = r_3block(1:nb, 1:na) + &
(deltaR(i, iatom) - deltaR(i, katom))* &
MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, j)))
r_3block(1:nb, 1:na) = r_3block(1:nb, 1:na) + &
(deltaR(i, jatom) - deltaR(i, katom))* &
MATMUL(bcint(1:nb, 1:np, j), TRANSPOSE(achint(1:na, 1:np, 1)))
END IF
END IF
EXIT ! We have found a match and there can be only one single match
END IF
END DO

11
src/grid/common/grid_constants.h
View file

@ -33,7 +33,16 @@ enum grid_func {
GRID_FUNC_DZDX = 703,
GRID_FUNC_DXDX = 801,
GRID_FUNC_DYDY = 802,
GRID_FUNC_DZDZ = 803
GRID_FUNC_DZDZ = 803,
GRID_FUNC_DAB_X = 901,
GRID_FUNC_DAB_Y = 902,
GRID_FUNC_DAB_Z = 903,
GRID_FUNC_ADB_X = 904,
GRID_FUNC_ADB_Y = 905,
GRID_FUNC_ADB_Z = 906,
GRID_FUNC_CORE_X = 1001,
GRID_FUNC_CORE_Y = 1002,
GRID_FUNC_CORE_Z = 1003,
};
enum grid_backend {

99
src/grid/common/grid_prepare_pab.h
View file

@ -131,6 +131,69 @@ GRID_DEVICE static void prepare_pab_DABpADB(const int idir, const orbital a,
prep_term(up(idir, a), b, -2.0 * zeta * pab_val, n, cab);
}
/*******************************************************************************
* \brief Implementation of function GRID_FUNC_DAB_{X,Y,Z}.
* This function takes the derivates with respect to nuclear positions
* which results in a change of signs compared to prepare_pab_DABpADB.
* Only the derivative with respect to the primitive on the left.
* \author Edward Ditler
******************************************************************************/
GRID_DEVICE static void prepare_pab_DAB(const int idir, const orbital a,
const orbital b, const double zeta,
const double pab_val, const int n,
double *cab) {
// creates cab such that mapping it with pgf_a pgf_b
// is equivalent to mapping pab with (nabla_{idir} pgf_a) pgf_b
// (ddX pgf_a)( pgf_b ) =
// (-lax pgf_{a-1x} - 2*zeta*pgf_{a+1x}) pgf_b
prep_term(down(idir, a), b, -a.l[idir] * pab_val, n, cab);
prep_term(up(idir, a), b, +2.0 * zeta * pab_val, n, cab);
}
/*******************************************************************************
* \brief Implementation of function GRID_FUNC_ADB_{X,Y,Z}.
* This function takes the derivates with respect to nuclear positions
* which results in a change of signs compared to prepare_pab_DABpADB.
* Only the derivative with respect to the primitive on the right.
* \author Edward Ditler
******************************************************************************/
GRID_DEVICE static void prepare_pab_ADB(const int idir, const orbital a,
const orbital b, const double zetb,
const double pab_val, const int n,
double *cab) {
// creates cab such that mapping it with pgf_a pgf_b
// is equivalent to mapping pab with
// pgf_a (nabla_{idir} pgf_b) + (nabla_{idir} pgf_a) pgf_b
// ( pgf_a ) (ddX pgf_b) =
// pgf_a *(-lbx pgf_{b-1x} - 2*zetb*pgf_{b+1x})
prep_term(a, down(idir, b), -b.l[idir] * pab_val, n, cab);
prep_term(a, up(idir, b), +2.0 * zetb * pab_val, n, cab);
}
/*******************************************************************************
* \brief Implementation of function GRID_FUNC_CORE_{X,Y,Z}.
* This function takes the derivates with respect to nuclear positions.
* \author Edward Ditler
******************************************************************************/
GRID_DEVICE static void prepare_pab_core(const int idir, const orbital a,
const orbital b, const double zeta,
const double pab_val, const int n,
double *cab) {
// creates cab such that mapping it with pgf_a pgf_b
// is equivalent to mapping pab with (nabla_{idir} pgf_a) pgf_b
// (ddX pgf_a)( pgf_b ) = 2*zeta*pgf_{a+1x}) pgf_b
prep_term(a, down(idir, b), 0.0, n, cab);
prep_term(a, up(idir, b), 0.0, n, cab);
prep_term(down(idir, a), b, 0.0, n, cab);
prep_term(up(idir, a), b, +2.0 * zeta * pab_val, n, cab);
}
/*******************************************************************************
* \brief Implementation of function GRID_FUNC_{DX,DY,DZ}.
* \author Ole Schuett
@ -319,6 +382,33 @@ GRID_DEVICE static void prepare_pab(const enum grid_func func, const orbital a,
case GRID_FUNC_DABpADB_Z:
prepare_pab_DABpADB(2, a, b, zeta, zetb, pab_val, n, cab);
break;
case GRID_FUNC_DAB_X:
prepare_pab_DAB(0, a, b, zeta, pab_val, n, cab);
break;
case GRID_FUNC_DAB_Y:
prepare_pab_DAB(1, a, b, zeta, pab_val, n, cab);
break;
case GRID_FUNC_DAB_Z:
prepare_pab_DAB(2, a, b, zeta, pab_val, n, cab);
break;
case GRID_FUNC_ADB_X:
prepare_pab_ADB(0, a, b, zetb, pab_val, n, cab);
break;
case GRID_FUNC_ADB_Y:
prepare_pab_ADB(1, a, b, zetb, pab_val, n, cab);
break;
case GRID_FUNC_ADB_Z:
prepare_pab_ADB(2, a, b, zetb, pab_val, n, cab);
break;
case GRID_FUNC_CORE_X:
prepare_pab_core(0, a, b, zeta, pab_val, n, cab);
break;
case GRID_FUNC_CORE_Y:
prepare_pab_core(1, a, b, zeta, pab_val, n, cab);
break;
case GRID_FUNC_CORE_Z:
prepare_pab_core(2, a, b, zeta, pab_val, n, cab);
break;
case GRID_FUNC_DX:
prepare_pab_Di(0, a, b, zeta, zetb, pab_val, n, cab);
break;
@ -383,6 +473,15 @@ static prepare_ldiffs prepare_get_ldiffs(const enum grid_func func) {
case GRID_FUNC_DABpADB_X:
case GRID_FUNC_DABpADB_Y:
case GRID_FUNC_DABpADB_Z:
case GRID_FUNC_DAB_X:
case GRID_FUNC_DAB_Y:
case GRID_FUNC_DAB_Z:
case GRID_FUNC_ADB_X:
case GRID_FUNC_ADB_Y:
case GRID_FUNC_ADB_Z:
case GRID_FUNC_CORE_X:
case GRID_FUNC_CORE_Y:
case GRID_FUNC_CORE_Z:
ldiffs.la_max_diff = +1;
ldiffs.la_min_diff = -1;
ldiffs.lb_max_diff = +1;

32
src/grid/grid_api.F
View file

@ -55,6 +55,16 @@ MODULE grid_api
INTEGER, PARAMETER, PUBLIC :: GRID_FUNC_DXDX = 801
INTEGER, PARAMETER, PUBLIC :: GRID_FUNC_DYDY = 802
INTEGER, PARAMETER, PUBLIC :: GRID_FUNC_DZDZ = 803
INTEGER, PARAMETER, PUBLIC :: GRID_FUNC_DAB_X = 901
INTEGER, PARAMETER, PUBLIC :: GRID_FUNC_DAB_Y = 902
INTEGER, PARAMETER, PUBLIC :: GRID_FUNC_DAB_Z = 903
INTEGER, PARAMETER, PUBLIC :: GRID_FUNC_ADB_X = 904
INTEGER, PARAMETER, PUBLIC :: GRID_FUNC_ADB_Y = 905
INTEGER, PARAMETER, PUBLIC :: GRID_FUNC_ADB_Z = 906
INTEGER, PARAMETER, PUBLIC :: GRID_FUNC_CORE_X = 1001
INTEGER, PARAMETER, PUBLIC :: GRID_FUNC_CORE_Y = 1002
INTEGER, PARAMETER, PUBLIC :: GRID_FUNC_CORE_Z = 1003
INTEGER, PARAMETER, PUBLIC :: GRID_BACKEND_AUTO = 10
INTEGER, PARAMETER, PUBLIC :: GRID_BACKEND_REF = 11
@ -262,7 +272,8 @@ CONTAINS
!> \param use_virial ...
!> \param my_virial_a ...
!> \param my_virial_b ...
!> \param hdab ...
!> \param hdab Derivative with respect to the primitive on the left.
!> \param hadb Derivative with respect to the primitive on the right.
!> \param a_hdab ...
!> \param use_subpatch ...
!> \param subpatch_pattern ...
@ -275,7 +286,7 @@ CONTAINS
calculate_forces, force_a, force_b, &
compute_tau, &
use_virial, my_virial_a, &
my_virial_b, hdab, a_hdab, use_subpatch, subpatch_pattern)
my_virial_b, hdab, hadb, a_hdab, use_subpatch, subpatch_pattern)
INTEGER, INTENT(IN) :: la_max
REAL(KIND=dp), INTENT(IN) :: zeta
@ -296,7 +307,7 @@ CONTAINS
LOGICAL, INTENT(IN), OPTIONAL :: compute_tau, use_virial
REAL(KIND=dp), DIMENSION(3, 3), OPTIONAL :: my_virial_a, my_virial_b
REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, &
POINTER :: hdab
POINTER :: hdab, hadb
REAL(KIND=dp), DIMENSION(:, :, :, :), OPTIONAL, &
POINTER :: a_hdab
LOGICAL, OPTIONAL :: use_subpatch
@ -310,8 +321,8 @@ CONTAINS
REAL(KIND=dp), DIMENSION(3, 2), TARGET :: forces
REAL(KIND=dp), DIMENSION(3, 3, 2), TARGET :: virials
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: grid
TYPE(C_PTR) :: a_hdab_cptr, forces_cptr, hdab_cptr, &
pab_cptr, virials_cptr
TYPE(C_PTR) :: a_hdab_cptr, forces_cptr, hadb_cptr, &
hdab_cptr, pab_cptr, virials_cptr
INTERFACE
SUBROUTINE grid_ref_integrate_pgf_product_c(orthorhombic, compute_tau, &
border_mask, &
@ -319,7 +330,7 @@ CONTAINS
zeta, zetb, dh, dh_inv, ra, rab, &
npts_global, npts_local, shift_local, border_width, &
radius, o1, o2, n1, n2, grid, hab, pab, &
forces, virials, hdab, a_hdab) &
forces, virials, hdab, hadb, a_hdab) &
BIND(C, name="grid_ref_integrate_pgf_product")
IMPORT :: C_PTR, C_INT, C_DOUBLE, C_BOOL
LOGICAL(KIND=C_BOOL), VALUE :: orthorhombic
@ -350,6 +361,7 @@ CONTAINS
TYPE(C_PTR), VALUE :: forces
TYPE(C_PTR), VALUE :: virials
TYPE(C_PTR), VALUE :: hdab
TYPE(C_PTR), VALUE :: hadb
TYPE(C_PTR), VALUE :: a_hdab
END SUBROUTINE grid_ref_integrate_pgf_product_c
END INTERFACE
@ -400,12 +412,17 @@ CONTAINS
END IF
IF (calculate_forces .AND. PRESENT(hdab)) THEN
hdab(:, :, :) = 0.0_dp
hdab_cptr = C_LOC(hdab(1, 1, 1))
ELSE
hdab_cptr = C_NULL_PTR
END IF
IF (calculate_forces .AND. PRESENT(hadb)) THEN
hadb_cptr = C_LOC(hadb(1, 1, 1))
ELSE
hadb_cptr = C_NULL_PTR
END IF
IF (calculate_forces .AND. my_use_virial .AND. PRESENT(a_hdab)) THEN
a_hdab_cptr = C_LOC(a_hdab(1, 1, 1, 1))
ELSE
@ -473,6 +490,7 @@ CONTAINS
forces=forces_cptr, &
virials=virials_cptr, &
hdab=hdab_cptr, &
hadb=hadb_cptr, &
a_hdab=a_hdab_cptr)
IF (PRESENT(force_a) .AND. C_ASSOCIATED(forces_cptr)) &

2
src/grid/grid_replay.c
View file

@ -405,7 +405,7 @@ double grid_replay(const char *filename, const int cycles, const bool collocate,
lb_min, zeta, zetb, dh, dh_inv, ra, rab, npts_global, npts_local,
shift_local, border_width, radius, o1, o2, n1, n2,
grid_ref->host_buffer, hab_test, pab, forces_test, virials_test,
NULL, NULL);
NULL, NULL, NULL);
}
for (int i = 0; i < 3; i++) {
for (int j = 0; j < 3; j++) {

2
src/grid/ref/grid_ref_collocate.c
View file

@ -155,7 +155,7 @@ static void record_collocate(
la_min, lb_max, lb_min, zeta, zetb, dh, dh_inv,
ra, rab, npts_global, npts_local, shift_local,
border_width, radius, o1, o2, n1, n2, grid,
hab, pab, forces, virials, NULL, NULL);
hab, pab, forces, virials, NULL, NULL, NULL);
for (int i = o2; i < ncoset(lb_max) + o2; i++) {
for (int j = o1; j < ncoset(la_max) + o1; j++) {

12
src/grid/ref/grid_ref_integrate.c
View file

@ -31,7 +31,8 @@ void grid_ref_integrate_pgf_product(
const int border_width[3], const double radius, const int o1, const int o2,
const int n1, const int n2, const double *grid, double hab[n2][n1],
const double pab[n2][n1], double forces[2][3], double virials[2][3][3],
double hdab[n2][n1][3], double a_hdab[n2][n1][3][3]) {
double hdab[n2][n1][3], double hadb[n2][n1][3],
double a_hdab[n2][n1][3][3]) {
const bool calculate_forces = (forces != NULL || hdab != NULL);
const bool calculate_virial = (virials != NULL || a_hdab != NULL);
@ -97,7 +98,7 @@ void grid_ref_integrate_pgf_product(
}
}
// Update hdab and a_hdab (not used in batch mode).
// Update hdab, hadb, and a_hdab (not used in batch mode).
if (hdab != NULL) {
assert(!compute_tau);
for (int i = 0; i < 3; i++) {
@ -105,6 +106,13 @@ void grid_ref_integrate_pgf_product(
get_force_a(a, b, i, zeta, zetb, m1, cab, false);
}
}
if (hadb != NULL) {
assert(!compute_tau);
for (int i = 0; i < 3; i++) {
hadb[o2 + idx(b)][o1 + idx(a)][i] +=
get_force_b(a, b, i, zeta, zetb, rab, m1, cab, false);
}
}
if (a_hdab != NULL) {
assert(!compute_tau);
for (int i = 0; i < 3; i++) {

3
src/grid/ref/grid_ref_integrate.h
View file

@ -36,7 +36,8 @@ void grid_ref_integrate_pgf_product(
const int border_width[3], const double radius, const int o1, const int o2,
const int n1, const int n2, const double *grid, double hab[n2][n1],
const double pab[n2][n1], double forces[2][3], double virials[2][3][3],
double hdab[n2][n1][3], double a_hdab[n2][n1][3][3]);
double hdab[n2][n1][3], double hadb[n2][n1][3],
double a_hdab[n2][n1][3][3]);
#endif
// EOF

1
src/grid/ref/grid_ref_task_list.c
View file

@ -567,6 +567,7 @@ static void integrate_one_grid_level(
/*forces=*/(forces != NULL) ? my_forces : NULL,
/*virials=*/(virial != NULL) ? my_virials : NULL,
/*hdab=*/NULL,
/*hadb=*/NULL,
/*a_hdab=*/NULL);
} // end of task loop

113
src/input_cp2k_properties_dft.F
View file

@ -251,6 +251,10 @@ CONTAINS
CALL section_add_subsection(section, subsection)
CALL section_release(subsection)
CALL create_dcdr_section(subsection)
CALL section_add_subsection(section, subsection)
CALL section_release(subsection)
CALL section_create(subsection, __LOCATION__, name="PRINT", &
description="printing of information during the linear response calculation", &
repeats=.FALSE.)
@ -278,6 +282,115 @@ CONTAINS
END SUBROUTINE create_linres_section
! **************************************************************************************************
!> \brief creates the input structure used to activate
!> calculation of position perturbation DFPT
!> \param section ...
!> \author Sandra Luber, Edward Ditler
! **************************************************************************************************
SUBROUTINE create_dcdr_section(section)
TYPE(section_type), POINTER :: section
CHARACTER(len=*), PARAMETER :: routineN = 'create_dcdr_section', &
routineP = moduleN//':'//routineN
LOGICAL :: failure
TYPE(keyword_type), POINTER :: keyword
TYPE(section_type), POINTER :: print_key, subsection
failure = .FALSE.
NULLIFY (keyword, print_key, subsection)
CPASSERT(.NOT. ASSOCIATED(section))
IF (.NOT. failure) THEN
CALL section_create(section, __LOCATION__, name="DCDR", &
description="Compute analytical gradients the dipole moments.", &
n_keywords=50, n_subsections=1, repeats=.FALSE.)
CALL keyword_create(keyword, __LOCATION__, name="_SECTION_PARAMETERS_", &
description="controls the activation of the APT calculation", &
usage="&DCDR 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="LIST_OF_ATOMS", &
description="Specifies a list of atoms.", &
usage="LIST {integer} {integer} .. {integer}", repeats=.TRUE., &
n_var=-1, type_of_var=integer_t)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="DISTRIBUTED_ORIGIN", &
variants=(/"DO_GAUGE"/), &
description="Use the distributed origin (DO) gauge?", &
usage="DISTRIBUTED_ORIGIN 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="ORBITAL_CENTER", &
description="The orbital center.", &
usage="ORBITAL_CENTER WANNIER", &
default_i_val=current_orb_center_wannier, &
enum_c_vals=s2a("WANNIER", "COMMON", "ATOM", "BOX"), &
enum_desc=s2a("Use the Wannier centers.", &
"Use a common center (works only for an isolate molecule).", &
"Use the atoms as center.", &
"Boxing."), &
enum_i_vals=(/current_orb_center_wannier, current_orb_center_common, &
current_orb_center_atom, current_orb_center_box/))
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="REFERENCE", &
description="Gauge origin of the velocity gauge factor.", &
enum_c_vals=s2a("COM", "COAC", "USER_DEFINED", "ZERO"), &
enum_desc=s2a("Use Center of Mass", &
"Use Center of Atomic Charges", &
"Use User-defined Point", &
"Use Origin of Coordinate System"), &
enum_i_vals=(/use_mom_ref_com, &
use_mom_ref_coac, &
use_mom_ref_user, &
use_mom_ref_zero/), &
default_i_val=use_mom_ref_zero)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="REFERENCE_POINT", &
description="User-defined reference point of the velocity gauge factor.", &
usage="REFERENCE_POINT x y z", &
repeats=.FALSE., n_var=3, type_of_var=real_t, unit_str='bohr')
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
NULLIFY (subsection)
CALL section_create(subsection, __LOCATION__, name="PRINT", &
description="print results of the magnetic dipole moment calculation", &
repeats=.FALSE.)
CALL cp_print_key_section_create(print_key, __LOCATION__, "APT", &
description="Controls the printing of the electric dipole gradient", &
print_level=low_print_level, add_last=add_last_numeric, filename="")
CALL section_add_subsection(subsection, print_key)
CALL section_release(print_key)
CALL section_add_subsection(section, subsection)
CALL section_release(subsection)
NULLIFY (subsection)
CALL create_interp_section(subsection)
CALL section_add_subsection(section, subsection)
CALL section_release(subsection)
END IF
END SUBROUTINE create_dcdr_section
! **************************************************************************************************
!> \brief creates the input structure used to activate
!> calculation of induced current DFPT

556
src/qs_collocate_density.F
View file

@ -29,6 +29,7 @@
!> - introduced map_consistent (Joost 02.04)
!> - Addition of the subroutine calculate_atomic_charge_density (TdK, 08.05)
!> - rewrite of the collocate/integrate kernels (Joost VandeVondele, 03.07)
!> - Extended by the derivatives for DFPT [Sandra Luber, Edward Ditler, 2021]
!> \author Matthias Krack (03.04.2001)
!> 1) Joost VandeVondele (01.2002)
!> Thomas D. Kuehne (04.08.2005)
@ -58,10 +59,11 @@ MODULE qs_collocate_density
USE gaussian_gridlevels, ONLY: gaussian_gridlevel,&
gridlevel_info_type
USE grid_api, ONLY: &
GRID_FUNC_AB, GRID_FUNC_DABpADB_X, GRID_FUNC_DABpADB_Y, GRID_FUNC_DABpADB_Z, &
GRID_FUNC_DADB, GRID_FUNC_DX, GRID_FUNC_DXDX, GRID_FUNC_DXDY, GRID_FUNC_DY, &
GRID_FUNC_DYDY, GRID_FUNC_DYDZ, GRID_FUNC_DZ, GRID_FUNC_DZDX, GRID_FUNC_DZDZ, &
collocate_pgf_product, grid_collocate_task_list
GRID_FUNC_AB, GRID_FUNC_CORE_X, GRID_FUNC_CORE_Y, GRID_FUNC_CORE_Z, GRID_FUNC_DAB_X, &
GRID_FUNC_DAB_Y, GRID_FUNC_DAB_Z, GRID_FUNC_DABpADB_X, GRID_FUNC_DABpADB_Y, &
GRID_FUNC_DABpADB_Z, GRID_FUNC_DADB, GRID_FUNC_DX, GRID_FUNC_DXDX, GRID_FUNC_DXDY, &
GRID_FUNC_DY, GRID_FUNC_DYDY, GRID_FUNC_DYDZ, GRID_FUNC_DZ, GRID_FUNC_DZDX, &
GRID_FUNC_DZDZ, collocate_pgf_product, grid_collocate_task_list
USE input_constants, ONLY: &
orb_dx2, orb_dxy, orb_dy2, orb_dyz, orb_dz2, orb_dzx, orb_px, orb_py, orb_pz, orb_s
USE kinds, ONLY: default_string_length,&
@ -134,7 +136,9 @@ MODULE qs_collocate_density
calculate_rho_elec, &
calculate_drho_elec, &
calculate_wavefunction, &
calculate_rho_nlcc
calculate_rho_nlcc, &
calculate_drho_elec_dR, &
calculate_drho_core
INTEGER :: debug_count = 0
@ -972,6 +976,145 @@ CONTAINS
END SUBROUTINE calculate_rho_core
! *****************************************************************************
!> \brief Computes the derivative of the density of the core charges with
!> respect to the nuclear coordinates on the grid.
!> \param drho_core The resulting density derivative
!> \param qs_env ...
!> \param beta Derivative direction
!> \param lambda Atom index
!> \note SL November 2014, ED 2021
! **************************************************************************************************
SUBROUTINE calculate_drho_core(drho_core, qs_env, beta, lambda)
TYPE(pw_p_type), INTENT(INOUT) :: drho_core
TYPE(qs_environment_type), POINTER :: qs_env
INTEGER, INTENT(IN) :: beta, lambda
CHARACTER(len=*), PARAMETER :: routineN = 'calculate_drho_core', &
routineP = moduleN//':'//routineN
INTEGER :: atom_a, dabqadb_func, handle, iatom, &
ikind, ithread, j, natom, npme, &
nthread, subpatch_pattern
INTEGER, DIMENSION(:), POINTER :: atom_list, cores
REAL(KIND=dp) :: alpha, eps_rho_rspace, radius
REAL(KIND=dp), DIMENSION(3) :: ra
REAL(KIND=dp), DIMENSION(:, :), POINTER :: pab
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(cell_type), POINTER :: cell
TYPE(cube_info_type) :: cube_info
TYPE(dft_control_type), POINTER :: dft_control
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(pw_env_type), POINTER :: pw_env
TYPE(pw_p_type) :: rhoc_r
TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(realspace_grid_type), POINTER :: rs_rho
CALL timeset(routineN, handle)
NULLIFY (cell, dft_control, pab, atomic_kind_set, qs_kind_set, particle_set, &
atom_list, pw_env, rs_rho, auxbas_pw_pool, cores)
ALLOCATE (pab(1, 1))
CALL get_qs_env(qs_env=qs_env, &
atomic_kind_set=atomic_kind_set, &
qs_kind_set=qs_kind_set, &
cell=cell, &
dft_control=dft_control, &
particle_set=particle_set, &
pw_env=pw_env)
CALL pw_env_get(pw_env, auxbas_rs_grid=rs_rho, &
auxbas_pw_pool=auxbas_pw_pool)
cube_info = pw_env%cube_info(1)
! be careful in parallel nsmax is chosen with multigrid in mind!
CALL rs_grid_retain(rs_rho)
CALL rs_grid_zero(rs_rho)
eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
SELECT CASE (beta)
CASE (1)
dabqadb_func = GRID_FUNC_CORE_X
CASE (2)
dabqadb_func = GRID_FUNC_CORE_Y
CASE (3)
dabqadb_func = GRID_FUNC_CORE_Z
CASE DEFAULT
PRINT *, 'beta', beta
CPABORT("invalid beta")
END SELECT
DO ikind = 1, SIZE(atomic_kind_set)
CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=atom_list)
CALL get_qs_kind(qs_kind_set(ikind), &
alpha_core_charge=alpha, ccore_charge=pab(1, 1))
IF (alpha == 0.0_dp .OR. pab(1, 1) == 0.0_dp) CYCLE
nthread = 1
ithread = 0
CALL reallocate(cores, 1, natom)
npme = 0
cores = 0
DO iatom = 1, natom
IF (rs_rho%desc%parallel .AND. .NOT. rs_rho%desc%distributed) THEN
! replicated realspace grid, split the atoms up between procs
IF (MODULO(iatom, rs_rho%desc%group_size) == rs_rho%desc%my_pos) THEN
npme = npme + 1
cores(npme) = iatom
END IF
ELSE
npme = npme + 1
cores(npme) = iatom
END IF
END DO
IF (npme .GT. 0) THEN
DO j = 1, npme
iatom = cores(j)
atom_a = atom_list(iatom)
IF (atom_a /= lambda) CYCLE
ra(:) = pbc(particle_set(atom_a)%r, cell)
subpatch_pattern = 0
radius = exp_radius_very_extended(la_min=0, la_max=0, &
lb_min=0, lb_max=0, &
ra=ra, rb=ra, rp=ra, &
zetp=alpha, eps=eps_rho_rspace, &
pab=pab, o1=0, o2=0, & ! without map_consistent
prefactor=-1.0_dp, cutoff=0.0_dp)
CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, &
(/0.0_dp, 0.0_dp, 0.0_dp/), -1.0_dp, pab, 0, 0, rs_rho, &
cell, cube_info, radius=radius, ga_gb_function=dabqadb_func, &
use_subpatch=.TRUE., subpatch_pattern=subpatch_pattern)
END DO
END IF
END DO
IF (ASSOCIATED(cores)) THEN
DEALLOCATE (cores)
END IF
DEALLOCATE (pab)
CALL pw_pool_create_pw(auxbas_pw_pool, rhoc_r%pw, &
use_data=REALDATA3D, in_space=REALSPACE)
CALL rs_pw_transfer(rs_rho, rhoc_r%pw, rs2pw)
CALL rs_grid_release(rs_rho)
CALL pw_transfer(rhoc_r%pw, drho_core%pw)
CALL pw_pool_give_back_pw(auxbas_pw_pool, rhoc_r%pw)
CALL timestop(handle)
END SUBROUTINE calculate_drho_core
! **************************************************************************************************
!> \brief collocate a single Gaussian on the grid
!> \param rho_gb charge density generated by a single gaussian
@ -1938,6 +2081,409 @@ CONTAINS
END SUBROUTINE calculate_drho_elec
! **************************************************************************************************
!> \brief Computes the gradient wrt. nuclear coordinates of a density on the grid
!> The density is given in terms of the density matrix_p
!> \param matrix_p Density matrix
!> \param matrix_p_kp ...
!> \param drho Density gradient on the grid
!> \param drho_gspace Density gradient on the reciprocal grid
!> \param qs_env ...
!> \param soft_valid ...
!> \param basis_type ...
!> \param beta Derivative direction
!> \param lambda Atom index
!> \note SL, ED 2021
!> Adapted from calculate_drho_elec
! **************************************************************************************************
SUBROUTINE calculate_drho_elec_dR(matrix_p, matrix_p_kp, drho, drho_gspace, qs_env, &
soft_valid, basis_type, beta, lambda)
TYPE(dbcsr_type), OPTIONAL, POINTER :: matrix_p
TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
POINTER :: matrix_p_kp
TYPE(pw_p_type), INTENT(INOUT) :: drho, drho_gspace
TYPE(qs_environment_type), POINTER :: qs_env
LOGICAL, INTENT(IN), OPTIONAL :: soft_valid
CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: basis_type
INTEGER, INTENT(IN) :: beta, lambda
CHARACTER(len=*), PARAMETER :: routineN = 'calculate_drho_elec_dR'
CHARACTER(LEN=default_string_length) :: my_basis_type
INTEGER :: bcol, brow, dabqadb_func, handle, i, iatom, iatom_old, igrid_level, ikind, &
ikind_old, img, img_old, ipgf, iset, iset_old, itask, ithread, jatom, jatom_old, jkind, &
jkind_old, jpgf, jset, jset_old, maxco, maxsgf_set, na1, na2, natoms, nb1, nb2, ncoa, &
ncob, nimages, nseta, nsetb, ntasks, nthread, sgfa, sgfb
INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
npgfb, nsgfa, nsgfb
INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
LOGICAL :: atom_pair_changed, distributed_rs_grids, &
do_kp, found, my_soft, use_subpatch
REAL(KIND=dp) :: eps_rho_rspace, f, prefactor, radius, &
scale, zetp
REAL(KIND=dp), DIMENSION(3) :: ra, rab, rab_inv, rb, rp
REAL(KIND=dp), DIMENSION(:, :), POINTER :: p_block, pab, sphi_a, sphi_b, work, &
zeta, zetb
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: pabt, workt
TYPE(atom_pair_type), DIMENSION(:), POINTER :: atom_pair_recv, atom_pair_send
TYPE(cell_type), POINTER :: cell
TYPE(cube_info_type), DIMENSION(:), POINTER :: cube_info
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: deltap
TYPE(dft_control_type), POINTER :: dft_control
TYPE(gridlevel_info_type), POINTER :: gridlevel_info
TYPE(gto_basis_set_type), POINTER :: orb_basis_set
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(pw_env_type), POINTER :: pw_env
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
POINTER :: rs_descs
TYPE(realspace_grid_p_type), DIMENSION(:), POINTER :: rs_rho
TYPE(task_list_type), POINTER :: task_list, task_list_soft
TYPE(task_type), DIMENSION(:), POINTER :: tasks
CALL timeset(routineN, handle)
CPASSERT(PRESENT(matrix_p) .OR. PRESENT(matrix_p_kp))
do_kp = PRESENT(matrix_p_kp)
NULLIFY (cell, dft_control, orb_basis_set, deltap, qs_kind_set, &
particle_set, rs_rho, pw_env, rs_descs, la_max, la_min, lb_max, &
lb_min, npgfa, npgfb, nsgfa, nsgfb, p_block, sphi_a, sphi_b, &
zeta, zetb, first_sgfa, first_sgfb, tasks, pabt, workt)
debug_count = debug_count + 1
! by default, the full density is calculated
my_soft = .FALSE.
IF (PRESENT(soft_valid)) my_soft = soft_valid
IF (PRESENT(basis_type)) THEN
my_basis_type = basis_type
ELSE
my_basis_type = "ORB"
END IF
CALL get_qs_env(qs_env=qs_env, &
qs_kind_set=qs_kind_set, &
cell=cell, &
dft_control=dft_control, &
particle_set=particle_set, &
pw_env=pw_env)
SELECT CASE (my_basis_type)
CASE ("ORB")
CALL get_qs_env(qs_env=qs_env, &
task_list=task_list, &
task_list_soft=task_list_soft)
CASE ("AUX_FIT")
CALL get_qs_env(qs_env=qs_env, &
task_list_aux_fit=task_list, &
task_list_soft=task_list_soft)
END SELECT
! *** assign from pw_env
gridlevel_info => pw_env%gridlevel_info
cube_info => pw_env%cube_info
! *** Allocate work storage ***
nthread = 1
CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
maxco=maxco, &
maxsgf_set=maxsgf_set, &
basis_type=my_basis_type)
CALL reallocate(pabt, 1, maxco, 1, maxco, 0, nthread - 1)
CALL reallocate(workt, 1, maxco, 1, maxsgf_set, 0, nthread - 1)
! find maximum numbers
nimages = dft_control%nimages
CPASSERT(nimages == 1 .OR. do_kp)
natoms = SIZE(particle_set)
! get the task lists
IF (my_soft) task_list => task_list_soft
CPASSERT(ASSOCIATED(task_list))
tasks => task_list%tasks
atom_pair_send => task_list%atom_pair_send
atom_pair_recv => task_list%atom_pair_recv
ntasks = task_list%ntasks
! *** set up the rs multi-grids
CPASSERT(ASSOCIATED(pw_env))
CALL pw_env_get(pw_env, rs_descs=rs_descs, rs_grids=rs_rho)
DO igrid_level = 1, gridlevel_info%ngrid_levels
CALL rs_grid_retain(rs_rho(igrid_level)%rs_grid)
distributed_rs_grids = rs_rho(igrid_level)%rs_grid%desc%distributed
END DO
eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
! *** Initialize working density matrix ***
! distributed rs grids require a matrix that will be changed
! whereas this is not the case for replicated grids
ALLOCATE (deltap(nimages))
IF (distributed_rs_grids) THEN
DO img = 1, nimages
END DO
! this matrix has no strict sparsity pattern in parallel
! deltap%sparsity_id=-1
IF (do_kp) THEN
DO img = 1, nimages
CALL dbcsr_copy(deltap(img)%matrix, matrix_p_kp(img)%matrix, &
name="DeltaP")
END DO
ELSE
CALL dbcsr_copy(deltap(1)%matrix, matrix_p, name="DeltaP")
END IF
ELSE
IF (do_kp) THEN
DO img = 1, nimages
deltap(img)%matrix => matrix_p_kp(img)%matrix
END DO
ELSE
deltap(1)%matrix => matrix_p
END IF
END IF
! distribute the matrix
IF (distributed_rs_grids) THEN
CALL rs_distribute_matrix(rs_descs=rs_descs, pmats=deltap, &
atom_pair_send=atom_pair_send, atom_pair_recv=atom_pair_recv, &
nimages=nimages, scatter=.TRUE.)
END IF
! map all tasks on the grids
ithread = 0
pab => pabt(:, :, ithread)
work => workt(:, :, ithread)
DO igrid_level = 1, gridlevel_info%ngrid_levels
CALL rs_grid_zero(rs_rho(igrid_level)%rs_grid)
END DO
iatom_old = -1; jatom_old = -1; iset_old = -1; jset_old = -1
ikind_old = -1; jkind_old = -1; img_old = -1
loop_tasks: DO itask = 1, ntasks
!decode the atom pair and basis info
igrid_level = tasks(itask)%grid_level
img = tasks(itask)%image
iatom = tasks(itask)%iatom
jatom = tasks(itask)%jatom
iset = tasks(itask)%iset
jset = tasks(itask)%jset
ipgf = tasks(itask)%ipgf
jpgf = tasks(itask)%jpgf
ikind = particle_set(iatom)%atomic_kind%kind_number
jkind = particle_set(jatom)%atomic_kind%kind_number
IF (iatom .NE. iatom_old .OR. jatom .NE. jatom_old .OR. img .NE. img_old) THEN
IF (iatom .NE. iatom_old) ra(:) = pbc(particle_set(iatom)%r, cell)
IF (iatom <= jatom) THEN
brow = iatom
bcol = jatom
ELSE
brow = jatom
bcol = iatom
END IF
IF (ikind .NE. ikind_old) THEN
CALL get_qs_kind(qs_kind_set(ikind), softb=my_soft, &
basis_set=orb_basis_set, basis_type=my_basis_type)
CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
first_sgf=first_sgfa, &
lmax=la_max, &
lmin=la_min, &
npgf=npgfa, &
nset=nseta, &
nsgf_set=nsgfa, &
sphi=sphi_a, &
zet=zeta)
END IF
IF (jkind .NE. jkind_old) THEN
CALL get_qs_kind(qs_kind_set(jkind), softb=my_soft, &
basis_set=orb_basis_set, basis_type=my_basis_type)
CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
first_sgf=first_sgfb, &
lmax=lb_max, &
lmin=lb_min, &
npgf=npgfb, &
nset=nsetb, &
nsgf_set=nsgfb, &
sphi=sphi_b, &
zet=zetb)
END IF
CALL dbcsr_get_block_p(matrix=deltap(img)%matrix, &
row=brow, col=bcol, BLOCK=p_block, found=found)
CPASSERT(found)
iatom_old = iatom
jatom_old = jatom
ikind_old = ikind
jkind_old = jkind
img_old = img
atom_pair_changed = .TRUE.
ELSE
atom_pair_changed = .FALSE.
END IF
IF (atom_pair_changed .OR. iset_old .NE. iset .OR. jset_old .NE. jset) THEN
ncoa = npgfa(iset)*ncoset(la_max(iset))
sgfa = first_sgfa(1, iset)
ncob = npgfb(jset)*ncoset(lb_max(jset))
sgfb = first_sgfb(1, jset)
IF (iatom <= jatom) THEN
CALL dgemm("N", "N", ncoa, nsgfb(jset), nsgfa(iset), &
1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
p_block(sgfa, sgfb), SIZE(p_block, 1), &
0.0_dp, work(1, 1), maxco)
CALL dgemm("N", "T", ncoa, ncob, nsgfb(jset), &
1.0_dp, work(1, 1), maxco, &
sphi_b(1, sgfb), SIZE(sphi_b, 1), &
0.0_dp, pab(1, 1), maxco)
ELSE
CALL dgemm("N", "N", ncob, nsgfa(iset), nsgfb(jset), &
1.0_dp, sphi_b(1, sgfb), SIZE(sphi_b, 1), &
p_block(sgfb, sgfa), SIZE(p_block, 1), &
0.0_dp, work(1, 1), maxco)
CALL dgemm("N", "T", ncob, ncoa, nsgfa(iset), &
1.0_dp, work(1, 1), maxco, &
sphi_a(1, sgfa), SIZE(sphi_a, 1), &
0.0_dp, pab(1, 1), maxco)
END IF
iset_old = iset
jset_old = jset
END IF
rab(:) = tasks(itask)%rab
rb(:) = ra(:) + rab(:)
zetp = zeta(ipgf, iset) + zetb(jpgf, jset)
f = zetb(jpgf, jset)/zetp
rp(:) = ra(:) + f*rab(:)
prefactor = EXP(-zeta(ipgf, iset)*f*DOT_PRODUCT(rab, rab))
radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
lb_min=lb_min(jset), lb_max=lb_max(jset), &
ra=ra, rb=rb, rp=rp, &
zetp=zetp, eps=eps_rho_rspace, &
prefactor=prefactor, cutoff=1.0_dp)
na1 = (ipgf - 1)*ncoset(la_max(iset)) + 1
na2 = ipgf*ncoset(la_max(iset))
nb1 = (jpgf - 1)*ncoset(lb_max(jset)) + 1
nb2 = jpgf*ncoset(lb_max(jset))
! takes the density matrix symmetry in account, i.e. off-diagonal blocks need to be mapped 'twice'
IF (iatom == jatom .AND. img == 1) THEN
scale = 1.0_dp
ELSE
scale = 2.0_dp
END IF
! check whether we need to use fawzi's generalised collocation scheme
IF (rs_rho(igrid_level)%rs_grid%desc%distributed) THEN
!tasks(4,:) is 0 for replicated, 1 for distributed 2 for exceptional distributed tasks
IF (tasks(itask)%dist_type .EQ. 2) THEN
use_subpatch = .TRUE.
ELSE
use_subpatch = .FALSE.
END IF
ELSE
use_subpatch = .FALSE.
END IF
SELECT CASE (beta)
CASE (1)
dabqadb_func = GRID_FUNC_DAB_X
CASE (2)
dabqadb_func = GRID_FUNC_DAB_Y
CASE (3)
dabqadb_func = GRID_FUNC_DAB_Z
CASE DEFAULT
CPABORT("invalid beta")
END SELECT
IF (iatom <= jatom) THEN
IF (iatom == lambda) &
CALL collocate_pgf_product( &
la_max(iset), zeta(ipgf, iset), la_min(iset), &
lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
ra, rab, scale, pab, na1 - 1, nb1 - 1, &
rsgrid=rs_rho(igrid_level)%rs_grid, cell=cell, cube_info=cube_info(igrid_level), &
ga_gb_function=dabqadb_func, radius=radius, &
use_subpatch=use_subpatch, &
subpatch_pattern=tasks(itask)%subpatch_pattern)
IF (jatom == lambda) &
CALL collocate_pgf_product( &
la_max(iset), zeta(ipgf, iset), la_min(iset), &
lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
ra, rab, scale, pab, na1 - 1, nb1 - 1, &
rsgrid=rs_rho(igrid_level)%rs_grid, cell=cell, cube_info=cube_info(igrid_level), &
ga_gb_function=dabqadb_func + 3, radius=radius, &
use_subpatch=use_subpatch, &
subpatch_pattern=tasks(itask)%subpatch_pattern)
ELSE
rab_inv = -rab
IF (jatom == lambda) &
CALL collocate_pgf_product( &
lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
la_max(iset), zeta(ipgf, iset), la_min(iset), &
rb, rab_inv, scale, pab, nb1 - 1, na1 - 1, &
rs_rho(igrid_level)%rs_grid, cell, cube_info(igrid_level), &
ga_gb_function=dabqadb_func, radius=radius, &
use_subpatch=use_subpatch, &
subpatch_pattern=tasks(itask)%subpatch_pattern)
IF (iatom == lambda) &
CALL collocate_pgf_product( &
lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
la_max(iset), zeta(ipgf, iset), la_min(iset), &
rb, rab_inv, scale, pab, nb1 - 1, na1 - 1, &
rs_rho(igrid_level)%rs_grid, cell, cube_info(igrid_level), &
ga_gb_function=dabqadb_func + 3, radius=radius, &
use_subpatch=use_subpatch, &
subpatch_pattern=tasks(itask)%subpatch_pattern)
END IF
END DO loop_tasks
CALL density_rs2pw_basic(pw_env, rs_rho, drho, drho_gspace)
! *** Release work storage ***
IF (ASSOCIATED(rs_rho)) THEN
DO i = 1, SIZE(rs_rho)
CALL rs_grid_release(rs_rho(i)%rs_grid)
END DO
END IF
IF (distributed_rs_grids) THEN
CALL dbcsr_deallocate_matrix_set(deltap)
ELSE
DO img = 1, nimages
NULLIFY (deltap(img)%matrix)
END DO
DEALLOCATE (deltap)
END IF
DEALLOCATE (pabt, workt)
CALL timestop(handle)
END SUBROUTINE calculate_drho_elec_dR
! **************************************************************************************************
!> \brief maps a given wavefunction on the grid
!> \param mo_vectors ...

669
src/qs_dcdr.F Normal file
View file

@ -0,0 +1,669 @@
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright 2000-2021 CP2K developers group <https://cp2k.org> !
! !
! SPDX-License-Identifier: GPL-2.0-or-later !
!--------------------------------------------------------------------------------------------------!
! **************************************************************************************************
!> \brief Calculate the derivatives of the MO coefficients wrt nuclear coordinates
!> \author Sandra Luber, Edward Ditler
! **************************************************************************************************
MODULE qs_dcdr
!#include "./common/cp_common_uses.f90"
USE atomic_kind_types, ONLY: get_atomic_kind
USE cell_types, ONLY: cell_type,&
pbc
USE cp_array_utils, ONLY: cp_2d_r_p_type
USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply,&
dbcsr_allocate_matrix_set,&
dbcsr_deallocate_matrix_set
USE cp_fm_basic_linalg, ONLY: cp_fm_gemm,&
cp_fm_scale,&
cp_fm_scale_and_add,&
cp_fm_trace
USE cp_fm_types, ONLY: cp_fm_create,&
cp_fm_get_diag,&
cp_fm_p_type,&
cp_fm_release,&
cp_fm_set_all,&
cp_fm_to_fm,&
cp_fm_type
USE cp_gemm_interface, ONLY: cp_gemm
USE cp_log_handling, ONLY: cp_get_default_logger,&
cp_logger_type
USE cp_output_handling, ONLY: cp_print_key_finished_output,&
cp_print_key_unit_nr
USE dbcsr_api, ONLY: dbcsr_add,&
dbcsr_copy,&
dbcsr_desymmetrize,&
dbcsr_p_type,&
dbcsr_set
USE input_section_types, ONLY: section_vals_get_subs_vals,&
section_vals_type
USE kinds, ONLY: dp
USE molecule_types, ONLY: molecule_of_atom,&
molecule_type
USE particle_types, ONLY: particle_type
USE qs_dcdr_ao, ONLY: apply_op_constant_term,&
core_dR,&
d_core_charge_density_dR,&
d_vhxc_dR,&
hr_mult_by_delta_1d,&
vhxc_R_perturbed_basis_functions
USE qs_dcdr_utils, ONLY: dcdr_read_restart,&
dcdr_write_restart,&
multiply_localization,&
shift_wannier_into_cell
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
USE qs_kind_types, ONLY: get_qs_kind,&
qs_kind_type
USE qs_linres_methods, ONLY: linres_solver
USE qs_linres_types, ONLY: dcdr_env_type,&
linres_control_type
USE qs_moments, ONLY: build_local_moment_matrix,&
dipole_deriv_ao
USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
USE qs_p_env_types, ONLY: qs_p_env_type
#include "./base/base_uses.f90"
IMPLICIT NONE
PRIVATE
PUBLIC :: prepare_per_atom, dcdr_response_dR, dcdr_build_op_dR, apt_dR, apt_dR_localization
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_dcdr'
CONTAINS
! **************************************************************************************************
!> \brief Prepare the environment for a choice of lambda
!> \param dcdr_env ...
!> \param qs_env ...
!> \author Edward Ditler
! **************************************************************************************************
SUBROUTINE prepare_per_atom(dcdr_env, qs_env)
TYPE(dcdr_env_type) :: dcdr_env
TYPE(qs_environment_type), POINTER :: qs_env
CHARACTER(LEN=*), PARAMETER :: routineN = 'prepare_per_atom', &
routineP = moduleN//':'//routineN
INTEGER :: handle, i, ispin, j, natom
TYPE(neighbor_list_set_p_type), DIMENSION(:), &
POINTER :: sab_all
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
CALL timeset(routineN, handle)
NULLIFY (sab_all, qs_kind_set, particle_set)
CALL get_qs_env(qs_env=qs_env, &
sab_all=sab_all, &
qs_kind_set=qs_kind_set, &
particle_set=particle_set)
natom = SIZE(particle_set)
IF (dcdr_env%distributed_origin) dcdr_env%ref_point(:) = particle_set(dcdr_env%lambda)%r(:)
dcdr_env%delta_basis_function = 0._dp
dcdr_env%delta_basis_function(:, dcdr_env%lambda) = 1._dp
! S matrix
! S1 = - < da/dr | b > * delta_a - < a | db/dr > * delta_b
! matrix_s(2:4) are anti-symmetric matrices and contain derivatives wrt. to < a |
! = < da/dR | b > = - < da/dr | b > = < a | db/dr >
! matrix_s1(2:4) = d/dR < a | b >
! and it's built as
! = - matrix_s * delta_b + matrix_s * delta_a
! = - < da/dR | b > * delta_b + < da/dR | b > * delta_a
! = + < da/dr | b > * delta_b - < da/dr | b > * delta_a
! = - < a | db/dr > * delta_b - < da/dr | b > * delta_a
DO i = 1, 3
! S matrix
CALL dbcsr_set(dcdr_env%matrix_nosym_temp(i)%matrix, 0._dp)
CALL dbcsr_desymmetrize(dcdr_env%matrix_s(1 + i)%matrix, dcdr_env%matrix_s1(1 + i)%matrix)
CALL dbcsr_desymmetrize(dcdr_env%matrix_s(1 + i)%matrix, dcdr_env%matrix_nosym_temp(i)%matrix)
CALL hr_mult_by_delta_1d(dcdr_env%matrix_s1(1 + i)%matrix, qs_kind_set, "ORB", &
sab_all, dcdr_env%lambda, direction_Or=.TRUE.)
CALL hr_mult_by_delta_1d(dcdr_env%matrix_nosym_temp(i)%matrix, qs_kind_set, "ORB", &
sab_all, dcdr_env%lambda, direction_Or=.FALSE.)
CALL dbcsr_add(dcdr_env%matrix_s1(1 + i)%matrix, dcdr_env%matrix_nosym_temp(i)%matrix, -1._dp, +1._dp)
CALL dbcsr_set(dcdr_env%matrix_nosym_temp(i)%matrix, 0._dp)
! T matrix
CALL dbcsr_set(dcdr_env%matrix_nosym_temp(i)%matrix, 0._dp)
CALL dbcsr_desymmetrize(dcdr_env%matrix_t(1 + i)%matrix, dcdr_env%matrix_t1(1 + i)%matrix)
CALL dbcsr_desymmetrize(dcdr_env%matrix_t(1 + i)%matrix, dcdr_env%matrix_nosym_temp(i)%matrix)
CALL hr_mult_by_delta_1d(dcdr_env%matrix_t1(1 + i)%matrix, qs_kind_set, "ORB", &
sab_all, dcdr_env%lambda, direction_Or=.TRUE.)
CALL hr_mult_by_delta_1d(dcdr_env%matrix_nosym_temp(i)%matrix, qs_kind_set, "ORB", &
sab_all, dcdr_env%lambda, direction_Or=.FALSE.)
CALL dbcsr_add(dcdr_env%matrix_t1(1 + i)%matrix, dcdr_env%matrix_nosym_temp(i)%matrix, -1._dp, +1._dp)
CALL dbcsr_set(dcdr_env%matrix_nosym_temp(i)%matrix, 0._dp)
END DO
! Operator:
ispin = 1
DO i = 1, 3
CALL dbcsr_set(dcdr_env%matrix_ppnl_1(i)%matrix, 0.0_dp)
CALL dbcsr_set(dcdr_env%matrix_hc(i)%matrix, 0.0_dp)
CALL dbcsr_set(dcdr_env%matrix_vhxc_perturbed_basis(ispin, i)%matrix, 0.0_dp)
CALL dbcsr_set(dcdr_env%matrix_vhxc_perturbed_basis(ispin, i + 3)%matrix, 0.0_dp)
CALL dbcsr_set(dcdr_env%matrix_d_vhxc_dR(i)%matrix, 0.0_dp)
CALL dbcsr_set(dcdr_env%matrix_core_charge_1(i)%matrix, 0.0_dp)
END DO
CALL core_dR(qs_env, dcdr_env) ! dcdr_env%matrix_ppnl_1, hc
CALL d_vhxc_dR(qs_env, dcdr_env) ! dcdr_env%matrix_d_vhxc_dR
CALL d_core_charge_density_dR(qs_env, dcdr_env) ! dcdr_env%matrix_core_charge_1
CALL vhxc_R_perturbed_basis_functions(qs_env, dcdr_env) ! dcdr_env%matrix_vhxc_perturbed_basis
! APT:
DO i = 1, 3
DO j = 1, 3
CALL dbcsr_set(dcdr_env%matrix_difdip(i, j)%matrix, 0._dp)
END DO
END DO
CALL dipole_deriv_ao(qs_env, dcdr_env%matrix_difdip, dcdr_env%delta_basis_function, 1, dcdr_env%ref_point)
CALL timestop(handle)
END SUBROUTINE prepare_per_atom
! **************************************************************************************************
!> \brief Build the operator for the position perturbation
!> \param dcdr_env ...
!> \param qs_env ...
!> \authors SL, ED
! **************************************************************************************************
SUBROUTINE dcdr_build_op_dR(dcdr_env, qs_env)
TYPE(dcdr_env_type) :: dcdr_env
TYPE(qs_environment_type), POINTER :: qs_env
CHARACTER(LEN=*), PARAMETER :: routineN = 'dcdr_build_op_dR', &
routineP = moduleN//':'//routineN
REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
INTEGER :: handle, nao, nmo
TYPE(cp_fm_type), POINTER :: buf
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: opdr_sym
CALL timeset(routineN, handle)
nao = dcdr_env%nao
nmo = dcdr_env%nmo
! allocate matrix for the sum of the perturbation terms of the operator (dbcsr matrix)
NULLIFY (opdr_sym)
CALL dbcsr_allocate_matrix_set(opdr_sym, 1)
ALLOCATE (opdr_sym(1)%matrix)
CALL dbcsr_copy(opdr_sym(1)%matrix, dcdr_env%matrix_s1(1)%matrix) ! symmetric
CALL dbcsr_set(opdr_sym(1)%matrix, 0.0_dp)
CALL apply_op_constant_term(qs_env, dcdr_env) ! dcdr_env%matrix_apply_op_constant
! Hartree and Exchange-Correlation contributions
CALL dbcsr_add(opdr_sym(1)%matrix, dcdr_env%matrix_core_charge_1(dcdr_env%beta)%matrix, zero, one)
CALL dbcsr_add(opdr_sym(1)%matrix, dcdr_env%matrix_d_vhxc_dR(dcdr_env%beta)%matrix, one, one)
CALL dbcsr_add(opdr_sym(1)%matrix, dcdr_env%matrix_vhxc_perturbed_basis(1, dcdr_env%beta)%matrix, one, one)
! Core Hamiltonian contributions
CALL dbcsr_add(opdr_sym(1)%matrix, dcdr_env%matrix_hc(dcdr_env%beta)%matrix, one, one)
CALL dbcsr_add(opdr_sym(1)%matrix, dcdr_env%matrix_ppnl_1(dcdr_env%beta)%matrix, one, one)
CALL dbcsr_add(opdr_sym(1)%matrix, dcdr_env%matrix_apply_op_constant(1)%matrix, one, one)
CALL dbcsr_desymmetrize(opdr_sym(1)%matrix, dcdr_env%hamiltonian1(1)%matrix)
CALL dbcsr_add(dcdr_env%hamiltonian1(1)%matrix, dcdr_env%matrix_t1(dcdr_env%beta + 1)%matrix, one, one)
CALL cp_dbcsr_sm_fm_multiply(dcdr_env%hamiltonian1(1)%matrix, dcdr_env%mo_coeff(1)%matrix, &
dcdr_env%op_dR(1)%matrix, ncol=nmo)
! The overlap derivative terms for the Sternheimer equation
! buf = mo * (-mo * matrix_ks * mo)
CALL cp_fm_create(buf, dcdr_env%likemos_fm_struct)
CALL cp_gemm('N', 'N', nao, nmo, nmo, &
-1.0_dp, dcdr_env%mo_coeff(1)%matrix, dcdr_env%chc(1)%matrix, &
0.0_dp, buf)
CALL cp_dbcsr_sm_fm_multiply(dcdr_env%matrix_s1(dcdr_env%beta + 1)%matrix, buf, dcdr_env%op_dR(1)%matrix, &
nmo, alpha=1.0_dp, beta=1.0_dp)
CALL cp_fm_release(buf)
! SL multiply by -1 for response solver (H-S<H> C + dR_coupled= - (op_dR)
CALL cp_fm_scale(-1.0_dp, dcdr_env%op_dR(1)%matrix)
CALL dbcsr_deallocate_matrix_set(opdr_sym)
CALL timestop(handle)
END SUBROUTINE dcdr_build_op_dR
! **************************************************************************************************
!> \brief Get the dC/dR by solving the Sternheimer equation, using the op_dR matrix
!> \param dcdr_env ...
!> \param p_env ...
!> \param qs_env ...
!> \authors SL, ED
! **************************************************************************************************
SUBROUTINE dcdr_response_dR(dcdr_env, p_env, qs_env)
TYPE(dcdr_env_type) :: dcdr_env
TYPE(qs_p_env_type), POINTER :: p_env
TYPE(qs_environment_type), POINTER :: qs_env
CHARACTER(LEN=*), PARAMETER :: routineN = 'dcdr_response_dR', &
routineP = moduleN//':'//routineN
INTEGER :: handle, output_unit
LOGICAL :: should_stop
TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: h1_psi0, psi0_order, psi1
TYPE(cp_logger_type), POINTER :: logger
TYPE(linres_control_type), POINTER :: linres_control
TYPE(section_vals_type), POINTER :: lr_section
CALL timeset(routineN, handle)
NULLIFY (linres_control, lr_section, logger, psi1, h1_psi0)
CALL get_qs_env(qs_env=qs_env, &
linres_control=linres_control)
logger => cp_get_default_logger()
lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES")
output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
extension=".linresLog")
IF (output_unit > 0) THEN
WRITE (UNIT=output_unit, FMT="(T10,A,/)") &
"*** Self consistent optimization of the response wavefunction ***"
END IF
! allocate the vectors
ALLOCATE (psi0_order(1), psi1(1), h1_psi0(1))
NULLIFY (psi1(1)%matrix, h1_psi0(1)%matrix)
psi0_order(1)%matrix => dcdr_env%mo_coeff(1)%matrix
CALL cp_fm_create(psi1(1)%matrix, dcdr_env%likemos_fm_struct)
CALL cp_fm_create(h1_psi0(1)%matrix, dcdr_env%likemos_fm_struct)
! Restart
IF (linres_control%linres_restart) THEN
CALL dcdr_read_restart(qs_env, lr_section, psi1, dcdr_env%lambda, dcdr_env%beta, "dCdR")
ELSE
CALL cp_fm_set_all(psi1(1)%matrix, 0.0_dp)
END IF
IF (output_unit > 0) THEN
WRITE (output_unit, "(T10,A,I4,A)") &
"Response to the perturbation operator referring to atom ", dcdr_env%lambda, &
" displaced in "//ACHAR(dcdr_env%beta + 119)
END IF
CALL cp_fm_set_all(dcdr_env%dCR(1)%matrix, 0.0_dp)
CALL cp_fm_to_fm(dcdr_env%op_dR(1)%matrix, h1_psi0(1)%matrix)
linres_control%lr_triplet = .FALSE. ! we do singlet response
linres_control%do_kernel = .TRUE.
linres_control%converged = .FALSE.
! Position perturbation to get dCR
! (H0-E0) psi1 = (H1-E1) psi0
! psi1 = the perturbed wavefunction
! h1_psi0 = (H1-E1-S1*\varepsilon)
! psi0_order = the unperturbed wavefunction
CALL linres_solver(p_env, qs_env, psi1, h1_psi0, psi0_order, &
output_unit, should_stop)
CALL cp_fm_to_fm(psi1(1)%matrix, dcdr_env%dCR(1)%matrix)
! Write the new result to the restart file
IF (linres_control%linres_restart) THEN
CALL dcdr_write_restart(qs_env, lr_section, psi1, dcdr_env%lambda, dcdr_env%beta, "dCdR")
END IF
! clean up
CALL cp_fm_release(psi1(1)%matrix)
CALL cp_fm_release(h1_psi0(1)%matrix)
DEALLOCATE (psi1, h1_psi0, psi0_order)
CALL cp_print_key_finished_output(output_unit, logger, lr_section, &
"PRINT%PROGRAM_RUN_INFO")
CALL timestop(handle)
END SUBROUTINE dcdr_response_dR
! **************************************************************************************************
!> \brief Calculate atomic polar tensor
!> \param qs_env ...
!> \param dcdr_env ...
!> \author Edward Ditler
! **************************************************************************************************
SUBROUTINE apt_dR(qs_env, dcdr_env)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(dcdr_env_type) :: dcdr_env
CHARACTER(LEN=*), PARAMETER :: routineN = 'apt_dR', routineP = moduleN//':'//routineN
INTEGER :: alpha, handle, ikind, nao, nmo
LOGICAL :: ghost
REAL(dp) :: apt_basis_derivative, &
apt_coeff_derivative, charge, f_spin
REAL(dp), DIMENSION(:, :, :), POINTER :: apt_el, apt_nuc
TYPE(cp_fm_p_type), POINTER :: overlap1_MO, tmp_fm_like_mos
TYPE(cp_fm_type), POINTER :: mo_coeff
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
apt_basis_derivative = 0._dp
apt_coeff_derivative = 0._dp
CALL timeset(routineN, handle)
NULLIFY (qs_kind_set, particle_set)
CALL get_qs_env(qs_env=qs_env, &
qs_kind_set=qs_kind_set, &
particle_set=particle_set)
nao = dcdr_env%nao
nmo = dcdr_env%nmo
mo_coeff => dcdr_env%mo_coeff(1)%matrix
apt_el => dcdr_env%apt_el_dcdr
apt_nuc => dcdr_env%apt_nuc_dcdr
f_spin = 2._dp
ALLOCATE (tmp_fm_like_mos)
ALLOCATE (overlap1_MO)
CALL cp_fm_create(tmp_fm_like_mos%matrix, dcdr_env%likemos_fm_struct)
CALL cp_fm_create(overlap1_MO%matrix, dcdr_env%momo_fm_struct)
! Compute S^(1,R)_(ij)
CALL cp_fm_set_all(tmp_fm_like_mos%matrix, 0.0_dp)
CALL cp_fm_scale_and_add(0._dp, dcdr_env%dCR_prime(1)%matrix, 1._dp, dcdr_env%dCR(1)%matrix)
CALL cp_dbcsr_sm_fm_multiply(dcdr_env%matrix_s1(dcdr_env%beta + 1)%matrix, mo_coeff, &
tmp_fm_like_mos%matrix, ncol=nmo)
CALL cp_fm_gemm("T", "N", nmo, nmo, nao, &
1.0_dp, mo_coeff, tmp_fm_like_mos%matrix, &
0.0_dp, overlap1_MO%matrix)
! C^1 <- -dCR - 0.5 * mo_coeff @ S1_ij
! We get the negative of the coefficients out of the linres solver
! And apply the constant correction due to the overlap derivative.
CALL cp_fm_gemm("N", "N", nao, nmo, nmo, &
-0.5_dp, mo_coeff, overlap1_MO%matrix, &
-1.0_dp, dcdr_env%dCR_prime(1)%matrix)
DO alpha = 1, 3
! FIRST CONTRIBUTION: dCR * moments * mo
CALL cp_fm_set_all(tmp_fm_like_mos%matrix, 0._dp)
CALL dbcsr_desymmetrize(dcdr_env%matrix_s1(1)%matrix, dcdr_env%matrix_nosym_temp(1)%matrix)
CALL dbcsr_desymmetrize(dcdr_env%moments(alpha)%matrix, dcdr_env%matrix_nosym_temp(2)%matrix)
CALL dbcsr_add(dcdr_env%matrix_nosym_temp(1)%matrix, dcdr_env%matrix_nosym_temp(2)%matrix, &
-dcdr_env%ref_point(alpha), 1._dp)
CALL cp_dbcsr_sm_fm_multiply(dcdr_env%matrix_nosym_temp(1)%matrix, dcdr_env%dCR_prime(1)%matrix, &
tmp_fm_like_mos%matrix, ncol=nmo)
CALL cp_fm_trace(mo_coeff, tmp_fm_like_mos%matrix, apt_coeff_derivative)
apt_coeff_derivative = (-2._dp)*f_spin*apt_coeff_derivative
apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) &
= apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) + apt_coeff_derivative
END DO ! alpha
DO alpha = 1, 3
! SECOND CONTRIBUTION: We assemble all combinations of r_i, d(chi)/d(idir)
! difdip contains derivatives with respect to atom dcdr_env%lambda
! difdip(alpha, beta): < a | r_alpha | db/dR_beta >
! Multiply by the MO coefficients
CALL cp_fm_set_all(tmp_fm_like_mos%matrix, 0.0_dp)
CALL cp_dbcsr_sm_fm_multiply(dcdr_env%matrix_difdip(alpha, dcdr_env%beta)%matrix, mo_coeff, &
tmp_fm_like_mos%matrix, ncol=nmo)
CALL cp_fm_trace(mo_coeff, tmp_fm_like_mos%matrix, apt_basis_derivative)
! The negative sign is due to dipole_deriv_ao computing the derivatives with respect to nuclear coordinates.
apt_basis_derivative = -f_spin*apt_basis_derivative
apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) = &
apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) + apt_basis_derivative
END DO ! alpha
! Finally the nuclear contribution: nuclear charge * Kronecker_delta_{dcdr_env%beta,i}
CALL get_atomic_kind(particle_set(dcdr_env%lambda)%atomic_kind, kind_number=ikind)
CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge, ghost=ghost)
IF (.NOT. ghost) THEN
apt_nuc(dcdr_env%beta, dcdr_env%beta, dcdr_env%lambda) = &
apt_nuc(dcdr_env%beta, dcdr_env%beta, dcdr_env%lambda) + charge
END IF
! And deallocate all the things!
CALL cp_fm_release(tmp_fm_like_mos%matrix)
CALL cp_fm_release(overlap1_MO%matrix)
DEALLOCATE (tmp_fm_like_mos)
DEALLOCATE (overlap1_MO)
CALL timestop(handle)
END SUBROUTINE apt_dR
! **************************************************************************************************
!> \brief Calculate atomic polar tensor using the localized dipole operator
!> \param qs_env ...
!> \param dcdr_env ...
!> \author Edward Ditler
! **************************************************************************************************
SUBROUTINE apt_dR_localization(qs_env, dcdr_env)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(dcdr_env_type) :: dcdr_env
CHARACTER(LEN=*), PARAMETER :: routineN = 'apt_dR_localization', &
routineP = moduleN//':'//routineN
INTEGER :: alpha, handle, i, icenter, ikind, &
map_atom, map_molecule, nao, natom, &
nmo, nsubset
INTEGER, ALLOCATABLE, DIMENSION(:) :: mapping_atom_molecule, &
mapping_wannier_atom
LOGICAL :: ghost
REAL(dp) :: apt_basis_derivative, &
apt_coeff_derivative, charge, f_spin, &
smallest_r, this_factor, tmp_aptcontr, &
tmp_r
REAL(dp), ALLOCATABLE, DIMENSION(:) :: diagonal_elements
REAL(dp), DIMENSION(3) :: distance, r_shifted
REAL(dp), DIMENSION(:, :, :), POINTER :: apt_el, apt_nuc
REAL(dp), DIMENSION(:, :, :, :), POINTER :: apt_center, apt_subset
TYPE(cell_type), POINTER :: cell
TYPE(cp_2d_r_p_type), DIMENSION(:), POINTER :: centers_set
TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: tmp_fm_like_mos, tmp_fm_momo
TYPE(cp_fm_p_type), POINTER :: tmp_fm
TYPE(cp_fm_type), POINTER :: mo_coeff
TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
CALL timeset(routineN, handle)
NULLIFY (qs_kind_set, particle_set, molecule_set, cell)
CALL get_qs_env(qs_env=qs_env, &
qs_kind_set=qs_kind_set, &
particle_set=particle_set, &
molecule_set=molecule_set, &
cell=cell)
nsubset = SIZE(molecule_set)
natom = SIZE(particle_set)
apt_el => dcdr_env%apt_el_dcdr
apt_nuc => dcdr_env%apt_nuc_dcdr
apt_subset => dcdr_env%apt_el_dcdr_per_subset
apt_center => dcdr_env%apt_el_dcdr_per_center
! Map wannier functions to atoms
centers_set => dcdr_env%centers_set
ALLOCATE (mapping_wannier_atom(dcdr_env%nbr_center(1)))
DO icenter = 1, dcdr_env%nbr_center(1)
! For every center we check which atom is closest
CALL shift_wannier_into_cell(r=centers_set(1)%array(1:3, icenter), &
cell=cell, &
r_shifted=r_shifted)
smallest_r = HUGE(0._dp)
DO i = 1, natom
distance = pbc(r_shifted, particle_set(i)%r(1:3), cell)
tmp_r = SUM(distance**2)
IF (tmp_r < smallest_r) THEN
mapping_wannier_atom(icenter) = i
smallest_r = tmp_r
END IF
END DO
END DO
! Map atoms to molecules
ALLOCATE (mapping_atom_molecule(natom))
CALL molecule_of_atom(molecule_set, atom_to_mol=mapping_atom_molecule)
IF (dcdr_env%lambda == 1 .AND. dcdr_env%beta == 1) THEN
! PRINT *, 'Mapping of atoms to subsets'
DO icenter = 1, dcdr_env%nbr_center(1)
map_atom = mapping_wannier_atom(icenter)
map_molecule = mapping_atom_molecule(map_atom)
! PRINT *, 'Center', icenter, 'map_atom', map_atom
END DO
END IF
nao = dcdr_env%nao
nmo = dcdr_env%nmo
mo_coeff => dcdr_env%mo_coeff(1)%matrix
f_spin = 2._dp
ALLOCATE (tmp_fm, tmp_fm_momo(3), tmp_fm_like_mos(3))
ALLOCATE (diagonal_elements(nmo))
! Allocate temporary matrices
CALL cp_fm_create(tmp_fm%matrix, dcdr_env%likemos_fm_struct)
DO i = 1, 3
CALL cp_fm_create(tmp_fm_momo(i)%matrix, dcdr_env%momo_fm_struct)
CALL cp_fm_create(tmp_fm_like_mos(i)%matrix, dcdr_env%likemos_fm_struct)
END DO
! Build the full coefficient derivatives.
CALL cp_fm_set_all(tmp_fm%matrix, 0.0_dp)
CALL cp_fm_scale_and_add(0._dp, dcdr_env%dCR_prime(1)%matrix, 1._dp, dcdr_env%dCR(1)%matrix)
CALL cp_dbcsr_sm_fm_multiply(dcdr_env%matrix_s1(dcdr_env%beta + 1)%matrix, mo_coeff, &
tmp_fm%matrix, ncol=nmo)
CALL cp_fm_gemm("T", "N", nmo, nmo, nao, &
1.0_dp, mo_coeff, tmp_fm%matrix, &
0.0_dp, tmp_fm_momo(1)%matrix)
! C^1 <- -dCR - 0.5 * mo_coeff @ S1_ij
CALL cp_fm_gemm("N", "N", nao, nmo, nmo, &
-0.5_dp, mo_coeff, tmp_fm_momo(1)%matrix, &
-1.0_dp, dcdr_env%dCR_prime(1)%matrix)
! FIRST CONTRIBUTION: dCR * moments * mo
this_factor = -2._dp*f_spin
DO alpha = 1, 3
DO icenter = 1, dcdr_env%nbr_center(1)
CALL build_local_moment_matrix(qs_env, dcdr_env%moments, 1, &
ref_point=centers_set(1)%array(1:3, icenter))
CALL multiply_localization(ao_matrix=dcdr_env%moments(alpha)%matrix, &
mo_coeff=dcdr_env%dCR_prime(1)%matrix, work=tmp_fm%matrix, nmo=nmo, &
icenter=icenter, &
res=tmp_fm_like_mos(alpha)%matrix)
CALL dbcsr_set(dcdr_env%moments(alpha)%matrix, 0.0_dp)
END DO
CALL cp_fm_gemm("T", "N", nmo, nmo, nao, &
1.0_dp, mo_coeff, tmp_fm_like_mos(alpha)%matrix, &
0.0_dp, tmp_fm_momo(alpha)%matrix)
CALL cp_fm_get_diag(tmp_fm_momo(alpha)%matrix, diagonal_elements)
DO icenter = 1, dcdr_env%nbr_center(1)
map_atom = mapping_wannier_atom(icenter)
map_molecule = mapping_atom_molecule(map_atom)
tmp_aptcontr = this_factor*diagonal_elements(icenter)
apt_subset(dcdr_env%beta, alpha, dcdr_env%lambda, map_molecule) &
= apt_subset(dcdr_env%beta, alpha, dcdr_env%lambda, map_molecule) + tmp_aptcontr
apt_center(dcdr_env%beta, alpha, dcdr_env%lambda, icenter) &
= apt_center(dcdr_env%beta, alpha, dcdr_env%lambda, icenter) + tmp_aptcontr
END DO
apt_coeff_derivative = this_factor*SUM(diagonal_elements)
apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) &
= apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) + apt_coeff_derivative
END DO
! SECOND CONTRIBUTION: We assemble all combinations of r_i, dphi/d(idir)
! build part with AOs differentiated with respect to nuclear coordinates
! difdip contains derivatives with respect to atom dcdr_env%lambda
! difdip(alpha, beta): < a | r_alpha | d b/dR_beta >
this_factor = -f_spin
DO alpha = 1, 3
DO icenter = 1, dcdr_env%nbr_center(1)
! Build the AO matrix with the right wannier center as reference point
CALL dbcsr_set(dcdr_env%matrix_difdip(1, dcdr_env%beta)%matrix, 0._dp)
CALL dbcsr_set(dcdr_env%matrix_difdip(2, dcdr_env%beta)%matrix, 0._dp)
CALL dbcsr_set(dcdr_env%matrix_difdip(3, dcdr_env%beta)%matrix, 0._dp)
CALL dipole_deriv_ao(qs_env, dcdr_env%matrix_difdip, dcdr_env%delta_basis_function, &
1, centers_set(1)%array(1:3, icenter))
CALL multiply_localization(ao_matrix=dcdr_env%matrix_difdip(alpha, dcdr_env%beta)%matrix, &
mo_coeff=mo_coeff, work=tmp_fm%matrix, nmo=nmo, &
icenter=icenter, &
res=tmp_fm_like_mos(alpha)%matrix)
END DO ! icenter
CALL cp_fm_gemm("T", "N", nmo, nmo, nao, &
1.0_dp, mo_coeff, tmp_fm_like_mos(alpha)%matrix, &
0.0_dp, tmp_fm_momo(alpha)%matrix)
CALL cp_fm_get_diag(tmp_fm_momo(alpha)%matrix, diagonal_elements)
DO icenter = 1, dcdr_env%nbr_center(1)
map_atom = mapping_wannier_atom(icenter)
map_molecule = mapping_atom_molecule(map_atom)
tmp_aptcontr = this_factor*diagonal_elements(icenter)
apt_subset(dcdr_env%beta, alpha, dcdr_env%lambda, map_molecule) &
= apt_subset(dcdr_env%beta, alpha, dcdr_env%lambda, map_molecule) + tmp_aptcontr
apt_center(dcdr_env%beta, alpha, dcdr_env%lambda, icenter) &
= apt_center(dcdr_env%beta, alpha, dcdr_env%lambda, icenter) + tmp_aptcontr
END DO
! The negative sign is due to dipole_deriv_ao computing the derivatives with respect to nuclear coordinates.
apt_basis_derivative = this_factor*SUM(diagonal_elements)
apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) &
= apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) + apt_basis_derivative
END DO ! alpha
! Finally the nuclear contribution: nuclear charge * Kronecker_delta_{dcdr_env%beta,i}
CALL get_atomic_kind(particle_set(dcdr_env%lambda)%atomic_kind, kind_number=ikind)
CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge, ghost=ghost)
IF (.NOT. ghost) THEN ! Those come from the pseudopotential, right?
apt_nuc(dcdr_env%beta, dcdr_env%beta, dcdr_env%lambda) = &
apt_nuc(dcdr_env%beta, dcdr_env%beta, dcdr_env%lambda) + charge
map_molecule = mapping_atom_molecule(dcdr_env%lambda)
apt_subset(dcdr_env%beta, dcdr_env%beta, dcdr_env%lambda, map_molecule) &
= apt_subset(dcdr_env%beta, dcdr_env%beta, dcdr_env%lambda, map_molecule) + charge
END IF
! And deallocate all the things!
CALL cp_fm_release(tmp_fm%matrix)
DO i = 1, 3
CALL cp_fm_release(tmp_fm_like_mos(i)%matrix)
CALL cp_fm_release(tmp_fm_momo(i)%matrix)
END DO
DEALLOCATE (tmp_fm, tmp_fm_like_mos, tmp_fm_momo)
DEALLOCATE (diagonal_elements)
CALL timestop(handle)
END SUBROUTINE apt_dR_localization
END MODULE qs_dcdr

789
src/qs_dcdr_ao.F Normal file
View file

@ -0,0 +1,789 @@
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright 2000-2021 CP2K developers group <https://cp2k.org> !
! !
! SPDX-License-Identifier: GPL-2.0-or-later !
!--------------------------------------------------------------------------------------------------!
! **************************************************************************************************
!> \brief Calculate the derivatives of the MO coefficients wrt nuclear coordinates
!> \author Sandra Luber, Edward Ditler
! **************************************************************************************************
MODULE qs_dcdr_ao
USE atomic_kind_types, ONLY: atomic_kind_type
USE basis_set_types, ONLY: gto_basis_set_p_type,&
gto_basis_set_type
USE core_ppl, ONLY: build_core_ppl
USE core_ppnl, ONLY: build_core_ppnl
USE cp_control_types, ONLY: dft_control_type
USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
copy_fm_to_dbcsr
USE cp_fm_types, ONLY: cp_fm_create,&
cp_fm_release,&
cp_fm_type
USE cp_gemm_interface, ONLY: cp_gemm
USE cp_log_handling, ONLY: cp_get_default_logger,&
cp_logger_type
USE dbcsr_api, ONLY: dbcsr_copy,&
dbcsr_get_block_p,&
dbcsr_p_type,&
dbcsr_set,&
dbcsr_type
USE input_constants, ONLY: do_ppl_analytic
USE input_section_types, ONLY: section_vals_get_subs_vals,&
section_vals_type
USE kinds, ONLY: default_string_length,&
dp
USE orbital_pointers, ONLY: ncoset
USE particle_types, ONLY: particle_type
USE pw_env_types, ONLY: pw_env_get,&
pw_env_type
USE pw_methods, ONLY: pw_axpy,&
pw_copy,&
pw_scale,&
pw_transfer,&
pw_zero
USE pw_poisson_methods, ONLY: pw_poisson_solve
USE pw_poisson_types, ONLY: pw_poisson_type
USE pw_pool_types, ONLY: pw_pool_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_collocate_density, ONLY: calculate_drho_core,&
calculate_drho_elec_dR
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_integral_utils, ONLY: basis_set_list_setup,&
get_memory_usage
USE qs_integrate_potential, ONLY: integrate_v_dbasis,&
integrate_v_rspace
USE qs_kind_types, ONLY: qs_kind_type
USE qs_ks_types, ONLY: get_ks_env,&
qs_ks_env_type
USE qs_linres_types, ONLY: dcdr_env_type
USE qs_neighbor_list_types, ONLY: get_iterator_info,&
get_neighbor_list_set_p,&
neighbor_list_iterate,&
neighbor_list_iterator_create,&
neighbor_list_iterator_p_type,&
neighbor_list_iterator_release,&
neighbor_list_set_p_type
USE qs_rho_methods, ONLY: qs_rho_rebuild,&
qs_rho_update_rho
USE qs_rho_types, ONLY: qs_rho_create,&
qs_rho_get,&
qs_rho_release,&
qs_rho_type
USE qs_vxc, ONLY: qs_vxc_create
USE virial_types, ONLY: virial_type
USE xc, ONLY: xc_calc_2nd_deriv,&
xc_prep_2nd_deriv
USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
xc_dset_release
USE xc_rho_set_types, ONLY: xc_rho_set_release,&
xc_rho_set_type
!$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_get_num_threads
!$ USE OMP_LIB, ONLY: omp_lock_kind, &
!$ omp_init_lock, omp_set_lock, &
!$ omp_unset_lock, omp_destroy_lock
#include "./base/base_uses.f90"
IMPLICIT NONE
PRIVATE
PUBLIC :: core_dR, d_vhxc_dR, d_core_charge_density_dR, apply_op_constant_term
PUBLIC :: vhxc_R_perturbed_basis_functions
PUBLIC :: hr_mult_by_delta_1d
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_dcdr_ao'
CONTAINS
! **************************************************************************************************
!> \brief Build the perturbed density matrix correction depending on the overlap derivative
!> \param qs_env ...
!> \param dcdr_env ...
!> \param overlap1 Overlap derivative in AO basis
!> \author Edward Ditler
! **************************************************************************************************
SUBROUTINE apply_op_constant_term(qs_env, dcdr_env, overlap1)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(dcdr_env_type) :: dcdr_env
TYPE(dbcsr_p_type), OPTIONAL :: overlap1
CHARACTER(len=*), PARAMETER :: routineN = 'apply_op_constant_term', &
routineP = moduleN//':'//routineN
INTEGER :: handle
REAL(KIND=dp) :: energy_hartree
TYPE(cp_fm_type), POINTER :: rho_ao_fm, rho_ao_s1, rho_ao_s1_rho_ao, &
s1_ao
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho1_ao, rho_ao
TYPE(pw_env_type), POINTER :: pw_env
TYPE(pw_p_type) :: rho1_tot_gspace, v_hartree_gspace, &
v_hartree_rspace
TYPE(pw_p_type), DIMENSION(:), POINTER :: rho1_g, rho1_g_pw, rho1_r, rho_r, &
v_rspace_new, v_xc
TYPE(pw_poisson_type), POINTER :: poisson_env
TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
TYPE(qs_rho_type), POINTER :: perturbed_density, rho
TYPE(section_vals_type), POINTER :: input, xc_section
TYPE(xc_derivative_set_type), POINTER :: deriv_set
TYPE(xc_rho_set_type), POINTER :: rho1_set, rho_set
! Build the perturbed density matrix correction depending on the overlap derivative
! P1 = C0 C1 + C1 C0
! - C0_(mu j) S1_(jk) C0_(k nu)
! This routine is adapted from apply_op_2_dft. There, build_dm_response builds
! C0 * dCR + dCR * C0.
! build_dm_response is computing $-1 * (C^0 C^1 + C^1 C^0)$ and later on in the
! integration the factor 2 is applied to account for the occupancy.
! The sign is negative because the kernel is on the RHS of the Sternheimer equation.
!
! The correction factor in this routine needs to have
! the opposite sign mathematically as (C0 C1 + C1 C0)
! so the same sign in the code because of the $-1$ in dCR
! so the opposite sign in the code because we are on the LHS of the Sternheimer equation.
!
! This term must not go into the kernel applied by the linear response solver, because
! for the (P)CG algorithm, all constant terms have to be on one side of the equations
! and all solution dependent terms must be on the other side.
CALL timeset(routineN, handle)
NULLIFY (auxbas_pw_pool, pw_env, v_rspace_new, rho1_r, rho1_g_pw, &
v_xc, rho1_set, poisson_env, input, rho, rho1_g)
CALL dbcsr_set(dcdr_env%perturbed_dm_correction, 0._dp)
CALL cp_fm_create(rho_ao_fm, dcdr_env%aoao_fm_struct)
CALL cp_fm_create(rho_ao_s1, dcdr_env%aoao_fm_struct)
CALL cp_fm_create(rho_ao_s1_rho_ao, dcdr_env%aoao_fm_struct)
CALL cp_fm_create(s1_ao, dcdr_env%aoao_fm_struct)
IF (PRESENT(overlap1)) THEN
CALL copy_dbcsr_to_fm(overlap1%matrix, s1_ao)
ELSE
CALL copy_dbcsr_to_fm(dcdr_env%matrix_s1(dcdr_env%beta + 1)%matrix, s1_ao)
END IF
CALL cp_gemm('N', 'T', dcdr_env%nao, dcdr_env%nao, dcdr_env%nmo, &
1.0_dp, dcdr_env%mo_coeff(1)%matrix, dcdr_env%mo_coeff(1)%matrix, &
0.0_dp, rho_ao_fm)
CALL cp_gemm('N', 'N', dcdr_env%nao, dcdr_env%nao, dcdr_env%nao, &
1.0_dp, rho_ao_fm, s1_ao, &
0.0_dp, rho_ao_s1)
CALL cp_gemm('N', 'N', dcdr_env%nao, dcdr_env%nao, dcdr_env%nao, &
-1._dp, rho_ao_s1, rho_ao_fm, & ! this is the sign mentioned above.
0.0_dp, rho_ao_s1_rho_ao)
CALL copy_fm_to_dbcsr(rho_ao_s1_rho_ao, dcdr_env%perturbed_dm_correction)
CALL cp_fm_release(rho_ao_fm)
CALL cp_fm_release(rho_ao_s1)
CALL cp_fm_release(rho_ao_s1_rho_ao)
CALL cp_fm_release(s1_ao)
! Done building the density matrix correction
CALL dbcsr_set(dcdr_env%matrix_apply_op_constant(1)%matrix, 0.0_dp)
! Build the density struct from the environment
NULLIFY (perturbed_density)
CALL qs_rho_create(perturbed_density)
CALL qs_rho_rebuild(perturbed_density, qs_env=qs_env)
! ... set the density matrix to be the perturbed density matrix
CALL qs_rho_get(perturbed_density, rho_ao=rho1_ao)
CALL dbcsr_copy(rho1_ao(1)%matrix, dcdr_env%perturbed_dm_correction)
! ... updates rho_r and rho_g to the rho%rho_ao.
CALL qs_rho_update_rho(rho_struct=perturbed_density, &
qs_env=qs_env)
! Also update the qs_env%rho
CALL get_qs_env(qs_env, rho=rho)
CALL qs_rho_update_rho(rho, qs_env=qs_env)
CALL qs_rho_get(rho, rho_ao=rho_ao, rho_r=rho_r)
energy_hartree = 0.0_dp
CALL get_qs_env(qs_env=qs_env, &
pw_env=pw_env, &
input=input)
! Create the temporary grids
CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
poisson_env=poisson_env)
! Allocate deriv_set and rho_set
NULLIFY (deriv_set, rho_set)
xc_section => section_vals_get_subs_vals(input, "DFT%XC")
CALL xc_prep_2nd_deriv(deriv_set, rho_set, &
rho_r, auxbas_pw_pool, &
xc_section=xc_section)
! Done with deriv_set and rho_set
ALLOCATE (v_rspace_new(1))
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, v_hartree_rspace%pw, &
use_data=REALDATA3D, &
in_space=REALSPACE)
! Calculate the Hartree potential on the total density
CALL pw_pool_create_pw(auxbas_pw_pool, rho1_tot_gspace%pw, &
use_data=COMPLEXDATA1D, &
in_space=RECIPROCALSPACE)
CALL qs_rho_get(perturbed_density, rho_g=rho1_g, rho_r=rho1_r)
CALL pw_copy(rho1_g(1)%pw, rho1_tot_gspace%pw)
CALL pw_poisson_solve(poisson_env, rho1_tot_gspace%pw, &
energy_hartree, &
v_hartree_gspace%pw)
CALL pw_transfer(v_hartree_gspace%pw, v_hartree_rspace%pw)
CALL pw_pool_give_back_pw(auxbas_pw_pool, rho1_tot_gspace%pw)
! Calculate the second derivative of the exchange-correlation potential
CALL xc_calc_2nd_deriv(v_xc, deriv_set, rho_set, &
rho1_r, rho1_g_pw, auxbas_pw_pool, xc_section, gapw=.FALSE.)
v_rspace_new(1)%pw => v_xc(1)%pw
DEALLOCATE (v_xc)
CALL xc_rho_set_release(rho1_set)
! Done calculating the potentials
!-------------------------------!
! Add both hartree and xc terms !
!-------------------------------!
CALL pw_scale(v_rspace_new(1)%pw, 2._dp*v_rspace_new(1)%pw%pw_grid%dvol)
CALL pw_scale(v_hartree_rspace%pw, 2._dp*v_hartree_rspace%pw%pw_grid%dvol)
! pw2 = 1.*pw1 + pw2
CALL pw_axpy(v_hartree_rspace%pw, v_rspace_new(1)%pw, 1._dp)
CALL dbcsr_set(dcdr_env%matrix_apply_op_constant(1)%matrix, 0.0_dp)
CALL integrate_v_rspace(v_rspace=v_rspace_new(1), &
hmat=dcdr_env%matrix_apply_op_constant(1), &
qs_env=qs_env, &
calculate_forces=.FALSE.)
CALL pw_pool_give_back_pw(auxbas_pw_pool, v_hartree_gspace%pw)
CALL pw_pool_give_back_pw(auxbas_pw_pool, v_hartree_rspace%pw)
CALL pw_pool_give_back_pw(auxbas_pw_pool, v_rspace_new(1)%pw)
DEALLOCATE (v_rspace_new)
CALL qs_rho_release(perturbed_density)
CALL xc_rho_set_release(rho_set, auxbas_pw_pool)
CALL xc_dset_release(deriv_set)
CALL timestop(handle)
END SUBROUTINE apply_op_constant_term
! **************************************************************************************************
!> \brief Calculate the derivative of the Hartree term due to the core charge density
!> \param qs_env ...
!> \param dcdr_env ...
!> \author Edward Ditler
! **************************************************************************************************
SUBROUTINE d_core_charge_density_dR(qs_env, dcdr_env)
! drho_core contribution
! sum over all directions
! output in ao x ao
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(dcdr_env_type) :: dcdr_env
CHARACTER(len=*), PARAMETER :: routineN = 'd_core_charge_density_dR', &
routineP = moduleN//':'//routineN
INTEGER :: beta, handle
TYPE(cp_logger_type), POINTER :: logger
TYPE(dft_control_type), POINTER :: dft_control
TYPE(pw_env_type), POINTER :: pw_env
TYPE(pw_p_type) :: drho_g, v_hartree_gspace, &
v_hartree_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_rho_type), POINTER :: rho
CALL timeset(routineN, handle)
logger => cp_get_default_logger()
NULLIFY (pw_env, auxbas_pw_pool, pw_pools, poisson_env, dft_control, &
v_hartree_gspace%pw, v_hartree_rspace%pw, rho)
CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, rho=rho, &
dft_control=dft_control)
CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, poisson_env=poisson_env, &
pw_pools=pw_pools)
! Create the Hartree potential grids in real and reciprocal space.
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, &
v_hartree_rspace%pw, &
use_data=REALDATA3D, &
in_space=REALSPACE)
! Create the grid for the derivative of the core potential
CALL pw_pool_create_pw(auxbas_pw_pool, drho_g%pw, &
use_data=COMPLEXDATA1D, in_space=RECIPROCALSPACE)
DO beta = 1, 3
CALL pw_zero(v_hartree_gspace%pw)
CALL pw_zero(v_hartree_rspace%pw)
CALL pw_zero(drho_g%pw)
! Calculate the Hartree potential on the perturbed density and Poisson solve it
CALL calculate_drho_core(drho_core=drho_g, qs_env=qs_env, &
beta=beta, lambda=dcdr_env%lambda)
CALL pw_poisson_solve(poisson_env, drho_g%pw, &
vhartree=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)
! Calculate the integrals
CALL integrate_v_rspace(v_rspace=v_hartree_rspace, &
hmat=dcdr_env%matrix_core_charge_1(beta), &
qs_env=qs_env, &
calculate_forces=.FALSE.)
END DO
CALL pw_pool_give_back_pw(auxbas_pw_pool, drho_g%pw)
CALL pw_pool_give_back_pw(auxbas_pw_pool, v_hartree_rspace%pw)
CALL pw_pool_give_back_pw(auxbas_pw_pool, v_hartree_gspace%pw)
CALL timestop(handle)
END SUBROUTINE d_core_charge_density_dR
! **************************************************************************************************
!> \brief Core Hamiltonian contributions to the operator (the pseudopotentials)
!> \param qs_env ...
!> \param dcdr_env ..
!> \author Edward Ditler
! **************************************************************************************************
SUBROUTINE core_dR(qs_env, dcdr_env)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(dcdr_env_type) :: dcdr_env
CHARACTER(LEN=*), PARAMETER :: routineN = 'core_dR', routineP = moduleN//':'//routineN
CHARACTER(LEN=default_string_length) :: my_basis_type
INTEGER :: handle, nder
INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
LOGICAL :: calculate_forces, failure, ppl_present, &
ppnl_present, use_virial
REAL(KIND=dp) :: eps_ppnl
REAL(KIND=dp), DIMENSION(:, :), POINTER :: deltaR
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_hc_pass, matrix_p_pass, &
matrix_ppnl_1_pass
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(qs_rho_type), POINTER :: rho
TYPE(virial_type), POINTER :: virial
CALL timeset(routineN, handle)
failure = .FALSE.
NULLIFY (atomic_kind_set, qs_kind_set, ks_env, dft_control, particle_set, sab_orb, sac_ppl, sap_ppnl, virial, rho, rho_ao)
CALL get_qs_env(qs_env=qs_env, &
atomic_kind_set=atomic_kind_set, &
qs_kind_set=qs_kind_set, &
ks_env=ks_env, &
dft_control=dft_control, &
particle_set=particle_set, &
sab_orb=sab_orb, &
sac_ppl=sac_ppl, &
sap_ppnl=sap_ppnl, &
virial=virial)
CALL get_ks_env(ks_env=ks_env, rho=rho)
CALL qs_rho_get(rho, rho_ao=rho_ao)
deltaR => dcdr_env%delta_basis_function
nder = 1
calculate_forces = .FALSE.
my_basis_type = "ORB"
! *** compute the ppl contribution to the core hamiltonian ***
ppl_present = ASSOCIATED(sac_ppl)
IF (ppl_present) THEN
IF (dft_control%qs_control%do_ppl_method == do_ppl_analytic) THEN
matrix_hc_pass(1:3, 1:1) => dcdr_env%matrix_hc(1:3)
matrix_p_pass(1:1, 1:1) => rho_ao(1:1)
CALL build_core_ppl(matrix_h=matrix_hc_pass, matrix_p=matrix_p_pass, &
force=force, virial=virial, calculate_forces=calculate_forces, &
use_virial=use_virial, nder=nder, qs_kind_set=qs_kind_set, &
atomic_kind_set=atomic_kind_set, particle_set=particle_set, &
sab_orb=sab_orb, sac_ppl=sac_ppl, basis_type=my_basis_type, &
nimages=1, cell_to_index=cell_to_index, deltaR=deltaR)
END IF ! ppl_analytic
END IF ! ppl_present
! *** compute the ppnl contribution to the core hamiltonian ***
eps_ppnl = dft_control%qs_control%eps_ppnl
ppnl_present = ASSOCIATED(sap_ppnl)
IF (ppnl_present) THEN
matrix_ppnl_1_pass(1:3, 1:1) => dcdr_env%matrix_ppnl_1(1:3)
CALL build_core_ppnl(matrix_h=matrix_ppnl_1_pass, matrix_p=matrix_p_pass, force=force, virial=virial, &
calculate_forces=calculate_forces, use_virial=use_virial, nder=nder, &
qs_kind_set=qs_kind_set, atomic_kind_set=atomic_kind_set, &
particle_set=particle_set, sab_orb=sab_orb, sap_ppnl=sap_ppnl, &
eps_ppnl=eps_ppnl, nimages=1, cell_to_index=cell_to_index, &
basis_type=my_basis_type, deltaR=deltaR)
END IF
CALL timestop(handle)
END SUBROUTINE core_dR
! **************************************************************************************************
!> \brief The derivatives of the basis functions going into the HXC potential wrt nuclear positions
!> \param qs_env ...
!> \param dcdr_env ...
!> \author Edward Ditler
! **************************************************************************************************
SUBROUTINE d_vhxc_dR(qs_env, dcdr_env)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(dcdr_env_type) :: dcdr_env
CHARACTER(len=*), PARAMETER :: routineN = 'd_vhxc_dR', routineP = moduleN//':'//routineN
INTEGER :: handle, idir
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao
TYPE(pw_env_type), POINTER :: pw_env
TYPE(pw_p_type) :: v_hartree_gspace, v_hartree_rspace
TYPE(pw_p_type), DIMENSION(:), POINTER :: drho_g, drho_r, rho_r, v_xc
TYPE(pw_poisson_type), POINTER :: poisson_env
TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
TYPE(qs_rho_type), POINTER :: rho
TYPE(section_vals_type), POINTER :: input, xc_section
TYPE(xc_derivative_set_type), POINTER :: my_deriv_set
TYPE(xc_rho_set_type), POINTER :: my_rho_set
CALL timeset(routineN, handle)
CALL get_qs_env(qs_env=qs_env, &
pw_env=pw_env, &
input=input, &
rho=rho)
CALL qs_rho_get(rho, rho_ao=rho_ao, rho_r=rho_r)
xc_section => section_vals_get_subs_vals(input, "DFT%XC")
! get the tmp grids
ALLOCATE (drho_r(1))
ALLOCATE (drho_g(1))
CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
pw_pools=pw_pools, poisson_env=poisson_env)
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, v_hartree_rspace%pw, &
use_data=REALDATA3D, &
in_space=REALSPACE)
CALL pw_pool_create_pw(auxbas_pw_pool, drho_r(1)%pw, &
use_data=REALDATA3D, in_space=REALSPACE)
CALL pw_pool_create_pw(auxbas_pw_pool, drho_g(1)%pw, &
use_data=COMPLEXDATA1D, in_space=RECIPROCALSPACE)
DO idir = 1, 3
NULLIFY (my_deriv_set, my_rho_set)
CALL pw_zero(v_hartree_gspace%pw)
CALL pw_zero(v_hartree_rspace%pw)
CALL pw_zero(drho_r(1)%pw)
CALL pw_zero(drho_g(1)%pw)
! Get the density
CALL calculate_drho_elec_dR(matrix_p=rho_ao(1)%matrix, &
drho=drho_r(1), &
drho_gspace=drho_g(1), &
qs_env=qs_env, &
beta=idir, lambda=dcdr_env%lambda)
! Get the Hartree potential corresponding to the perturbed density
CALL pw_poisson_solve(poisson_env, drho_g(1)%pw, &
vhartree=v_hartree_gspace%pw)
CALL pw_transfer(v_hartree_gspace%pw, v_hartree_rspace%pw)
! Get the XC potential corresponding to the perturbed density
CALL xc_prep_2nd_deriv(my_deriv_set, my_rho_set, &
rho_r, auxbas_pw_pool, &
xc_section=xc_section)
CALL xc_calc_2nd_deriv(v_xc, my_deriv_set, my_rho_set, &
drho_r, drho_g, auxbas_pw_pool, xc_section, gapw=.FALSE.)
CALL xc_dset_release(my_deriv_set)
CALL xc_rho_set_release(my_rho_set)
!-------------------------------!
! Add both hartree and xc terms !
!-------------------------------!
! Can the dvol be different?
CALL pw_scale(v_xc(1)%pw, v_xc(1)%pw%pw_grid%dvol)
CALL pw_axpy(v_hartree_rspace%pw, v_xc(1)%pw, v_hartree_rspace%pw%pw_grid%dvol)
CALL integrate_v_rspace(v_rspace=v_xc(1), &
hmat=dcdr_env%matrix_d_vhxc_dR(idir), &
qs_env=qs_env, &
calculate_forces=.FALSE.)
! v_xc gets allocated again in xc_calc_2nd_deriv
CALL pw_pool_give_back_pw(auxbas_pw_pool, v_xc(1)%pw)
DEALLOCATE (v_xc)
END DO ! idir
CALL pw_pool_give_back_pw(auxbas_pw_pool, v_hartree_gspace%pw)
CALL pw_pool_give_back_pw(auxbas_pw_pool, v_hartree_rspace%pw)
CALL pw_pool_give_back_pw(auxbas_pw_pool, drho_g(1)%pw)
CALL pw_pool_give_back_pw(auxbas_pw_pool, drho_r(1)%pw)
DEALLOCATE (drho_g)
DEALLOCATE (drho_r)
NULLIFY (my_deriv_set)
CALL timestop(handle)
END SUBROUTINE d_vhxc_dR
! **************************************************************************************************
!> \brief The derivatives of the basis functions over which the HXC potential is integrated,
!> so < da/dR | Vhxc | b >
!> \param qs_env ...
!> \param dcdr_env ...
!> \author Edward Ditler
! **************************************************************************************************
SUBROUTINE vhxc_R_perturbed_basis_functions(qs_env, dcdr_env)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(dcdr_env_type) :: dcdr_env
CHARACTER(LEN=*), PARAMETER :: routineN = 'vhxc_R_perturbed_basis_functions', &
routineP = moduleN//':'//routineN
INTEGER :: handle
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_vhxc_dbasis
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p
TYPE(pw_env_type), POINTER :: pw_env
TYPE(pw_p_type) :: v_hartree_r
TYPE(pw_p_type), DIMENSION(:), POINTER :: v_hxc_r, v_tau_rspace
TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
TYPE(qs_energy_type), POINTER :: energy
TYPE(qs_ks_env_type), POINTER :: ks_env
TYPE(qs_rho_type), POINTER :: rho_struct
TYPE(section_vals_type), POINTER :: input, xc_section
CALL timeset(routineN, handle)
NULLIFY (rho_struct, energy, input, ks_env, pw_env, matrix_p)
CALL get_qs_env(qs_env, &
rho=rho_struct, &
energy=energy, &
input=input, &
ks_env=ks_env, &
pw_env=pw_env, &
v_hartree_rspace=v_hartree_r%pw)
CALL qs_rho_get(rho_struct, rho_ao_kp=matrix_p)
xc_section => section_vals_get_subs_vals(input, "DFT%XC")
NULLIFY (auxbas_pw_pool)
CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
! *** calculate the xc potential on the pw density ***
! *** associates v_hxc_r if the xc potential needs to be computed.
! If we do wavefunction fitting, we need the vxc_potential in the auxiliary basis set
NULLIFY (v_hxc_r, v_tau_rspace)
CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho_struct, xc_section=xc_section, &
vxc_rho=v_hxc_r, vxc_tau=v_tau_rspace, exc=energy%exc)
CALL pw_scale(v_hxc_r(1)%pw, v_hxc_r(1)%pw%pw_grid%dvol)
! sum up potentials and integrate
CALL pw_axpy(v_hartree_r%pw, v_hxc_r(1)%pw, 1._dp)
matrix_vhxc_dbasis => dcdr_env%matrix_vhxc_perturbed_basis(1, :)
CALL integrate_v_dbasis(v_rspace=v_hxc_r(1), &
matrix_p=matrix_p(1, 1)%matrix, &
matrix_vhxc_dbasis=matrix_vhxc_dbasis, &
qs_env=qs_env, &
lambda=dcdr_env%lambda)
CALL pw_pool_give_back_pw(auxbas_pw_pool, v_hxc_r(1)%pw)
DEALLOCATE (v_hxc_r)
CALL timestop(handle)
END SUBROUTINE vhxc_R_perturbed_basis_functions
! **************************************************************************************************
!> \brief Enforce that one of the basis functions in < a | O | b > is centered on atom lambda.
!> \param matrix ...
!> \param qs_kind_set ...
!> \param basis_type ...
!> \param sab_nl ...
!> \param lambda Atom index
!> \param direction_Or True: < a | O | b==lambda >, False: < a==lambda | O | b >
! **************************************************************************************************
SUBROUTINE hr_mult_by_delta_1d(matrix, qs_kind_set, basis_type, sab_nl, lambda, direction_Or)
TYPE(dbcsr_type), POINTER :: matrix
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
CHARACTER(LEN=*), INTENT(IN) :: basis_type
TYPE(neighbor_list_set_p_type), DIMENSION(:), &
POINTER :: sab_nl
INTEGER, INTENT(IN) :: lambda
LOGICAL, INTENT(IN) :: direction_Or
CHARACTER(len=*), PARAMETER :: routineN = 'hr_mult_by_delta_1d', &
routineP = moduleN//':'//routineN
INTEGER :: handle, iatom, icol, ikind, irow, jatom, &
jkind, ldsab, mepos, nkind, nseta, &
nsetb, nthread
INTEGER, DIMENSION(3) :: cell
INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
npgfb, nsgfa, nsgfb
INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
LOGICAL :: do_symmetric, found
REAL(KIND=dp), DIMENSION(3) :: rab
REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
REAL(KIND=dp), DIMENSION(:, :), POINTER :: k_block, rpgfa, rpgfb, scon_a, scon_b, &
zeta, zetb
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
CALL timeset(routineN, handle)
nkind = SIZE(qs_kind_set)
! check for symmetry
CPASSERT(SIZE(sab_nl) > 0)
CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
! prepare basis set
ALLOCATE (basis_set_list(nkind))
CALL basis_set_list_setup(basis_set_list, basis_type, qs_kind_set)
! *** Allocate work storage ***
ldsab = get_memory_usage(qs_kind_set, basis_type)
nthread = 1
!$ nthread = omp_get_max_threads()
! Iterate of neighbor list
CALL neighbor_list_iterator_create(nl_iterator, sab_nl, nthread=nthread)
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP SHARED (nthread,ldsab,nl_iterator, do_symmetric) &
!$OMP SHARED (ncoset,matrix,basis_set_list) &
!$OMP SHARED (direction_or, lambda) &
!$OMP PRIVATE (k_block,mepos,ikind,jkind,iatom,jatom,rab,cell) &
!$OMP PRIVATE (basis_set_a,basis_set_b) &
!$OMP PRIVATE (first_sgfa, la_max, la_min, npgfa, nsgfa, nseta, rpgfa, set_radius_a) &
!$OMP PRIVATE (zeta, first_sgfb, lb_max, lb_min, npgfb, nsetb, rpgfb, set_radius_b, nsgfb) &
!$OMP PRIVATE (zetb, scon_a, scon_b, irow, icol, found)
mepos = 0
!$ mepos = omp_get_thread_num()
DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0)
CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, &
iatom=iatom, jatom=jatom, r=rab, cell=cell)
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
scon_a => basis_set_a%scon
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
scon_b => basis_set_b%scon
zetb => basis_set_b%zet
IF (do_symmetric) THEN
IF (iatom <= jatom) THEN
irow = iatom
icol = jatom
ELSE
irow = jatom
icol = iatom
END IF
ELSE
irow = iatom
icol = jatom
END IF
NULLIFY (k_block)
CALL dbcsr_get_block_p(matrix, irow, icol, k_block, found)
CPASSERT(found)
IF (direction_Or) THEN
IF (jatom /= lambda) k_block(:, :) = 0._dp
ELSE IF (.NOT. direction_Or) THEN
IF (iatom /= lambda) k_block(:, :) = 0._dp
END IF
END DO
!$OMP END PARALLEL
CALL neighbor_list_iterator_release(nl_iterator)
! Release work storage
DEALLOCATE (basis_set_list)
CALL timestop(handle)
END SUBROUTINE hr_mult_by_delta_1d
END MODULE qs_dcdr_ao

1077
src/qs_dcdr_utils.F Normal file
View file

File diff suppressed because it is too large Load diff

6
src/qs_integrate_potential.F
View file

@ -21,7 +21,8 @@
! **************************************************************************************************
MODULE qs_integrate_potential
USE grid_api, ONLY: integrate_pgf_product
USE qs_integrate_potential_product, ONLY: integrate_v_rspace
USE qs_integrate_potential_product, ONLY: integrate_v_dbasis,&
integrate_v_rspace
USE qs_integrate_potential_single, ONLY: integrate_ppl_rspace,&
integrate_rho_nlcc,&
integrate_v_core_rspace,&
@ -38,7 +39,8 @@ MODULE qs_integrate_potential
! *** Public subroutines ***
! included from qs_integrate_potential_product
PUBLIC :: integrate_v_rspace
PUBLIC :: integrate_v_rspace, &
integrate_v_dbasis
! included from qs_integrate_potential_single
PUBLIC :: integrate_v_rspace_one_center, &

492
src/qs_integrate_potential_product.F
View file

@ -30,17 +30,29 @@
! **************************************************************************************************
MODULE qs_integrate_potential_product
USE admm_types, ONLY: admm_type
USE ao_util, ONLY: exp_radius_very_extended
USE atomic_kind_types, ONLY: atomic_kind_type,&
get_atomic_kind_set
USE cell_types, ONLY: cell_type
USE basis_set_types, ONLY: get_gto_basis_set,&
gto_basis_set_type
USE block_p_types, ONLY: block_p_type
USE cell_types, ONLY: cell_type,&
pbc
USE cp_control_types, ONLY: dft_control_type
USE cp_dbcsr_operations, ONLY: dbcsr_deallocate_matrix_set
USE cube_utils, ONLY: cube_info_type
USE dbcsr_api, ONLY: dbcsr_p_type
USE dbcsr_api, ONLY: dbcsr_copy,&
dbcsr_finalize,&
dbcsr_get_block_p,&
dbcsr_p_type,&
dbcsr_type
USE gaussian_gridlevels, ONLY: gridlevel_info_type
USE grid_api, ONLY: grid_integrate_task_list
USE grid_api, ONLY: grid_integrate_task_list,&
integrate_pgf_product
USE input_constants, ONLY: do_admm_exch_scaling_merlot
USE kinds, ONLY: default_string_length,&
dp
USE orbital_pointers, ONLY: ncoset
USE particle_types, ONLY: particle_type
USE pw_env_types, ONLY: pw_env_get,&
pw_env_type
@ -48,21 +60,28 @@ MODULE qs_integrate_potential_product
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_set,&
USE qs_kind_types, ONLY: get_qs_kind,&
get_qs_kind_set,&
qs_kind_type
USE realspace_grid_types, ONLY: realspace_grid_p_type,&
USE realspace_grid_types, ONLY: realspace_grid_desc_p_type,&
realspace_grid_p_type,&
rs_grid_release,&
rs_grid_retain
USE rs_pw_interface, ONLY: potential_pw2rs
USE task_list_methods, ONLY: rs_copy_to_buffer,&
rs_copy_to_matrices,&
rs_distribute_matrix,&
rs_gather_matrices,&
rs_scatter_matrices
USE task_list_types, ONLY: task_list_type
USE task_list_types, ONLY: atom_pair_type,&
task_list_type,&
task_type
USE virial_types, ONLY: virial_type
!$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_get_num_threads
!$ USE OMP_LIB, ONLY: omp_lock_kind, &
!$ omp_init_lock, omp_set_lock, &
!$ omp_unset_lock, omp_destroy_lock
#include "./base/base_uses.f90"
IMPLICIT NONE
@ -76,9 +95,468 @@ MODULE qs_integrate_potential_product
! *** qs_integrate_potential
PUBLIC :: integrate_v_rspace
PUBLIC :: integrate_v_dbasis
CONTAINS
! **************************************************************************************************
!> \brief Integrate a potential v_rspace over the derivatives of the basis functions
!> < da/dR | V | b > + < a | V | db/dR >
!> Adapted from the old version of integrate_v_rspace (ED)
!> \param v_rspace ...
!> \param matrix_vhxc_dbasis ...
!> \param matrix_p ...
!> \param qs_env ...
!> \param lambda The atom index.
! **************************************************************************************************
SUBROUTINE integrate_v_dbasis(v_rspace, matrix_vhxc_dbasis, matrix_p, qs_env, lambda)
TYPE(pw_p_type) :: v_rspace
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_vhxc_dbasis
TYPE(dbcsr_type), POINTER :: matrix_p
TYPE(qs_environment_type), POINTER :: qs_env
INTEGER, INTENT(IN) :: lambda
CHARACTER(len=*), PARAMETER :: routineN = 'integrate_v_dbasis'
INTEGER :: bcol, brow, group, handle, i, iatom, igrid_level, ikind, ikind_old, ilevel, img, &
ipair, ipgf, ipgf_new, iset, iset_new, iset_old, itask, ithread, jatom, jkind, jkind_old, &
jpgf, jpgf_new, jset, jset_new, jset_old, maxco, maxpgf, maxset, maxsgf_set, na1, na2, &
natom, nb1, nb2, ncoa, ncob, nimages, nkind, nseta, nsetb, nthread, sgfa, sgfb
INTEGER, ALLOCATABLE, DIMENSION(:, :) :: block_touched
INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
npgfb, nsgfa, nsgfb
INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
LOGICAL :: atom_pair_changed, atom_pair_done, dh_duplicated, distributed_grids, found, &
my_compute_tau, my_gapw, new_set_pair_coming, pab_required, scatter, use_subpatch
REAL(KIND=dp) :: eps_rho_rspace, f, prefactor, radius, &
scalef, zetp
REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, ra, rab, rab_inv, rb, &
rp
REAL(KIND=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b
REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
REAL(KIND=dp), DIMENSION(:, :), POINTER :: h_block, hab, p_block, pab, rpgfa, &
rpgfb, sphi_a, sphi_b, work, zeta, zetb
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: habt, hadb, hdab, pabt, workt
REAL(kind=dp), DIMENSION(:, :, :, :), POINTER :: hadbt, hdabt
TYPE(atom_pair_type), DIMENSION(:), POINTER :: atom_pair_recv, atom_pair_send
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: vhxc_block
TYPE(cell_type), POINTER :: cell
TYPE(cube_info_type), DIMENSION(:), POINTER :: cube_info
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: deltap
TYPE(dft_control_type), POINTER :: dft_control
TYPE(gridlevel_info_type), POINTER :: gridlevel_info
TYPE(gto_basis_set_type), POINTER :: orb_basis_set
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(pw_env_type), POINTER :: pw_env
TYPE(qs_force_type), DIMENSION(:), POINTER :: force
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
POINTER :: rs_descs
TYPE(realspace_grid_p_type), DIMENSION(:), POINTER :: rs_rho
TYPE(task_list_type), POINTER :: task_list, task_list_soft
TYPE(task_type), DIMENSION(:), POINTER :: tasks
CALL timeset(routineN, handle)
NULLIFY (pw_env)
! get the task lists
my_gapw = .FALSE.
CALL get_qs_env(qs_env=qs_env, &
task_list=task_list, &
task_list_soft=task_list_soft)
CPASSERT(ASSOCIATED(task_list))
! the information on the grids is provided through pw_env
! pw_env has to be the parent env for the potential grid (input)
! there is an option to provide an external grid
CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
! *** assign from pw_env
gridlevel_info => pw_env%gridlevel_info
cube_info => pw_env%cube_info
! get all the general information on the system we are working on
CALL get_qs_env(qs_env=qs_env, &
atomic_kind_set=atomic_kind_set, &
qs_kind_set=qs_kind_set, &
cell=cell, &
natom=natom, &
dft_control=dft_control, &
particle_set=particle_set)
! *** set up the rs multi-grids
CPASSERT(ASSOCIATED(pw_env))
CALL pw_env_get(pw_env, rs_descs=rs_descs, rs_grids=rs_rho)
DO igrid_level = 1, gridlevel_info%ngrid_levels
CALL rs_grid_retain(rs_rho(igrid_level)%rs_grid)
distributed_grids = rs_rho(igrid_level)%rs_grid%desc%distributed
END DO
! get mpi group from rs_rho
group = rs_rho(1)%rs_grid%desc%group
! transform the potential on the rs_multigrids
CALL potential_pw2rs(rs_rho, v_rspace, pw_env)
nkind = SIZE(qs_kind_set)
CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
maxco=maxco, &
maxsgf_set=maxsgf_set, &
basis_type="ORB")
! short cuts to task list variables
tasks => task_list%tasks
atom_pair_send => task_list%atom_pair_send
atom_pair_recv => task_list%atom_pair_recv
! needs to be consistent with rho_rspace
eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
! *** Initialize working density matrix ***
! distributed rs grids require a matrix that will be changed
! whereas this is not the case for replicated grids
ALLOCATE (deltap(dft_control%nimages))
IF (distributed_grids) THEN
! this matrix has no strict sparsity pattern in parallel
! deltap%sparsity_id=-1
CALL dbcsr_copy(deltap(1)%matrix, matrix_p, name="DeltaP")
ELSE
deltap(1)%matrix => matrix_p
END IF
nthread = 1
!$ nthread = omp_get_max_threads()
! *** Allocate work storage ***
NULLIFY (pabt, habt, workt)
ALLOCATE (habt(maxco, maxco, 0:nthread))
ALLOCATE (workt(maxco, maxsgf_set, 0:nthread))
ALLOCATE (hdabt(3, maxco, maxco, 0:nthread))
ALLOCATE (hadbt(3, maxco, maxco, 0:nthread))
ALLOCATE (pabt(maxco, maxco, 0:nthread))
IF (distributed_grids) THEN
CALL rs_distribute_matrix(rs_descs, deltap, atom_pair_send, atom_pair_recv, &
nimages, scatter=.TRUE.)
END IF
!$OMP PARALLEL DEFAULT(NONE), &
!$OMP SHARED(workt,habt,hdabt,hadbt,pabt,tasks,particle_set,natom,maxset), &
!$OMP SHARED(maxpgf,my_gapw,matrix_vhxc_dbasis,deltap), &
!$OMP SHARED(pab_required,ncoset,rs_rho,cube_info,my_compute_tau), &
!$OMP SHARED(eps_rho_rspace,force,cell), &
!$OMP SHARED(gridlevel_info,task_list,block_touched,nthread,qs_kind_set), &
!$OMP SHARED(nimages,lambda, dh_duplicated), &
!$OMP PRIVATE(ithread,work,hab,hdab,hadb,pab,iset_old,jset_old), &
!$OMP PRIVATE(ikind_old,jkind_old,iatom,jatom,iset,jset,ikind,jkind,ilevel,ipgf,jpgf), &
!$OMP PRIVATE(img,brow,bcol,orb_basis_set,first_sgfa,la_max,la_min,npgfa,nseta,nsgfa), &
!$OMP PRIVATE(rpgfa,set_radius_a,sphi_a,zeta,first_sgfb,lb_max,lb_min,npgfb), &
!$OMP PRIVATE(nsetb,nsgfb,rpgfb,set_radius_b,sphi_b,zetb,found), &
!$OMP PRIVATE(force_a,force_b,my_virial_a,my_virial_b,atom_pair_changed,h_block, vhxc_block), &
!$OMP PRIVATE(p_block,ncoa,sgfa,ncob,sgfb,rab,ra,rb,rp,zetp,f,prefactor,radius,igrid_level), &
!$OMP PRIVATE(na1,na2,nb1,nb2,use_subpatch,rab_inv,new_set_pair_coming,atom_pair_done), &
!$OMP PRIVATE(iset_new,jset_new,ipgf_new,jpgf_new,scalef), &
!$OMP PRIVATE(itask)
IF (.NOT. ALLOCATED(vhxc_block)) ALLOCATE (vhxc_block(3))
ithread = 0
!$ ithread = omp_get_thread_num()
work => workt(:, :, ithread)
hab => habt(:, :, ithread)
pab => pabt(:, :, ithread)
hdab => hdabt(:, :, :, ithread)
hadb => hadbt(:, :, :, ithread)
iset_old = -1; jset_old = -1
ikind_old = -1; jkind_old = -1
! Here we loop over gridlevels first, finalising the matrix after each grid level is
! completed. On each grid level, we loop over atom pairs, which will only access
! a single block of each matrix, so with OpenMP, each matrix block is only touched
! by a single thread for each grid level
loop_gridlevels: DO igrid_level = 1, gridlevel_info%ngrid_levels
!$OMP BARRIER
!$OMP DO schedule (dynamic, MAX(1,task_list%npairs(igrid_level)/(nthread*50)))
loop_pairs: DO ipair = 1, task_list%npairs(igrid_level)
loop_tasks: DO itask = task_list%taskstart(ipair, igrid_level), task_list%taskstop(ipair, igrid_level)
ilevel = tasks(itask)%grid_level
img = tasks(itask)%image
iatom = tasks(itask)%iatom
jatom = tasks(itask)%jatom
iset = tasks(itask)%iset
jset = tasks(itask)%jset
ipgf = tasks(itask)%ipgf
jpgf = tasks(itask)%jpgf
! At the start of a block of tasks, get atom data (and kind data, if needed)
IF (itask .EQ. task_list%taskstart(ipair, igrid_level)) THEN
ikind = particle_set(iatom)%atomic_kind%kind_number
jkind = particle_set(jatom)%atomic_kind%kind_number
ra(:) = pbc(particle_set(iatom)%r, cell)
IF (iatom <= jatom) THEN
brow = iatom
bcol = jatom
ELSE
brow = jatom
bcol = iatom
END IF
IF (ikind .NE. ikind_old) THEN
CALL get_qs_kind(qs_kind_set(ikind), &
softb=my_gapw, &
basis_set=orb_basis_set, basis_type="ORB")
CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
first_sgf=first_sgfa, &
lmax=la_max, &
lmin=la_min, &
npgf=npgfa, &
nset=nseta, &
nsgf_set=nsgfa, &
pgf_radius=rpgfa, &
set_radius=set_radius_a, &
sphi=sphi_a, &
zet=zeta)
END IF
IF (jkind .NE. jkind_old) THEN
CALL get_qs_kind(qs_kind_set(jkind), &
softb=my_gapw, &
basis_set=orb_basis_set, basis_type="ORB")
CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
first_sgf=first_sgfb, &
lmax=lb_max, &
lmin=lb_min, &
npgf=npgfb, &
nset=nsetb, &
nsgf_set=nsgfb, &
pgf_radius=rpgfb, &
set_radius=set_radius_b, &
sphi=sphi_b, &
zet=zetb)
END IF
DO i = 1, 3
NULLIFY (vhxc_block(i)%block)
CALL dbcsr_get_block_p(matrix_vhxc_dbasis(i)%matrix, brow, bcol, vhxc_block(i)%block, found)
CPASSERT(found)
END DO
CALL dbcsr_get_block_p(matrix=deltap(img)%matrix, &
row=brow, col=bcol, BLOCK=p_block, found=found)
CPASSERT(found)
ikind_old = ikind
jkind_old = jkind
atom_pair_changed = .TRUE.
ELSE
atom_pair_changed = .FALSE.
END IF
IF (atom_pair_changed .OR. iset_old .NE. iset .OR. jset_old .NE. jset) THEN
ncoa = npgfa(iset)*ncoset(la_max(iset))
sgfa = first_sgfa(1, iset)
ncob = npgfb(jset)*ncoset(lb_max(jset))
sgfb = first_sgfb(1, jset)
IF (iatom <= jatom) THEN
work(1:ncoa, 1:nsgfb(jset)) = MATMUL(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1), &
p_block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1))
pab(1:ncoa, 1:ncob) = MATMUL(work(1:ncoa, 1:nsgfb(jset)), TRANSPOSE(sphi_b(1:ncob, sgfb:sgfb + nsgfb(jset) - 1)))
ELSE
work(1:ncob, 1:nsgfa(iset)) = MATMUL(sphi_b(1:ncob, sgfb:sgfb + nsgfb(jset) - 1), &
p_block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1))
pab(1:ncob, 1:ncoa) = MATMUL(work(1:ncob, 1:nsgfa(iset)), TRANSPOSE(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1)))
END IF
IF (iatom <= jatom) THEN
hab(1:ncoa, 1:ncob) = 0._dp
hdab(:, 1:ncoa, 1:ncob) = 0._dp
hadb(:, 1:ncoa, 1:ncob) = 0._dp
ELSE
hab(1:ncob, 1:ncoa) = 0._dp
hdab(:, 1:ncob, 1:ncoa) = 0._dp
hadb(:, 1:ncob, 1:ncoa) = 0._dp
END IF
iset_old = iset
jset_old = jset
END IF
rab = tasks(itask)%rab
rb(:) = ra(:) + rab(:)
zetp = zeta(ipgf, iset) + zetb(jpgf, jset)
f = zetb(jpgf, jset)/zetp
rp(:) = ra(:) + f*rab(:)
prefactor = EXP(-zeta(ipgf, iset)*f*DOT_PRODUCT(rab, rab))
radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
lb_min=lb_min(jset), lb_max=lb_max(jset), &
ra=ra, rb=rb, rp=rp, &
zetp=zetp, eps=eps_rho_rspace, &
prefactor=prefactor, cutoff=1.0_dp)
na1 = (ipgf - 1)*ncoset(la_max(iset)) + 1
na2 = ipgf*ncoset(la_max(iset))
nb1 = (jpgf - 1)*ncoset(lb_max(jset)) + 1
nb2 = jpgf*ncoset(lb_max(jset))
! check whether we need to use fawzi's generalised collocation scheme
IF (rs_rho(igrid_level)%rs_grid%desc%distributed) THEN
!tasks(4,:) is 0 for replicated, 1 for distributed 2 for exceptional distributed tasks
IF (tasks(itask)%dist_type .EQ. 2) THEN
use_subpatch = .TRUE.
ELSE
use_subpatch = .FALSE.
END IF
ELSE
use_subpatch = .FALSE.
END IF
IF (iatom <= jatom) THEN
IF (iatom == lambda) &
CALL integrate_pgf_product( &
la_max(iset), zeta(ipgf, iset), la_min(iset), &
lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
ra, rab, rs_rho(igrid_level)%rs_grid, cell, &
cube_info(igrid_level), &
hab, o1=na1 - 1, o2=nb1 - 1, &
radius=radius, &
calculate_forces=.TRUE., &
compute_tau=.FALSE., &
use_subpatch=use_subpatch, subpatch_pattern=tasks(itask)%subpatch_pattern, &
hdab=hdab, pab=pab)
IF (jatom == lambda) &
CALL integrate_pgf_product( &
la_max(iset), zeta(ipgf, iset), la_min(iset), &
lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
ra, rab, rs_rho(igrid_level)%rs_grid, cell, &
cube_info(igrid_level), &
hab, o1=na1 - 1, o2=nb1 - 1, &
radius=radius, &
calculate_forces=.TRUE., &
compute_tau=.FALSE., &
use_subpatch=use_subpatch, subpatch_pattern=tasks(itask)%subpatch_pattern, &
hadb=hadb, pab=pab)
ELSE
rab_inv = -rab
IF (iatom == lambda) &
CALL integrate_pgf_product( &
lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
la_max(iset), zeta(ipgf, iset), la_min(iset), &
rb, rab_inv, rs_rho(igrid_level)%rs_grid, cell, &
cube_info(igrid_level), &
hab, o1=nb1 - 1, o2=na1 - 1, &
radius=radius, &
calculate_forces=.TRUE., &
force_a=force_b, force_b=force_a, &
compute_tau=.FALSE., &
use_subpatch=use_subpatch, subpatch_pattern=tasks(itask)%subpatch_pattern, &
hadb=hadb, pab=pab)
IF (jatom == lambda) &
CALL integrate_pgf_product( &
lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
la_max(iset), zeta(ipgf, iset), la_min(iset), &
rb, rab_inv, rs_rho(igrid_level)%rs_grid, cell, &
cube_info(igrid_level), &
hab, o1=nb1 - 1, o2=na1 - 1, &
radius=radius, &
calculate_forces=.TRUE., &
force_a=force_b, force_b=force_a, &
compute_tau=.FALSE., &
use_subpatch=use_subpatch, subpatch_pattern=tasks(itask)%subpatch_pattern, &
hdab=hdab, pab=pab)
END IF
new_set_pair_coming = .FALSE.
atom_pair_done = .FALSE.
IF (itask < task_list%taskstop(ipair, igrid_level)) THEN
ilevel = tasks(itask + 1)%grid_level
img = tasks(itask + 1)%image
iatom = tasks(itask + 1)%iatom
jatom = tasks(itask + 1)%jatom
iset_new = tasks(itask + 1)%iset
jset_new = tasks(itask + 1)%jset
ipgf_new = tasks(itask + 1)%ipgf
jpgf_new = tasks(itask + 1)%jpgf
IF (iset_new .NE. iset .OR. jset_new .NE. jset) THEN
new_set_pair_coming = .TRUE.
END IF
ELSE
! do not forget the last block
new_set_pair_coming = .TRUE.
atom_pair_done = .TRUE.
END IF
IF (new_set_pair_coming) THEN
DO i = 1, 3
hdab(i, :, :) = hdab(i, :, :) + hadb(i, :, :)
IF (iatom <= jatom) THEN
work(1:ncoa, 1:nsgfb(jset)) = MATMUL(hdab(i, 1:ncoa, 1:ncob), sphi_b(1:ncob, sgfb:sgfb + nsgfb(jset) - 1))
vhxc_block(i)%block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1) = &
vhxc_block(i)%block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1) + &
MATMUL(TRANSPOSE(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1)), work(1:ncoa, 1:nsgfb(jset)))
ELSE
work(1:ncob, 1:nsgfa(iset)) = MATMUL(hdab(i, 1:ncob, 1:ncoa), sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1))
vhxc_block(i)%block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) = &
vhxc_block(i)%block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) + &
MATMUL(TRANSPOSE(sphi_b(1:ncob, sgfb:sgfb + nsgfb(jset) - 1)), work(1:ncob, 1:nsgfa(iset)))
END IF
END DO
END IF ! new_set_pair_coming
END DO loop_tasks
END DO loop_pairs
!$OMP END DO
DO i = 1, 3
CALL dbcsr_finalize(matrix_vhxc_dbasis(i)%matrix)
END DO
END DO loop_gridlevels
!$OMP END PARALLEL
IF (distributed_grids) THEN
! Reconstruct H matrix if using distributed RS grids
! note send and recv direction reversed WRT collocate
scatter = .FALSE.
CALL rs_distribute_matrix(rs_descs, matrix_vhxc_dbasis, atom_pair_recv, atom_pair_send, &
dft_control%nimages, scatter=.FALSE.)
END IF
IF (distributed_grids) THEN
CALL dbcsr_deallocate_matrix_set(deltap)
ELSE
DO img = 1, dft_control%nimages
NULLIFY (deltap(img)%matrix)
END DO
DEALLOCATE (deltap)
END IF
DEALLOCATE (pabt, habt, workt, hdabt, hadbt)
IF (ASSOCIATED(rs_rho)) THEN
DO i = 1, SIZE(rs_rho)
CALL rs_grid_release(rs_rho(i)%rs_grid)
END DO
END IF
CALL timestop(handle)
END SUBROUTINE integrate_v_dbasis
! **************************************************************************************************
!> \brief computes matrix elements corresponding to a given potential
!> \param v_rspace ...

102
src/qs_kinetic.F
View file

@ -22,6 +22,7 @@ MODULE qs_kinetic
get_atomic_kind_set
USE basis_set_types, ONLY: gto_basis_set_p_type,&
gto_basis_set_type
USE block_p_types, ONLY: block_p_type
USE cp_control_types, ONLY: dft_control_type
USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set
USE dbcsr_api, ONLY: dbcsr_filter,&
@ -64,6 +65,11 @@ MODULE qs_kinetic
PUBLIC :: build_kinetic_matrix
INTEGER, DIMENSION(1:56), SAVE :: ndod = (/0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, &
1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, &
0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, &
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1/)
CONTAINS
! **************************************************************************************************
@ -78,18 +84,20 @@ CONTAINS
!> \param matrix_p density matrix for force calculation (optional)
!> \param matrixkp_p density matrix for force calculation with kpoints (optional)
!> \param eps_filter Filter final matrix (optional)
!> \param nderivative The number of calculated derivatives
!> \date 11.10.2010
!> \par History
!> Ported from qs_overlap, replaces code in build_core_hamiltonian
!> Refactoring [07.2014] JGH
!> Simplify options and use new kinetic energy integral routine
!> kpoints [08.2014] JGH
!> Include the derivatives [2021] SL, ED
!> \author JGH
!> \version 1.0
! **************************************************************************************************
SUBROUTINE build_kinetic_matrix(ks_env, matrix_t, matrixkp_t, matrix_name, &
basis_type, sab_nl, calculate_forces, matrix_p, matrixkp_p, &
eps_filter)
eps_filter, nderivative)
TYPE(qs_ks_env_type), POINTER :: ks_env
TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
@ -105,13 +113,15 @@ CONTAINS
TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
POINTER :: matrixkp_p
REAL(KIND=dp), INTENT(IN), OPTIONAL :: eps_filter
INTEGER, INTENT(IN), OPTIONAL :: nderivative
INTEGER :: natom
CALL get_ks_env(ks_env, natom=natom)
CALL build_kinetic_matrix_low(ks_env, matrix_t, matrixkp_t, matrix_name, basis_type, &
sab_nl, calculate_forces, matrix_p, matrixkp_p, eps_filter, natom)
sab_nl, calculate_forces, matrix_p, matrixkp_p, eps_filter, natom, &
nderivative)
END SUBROUTINE build_kinetic_matrix
@ -129,9 +139,11 @@ CONTAINS
!> \param matrixkp_p ...
!> \param eps_filter ...
!> \param natom ...
!> \param nderivative ...
! **************************************************************************************************
SUBROUTINE build_kinetic_matrix_low(ks_env, matrix_t, matrixkp_t, matrix_name, basis_type, &
sab_nl, calculate_forces, matrix_p, matrixkp_p, eps_filter, natom)
sab_nl, calculate_forces, matrix_p, matrixkp_p, eps_filter, natom, &
nderivative)
TYPE(qs_ks_env_type), POINTER :: ks_env
TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
@ -148,13 +160,12 @@ CONTAINS
POINTER :: matrixkp_p
REAL(KIND=dp), INTENT(IN), OPTIONAL :: eps_filter
INTEGER, INTENT(IN) :: natom
INTEGER, INTENT(IN), OPTIONAL :: nderivative
CHARACTER(len=*), PARAMETER :: routineN = 'build_kinetic_matrix_low'
INTEGER :: atom_a, handle, iatom, ic, icol, ikind, &
irow, iset, jatom, jkind, jset, ldsab, &
ncoa, ncob, nimg, nkind, nseta, nsetb, &
sgfa, sgfb, slot
INTEGER :: atom_a, handle, i, iatom, ic, icol, ikind, irow, iset, jatom, jkind, jset, ldsab, &
maxder, ncoa, ncob, nder, nimg, nkind, nseta, nsetb, sgfa, sgfb, slot
INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
INTEGER, DIMENSION(3) :: cell
INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
@ -163,16 +174,18 @@ CONTAINS
INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
LOGICAL :: do_forces, do_symmetric, dokp, found, &
trans, use_cell_mapping, use_virial
REAL(KIND=dp) :: f0, ff, rab2, tab
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: kab, pab, qab
REAL(KIND=dp) :: f, f0, ff, rab2, tab
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: pab, qab
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: kab
REAL(KIND=dp), DIMENSION(3) :: force_a, rab
REAL(KIND=dp), DIMENSION(3, 3) :: pv_thread
REAL(KIND=dp), DIMENSION(3, natom) :: force_thread
REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
REAL(KIND=dp), DIMENSION(:, :), POINTER :: k_block, p_block, rpgfa, rpgfb, scon_a, &
scon_b, zeta, zetb
REAL(KIND=dp), DIMENSION(:, :), POINTER :: p_block, rpgfa, rpgfb, scon_a, scon_b, &
zeta, zetb
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: dab
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: k_block
TYPE(dft_control_type), POINTER :: dft_control
TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
@ -216,6 +229,10 @@ CONTAINS
do_forces = .FALSE.
IF (PRESENT(calculate_forces)) do_forces = calculate_forces
nder = 0
IF (PRESENT(nderivative)) nder = nderivative
maxder = ncoset(nder)
! check for symmetry
CPASSERT(SIZE(sab_nl) > 0)
CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
@ -229,7 +246,7 @@ CONTAINS
CALL create_sab_matrix(ks_env, matrixkp_t, matrix_name, basis_set_list, basis_set_list, &
sab_nl, do_symmetric)
ELSE
CALL dbcsr_allocate_matrix_set(matrix_t, 1)
CALL dbcsr_allocate_matrix_set(matrix_t, maxder)
CALL create_sab_matrix(ks_env, matrix_t, matrix_name, basis_set_list, basis_set_list, &
sab_nl, do_symmetric)
END IF
@ -255,12 +272,12 @@ CONTAINS
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP SHARED (do_forces, ldsab, use_cell_mapping, do_symmetric, dokp,&
!$OMP sab_nl, ncoset, use_virial, matrix_t, matrixkp_t,&
!$OMP sab_nl, ncoset, maxder, nder, ndod, use_virial, matrix_t, matrixkp_t,&
!$OMP matrix_p, basis_set_list, atom_of_kind, cell_to_index, matrixkp_p, locks, natom) &
!$OMP PRIVATE (k_block, kab, qab, pab, ikind, jkind, iatom, jatom, rab, cell, basis_set_a, basis_set_b,&
!$OMP PRIVATE (k_block, kab, qab, pab, ikind, jkind, iatom, jatom, rab, cell, basis_set_a, basis_set_b, f, &
!$OMP first_sgfa, la_max, la_min, npgfa, nsgfa, nseta, rpgfa, set_radius_a, ncoa, ncob, force_a, &
!$OMP zeta, first_sgfb, lb_max, lb_min, npgfb, nsetb, rpgfb, set_radius_b, nsgfb, p_block, dab, tab, &
!$OMP slot, zetb, scon_a, scon_b, ic, irow, icol, f0, ff, found, trans, rab2, sgfa, sgfb, iset, jset, &
!$OMP slot, zetb, scon_a, scon_b, i, ic, irow, icol, f0, ff, found, trans, rab2, sgfa, sgfb, iset, jset, &
!$OMP hash, hash1, hash2, iatom8) &
!$OMP REDUCTION (+ : pv_thread, force_thread )
@ -274,11 +291,16 @@ CONTAINS
!$ END DO
!$OMP END DO
ALLOCATE (kab(ldsab, ldsab), qab(ldsab, ldsab))
ALLOCATE (kab(ldsab, ldsab, maxder), qab(ldsab, ldsab))
IF (do_forces) THEN
ALLOCATE (dab(ldsab, ldsab, 3), pab(ldsab, ldsab))
END IF
ALLOCATE (k_block(maxder))
DO i = 1, maxder
NULLIFY (k_block(i)%block)
END DO
!$OMP DO SCHEDULE(GUIDED)
DO slot = 1, sab_nl(1)%nl_size
@ -344,15 +366,17 @@ CONTAINS
f0 = 1.0_dp
ff = 1.0_dp
END IF
NULLIFY (k_block)
IF (dokp) THEN
CALL dbcsr_get_block_p(matrix=matrixkp_t(1, ic)%matrix, &
row=irow, col=icol, BLOCK=k_block, found=found)
row=irow, col=icol, BLOCK=k_block(1)%block, found=found)
CPASSERT(found)
ELSE
CALL dbcsr_get_block_p(matrix=matrix_t(1)%matrix, &
row=irow, col=icol, BLOCK=k_block, found=found)
CPASSERT(found)
DO i = 1, maxder
NULLIFY (k_block(i)%block)
CALL dbcsr_get_block_p(matrix=matrix_t(i)%matrix, &
row=irow, col=icol, BLOCK=k_block(i)%block, found=found)
CPASSERT(found)
END DO
END IF
IF (do_forces) THEN
@ -390,13 +414,13 @@ CONTAINS
IF (do_forces .AND. ASSOCIATED(p_block) .AND. ((iatom /= jatom) .OR. use_virial)) THEN
! Decontract P matrix block
kab = 0.0_dp
CALL block_add("OUT", kab, nsgfa(iset), nsgfb(jset), p_block, sgfa, sgfb, trans=trans)
CALL decontraction(kab, pab, scon_a(:, sgfa:), ncoa, nsgfa(iset), scon_b(:, sgfb:), ncob, nsgfb(jset), &
CALL block_add("OUT", kab(:, :, 1), nsgfa(iset), nsgfb(jset), p_block, sgfa, sgfb, trans=trans)
CALL decontraction(kab(:, :, 1), pab, scon_a(:, sgfa:), ncoa, nsgfa(iset), scon_b(:, sgfb:), ncob, nsgfb(jset), &
trans=trans)
! calculate integrals and derivatives
CALL kinetic(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
rab, kab, dab)
rab, kab(:, :, 1), dab)
CALL force_trace(force_a, dab, pab, ncoa, ncob, 3)
force_thread(:, iatom) = force_thread(:, iatom) + ff*force_a(:)
force_thread(:, jatom) = force_thread(:, jatom) - ff*force_a(:)
@ -405,17 +429,27 @@ CONTAINS
END IF
ELSE
! calclulate integrals
CALL kinetic(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
rab, kab)
IF (nder == 0) THEN
CALL kinetic(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
rab, kab=kab(:, :, 1))
ELSE IF (nder == 1) THEN
CALL kinetic(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
rab, kab=kab(:, :, 1), dab=kab(:, :, 2:4))
END IF
END IF
! Contraction step
CALL contraction(kab, qab, ca=scon_a(:, sgfa:), na=ncoa, ma=nsgfa(iset), &
cb=scon_b(:, sgfb:), nb=ncob, mb=nsgfb(jset), &
trans=trans)
!$ CALL omp_set_lock(locks(hash))
CALL block_add("IN", qab, nsgfa(iset), nsgfb(jset), k_block, sgfa, sgfb, trans=trans)
!$ CALL omp_unset_lock(locks(hash))
DO i = 1, maxder
f = 1.0_dp
IF (ndod(i) == 1 .AND. trans) f = -1.0_dp
! Contraction step
CALL contraction(kab(:, :, i), qab, ca=scon_a(:, sgfa:), na=ncoa, ma=nsgfa(iset), &
cb=scon_b(:, sgfb:), nb=ncob, mb=nsgfb(jset), fscale=f, &
trans=trans)
!$ CALL omp_set_lock(locks(hash))
CALL block_add("IN", qab, nsgfa(iset), nsgfb(jset), k_block(i)%block, sgfa, sgfb, trans=trans)
!$ CALL omp_unset_lock(locks(hash))
END DO
END DO
END DO

72
src/qs_linres_module.F
View file

@ -16,7 +16,8 @@
!> \author MI
! **************************************************************************************************
MODULE qs_linres_module
USE bibliography, ONLY: Weber2009,&
USE bibliography, ONLY: Ditler2021,&
Weber2009,&
cite_reference
USE cp_control_types, ONLY: dft_control_type
USE cp_log_handling, ONLY: cp_get_default_logger,&
@ -40,6 +41,15 @@ MODULE qs_linres_module
section_vals_get_subs_vals,&
section_vals_type,&
section_vals_val_get
USE kinds, ONLY: dp
USE qs_dcdr, ONLY: apt_dR,&
apt_dR_localization,&
dcdr_build_op_dR,&
dcdr_response_dR,&
prepare_per_atom
USE qs_dcdr_utils, ONLY: dcdr_env_cleanup,&
dcdr_env_init,&
dcdr_print
USE qs_density_matrices, ONLY: calculate_density_matrix
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type,&
@ -75,6 +85,7 @@ MODULE qs_linres_module
polar_print,&
polar_response
USE qs_linres_types, ONLY: current_env_type,&
dcdr_env_type,&
epr_env_type,&
issc_env_type,&
linres_control_create,&
@ -102,6 +113,54 @@ MODULE qs_linres_module
CONTAINS
! *****************************************************************************
!> \brief Calculates the derivatives of the MO coefficients dC/dR^lambda_beta
!> wrt to nuclear coordinates. The derivative is index by `beta`, the
!> electric dipole operator by `alpha`.
!> Also calculates the APT
!> P^lambda_alpha,beta = d< mu_alpha >/dR^lambda_beta
!> and calculates the sum rules for the APT elements.
!> \param qs_env ...
!> \param p_env ...
! **************************************************************************************************
SUBROUTINE dcdr_linres(qs_env, p_env)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(qs_p_env_type), POINTER :: p_env
CHARACTER(len=*), PARAMETER :: routineN = 'dcdr_linres', routineP = moduleN//':'//routineN
INTEGER :: beta, latom
TYPE(dcdr_env_type) :: dcdr_env
CALL cite_reference(Ditler2021)
CALL dcdr_env_init(dcdr_env, qs_env)
DO latom = 1, SIZE(dcdr_env%list_of_atoms)
dcdr_env%lambda = dcdr_env%list_of_atoms(latom)
CALL prepare_per_atom(dcdr_env, qs_env)
DO beta = 1, 3 ! in every direction
dcdr_env%beta = beta
dcdr_env%deltaR(dcdr_env%beta, dcdr_env%lambda) = 1._dp
CALL dcdr_build_op_dR(dcdr_env, qs_env)
CALL dcdr_response_dR(dcdr_env, p_env, qs_env)
IF (.NOT. dcdr_env%localized_psi0) THEN
CALL apt_dR(qs_env, dcdr_env)
ELSE IF (dcdr_env%localized_psi0) THEN
CALL apt_dR_localization(qs_env, dcdr_env)
END IF
END DO !beta
dcdr_env%apt_total_dcdr(:, :, dcdr_env%lambda) = &
dcdr_env%apt_el_dcdr(:, :, dcdr_env%lambda) + dcdr_env%apt_nuc_dcdr(:, :, dcdr_env%lambda)
END DO !lambda
CALL dcdr_print(dcdr_env, qs_env)
CALL dcdr_env_cleanup(qs_env, dcdr_env)
END SUBROUTINE dcdr_linres
! **************************************************************************************************
!> \brief Driver for the linear response calculatios
!> \param force_env ...
@ -160,7 +219,7 @@ CONTAINS
CHARACTER(LEN=*), PARAMETER :: routineN = 'linres_calculation_low'
INTEGER :: handle, iounit
LOGICAL :: epr_present, issc_present, &
LOGICAL :: dcdr_present, epr_present, issc_present, &
lr_calculation, nmr_present, &
polar_present
TYPE(cp_logger_type), POINTER :: logger
@ -176,6 +235,7 @@ CONTAINS
epr_present = .FALSE.
issc_present = .FALSE.
polar_present = .FALSE.
dcdr_present = .FALSE.
NULLIFY (dft_control, p_env, linres_control, logger, prop_section, lr_section)
logger => cp_get_default_logger()
@ -219,6 +279,14 @@ CONTAINS
CALL polar_linres(qs_env, p_env)
END IF
! Nuclear Position Perturbation
prop_section => section_vals_get_subs_vals(lr_section, "dcdr")
CALL section_vals_get(prop_section, explicit=dcdr_present)
IF (dcdr_present) THEN
CALL dcdr_linres(qs_env, p_env)
END IF
! Other possible LR calculations can be introduced here
CALL p_env_release(p_env)

41
src/qs_linres_types.F
View file

@ -17,9 +17,11 @@ MODULE qs_linres_types
gto_basis_set_type
USE cp_array_utils, ONLY: cp_2d_i_p_type,&
cp_2d_r_p_type
USE cp_fm_struct, ONLY: cp_fm_struct_type
USE cp_fm_types, ONLY: cp_fm_p_type,&
cp_fm_release
USE dbcsr_api, ONLY: dbcsr_p_type
USE dbcsr_api, ONLY: dbcsr_p_type,&
dbcsr_type
USE kinds, ONLY: dp
USE qs_grid_atom, ONLY: grid_atom_type
USE qs_harmonics_atom, ONLY: harmonics_atom_type
@ -210,13 +212,48 @@ MODULE qs_linres_types
POINTER :: jrho_b_h_iii, jrho_b_s_iii
END TYPE jrho_atom_type
! \param type for dC/dR calculation
TYPE dcdr_env_type
INTEGER :: nao, nmo, orb_center, beta, lambda, output_unit
TYPE(dbcsr_type), POINTER :: perturbed_dm_correction
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_hc, matrix_s1, matrix_t1, matrix_s, matrix_t, &
matrix_ppnl_1, matrix_d_vhxc_dR, matrix_core_charge_1, &
matrix_nosym_temp, &
moments, &
matrix_apply_op_constant, &
hamiltonian1
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_vhxc_perturbed_basis, &
matrix_difdip
REAL(dp), DIMENSION(:, :), POINTER :: deltaR, delta_basis_function
REAL(dp), DIMENSION(:, :, :, :), POINTER :: apt_subset, apt_at_dcdr_per_center
TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: mo_coeff, dCR, dCR_prime, op_dR, chc
CHARACTER(LEN=30) :: orb_center_name
TYPE(cp_2d_i_p_type), DIMENSION(:), POINTER :: center_list
TYPE(cp_2d_r_p_type), DIMENSION(:), POINTER :: centers_set
INTEGER, DIMENSION(2) :: nbr_center, nstates
REAL(dp), DIMENSION(3) :: ref_point, dipole_pos
LOGICAL :: localized_psi0
INTEGER, POINTER :: list_of_atoms(:)
LOGICAL :: distributed_origin
TYPE(cp_fm_struct_type), POINTER :: aoao_fm_struct
TYPE(cp_fm_struct_type), POINTER :: momo_fm_struct
TYPE(cp_fm_struct_type), POINTER :: homohomo_fm_struct
TYPE(cp_fm_struct_type), POINTER :: likemos_fm_struct
REAL(dp), DIMENSION(:, :, :), POINTER :: apt_el_dcdr, apt_nuc_dcdr, apt_total_dcdr
REAL(dp), DIMENSION(:, :, :, :), POINTER :: apt_el_dcdr_per_center, apt_el_dcdr_per_subset
END TYPE dcdr_env_type
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_linres_types'
! *** Public data types ***
PUBLIC :: linres_control_type, &
nmr_env_type, issc_env_type, jrho_atom_type, &
epr_env_type, &
epr_env_type, dcdr_env_type, &
nablavks_atom_type, current_env_type, &
realspaces_grid_p_type, polar_env_type

219
src/qs_moments.F
View file

@ -16,6 +16,7 @@ MODULE qs_moments
USE ai_moments, ONLY: contract_cossin,&
cossin,&
diff_momop,&
diff_momop2,&
moment
USE atomic_kind_types, ONLY: atomic_kind_type,&
get_atomic_kind
@ -103,6 +104,7 @@ MODULE qs_moments
PUBLIC :: build_berry_moment_matrix, build_local_moment_matrix
PUBLIC :: build_berry_kpoint_matrix
PUBLIC :: qs_moment_berry_phase, qs_moment_locop
PUBLIC :: dipole_deriv_ao
CONTAINS
@ -2498,4 +2500,221 @@ CONTAINS
CALL timestop(handle)
END SUBROUTINE calculate_commutator_nl_terms
! *****************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param difdip ...
!> \param deltaR ...
!> \param order ...
!> \param rcc ...
!> \note calculate matrix elements <a|r_beta |db/dR_alpha > + <da/dR_alpha | r_beta | b >
!> be aware: < a | r_beta| db/dR_alpha > = - < da/dR_alpha | r_beta | b > only valid
!> if alpha .neq.beta
!> if alpha=beta: < a | r_beta| db/dR_alpha > = - < da/dR_alpha | r_beta | b > - < a | b >
!> modified from qs_efield_mo_derivatives
!> SL July 2015
! **************************************************************************************************
SUBROUTINE dipole_deriv_ao(qs_env, difdip, deltaR, order, rcc)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(dbcsr_p_type), DIMENSION(:, :), &
INTENT(INOUT), POINTER :: difdip
REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
POINTER :: deltaR
INTEGER, INTENT(IN) :: order
REAL(KIND=dp), DIMENSION(3), OPTIONAL :: rcc
CHARACTER(LEN=*), PARAMETER :: routineN = 'dipole_deriv_ao', &
routineP = moduleN//':'//routineN
INTEGER :: handle, i, iatom, icol, idir, ikind, inode, irow, iset, j, jatom, jkind, jset, &
last_jatom, lda, ldab, ldb, M_dim, maxsgf, natom, ncoa, ncob, nkind, nseta, nsetb, sgfa, &
sgfb
INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
npgfb, nsgfa, nsgfb
INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
LOGICAL :: found
REAL(dp) :: dab
REAL(dp), DIMENSION(3) :: ra, rab, rac, rb, rbc, rc
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: work
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: difmab
REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: mab
REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: difmab2
TYPE(block_p_type), ALLOCATABLE, DIMENSION(:, :) :: mint, mint2
TYPE(cell_type), POINTER :: cell
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_all, sab_orb
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(qs_kind_type), POINTER :: qs_kind
CALL timeset(routineN, handle)
NULLIFY (cell, particle_set, qs_kind_set, sab_orb, sab_all)
CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set, &
qs_kind_set=qs_kind_set, sab_orb=sab_orb, sab_all=sab_all)
CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
maxco=ldab, maxsgf=maxsgf)
nkind = SIZE(qs_kind_set)
natom = SIZE(particle_set)
M_dim = ncoset(order) - 1
IF (PRESENT(rcc)) THEN
rc = rcc
ELSE
rc = 0._dp
END IF
ALLOCATE (basis_set_list(nkind))
ALLOCATE (mab(ldab, ldab, M_dim))
ALLOCATE (difmab2(ldab, ldab, M_dim, 3))
ALLOCATE (work(ldab, maxsgf))
ALLOCATE (mint(3, 3))
ALLOCATE (mint2(3, 3))
mab(1:ldab, 1:ldab, 1:M_dim) = 0.0_dp
difmab2(1:ldab, 1:ldab, 1:M_dim, 1:3) = 0.0_dp
work(1:ldab, 1:maxsgf) = 0.0_dp
DO i = 1, 3
DO j = 1, 3
NULLIFY (mint(i, j)%block)
NULLIFY (mint2(i, j)%block)
END DO
END DO
! Set the basis_set_list(nkind) to point to the corresponding basis sets
DO ikind = 1, nkind
qs_kind => qs_kind_set(ikind)
CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a)
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_all)
DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
iatom=iatom, jatom=jatom, r=rab)
basis_set_a => basis_set_list(ikind)%gto_basis_set
basis_set_b => basis_set_list(jkind)%gto_basis_set
IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
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
IF (inode == 1) last_jatom = 0
! this guarentees minimum image convention
! anything else would not make sense
IF (jatom == last_jatom) THEN
CYCLE
END IF
last_jatom = jatom
irow = iatom
icol = jatom
DO i = 1, 3
DO j = 1, 3
NULLIFY (mint(i, j)%block)
CALL dbcsr_get_block_p(matrix=difdip(i, j)%matrix, &
row=irow, col=icol, BLOCK=mint(i, j)%block, &
found=found)
CPASSERT(found)
END DO
END DO
ra(:) = particle_set(iatom)%r(:)
rb(:) = particle_set(jatom)%r(:)
rab(:) = pbc(rb, ra, cell)
rac(:) = pbc(ra - rc, cell)
rbc(:) = pbc(rb - rc, cell)
dab = SQRT(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3))
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)
ldab = MAX(ncoa, ncob)
lda = ncoset(la_max(iset))*npgfa(iset)
ldb = ncoset(lb_max(jset))*npgfb(jset)
ALLOCATE (difmab(lda, ldb, M_dim, 3))
! Calculate integral (da|r|b)
CALL diff_momop2(la_max(iset), npgfa(iset), zeta(:, iset), &
rpgfa(:, iset), la_min(iset), lb_max(jset), npgfb(jset), &
zetb(:, jset), rpgfb(:, jset), lb_min(jset), order, rac, rbc, &
difmab, deltaR=deltaR, iatom=iatom, jatom=jatom)
! *** Contraction step ***
DO idir = 1, 3 ! derivative of AO function
DO j = 1, 3 ! position operator r_j
CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
1.0_dp, difmab(1, 1, j, idir), SIZE(difmab, 1), &
sphi_b(1, sgfb), SIZE(sphi_b, 1), &
0.0_dp, work(1, 1), SIZE(work, 1))
CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
work(1, 1), SIZE(work, 1), &
1.0_dp, mint(j, idir)%block(sgfa, sgfb), &
SIZE(mint(j, idir)%block, 1))
END DO !j
END DO !idir
DEALLOCATE (difmab)
END DO !jset
END DO !iset
END DO!iterator
CALL neighbor_list_iterator_release(nl_iterator)
DO i = 1, 3
DO j = 1, 3
NULLIFY (mint(i, j)%block)
END DO
END DO
DEALLOCATE (mab, difmab2, basis_set_list, work, mint, mint2)
CALL timestop(handle)
END SUBROUTINE dipole_deriv_ao
END MODULE qs_moments

58
src/qs_operators_ao.F
View file

@ -25,6 +25,7 @@ MODULE qs_operators_ao
USE cp_para_types, ONLY: cp_para_env_type
USE dbcsr_api, ONLY: dbcsr_get_block_p,&
dbcsr_get_matrix_type,&
dbcsr_has_symmetry,&
dbcsr_p_type,&
dbcsr_type_antisymmetric,&
dbcsr_type_no_symmetry
@ -81,7 +82,7 @@ CONTAINS
INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
npgfb, nsgfa, nsgfb
INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
LOGICAL :: found, new_atom_b
LOGICAL :: do_symmetric, found, new_atom_b
REAL(KIND=dp) :: dab, rab2
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: work
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: intab, rr_work
@ -97,14 +98,14 @@ CONTAINS
TYPE(neighbor_list_iterator_p_type), &
DIMENSION(:), POINTER :: nl_iterator
TYPE(neighbor_list_set_p_type), DIMENSION(:), &
POINTER :: sab_orb
POINTER :: sab_nl
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(qs_kind_type), POINTER :: qs_kind
CALL timeset(routineN, handle)
NULLIFY (cell, sab_orb, qs_kind_set, particle_set, para_env)
NULLIFY (cell, sab_nl, qs_kind_set, particle_set, para_env)
NULLIFY (logger)
logger => cp_get_default_logger()
@ -114,12 +115,18 @@ CONTAINS
particle_set=particle_set, &
neighbor_list_id=neighbor_list_id, &
para_env=para_env, &
sab_orb=sab_orb, &
cell=cell)
nkind = SIZE(qs_kind_set)
natom = SIZE(particle_set)
! Take into account the symmetry of the input matrix
do_symmetric = dbcsr_has_symmetry(matrix(1)%matrix)
IF (do_symmetric) THEN
CALL get_qs_env(qs_env=qs_env, sab_orb=sab_nl)
ELSE
CALL get_qs_env(qs_env=qs_env, sab_all=sab_nl)
END IF
! *** Allocate work storage ***
CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
@ -145,7 +152,7 @@ CONTAINS
NULLIFY (basis_set_list(ikind)%gto_basis_set)
END IF
END DO
CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
CALL neighbor_list_iterator_create(nl_iterator, sab_nl)
DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
iatom=iatom, jatom=jatom, r=rab)
@ -185,12 +192,17 @@ CONTAINS
END IF
IF (new_atom_b) THEN
IF (iatom <= jatom) THEN
IF (do_symmetric) THEN
IF (iatom <= jatom) THEN
irow = iatom
icol = jatom
ELSE
irow = jatom
icol = iatom
END IF
ELSE
irow = iatom
icol = jatom
ELSE
irow = jatom
icol = iatom
END IF
DO i = 1, 3
@ -233,22 +245,30 @@ CONTAINS
sphi_b(1, sgfb), SIZE(sphi_b, 1), &
0.0_dp, work(1, 1), SIZE(work, 1))
IF (iatom <= jatom) THEN
IF (do_symmetric) THEN
IF (iatom <= jatom) THEN
CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
work(1, 1), SIZE(work, 1), &
1.0_dp, integral(i)%block(sgfa, sgfb), &
SIZE(integral(i)%block, 1))
ELSE
CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
-1.0_dp, work(1, 1), SIZE(work, 1), &
sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1.0_dp, integral(i)%block(sgfb, sgfa), &
SIZE(integral(i)%block, 1))
END IF
ELSE
CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
work(1, 1), SIZE(work, 1), &
1.0_dp, integral(i)%block(sgfa, sgfb), &
SIZE(integral(i)%block, 1))
ELSE
CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
-1.0_dp, work(1, 1), SIZE(work, 1), &
sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1.0_dp, integral(i)%block(sgfb, sgfa), &
SIZE(integral(i)%block, 1))
END IF
END DO

10
tests/QS/regtest-dcdr/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_apt.inp 96 1e-06 -0.879586
h2o_apt_loc.inp 96 1e-06 -0.879586
h2o_apt_pbc.inp 96 1e-06 -0.938707
h2o_apt_pbc_loc.inp 96 1e-06 -0.938770
#EOF

85
tests/QS/regtest-dcdr/h2o_apt.inp Normal file
View file

@ -0,0 +1,85 @@
###################################
@SET RUN_TYPE ENERGY_FORCE
@SET CUTOFF 200
@SET FUNCTIONAL LDA
@SET PRINT_LEVEL MEDIUM
@SET BASIS_SET_FILE_NAME GTH_BASIS_SETS
@SET BASIS_SET SZV-GTH
@SET EPS_SCF 1.08E-5
@SET EPS_LINRES 5.0E-5
###################################
&GLOBAL
PROJECT second
RUN_TYPE $RUN_TYPE
PRINT_LEVEL $PRINT_LEVEL
&END GLOBAL
&FORCE_EVAL
&PROPERTIES
&LINRES
MAX_ITER 1000
PRECONDITIONER FULL_SINGLE_INVERSE
EPS $EPS_LINRES
&PRINT
&PROGRAM_RUN_INFO
&END
&END PRINT
&DCDR
&PRINT
&APT
FILENAME __STD_OUT__
&END
&END PRINT
&END DCDR
&END LINRES
&END PROPERTIES
METHOD Quickstep
&DFT
CHARGE 0
BASIS_SET_FILE_NAME $BASIS_SET_FILE_NAME
POTENTIAL_FILE_NAME POTENTIAL
&MGRID
NGRIDS 1
CUTOFF $CUTOFF
&END MGRID
&QS
EXTRAPOLATION ASPC
EXTRAPOLATION_ORDER 3
METHOD GPW
&END QS
&SCF
SCF_GUESS ATOMIC
EPS_SCF $EPS_SCF
&OT
PRECONDITIONER FULL_SINGLE_INVERSE
&END OT
&END SCF
&PRINT
&MOMENTS
PERIODIC FALSE
&END MOMENTS
&END PRINT
&POISSON
POISSON_SOLVER ANALYTIC
PERIODIC NONE
&END
&XC
&XC_FUNCTIONAL $FUNCTIONAL
&END
&END XC
&END DFT
&SUBSYS
&CELL
PERIODIC NONE
ABC [angstrom] 5.0 5.0 5.0
&END CELL
&KIND DEFAULT
BASIS_SET $BASIS_SET
POTENTIAL GTH-$FUNCTIONAL
&END KIND
&COORD
O 0.000000 0.000000 0.000000
H 0.000000 0.769665 -0.591648
H 0.000000 -0.769665 -0.591648
&END COORD
&END SUBSYS
&END FORCE_EVAL

88
tests/QS/regtest-dcdr/h2o_apt_loc.inp Normal file
View file

@ -0,0 +1,88 @@
###################################
@SET RUN_TYPE ENERGY_FORCE
@SET CUTOFF 200
@SET FUNCTIONAL LDA
@SET PRINT_LEVEL MEDIUM
@SET BASIS_SET_FILE_NAME GTH_BASIS_SETS
@SET BASIS_SET SZV-GTH
@SET EPS_SCF 1.08E-5
@SET EPS_LINRES 5.0E-5
###################################
&GLOBAL
PROJECT second
RUN_TYPE $RUN_TYPE
PRINT_LEVEL $PRINT_LEVEL
&END GLOBAL
&FORCE_EVAL
&PROPERTIES
&LINRES
&LOCALIZE
EPS_LOCALIZATION 1.E-9
&END LOCALIZE
MAX_ITER 1000
PRECONDITIONER FULL_SINGLE_INVERSE
EPS $EPS_LINRES
&PRINT
&PROGRAM_RUN_INFO
&END
&END PRINT
&DCDR
&PRINT
&APT
FILENAME __STD_OUT__
&END
&END PRINT
&END DCDR
&END LINRES
&END PROPERTIES
METHOD Quickstep
&DFT
CHARGE 0
BASIS_SET_FILE_NAME $BASIS_SET_FILE_NAME
POTENTIAL_FILE_NAME POTENTIAL
&MGRID
NGRIDS 1
CUTOFF $CUTOFF
&END MGRID
&QS
EXTRAPOLATION ASPC
EXTRAPOLATION_ORDER 3
METHOD GPW
&END QS
&SCF
SCF_GUESS ATOMIC
EPS_SCF $EPS_SCF
&OT
PRECONDITIONER FULL_SINGLE_INVERSE
&END OT
&END SCF
&PRINT
&MOMENTS
PERIODIC FALSE
&END MOMENTS
&END PRINT
&POISSON
POISSON_SOLVER ANALYTIC
PERIODIC NONE
&END
&XC
&XC_FUNCTIONAL $FUNCTIONAL
&END
&END XC
&END DFT
&SUBSYS
&CELL
PERIODIC NONE
ABC [angstrom] 5.0 5.0 5.0
&END CELL
&KIND DEFAULT
BASIS_SET $BASIS_SET
POTENTIAL GTH-$FUNCTIONAL
&END KIND
&COORD
O 0.000000 0.000000 0.000000
H 0.000000 0.769665 -0.591648
H 0.000000 -0.769665 -0.591648
&END COORD
&END SUBSYS
&END FORCE_EVAL

85
tests/QS/regtest-dcdr/h2o_apt_pbc.inp Normal file
View file

@ -0,0 +1,85 @@
###################################
@SET RUN_TYPE ENERGY_FORCE
@SET CUTOFF 200
@SET FUNCTIONAL LDA
@SET PRINT_LEVEL MEDIUM
@SET BASIS_SET_FILE_NAME GTH_BASIS_SETS
@SET BASIS_SET SZV-GTH
@SET EPS_SCF 1.08E-5
@SET EPS_LINRES 5.0E-5
###################################
&GLOBAL
PROJECT second
RUN_TYPE $RUN_TYPE
PRINT_LEVEL $PRINT_LEVEL
&END GLOBAL
&FORCE_EVAL
&PROPERTIES
&LINRES
MAX_ITER 1000
PRECONDITIONER FULL_SINGLE_INVERSE
EPS $EPS_LINRES
&PRINT
&PROGRAM_RUN_INFO
&END
&END PRINT
&DCDR
&PRINT
&APT
FILENAME __STD_OUT__
&END
&END PRINT
&END DCDR
&END LINRES
&END PROPERTIES
METHOD Quickstep
&DFT
CHARGE 0
BASIS_SET_FILE_NAME $BASIS_SET_FILE_NAME
POTENTIAL_FILE_NAME POTENTIAL
&MGRID
NGRIDS 1
CUTOFF $CUTOFF
&END MGRID
&QS
EXTRAPOLATION ASPC
EXTRAPOLATION_ORDER 3
METHOD GPW
&END QS
&SCF
SCF_GUESS ATOMIC
EPS_SCF $EPS_SCF
&OT
PRECONDITIONER FULL_SINGLE_INVERSE
&END OT
&END SCF
&PRINT
&MOMENTS
PERIODIC TRUE
&END MOMENTS
&END PRINT
&POISSON
POISSON_SOLVER PERIODIC
PERIODIC XYZ
&END
&XC
&XC_FUNCTIONAL $FUNCTIONAL
&END
&END XC
&END DFT
&SUBSYS
&CELL
PERIODIC XYZ
ABC [angstrom] 5.0 5.0 5.0
&END CELL
&KIND DEFAULT
BASIS_SET $BASIS_SET
POTENTIAL GTH-$FUNCTIONAL
&END KIND
&COORD
O 0.000000 0.000000 0.000000
H 0.000000 0.769665 -0.591648
H 0.000000 -0.769665 -0.591648
&END COORD
&END SUBSYS
&END FORCE_EVAL

88
tests/QS/regtest-dcdr/h2o_apt_pbc_loc.inp Normal file
View file

@ -0,0 +1,88 @@
###################################
@SET RUN_TYPE ENERGY_FORCE
@SET CUTOFF 200
@SET FUNCTIONAL LDA
@SET PRINT_LEVEL MEDIUM
@SET BASIS_SET_FILE_NAME GTH_BASIS_SETS
@SET BASIS_SET SZV-GTH
@SET EPS_SCF 1.08E-5
@SET EPS_LINRES 5.0E-5
###################################
&GLOBAL
PROJECT second
RUN_TYPE $RUN_TYPE
PRINT_LEVEL $PRINT_LEVEL
&END GLOBAL
&FORCE_EVAL
&PROPERTIES
&LINRES
&LOCALIZE
EPS_LOCALIZATION 1.E-9
&END LOCALIZE
MAX_ITER 1000
PRECONDITIONER FULL_SINGLE_INVERSE
EPS $EPS_LINRES
&PRINT
&PROGRAM_RUN_INFO
&END
&END PRINT
&DCDR
&PRINT
&APT
FILENAME __STD_OUT__
&END
&END PRINT
&END DCDR
&END LINRES
&END PROPERTIES
METHOD Quickstep
&DFT
CHARGE 0
BASIS_SET_FILE_NAME $BASIS_SET_FILE_NAME
POTENTIAL_FILE_NAME POTENTIAL
&MGRID
NGRIDS 1
CUTOFF $CUTOFF
&END MGRID
&QS
EXTRAPOLATION ASPC
EXTRAPOLATION_ORDER 3
METHOD GPW
&END QS
&SCF
SCF_GUESS ATOMIC
EPS_SCF $EPS_SCF
&OT
PRECONDITIONER FULL_SINGLE_INVERSE
&END OT
&END SCF
&PRINT
&MOMENTS
PERIODIC TRUE
&END MOMENTS
&END PRINT
&POISSON
POISSON_SOLVER PERIODIC
PERIODIC XYZ
&END
&XC
&XC_FUNCTIONAL $FUNCTIONAL
&END
&END XC
&END DFT
&SUBSYS
&CELL
PERIODIC XYZ
ABC [angstrom] 5.0 5.0 5.0
&END CELL
&KIND DEFAULT
BASIS_SET $BASIS_SET
POTENTIAL GTH-$FUNCTIONAL
&END KIND
&COORD
O 0.000000 0.000000 0.000000
H 0.000000 0.769665 -0.591648
H 0.000000 -0.769665 -0.591648
&END COORD
&END SUBSYS
&END FORCE_EVAL

1
tests/TEST_DIRS
View file

@ -170,6 +170,7 @@ QS/regtest-ot-2
Fist/regtest-7-2
QS/regtest-sparsity libint
QS/regtest-polar
QS/regtest-dcdr
SE/regtest-4
QS/regtest-nonortho
QS/regtest-hybrid-3 libint

3
tests/TEST_TYPES
View file

@ -1,4 +1,4 @@
95
96
Total energy:!3
MD| Potential energy!5
Total energy \[eV\]:!4
@ -94,6 +94,7 @@ FCIDUMP| Checksum: !3
SPGR| SPACE GROUP NUMBER: !5
KS CSR write| !4
Fermi energy: !3
APT | 1 2 !7
#
# these are the tests the can be selected for regtesting.
# do regtest will grep for test_grep (first column) and look if the numeric value