41 #include "./base/base_uses.f90"
46 LOGICAL,
PRIVATE,
PARAMETER :: debug_this_module = .true.
47 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'ewald_spline_util'
72 param_section, tag, print_section, para_env)
73 TYPE(pw_grid_type),
POINTER :: pw_grid
74 TYPE(pw_pool_type),
POINTER :: pw_pool
75 TYPE(pw_r3d_rs_type),
POINTER :: coeff
76 REAL(kind=
dp),
DIMENSION(:),
POINTER :: lg, gx, gy, gz
77 REAL(kind=
dp),
INTENT(IN) :: hmat(3, 3)
78 INTEGER,
INTENT(IN) :: npts(3)
79 TYPE(section_vals_type),
POINTER :: param_section
80 CHARACTER(LEN=*),
INTENT(IN) :: tag
81 TYPE(section_vals_type),
POINTER :: print_section
82 TYPE(mp_para_env_type),
POINTER :: para_env
84 INTEGER :: bo(2, 3), iounit
85 TYPE(cell_type),
POINTER :: cell
86 TYPE(cp_logger_type),
POINTER :: logger
87 TYPE(pw_r3d_rs_type) :: pw
93 cpassert(.NOT.
ASSOCIATED(pw_grid))
94 cpassert(.NOT.
ASSOCIATED(pw_pool))
95 cpassert(.NOT.
ASSOCIATED(coeff))
98 CALL cell_create(cell, hmat=hmat, periodic=(/1, 1, 1/))
104 bo(2, 1:3) = npts(1:3) - 1
112 CALL pw_pool%create_pw(pw)
113 CALL pw_pool%create_pw(coeff)
115 CALL eval_pw_tablr(pw, pw_pool, coeff, lg, gx, gy, gz, hmat_mm=hmat, &
116 param_section=param_section, tag=tag)
117 CALL pw_pool%give_back_pw(pw)
139 SUBROUTINE eval_pw_tablr(grid, pw_pool, TabLR, Lg, gx, gy, gz, hmat_mm, &
141 TYPE(pw_r3d_rs_type),
INTENT(INOUT) :: grid
142 TYPE(pw_pool_type),
POINTER :: pw_pool
143 TYPE(pw_r3d_rs_type),
POINTER :: tablr
144 REAL(kind=
dp),
DIMENSION(:),
POINTER :: lg, gx, gy, gz
145 REAL(kind=
dp),
DIMENSION(3, 3) :: hmat_mm
146 TYPE(section_vals_type),
POINTER :: param_section
147 CHARACTER(LEN=*),
INTENT(IN) :: tag
149 CHARACTER(len=*),
PARAMETER :: routinen =
'eval_pw_TabLR'
151 INTEGER :: act_nx, act_ny, aint_precond, handle, i, iii, is, j, js, k, ks, lg_loc_max, &
152 lg_loc_min, max_iter, my_i, my_j, my_k, n1, n2, n3, n_extra, nlg_loc, nxlim, nylim, &
154 INTEGER,
DIMENSION(2, 3) :: gbo
156 REAL(kind=
dp) :: dr1, dr2, dr3, eps_r, eps_x, xs1, xs2, &
158 REAL(kind=
dp),
ALLOCATABLE :: cos_gx(:, :), cos_gy(:, :), &
159 cos_gz(:, :), lhs(:, :), rhs(:, :), &
160 sin_gx(:, :), sin_gy(:, :), &
162 TYPE(pw_spline_precond_type) :: precond
163 TYPE(section_vals_type),
POINTER :: interp_section
170 CALL timeset(routinen, handle)
171 n1 = grid%pw_grid%npts(1)
172 n2 = grid%pw_grid%npts(2)
173 n3 = grid%pw_grid%npts(3)
174 dr1 = grid%pw_grid%dr(1)
175 dr2 = grid%pw_grid%dr(2)
176 dr3 = grid%pw_grid%dr(3)
177 gbo = grid%pw_grid%bounds
178 nxlim = floor(real(n1, kind=
dp)/2.0_dp)
179 nylim = floor(real(n2, kind=
dp)/2.0_dp)
180 nzlim = floor(real(n3, kind=
dp)/2.0_dp)
184 IF (2*nxlim /= n1) is = 1
185 IF (2*nylim /= n2) js = 1
186 IF (2*nzlim /= n3) ks = 1
193 ALLOCATE (cos_gx(
SIZE(lg), gbo(1, 1):gbo(2, 1)))
194 ALLOCATE (sin_gx(
SIZE(lg), gbo(1, 1):gbo(2, 1)))
195 ALLOCATE (cos_gy(
SIZE(lg), gbo(1, 2):gbo(2, 2)))
196 ALLOCATE (sin_gy(
SIZE(lg), gbo(1, 2):gbo(2, 2)))
197 ALLOCATE (cos_gz(
SIZE(lg), gbo(1, 3):gbo(2, 3)))
198 ALLOCATE (sin_gz(
SIZE(lg), gbo(1, 3):gbo(2, 3)))
200 DO k = gbo(1, 3), gbo(2, 3)
202 xs3 = real(my_k,
dp)*dr3
204 cos_gz(1:
SIZE(lg), k) = cos(gz(1:
SIZE(lg))*xs3)
205 sin_gz(1:
SIZE(lg), k) = sin(gz(1:
SIZE(lg))*xs3)
208 DO j = gbo(1, 2), gbo(2, 2)
210 cos_gy(1:
SIZE(lg), j) = cos(gy(1:
SIZE(lg))*xs2)
211 sin_gy(1:
SIZE(lg), j) = sin(gy(1:
SIZE(lg))*xs2)
215 DO i = gbo(1, 1), gbo(2, 1)
217 cos_gx(1:
SIZE(lg), i) = cos(gx(1:
SIZE(lg))*xs1)
218 sin_gx(1:
SIZE(lg), i) = sin(gx(1:
SIZE(lg))*xs1)
224 nlg_loc =
SIZE(lg)/grid%pw_grid%para%group_size
226 n_extra = mod(
SIZE(lg), grid%pw_grid%para%group_size)
228 IF (grid%pw_grid%para%my_pos < n_extra)
THEN
229 lg_loc_min = (nlg_loc + 1)*grid%pw_grid%para%my_pos + 1
230 lg_loc_max = lg_loc_min + (nlg_loc + 1) - 1
232 lg_loc_min = (nlg_loc + 1)*n_extra + nlg_loc*(grid%pw_grid%para%my_pos - n_extra) + 1
233 lg_loc_max = lg_loc_min + nlg_loc - 1
236 lg_loc_max = min(
SIZE(lg), lg_loc_max)
237 nlg_loc = lg_loc_max - lg_loc_min + 1
239 IF (nlg_loc > 0)
THEN
241 act_nx = min(gbo(2, 1), nxlim) - gbo(1, 1) + 1
242 act_ny = min(gbo(2, 2), nylim) - gbo(1, 2) + 1
244 ALLOCATE (lhs(act_nx, nlg_loc))
245 ALLOCATE (rhs(act_ny, nlg_loc))
248 DO i = gbo(1, 1), gbo(2, 1)
250 lhs(i - gbo(1, 1) + 1, 1:nlg_loc) = lg(lg_loc_min:lg_loc_max)*cos_gx(lg_loc_min:lg_loc_max, i)
252 DO k = gbo(1, 3), gbo(2, 3)
254 DO j = gbo(1, 2), gbo(2, 2)
256 rhs(j - gbo(1, 2) + 1, 1:nlg_loc) = cos_gy(lg_loc_min:lg_loc_max, j)*cos_gz(lg_loc_min:lg_loc_max, k) - &
257 sin_gy(lg_loc_min:lg_loc_max, j)*sin_gz(lg_loc_min:lg_loc_max, k)
259 CALL dgemm(
'N',
'T', act_nx, act_ny, nlg_loc, 1.0d0, lhs(1, 1), act_nx, rhs(1, 1), act_ny, 0.0d0, &
260 grid%array(gbo(1, 1), gbo(1, 2), k),
SIZE(grid%array, 1))
264 DO i = gbo(1, 1), gbo(2, 1)
266 lhs(i - gbo(1, 1) + 1, 1:nlg_loc) = -lg(lg_loc_min:lg_loc_max)*sin_gx(lg_loc_min:lg_loc_max, i)
268 DO k = gbo(1, 3), gbo(2, 3)
270 DO j = gbo(1, 2), gbo(2, 2)
272 rhs(j - gbo(1, 2) + 1, 1:nlg_loc) = cos_gy(lg_loc_min:lg_loc_max, j)*sin_gz(lg_loc_min:lg_loc_max, k) + &
273 sin_gy(lg_loc_min:lg_loc_max, j)*cos_gz(lg_loc_min:lg_loc_max, k)
275 CALL dgemm(
'N',
'T', act_nx, act_ny, nlg_loc, 1.0d0, lhs(1, 1), act_nx, rhs(1, 1), act_ny, 1.0d0, &
276 grid%array(gbo(1, 1), gbo(1, 2), k),
SIZE(grid%array, 1))
291 grid%array(gbo(1, 1):nxlim, gbo(1, 2):nylim, gbo(1, 3):nzlim) = 0.0_dp
294 CALL grid%pw_grid%para%group%sum(grid%array(gbo(1, 1):nxlim, gbo(1, 2):nylim, gbo(1, 3):nzlim))
296 fake_loopongrid:
DO k = gbo(1, 3), gbo(2, 3)
298 IF (k > nzlim) my_k = nzlim - abs(nzlim - k) + ks
299 DO j = gbo(1, 2), gbo(2, 2)
301 IF (j > nylim) my_j = nylim - abs(nylim - j) + js
302 DO i = gbo(1, 1), gbo(2, 1)
304 IF (i > nxlim) my_i = nxlim - abs(nxlim - i) + is
305 grid%array(i, j, k) = grid%array(my_i, my_j, my_k)
308 END DO fake_loopongrid
322 pool=pw_pool,
pbc=.true., transpose=.false.)
327 eps_r=eps_r, eps_x=eps_x, &
334 CALL check_spline_interp_tablr(hmat_mm, lg, gx, gy, gz, tablr, param_section, &
336 CALL timestop(handle)
337 END SUBROUTINE eval_pw_tablr
353 SUBROUTINE check_spline_interp_tablr(hmat_mm, Lg, gx, gy, gz, TabLR, &
355 REAL(kind=
dp),
DIMENSION(3, 3) :: hmat_mm
356 REAL(kind=
dp),
DIMENSION(:),
POINTER :: lg, gx, gy, gz
357 TYPE(pw_r3d_rs_type),
POINTER :: tablr
358 TYPE(section_vals_type),
POINTER :: param_section
359 CHARACTER(LEN=*),
INTENT(IN) :: tag
361 CHARACTER(len=*),
PARAMETER :: routinen =
'check_spline_interp_TabLR'
363 INTEGER :: handle, i, iw, kg, npoints
364 REAL(kind=
dp) :: dn(3), dr1, dr2, dr3, dxterm, dyterm, &
365 dzterm, errd, errf, fterm, maxerrord, &
366 maxerrorf, na, nn, term, tmp1, tmp2, &
367 vec(3), xs1, xs2, xs3
368 TYPE(cp_logger_type),
POINTER :: logger
373 extension=
"."//trim(tag)//
"Log")
374 CALL timeset(routinen, handle)
381 dr1 = hmat_mm(1, 1)/real(npoints, kind=
dp)
382 dr2 = hmat_mm(2, 2)/real(npoints, kind=
dp)
383 dr3 = hmat_mm(3, 3)/real(npoints, kind=
dp)
387 WRITE (iw,
'(A,T5,A15,4X,A17,T50,4X,A,5X,A,T80,A,T85,A15,4X,A17,T130,4X,A,5X,A)') &
388 "#",
"Analytical Term",
"Interpolated Term",
"Error",
"MaxError", &
389 "*",
" Analyt Deriv ",
"Interp Deriv Mod ",
"Error",
"MaxError"
390 DO i = 1, npoints + 1
397 vec = (/real(gx(kg), kind=
dp), real(gy(kg), kind=
dp), real(gz(kg), kind=
dp)/)
398 term = term + lg(kg)*cos(vec(1)*xs1 + vec(2)*xs2 + vec(3)*xs3)
399 dxterm = dxterm - lg(kg)*sin(vec(1)*xs1 + vec(2)*xs2 + vec(3)*xs3)*vec(1)
400 dyterm = dyterm - lg(kg)*sin(vec(1)*xs1 + vec(2)*xs2 + vec(3)*xs3)*vec(2)
401 dzterm = dzterm - lg(kg)*sin(vec(1)*xs1 + vec(2)*xs2 + vec(3)*xs3)*vec(3)
403 na = sqrt(dxterm*dxterm + dyterm*dyterm + dzterm*dzterm)
405 nn = sqrt(dot_product(dn, dn))
407 tmp1 = abs(term - fterm)
408 tmp2 = sqrt(dot_product(dn - (/dxterm, dyterm, dzterm/), dn - (/dxterm, dyterm, dzterm/)))
410 maxerrorf = max(maxerrorf, tmp1)
412 maxerrord = max(maxerrord, tmp2)
413 WRITE (iw,
'(T5,F15.10,5X,F15.10,T50,2F12.9,T80,A,T85,F15.10,5X,F15.10,T130,2F12.9)') &
414 term, fterm, tmp1, maxerrorf,
"*", na, nn, tmp2, maxerrord
419 WRITE (iw,
'(A,T5,A,T50,F12.9,T130,F12.9)')
"#",
"Averages", errf/real(npoints, kind=
dp), &
420 errd/real(npoints, kind=
dp)
422 CALL timestop(handle)
425 END SUBROUTINE check_spline_interp_tablr
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Handles all functions related to the CELL.
subroutine, public cell_create(cell, hmat, periodic, tag)
allocates and initializes a cell
Handles all functions related to the CELL.
subroutine, public cell_release(cell)
releases the given cell (see doc/ReferenceCounting.html)
various routines to log and control the output. The idea is that decisions about where to log should ...
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...
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
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
Interface to the message passing library MPI.
computes preconditioners, and implements methods to apply them currently used in qs_ot
integer, parameter, public halfspace
This module defines the grid data type and some basic operations on it.
subroutine, public pw_grid_setup(cell_hmat, pw_grid, grid_span, cutoff, bounds, bounds_local, npts, spherical, odd, fft_usage, ncommensurate, icommensurate, blocked, ref_grid, rs_dims, iounit)
sets up a pw_grid
subroutine, public pw_grid_create(pw_grid, pe_group, local)
Initialize a PW grid with all defaults.
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
subroutine, public pw_pool_create(pool, pw_grid, max_cache)
creates a pool for pw
different utils that are useful to manipulate splines on the regular grid of a pw
logical function, public find_coeffs(values, coeffs, linOp, preconditioner, pool, eps_r, eps_x, max_iter, sumtype)
solves iteratively (CG) a systmes of linear equations linOp(coeffs)=values (for example those needed ...
subroutine, public pw_spline_precond_release(preconditioner)
releases the preconditioner
subroutine, public pw_spline_precond_create(preconditioner, precond_kind, pool, pbc, transpose)
...
subroutine, public pw_spline_do_precond(preconditioner, in_v, out_v)
applies the preconditioner to the system of equations to find the coefficients of the spline
subroutine, public pw_spline_precond_set_kind(preconditioner, precond_kind, pbc, transpose)
switches the types of precoditioner to use
real(kind=dp) function, dimension(3), public eval_d_interp_spl3_pbc(vec, pw)
Evaluates the derivatives of the PBC interpolated Spline (pw) function on the generic input vector (v...
real(kind=dp) function, public eval_interp_spl3_pbc(vec, pw)
Evaluates the PBC interpolated Spline (pw) function on the generic input vector (vec)
subroutine, public spl3_pbc(pw_in, pw_out)
...