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 : !> \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
|
