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 : !> \brief Setting up the potential for QM/MM periodic boundary conditions calculations
10 : !> \par History
11 : !> 07.2005 created [tlaino]
12 : !> \author Teodoro Laino
13 : ! **************************************************************************************************
14 : MODULE qmmm_per_elpot
15 : USE ao_util, ONLY: exp_radius
16 : USE cell_types, ONLY: cell_type
17 : USE cp_log_handling, ONLY: cp_get_default_logger,&
18 : cp_logger_get_default_io_unit,&
19 : cp_logger_type
20 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
21 : cp_print_key_unit_nr
22 : USE ewald_environment_types, ONLY: ewald_env_create,&
23 : ewald_env_get,&
24 : ewald_env_set,&
25 : ewald_environment_type,&
26 : read_ewald_section
27 : USE ewald_pw_types, ONLY: ewald_pw_create,&
28 : ewald_pw_type
29 : USE ewald_spline_util, ONLY: Setup_Ewald_Spline
30 : USE input_constants, ONLY: do_qmmm_coulomb,&
31 : do_qmmm_gauss,&
32 : do_qmmm_pcharge,&
33 : do_qmmm_swave
34 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
35 : section_vals_type,&
36 : section_vals_val_get
37 : USE kinds, ONLY: dp
38 : USE mathconstants, ONLY: pi
39 : USE message_passing, ONLY: mp_para_env_type
40 : USE pw_poisson_types, ONLY: do_ewald_ewald,&
41 : do_ewald_none,&
42 : do_ewald_pme,&
43 : do_ewald_spme
44 : USE qmmm_gaussian_types, ONLY: qmmm_gaussian_p_type,&
45 : qmmm_gaussian_type
46 : USE qmmm_types_low, ONLY: qmmm_per_pot_p_type,&
47 : qmmm_per_pot_type,&
48 : qmmm_pot_p_type
49 : #include "./base/base_uses.f90"
50 :
51 : IMPLICIT NONE
52 : PRIVATE
53 :
54 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
55 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qmmm_per_elpot'
56 : PUBLIC :: qmmm_per_potential_init, qmmm_ewald_potential_init
57 :
58 : CONTAINS
59 :
60 : ! **************************************************************************************************
61 : !> \brief Initialize the QMMM potential stored on vector,
62 : !> according the qmmm_coupl_type
63 : !> \param qmmm_coupl_type ...
64 : !> \param per_potentials ...
65 : !> \param potentials ...
66 : !> \param pgfs ...
67 : !> \param qm_cell_small ...
68 : !> \param mm_cell ...
69 : !> \param compatibility ...
70 : !> \param qmmm_periodic ...
71 : !> \param print_section ...
72 : !> \param eps_mm_rspace ...
73 : !> \param maxchrg ...
74 : !> \param ncp ...
75 : !> \param ncpl ...
76 : !> \par History
77 : !> 09.2004 created [tlaino]
78 : !> \author Teodoro Laino
79 : ! **************************************************************************************************
80 40 : SUBROUTINE qmmm_per_potential_init(qmmm_coupl_type, per_potentials, potentials, &
81 : pgfs, qm_cell_small, mm_cell, compatibility, qmmm_periodic, print_section, &
82 : eps_mm_rspace, maxchrg, ncp, ncpl)
83 : INTEGER, INTENT(IN) :: qmmm_coupl_type
84 : TYPE(qmmm_per_pot_p_type), DIMENSION(:), POINTER :: per_potentials
85 : TYPE(qmmm_pot_p_type), DIMENSION(:), POINTER :: potentials
86 : TYPE(qmmm_gaussian_p_type), DIMENSION(:), POINTER :: pgfs
87 : TYPE(cell_type), POINTER :: qm_cell_small, mm_cell
88 : LOGICAL, INTENT(IN) :: compatibility
89 : TYPE(section_vals_type), POINTER :: qmmm_periodic, print_section
90 : REAL(KIND=dp), INTENT(IN) :: eps_mm_rspace, maxchrg
91 : INTEGER, INTENT(IN) :: ncp(3), ncpl(3)
92 :
93 : INTEGER :: I, idim, ig, ig_start, iw, ix, iy, iz, &
94 : K, Kmax(3), n_rep_real(3), &
95 : n_rep_real_val, ncoarsel, ncoarset, &
96 : Ndim, output_unit
97 40 : INTEGER, DIMENSION(:), POINTER :: mm_atom_index
98 : REAL(KIND=dp) :: Ak, alpha, box(3), Fac(3), fs, g, g2, &
99 : Gk, Gmax, mymaxradius, npl, npt, &
100 : Prefactor, rc, rc2, Rmax, tmp, vec(3), &
101 : vol
102 40 : REAL(KIND=dp), DIMENSION(:), POINTER :: gx, gy, gz, Lg
103 : TYPE(cp_logger_type), POINTER :: logger
104 : TYPE(qmmm_gaussian_type), POINTER :: pgf
105 :
106 40 : NULLIFY (Lg, gx, gy, gz)
107 160 : ncoarset = PRODUCT(ncp)
108 160 : ncoarsel = PRODUCT(ncpl)
109 40 : logger => cp_get_default_logger()
110 40 : output_unit = cp_logger_get_default_io_unit(logger)
111 : Rmax = SQRT(mm_cell%hmat(1, 1)**2 + &
112 : mm_cell%hmat(2, 2)**2 + &
113 : mm_cell%hmat(3, 3)**2)
114 40 : CALL section_vals_val_get(qmmm_periodic, "GMAX", r_val=Gmax)
115 40 : CALL section_vals_val_get(qmmm_periodic, "REPLICA", i_val=n_rep_real_val)
116 160 : fac = 2.0e0_dp*Pi/(/mm_cell%hmat(1, 1), mm_cell%hmat(2, 2), mm_cell%hmat(3, 3)/)
117 160 : Kmax = CEILING(Gmax/Fac)
118 : Vol = mm_cell%hmat(1, 1)* &
119 : mm_cell%hmat(2, 2)* &
120 40 : mm_cell%hmat(3, 3)
121 40 : Ndim = (Kmax(1) + 1)*(2*Kmax(2) + 1)*(2*Kmax(3) + 1)
122 40 : ig_start = 1
123 160 : n_rep_real = n_rep_real_val
124 40 : IF (compatibility .AND. (qmmm_coupl_type == do_qmmm_gauss)) ig_start = 2
125 :
126 40 : CPASSERT(.NOT. ASSOCIATED(per_potentials))
127 186 : ALLOCATE (per_potentials(SIZE(pgfs)))
128 40 : CPASSERT(SIZE(pgfs) == SIZE(potentials))
129 106 : Potential_Type: DO K = 1, SIZE(pgfs)
130 :
131 66 : rc = pgfs(K)%pgf%Elp_Radius
132 660 : ALLOCATE (per_potentials(K)%Pot)
133 66 : SELECT CASE (qmmm_coupl_type)
134 : CASE (do_qmmm_coulomb, do_qmmm_pcharge)
135 : ! Not yet implemented for this case
136 0 : CPABORT("")
137 : CASE (do_qmmm_gauss, do_qmmm_swave)
138 198 : ALLOCATE (Lg(Ndim))
139 132 : ALLOCATE (gx(Ndim))
140 132 : ALLOCATE (gy(Ndim))
141 198 : ALLOCATE (gz(Ndim))
142 : END SELECT
143 :
144 62796 : LG = 0.0_dp
145 62796 : gx = 0.0_dp
146 62796 : gy = 0.0_dp
147 62796 : gz = 0.0_dp
148 :
149 0 : SELECT CASE (qmmm_coupl_type)
150 : CASE (do_qmmm_coulomb, do_qmmm_pcharge)
151 : ! Not yet implemented for this case
152 0 : CPABORT("")
153 : CASE (do_qmmm_gauss, do_qmmm_swave)
154 66 : pgf => pgfs(K)%pgf
155 66 : idim = 0
156 412 : DO ix = 0, kmax(1)
157 4654 : DO iy = -kmax(2), kmax(2)
158 67318 : DO iz = -kmax(3), kmax(3)
159 62730 : idim = idim + 1
160 62730 : IF (ix == 0 .AND. iy == 0 .AND. iz == 0) THEN
161 418 : DO Ig = ig_start, pgf%number_of_gaussians
162 352 : Gk = pgf%Gk(Ig)
163 352 : Ak = pgf%Ak(Ig)*Pi**(3.0_dp/2.0_dp)*Gk**3.0_dp
164 418 : LG(idim) = LG(idim) - Ak
165 : END DO
166 : ELSE
167 62664 : fs = 2.0_dp; IF (ix == 0) fs = 1.0_dp
168 250656 : vec = fac*(/REAL(ix, KIND=dp), REAL(iy, KIND=dp), REAL(iz, KIND=dp)/)
169 250656 : g2 = DOT_PRODUCT(vec, vec)
170 62664 : rc2 = rc*rc
171 62664 : g = SQRT(g2)
172 62664 : IF (qmmm_coupl_type == do_qmmm_gauss) THEN
173 62664 : LG(idim) = 4.0_dp*Pi/g2*EXP(-(g2*rc2)/4.0_dp)
174 0 : ELSEIF (qmmm_coupl_type == do_qmmm_swave) THEN
175 0 : tmp = 4.0_dp/rc2
176 0 : LG(idim) = 4.0_dp*Pi*tmp**2/(g2*(g2 + tmp)**2)
177 : END IF
178 386924 : DO Ig = ig_start, pgf%number_of_gaussians
179 324260 : Gk = pgf%Gk(Ig)
180 324260 : Ak = pgf%Ak(Ig)*Pi**(3.0_dp/2.0_dp)*Gk**3.0_dp
181 386924 : LG(idim) = LG(idim) - Ak*EXP(-(g*Gk)**2.0_dp/4.0_dp)
182 : END DO
183 : END IF
184 62730 : LG(idim) = fs*LG(idim)*1.0_dp/Vol
185 62730 : gx(idim) = fac(1)*REAL(ix, KIND=dp)
186 62730 : gy(idim) = fac(2)*REAL(iy, KIND=dp)
187 66972 : gz(idim) = fac(3)*REAL(iz, KIND=dp)
188 : END DO
189 : END DO
190 : END DO
191 :
192 186 : IF (ALL(n_rep_real == -1)) THEN
193 40 : mymaxradius = 0.0_dp
194 276 : DO I = 1, pgf%number_of_gaussians
195 276 : IF (pgf%Gk(I) /= 0.0_dp) THEN
196 236 : alpha = 1.0_dp/pgf%Gk(I)
197 236 : alpha = alpha*alpha
198 236 : Prefactor = pgf%Ak(I)*maxchrg
199 236 : mymaxradius = MAX(mymaxradius, exp_radius(0, alpha, eps_mm_rspace, Prefactor, rlow=mymaxradius))
200 : END IF
201 : END DO
202 40 : box(1) = (qm_cell_small%hmat(1, 1) - mm_cell%hmat(1, 1))/2.0_dp
203 40 : box(2) = (qm_cell_small%hmat(2, 2) - mm_cell%hmat(2, 2))/2.0_dp
204 40 : box(3) = (qm_cell_small%hmat(3, 3) - mm_cell%hmat(3, 3))/2.0_dp
205 160 : IF (ANY(box > 0.0_dp)) THEN
206 0 : CPABORT("")
207 : END IF
208 40 : n_rep_real(1) = CEILING((box(1) + mymaxradius)/mm_cell%hmat(1, 1))
209 40 : n_rep_real(2) = CEILING((box(2) + mymaxradius)/mm_cell%hmat(2, 2))
210 40 : n_rep_real(3) = CEILING((box(3) + mymaxradius)/mm_cell%hmat(3, 3))
211 : END IF
212 :
213 : CASE DEFAULT
214 0 : DEALLOCATE (per_potentials(K)%Pot)
215 0 : IF (output_unit > 0) WRITE (output_unit, '(A)') " QMMM Periodic Potential - not Initialized!"
216 66 : CYCLE Potential_Type
217 : END SELECT
218 :
219 : NULLIFY (mm_atom_index)
220 198 : ALLOCATE (mm_atom_index(SIZE(potentials(K)%pot%mm_atom_index)))
221 1606 : mm_atom_index = potentials(K)%pot%mm_atom_index
222 :
223 66 : NULLIFY (per_potentials(K)%Pot%LG, per_potentials(K)%Pot%mm_atom_index, &
224 66 : per_potentials(K)%Pot%gx, per_potentials(K)%Pot%gy, per_potentials(K)%Pot%gz)
225 : CALL qmmm_per_pot_type_create(per_potentials(K)%Pot, LG=LG, gx=gx, gy=gy, gz=gz, &
226 : Gmax=Gmax, Kmax=Kmax, n_rep_real=n_rep_real, &
227 : Fac=Fac, mm_atom_index=mm_atom_index, &
228 : mm_cell=mm_cell, &
229 66 : qmmm_per_section=qmmm_periodic, print_section=print_section)
230 :
231 : iw = cp_print_key_unit_nr(logger, print_section, "PERIODIC_INFO", &
232 66 : extension=".log")
233 66 : IF (iw > 0) THEN
234 37 : npt = REAL(ncoarset, KIND=dp)*REAL(ndim, KIND=dp)*REAL(SIZE(mm_atom_index), KIND=dp)
235 37 : npl = REAL(ncoarsel, KIND=dp)*REAL(ndim, KIND=dp)*REAL(SIZE(mm_atom_index), KIND=dp)
236 37 : WRITE (UNIT=iw, FMT="(/,T2,A)") REPEAT("-", 79)
237 37 : WRITE (UNIT=iw, FMT="(T2,A,T20,A,T80,A)") "-", "QMMM PERIODIC BOUNDARY CONDITION INFO", "-"
238 37 : WRITE (UNIT=iw, FMT="(T2,A)") REPEAT("-", 79)
239 37 : WRITE (UNIT=iw, FMT="(T2,A,T10,A,F15.6,T50,A,3I5,T80,A)") "-", "RADIUS =", rc, "REPLICA =", n_rep_real, "-"
240 34343 : WRITE (UNIT=iw, FMT="(T2,A,T10,A,F15.6,T50,A,I15,T80,A)") "-", "MINGVAL =", MINVAL(ABS(Lg)), &
241 74 : "GPOINTS =", ndim, "-"
242 37 : WRITE (UNIT=iw, FMT="(T2,A,T10,A,3I5,T50,A,3I5,T80,A)") "-", "NCOARSL =", ncpl, &
243 74 : "NCOARST =", ncp, "-"
244 37 : WRITE (UNIT=iw, FMT="(T2,A,T10,A,F15.0,T50,A,F15.0,T80,A)") "-", "NFLOP-L ~", npl, &
245 74 : "NFLOP-T ~", npt, "-"
246 37 : WRITE (UNIT=iw, FMT="(T2,A)") REPEAT("-", 79)
247 : END IF
248 : CALL cp_print_key_finished_output(iw, logger, print_section, &
249 106 : "PERIODIC_INFO")
250 :
251 : END DO Potential_Type
252 :
253 80 : END SUBROUTINE qmmm_per_potential_init
254 :
255 : ! **************************************************************************************************
256 : !> \brief Creates the qmmm_pot_type structure
257 : !> \param Pot ...
258 : !> \param LG ...
259 : !> \param gx ...
260 : !> \param gy ...
261 : !> \param gz ...
262 : !> \param GMax ...
263 : !> \param Kmax ...
264 : !> \param n_rep_real ...
265 : !> \param Fac ...
266 : !> \param mm_atom_index ...
267 : !> \param mm_cell ...
268 : !> \param qmmm_per_section ...
269 : !> \param print_section ...
270 : !> \par History
271 : !> 09.2004 created [tlaino]
272 : !> \author Teodoro Laino
273 : ! **************************************************************************************************
274 66 : SUBROUTINE qmmm_per_pot_type_create(Pot, LG, gx, gy, gz, GMax, Kmax, n_rep_real, &
275 : Fac, mm_atom_index, mm_cell, qmmm_per_section, print_section)
276 : TYPE(qmmm_per_pot_type), POINTER :: Pot
277 : REAL(KIND=dp), DIMENSION(:), POINTER :: LG, gx, gy, gz
278 : REAL(KIND=dp), INTENT(IN) :: Gmax
279 : INTEGER, INTENT(IN) :: Kmax(3), n_rep_real(3)
280 : REAL(KIND=dp), INTENT(IN) :: Fac(3)
281 : INTEGER, DIMENSION(:), POINTER :: mm_atom_index
282 : TYPE(cell_type), POINTER :: mm_cell
283 : TYPE(section_vals_type), POINTER :: qmmm_per_section, print_section
284 :
285 : INTEGER :: npts(3)
286 66 : INTEGER, DIMENSION(:), POINTER :: ngrids
287 : REAL(KIND=dp) :: hmat(3, 3)
288 : TYPE(section_vals_type), POINTER :: grid_print_section
289 :
290 66 : Pot%LG => LG
291 66 : Pot%gx => gx
292 66 : Pot%gy => gy
293 66 : Pot%gz => gz
294 66 : Pot%mm_atom_index => mm_atom_index
295 66 : Pot%Gmax = Gmax
296 264 : Pot%Kmax = Kmax
297 264 : Pot%n_rep_real = n_rep_real
298 264 : Pot%Fac = Fac
299 : !
300 : ! Setting Up Fit Procedure
301 : !
302 66 : NULLIFY (Pot%pw_grid)
303 66 : NULLIFY (Pot%pw_pool)
304 66 : NULLIFY (Pot%TabLR, ngrids)
305 66 : CALL section_vals_val_get(qmmm_per_section, "ngrids", i_vals=ngrids)
306 264 : npts = ngrids
307 858 : hmat = mm_cell%hmat
308 :
309 66 : grid_print_section => section_vals_get_subs_vals(print_section, "GRID_INFORMATION")
310 : CALL Setup_Ewald_Spline(pw_grid=Pot%pw_grid, pw_pool=Pot%pw_pool, coeff=Pot%TabLR, &
311 : LG=LG, gx=gx, gy=gy, gz=gz, hmat=hmat, npts=npts, param_section=qmmm_per_section, &
312 66 : tag="qmmm", print_section=grid_print_section)
313 :
314 66 : END SUBROUTINE qmmm_per_pot_type_create
315 :
316 : ! **************************************************************************************************
317 : !> \brief Initialize the QMMM Ewald potential needed for QM-QM Coupling using
318 : !> point charges
319 : !> \param ewald_env ...
320 : !> \param ewald_pw ...
321 : !> \param qmmm_coupl_type ...
322 : !> \param mm_cell ...
323 : !> \param para_env ...
324 : !> \param qmmm_periodic ...
325 : !> \param print_section ...
326 : !> \par History
327 : !> 10.2014 created [JGH]
328 : !> \author JGH
329 : ! **************************************************************************************************
330 32 : SUBROUTINE qmmm_ewald_potential_init(ewald_env, ewald_pw, qmmm_coupl_type, mm_cell, para_env, &
331 : qmmm_periodic, print_section)
332 : TYPE(ewald_environment_type), POINTER :: ewald_env
333 : TYPE(ewald_pw_type), POINTER :: ewald_pw
334 : INTEGER, INTENT(IN) :: qmmm_coupl_type
335 : TYPE(cell_type), POINTER :: mm_cell
336 : TYPE(mp_para_env_type), POINTER :: para_env
337 : TYPE(section_vals_type), POINTER :: qmmm_periodic, print_section
338 :
339 : INTEGER :: ewald_type, gmax(3), iw, o_spline, ounit
340 : LOGICAL :: do_multipoles
341 : REAL(KIND=dp) :: alpha, rcut
342 : TYPE(cp_logger_type), POINTER :: logger
343 : TYPE(section_vals_type), POINTER :: ewald_print_section, ewald_section, &
344 : poisson_section
345 :
346 8 : logger => cp_get_default_logger()
347 8 : ounit = cp_logger_get_default_io_unit(logger)
348 8 : CPASSERT(.NOT. ASSOCIATED(ewald_env))
349 8 : CPASSERT(.NOT. ASSOCIATED(ewald_pw))
350 :
351 : ! Create Ewald environments
352 8 : poisson_section => section_vals_get_subs_vals(qmmm_periodic, "POISSON")
353 144 : ALLOCATE (ewald_env)
354 8 : CALL ewald_env_create(ewald_env, para_env)
355 8 : CALL ewald_env_set(ewald_env, poisson_section=poisson_section)
356 8 : ewald_section => section_vals_get_subs_vals(poisson_section, "EWALD")
357 8 : CALL read_ewald_section(ewald_env, ewald_section)
358 8 : ewald_print_section => section_vals_get_subs_vals(print_section, "GRID_INFORMATION")
359 8 : ALLOCATE (ewald_pw)
360 8 : CALL ewald_pw_create(ewald_pw, ewald_env, mm_cell, mm_cell, print_section=ewald_print_section)
361 :
362 : CALL ewald_env_get(ewald_env, ewald_type=ewald_type, do_multipoles=do_multipoles, &
363 8 : gmax=gmax, o_spline=o_spline, alpha=alpha, rcut=rcut)
364 8 : IF (do_multipoles) &
365 0 : CPABORT("No multipole force fields allowed in QM-QM Ewald long range correction")
366 :
367 8 : SELECT CASE (qmmm_coupl_type)
368 : CASE (do_qmmm_coulomb)
369 0 : CPABORT("QM-QM long range correction not possible with COULOMB coupling")
370 : CASE (do_qmmm_pcharge)
371 : ! OK
372 : CASE (do_qmmm_gauss, do_qmmm_swave)
373 0 : CPABORT("QM-QM long range correction not possible with GAUSS/SWAVE coupling")
374 : CASE DEFAULT
375 : ! We should never get to this point
376 8 : CPABORT("")
377 : END SELECT
378 :
379 8 : iw = cp_print_key_unit_nr(logger, print_section, "PERIODIC_INFO", extension=".log")
380 8 : IF (iw > 0) THEN
381 2 : WRITE (UNIT=iw, FMT="(/,T2,A)") REPEAT("-", 79)
382 2 : WRITE (UNIT=iw, FMT="(T2,A,T20,A,T80,A)") "-", "QMMM PERIODIC BOUNDARY CONDITION INFO", "-"
383 2 : WRITE (UNIT=iw, FMT="(T2,A,T25,A,T80,A)") "-", "QM-QM Long Range Correction", "-"
384 2 : WRITE (UNIT=iw, FMT="(T2,A)") REPEAT("-", 79)
385 0 : SELECT CASE (ewald_type)
386 : CASE (do_ewald_none)
387 0 : CPABORT("QM-QM long range correction not compatible with Ewald=NONE")
388 : CASE (do_ewald_pme)
389 0 : CPABORT("QM-QM long range correction not possible with Ewald=PME")
390 : CASE (do_ewald_ewald)
391 0 : CPABORT("QM-QM long range correction not possible with Ewald method")
392 : CASE (do_ewald_spme)
393 2 : WRITE (UNIT=iw, FMT="(T2,A,T35,A,T75,A,T80,A)") "-", "Ewald type", "SPME", "-"
394 2 : WRITE (UNIT=iw, FMT="(T2,A,T35,A,T61,3I6,T80,A)") "-", "GMAX values", gmax, "-"
395 2 : WRITE (UNIT=iw, FMT="(T2,A,T35,A,T73,I6,T80,A)") "-", "Spline order", o_spline, "-"
396 2 : WRITE (UNIT=iw, FMT="(T2,A,T35,A,T67,F12.4,T80,A)") "-", "Alpha value", alpha, "-"
397 4 : WRITE (UNIT=iw, FMT="(T2,A,T35,A,T67,F12.4,T80,A)") "-", "Real space cutoff value", rcut, "-"
398 : END SELECT
399 2 : WRITE (UNIT=iw, FMT="(T2,A)") REPEAT("-", 79)
400 : END IF
401 8 : CALL cp_print_key_finished_output(iw, logger, print_section, "PERIODIC_INFO")
402 :
403 8 : END SUBROUTINE qmmm_ewald_potential_init
404 :
405 : END MODULE qmmm_per_elpot
|
