1!--------------------------------------------------------------------------------------------------!
2!   CP2K: A general program to perform molecular dynamics simulations                              !
3!   Copyright (C) 2000 - 2019  CP2K developers group                                               !
4!--------------------------------------------------------------------------------------------------!
5
6! **************************************************************************************************
7!> \brief Initialize the use of the gaussians to treat the QMMM
8!>      coupling potential
9!> \par History
10!>      6.2004 created [tlaino]
11!> \author Teodoro Laino
12! **************************************************************************************************
13MODULE qmmm_gaussian_init
14   USE ao_util,                         ONLY: exp_radius
15   USE cp_log_handling,                 ONLY: cp_get_default_logger,&
16                                              cp_logger_type
17   USE cp_output_handling,              ONLY: cp_print_key_finished_output,&
18                                              cp_print_key_unit_nr
19   USE cp_para_types,                   ONLY: cp_para_env_type
20   USE gaussian_gridlevels,             ONLY: gaussian_gridlevel,&
21                                              gridlevel_info_type
22   USE input_constants,                 ONLY: do_qmmm_gauss,&
23                                              do_qmmm_swave
24   USE input_section_types,             ONLY: section_vals_type,&
25                                              section_vals_val_get
26   USE kinds,                           ONLY: dp
27   USE memory_utilities,                ONLY: reallocate
28   USE pw_env_types,                    ONLY: pw_env_get,&
29                                              pw_env_type
30   USE pw_pool_types,                   ONLY: pw_pool_p_type
31   USE qmmm_gaussian_data,              ONLY: max_geep_lib_gauss,&
32                                              min_geep_lib_gauss
33   USE qmmm_gaussian_input,             ONLY: read_mm_potential,&
34                                              set_mm_potential_erf,&
35                                              set_mm_potential_swave
36   USE qmmm_gaussian_types,             ONLY: qmmm_gaussian_p_type,&
37                                              qmmm_gaussian_type
38#include "./base/base_uses.f90"
39
40   IMPLICIT NONE
41   PRIVATE
42
43   LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
44   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qmmm_gaussian_init'
45   PUBLIC :: qmmm_gaussian_initialize
46
47CONTAINS
48
49! **************************************************************************************************
50!> \brief Initialize the Gaussian QMMM Environment
51!> \param qmmm_gaussian_fns ...
52!> \param para_env ...
53!> \param pw_env ...
54!> \param mm_el_pot_radius ...
55!> \param mm_el_pot_radius_corr ...
56!> \param qmmm_coupl_type ...
57!> \param eps_mm_rspace ...
58!> \param maxradius ...
59!> \param maxchrg ...
60!> \param compatibility ...
61!> \param print_section ...
62!> \param qmmm_section ...
63!> \par History
64!>      06.2004 created [tlaino]
65!> \author Teodoro Laino
66! **************************************************************************************************
67   SUBROUTINE qmmm_gaussian_initialize(qmmm_gaussian_fns, para_env, pw_env, &
68                                       mm_el_pot_radius, mm_el_pot_radius_corr, &
69                                       qmmm_coupl_type, eps_mm_rspace, maxradius, maxchrg, compatibility, &
70                                       print_section, qmmm_section)
71      TYPE(qmmm_gaussian_p_type), DIMENSION(:), POINTER  :: qmmm_gaussian_fns
72      TYPE(cp_para_env_type), POINTER                    :: para_env
73      TYPE(pw_env_type), POINTER                         :: pw_env
74      REAL(KIND=dp), DIMENSION(:), POINTER               :: mm_el_pot_radius, mm_el_pot_radius_corr
75      INTEGER, INTENT(IN)                                :: qmmm_coupl_type
76      REAL(KIND=dp), INTENT(IN)                          :: eps_mm_rspace
77      REAL(KIND=dp), DIMENSION(:), POINTER               :: maxradius
78      REAL(KIND=dp), INTENT(IN)                          :: maxchrg
79      LOGICAL, INTENT(IN)                                :: compatibility
80      TYPE(section_vals_type), POINTER                   :: print_section, qmmm_section
81
82      CHARACTER(len=*), PARAMETER :: routineN = 'qmmm_gaussian_initialize', &
83         routineP = moduleN//':'//routineN
84
85      INTEGER                                            :: i, ilevel, j, Ndim, num_geep_gauss, &
86                                                            output_unit
87      LOGICAL                                            :: Found, use_geep_lib
88      REAL(KIND=dp)                                      :: alpha, mymaxradius, Prefactor
89      REAL(KIND=dp), DIMENSION(:), POINTER               :: c_radius, radius
90      TYPE(cp_logger_type), POINTER                      :: logger
91      TYPE(gridlevel_info_type), POINTER                 :: gridlevel_info
92      TYPE(pw_pool_p_type), DIMENSION(:), POINTER        :: pools
93      TYPE(qmmm_gaussian_type), POINTER                  :: mypgf
94
95! Statements
96
97      NULLIFY (mypgf, gridlevel_info, radius, c_radius, logger)
98      logger => cp_get_default_logger()
99      CALL section_vals_val_get(qmmm_section, "USE_GEEP_LIB", i_val=num_geep_gauss)
100      IF (num_geep_gauss == 0) THEN
101         use_geep_lib = .FALSE.
102      ELSE
103         use_geep_lib = .TRUE.
104         CPASSERT(num_geep_gauss >= min_geep_lib_gauss)
105         CPASSERT(num_geep_gauss <= max_geep_lib_gauss)
106      END IF
107      SELECT CASE (qmmm_coupl_type)
108      CASE (do_qmmm_gauss, do_qmmm_swave)
109         !
110         ! Preprocessing...
111         !
112         ALLOCATE (radius(1))
113         ALLOCATE (c_radius(1))
114         Ndim = SIZE(radius)
115         Loop_on_all_values: DO I = 1, SIZE(mm_el_pot_radius)
116            Found = .FALSE.
117            Loop_on_found_values: DO J = 1, SIZE(radius) - 1
118               IF (mm_el_pot_radius(i) .EQ. radius(j)) THEN
119                  Found = .TRUE.
120                  EXIT Loop_on_found_values
121               END IF
122            END DO Loop_on_found_values
123            IF (.NOT. Found) THEN
124               Ndim = SIZE(radius)
125               radius(Ndim) = mm_el_pot_radius(i)
126               c_radius(Ndim) = mm_el_pot_radius_corr(i)
127               Ndim = Ndim + 1
128               CALL REALLOCATE(radius, 1, Ndim)
129               CALL REALLOCATE(c_radius, 1, Ndim)
130            END IF
131         END DO Loop_on_all_values
132         !
133         IF (Ndim - 1 > 0) THEN
134            CALL REALLOCATE(radius, 1, Ndim - 1)
135            CALL REALLOCATE(c_radius, 1, Ndim - 1)
136         ELSE IF (Ndim - 1 == 0) THEN
137            DEALLOCATE (radius)
138            DEALLOCATE (c_radius)
139         ELSE
140            CPABORT("")
141         END IF
142         !
143         ALLOCATE (qmmm_gaussian_fns(Ndim - 1))
144         DO I = 1, Ndim - 1
145            NULLIFY (qmmm_gaussian_fns(I)%pgf)
146            ALLOCATE (qmmm_gaussian_fns(I)%pgf)
147            NULLIFY (qmmm_gaussian_fns(I)%pgf%Ak)
148            NULLIFY (qmmm_gaussian_fns(I)%pgf%Gk)
149            NULLIFY (qmmm_gaussian_fns(I)%pgf%grid_level)
150            !
151            ! Default Values
152            !
153            qmmm_gaussian_fns(I)%pgf%Elp_Radius = radius(I)
154            qmmm_gaussian_fns(I)%pgf%Elp_Radius_corr = c_radius(I)
155         END DO
156         IF (ASSOCIATED(radius)) THEN
157            DEALLOCATE (radius)
158         END IF
159         IF (ASSOCIATED(c_radius)) THEN
160            DEALLOCATE (c_radius)
161         END IF
162         !
163         IF (use_geep_lib) THEN
164            IF (qmmm_coupl_type == do_qmmm_gauss) THEN
165               CALL set_mm_potential_erf(qmmm_gaussian_fns, &
166                                         compatibility, num_geep_gauss)
167            ELSEIF (qmmm_coupl_type == do_qmmm_swave) THEN
168               CALL set_mm_potential_swave(qmmm_gaussian_fns, &
169                                           num_geep_gauss)
170            END IF
171         ELSE
172            CALL read_mm_potential(para_env, qmmm_gaussian_fns, &
173                                   (compatibility .AND. (qmmm_coupl_type == do_qmmm_gauss)), qmmm_section)
174         END IF
175         !
176         CALL pw_env_get(pw_env, pw_pools=pools, gridlevel_info=gridlevel_info)
177         ALLOCATE (maxradius(SIZE(pools)))
178         maxradius = 0.0_dp
179         DO J = 1, SIZE(qmmm_gaussian_fns)
180            mypgf => qmmm_gaussian_fns(J)%pgf
181            ALLOCATE (mypgf%grid_level(SIZE(mypgf%Ak)))
182            mypgf%grid_level = 0
183            mymaxradius = 0.0_dp
184            DO I = 1, mypgf%number_of_gaussians
185               IF (mypgf%Gk(I) /= 0.0_dp) THEN
186                  alpha = 1.0_dp/mypgf%Gk(I)
187                  alpha = alpha*alpha
188                  ilevel = gaussian_gridlevel(gridlevel_info, alpha)
189                  Prefactor = mypgf%Ak(I)*maxchrg
190                  mymaxradius = exp_radius(0, alpha, eps_mm_rspace, Prefactor)
191                  maxradius(ilevel) = MAX(maxradius(ilevel), mymaxradius)
192                  mypgf%grid_level(i) = ilevel
193               END IF
194            END DO
195         END DO
196         !
197         ! End of gaussian initialization...
198      CASE DEFAULT
199         output_unit = cp_print_key_unit_nr(logger, print_section, "PROGRAM_RUN_INFO", &
200                                            extension=".qmmmLog")
201         IF (output_unit > 0) WRITE (output_unit, '(A)') " QMMM Gaussian Data Not Initialized!"
202         CALL cp_print_key_finished_output(output_unit, logger, print_section, "PROGRAM_RUN_INFO")
203      END SELECT
204   END SUBROUTINE qmmm_gaussian_initialize
205
206END MODULE qmmm_gaussian_init
207