40 #include "./base/base_uses.f90"
44 LOGICAL,
PRIVATE,
PARAMETER :: debug_this_module = .true.
45 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'cp_ddapc_types'
54 REAL(KIND=
dp) :: c0 = 0.0_dp
55 REAL(KIND=
dp),
DIMENSION(:, :),
POINTER :: ami => null()
56 REAL(KIND=
dp),
DIMENSION(:, :),
POINTER :: md => null()
57 REAL(KIND=
dp),
DIMENSION(:, :),
POINTER :: mr => null()
58 REAL(KIND=
dp),
DIMENSION(:, :),
POINTER :: mt => null()
59 REAL(KIND=
dp),
DIMENSION(:, :),
POINTER :: ms => null()
60 REAL(KIND=
dp),
POINTER,
DIMENSION(:, :) :: gfunc => null()
61 REAL(KIND=
dp),
POINTER,
DIMENSION(:) :: w => null()
62 END TYPE cp_ddapc_type
65 TYPE cp_ddapc_ewald_type
66 LOGICAL :: do_decoupling = .false.
67 LOGICAL :: do_qmmm_periodic_decpl = .false.
68 LOGICAL :: do_solvation = .false.
69 LOGICAL :: do_property = .false.
70 LOGICAL :: do_restraint = .false.
71 TYPE(section_vals_type),
POINTER :: ewald_section => null()
72 TYPE(pw_pool_type),
POINTER :: pw_pool_qm => null(), pw_pool_mm => null()
73 TYPE(pw_grid_type),
POINTER :: pw_grid_qm => null(), pw_grid_mm => null()
74 TYPE(pw_r3d_rs_type),
POINTER :: coeff_qm => null(), coeff_mm => null()
75 END TYPE cp_ddapc_ewald_type
97 particle_set, radii, cell, super_cell, rho_tot_g, gcut, iw2, Vol, &
99 TYPE(mp_para_env_type),
POINTER :: cp_para_env
100 TYPE(cp_ddapc_type),
INTENT(OUT) :: cp_ddapc_env
101 TYPE(cp_ddapc_ewald_type),
POINTER :: cp_ddapc_ewald
102 TYPE(particle_type),
DIMENSION(:),
POINTER :: particle_set
103 REAL(kind=
dp),
DIMENSION(:),
POINTER :: radii
104 TYPE(cell_type),
POINTER :: cell, super_cell
105 TYPE(pw_c1d_gs_type),
INTENT(IN) :: rho_tot_g
106 REAL(kind=
dp),
INTENT(IN) :: gcut
107 INTEGER,
INTENT(IN) :: iw2
108 REAL(kind=
dp),
INTENT(IN) :: vol
109 TYPE(section_vals_type),
POINTER :: force_env_section
111 CHARACTER(len=*),
PARAMETER :: routinen =
'cp_ddapc_create'
114 TYPE(section_vals_type),
POINTER :: param_section, solvation_section
116 CALL timeset(routinen, handle)
117 NULLIFY (cp_ddapc_env%AmI, &
122 cp_ddapc_env%gfunc, &
125 CALL ddapc_eval_gfunc(cp_ddapc_env%gfunc, cp_ddapc_env%w, gcut, rho_tot_g, radii)
128 cp_ddapc_env%gfunc, &
136 IF (cp_ddapc_ewald%do_qmmm_periodic_decpl .OR. &
137 cp_ddapc_ewald%do_decoupling)
THEN
141 param_section => cp_ddapc_ewald%ewald_section
150 IF (cp_ddapc_ewald%do_qmmm_periodic_decpl .OR. cp_ddapc_ewald%do_decoupling)
THEN
151 ALLOCATE (cp_ddapc_env%Mt(
SIZE(cp_ddapc_env%Md, 1),
SIZE(cp_ddapc_env%Md, 2)))
152 IF (cp_ddapc_ewald%do_decoupling)
THEN
154 cp_ddapc_env%Mt = cp_ddapc_env%Md
158 CALL ewald_ddapc_pot(cp_para_env, cp_ddapc_ewald%coeff_mm, -1.0_dp, super_cell, param_section, &
159 particle_set, cp_ddapc_env%Mr, radii)
160 cp_ddapc_env%Mt = cp_ddapc_env%Md + cp_ddapc_env%Mr
164 IF (cp_ddapc_ewald%do_solvation)
THEN
168 particle_set, cp_ddapc_env%Ms, radii)
170 CALL timestop(handle)
181 TYPE(cp_ddapc_type),
INTENT(INOUT) :: cp_ddapc_env
183 IF (
ASSOCIATED(cp_ddapc_env%AmI))
THEN
184 DEALLOCATE (cp_ddapc_env%AmI)
186 IF (
ASSOCIATED(cp_ddapc_env%Mt))
THEN
187 DEALLOCATE (cp_ddapc_env%Mt)
189 IF (
ASSOCIATED(cp_ddapc_env%Md))
THEN
190 DEALLOCATE (cp_ddapc_env%Md)
192 IF (
ASSOCIATED(cp_ddapc_env%Mr))
THEN
193 DEALLOCATE (cp_ddapc_env%Mr)
195 IF (
ASSOCIATED(cp_ddapc_env%Ms))
THEN
196 DEALLOCATE (cp_ddapc_env%Ms)
198 IF (
ASSOCIATED(cp_ddapc_env%gfunc))
THEN
199 DEALLOCATE (cp_ddapc_env%gfunc)
201 IF (
ASSOCIATED(cp_ddapc_env%w))
THEN
202 DEALLOCATE (cp_ddapc_env%w)
220 force_env_section, subsys_section, para_env)
221 TYPE(cp_ddapc_ewald_type),
POINTER :: cp_ddapc_ewald
222 LOGICAL,
INTENT(IN) :: qmmm_decoupl
223 TYPE(cell_type),
POINTER :: qm_cell
224 TYPE(section_vals_type),
POINTER :: force_env_section, subsys_section
225 TYPE(mp_para_env_type),
POINTER :: para_env
227 INTEGER :: my_val, npts(3), output_unit
228 INTEGER,
DIMENSION(:),
POINTER :: ngrids
229 LOGICAL :: analyt, decoupling, &
230 do_qmmm_periodic_decpl, do_restraint, &
231 do_restraintb, do_solvation
232 REAL(kind=
dp) :: hmat(3, 3)
233 REAL(kind=
dp),
DIMENSION(:),
POINTER :: gx, gy, gz, lg
234 TYPE(cell_type),
POINTER :: dummy_cell, mm_cell
235 TYPE(cp_logger_type),
POINTER :: logger
236 TYPE(section_vals_type),
POINTER :: cell_section, grid_print_section, multipole_section, &
237 poisson_section, printc_section, qmmm_per_section, restraint_section, restraint_sectionb, &
242 cpassert(.NOT.
ASSOCIATED(cp_ddapc_ewald))
243 ALLOCATE (cp_ddapc_ewald)
244 NULLIFY (cp_ddapc_ewald%pw_grid_mm, &
245 cp_ddapc_ewald%pw_grid_qm, &
246 cp_ddapc_ewald%ewald_section, &
247 cp_ddapc_ewald%pw_pool_mm, &
248 cp_ddapc_ewald%pw_pool_qm, &
249 cp_ddapc_ewald%coeff_mm, &
250 cp_ddapc_ewald%coeff_qm)
251 NULLIFY (multipole_section)
259 "PROPERTIES%ET_COUPLING%DDAPC_RESTRAINT_A")
264 do_qmmm_periodic_decpl = qmmm_decoupl
265 cp_ddapc_ewald%do_solvation = do_solvation
266 cp_ddapc_ewald%do_qmmm_periodic_decpl = do_qmmm_periodic_decpl
268 cp_ddapc_ewald%do_restraint = do_restraint .OR. do_restraintb
270 IF (do_qmmm_periodic_decpl .AND. decoupling)
THEN
272 IF (output_unit > 0) &
273 WRITE (output_unit,
'(T2,"WARNING",A)') &
274 "A calculation with the QMMM periodic model has been requested.", &
275 "The explicit POISSON section in DFT section will be IGNORED.", &
276 "QM Electrostatic controlled only by the PERIODIC section in QMMM section"
289 cp_ddapc_ewald%do_decoupling = decoupling
290 IF (cp_ddapc_ewald%do_qmmm_periodic_decpl)
THEN
294 cp_ddapc_ewald%ewald_section => multipole_section
295 IF (cp_ddapc_ewald%do_decoupling .OR. cp_ddapc_ewald%do_qmmm_periodic_decpl)
THEN
298 IF (.NOT. analyt)
THEN
302 NULLIFY (lg, gx, gy, gz)
304 CALL eval_lg(multipole_section, hmat, qm_cell%deth, lg, gx, gy, gz)
306 CALL setup_ewald_spline(pw_grid=cp_ddapc_ewald%pw_grid_qm, pw_pool=cp_ddapc_ewald%pw_pool_qm, &
307 coeff=cp_ddapc_ewald%coeff_qm, lg=lg, gx=gx, gy=gy, gz=gz, hmat=hmat, npts=npts, &
308 param_section=multipole_section, tag=
"ddapc", &
309 para_env=para_env, print_section=grid_print_section)
314 IF (cp_ddapc_ewald%do_qmmm_periodic_decpl)
THEN
315 NULLIFY (mm_cell, dummy_cell)
317 CALL read_cell(mm_cell, dummy_cell, cell_section=cell_section, para_env=para_env)
319 CALL eval_lg(multipole_section, hmat, mm_cell%deth, lg, gx, gy, gz)
321 CALL setup_ewald_spline(pw_grid=cp_ddapc_ewald%pw_grid_mm, pw_pool=cp_ddapc_ewald%pw_pool_mm, &
322 coeff=cp_ddapc_ewald%coeff_mm, lg=lg, gx=gx, gy=gy, gz=gz, hmat=hmat, npts=npts, &
323 param_section=multipole_section, tag=
"ddapc", para_env=para_env, &
324 print_section=grid_print_section)
349 SUBROUTINE eval_lg(multipole_section, hmat, deth, LG, gx, gy, gz)
350 TYPE(section_vals_type),
POINTER :: multipole_section
351 REAL(kind=
dp),
INTENT(IN) :: hmat(3, 3), deth
352 REAL(kind=
dp),
DIMENSION(:),
POINTER :: lg, gx, gy, gz
354 INTEGER :: i, k1, k2, k3, n_rep, ndim, nmax1, &
356 REAL(kind=
dp) :: alpha, eps,
fac, fs, fvec(3), galpha, &
357 gsq, gsqi, rcut, tol, tol1
359 rcut = min(hmat(1, 1), hmat(2, 2), hmat(3, 3))/2.0_dp
363 eps = min(abs(eps), 0.5_dp)
364 tol = sqrt(abs(log(eps*rcut)))
365 alpha = sqrt(abs(log(eps*rcut*tol)))/rcut
366 galpha = 1.0_dp/(4.0_dp*alpha*alpha)
367 tol1 = sqrt(-log(eps*rcut*(2.0_dp*tol*alpha)**2))
368 nmax1 = nint(0.25_dp + hmat(1, 1)*alpha*tol1/
pi)
369 nmax2 = nint(0.25_dp + hmat(2, 2)*alpha*tol1/
pi)
370 nmax3 = nint(0.25_dp + hmat(3, 3)*alpha*tol1/
pi)
372 fvec = 2.0_dp*
pi/(/hmat(1, 1), hmat(2, 2), hmat(3, 3)/)
373 ndim = (nmax1 + 1)*(2*nmax2 + 1)*(2*nmax3 + 1) - 1
381 DO k2 = -nmax2, nmax2
382 DO k3 = -nmax3, nmax3
383 IF (k1 == 0 .AND. k2 == 0 .AND. k3 == 0) cycle
385 fs = 2.0_dp;
IF (k1 == 0) fs = 1.0_dp
386 gx(i) = fvec(1)*real(k1, kind=
dp)
387 gy(i) = fvec(2)*real(k2, kind=
dp)
388 gz(i) = fvec(3)*real(k3, kind=
dp)
389 gsq = gx(i)*gx(i) + gy(i)*gy(i) + gz(i)*gz(i)
391 lg(i) =
fac*gsqi*exp(-galpha*gsq)
396 END SUBROUTINE eval_lg
406 TYPE(cp_ddapc_ewald_type),
POINTER :: cp_ddapc_ewald
408 IF (
ASSOCIATED(cp_ddapc_ewald))
THEN
409 IF (
ASSOCIATED(cp_ddapc_ewald%coeff_qm))
THEN
410 CALL cp_ddapc_ewald%pw_pool_qm%give_back_pw(cp_ddapc_ewald%coeff_qm)
411 DEALLOCATE (cp_ddapc_ewald%coeff_qm)
413 IF (
ASSOCIATED(cp_ddapc_ewald%coeff_mm))
THEN
414 CALL cp_ddapc_ewald%pw_pool_mm%give_back_pw(cp_ddapc_ewald%coeff_mm)
415 DEALLOCATE (cp_ddapc_ewald%coeff_mm)
417 IF (
ASSOCIATED(cp_ddapc_ewald%pw_pool_qm))
THEN
419 cpassert(.NOT.
ASSOCIATED(cp_ddapc_ewald%pw_pool_qm))
421 IF (
ASSOCIATED(cp_ddapc_ewald%pw_pool_mm))
THEN
423 cpassert(.NOT.
ASSOCIATED(cp_ddapc_ewald%pw_pool_mm))
425 IF (
ASSOCIATED(cp_ddapc_ewald%pw_grid_qm))
THEN
427 cpassert(.NOT.
ASSOCIATED(cp_ddapc_ewald%pw_grid_qm))
429 IF (
ASSOCIATED(cp_ddapc_ewald%pw_grid_mm))
THEN
431 cpassert(.NOT.
ASSOCIATED(cp_ddapc_ewald%pw_grid_mm))
433 DEALLOCATE (cp_ddapc_ewald)
static GRID_HOST_DEVICE double fac(const int i)
Factorial function, e.g. fac(5) = 5! = 120.
Handles all functions related to the CELL.
recursive subroutine, public read_cell(cell, cell_ref, use_ref_cell, cell_section, check_for_ref, para_env)
...
Handles all functions related to the CELL.
subroutine, public cell_release(cell)
releases the given cell (see doc/ReferenceCounting.html)
contains information regarding the decoupling/recoupling method of Bloechl
subroutine, public ddapc_eval_ami(GAmI, c0, gfunc, w, particle_set, gcut, rho_tot_g, radii, iw, Vol)
Computes the inverse AmI of the Am matrix.
recursive subroutine, public ewald_ddapc_pot(cp_para_env, coeff, factor, cell, multipole_section, particle_set, M, radii)
Evaluates the Ewald term E2 and E3 energy term for the decoupling/coupling of periodic images.
subroutine, public solvation_ddapc_pot(solvation_section, particle_set, M, radii)
Evaluates the electrostatic potential due to a simple solvation model Spherical cavity in a dieletric...
subroutine, public ddapc_eval_gfunc(gfunc, w, gcut, rho_tot_g, radii)
...
contains information regarding the decoupling/recoupling method of Bloechl
subroutine, public cp_ddapc_release(cp_ddapc_env)
...
subroutine, public cp_ddapc_create(cp_para_env, cp_ddapc_env, cp_ddapc_ewald, particle_set, radii, cell, super_cell, rho_tot_g, gcut, iw2, Vol, force_env_section)
...
subroutine, public cp_ddapc_ewald_create(cp_ddapc_ewald, qmmm_decoupl, qm_cell, force_env_section, subsys_section, para_env)
...
subroutine, public cp_ddapc_ewald_release(cp_ddapc_ewald)
...
various routines to log and control the output. The idea is that decisions about where to log should ...
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
logical function, public cp_printkey_is_on(iteration_info, print_key)
returns true if the printlevel activates this printkey does not look if this iteration it should be p...
Setting up the Spline coefficients used to Interpolate the G-Term in Ewald sums.
subroutine, public setup_ewald_spline(pw_grid, pw_pool, coeff, LG, gx, gy, gz, hmat, npts, param_section, tag, print_section, para_env)
Setup of the G-space Ewald Term Spline Coefficients.
Defines the basic variable types.
integer, parameter, public dp
Definition of mathematical constants and functions.
real(kind=dp), parameter, public pi
Interface to the message passing library MPI.
Define the data structure for the particle information.
This module defines the grid data type and some basic operations on it.
subroutine, public pw_grid_release(pw_grid)
releases the given pw grid
functions related to the poisson solver on regular grids
integer, parameter, public pw_poisson_multipole
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
subroutine, public pw_pool_release(pool)
releases the given pool (see cp2k/doc/ReferenceCounting.html)