1!--------------------------------------------------------------------------------------------------!
2!   CP2K: A general program to perform molecular dynamics simulations                              !
3!   Copyright (C) 2000 - 2019  CP2K developers group                                               !
4!--------------------------------------------------------------------------------------------------!
5
6! **************************************************************************************************
7!> \par History
8!>      09.2004 created [tlaino]
9!> \author Teodoro Laino
10! **************************************************************************************************
11MODULE qmmm_elpot
12   USE cell_types,                      ONLY: cell_type
13   USE cp_log_handling,                 ONLY: cp_get_default_logger,&
14                                              cp_logger_get_default_io_unit,&
15                                              cp_logger_type
16   USE cp_output_handling,              ONLY: cp_p_file,&
17                                              cp_print_key_finished_output,&
18                                              cp_print_key_should_output,&
19                                              cp_print_key_unit_nr
20   USE input_constants,                 ONLY: do_qmmm_coulomb,&
21                                              do_qmmm_gauss,&
22                                              do_qmmm_pcharge,&
23                                              do_qmmm_swave
24   USE input_section_types,             ONLY: section_vals_type
25   USE kinds,                           ONLY: default_path_length,&
26                                              default_string_length,&
27                                              dp
28   USE mathconstants,                   ONLY: rootpi
29   USE memory_utilities,                ONLY: reallocate
30   USE qmmm_gaussian_types,             ONLY: qmmm_gaussian_p_type,&
31                                              qmmm_gaussian_type
32   USE qmmm_types_low,                  ONLY: qmmm_Pot_Type,&
33                                              qmmm_pot_p_type
34#include "./base/base_uses.f90"
35
36   IMPLICIT NONE
37   PRIVATE
38
39   LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
40   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qmmm_elpot'
41   PUBLIC :: qmmm_potential_init
42
43CONTAINS
44
45! **************************************************************************************************
46!> \brief Initialize the QMMM potential stored on vector,
47!>      according the qmmm_coupl_type
48!> \param qmmm_coupl_type ...
49!> \param mm_el_pot_radius ...
50!> \param potentials ...
51!> \param pgfs ...
52!> \param mm_cell ...
53!> \param compatibility ...
54!> \param print_section ...
55!> \par History
56!>      09.2004 created [tlaino]
57!> \author Teodoro Laino
58! **************************************************************************************************
59   SUBROUTINE qmmm_potential_init(qmmm_coupl_type, mm_el_pot_radius, potentials, &
60                                  pgfs, mm_cell, compatibility, print_section)
61      INTEGER, INTENT(IN)                                :: qmmm_coupl_type
62      REAL(KIND=dp), DIMENSION(:), POINTER               :: mm_el_pot_radius
63      TYPE(qmmm_pot_p_type), DIMENSION(:), POINTER       :: potentials
64      TYPE(qmmm_gaussian_p_type), DIMENSION(:), POINTER  :: pgfs
65      TYPE(cell_type), POINTER                           :: mm_cell
66      LOGICAL, INTENT(IN)                                :: compatibility
67      TYPE(section_vals_type), POINTER                   :: print_section
68
69      CHARACTER(len=*), PARAMETER :: routineN = 'qmmm_potential_init', &
70         routineP = moduleN//':'//routineN
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      INTEGER, DIMENSION(:), POINTER                     :: mm_atom_index
78      LOGICAL                                            :: found
79      REAL(KIND=dp)                                      :: A, G, rc, Rmax, Rmin, t, t1, t2, x
80      REAL(KIND=dp), DIMENSION(:), POINTER               :: radius
81      REAL(KIND=dp), DIMENSION(:, :), POINTER            :: pot0_2
82      TYPE(cp_logger_type), POINTER                      :: logger
83      TYPE(qmmm_gaussian_type), POINTER                  :: pgf
84
85      logger => cp_get_default_logger()
86      output_unit = cp_logger_get_default_io_unit(logger)
87      Rmin = 0.0_dp
88      Rmax = SQRT(mm_cell%hmat(1, 1)**2 + &
89                  mm_cell%hmat(2, 2)**2 + &
90                  mm_cell%hmat(3, 3)**2)
91      np = CEILING(rmax/dx) + 1
92      !
93      ! Preprocessing
94      !
95      IF (SIZE(mm_el_pot_radius) /= 0) THEN
96         ALLOCATE (radius(1))
97         radius(1) = mm_el_pot_radius(1)
98      ELSE
99         ALLOCATE (radius(0))
100      END IF
101      Loop_on_all_values: DO I = 2, SIZE(mm_el_pot_radius)
102         Found = .FALSE.
103         Loop_on_found_values: DO J = 1, SIZE(radius)
104            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         IF (.NOT. Found) THEN
110            Ndim = SIZE(radius)
111            Ndim = Ndim + 1
112            CALL REALLOCATE(radius, 1, Ndim)
113            radius(Ndim) = mm_el_pot_radius(i)
114         END IF
115      END DO Loop_on_all_values
116      !
117      CPASSERT(.NOT. ASSOCIATED(potentials))
118      ALLOCATE (potentials(SIZE(radius)))
119
120      Potential_Type: DO K = 1, SIZE(radius)
121
122         rc = radius(K)
123         ALLOCATE (potentials(K)%Pot)
124         SELECT CASE (qmmm_coupl_type)
125         CASE (do_qmmm_coulomb)
126            NULLIFY (pot0_2)
127         CASE (do_qmmm_pcharge)
128            NULLIFY (pot0_2)
129         CASE (do_qmmm_gauss, do_qmmm_swave)
130            ALLOCATE (pot0_2(2, np))
131         END SELECT
132
133         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            IF (qmmm_coupl_type == do_qmmm_gauss) THEN
138               ! Smooth Coulomb Potential ::  Erf(x/rc)/x
139               pot0_2(1, 1) = 2.0_dp/(rootpi*rc)
140               pot0_2(2, 1) = 0.0_dp
141               x = 0.0_dp
142               DO i = 2, np
143                  x = x + dx
144                  pot0_2(1, i) = erf(x/rc)/x
145                  t = 2._dp/(rootpi*x*rc)*EXP(-(x/rc)**2)
146                  pot0_2(2, i) = (t - pot0_2(1, i)/x)*dx
147               END DO
148            ELSEIF (qmmm_coupl_type == do_qmmm_swave) THEN
149               ! S-wave expansion :: 1/x - exp(-2*x/rc) * ( 1/x - 1/rc )
150               pot0_2(1, 1) = 1.0_dp/rc
151               pot0_2(2, 1) = 0.0_dp
152               x = 0.0_dp
153               DO i = 2, np
154                  x = x + dx
155                  t = EXP(-2.0_dp*x/rc)/rc
156                  pot0_2(1, i) = (1.0_dp - t*(rc + x))/x
157                  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            pgf => pgfs(K)%pgf
161            CPASSERT(pgf%Elp_Radius == rc)
162            ig_start = 1
163            IF (compatibility .AND. (qmmm_coupl_type == do_qmmm_gauss)) ig_start = 2
164            DO Ig = ig_start, pgf%number_of_gaussians
165               A = pgf%Ak(Ig)
166               G = pgf%Gk(Ig)
167               pot0_2(1, 1) = pot0_2(1, 1) - A
168               x = 0.0_dp
169               DO i = 2, np
170                  x = x + dx
171                  t = EXP(-(x/G)**2)*A
172                  t1 = 1/G**2
173                  t2 = t1*t
174                  pot0_2(1, i) = pot0_2(1, i) - t
175                  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            IF (BTEST(cp_print_key_should_output(logger%iter_info, print_section, "MM_POTENTIAL") &
181                      , cp_p_file)) THEN
182               WRITE (rc_s, '(F6.3)') rc
183               unit_nr = cp_print_key_unit_nr(logger, print_section, "MM_POTENTIAL", &
184                                              extension="_rc="//TRIM(ADJUSTL(rc_s))//".data")
185               IF (unit_nr > 0) THEN
186                  WRITE (unit_nr, '(A)') "# MM ELECTROSTATIC POTENTIAL - UNIDIMENSIONAL - ATOMIC UNITS"
187                  WRITE (unit_nr, '(A,I5)') "# MM ELECTROSTATIC POTENTIAL - Nr. of Gaussians:", pgf%number_of_gaussians
188                  WRITE (unit_nr, '(A,T10,A,T30,A,T300,A)') "#", "Xval", "Gaussians", "LongRange"
189                  myFormat = "T10,F15.9,T30,"
190                  DO Ig = 1, pgf%number_of_gaussians
191                     myind = INDEX(myFormat, " ")
192                     WRITE (myFormat(myind:), '(A6)') "F12.9,"
193                  END DO
194                  myind = INDEX(myFormat, " ") - 1
195                  myFormat = myFormat(1:myind)//"T300,F15.9"
196                  myind = INDEX(myFormat, " ") - 1
197                  x = 0.0_dp
198                  DO i = 1, np
199                     WRITE (unit_nr, '('//myFormat(1:myind)//')') &
200                        x, (EXP(-(x/pgf%Gk(Ig))**2)*pgf%Ak(Ig), Ig=1, pgf%number_of_gaussians), pot0_2(1, i)
201                     x = x + dx
202                  END DO
203               END IF
204               CALL cp_print_key_finished_output(unit_nr, logger, print_section, &
205                                                 "MM_POTENTIAL")
206            END IF
207         CASE DEFAULT
208            DEALLOCATE (potentials(K)%Pot)
209            IF (output_unit > 0) WRITE (output_unit, '(A)') " QMMM Potential - Spline Interpolation - not Initialized!"
210            CYCLE Potential_Type
211         END SELECT
212         NULLIFY (mm_atom_index)
213         ALLOCATE (mm_atom_index(1))
214         ! Build mm_atom_index List
215         DO J = 1, SIZE(mm_el_pot_radius)
216            IF (rc .EQ. mm_el_pot_radius(J)) THEN
217               Ndim = SIZE(mm_atom_index)
218               mm_atom_index(Ndim) = J
219               CALL reallocate(mm_atom_index, 1, Ndim + 1)
220            ENDIF
221         END DO
222         CALL reallocate(mm_atom_index, 1, Ndim)
223
224         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                                   mm_atom_index=mm_atom_index)
228
229      END DO Potential_Type
230      DEALLOCATE (radius)
231   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   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      Pot%pot0_2 => pot0_2
257      Pot%mm_atom_index => mm_atom_index
258      Pot%Rmax = Rmax
259      Pot%Rmin = Rmin
260      Pot%Rc = rc
261      Pot%dx = dx
262      Pot%npts = npts
263
264   END SUBROUTINE qmmm_pot_type_create
265
266END MODULE qmmm_elpot
267