semi_empirical_mpole_methods.F

Go to the documentation of this file.

!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2026 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!
 
! **************************************************************************************************
!> \brief Setup and Methods for semi-empirical multipole types
!> \author Teodoro Laino [tlaino] - 08.2008 Zurich University
! **************************************************************************************************
MODULE semi_empirical_mpole_methods
 
   USE input_constants,                 ONLY: do_method_pnnl
   USE kinds,                           ONLY: dp
   USE mathconstants,                   ONLY: sqrt3
   USE semi_empirical_int_arrays,       ONLY: alm,&
                                              indexa,&
                                              indexb,&
                                              se_map_alm
   USE semi_empirical_mpole_types,      ONLY: nddo_mpole_create,&
                                              nddo_mpole_release,&
                                              nddo_mpole_type,&
                                              semi_empirical_mpole_p_create,&
                                              semi_empirical_mpole_p_type,&
                                              semi_empirical_mpole_type
   USE semi_empirical_par_utils,        ONLY: amn_l
   USE semi_empirical_types,            ONLY: semi_empirical_type
#include "./base/base_uses.f90"
 
   IMPLICIT NONE
 
   PRIVATE
 
! *** Global parameters ***
   LOGICAL, PARAMETER, PRIVATE          :: debug_this_module = .false.
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'semi_empirical_mpole_methods'
 
   PUBLIC :: semi_empirical_mpole_p_setup, &
             nddo_mpole_setup, &
             quadrupole_sph_to_cart
 
CONTAINS
 
! **************************************************************************************************
!> \brief Setup semi-empirical mpole type
!>        This function setup for each semi-empirical type a structure containing
!>        the multipolar expansion for all possible combination on-site of atomic
!>        orbitals ( \mu \nu |
!> \param mpoles ...
!> \param se_parameter ...
!> \param method ...
!> \date   09.2008
!> \author Teodoro Laino [tlaino] - University of Zurich
! **************************************************************************************************
   SUBROUTINE semi_empirical_mpole_p_setup(mpoles, se_parameter, method)
      TYPE(semi_empirical_mpole_p_type), DIMENSION(:), &
         POINTER                                         :: mpoles
      TYPE(semi_empirical_type), POINTER                 :: se_parameter
      INTEGER, INTENT(IN)                                :: method
 
      CHARACTER(LEN=3), DIMENSION(9), PARAMETER :: &
         label_print_orb = ["  s", " px", " py", " pz", "dx2", "dzx", "dz2", "dzy", "dxy"]
      INTEGER, DIMENSION(9), PARAMETER :: loc_index = [1, 2, 2, 2, 3, 3, 3, 3, 3]
 
      INTEGER                                            :: a, b, i, ind1, ind2, j, k, k1, k2, mu, &
                                                            natorb, ndim, nr
      REAL(kind=dp)                                      :: dlm, tmp, wp, ws, zb, zp, zs, zt
      REAL(kind=dp), DIMENSION(3, 3, 45)                 :: m2
      REAL(kind=dp), DIMENSION(3, 45)                    :: m1
      REAL(kind=dp), DIMENSION(45)                       :: m0
      REAL(kind=dp), DIMENSION(6, 0:2)                   :: amn
      TYPE(semi_empirical_mpole_type), POINTER           :: mpole
 
      cpassert(.NOT. ASSOCIATED(mpoles))
      ! If there are atomic orbitals proceed with the expansion in multipoles
      natorb = se_parameter%natorb
      IF (natorb /= 0) THEN
         ndim = natorb*(natorb + 1)/2
         CALL semi_empirical_mpole_p_create(mpoles, ndim)
 
         ! Select method for multipolar expansion
 
         ! Fill in information on multipole expansion due to atomic orbitals charge
         ! distribution
         NULLIFY (mpole)
         CALL amn_l(se_parameter, amn)
         DO i = 1, natorb
            DO j = 1, i
               ind1 = indexa(se_map_alm(i), se_map_alm(j))
               ind2 = indexb(i, j)
               ! the order in the mpoles structure is like the standard one for the
               ! integrals: s px py pz dx2-y2 dzx dz2 dzy dxy (lower triangular)
               ! which differs from the order of the Hamiltonian in CP2K. But I
               ! preferred to keep this order for consistency with the integrals
               mpole => mpoles(ind2)%mpole
               mpole%indi = i
               mpole%indj = j
               a = loc_index(i)
               b = loc_index(j)
               mpole%c = huge(0.0_dp)
               mpole%d = huge(0.0_dp)
               mpole%qs = huge(0.0_dp)
               mpole%qc = huge(0.0_dp)
 
               ! Charge
               IF (alm(ind1, 0, 0) /= 0.0_dp) THEN
                  dlm = 1.0_dp/sqrt(real((2*0 + 1), kind=dp))
                  tmp = -dlm*amn(indexb(a, b), 0)
                  mpole%c = tmp*alm(ind1, 0, 0)
                  mpole%task(1) = .true.
               END IF
 
               ! Dipole
               IF (any(alm(ind1, 1, -1:1) /= 0.0_dp)) THEN
                  dlm = 1.0_dp/sqrt(real((2*1 + 1), kind=dp))
                  tmp = -dlm*amn(indexb(a, b), 1)
                  mpole%d(1) = tmp*alm(ind1, 1, 1)
                  mpole%d(2) = tmp*alm(ind1, 1, -1)
                  mpole%d(3) = tmp*alm(ind1, 1, 0)
                  mpole%task(2) = .true.
               END IF
 
               ! Quadrupole
               IF (any(alm(ind1, 2, -2:2) /= 0.0_dp)) THEN
                  dlm = 1.0_dp/sqrt(real((2*2 + 1), kind=dp))
                  tmp = -dlm*amn(indexb(a, b), 2)
 
                  ! Spherical components
                  mpole%qs(1) = tmp*alm(ind1, 2, 0) ! d3z2-r2
                  mpole%qs(2) = tmp*alm(ind1, 2, 1) ! dzx
                  mpole%qs(3) = tmp*alm(ind1, 2, -1) ! dzy
                  mpole%qs(4) = tmp*alm(ind1, 2, 2) ! dx2-y2
                  mpole%qs(5) = tmp*alm(ind1, 2, -2) ! dxy
 
                  ! Convert into cartesian components
                  CALL quadrupole_sph_to_cart(mpole%qc, mpole%qs)
                  mpole%task(3) = .true.
               END IF
 
               IF (debug_this_module) THEN
                  WRITE (*, '(A,2I6,A)') "Orbitals ", i, j, &
                     " ("//label_print_orb(i)//","//label_print_orb(j)//")"
                  IF (mpole%task(1)) WRITE (*, '(9F12.6)') mpole%c
                  IF (mpole%task(2)) WRITE (*, '(9F12.6)') mpole%d
                  IF (mpole%task(3)) WRITE (*, '(9F12.6)') mpole%qc
                  WRITE (*, *)
               END IF
            END DO
         END DO
 
         IF (method == do_method_pnnl) THEN
            ! No d-function for Schenter type integrals
            cpassert(natorb <= 4)
 
            m0 = 0.0_dp
            m1 = 0.0_dp
            m2 = 0.0_dp
 
            DO mu = 1, se_parameter%natorb
               m0(indexb(mu, mu)) = 1.0_dp
            END DO
 
            zs = se_parameter%sto_exponents(0)
            zp = se_parameter%sto_exponents(1)
            nr = se_parameter%nr
 
            ws = real((2*nr + 2)*(2*nr + 1), dp)/(24.0_dp*zs**2)
            DO k = 1, 3
               m2(k, k, indexb(1, 1)) = ws
            END DO
 
            IF (zp > 0._dp) THEN
               zt = sqrt(zs*zp)
               zb = 0.5_dp*(zs + zp)
               DO k = 1, 3
                  m1(k, indexb(1, 1 + k)) = (zt/zb)**(2*nr + 1)*real(2*nr + 1, dp)/(2.0_dp*zb*sqrt3)
               END DO
 
               wp = real((2*nr + 2)*(2*nr + 1), dp)/(40.0_dp*zp**2)
               DO k1 = 1, 3
                  DO k2 = 1, 3
                     IF (k1 == k2) THEN
                        m2(k2, k2, indexb(1 + k1, 1 + k1)) = 3.0_dp*wp
                     ELSE
                        m2(k2, k2, indexb(1 + k1, 1 + k1)) = wp
                     END IF
                  END DO
               END DO
               m2(1, 2, indexb(1 + 1, 1 + 2)) = wp
               m2(2, 1, indexb(1 + 1, 1 + 2)) = wp
               m2(2, 3, indexb(1 + 2, 1 + 3)) = wp
               m2(3, 2, indexb(1 + 2, 1 + 3)) = wp
               m2(3, 1, indexb(1 + 3, 1 + 1)) = wp
               m2(1, 3, indexb(1 + 3, 1 + 1)) = wp
            END IF
 
            DO i = 1, natorb
               DO j = 1, i
                  ind1 = indexa(se_map_alm(i), se_map_alm(j))
                  ind2 = indexb(i, j)
                  mpole => mpoles(ind2)%mpole
                  mpole%indi = i
                  mpole%indj = j
                  ! Charge
                  mpole%cs = -m0(indexb(i, j))
                  ! Dipole
                  mpole%ds = -m1(1:3, indexb(i, j))
                  ! Quadrupole
                  mpole%qq = -3._dp*m2(1:3, 1:3, indexb(i, j))
                  IF (debug_this_module) THEN
                     WRITE (*, '(A,2I6,A)') "Orbitals ", i, j, &
                        " ("//label_print_orb(i)//","//label_print_orb(j)//")"
                     WRITE (*, '(9F12.6)') mpole%cs
                     WRITE (*, '(9F12.6)') mpole%ds
                     WRITE (*, '(9F12.6)') mpole%qq
                     WRITE (*, *)
                  END IF
               END DO
            END DO
         ELSE
            mpole%cs = mpole%c
            mpole%ds = mpole%d
            mpole%qq = mpole%qc
         END IF
      END IF
 
   END SUBROUTINE semi_empirical_mpole_p_setup
 
! **************************************************************************************************
!> \brief  Transforms the quadrupole components from sphericals to cartesians
!> \param qcart ...
!> \param qsph ...
!> \date   09.2008
!> \author Teodoro Laino [tlaino] - University of Zurich
! **************************************************************************************************
   SUBROUTINE quadrupole_sph_to_cart(qcart, qsph)
      REAL(kind=dp), DIMENSION(3, 3), INTENT(OUT)        :: qcart
      REAL(kind=dp), DIMENSION(5), INTENT(IN)            :: qsph
 
! Notation
!          qs(1) - d3z2-r2
!          qs(2) - dzx
!          qs(3) - dzy
!          qs(4) - dx2-y2
!          qs(5) - dxy
! Cartesian components
 
      qcart(1, 1) = (qsph(4) - qsph(1)/sqrt(3.0_dp))*sqrt(3.0_dp)/2.0_dp
      qcart(2, 1) = qsph(5)*sqrt(3.0_dp)/2.0_dp
      qcart(3, 1) = qsph(2)*sqrt(3.0_dp)/2.0_dp
      qcart(2, 2) = -(qsph(4) + qsph(1)/sqrt(3.0_dp))*sqrt(3.0_dp)/2.0_dp
      qcart(3, 2) = qsph(3)*sqrt(3.0_dp)/2.0_dp
      qcart(3, 3) = qsph(1)
      ! Symmetrize tensor
      qcart(1, 2) = qcart(2, 1)
      qcart(1, 3) = qcart(3, 1)
      qcart(2, 3) = qcart(3, 2)
 
   END SUBROUTINE quadrupole_sph_to_cart
 
! **************************************************************************************************
!> \brief Setup NDDO multipole type
!> \param nddo_mpole ...
!> \param natom ...
!> \date   09.2008
!> \author Teodoro Laino [tlaino] - University of Zurich
! **************************************************************************************************
   SUBROUTINE nddo_mpole_setup(nddo_mpole, natom)
      TYPE(nddo_mpole_type), POINTER                     :: nddo_mpole
      INTEGER, INTENT(IN)                                :: natom
 
      CHARACTER(len=*), PARAMETER                        :: routinen = 'nddo_mpole_setup'
 
      INTEGER                                            :: handle
 
      CALL timeset(routinen, handle)
 
      IF (ASSOCIATED(nddo_mpole)) THEN
         CALL nddo_mpole_release(nddo_mpole)
      END IF
      CALL nddo_mpole_create(nddo_mpole)
      ! Allocate Global Arrays
      ALLOCATE (nddo_mpole%charge(natom))
      ALLOCATE (nddo_mpole%dipole(3, natom))
      ALLOCATE (nddo_mpole%quadrupole(3, 3, natom))
 
      ALLOCATE (nddo_mpole%efield0(natom))
      ALLOCATE (nddo_mpole%efield1(3, natom))
      ALLOCATE (nddo_mpole%efield2(9, natom))
 
      CALL timestop(handle)
 
   END SUBROUTINE nddo_mpole_setup
 
END MODULE semi_empirical_mpole_methods