qmmm_per_elpot.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 Setting up the potential for QM/MM periodic boundary conditions calculations
!> \par History
!>      07.2005 created [tlaino]
!> \author Teodoro Laino
! **************************************************************************************************
MODULE qmmm_per_elpot
   USE ao_util,                         ONLY: exp_radius
   USE cell_types,                      ONLY: cell_type
   USE cp_log_handling,                 ONLY: cp_get_default_logger,&
                                              cp_logger_get_default_io_unit,&
                                              cp_logger_type
   USE cp_output_handling,              ONLY: cp_print_key_finished_output,&
                                              cp_print_key_unit_nr
   USE ewald_environment_types,         ONLY: ewald_env_create,&
                                              ewald_env_get,&
                                              ewald_env_set,&
                                              ewald_environment_type,&
                                              read_ewald_section
   USE ewald_pw_types,                  ONLY: ewald_pw_create,&
                                              ewald_pw_type
   USE ewald_spline_util,               ONLY: setup_ewald_spline
   USE input_constants,                 ONLY: do_qmmm_coulomb,&
                                              do_qmmm_gauss,&
                                              do_qmmm_pcharge,&
                                              do_qmmm_swave
   USE input_section_types,             ONLY: section_vals_get_subs_vals,&
                                              section_vals_type,&
                                              section_vals_val_get
   USE kinds,                           ONLY: dp
   USE mathconstants,                   ONLY: pi
   USE message_passing,                 ONLY: mp_para_env_type
   USE pw_poisson_types,                ONLY: do_ewald_ewald,&
                                              do_ewald_none,&
                                              do_ewald_pme,&
                                              do_ewald_spme
   USE qmmm_gaussian_types,             ONLY: qmmm_gaussian_p_type,&
                                              qmmm_gaussian_type
   USE qmmm_types_low,                  ONLY: qmmm_per_pot_p_type,&
                                              qmmm_per_pot_type,&
                                              qmmm_pot_p_type
#include "./base/base_uses.f90"
 
   IMPLICIT NONE
   PRIVATE
 
   LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .true.
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qmmm_per_elpot'
   PUBLIC :: qmmm_per_potential_init, qmmm_ewald_potential_init
 
CONTAINS
 
! **************************************************************************************************
!> \brief Initialize the QMMM potential stored on vector,
!>      according the qmmm_coupl_type
!> \param qmmm_coupl_type ...
!> \param per_potentials ...
!> \param potentials ...
!> \param pgfs ...
!> \param qm_cell_small ...
!> \param mm_cell ...
!> \param compatibility ...
!> \param qmmm_periodic ...
!> \param print_section ...
!> \param eps_mm_rspace ...
!> \param maxchrg ...
!> \param ncp ...
!> \param ncpl ...
!> \par History
!>      09.2004 created [tlaino]
!> \author Teodoro Laino
! **************************************************************************************************
   SUBROUTINE qmmm_per_potential_init(qmmm_coupl_type, per_potentials, potentials, &
                                      pgfs, qm_cell_small, mm_cell, compatibility, qmmm_periodic, print_section, &
                                      eps_mm_rspace, maxchrg, ncp, ncpl)
      INTEGER, INTENT(IN)                                :: qmmm_coupl_type
      TYPE(qmmm_per_pot_p_type), DIMENSION(:), POINTER   :: per_potentials
      TYPE(qmmm_pot_p_type), DIMENSION(:), POINTER       :: potentials
      TYPE(qmmm_gaussian_p_type), DIMENSION(:), POINTER  :: pgfs
      TYPE(cell_type), POINTER                           :: qm_cell_small, mm_cell
      LOGICAL, INTENT(IN)                                :: compatibility
      TYPE(section_vals_type), POINTER                   :: qmmm_periodic, print_section
      REAL(kind=dp), INTENT(IN)                          :: eps_mm_rspace, maxchrg
      INTEGER, INTENT(IN)                                :: ncp(3), ncpl(3)
 
      INTEGER                                            :: i, idim, ig, ig_start, iw, ix, iy, iz, &
                                                            k, kmax(3), n_rep_real(3), &
                                                            n_rep_real_val, ncoarsel, ncoarset, &
                                                            ndim, output_unit
      INTEGER, DIMENSION(:), POINTER                     :: mm_atom_index
      REAL(kind=dp)                                      :: ak, alpha, box(3), fac(3), fs, g, g2, &
                                                            gk, gmax, mymaxradius, npl, npt, &
                                                            prefactor, rc, rc2, rmax, tmp, vec(3), &
                                                            vol
      REAL(kind=dp), DIMENSION(:), POINTER               :: gx, gy, gz, lg
      TYPE(cp_logger_type), POINTER                      :: logger
      TYPE(qmmm_gaussian_type), POINTER                  :: pgf
 
      NULLIFY (lg, gx, gy, gz)
      ncoarset = product(ncp)
      ncoarsel = product(ncpl)
      logger => cp_get_default_logger()
      output_unit = cp_logger_get_default_io_unit(logger)
      rmax = sqrt(mm_cell%hmat(1, 1)**2 + &
                  mm_cell%hmat(2, 2)**2 + &
                  mm_cell%hmat(3, 3)**2)
      CALL section_vals_val_get(qmmm_periodic, "GMAX", r_val=gmax)
      CALL section_vals_val_get(qmmm_periodic, "REPLICA", i_val=n_rep_real_val)
      fac = 2.0e0_dp*pi/[mm_cell%hmat(1, 1), mm_cell%hmat(2, 2), mm_cell%hmat(3, 3)]
      kmax = ceiling(gmax/fac)
      vol = mm_cell%hmat(1, 1)* &
            mm_cell%hmat(2, 2)* &
            mm_cell%hmat(3, 3)
      ndim = (kmax(1) + 1)*(2*kmax(2) + 1)*(2*kmax(3) + 1)
      ig_start = 1
      n_rep_real = n_rep_real_val
      IF (compatibility .AND. (qmmm_coupl_type == do_qmmm_gauss)) ig_start = 2
 
      cpassert(.NOT. ASSOCIATED(per_potentials))
      ALLOCATE (per_potentials(SIZE(pgfs)))
      cpassert(SIZE(pgfs) == SIZE(potentials))
      potential_type: DO k = 1, SIZE(pgfs)
 
         rc = pgfs(k)%pgf%Elp_Radius
         ALLOCATE (per_potentials(k)%Pot)
         SELECT CASE (qmmm_coupl_type)
         CASE (do_qmmm_coulomb, do_qmmm_pcharge)
            CALL cp_abort(__location__, &
                          "QM/MM coupling type Coulomb or Point_Charge not "// &
                          "yet implemented for qmmm_per_potential_init")
         CASE (do_qmmm_gauss, do_qmmm_swave)
            ALLOCATE (lg(ndim))
            ALLOCATE (gx(ndim))
            ALLOCATE (gy(ndim))
            ALLOCATE (gz(ndim))
         END SELECT
 
         lg = 0.0_dp
         gx = 0.0_dp
         gy = 0.0_dp
         gz = 0.0_dp
 
         SELECT CASE (qmmm_coupl_type)
         CASE (do_qmmm_coulomb, do_qmmm_pcharge)
            CALL cp_abort(__location__, &
                          "QM/MM coupling type Coulomb or Point_Charge not "// &
                          "yet implemented for qmmm_per_potential_init")
         CASE (do_qmmm_gauss, do_qmmm_swave)
            pgf => pgfs(k)%pgf
            idim = 0
            DO ix = 0, kmax(1)
               DO iy = -kmax(2), kmax(2)
                  DO iz = -kmax(3), kmax(3)
                     idim = idim + 1
                     IF (ix == 0 .AND. iy == 0 .AND. iz == 0) THEN
                        DO ig = ig_start, pgf%number_of_gaussians
                           gk = pgf%Gk(ig)
                           ak = pgf%Ak(ig)*pi**(3.0_dp/2.0_dp)*gk**3.0_dp
                           lg(idim) = lg(idim) - ak
                        END DO
                     ELSE
                        fs = 2.0_dp; IF (ix == 0) fs = 1.0_dp
                        vec = fac*[real(ix, kind=dp), real(iy, kind=dp), real(iz, kind=dp)]
                        g2 = dot_product(vec, vec)
                        rc2 = rc*rc
                        g = sqrt(g2)
                        IF (qmmm_coupl_type == do_qmmm_gauss) THEN
                           lg(idim) = 4.0_dp*pi/g2*exp(-(g2*rc2)/4.0_dp)
                        ELSEIF (qmmm_coupl_type == do_qmmm_swave) THEN
                           tmp = 4.0_dp/rc2
                           lg(idim) = 4.0_dp*pi*tmp**2/(g2*(g2 + tmp)**2)
                        END IF
                        DO ig = ig_start, pgf%number_of_gaussians
                           gk = pgf%Gk(ig)
                           ak = pgf%Ak(ig)*pi**(3.0_dp/2.0_dp)*gk**3.0_dp
                           lg(idim) = lg(idim) - ak*exp(-(g*gk)**2.0_dp/4.0_dp)
                        END DO
                     END IF
                     lg(idim) = fs*lg(idim)*1.0_dp/vol
                     gx(idim) = fac(1)*real(ix, kind=dp)
                     gy(idim) = fac(2)*real(iy, kind=dp)
                     gz(idim) = fac(3)*real(iz, kind=dp)
                  END DO
               END DO
            END DO
 
            IF (all(n_rep_real == -1)) THEN
               mymaxradius = 0.0_dp
               DO i = 1, pgf%number_of_gaussians
                  IF (pgf%Gk(i) /= 0.0_dp) THEN
                     alpha = 1.0_dp/pgf%Gk(i)
                     alpha = alpha*alpha
                     prefactor = pgf%Ak(i)*maxchrg
                     mymaxradius = max(mymaxradius, exp_radius(0, alpha, eps_mm_rspace, prefactor, rlow=mymaxradius))
                  END IF
               END DO
               box(1) = (qm_cell_small%hmat(1, 1) - mm_cell%hmat(1, 1))/2.0_dp
               box(2) = (qm_cell_small%hmat(2, 2) - mm_cell%hmat(2, 2))/2.0_dp
               box(3) = (qm_cell_small%hmat(3, 3) - mm_cell%hmat(3, 3))/2.0_dp
               IF (any(box > 0.0_dp)) THEN
                  cpabort("QM cell and MM cell differ in size")
               END IF
               n_rep_real(1) = ceiling((box(1) + mymaxradius)/mm_cell%hmat(1, 1))
               n_rep_real(2) = ceiling((box(2) + mymaxradius)/mm_cell%hmat(2, 2))
               n_rep_real(3) = ceiling((box(3) + mymaxradius)/mm_cell%hmat(3, 3))
            END IF
 
         CASE DEFAULT
            DEALLOCATE (per_potentials(k)%Pot)
            IF (output_unit > 0) WRITE (output_unit, '(A)') " QMMM Periodic Potential - not Initialized!"
            cycle potential_type
         END SELECT
 
         NULLIFY (mm_atom_index)
         ALLOCATE (mm_atom_index(SIZE(potentials(k)%pot%mm_atom_index)))
         mm_atom_index = potentials(k)%pot%mm_atom_index
 
         NULLIFY (per_potentials(k)%Pot%LG, per_potentials(k)%Pot%mm_atom_index, &
                  per_potentials(k)%Pot%gx, per_potentials(k)%Pot%gy, per_potentials(k)%Pot%gz)
         CALL qmmm_per_pot_type_create(per_potentials(k)%Pot, lg=lg, gx=gx, gy=gy, gz=gz, &
                                       gmax=gmax, kmax=kmax, n_rep_real=n_rep_real, &
                                       fac=fac, mm_atom_index=mm_atom_index, &
                                       mm_cell=mm_cell, &
                                       qmmm_per_section=qmmm_periodic, print_section=print_section)
 
         iw = cp_print_key_unit_nr(logger, print_section, "PERIODIC_INFO", &
                                   extension=".log")
         IF (iw > 0) THEN
            npt = real(ncoarset, kind=dp)*real(ndim, kind=dp)*real(SIZE(mm_atom_index), kind=dp)
            npl = real(ncoarsel, kind=dp)*real(ndim, kind=dp)*real(SIZE(mm_atom_index), kind=dp)
            WRITE (unit=iw, fmt="(/,T2,A)") repeat("-", 79)
            WRITE (unit=iw, fmt="(T2,A,T20,A,T80,A)") "-", "QMMM PERIODIC BOUNDARY CONDITION INFO", "-"
            WRITE (unit=iw, fmt="(T2,A)") repeat("-", 79)
            WRITE (unit=iw, fmt="(T2,A,T10,A,F15.6,T50,A,3I5,T80,A)") "-", "RADIUS  =", rc, "REPLICA =", n_rep_real, "-"
            WRITE (unit=iw, fmt="(T2,A,T10,A,F15.6,T50,A,I15,T80,A)") "-", "MINGVAL =", minval(abs(lg)), &
               "GPOINTS =", ndim, "-"
            WRITE (unit=iw, fmt="(T2,A,T10,A,3I5,T50,A,3I5,T80,A)") "-", "NCOARSL =", ncpl, &
               "NCOARST =", ncp, "-"
            WRITE (unit=iw, fmt="(T2,A,T10,A,F15.0,T50,A,F15.0,T80,A)") "-", "NFLOP-L ~", npl, &
               "NFLOP-T ~", npt, "-"
            WRITE (unit=iw, fmt="(T2,A)") repeat("-", 79)
         END IF
         CALL cp_print_key_finished_output(iw, logger, print_section, &
                                           "PERIODIC_INFO")
 
      END DO potential_type
 
   END SUBROUTINE qmmm_per_potential_init
 
! **************************************************************************************************
!> \brief Creates the qmmm_pot_type structure
!> \param Pot ...
!> \param LG ...
!> \param gx ...
!> \param gy ...
!> \param gz ...
!> \param GMax ...
!> \param Kmax ...
!> \param n_rep_real ...
!> \param Fac ...
!> \param mm_atom_index ...
!> \param mm_cell ...
!> \param qmmm_per_section ...
!> \param print_section ...
!> \par History
!>      09.2004 created [tlaino]
!> \author Teodoro Laino
! **************************************************************************************************
   SUBROUTINE qmmm_per_pot_type_create(Pot, LG, gx, gy, gz, GMax, Kmax, n_rep_real, &
                                       Fac, mm_atom_index, mm_cell, qmmm_per_section, print_section)
      TYPE(qmmm_per_pot_type), POINTER                   :: pot
      REAL(kind=dp), DIMENSION(:), POINTER               :: lg, gx, gy, gz
      REAL(kind=dp), INTENT(IN)                          :: gmax
      INTEGER, INTENT(IN)                                :: kmax(3), n_rep_real(3)
      REAL(kind=dp), INTENT(IN)                          :: fac(3)
      INTEGER, DIMENSION(:), POINTER                     :: mm_atom_index
      TYPE(cell_type), POINTER                           :: mm_cell
      TYPE(section_vals_type), POINTER                   :: qmmm_per_section, print_section
 
      INTEGER                                            :: npts(3)
      INTEGER, DIMENSION(:), POINTER                     :: ngrids
      REAL(kind=dp)                                      :: hmat(3, 3)
      TYPE(section_vals_type), POINTER                   :: grid_print_section
 
      pot%LG => lg
      pot%gx => gx
      pot%gy => gy
      pot%gz => gz
      pot%mm_atom_index => mm_atom_index
      pot%Gmax = gmax
      pot%Kmax = kmax
      pot%n_rep_real = n_rep_real
      pot%Fac = fac
      !
      ! Setting Up Fit Procedure
      !
      NULLIFY (pot%pw_grid)
      NULLIFY (pot%pw_pool)
      NULLIFY (pot%TabLR, ngrids)
      CALL section_vals_val_get(qmmm_per_section, "ngrids", i_vals=ngrids)
      npts = ngrids
      hmat = mm_cell%hmat
 
      grid_print_section => section_vals_get_subs_vals(print_section, "GRID_INFORMATION")
      CALL setup_ewald_spline(pw_grid=pot%pw_grid, pw_pool=pot%pw_pool, coeff=pot%TabLR, &
                              lg=lg, gx=gx, gy=gy, gz=gz, hmat=hmat, npts=npts, param_section=qmmm_per_section, &
                              tag="qmmm", print_section=grid_print_section)
 
   END SUBROUTINE qmmm_per_pot_type_create
 
! **************************************************************************************************
!> \brief Initialize the QMMM Ewald potential needed for QM-QM Coupling using
!>      point charges
!> \param ewald_env ...
!> \param ewald_pw ...
!> \param qmmm_coupl_type ...
!> \param mm_cell ...
!> \param para_env ...
!> \param qmmm_periodic ...
!> \param print_section ...
!> \par History
!>      10.2014 created [JGH]
!> \author JGH
! **************************************************************************************************
   SUBROUTINE qmmm_ewald_potential_init(ewald_env, ewald_pw, qmmm_coupl_type, mm_cell, para_env, &
                                        qmmm_periodic, print_section)
      TYPE(ewald_environment_type), POINTER              :: ewald_env
      TYPE(ewald_pw_type), POINTER                       :: ewald_pw
      INTEGER, INTENT(IN)                                :: qmmm_coupl_type
      TYPE(cell_type), POINTER                           :: mm_cell
      TYPE(mp_para_env_type), POINTER                    :: para_env
      TYPE(section_vals_type), POINTER                   :: qmmm_periodic, print_section
 
      INTEGER                                            :: ewald_type, gmax(3), iw, o_spline, ounit
      LOGICAL                                            :: do_multipoles
      REAL(kind=dp)                                      :: alpha, rcut
      TYPE(cp_logger_type), POINTER                      :: logger
      TYPE(section_vals_type), POINTER                   :: ewald_print_section, ewald_section, &
                                                            poisson_section
 
      logger => cp_get_default_logger()
      ounit = cp_logger_get_default_io_unit(logger)
      cpassert(.NOT. ASSOCIATED(ewald_env))
      cpassert(.NOT. ASSOCIATED(ewald_pw))
 
      ! Create Ewald environments
      poisson_section => section_vals_get_subs_vals(qmmm_periodic, "POISSON")
      ALLOCATE (ewald_env)
      CALL ewald_env_create(ewald_env, para_env)
      CALL ewald_env_set(ewald_env, poisson_section=poisson_section)
      ewald_section => section_vals_get_subs_vals(poisson_section, "EWALD")
      CALL read_ewald_section(ewald_env, ewald_section)
      ewald_print_section => section_vals_get_subs_vals(print_section, "GRID_INFORMATION")
      ALLOCATE (ewald_pw)
      CALL ewald_pw_create(ewald_pw, ewald_env, mm_cell, mm_cell, print_section=ewald_print_section)
 
      CALL ewald_env_get(ewald_env, ewald_type=ewald_type, do_multipoles=do_multipoles, &
                         gmax=gmax, o_spline=o_spline, alpha=alpha, rcut=rcut)
      IF (do_multipoles) &
         cpabort("No multipole force fields allowed in QM-QM Ewald long range correction")
 
      SELECT CASE (qmmm_coupl_type)
      CASE (do_qmmm_coulomb)
         cpabort("QM-QM long range correction not possible with COULOMB coupling")
      CASE (do_qmmm_pcharge)
         ! OK
      CASE (do_qmmm_gauss, do_qmmm_swave)
         cpabort("QM-QM long range correction not possible with GAUSS/SWAVE coupling")
      CASE DEFAULT
         ! We should never get to this point
         cpabort("Unknown QM/MM coupling type")
      END SELECT
 
      iw = cp_print_key_unit_nr(logger, print_section, "PERIODIC_INFO", extension=".log")
      IF (iw > 0) THEN
         WRITE (unit=iw, fmt="(/,T2,A)") repeat("-", 79)
         WRITE (unit=iw, fmt="(T2,A,T20,A,T80,A)") "-", "QMMM PERIODIC BOUNDARY CONDITION INFO", "-"
         WRITE (unit=iw, fmt="(T2,A,T25,A,T80,A)") "-", "QM-QM Long Range Correction", "-"
         WRITE (unit=iw, fmt="(T2,A)") repeat("-", 79)
         SELECT CASE (ewald_type)
         CASE (do_ewald_none)
            cpabort("QM-QM long range correction not compatible with Ewald=NONE")
         CASE (do_ewald_pme)
            cpabort("QM-QM long range correction not possible with Ewald=PME")
         CASE (do_ewald_ewald)
            cpabort("QM-QM long range correction not possible with Ewald method")
         CASE (do_ewald_spme)
            WRITE (unit=iw, fmt="(T2,A,T35,A,T75,A,T80,A)") "-", "Ewald type", "SPME", "-"
            WRITE (unit=iw, fmt="(T2,A,T35,A,T61,3I6,T80,A)") "-", "GMAX values", gmax, "-"
            WRITE (unit=iw, fmt="(T2,A,T35,A,T73,I6,T80,A)") "-", "Spline order", o_spline, "-"
            WRITE (unit=iw, fmt="(T2,A,T35,A,T67,F12.4,T80,A)") "-", "Alpha value", alpha, "-"
            WRITE (unit=iw, fmt="(T2,A,T35,A,T67,F12.4,T80,A)") "-", "Real space cutoff value", rcut, "-"
         END SELECT
         WRITE (unit=iw, fmt="(T2,A)") repeat("-", 79)
      END IF
      CALL cp_print_key_finished_output(iw, logger, print_section, "PERIODIC_INFO")
 
   END SUBROUTINE qmmm_ewald_potential_init
 
END MODULE qmmm_per_elpot