qmmm_gaussian_init.F

Go to the documentation of this file.

!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!
 
! **************************************************************************************************
!> \brief Initialize the use of the gaussians to treat the QMMM
!>      coupling potential
!> \par History
!>      6.2004 created [tlaino]
!> \author Teodoro Laino
! **************************************************************************************************
MODULE qmmm_gaussian_init
   USE ao_util,                         ONLY: exp_radius
   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 gaussian_gridlevels,             ONLY: gaussian_gridlevel,&
                                              gridlevel_info_type
   USE input_constants,                 ONLY: do_qmmm_gauss,&
                                              do_qmmm_swave
   USE input_section_types,             ONLY: section_vals_type,&
                                              section_vals_val_get
   USE kinds,                           ONLY: dp
   USE memory_utilities,                ONLY: reallocate
   USE message_passing,                 ONLY: mp_para_env_type
   USE pw_env_types,                    ONLY: pw_env_get,&
                                              pw_env_type
   USE pw_pool_types,                   ONLY: pw_pool_p_type
   USE qmmm_gaussian_data,              ONLY: max_geep_lib_gauss,&
                                              min_geep_lib_gauss
   USE qmmm_gaussian_input,             ONLY: read_mm_potential,&
                                              set_mm_potential_erf,&
                                              set_mm_potential_swave
   USE qmmm_gaussian_types,             ONLY: qmmm_gaussian_p_type,&
                                              qmmm_gaussian_type
#include "./base/base_uses.f90"
 
   IMPLICIT NONE
   PRIVATE
 
   LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .true.
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qmmm_gaussian_init'
   PUBLIC :: qmmm_gaussian_initialize
 
CONTAINS
 
! **************************************************************************************************
!> \brief Initialize the Gaussian QMMM Environment
!> \param qmmm_gaussian_fns ...
!> \param para_env ...
!> \param pw_env ...
!> \param mm_el_pot_radius ...
!> \param mm_el_pot_radius_corr ...
!> \param qmmm_coupl_type ...
!> \param eps_mm_rspace ...
!> \param maxradius ...
!> \param maxchrg ...
!> \param compatibility ...
!> \param print_section ...
!> \param qmmm_section ...
!> \par History
!>      06.2004 created [tlaino]
!> \author Teodoro Laino
! **************************************************************************************************
   SUBROUTINE qmmm_gaussian_initialize(qmmm_gaussian_fns, para_env, pw_env, &
                                       mm_el_pot_radius, mm_el_pot_radius_corr, &
                                       qmmm_coupl_type, eps_mm_rspace, maxradius, maxchrg, compatibility, &
                                       print_section, qmmm_section)
      TYPE(qmmm_gaussian_p_type), DIMENSION(:), POINTER  :: qmmm_gaussian_fns
      TYPE(mp_para_env_type), POINTER                    :: para_env
      TYPE(pw_env_type), POINTER                         :: pw_env
      REAL(kind=dp), DIMENSION(:), POINTER               :: mm_el_pot_radius, mm_el_pot_radius_corr
      INTEGER, INTENT(IN)                                :: qmmm_coupl_type
      REAL(kind=dp), INTENT(IN)                          :: eps_mm_rspace
      REAL(kind=dp), DIMENSION(:), POINTER               :: maxradius
      REAL(kind=dp), INTENT(IN)                          :: maxchrg
      LOGICAL, INTENT(IN)                                :: compatibility
      TYPE(section_vals_type), POINTER                   :: print_section, qmmm_section
 
      INTEGER                                            :: i, ilevel, j, ndim, num_geep_gauss, &
                                                            output_unit
      LOGICAL                                            :: found, use_geep_lib
      REAL(kind=dp)                                      :: alpha, mymaxradius, prefactor
      REAL(kind=dp), DIMENSION(:), POINTER               :: c_radius, radius
      TYPE(cp_logger_type), POINTER                      :: logger
      TYPE(gridlevel_info_type), POINTER                 :: gridlevel_info
      TYPE(pw_pool_p_type), DIMENSION(:), POINTER        :: pools
      TYPE(qmmm_gaussian_type), POINTER                  :: mypgf
 
! Statements
 
      NULLIFY (mypgf, gridlevel_info, radius, c_radius, logger)
      logger => cp_get_default_logger()
      CALL section_vals_val_get(qmmm_section, "USE_GEEP_LIB", i_val=num_geep_gauss)
      IF (num_geep_gauss == 0) THEN
         use_geep_lib = .false.
      ELSE
         use_geep_lib = .true.
         cpassert(num_geep_gauss >= min_geep_lib_gauss)
         cpassert(num_geep_gauss <= max_geep_lib_gauss)
      END IF
      SELECT CASE (qmmm_coupl_type)
      CASE (do_qmmm_gauss, do_qmmm_swave)
         !
         ! Preprocessing...
         !
         ALLOCATE (radius(1))
         ALLOCATE (c_radius(1))
         ndim = SIZE(radius)
         loop_on_all_values: DO i = 1, SIZE(mm_el_pot_radius)
            found = .false.
            loop_on_found_values: DO j = 1, SIZE(radius) - 1
               IF (mm_el_pot_radius(i) .EQ. radius(j)) THEN
                  found = .true.
                  EXIT loop_on_found_values
               END IF
            END DO loop_on_found_values
            IF (.NOT. found) THEN
               ndim = SIZE(radius)
               radius(ndim) = mm_el_pot_radius(i)
               c_radius(ndim) = mm_el_pot_radius_corr(i)
               ndim = ndim + 1
               CALL reallocate(radius, 1, ndim)
               CALL reallocate(c_radius, 1, ndim)
            END IF
         END DO loop_on_all_values
         !
         IF (ndim - 1 > 0) THEN
            CALL reallocate(radius, 1, ndim - 1)
            CALL reallocate(c_radius, 1, ndim - 1)
         ELSE IF (ndim - 1 == 0) THEN
            DEALLOCATE (radius)
            DEALLOCATE (c_radius)
         ELSE
            cpabort("")
         END IF
         !
         ALLOCATE (qmmm_gaussian_fns(ndim - 1))
         DO i = 1, ndim - 1
            NULLIFY (qmmm_gaussian_fns(i)%pgf)
            ALLOCATE (qmmm_gaussian_fns(i)%pgf)
            NULLIFY (qmmm_gaussian_fns(i)%pgf%Ak)
            NULLIFY (qmmm_gaussian_fns(i)%pgf%Gk)
            NULLIFY (qmmm_gaussian_fns(i)%pgf%grid_level)
            !
            ! Default Values
            !
            qmmm_gaussian_fns(i)%pgf%Elp_Radius = radius(i)
            qmmm_gaussian_fns(i)%pgf%Elp_Radius_corr = c_radius(i)
         END DO
         IF (ASSOCIATED(radius)) THEN
            DEALLOCATE (radius)
         END IF
         IF (ASSOCIATED(c_radius)) THEN
            DEALLOCATE (c_radius)
         END IF
         !
         IF (use_geep_lib) THEN
            IF (qmmm_coupl_type == do_qmmm_gauss) THEN
               CALL set_mm_potential_erf(qmmm_gaussian_fns, &
                                         compatibility, num_geep_gauss)
            ELSEIF (qmmm_coupl_type == do_qmmm_swave) THEN
               CALL set_mm_potential_swave(qmmm_gaussian_fns, &
                                           num_geep_gauss)
            END IF
         ELSE
            CALL read_mm_potential(para_env, qmmm_gaussian_fns, &
                                   (compatibility .AND. (qmmm_coupl_type == do_qmmm_gauss)), qmmm_section)
         END IF
         !
         CALL pw_env_get(pw_env, pw_pools=pools, gridlevel_info=gridlevel_info)
         ALLOCATE (maxradius(SIZE(pools)))
         maxradius = 0.0_dp
         DO j = 1, SIZE(qmmm_gaussian_fns)
            mypgf => qmmm_gaussian_fns(j)%pgf
            ALLOCATE (mypgf%grid_level(SIZE(mypgf%Ak)))
            mypgf%grid_level = 0
            mymaxradius = 0.0_dp
            DO i = 1, mypgf%number_of_gaussians
               IF (mypgf%Gk(i) /= 0.0_dp) THEN
                  alpha = 1.0_dp/mypgf%Gk(i)
                  alpha = alpha*alpha
                  ilevel = gaussian_gridlevel(gridlevel_info, alpha)
                  prefactor = mypgf%Ak(i)*maxchrg
                  mymaxradius = exp_radius(0, alpha, eps_mm_rspace, prefactor, rlow=mymaxradius)
                  maxradius(ilevel) = max(maxradius(ilevel), mymaxradius)
                  mypgf%grid_level(i) = ilevel
               END IF
            END DO
         END DO
         !
         ! End of gaussian initialization...
      CASE DEFAULT
         output_unit = cp_print_key_unit_nr(logger, print_section, "PROGRAM_RUN_INFO", &
                                            extension=".qmmmLog")
         IF (output_unit > 0) WRITE (output_unit, '(A)') " QMMM Gaussian Data Not Initialized!"
         CALL cp_print_key_finished_output(output_unit, logger, print_section, "PROGRAM_RUN_INFO")
      END SELECT
   END SUBROUTINE qmmm_gaussian_initialize
 
END MODULE qmmm_gaussian_init