From 67208cef2f55ffbaf6de46701c20b6094b782f86 Mon Sep 17 00:00:00 2001 From: Edward Ditler Date: Wed, 17 Nov 2021 15:21:58 +0100 Subject: [PATCH] Add analytical derivatives of MO coefficients wrt nuclear coordinates --- src/aobasis/ai_moments.F | 105 ++ src/aobasis/ai_oneelectron.F | 60 +- src/aobasis/ai_overlap_ppl.F | 21 +- src/common/bibliography.F | 18 +- src/core_ppl.F | 115 ++- src/core_ppnl.F | 103 +- src/grid/common/grid_constants.h | 11 +- src/grid/common/grid_prepare_pab.h | 99 ++ src/grid/grid_api.F | 32 +- src/grid/grid_replay.c | 2 +- src/grid/ref/grid_ref_collocate.c | 2 +- src/grid/ref/grid_ref_integrate.c | 12 +- src/grid/ref/grid_ref_integrate.h | 3 +- src/grid/ref/grid_ref_task_list.c | 1 + src/input_cp2k_properties_dft.F | 113 +++ src/qs_collocate_density.F | 556 ++++++++++- src/qs_dcdr.F | 669 +++++++++++++ src/qs_dcdr_ao.F | 789 +++++++++++++++ src/qs_dcdr_utils.F | 1077 +++++++++++++++++++++ src/qs_integrate_potential.F | 6 +- src/qs_integrate_potential_product.F | 492 +++++++++- src/qs_kinetic.F | 102 +- src/qs_linres_module.F | 72 +- src/qs_linres_types.F | 41 +- src/qs_moments.F | 219 +++++ src/qs_operators_ao.F | 58 +- tests/QS/regtest-dcdr/TEST_FILES | 10 + tests/QS/regtest-dcdr/h2o_apt.inp | 85 ++ tests/QS/regtest-dcdr/h2o_apt_loc.inp | 88 ++ tests/QS/regtest-dcdr/h2o_apt_pbc.inp | 85 ++ tests/QS/regtest-dcdr/h2o_apt_pbc_loc.inp | 88 ++ tests/TEST_DIRS | 1 + tests/TEST_TYPES | 3 +- 33 files changed, 5023 insertions(+), 115 deletions(-) create mode 100644 src/qs_dcdr.F create mode 100644 src/qs_dcdr_ao.F create mode 100644 src/qs_dcdr_utils.F create mode 100644 tests/QS/regtest-dcdr/TEST_FILES create mode 100644 tests/QS/regtest-dcdr/h2o_apt.inp create mode 100644 tests/QS/regtest-dcdr/h2o_apt_loc.inp create mode 100644 tests/QS/regtest-dcdr/h2o_apt_pbc.inp create mode 100644 tests/QS/regtest-dcdr/h2o_apt_pbc_loc.inp diff --git a/src/aobasis/ai_moments.F b/src/aobasis/ai_moments.F index bcbdd684a7..ff2e847c4c 100644 --- a/src/aobasis/ai_moments.F +++ b/src/aobasis/ai_moments.F @@ -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 ... diff --git a/src/aobasis/ai_oneelectron.F b/src/aobasis/ai_oneelectron.F index 14e658d598..ae8a2e90af 100644 --- a/src/aobasis/ai_oneelectron.F +++ b/src/aobasis/ai_oneelectron.F @@ -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 diff --git a/src/aobasis/ai_overlap_ppl.F b/src/aobasis/ai_overlap_ppl.F index 5790b91732..ba25f4d907 100644 --- a/src/aobasis/ai_overlap_ppl.F +++ b/src/aobasis/ai_overlap_ppl.F @@ -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) diff --git a/src/common/bibliography.F b/src/common/bibliography.F index 57f34d6b7f..63bece1a4d 100644 --- a/src/common/bibliography.F +++ b/src/common/bibliography.F @@ -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 diff --git a/src/core_ppl.F b/src/core_ppl.F index 2bf4663530..de2f192f25 100644 --- a/src/core_ppl.F +++ b/src/core_ppl.F @@ -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 diff --git a/src/core_ppnl.F b/src/core_ppnl.F index ee21a08ab4..90d853f491 100644 --- a/src/core_ppnl.F +++ b/src/core_ppnl.F @@ -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 diff --git a/src/grid/common/grid_constants.h b/src/grid/common/grid_constants.h index a0bdadbfe0..680bc32bea 100644 --- a/src/grid/common/grid_constants.h +++ b/src/grid/common/grid_constants.h @@ -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 { diff --git a/src/grid/common/grid_prepare_pab.h b/src/grid/common/grid_prepare_pab.h index 3aa29bfe7b..06ce32cc56 100644 --- a/src/grid/common/grid_prepare_pab.h +++ b/src/grid/common/grid_prepare_pab.h @@ -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; diff --git a/src/grid/grid_api.F b/src/grid/grid_api.F index eb7c39c294..118eb233c2 100644 --- a/src/grid/grid_api.F +++ b/src/grid/grid_api.F @@ -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)) & diff --git a/src/grid/grid_replay.c b/src/grid/grid_replay.c index 70407ba5c7..4bbe06661b 100644 --- a/src/grid/grid_replay.c +++ b/src/grid/grid_replay.c @@ -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++) { diff --git a/src/grid/ref/grid_ref_collocate.c b/src/grid/ref/grid_ref_collocate.c index 8a16fd09e8..e427987e59 100644 --- a/src/grid/ref/grid_ref_collocate.c +++ b/src/grid/ref/grid_ref_collocate.c @@ -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++) { diff --git a/src/grid/ref/grid_ref_integrate.c b/src/grid/ref/grid_ref_integrate.c index 66c191ec6a..59966e19e4 100644 --- a/src/grid/ref/grid_ref_integrate.c +++ b/src/grid/ref/grid_ref_integrate.c @@ -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++) { diff --git a/src/grid/ref/grid_ref_integrate.h b/src/grid/ref/grid_ref_integrate.h index b0e62fb877..36ea8dacb5 100644 --- a/src/grid/ref/grid_ref_integrate.h +++ b/src/grid/ref/grid_ref_integrate.h @@ -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 diff --git a/src/grid/ref/grid_ref_task_list.c b/src/grid/ref/grid_ref_task_list.c index 0d867476b1..cfa9e61a36 100644 --- a/src/grid/ref/grid_ref_task_list.c +++ b/src/grid/ref/grid_ref_task_list.c @@ -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 diff --git a/src/input_cp2k_properties_dft.F b/src/input_cp2k_properties_dft.F index cd807c0b39..b0d7eeb5e5 100644 --- a/src/input_cp2k_properties_dft.F +++ b/src/input_cp2k_properties_dft.F @@ -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 diff --git a/src/qs_collocate_density.F b/src/qs_collocate_density.F index 8624dd000e..eb83843080 100644 --- a/src/qs_collocate_density.F +++ b/src/qs_collocate_density.F @@ -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 ... diff --git a/src/qs_dcdr.F b/src/qs_dcdr.F new file mode 100644 index 0000000000..0553018af0 --- /dev/null +++ b/src/qs_dcdr.F @@ -0,0 +1,669 @@ +!--------------------------------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations ! +! Copyright 2000-2021 CP2K developers group ! +! ! +! 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 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 diff --git a/src/qs_dcdr_ao.F b/src/qs_dcdr_ao.F new file mode 100644 index 0000000000..63b323bbf1 --- /dev/null +++ b/src/qs_dcdr_ao.F @@ -0,0 +1,789 @@ +!--------------------------------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations ! +! Copyright 2000-2021 CP2K developers group ! +! ! +! 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 diff --git a/src/qs_dcdr_utils.F b/src/qs_dcdr_utils.F new file mode 100644 index 0000000000..efec4badd0 --- /dev/null +++ b/src/qs_dcdr_utils.F @@ -0,0 +1,1077 @@ +!--------------------------------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations ! +! Copyright 2000-2021 CP2K developers group ! +! ! +! 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_utils +!#include "./common/cp_common_uses.f90" + USE cell_types, ONLY: cell_type,& + get_cell + USE cp_control_types, ONLY: dft_control_type + USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl + USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply,& + dbcsr_allocate_matrix_set,& + dbcsr_deallocate_matrix_set + USE cp_files, ONLY: close_file,& + open_file + USE cp_fm_struct, ONLY: cp_fm_struct_create,& + cp_fm_struct_release + USE cp_fm_types, ONLY: & + cp_fm_create, cp_fm_get_info, cp_fm_get_submatrix, cp_fm_p_type, cp_fm_release, & + cp_fm_set_all, cp_fm_set_submatrix, 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_get_default_io_unit,& + cp_logger_type,& + cp_to_string + USE cp_output_handling, ONLY: cp_p_file,& + cp_print_key_finished_output,& + cp_print_key_generate_filename,& + cp_print_key_should_output,& + cp_print_key_unit_nr + USE cp_para_types, ONLY: cp_para_env_type + USE cp_result_methods, ONLY: get_results + USE cp_result_types, ONLY: cp_result_type + USE dbcsr_api, ONLY: & + dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_init_p, dbcsr_p_type, dbcsr_set, & + dbcsr_type, dbcsr_type_no_symmetry, dbcsr_type_symmetric + USE input_constants, ONLY: current_orb_center_wannier,& + use_mom_ref_user + USE input_section_types, ONLY: section_vals_get_subs_vals,& + section_vals_type,& + section_vals_val_get + USE kinds, ONLY: default_path_length,& + default_string_length,& + dp + USE memory_utilities, ONLY: reallocate + USE message_passing, ONLY: mp_bcast + USE molecule_types, ONLY: molecule_type + USE moments_utils, ONLY: get_reference_point + USE particle_types, ONLY: particle_type + USE qs_environment_types, ONLY: get_qs_env,& + qs_environment_type + USE qs_kinetic, ONLY: build_kinetic_matrix + USE qs_ks_types, ONLY: qs_ks_env_type + USE qs_linres_types, ONLY: dcdr_env_type,& + linres_control_type + USE qs_loc_types, ONLY: get_qs_loc_env,& + localized_wfn_control_type,& + qs_loc_env_new_type + USE qs_mo_types, ONLY: get_mo_set,& + mo_set_p_type + USE qs_moments, ONLY: build_local_moment_matrix + USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type + USE qs_overlap, ONLY: build_overlap_matrix + USE string_utilities, ONLY: xstring +#include "./base/base_uses.f90" + + IMPLICIT NONE + + PRIVATE + PUBLIC :: dcdr_env_cleanup, dcdr_env_init, dcdr_print, & + get_loc_setting, shift_wannier_into_cell, & + dcdr_write_restart, dcdr_read_restart, & + multiply_localization + + CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_dcdr_utils' + + REAL(dp), DIMENSION(3, 3, 3), PARAMETER :: Levi_Civita = RESHAPE((/ & + 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, -1.0_dp, 0.0_dp, 1.0_dp, 0.0_dp, & + 0.0_dp, 0.0_dp, 1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, -1.0_dp, 0.0_dp, 0.0_dp, & + 0.0_dp, -1.0_dp, 0.0_dp, 1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp/), & + (/3, 3, 3/)) + +CONTAINS + +! ************************************************************************************************** +!> \brief Multiply (ao_matrix @ mo_coeff) and store the column icenter in res +!> \param ao_matrix ... +!> \param mo_coeff ... +!> \param work Working space +!> \param nmo ... +!> \param icenter ... +!> \param res ... +!> \author Edward Ditler +! ************************************************************************************************** + SUBROUTINE multiply_localization(ao_matrix, mo_coeff, work, nmo, icenter, res) + TYPE(dbcsr_type), INTENT(IN), POINTER :: ao_matrix + TYPE(cp_fm_type), INTENT(IN), POINTER :: mo_coeff, work + INTEGER, INTENT(IN) :: nmo, icenter + TYPE(cp_fm_type), INTENT(IN), POINTER :: res + + CHARACTER(LEN=*), PARAMETER :: routineN = 'multiply_localization', & + routineP = moduleN//':'//routineN + + INTEGER :: handle + + CALL timeset(routineN, handle) + + ! Multiply by the MO coefficients + CALL cp_dbcsr_sm_fm_multiply(ao_matrix, mo_coeff, work, ncol=nmo) + + ! Only keep the icenter-th column + CALL cp_fm_to_fm(work, res, 1, icenter, icenter) + + ! Reset the matrices + CALL cp_fm_set_all(work, 0.0_dp) + + CALL timestop(handle) + END SUBROUTINE multiply_localization + +! ************************************************************************************************** +!> \brief Copied from linres_read_restart +!> \param qs_env ... +!> \param linres_section ... +!> \param vec ... +!> \param lambda ... +!> \param beta ... +!> \param tag ... +!> \note Adapted from linres_read_restart (ED) +!> Would be nice not to crash but to start from zero if the present file doesn't match. +! ************************************************************************************************** + SUBROUTINE dcdr_read_restart(qs_env, linres_section, vec, lambda, beta, tag) + TYPE(qs_environment_type), POINTER :: qs_env + TYPE(section_vals_type), POINTER :: linres_section + TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: vec + INTEGER, INTENT(IN) :: lambda, beta + CHARACTER(LEN=*) :: tag + + CHARACTER(LEN=*), PARAMETER :: routineN = 'dcdr_read_restart' + + CHARACTER(LEN=default_path_length) :: filename + CHARACTER(LEN=default_string_length) :: my_middle + INTEGER :: beta_tmp, group, handle, i, i_block, ia, ie, iostat, iounit, ispin, j, & + lambda_tmp, max_block, n_rep_val, nao, nao_tmp, nmo, nmo_tmp, nspins, nspins_tmp, & + rst_unit, source + LOGICAL :: file_exists + REAL(KIND=dp), DIMENSION(:, :), POINTER :: vecbuffer + TYPE(cp_fm_type), POINTER :: mo_coeff + TYPE(cp_logger_type), POINTER :: logger + TYPE(cp_para_env_type), POINTER :: para_env + TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos + TYPE(section_vals_type), POINTER :: print_key + + file_exists = .FALSE. + + CALL timeset(routineN, handle) + + NULLIFY (mos, para_env, logger, print_key, vecbuffer) + logger => cp_get_default_logger() + + iounit = cp_print_key_unit_nr(logger, linres_section, & + "PRINT%PROGRAM_RUN_INFO", extension=".Log") + + CALL get_qs_env(qs_env=qs_env, & + para_env=para_env, & + mos=mos) + + nspins = SIZE(mos) + group = para_env%group + source = para_env%source !ionode??? + + rst_unit = -1 + IF (para_env%ionode) THEN + CALL section_vals_val_get(linres_section, "WFN_RESTART_FILE_NAME", & + n_rep_val=n_rep_val) + + CALL XSTRING(tag, ia, ie) + my_middle = "RESTART-"//tag(ia:ie)//TRIM("-")//TRIM(ADJUSTL(cp_to_string(beta))) & + //TRIM("-")//TRIM(ADJUSTL(cp_to_string(lambda))) + + IF (n_rep_val > 0) THEN + CALL section_vals_val_get(linres_section, "WFN_RESTART_FILE_NAME", c_val=filename) + CALL xstring(filename, ia, ie) + filename = filename(ia:ie)//TRIM(my_middle)//".lr" + ELSE + ! try to read from the filename that is generated automatically from the printkey + print_key => section_vals_get_subs_vals(linres_section, "PRINT%RESTART") + filename = cp_print_key_generate_filename(logger, print_key, & + extension=".lr", middle_name=TRIM(my_middle), my_local=.FALSE.) + END IF + INQUIRE (FILE=filename, exist=file_exists) + ! + ! open file + IF (file_exists) THEN + CALL open_file(file_name=TRIM(filename), & + file_action="READ", & + file_form="UNFORMATTED", & + file_position="REWIND", & + file_status="OLD", & + unit_number=rst_unit) + + IF (iounit > 0) WRITE (iounit, "(T2,A)") & + "LINRES| Reading response wavefunctions from the restart file <"//TRIM(ADJUSTL(filename))//">" + ELSE + IF (iounit > 0) WRITE (iounit, "(T2,A)") & + "LINRES| Restart file <"//TRIM(ADJUSTL(filename))//"> not found" + END IF + END IF + + CALL mp_bcast(file_exists, source, group) + + IF (file_exists) THEN + + CALL get_mo_set(mos(1)%mo_set, mo_coeff=mo_coeff) + CALL cp_fm_get_info(mo_coeff, nrow_global=nao, ncol_block=max_block) + + ALLOCATE (vecbuffer(nao, max_block)) + ! + ! read headers + IF (rst_unit > 0) READ (rst_unit, IOSTAT=iostat) lambda_tmp, beta_tmp, nspins_tmp, nao_tmp + CALL mp_bcast(iostat, source, group) + + CALL mp_bcast(beta_tmp, source, group) + CALL mp_bcast(lambda_tmp, source, group) + CALL mp_bcast(nspins_tmp, source, group) + CALL mp_bcast(nao_tmp, source, group) + + ! check that the number nao, nmo and nspins are + ! the same as in the current mos + IF (nspins_tmp .NE. nspins) THEN + CPABORT("nspins not consistent") + END IF + IF (nao_tmp .NE. nao) CPABORT("nao not consistent") + ! check that it's the right file + ! the same as in the current mos + IF (lambda_tmp .NE. lambda) CPABORT("lambda not consistent") + IF (beta_tmp .NE. beta) CPABORT("beta not consistent") + ! + DO ispin = 1, nspins + CALL get_mo_set(mos(ispin)%mo_set, mo_coeff=mo_coeff) + CALL cp_fm_get_info(mo_coeff, ncol_global=nmo) + ! + IF (rst_unit > 0) READ (rst_unit) nmo_tmp + CALL mp_bcast(nmo_tmp, source, group) + IF (nmo_tmp .NE. nmo) CPABORT("nmo not consistent") + ! + ! read the response + DO i = 1, nmo, MAX(max_block, 1) + i_block = MIN(max_block, nmo - i + 1) + DO j = 1, i_block + IF (rst_unit > 0) READ (rst_unit) vecbuffer(1:nao, j) + END DO + CALL mp_bcast(vecbuffer, source, group) + CALL cp_fm_set_submatrix(vec(ispin)%matrix, vecbuffer, 1, i, nao, i_block) + END DO + END DO + + IF (iostat /= 0) THEN + IF (iounit > 0) WRITE (iounit, "(T2,A)") & + "LINRES| Restart file <"//TRIM(ADJUSTL(filename))//"> not found" + END IF + + DEALLOCATE (vecbuffer) + + END IF + + IF (para_env%ionode) THEN + IF (file_exists) CALL close_file(unit_number=rst_unit) + END IF + + CALL timestop(handle) + + END SUBROUTINE dcdr_read_restart + +! ************************************************************************************************** +!> \brief Copied from linres_write_restart +!> \param qs_env ... +!> \param linres_section ... +!> \param vec ... +!> \param lambda ... +!> \param beta ... +!> \param tag ... +!> \note Adapted from linres_read_restart (ED) +!> Would be nice not to crash but to start from zero if the present file doesn't match. +! ************************************************************************************************** + SUBROUTINE dcdr_write_restart(qs_env, linres_section, vec, lambda, beta, tag) + TYPE(qs_environment_type), POINTER :: qs_env + TYPE(section_vals_type), POINTER :: linres_section + TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: vec + INTEGER, INTENT(IN) :: lambda, beta + CHARACTER(LEN=*) :: tag + + CHARACTER(LEN=*), PARAMETER :: routineN = 'dcdr_write_restart' + + CHARACTER(LEN=default_path_length) :: filename + CHARACTER(LEN=default_string_length) :: my_middle, my_pos, my_status + INTEGER :: handle, i, i_block, ia, ie, iounit, & + ispin, j, max_block, nao, nmo, nspins, & + rst_unit + REAL(KIND=dp), DIMENSION(:, :), POINTER :: vecbuffer + TYPE(cp_fm_type), POINTER :: mo_coeff + TYPE(cp_logger_type), POINTER :: logger + TYPE(cp_para_env_type), POINTER :: para_env + TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos + TYPE(section_vals_type), POINTER :: print_key + + NULLIFY (logger, mo_coeff, mos, para_env, print_key, vecbuffer) + + CALL timeset(routineN, handle) + + logger => cp_get_default_logger() + + IF (BTEST(cp_print_key_should_output(logger%iter_info, linres_section, "PRINT%RESTART", & + used_print_key=print_key), & + cp_p_file)) THEN + + iounit = cp_print_key_unit_nr(logger, linres_section, & + "PRINT%PROGRAM_RUN_INFO", extension=".Log") + + CALL get_qs_env(qs_env=qs_env, & + mos=mos, & + para_env=para_env) + + nspins = SIZE(mos) + + my_status = "REPLACE" + my_pos = "REWIND" + CALL XSTRING(tag, ia, ie) + my_middle = "RESTART-"//tag(ia:ie)//TRIM("-")//TRIM(ADJUSTL(cp_to_string(beta))) & + //TRIM("-")//TRIM(ADJUSTL(cp_to_string(lambda))) + rst_unit = cp_print_key_unit_nr(logger, linres_section, "PRINT%RESTART", & + extension=".lr", middle_name=TRIM(my_middle), file_status=TRIM(my_status), & + file_position=TRIM(my_pos), file_action="WRITE", file_form="UNFORMATTED") + + filename = cp_print_key_generate_filename(logger, print_key, & + extension=".lr", middle_name=TRIM(my_middle), my_local=.FALSE.) + + IF (iounit > 0) THEN + WRITE (UNIT=iounit, FMT="(T2,A)") & + "LINRES| Writing response functions to the restart file <"//TRIM(ADJUSTL(filename))//">" + END IF + + ! + ! write data to file + ! use the scalapack block size as a default for buffering columns + CALL get_mo_set(mos(1)%mo_set, mo_coeff=mo_coeff) + CALL cp_fm_get_info(mo_coeff, nrow_global=nao, ncol_block=max_block) + ALLOCATE (vecbuffer(nao, max_block)) + + IF (rst_unit > 0) WRITE (rst_unit) lambda, beta, nspins, nao + + DO ispin = 1, nspins + CALL cp_fm_get_info(vec(ispin)%matrix, ncol_global=nmo) + + IF (rst_unit > 0) WRITE (rst_unit) nmo + + DO i = 1, nmo, MAX(max_block, 1) + i_block = MIN(max_block, nmo - i + 1) + CALL cp_fm_get_submatrix(vec(ispin)%matrix, vecbuffer, 1, i, nao, i_block) + ! doing this in one write would increase efficiency, but breaks RESTART compatibility. + ! to old ones, and in cases where max_block is different between runs, as might happen during + ! restarts with a different number of CPUs + DO j = 1, i_block + IF (rst_unit > 0) WRITE (rst_unit) vecbuffer(1:nao, j) + END DO + END DO + END DO + + DEALLOCATE (vecbuffer) + + CALL cp_print_key_finished_output(rst_unit, logger, linres_section, & + "PRINT%RESTART") + END IF + + CALL timestop(handle) + + END SUBROUTINE dcdr_write_restart + +! ************************************************************************************************** +!> \brief Print the APT and sum rules +!> \param dcdr_env ... +!> \param qs_env ... +!> \author Edward Ditler +! ************************************************************************************************** + SUBROUTINE dcdr_print(dcdr_env, qs_env) + TYPE(dcdr_env_type) :: dcdr_env + TYPE(qs_environment_type), POINTER :: qs_env + + CHARACTER(len=*), PARAMETER :: routineN = 'dcdr_print', routineP = moduleN//':'//routineN + + CHARACTER(LEN=default_string_length) :: description + INTEGER :: alpha, beta, delta, gamma, i, k, l, & + lambda, natom, nsubset, output_unit + REAL(dp), DIMENSION(:, :, :), POINTER :: apt_el_dcdr, apt_nuc_dcdr, apt_total_dcdr + REAL(dp), DIMENSION(:, :, :, :), POINTER :: apt_center_dcdr, apt_subset_dcdr + REAL(kind=dp), DIMENSION(3, 3) :: sum_rule_0, sum_rule_1, sum_rule_2 + TYPE(cp_logger_type), POINTER :: logger + TYPE(cp_result_type), POINTER :: results + TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set + TYPE(particle_type), DIMENSION(:), POINTER :: particle_set + TYPE(section_vals_type), POINTER :: dcdr_section + + NULLIFY (logger) + + logger => cp_get_default_logger() + output_unit = cp_logger_get_default_io_unit(logger) + + dcdr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES%DCDR") + + NULLIFY (particle_set) + CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, molecule_set=molecule_set) + natom = SIZE(particle_set) + nsubset = SIZE(molecule_set) + + apt_el_dcdr => dcdr_env%apt_el_dcdr + apt_nuc_dcdr => dcdr_env%apt_nuc_dcdr + apt_total_dcdr => dcdr_env%apt_total_dcdr + apt_subset_dcdr => dcdr_env%apt_el_dcdr_per_subset + apt_center_dcdr => dcdr_env%apt_el_dcdr_per_center + + IF (dcdr_env%localized_psi0) THEN + IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A)") 'APT | Write the final apt matrix per atom per subset' + DO k = 1, natom + DO l = 1, nsubset + IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, *) 'APT | Subset', l + DO i = 1, 3 + IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A,I3,F15.6,F15.6,F15.6)") & + 'APT | apt_subset ', i, apt_subset_dcdr(i, :, k, l) + END DO + END DO + END DO + END IF + + IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A)") & + 'APT | Write the final apt matrix per atom (Position perturbation)' + DO l = 1, natom + IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A,I3,A,F15.6)") & + 'APT | Atom', l, ' - GAPT ', & + (apt_total_dcdr(1, 1, l) & + + apt_total_dcdr(2, 2, l) & + + apt_total_dcdr(3, 3, l))/3._dp + DO i = 1, 3 + IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A,F15.6,F15.6,F15.6)") "APT | ", apt_total_dcdr(i, :, l) + END DO + END DO + + IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A)") 'APT | Write the total apt matrix' + DO i = 1, 3 + IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, & + "(A,F15.6,F15.6,F15.6)") "APT | ", SUM(apt_total_dcdr(i, :, :), dim=2) + END DO + IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A)") 'APT | End Write the final apt matrix' + + ! Get the dipole + CALL get_qs_env(qs_env, results=results) + description = "[DIPOLE]" + CALL get_results(results=results, description=description, values=dcdr_env%dipole_pos(1:3)) + + ! Sum rules [for all alpha, beta] + sum_rule_0 = 0._dp + sum_rule_1 = 0._dp + sum_rule_2 = 0._dp + + DO alpha = 1, 3 + DO beta = 1, 3 + ! 0: sum_lambda apt(alpha, beta, lambda) + DO lambda = 1, natom + sum_rule_0(alpha, beta) = sum_rule_0(alpha, beta) & + + apt_total_dcdr(alpha, beta, lambda) + END DO + + ! 1: sum_gamma epsilon_(alpha beta gamma) mu_gamma + DO gamma = 1, 3 + sum_rule_1(alpha, beta) = sum_rule_1(alpha, beta) & + + Levi_Civita(alpha, beta, gamma)*dcdr_env%dipole_pos(gamma) + END DO + + ! 2: sum_(lambda gamma delta) R^lambda_gamma apt(delta, alpha, lambda) + DO lambda = 1, natom + DO gamma = 1, 3 + DO delta = 1, 3 + sum_rule_2(alpha, beta) = sum_rule_2(alpha, beta) & + + Levi_Civita(beta, gamma, delta) & + *particle_set(lambda)%r(gamma) & + *apt_total_dcdr(delta, alpha, lambda) + END DO + END DO + END DO + + END DO ! beta + END DO ! alpha + + IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A)") "APT | Position perturbation sum rules" + IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A,T18,A,T34,A,T49,A)") & + "APT |", "Total APT", "Dipole", "R * APT" + DO alpha = 1, 3 + DO beta = 1, 3 + IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, & + "(A,I3,I3,F15.6,F15.6,F15.6)") & + "APT | ", & + alpha, beta, & + sum_rule_0(alpha, beta), & + sum_rule_1(alpha, beta), & + sum_rule_2(alpha, beta) + END DO + END DO + + END SUBROUTINE dcdr_print + +! ************************************************************************************************** +!> \brief ... +!> \param r ... +!> \param cell ... +!> \param r_shifted ... +! ************************************************************************************************** + SUBROUTINE shift_wannier_into_cell(r, cell, r_shifted) + REAL(dp), DIMENSION(3), INTENT(in) :: r + TYPE(cell_type), INTENT(in), POINTER :: cell + REAL(dp), DIMENSION(3), INTENT(out) :: r_shifted + + INTEGER :: i + REAL(kind=dp), DIMENSION(3) :: abc + + ! Only orthorombic cell for now + CALL get_cell(cell, abc=abc) + + DO i = 1, 3 + IF (r(i) < 0._dp) THEN + r_shifted(i) = r(i) + abc(i) + ELSE IF (r(i) > abc(i)) THEN + r_shifted(i) = r(i) - abc(i) + ELSE + r_shifted(i) = r(i) + END IF + END DO + END SUBROUTINE shift_wannier_into_cell + +! ************************************************************************************************** +!> \brief ... +!> \param dcdr_env ... +!> \param qs_env ... +! ************************************************************************************************** + SUBROUTINE get_loc_setting(dcdr_env, qs_env) + TYPE(dcdr_env_type) :: dcdr_env + TYPE(qs_environment_type), POINTER :: qs_env + + CHARACTER(LEN=*), PARAMETER :: routineN = 'get_loc_setting', & + routineP = moduleN//':'//routineN + + INTEGER :: handle, is, istate, max_states, nmo, & + nmoloc, nstate, nstate_list(2) + INTEGER, ALLOCATABLE, DIMENSION(:, :) :: state_list + REAL(dp), DIMENSION(:, :), POINTER :: center_array + TYPE(linres_control_type), POINTER :: linres_control + TYPE(localized_wfn_control_type), POINTER :: localized_wfn_control + TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos + TYPE(qs_loc_env_new_type), POINTER :: qs_loc_env + TYPE(section_vals_type), POINTER :: dcdr_section + + CALL timeset(routineN, handle) + + CALL get_qs_env(qs_env=qs_env, & + linres_control=linres_control, & + mos=mos) + + ! Some checks + max_states = 0 + CALL get_mo_set(mo_set=mos(1)%mo_set, nmo=nmo) + max_states = MAX(max_states, nmo) + + ! check that the number of localized states is equal to the number of states + nmoloc = SIZE(linres_control%qs_loc_env%localized_wfn_control%centers_set(1)%array, 2) + IF (nmoloc .NE. nmo) THEN + CPABORT("The number of localized functions is not equal to the number of states.") + END IF + + ! which center for the orbitals shall we use + dcdr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES%DCDR") + CALL section_vals_val_get(dcdr_section, "ORBITAL_CENTER", i_val=dcdr_env%orb_center) + SELECT CASE (dcdr_env%orb_center) + CASE (current_orb_center_wannier) + dcdr_env%orb_center_name = "WANNIER" + CASE DEFAULT + CPABORT(" ") + END SELECT + + qs_loc_env => linres_control%qs_loc_env + CALL get_qs_loc_env(qs_loc_env, localized_wfn_control=localized_wfn_control) + + ALLOCATE (dcdr_env%centers_set(1)) + ALLOCATE (dcdr_env%center_list(1)) + ALLOCATE (state_list(max_states, 1)) + state_list(:, :) = HUGE(0) + nstate_list(:) = HUGE(0) + + ! Build the state_list + center_array => localized_wfn_control%centers_set(1)%array + nstate = 0 + DO istate = 1, SIZE(center_array, 2) + nstate = nstate + 1 + state_list(nstate, 1) = istate + END DO + nstate_list(1) = nstate + + ! clustering the states + nstate = nstate_list(1) + dcdr_env%nstates(1) = nstate + + ALLOCATE (dcdr_env%center_list(1)%array(2, nstate + 1)) + ALLOCATE (dcdr_env%centers_set(1)%array(3, nstate)) + dcdr_env%center_list(1)%array(:, :) = HUGE(0) + dcdr_env%centers_set(1)%array(:, :) = HUGE(0.0_dp) + + center_array => localized_wfn_control%centers_set(1)%array + + ! point to the psi0 centers + SELECT CASE (dcdr_env%orb_center) + CASE (current_orb_center_wannier) + ! use the wannier center as -center- + dcdr_env%nbr_center(1) = nstate + DO is = 1, nstate + istate = state_list(is, 1) + dcdr_env%centers_set(1)%array(1:3, is) = center_array(1:3, istate) + dcdr_env%center_list(1)%array(1, is) = is + dcdr_env%center_list(1)%array(2, is) = istate + END DO + dcdr_env%center_list(1)%array(1, nstate + 1) = nstate + 1 + + CASE DEFAULT + CPABORT("Unknown orbital center...") + END SELECT + + ! This doesn't do anything right? + ! ! hack to make centers_set%array closest to chosen reference + ! nsubset = SIZE(molecule_set) + ! natoms = SIZE(particle_set) + ! ALLOCATE (map_rc(3, natoms)) + ! ALLOCATE (mapm_rc(3, nsubset)) + + ! ! Get reference point for this molecule + ! NULLIFY (ref_point) + ! reference = section_get_ival(section_vals=qs_env%input, & + ! keyword_name="DFT%PRINT%MOMENTS%REFERENCE") + ! CALL section_vals_val_get(qs_env%input, "DFT%PRINT%MOMENTS%REF_POINT", r_vals=ref_point) + + ! DO i = 1, nsubset + ! first_atom = molecule_set(i)%first_atom + ! last_atom = molecule_set(i)%last_atom + + ! CALL get_reference_point(rcc, qs_env=qs_env, reference=reference, & + ! ref_point=ref_point, & + ! ifirst=first_atom, ilast=last_atom) + ! mapm_rc(:, i) = rcc + ! DO j = first_atom, last_atom + ! map_rc(:, j) = rcc + ! END DO + ! END DO + + ! DEALLOCATE (state_list) + ! DEALLOCATE (map_rc) + ! DEALLOCATE (mapm_rc) + + CALL timestop(handle) + END SUBROUTINE get_loc_setting + +! ************************************************************************************************** +!> \brief Initialize the dcdr environment +!> \param dcdr_env ... +!> \param qs_env ... +! ************************************************************************************************** + SUBROUTINE dcdr_env_init(dcdr_env, qs_env) + TYPE(dcdr_env_type) :: dcdr_env + TYPE(qs_environment_type), POINTER :: qs_env + + CHARACTER(LEN=*), PARAMETER :: routineN = 'dcdr_env_init', routineP = moduleN//':'//routineN + + INTEGER :: handle, homo, i, isize, j, jg, n_rep, & + nao, natom, nmo, nspins, nsubset, & + output_unit, reference, unit_number + INTEGER, DIMENSION(:), POINTER :: tmplist + LOGICAL :: explicit + REAL(KIND=dp), DIMENSION(:), POINTER :: ref_point + TYPE(cp_fm_type), POINTER :: buf, mo_coeff + TYPE(cp_logger_type), POINTER :: logger + TYPE(cp_para_env_type), POINTER :: para_env + TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s + TYPE(dft_control_type), POINTER :: dft_control + TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos + TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set + TYPE(neighbor_list_set_p_type), DIMENSION(:), & + POINTER :: sab_all, sab_orb + TYPE(particle_type), DIMENSION(:), POINTER :: particle_set + TYPE(qs_ks_env_type), POINTER :: ks_env + TYPE(section_vals_type), POINTER :: dcdr_section, loc_section, lr_section + + CALL timeset(routineN, handle) + + ! Set up the logger + NULLIFY (logger, loc_section, dcdr_section, lr_section) + logger => cp_get_default_logger() + loc_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES%LOCALIZE") + dcdr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES%DCDR") + dcdr_env%output_unit = cp_print_key_unit_nr(logger, dcdr_section, "PRINT%APT", & + extension=".data", middle_name="dcdr", log_filename=.FALSE., & + file_position="REWIND", file_status="REPLACE") + + 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") + unit_number = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", extension=".linresLog") + + IF (output_unit > 0) THEN + WRITE (output_unit, "(/,T20,A,/)") "*** Start DCDR calculation ***" + END IF + + NULLIFY (ks_env, dft_control, sab_orb, sab_all, particle_set, molecule_set, matrix_s, matrix_ks, mos, para_env) + CALL get_qs_env(qs_env=qs_env, & + ks_env=ks_env, & + dft_control=dft_control, & + sab_orb=sab_orb, & + sab_all=sab_all, & + particle_set=particle_set, & + molecule_set=molecule_set, & + matrix_s=matrix_s, & + matrix_ks=matrix_ks, & + mos=mos, & + para_env=para_env) + + natom = SIZE(particle_set) + nsubset = SIZE(molecule_set) + nspins = dft_control%nspins + IF (nspins == 2) THEN + CPABORT("NYE") + END IF + + CALL build_overlap_matrix(ks_env, matrix_s=dcdr_env%matrix_s, & + matrix_name="OVERLAP MATRIX", & + nderivative=1, & + basis_type_a="ORB", & + basis_type_b="ORB", & + sab_nl=sab_orb) + + CALL build_kinetic_matrix(ks_env, matrix_t=dcdr_env%matrix_t, & + matrix_name="KINETIC ENERGY MATRIX", & + basis_type="ORB", & + sab_nl=sab_orb, nderivative=1, & + eps_filter=dft_control%qs_control%eps_filter_matrix) + + ! Get inputs + CALL section_vals_val_get(dcdr_section, "DISTRIBUTED_ORIGIN", l_val=dcdr_env%distributed_origin) + CALL section_vals_val_get(loc_section, "_SECTION_PARAMETERS_", l_val=dcdr_env%localized_psi0) + CALL section_vals_val_get(lr_section, "RESTART", l_val=qs_env%linres_control%linres_restart) + dcdr_env%ref_point = 0._dp + + ! List of atoms + NULLIFY (tmplist) + isize = 0 + CALL section_vals_val_get(dcdr_section, "LIST_OF_ATOMS", n_rep_val=n_rep) + IF (n_rep == 0) THEN + ALLOCATE (dcdr_env%list_of_atoms(natom)) + DO jg = 1, natom + dcdr_env%list_of_atoms(jg) = jg + END DO + ELSE + DO jg = 1, n_rep + ALLOCATE (dcdr_env%list_of_atoms(isize)) + CALL section_vals_val_get(dcdr_section, "LIST_OF_ATOMS", i_rep_val=jg, i_vals=tmplist) + CALL reallocate(dcdr_env%list_of_atoms, 1, isize + SIZE(tmplist)) + dcdr_env%list_of_atoms(isize + 1:isize + SIZE(tmplist)) = tmplist + isize = SIZE(dcdr_env%list_of_atoms) + END DO + END IF + + ! Reference point + IF (dcdr_env%localized_psi0) THEN + ! Get the Wannier localized wave functions and centers + CALL get_loc_setting(dcdr_env, qs_env) + ELSE + ! Get the reference point from the input + CALL section_vals_val_get(dcdr_section, "REFERENCE", i_val=reference) + CALL section_vals_val_get(dcdr_section, "REFERENCE_POINT", explicit=explicit) + IF (explicit) THEN + CALL section_vals_val_get(dcdr_section, "REFERENCE_POINT", r_vals=ref_point) + ELSE + IF (reference == use_mom_ref_user) & + CPABORT("User-defined reference point should be given explicitly") + END IF + + CALL get_reference_point(rpoint=dcdr_env%ref_point, qs_env=qs_env, & + reference=reference, & + ref_point=ref_point) + END IF + + ! Helper matrix structs + NULLIFY (dcdr_env%aoao_fm_struct, & + dcdr_env%momo_fm_struct, & + dcdr_env%likemos_fm_struct, & + dcdr_env%homohomo_fm_struct) + CALL get_mo_set(mo_set=mos(1)%mo_set, mo_coeff=mo_coeff, & + nao=nao, nmo=nmo, homo=homo) + CALL cp_fm_struct_create(dcdr_env%aoao_fm_struct, nrow_global=nao, & + ncol_global=nao, para_env=para_env, & + context=mo_coeff%matrix_struct%context) + CALL cp_fm_struct_create(dcdr_env%momo_fm_struct, nrow_global=nmo, & + ncol_global=nmo, para_env=para_env, & + context=mo_coeff%matrix_struct%context) + CALL cp_fm_struct_create(dcdr_env%homohomo_fm_struct, nrow_global=homo, & + ncol_global=homo, para_env=para_env, & + context=mo_coeff%matrix_struct%context) + CALL cp_fm_struct_create(dcdr_env%likemos_fm_struct, & + template_fmstruct=mo_coeff%matrix_struct) + dcdr_env%nao = nao + dcdr_env%nmo = nmo + + ! Fields of reals + ALLOCATE (dcdr_env%deltaR(3, natom)) + ALLOCATE (dcdr_env%delta_basis_function(3, natom)) + ALLOCATE (dcdr_env%apt_el_dcdr(3, 3, natom)) + ALLOCATE (dcdr_env%apt_nuc_dcdr(3, 3, natom)) + ALLOCATE (dcdr_env%apt_total_dcdr(3, 3, natom)) + + dcdr_env%apt_el_dcdr = 0._dp + dcdr_env%apt_nuc_dcdr = 0._dp + dcdr_env%apt_total_dcdr = 0._dp + + dcdr_env%deltaR = 0.0_dp + dcdr_env%delta_basis_function = 0._dp + + ! Localization + IF (dcdr_env%localized_psi0) THEN + ALLOCATE (dcdr_env%apt_el_dcdr_per_center(3, 3, natom, dcdr_env%nbr_center(1))) + ALLOCATE (dcdr_env%apt_el_dcdr_per_subset(3, 3, natom, nsubset)) + ALLOCATE (dcdr_env%apt_subset(3, 3, natom, nsubset)) + dcdr_env%apt_el_dcdr_per_center = 0._dp + dcdr_env%apt_el_dcdr_per_subset = 0._dp + dcdr_env%apt_subset = 0.0_dp + END IF + + ! Full matrices + ALLOCATE (dcdr_env%mo_coeff(nspins)) + ALLOCATE (dcdr_env%dCR(nspins)) + ALLOCATE (dcdr_env%dCR_prime(nspins)) + ALLOCATE (dcdr_env%chc(nspins)) + ALLOCATE (dcdr_env%op_dR(nspins)) + + NULLIFY (dcdr_env%mo_coeff(1)%matrix) + NULLIFY (dcdr_env%dCR(1)%matrix) + NULLIFY (dcdr_env%dCR_prime(1)%matrix) + NULLIFY (dcdr_env%chc(1)%matrix) + NULLIFY (dcdr_env%op_dR(1)%matrix) + + CALL cp_fm_create(dcdr_env%dCR(1)%matrix, dcdr_env%likemos_fm_struct) + CALL cp_fm_create(dcdr_env%dCR_prime(1)%matrix, dcdr_env%likemos_fm_struct) + CALL cp_fm_create(dcdr_env%mo_coeff(1)%matrix, dcdr_env%likemos_fm_struct) + CALL cp_fm_create(dcdr_env%chc(1)%matrix, dcdr_env%momo_fm_struct) + CALL cp_fm_create(dcdr_env%op_dR(1)%matrix, dcdr_env%likemos_fm_struct) + + CALL cp_fm_to_fm(mo_coeff, dcdr_env%mo_coeff(1)%matrix) + + ! DBCSR matrices + NULLIFY (dcdr_env%hamiltonian1) + NULLIFY (dcdr_env%moments) + NULLIFY (dcdr_env%matrix_difdip) + NULLIFY (dcdr_env%matrix_core_charge_1) + NULLIFY (dcdr_env%matrix_nosym_temp) + NULLIFY (dcdr_env%matrix_s1) + NULLIFY (dcdr_env%matrix_t1) + NULLIFY (dcdr_env%matrix_apply_op_constant) + NULLIFY (dcdr_env%matrix_d_vhxc_dR) + NULLIFY (dcdr_env%matrix_vhxc_perturbed_basis) + NULLIFY (dcdr_env%matrix_hc) + NULLIFY (dcdr_env%matrix_ppnl_1) + CALL dbcsr_allocate_matrix_set(dcdr_env%hamiltonian1, nspins) + CALL dbcsr_allocate_matrix_set(dcdr_env%moments, 3) + CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_difdip, 3, 3) + CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_core_charge_1, 3) + CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_nosym_temp, 3) + CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_s1, 4) + CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_t1, 4) + CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_apply_op_constant, 1) + CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_d_vhxc_dR, 3) + CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_vhxc_perturbed_basis, nspins, 6) + CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_hc, 3) + CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_ppnl_1, 3) + + CALL dbcsr_init_p(dcdr_env%perturbed_dm_correction) + CALL dbcsr_copy(dcdr_env%perturbed_dm_correction, matrix_ks(1)%matrix) + + ! temporary no_symmetry matrix: + DO i = 1, 3 + CALL dbcsr_init_p(dcdr_env%matrix_nosym_temp(i)%matrix) + CALL dbcsr_create(dcdr_env%matrix_nosym_temp(i)%matrix, template=matrix_ks(1)%matrix, & + matrix_type=dbcsr_type_no_symmetry) + CALL cp_dbcsr_alloc_block_from_nbl(dcdr_env%matrix_nosym_temp(i)%matrix, sab_all) + CALL dbcsr_set(dcdr_env%matrix_nosym_temp(i)%matrix, 0._dp) + + END DO + + ! moments carry the result of build_local_moment_matrix + DO i = 1, 3 + CALL dbcsr_init_p(dcdr_env%moments(i)%matrix) + CALL dbcsr_copy(dcdr_env%moments(i)%matrix, matrix_ks(1)%matrix, "dcdr_env%moments") + CALL dbcsr_set(dcdr_env%moments(i)%matrix, 0.0_dp) + END DO + CALL build_local_moment_matrix(qs_env, dcdr_env%moments, 1, ref_point=[0._dp, 0._dp, 0._dp]) + + DO i = 1, 3 + DO j = 1, 3 + CALL dbcsr_init_p(dcdr_env%matrix_difdip(i, j)%matrix) + CALL dbcsr_copy(dcdr_env%matrix_difdip(i, j)%matrix, dcdr_env%matrix_nosym_temp(1)%matrix) + CALL dbcsr_set(dcdr_env%matrix_difdip(i, j)%matrix, 0.0_dp) + END DO + END DO + + CALL dbcsr_init_p(dcdr_env%hamiltonian1(1)%matrix) + + ! overlap/kinetic matrix: s(1) normal overlap matrix; + ! s(2:4) derivatives wrt. nuclear coordinates + CALL dbcsr_init_p(dcdr_env%matrix_s1(1)%matrix) + CALL dbcsr_init_p(dcdr_env%matrix_t1(1)%matrix) + + CALL dbcsr_copy(dcdr_env%matrix_s1(1)%matrix, matrix_s(1)%matrix) + CALL dbcsr_copy(dcdr_env%matrix_t1(1)%matrix, dcdr_env%matrix_t(1)%matrix) + + DO i = 2, 4 + CALL dbcsr_init_p(dcdr_env%matrix_s1(i)%matrix) + CALL dbcsr_copy(dcdr_env%matrix_s1(i)%matrix, dcdr_env%matrix_nosym_temp(1)%matrix) + + CALL dbcsr_init_p(dcdr_env%matrix_t1(i)%matrix) + CALL dbcsr_copy(dcdr_env%matrix_t1(i)%matrix, dcdr_env%matrix_nosym_temp(1)%matrix) + END DO + + CALL dbcsr_init_p(dcdr_env%matrix_apply_op_constant(1)%matrix) + CALL dbcsr_copy(dcdr_env%matrix_apply_op_constant(1)%matrix, dcdr_env%matrix_s1(1)%matrix) + + ! j=1...3: derivative wrt nucleus A, 4...6: wrt nucleus B + DO j = 1, 6 + CALL dbcsr_init_p(dcdr_env%matrix_vhxc_perturbed_basis(1, j)%matrix) + CALL dbcsr_copy(dcdr_env%matrix_vhxc_perturbed_basis(1, j)%matrix, dcdr_env%matrix_s1(1)%matrix) + END DO + + DO i = 1, 3 + CALL dbcsr_init_p(dcdr_env%matrix_hc(i)%matrix) + CALL dbcsr_create(dcdr_env%matrix_hc(i)%matrix, template=matrix_ks(1)%matrix, & + matrix_type=dbcsr_type_symmetric) + CALL cp_dbcsr_alloc_block_from_nbl(dcdr_env%matrix_hc(i)%matrix, sab_orb) + CALL dbcsr_set(dcdr_env%matrix_hc(i)%matrix, 0.0_dp) + END DO + + DO i = 1, 3 + CALL dbcsr_init_p(dcdr_env%matrix_ppnl_1(i)%matrix) + CALL dbcsr_create(dcdr_env%matrix_ppnl_1(i)%matrix, template=matrix_ks(1)%matrix, & + matrix_type=dbcsr_type_symmetric) + CALL cp_dbcsr_alloc_block_from_nbl(dcdr_env%matrix_ppnl_1(i)%matrix, sab_orb) + CALL dbcsr_set(dcdr_env%matrix_ppnl_1(i)%matrix, 0.0_dp) + END DO + + DO i = 1, 3 + CALL dbcsr_init_p(dcdr_env%matrix_d_vhxc_dR(i)%matrix) + CALL dbcsr_copy(dcdr_env%matrix_d_vhxc_dR(i)%matrix, dcdr_env%matrix_s1(1)%matrix) + + CALL dbcsr_init_p(dcdr_env%matrix_core_charge_1(i)%matrix) + CALL dbcsr_copy(dcdr_env%matrix_core_charge_1(i)%matrix, dcdr_env%matrix_s1(1)%matrix) + CALL dbcsr_set(dcdr_env%matrix_core_charge_1(i)%matrix, 0.0_dp) + END DO + + ! CHC + CALL get_mo_set(mo_set=mos(1)%mo_set, mo_coeff=mo_coeff, nao=nao, nmo=nmo) + CALL cp_fm_create(buf, dcdr_env%likemos_fm_struct) + + CALL cp_dbcsr_sm_fm_multiply(matrix_ks(1)%matrix, mo_coeff, buf, nmo) + ! chc = mo * matrix_ks * mo + CALL cp_fm_get_info(buf) + CALL cp_fm_get_info(dcdr_env%chc(1)%matrix) + CALL cp_gemm('T', 'N', nmo, nmo, nao, & + 1.0_dp, mo_coeff, buf, & + 0.0_dp, dcdr_env%chc(1)%matrix) + + CALL cp_fm_release(buf) + + CALL cp_print_key_finished_output(output_unit, logger, lr_section, & + "PRINT%PROGRAM_RUN_INFO") + + IF (dcdr_env%output_unit > 0) & + WRITE (dcdr_env%output_unit, *) 'List of atoms: ', dcdr_env%list_of_atoms + + CALL timestop(handle) + END SUBROUTINE dcdr_env_init + +! ************************************************************************************************** +!> \brief Deallocate the dcdr environment +!> \param qs_env ... +!> \param dcdr_env ... +! ************************************************************************************************** + SUBROUTINE dcdr_env_cleanup(qs_env, dcdr_env) + + TYPE(qs_environment_type), POINTER :: qs_env + TYPE(dcdr_env_type) :: dcdr_env + + CHARACTER(LEN=*), PARAMETER :: routineN = 'dcdr_env_cleanup', & + routineP = moduleN//':'//routineN + + TYPE(cp_logger_type), POINTER :: logger + TYPE(section_vals_type), POINTER :: dcdr_section + + ! Destroy the logger + logger => cp_get_default_logger() + dcdr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES%DCDR") + CALL cp_print_key_finished_output(dcdr_env%output_unit, logger, dcdr_section, "PRINT%APT") + + DEALLOCATE (dcdr_env%list_of_atoms) + + CALL cp_fm_struct_release(dcdr_env%aoao_fm_struct) + CALL cp_fm_struct_release(dcdr_env%momo_fm_struct) + CALL cp_fm_struct_release(dcdr_env%homohomo_fm_struct) + CALL cp_fm_struct_release(dcdr_env%likemos_fm_struct) + + DEALLOCATE (dcdr_env%deltar) + DEALLOCATE (dcdr_env%delta_basis_function) + + IF (dcdr_env%localized_psi0) THEN + ! DEALLOCATE (dcdr_env%psi0_order) + DEALLOCATE (dcdr_env%centers_set(1)%array) + DEALLOCATE (dcdr_env%center_list(1)%array) + DEALLOCATE (dcdr_env%centers_set) + DEALLOCATE (dcdr_env%center_list) + DEALLOCATE (dcdr_env%apt_subset) + END IF + + DEALLOCATE (dcdr_env%apt_el_dcdr) + DEALLOCATE (dcdr_env%apt_nuc_dcdr) + DEALLOCATE (dcdr_env%apt_total_dcdr) + IF (dcdr_env%localized_psi0) THEN + DEALLOCATE (dcdr_env%apt_el_dcdr_per_center) + DEALLOCATE (dcdr_env%apt_el_dcdr_per_subset) + END IF + + ! Full matrices + CALL cp_fm_release(dcdr_env%dCR(1)%matrix) + CALL cp_fm_release(dcdr_env%dCR_prime(1)%matrix) + CALL cp_fm_release(dcdr_env%mo_coeff(1)%matrix) + CALL cp_fm_release(dcdr_env%chc(1)%matrix) + CALL cp_fm_release(dcdr_env%op_dR(1)%matrix) + + DEALLOCATE (dcdr_env%op_dR) + DEALLOCATE (dcdr_env%dCR) + DEALLOCATE (dcdr_env%dCR_prime) + DEALLOCATE (dcdr_env%mo_coeff) + DEALLOCATE (dcdr_env%chc) + + ! DBCSR matrices + CALL dbcsr_deallocate_matrix(dcdr_env%perturbed_dm_correction) + CALL dbcsr_deallocate_matrix_set(dcdr_env%hamiltonian1) + CALL dbcsr_deallocate_matrix_set(dcdr_env%moments) + CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_difdip) + CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_core_charge_1) + CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_nosym_temp) + CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_s) + CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_t) + CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_s1) + CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_t1) + CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_apply_op_constant) + CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_d_vhxc_dR) + CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_vhxc_perturbed_basis) + CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_hc) + CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_ppnl_1) + + END SUBROUTINE dcdr_env_cleanup + +END MODULE qs_dcdr_utils diff --git a/src/qs_integrate_potential.F b/src/qs_integrate_potential.F index a6d1ae9092..a13da2a48d 100644 --- a/src/qs_integrate_potential.F +++ b/src/qs_integrate_potential.F @@ -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, & diff --git a/src/qs_integrate_potential_product.F b/src/qs_integrate_potential_product.F index e9310a0b32..5ef5383209 100644 --- a/src/qs_integrate_potential_product.F +++ b/src/qs_integrate_potential_product.F @@ -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 ... diff --git a/src/qs_kinetic.F b/src/qs_kinetic.F index 89ad5e2d0d..00fb3561bb 100644 --- a/src/qs_kinetic.F +++ b/src/qs_kinetic.F @@ -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 diff --git a/src/qs_linres_module.F b/src/qs_linres_module.F index 99bfa0fad1..0f1967cf51 100644 --- a/src/qs_linres_module.F +++ b/src/qs_linres_module.F @@ -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) diff --git a/src/qs_linres_types.F b/src/qs_linres_types.F index 890e92eb64..7ea005e9ea 100644 --- a/src/qs_linres_types.F +++ b/src/qs_linres_types.F @@ -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 diff --git a/src/qs_moments.F b/src/qs_moments.F index a90c9b82eb..96bfd7e8e5 100644 --- a/src/qs_moments.F +++ b/src/qs_moments.F @@ -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 + +!> 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 diff --git a/src/qs_operators_ao.F b/src/qs_operators_ao.F index efaf918ef0..18256fedbe 100644 --- a/src/qs_operators_ao.F +++ b/src/qs_operators_ao.F @@ -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 diff --git a/tests/QS/regtest-dcdr/TEST_FILES b/tests/QS/regtest-dcdr/TEST_FILES new file mode 100644 index 0000000000..5d471bd374 --- /dev/null +++ b/tests/QS/regtest-dcdr/TEST_FILES @@ -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 diff --git a/tests/QS/regtest-dcdr/h2o_apt.inp b/tests/QS/regtest-dcdr/h2o_apt.inp new file mode 100644 index 0000000000..6e427682c5 --- /dev/null +++ b/tests/QS/regtest-dcdr/h2o_apt.inp @@ -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 diff --git a/tests/QS/regtest-dcdr/h2o_apt_loc.inp b/tests/QS/regtest-dcdr/h2o_apt_loc.inp new file mode 100644 index 0000000000..f7f63f8284 --- /dev/null +++ b/tests/QS/regtest-dcdr/h2o_apt_loc.inp @@ -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 diff --git a/tests/QS/regtest-dcdr/h2o_apt_pbc.inp b/tests/QS/regtest-dcdr/h2o_apt_pbc.inp new file mode 100644 index 0000000000..147b039825 --- /dev/null +++ b/tests/QS/regtest-dcdr/h2o_apt_pbc.inp @@ -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 diff --git a/tests/QS/regtest-dcdr/h2o_apt_pbc_loc.inp b/tests/QS/regtest-dcdr/h2o_apt_pbc_loc.inp new file mode 100644 index 0000000000..2e43afeef9 --- /dev/null +++ b/tests/QS/regtest-dcdr/h2o_apt_pbc_loc.inp @@ -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 diff --git a/tests/TEST_DIRS b/tests/TEST_DIRS index 8ffc9a9308..4e159d59c5 100644 --- a/tests/TEST_DIRS +++ b/tests/TEST_DIRS @@ -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 diff --git a/tests/TEST_TYPES b/tests/TEST_TYPES index 8f48bf05dd..8122aa36e6 100644 --- a/tests/TEST_TYPES +++ b/tests/TEST_TYPES @@ -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