Line data    Source code

       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \author T.Laino
      10             : ! **************************************************************************************************
      11             : MODULE qs_ks_qmmm_methods
      12             :    USE cp_control_types,                ONLY: dft_control_type
      13             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      14             :                                               cp_logger_type
      15             :    USE cp_output_handling,              ONLY: cp_print_key_finished_output,&
      16             :                                               cp_print_key_unit_nr
      17             :    USE cube_utils,                      ONLY: cube_info_type,&
      18             :                                               init_cube_info
      19             :    USE d3_poly,                         ONLY: init_d3_poly_module
      20             :    USE input_constants,                 ONLY: do_qmmm_gauss,&
      21             :                                               do_qmmm_swave
      22             :    USE input_section_types,             ONLY: section_vals_type
      23             :    USE kinds,                           ONLY: dp
      24             :    USE pw_env_types,                    ONLY: pw_env_get,&
      25             :                                               pw_env_retain
      26             :    USE pw_methods,                      ONLY: pw_axpy,&
      27             :                                               pw_integral_ab
      28             :    USE pw_pool_types,                   ONLY: pw_pool_p_type,&
      29             :                                               pw_pool_type
      30             :    USE pw_types,                        ONLY: pw_r3d_rs_type
      31             :    USE qmmm_types_low,                  ONLY: qmmm_env_qm_type
      32             :    USE qs_environment_types,            ONLY: get_qs_env,&
      33             :                                               qs_environment_type,&
      34             :                                               set_qs_env
      35             :    USE qs_ks_qmmm_types,                ONLY: qs_ks_qmmm_env_type
      36             : #include "./base/base_uses.f90"
      37             : 
      38             :    IMPLICIT NONE
      39             : 
      40             :    PRIVATE
      41             : 
      42             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_ks_qmmm_methods'
      43             : 
      44             :    PUBLIC :: ks_qmmm_env_rebuild, qmmm_calculate_energy, &
      45             :              qmmm_modify_hartree_pot
      46             : 
      47             : CONTAINS
      48             : 
      49             : ! **************************************************************************************************
      50             : !> \brief Initialize the ks_qmmm_env
      51             : !> \param qs_env ...
      52             : !> \param qmmm_env ...
      53             : !> \par History
      54             : !>      05.2004 created [tlaino]
      55             : !> \author Teodoro Laino
      56             : ! **************************************************************************************************
      57        3802 :    SUBROUTINE ks_qmmm_env_rebuild(qs_env, qmmm_env)
      58             :       TYPE(qs_environment_type), OPTIONAL, POINTER       :: qs_env
      59             :       TYPE(qmmm_env_qm_type), POINTER                    :: qmmm_env
      60             : 
      61             :       TYPE(qs_ks_qmmm_env_type), POINTER                 :: ks_qmmm_env
      62             : 
      63        3802 :       NULLIFY (ks_qmmm_env)
      64             :       CALL get_qs_env(qs_env=qs_env, &
      65        3802 :                       ks_qmmm_env=ks_qmmm_env)
      66             : 
      67             :       !   *** allocate the ks_qmmm env if not allocated yet!**
      68        3802 :       IF (.NOT. ASSOCIATED(ks_qmmm_env)) THEN
      69         378 :          ALLOCATE (ks_qmmm_env)
      70             :          CALL qs_ks_qmmm_create(ks_qmmm_env=ks_qmmm_env, qs_env=qs_env, &
      71         378 :                                 qmmm_env=qmmm_env)
      72         378 :          CALL set_qs_env(qs_env=qs_env, ks_qmmm_env=ks_qmmm_env)
      73             :       END IF
      74        3802 :    END SUBROUTINE ks_qmmm_env_rebuild
      75             : 
      76             : ! **************************************************************************************************
      77             : !> \brief allocates and initializes the given ks_qmmm_env.
      78             : !> \param ks_qmmm_env the ks_qmmm env to be initialized
      79             : !> \param qs_env the qs environment
      80             : !> \param qmmm_env ...
      81             : !> \par History
      82             : !>      05.2004 created [tlaino]
      83             : !> \author Teodoro Laino
      84             : ! **************************************************************************************************
      85         378 :    SUBROUTINE qs_ks_qmmm_create(ks_qmmm_env, qs_env, qmmm_env)
      86             :       TYPE(qs_ks_qmmm_env_type), INTENT(OUT)             :: ks_qmmm_env
      87             :       TYPE(qs_environment_type), POINTER                 :: qs_env
      88             :       TYPE(qmmm_env_qm_type), POINTER                    :: qmmm_env
      89             : 
      90             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'qs_ks_qmmm_create'
      91             : 
      92             :       INTEGER                                            :: handle, igrid
      93         378 :       TYPE(cube_info_type), DIMENSION(:), POINTER        :: cube_info
      94         378 :       TYPE(pw_pool_p_type), DIMENSION(:), POINTER        :: pools
      95             :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
      96             : 
      97         378 :       CALL timeset(routineN, handle)
      98             : 
      99             :       NULLIFY (ks_qmmm_env%pw_env, &
     100         378 :                ks_qmmm_env%cube_info)
     101         378 :       NULLIFY (auxbas_pw_pool)
     102             :       CALL get_qs_env(qs_env=qs_env, &
     103         378 :                       pw_env=ks_qmmm_env%pw_env)
     104         378 :       CALL pw_env_get(ks_qmmm_env%pw_env, auxbas_pw_pool=auxbas_pw_pool)
     105         378 :       CALL pw_env_retain(ks_qmmm_env%pw_env)
     106             : 
     107         378 :       ks_qmmm_env%pc_ener = 0.0_dp
     108         378 :       ks_qmmm_env%n_evals = 0
     109             : 
     110         378 :       CALL auxbas_pw_pool%create_pw(ks_qmmm_env%v_qmmm_rspace)
     111             : 
     112         378 :       IF ((qmmm_env%qmmm_coupl_type == do_qmmm_gauss) .OR. (qmmm_env%qmmm_coupl_type == do_qmmm_swave)) THEN
     113         236 :          CALL init_d3_poly_module() ! a fairly arbitrary but sufficient spot to do this
     114         236 :          CALL pw_env_get(ks_qmmm_env%pw_env, pw_pools=pools)
     115        8748 :          ALLOCATE (cube_info(SIZE(pools)))
     116        1196 :          DO igrid = 1, SIZE(pools)
     117             :             CALL init_cube_info(cube_info(igrid), &
     118             :                                 pools(igrid)%pool%pw_grid%dr(:), &
     119             :                                 pools(igrid)%pool%pw_grid%dh(:, :), &
     120             :                                 pools(igrid)%pool%pw_grid%dh_inv(:, :), &
     121             :                                 pools(igrid)%pool%pw_grid%orthorhombic, &
     122        1196 :                                 qmmm_env%maxRadius(igrid))
     123             :          END DO
     124         236 :          ks_qmmm_env%cube_info => cube_info
     125             :       END IF
     126         378 :       NULLIFY (ks_qmmm_env%matrix_h)
     127             :       !
     128             : 
     129         378 :       CALL timestop(handle)
     130             : 
     131         378 :    END SUBROUTINE qs_ks_qmmm_create
     132             : 
     133             : ! **************************************************************************************************
     134             : !> \brief Computes the contribution to the total energy of the QM/MM
     135             : !>      electrostatic coupling
     136             : !> \param qs_env ...
     137             : !> \param rho ...
     138             : !> \param v_qmmm ...
     139             : !> \param qmmm_energy ...
     140             : !> \par History
     141             : !>      05.2004 created [tlaino]
     142             : !> \author Teodoro Laino
     143             : ! **************************************************************************************************
     144        6298 :    SUBROUTINE qmmm_calculate_energy(qs_env, rho, v_qmmm, qmmm_energy)
     145             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     146             :       TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(IN)     :: rho
     147             :       TYPE(pw_r3d_rs_type), INTENT(IN)                   :: v_qmmm
     148             :       REAL(KIND=dp), INTENT(INOUT)                       :: qmmm_energy
     149             : 
     150             :       CHARACTER(len=*), PARAMETER :: routineN = 'qmmm_calculate_energy'
     151             : 
     152             :       INTEGER                                            :: handle, ispin, output_unit
     153             :       TYPE(cp_logger_type), POINTER                      :: logger
     154             :       TYPE(dft_control_type), POINTER                    :: dft_control
     155             :       TYPE(pw_r3d_rs_type), POINTER                      :: rho0_s_rs
     156             :       TYPE(section_vals_type), POINTER                   :: input
     157             : 
     158        6298 :       CALL timeset(routineN, handle)
     159        6298 :       NULLIFY (dft_control, input, logger)
     160        6298 :       logger => cp_get_default_logger()
     161             : 
     162             :       CALL get_qs_env(qs_env=qs_env, &
     163             :                       dft_control=dft_control, &
     164        6298 :                       input=input)
     165             : 
     166             :       output_unit = cp_print_key_unit_nr(logger, input, "QMMM%PRINT%PROGRAM_RUN_INFO", &
     167        6298 :                                          extension=".qmmmLog")
     168        6298 :       IF (output_unit > 0) WRITE (UNIT=output_unit, FMT="(T3,A)") &
     169        3145 :          "Adding QM/MM electrostatic potential to the Kohn-Sham potential."
     170             : 
     171        6298 :       qmmm_energy = 0.0_dp
     172       12826 :       DO ispin = 1, dft_control%nspins
     173       12826 :          qmmm_energy = qmmm_energy + pw_integral_ab(rho(ispin), v_qmmm)
     174             :       END DO
     175        6298 :       IF (dft_control%qs_control%gapw) THEN
     176             :          CALL get_qs_env(qs_env=qs_env, &
     177        1446 :                          rho0_s_rs=rho0_s_rs)
     178        1446 :          CPASSERT(ASSOCIATED(rho0_s_rs))
     179        1446 :          qmmm_energy = qmmm_energy + pw_integral_ab(rho0_s_rs, v_qmmm)
     180             :       END IF
     181             : 
     182             :       CALL cp_print_key_finished_output(output_unit, logger, input, &
     183        6298 :                                         "QMMM%PRINT%PROGRAM_RUN_INFO")
     184             : 
     185        6298 :       CALL timestop(handle)
     186        6298 :    END SUBROUTINE qmmm_calculate_energy
     187             : 
     188             : ! **************************************************************************************************
     189             : !> \brief Modify the hartree potential in order to include the QM/MM correction
     190             : !> \param v_hartree ...
     191             : !> \param v_qmmm ...
     192             : !> \param scale ...
     193             : !> \par History
     194             : !>      05.2004 created [tlaino]
     195             : !> \author Teodoro Laino
     196             : ! **************************************************************************************************
     197        6624 :    SUBROUTINE qmmm_modify_hartree_pot(v_hartree, v_qmmm, scale)
     198             :       TYPE(pw_r3d_rs_type), INTENT(INOUT)                :: v_hartree
     199             :       TYPE(pw_r3d_rs_type), INTENT(IN)                   :: v_qmmm
     200             :       REAL(KIND=dp), INTENT(IN)                          :: scale
     201             : 
     202        6624 :       CALL pw_axpy(v_qmmm, v_hartree, v_qmmm%pw_grid%dvol*scale)
     203             : 
     204        6624 :    END SUBROUTINE qmmm_modify_hartree_pot
     205             : 
     206             : END MODULE qs_ks_qmmm_methods