qmmm_elpot.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                                                      !
!--------------------------------------------------------------------------------------------------!
 
! **************************************************************************************************
!> \par History
!>      09.2004 created [tlaino]
!> \author Teodoro Laino
! **************************************************************************************************
MODULE qmmm_elpot
   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_p_file,&
                                              cp_print_key_finished_output,&
                                              cp_print_key_should_output,&
                                              cp_print_key_unit_nr
   USE input_constants,                 ONLY: do_qmmm_coulomb,&
                                              do_qmmm_gauss,&
                                              do_qmmm_pcharge,&
                                              do_qmmm_swave
   USE input_section_types,             ONLY: section_vals_type
   USE kinds,                           ONLY: default_path_length,&
                                              default_string_length,&
                                              dp
   USE mathconstants,                   ONLY: rootpi
   USE memory_utilities,                ONLY: reallocate
   USE qmmm_gaussian_types,             ONLY: qmmm_gaussian_p_type,&
                                              qmmm_gaussian_type
   USE qmmm_types_low,                  ONLY: qmmm_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_elpot'
   PUBLIC :: qmmm_potential_init
 
CONTAINS
 
! **************************************************************************************************
!> \brief Initialize the QMMM potential stored on vector,
!>      according the qmmm_coupl_type
!> \param qmmm_coupl_type ...
!> \param mm_el_pot_radius ...
!> \param potentials ...
!> \param pgfs ...
!> \param mm_cell ...
!> \param compatibility ...
!> \param print_section ...
!> \par History
!>      09.2004 created [tlaino]
!> \author Teodoro Laino
! **************************************************************************************************
   SUBROUTINE qmmm_potential_init(qmmm_coupl_type, mm_el_pot_radius, potentials, &
                                  pgfs, mm_cell, compatibility, print_section)
      INTEGER, INTENT(IN)                                :: qmmm_coupl_type
      REAL(kind=dp), DIMENSION(:), POINTER               :: mm_el_pot_radius
      TYPE(qmmm_pot_p_type), DIMENSION(:), POINTER       :: potentials
      TYPE(qmmm_gaussian_p_type), DIMENSION(:), POINTER  :: pgfs
      TYPE(cell_type), POINTER                           :: mm_cell
      LOGICAL, INTENT(IN)                                :: compatibility
      TYPE(section_vals_type), POINTER                   :: print_section
 
      REAL(kind=dp), PARAMETER                           :: dx = 0.05_dp
 
      CHARACTER(LEN=default_path_length)                 :: myformat
      CHARACTER(LEN=default_string_length)               :: rc_s
      INTEGER                                            :: i, ig, ig_start, j, k, myind, ndim, np, &
                                                            output_unit, unit_nr
      INTEGER, DIMENSION(:), POINTER                     :: mm_atom_index
      LOGICAL                                            :: found
      REAL(kind=dp)                                      :: a, g, rc, rmax, rmin, t, t1, t2, x
      REAL(kind=dp), DIMENSION(:), POINTER               :: radius
      REAL(kind=dp), DIMENSION(:, :), POINTER            :: pot0_2
      TYPE(cp_logger_type), POINTER                      :: logger
      TYPE(qmmm_gaussian_type), POINTER                  :: pgf
 
      logger => cp_get_default_logger()
      output_unit = cp_logger_get_default_io_unit(logger)
      rmin = 0.0_dp
      rmax = sqrt(mm_cell%hmat(1, 1)**2 + &
                  mm_cell%hmat(2, 2)**2 + &
                  mm_cell%hmat(3, 3)**2)
      np = ceiling(rmax/dx) + 1
      !
      ! Preprocessing
      !
      IF (SIZE(mm_el_pot_radius) /= 0) THEN
         ALLOCATE (radius(1))
         radius(1) = mm_el_pot_radius(1)
      ELSE
         ALLOCATE (radius(0))
      END IF
      loop_on_all_values: DO i = 2, SIZE(mm_el_pot_radius)
         found = .false.
         loop_on_found_values: DO j = 1, SIZE(radius)
            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)
            ndim = ndim + 1
            CALL reallocate(radius, 1, ndim)
            radius(ndim) = mm_el_pot_radius(i)
         END IF
      END DO loop_on_all_values
      !
      cpassert(.NOT. ASSOCIATED(potentials))
      ALLOCATE (potentials(SIZE(radius)))
 
      potential_type: DO k = 1, SIZE(radius)
 
         rc = radius(k)
         ALLOCATE (potentials(k)%Pot)
         SELECT CASE (qmmm_coupl_type)
         CASE (do_qmmm_coulomb)
            NULLIFY (pot0_2)
         CASE (do_qmmm_pcharge)
            NULLIFY (pot0_2)
         CASE (do_qmmm_gauss, do_qmmm_swave)
            ALLOCATE (pot0_2(2, np))
         END SELECT
 
         SELECT CASE (qmmm_coupl_type)
         CASE (do_qmmm_coulomb, do_qmmm_pcharge)
            ! Do Nothing
         CASE (do_qmmm_gauss, do_qmmm_swave)
            IF (qmmm_coupl_type == do_qmmm_gauss) THEN
               ! Smooth Coulomb Potential ::  Erf(x/rc)/x
               pot0_2(1, 1) = 2.0_dp/(rootpi*rc)
               pot0_2(2, 1) = 0.0_dp
               x = 0.0_dp
               DO i = 2, np
                  x = x + dx
                  pot0_2(1, i) = erf(x/rc)/x
                  t = 2._dp/(rootpi*x*rc)*exp(-(x/rc)**2)
                  pot0_2(2, i) = (t - pot0_2(1, i)/x)*dx
               END DO
            ELSEIF (qmmm_coupl_type == do_qmmm_swave) THEN
               ! S-wave expansion :: 1/x - exp(-2*x/rc) * ( 1/x - 1/rc )
               pot0_2(1, 1) = 1.0_dp/rc
               pot0_2(2, 1) = 0.0_dp
               x = 0.0_dp
               DO i = 2, np
                  x = x + dx
                  t = exp(-2.0_dp*x/rc)/rc
                  pot0_2(1, i) = (1.0_dp - t*(rc + x))/x
                  pot0_2(2, i) = ((t*(rc**2 + 2.0_dp*rc*x + 2.0_dp*x**2)/rc - 1.0_dp)/x**2)*dx
               END DO
            END IF
            pgf => pgfs(k)%pgf
            cpassert(pgf%Elp_Radius == rc)
            ig_start = 1
            IF (compatibility .AND. (qmmm_coupl_type == do_qmmm_gauss)) ig_start = 2
            DO ig = ig_start, pgf%number_of_gaussians
               a = pgf%Ak(ig)
               g = pgf%Gk(ig)
               pot0_2(1, 1) = pot0_2(1, 1) - a
               x = 0.0_dp
               DO i = 2, np
                  x = x + dx
                  t = exp(-(x/g)**2)*a
                  t1 = 1/g**2
                  t2 = t1*t
                  pot0_2(1, i) = pot0_2(1, i) - t
                  pot0_2(2, i) = pot0_2(2, i) + 2.0_dp*x*t2*dx
               END DO
            END DO
 
            ! Print info on the unidimensional MM electrostatic potential
            IF (btest(cp_print_key_should_output(logger%iter_info, print_section, "MM_POTENTIAL") &
                      , cp_p_file)) THEN
               WRITE (rc_s, '(F6.3)') rc
               unit_nr = cp_print_key_unit_nr(logger, print_section, "MM_POTENTIAL", &
                                              extension="_rc="//trim(adjustl(rc_s))//".data")
               IF (unit_nr > 0) THEN
                  WRITE (unit_nr, '(A)') "# MM ELECTROSTATIC POTENTIAL - UNIDIMENSIONAL - ATOMIC UNITS"
                  WRITE (unit_nr, '(A,I5)') "# MM ELECTROSTATIC POTENTIAL - Nr. of Gaussians:", pgf%number_of_gaussians
                  WRITE (unit_nr, '(A,T10,A,T30,A,T300,A)') "#", "Xval", "Gaussians", "LongRange"
                  myformat = "T10,F15.9,T30,"
                  DO ig = 1, pgf%number_of_gaussians
                     myind = index(myformat, " ")
                     WRITE (myformat(myind:), '(A6)') "F12.9,"
                  END DO
                  myind = index(myformat, " ") - 1
                  myformat = myformat(1:myind)//"T300,F15.9"
                  myind = index(myformat, " ") - 1
                  x = 0.0_dp
                  DO i = 1, np
                     WRITE (unit_nr, '('//myformat(1:myind)//')') &
                        x, (exp(-(x/pgf%Gk(ig))**2)*pgf%Ak(ig), ig=1, pgf%number_of_gaussians), pot0_2(1, i)
                     x = x + dx
                  END DO
               END IF
               CALL cp_print_key_finished_output(unit_nr, logger, print_section, &
                                                 "MM_POTENTIAL")
            END IF
         CASE DEFAULT
            DEALLOCATE (potentials(k)%Pot)
            IF (output_unit > 0) WRITE (output_unit, '(A)') " QMMM Potential - Spline Interpolation - not Initialized!"
            cycle potential_type
         END SELECT
         NULLIFY (mm_atom_index)
         ALLOCATE (mm_atom_index(1))
         ! Build mm_atom_index List
         DO j = 1, SIZE(mm_el_pot_radius)
            IF (rc .EQ. mm_el_pot_radius(j)) THEN
               ndim = SIZE(mm_atom_index)
               mm_atom_index(ndim) = j
               CALL reallocate(mm_atom_index, 1, ndim + 1)
            END IF
         END DO
         CALL reallocate(mm_atom_index, 1, ndim)
 
         NULLIFY (potentials(k)%Pot%pot0_2)
         CALL qmmm_pot_type_create(potentials(k)%Pot, pot0_2=pot0_2, &
                                   rmax=rmax, rmin=rmin, dx=dx, rc=rc, npts=np, &
                                   mm_atom_index=mm_atom_index)
 
      END DO potential_type
      DEALLOCATE (radius)
   END SUBROUTINE qmmm_potential_init
 
! **************************************************************************************************
!> \brief Creates the qmmm_pot_type structure
!> \param Pot ...
!> \param pot0_2 ...
!> \param Rmax ...
!> \param Rmin ...
!> \param dx ...
!> \param npts ...
!> \param rc ...
!> \param mm_atom_index ...
!> \par History
!>      09.2004 created [tlaino]
!> \author Teodoro Laino
! **************************************************************************************************
   SUBROUTINE qmmm_pot_type_create(Pot, pot0_2, Rmax, Rmin, dx, npts, rc, &
                                   mm_atom_index)
      TYPE(qmmm_pot_type), POINTER                       :: pot
      REAL(kind=dp), DIMENSION(:, :), POINTER            :: pot0_2
      REAL(kind=dp), INTENT(IN)                          :: rmax, rmin, dx
      INTEGER, INTENT(IN)                                :: npts
      REAL(kind=dp), INTENT(IN)                          :: rc
      INTEGER, DIMENSION(:), POINTER                     :: mm_atom_index
 
      pot%pot0_2 => pot0_2
      pot%mm_atom_index => mm_atom_index
      pot%Rmax = rmax
      pot%Rmin = rmin
      pot%Rc = rc
      pot%dx = dx
      pot%npts = npts
 
   END SUBROUTINE qmmm_pot_type_create
 
END MODULE qmmm_elpot