LCOV - code coverage report
Current view: top level - src - qmmm_elpot.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 100.0 % 108 108
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 2 2

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \par History
      10              : !>      09.2004 created [tlaino]
      11              : !> \author Teodoro Laino
      12              : ! **************************************************************************************************
      13              : MODULE qmmm_elpot
      14              :    USE cell_types,                      ONLY: cell_type
      15              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      16              :                                               cp_logger_get_default_io_unit,&
      17              :                                               cp_logger_type
      18              :    USE cp_output_handling,              ONLY: cp_p_file,&
      19              :                                               cp_print_key_finished_output,&
      20              :                                               cp_print_key_should_output,&
      21              :                                               cp_print_key_unit_nr
      22              :    USE input_constants,                 ONLY: do_qmmm_coulomb,&
      23              :                                               do_qmmm_gauss,&
      24              :                                               do_qmmm_pcharge,&
      25              :                                               do_qmmm_swave
      26              :    USE input_section_types,             ONLY: section_vals_type
      27              :    USE kinds,                           ONLY: default_path_length,&
      28              :                                               default_string_length,&
      29              :                                               dp
      30              :    USE mathconstants,                   ONLY: rootpi
      31              :    USE memory_utilities,                ONLY: reallocate
      32              :    USE qmmm_gaussian_types,             ONLY: qmmm_gaussian_p_type,&
      33              :                                               qmmm_gaussian_type
      34              :    USE qmmm_types_low,                  ONLY: qmmm_Pot_Type,&
      35              :                                               qmmm_pot_p_type
      36              : #include "./base/base_uses.f90"
      37              : 
      38              :    IMPLICIT NONE
      39              :    PRIVATE
      40              : 
      41              :    LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
      42              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qmmm_elpot'
      43              :    PUBLIC :: qmmm_potential_init
      44              : 
      45              : CONTAINS
      46              : 
      47              : ! **************************************************************************************************
      48              : !> \brief Initialize the QMMM potential stored on vector,
      49              : !>      according the qmmm_coupl_type
      50              : !> \param qmmm_coupl_type ...
      51              : !> \param mm_el_pot_radius ...
      52              : !> \param potentials ...
      53              : !> \param pgfs ...
      54              : !> \param mm_cell ...
      55              : !> \param compatibility ...
      56              : !> \param print_section ...
      57              : !> \par History
      58              : !>      09.2004 created [tlaino]
      59              : !> \author Teodoro Laino
      60              : ! **************************************************************************************************
      61          404 :    SUBROUTINE qmmm_potential_init(qmmm_coupl_type, mm_el_pot_radius, potentials, &
      62              :                                   pgfs, mm_cell, compatibility, print_section)
      63              :       INTEGER, INTENT(IN)                                :: qmmm_coupl_type
      64              :       REAL(KIND=dp), DIMENSION(:), POINTER               :: mm_el_pot_radius
      65              :       TYPE(qmmm_pot_p_type), DIMENSION(:), POINTER       :: potentials
      66              :       TYPE(qmmm_gaussian_p_type), DIMENSION(:), POINTER  :: pgfs
      67              :       TYPE(cell_type), POINTER                           :: mm_cell
      68              :       LOGICAL, INTENT(IN)                                :: compatibility
      69              :       TYPE(section_vals_type), POINTER                   :: print_section
      70              : 
      71              :       REAL(KIND=dp), PARAMETER                           :: dx = 0.05_dp
      72              : 
      73              :       CHARACTER(LEN=default_path_length)                 :: myFormat
      74              :       CHARACTER(LEN=default_string_length)               :: rc_s
      75              :       INTEGER                                            :: I, ig, ig_start, J, K, myind, ndim, Np, &
      76              :                                                             output_unit, unit_nr
      77          404 :       INTEGER, DIMENSION(:), POINTER                     :: mm_atom_index
      78              :       LOGICAL                                            :: found
      79              :       REAL(KIND=dp)                                      :: A, G, rc, Rmax, Rmin, t, t1, t2, x
      80          404 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: radius
      81          404 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: pot0_2
      82              :       TYPE(cp_logger_type), POINTER                      :: logger
      83              :       TYPE(qmmm_gaussian_type), POINTER                  :: pgf
      84              : 
      85          808 :       logger => cp_get_default_logger()
      86          404 :       output_unit = cp_logger_get_default_io_unit(logger)
      87          404 :       Rmin = 0.0_dp
      88              :       Rmax = SQRT(mm_cell%hmat(1, 1)**2 + &
      89              :                   mm_cell%hmat(2, 2)**2 + &
      90          404 :                   mm_cell%hmat(3, 3)**2)
      91          404 :       np = CEILING(rmax/dx) + 1
      92              :       !
      93              :       ! Preprocessing
      94              :       !
      95          404 :       IF (SIZE(mm_el_pot_radius) /= 0) THEN
      96          402 :          ALLOCATE (radius(1))
      97          402 :          radius(1) = mm_el_pot_radius(1)
      98              :       ELSE
      99            2 :          ALLOCATE (radius(0))
     100              :       END IF
     101       187866 :       Loop_on_all_values: DO I = 2, SIZE(mm_el_pot_radius)
     102       187462 :          Found = .FALSE.
     103       357796 :          Loop_on_found_values: DO J = 1, SIZE(radius)
     104       357796 :             IF (mm_el_pot_radius(i) .EQ. radius(j)) THEN
     105              :                Found = .TRUE.
     106              :                EXIT Loop_on_found_values
     107              :             END IF
     108              :          END DO Loop_on_found_values
     109       187866 :          IF (.NOT. Found) THEN
     110              :             Ndim = SIZE(radius)
     111          232 :             Ndim = Ndim + 1
     112          232 :             CALL REALLOCATE(radius, 1, Ndim)
     113          232 :             radius(Ndim) = mm_el_pot_radius(i)
     114              :          END IF
     115              :       END DO Loop_on_all_values
     116              :       !
     117          404 :       CPASSERT(.NOT. ASSOCIATED(potentials))
     118         1844 :       ALLOCATE (potentials(SIZE(radius)))
     119              : 
     120         1038 :       Potential_Type: DO K = 1, SIZE(radius)
     121              : 
     122          634 :          rc = radius(K)
     123          634 :          ALLOCATE (potentials(K)%Pot)
     124          732 :          SELECT CASE (qmmm_coupl_type)
     125              :          CASE (do_qmmm_coulomb)
     126           98 :             NULLIFY (pot0_2)
     127              :          CASE (do_qmmm_pcharge)
     128           56 :             NULLIFY (pot0_2)
     129              :          CASE (do_qmmm_gauss, do_qmmm_swave)
     130         1490 :             ALLOCATE (pot0_2(2, np))
     131              :          END SELECT
     132              : 
     133          428 :          SELECT CASE (qmmm_coupl_type)
     134              :          CASE (do_qmmm_coulomb, do_qmmm_pcharge)
     135              :             ! Do Nothing
     136              :          CASE (do_qmmm_gauss, do_qmmm_swave)
     137          428 :             IF (qmmm_coupl_type == do_qmmm_gauss) THEN
     138              :                ! Smooth Coulomb Potential ::  Erf(x/rc)/x
     139          360 :                pot0_2(1, 1) = 2.0_dp/(rootpi*rc)
     140          360 :                pot0_2(2, 1) = 0.0_dp
     141          360 :                x = 0.0_dp
     142       507998 :                DO i = 2, np
     143       507638 :                   x = x + dx
     144       507638 :                   pot0_2(1, i) = erf(x/rc)/x
     145       507638 :                   t = 2._dp/(rootpi*x*rc)*EXP(-(x/rc)**2)
     146       507998 :                   pot0_2(2, i) = (t - pot0_2(1, i)/x)*dx
     147              :                END DO
     148           68 :             ELSEIF (qmmm_coupl_type == do_qmmm_swave) THEN
     149              :                ! S-wave expansion :: 1/x - exp(-2*x/rc) * ( 1/x - 1/rc )
     150           68 :                pot0_2(1, 1) = 1.0_dp/rc
     151           68 :                pot0_2(2, 1) = 0.0_dp
     152           68 :                x = 0.0_dp
     153       111180 :                DO i = 2, np
     154       111112 :                   x = x + dx
     155       111112 :                   t = EXP(-2.0_dp*x/rc)/rc
     156       111112 :                   pot0_2(1, i) = (1.0_dp - t*(rc + x))/x
     157       111180 :                   pot0_2(2, i) = ((t*(rc**2 + 2.0_dp*rc*x + 2.0_dp*x**2)/rc - 1.0_dp)/x**2)*dx
     158              :                END DO
     159              :             END IF
     160          428 :             pgf => pgfs(K)%pgf
     161          428 :             CPASSERT(pgf%Elp_Radius == rc)
     162          428 :             ig_start = 1
     163          428 :             IF (compatibility .AND. (qmmm_coupl_type == do_qmmm_gauss)) ig_start = 2
     164         3662 :             DO Ig = ig_start, pgf%number_of_gaussians
     165         3234 :                A = pgf%Ak(Ig)
     166         3234 :                G = pgf%Gk(Ig)
     167         3234 :                pot0_2(1, 1) = pot0_2(1, 1) - A
     168         3234 :                x = 0.0_dp
     169      4760776 :                DO i = 2, np
     170      4757114 :                   x = x + dx
     171      4757114 :                   t = EXP(-(x/G)**2)*A
     172      4757114 :                   t1 = 1/G**2
     173      4757114 :                   t2 = t1*t
     174      4757114 :                   pot0_2(1, i) = pot0_2(1, i) - t
     175      4760348 :                   pot0_2(2, i) = pot0_2(2, i) + 2.0_dp*x*t2*dx
     176              :                END DO
     177              :             END DO
     178              : 
     179              :             ! Print info on the unidimensional MM electrostatic potential
     180          428 :             IF (BTEST(cp_print_key_should_output(logger%iter_info, print_section, "MM_POTENTIAL") &
     181           52 :                       , cp_p_file)) THEN
     182          178 :                WRITE (rc_s, '(F6.3)') rc
     183              :                unit_nr = cp_print_key_unit_nr(logger, print_section, "MM_POTENTIAL", &
     184          178 :                                               extension="_rc="//TRIM(ADJUSTL(rc_s))//".data")
     185          178 :                IF (unit_nr > 0) THEN
     186           93 :                   WRITE (unit_nr, '(A)') "# MM ELECTROSTATIC POTENTIAL - UNIDIMENSIONAL - ATOMIC UNITS"
     187           93 :                   WRITE (unit_nr, '(A,I5)') "# MM ELECTROSTATIC POTENTIAL - Nr. of Gaussians:", pgf%number_of_gaussians
     188           93 :                   WRITE (unit_nr, '(A,T10,A,T30,A,T300,A)') "#", "Xval", "Gaussians", "LongRange"
     189           93 :                   myFormat = "T10,F15.9,T30,"
     190          881 :                   DO Ig = 1, pgf%number_of_gaussians
     191          788 :                      myind = INDEX(myFormat, " ")
     192          881 :                      WRITE (myFormat(myind:), '(A6)') "F12.9,"
     193              :                   END DO
     194           93 :                   myind = INDEX(myFormat, " ") - 1
     195           93 :                   myFormat = myFormat(1:myind)//"T300,F15.9"
     196           93 :                   myind = INDEX(myFormat, " ") - 1
     197           93 :                   x = 0.0_dp
     198       137156 :                   DO i = 1, np
     199              :                      WRITE (unit_nr, '('//myFormat(1:myind)//')') &
     200      1362883 :                         x, (EXP(-(x/pgf%Gk(Ig))**2)*pgf%Ak(Ig), Ig=1, pgf%number_of_gaussians), pot0_2(1, i)
     201       137156 :                      x = x + dx
     202              :                   END DO
     203              :                END IF
     204              :                CALL cp_print_key_finished_output(unit_nr, logger, print_section, &
     205          178 :                                                  "MM_POTENTIAL")
     206              :             END IF
     207              :          CASE DEFAULT
     208           52 :             DEALLOCATE (potentials(K)%Pot)
     209           52 :             IF (output_unit > 0) WRITE (output_unit, '(A)') " QMMM Potential - Spline Interpolation - not Initialized!"
     210          634 :             CYCLE Potential_Type
     211              :          END SELECT
     212              :          NULLIFY (mm_atom_index)
     213          582 :          ALLOCATE (mm_atom_index(1))
     214              :          ! Build mm_atom_index List
     215       360602 :          DO J = 1, SIZE(mm_el_pot_radius)
     216       360602 :             IF (rc .EQ. mm_el_pot_radius(J)) THEN
     217       131110 :                Ndim = SIZE(mm_atom_index)
     218       131110 :                mm_atom_index(Ndim) = J
     219       131110 :                CALL reallocate(mm_atom_index, 1, Ndim + 1)
     220              :             END IF
     221              :          END DO
     222          582 :          CALL reallocate(mm_atom_index, 1, Ndim)
     223              : 
     224          582 :          NULLIFY (potentials(K)%Pot%pot0_2)
     225              :          CALL qmmm_pot_type_create(potentials(K)%Pot, pot0_2=pot0_2, &
     226              :                                    Rmax=Rmax, Rmin=Rmin, dx=dx, Rc=rc, npts=np, &
     227          986 :                                    mm_atom_index=mm_atom_index)
     228              : 
     229              :       END DO Potential_Type
     230          404 :       DEALLOCATE (radius)
     231          404 :    END SUBROUTINE qmmm_potential_init
     232              : 
     233              : ! **************************************************************************************************
     234              : !> \brief Creates the qmmm_pot_type structure
     235              : !> \param Pot ...
     236              : !> \param pot0_2 ...
     237              : !> \param Rmax ...
     238              : !> \param Rmin ...
     239              : !> \param dx ...
     240              : !> \param npts ...
     241              : !> \param rc ...
     242              : !> \param mm_atom_index ...
     243              : !> \par History
     244              : !>      09.2004 created [tlaino]
     245              : !> \author Teodoro Laino
     246              : ! **************************************************************************************************
     247          582 :    SUBROUTINE qmmm_pot_type_create(Pot, pot0_2, Rmax, Rmin, dx, npts, rc, &
     248              :                                    mm_atom_index)
     249              :       TYPE(qmmm_Pot_Type), POINTER                       :: Pot
     250              :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: pot0_2
     251              :       REAL(KIND=dp), INTENT(IN)                          :: Rmax, Rmin, dx
     252              :       INTEGER, INTENT(IN)                                :: npts
     253              :       REAL(KIND=dp), INTENT(IN)                          :: Rc
     254              :       INTEGER, DIMENSION(:), POINTER                     :: mm_atom_index
     255              : 
     256          582 :       Pot%pot0_2 => pot0_2
     257          582 :       Pot%mm_atom_index => mm_atom_index
     258          582 :       Pot%Rmax = Rmax
     259          582 :       Pot%Rmin = Rmin
     260          582 :       Pot%Rc = rc
     261          582 :       Pot%dx = dx
     262          582 :       Pot%npts = npts
     263              : 
     264          582 :    END SUBROUTINE qmmm_pot_type_create
     265              : 
     266              : END MODULE qmmm_elpot
        

Generated by: LCOV version 2.0-1