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