Add analytical derivatives of MO coefficients wrt nuclear coordinates

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

src/aobasis/ai_oneelectron.F

@@ -78,13 +78,21 @@ MODULE ai_oneelectron
 !> \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 @@ SUBROUTINE os_3center(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
          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 @@ SUBROUTINE os_3center(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
          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 @@ SUBROUTINE os_3center(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
                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 @@ SUBROUTINE os_3center(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
                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

src/aobasis/ai_overlap_ppl.F

@@ -85,13 +85,21 @@ MODULE ai_overlap_ppl
 !> \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 @@ SUBROUTINE ppl_integral(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
          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 @@ SUBROUTINE ppl_integral(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
          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 @@ SUBROUTINE ppl_integral(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
 
       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)
 

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 @@ SUBROUTINE add_all_references()
                          "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