|
CP2K 2.4 (Revision 12889)
|
00001 !-----------------------------------------------------------------------------! 00002 ! CP2K: A general program to perform molecular dynamics simulations ! 00003 ! Copyright (C) 2000 - 2013 CP2K developers group ! 00004 !-----------------------------------------------------------------------------! 00005 00006 ! ***************************************************************************** 00014 MODULE gauss_colloc 00015 USE d3_poly 00016 USE f77_blas 00017 USE kinds, ONLY: dp,& 00018 int_8 00019 USE lgrid_types, ONLY: lgrid_type 00020 00021 !$ USE OMP_LIB 00022 #include "cp_common_uses.h" 00023 00024 IMPLICIT NONE 00025 PRIVATE 00026 00027 PUBLIC :: collocGauss, collocGauss_safe,sphRad2,calc_max_r2,& 00028 calcBox , integrateGauss, calcRadius2, integrateGaussFull,& 00029 collocGaussFlat,integrateGaussFullFlat 00030 00031 #ifdef FD_DEBUG 00032 #define IF_CHECK(x,y) x 00033 #else 00034 #define IF_CHECK(x,y) y 00035 #endif 00036 00037 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gauss_colloc' 00038 00039 REAL(dp), PARAMETER :: small=TINY(1.0_dp) 00040 00041 ! keep prettify happy (does not see the include) 00042 INTEGER(KIND=int_8), PARAMETER :: unused_import_of_int_8=1 00043 00044 00045 CONTAINS 00046 00047 !**** radius calculation 00048 00049 ! ***************************************************************************** 00055 FUNCTION calc_max_r2(poly,poly_shift,alpha,epsilon,error,scale) RESULT(res) 00056 REAL(dp), DIMENSION(0:), INTENT(in) :: poly 00057 REAL(dp), DIMENSION(3) :: poly_shift 00058 REAL(dp), INTENT(in) :: alpha, epsilon 00059 TYPE(cp_error_type), INTENT(inout) :: error 00060 REAL(dp), INTENT(in), OPTIONAL :: scale 00061 REAL(dp) :: res 00062 00063 CHARACTER(len=*), PARAMETER :: routineN = 'calc_max_r2', 00064 routineP = moduleN//':'//routineN 00065 REAL(dp), DIMENSION(3, 3), PARAMETER :: 00066 identity_m = RESHAPE((/ 1,0,0, 0,1,0, 0,0,1 /),(/ 3,3 /)) 00067 00068 INTEGER :: grad, stat 00069 REAL(dp) :: eps 00070 REAL(dp), ALLOCATABLE, DIMENSION(:) :: poly_ijk 00071 00072 eps=epsilon 00073 IF (PRESENT(scale)) THEN 00074 IF (scale==0.0_dp) THEN 00075 res=0.0_dp 00076 RETURN 00077 END IF 00078 eps=eps/scale 00079 END IF 00080 grad=grad_size3(SIZE(poly)) 00081 ALLOCATE(poly_ijk(0:poly_size3(grad)-1),stat=stat) 00082 CPPostconditionNoFail(stat==0,cp_fatal_level,routineP,error) 00083 CALL poly_affine_t3(poly,identity_m,poly_shift,poly_ijk,error=error) 00084 res=calcRadius2(poly_ijk,alpha,eps) 00085 DEALLOCATE(poly_ijk,stat=stat) 00086 CPPostconditionNoFail(stat==0,cp_fatal_level,routineP,error) 00087 END FUNCTION 00088 00089 ! ***************************************************************************** 00107 FUNCTION calcRadius2(poly,alpha,epsilon) RESULT(res) 00108 REAL(dp), DIMENSION(0:), INTENT(in) :: poly 00109 REAL(dp), INTENT(in) :: alpha, epsilon 00110 REAL(dp) :: res 00111 00112 INTEGER :: from_i, grad, gradSize, i, 00113 igrad, p_size, to_i 00114 REAL(dp) :: coeffAtt, rAtt, residual, t 00115 REAL(dp), ALLOCATABLE, DIMENSION(:) :: rAtts 00116 00117 res=0.0_dp 00118 grad=grad_size3(SIZE(poly)) 00119 ALLOCATE(rAtts(0:grad)) 00120 from_i=0 00121 gradSize=1 00122 p_size=SIZE(poly) 00123 DO igrad=0,grad 00124 to_i=MIN(from_i+gradSize,p_size) 00125 !! igrad-norm 00126 IF (igrad==0) THEN 00127 coeffAtt=ABS(poly(from_i)) 00128 ELSE 00129 coeffAtt=0.0_dp 00130 DO i=from_i,to_i-1 00131 coeffAtt=coeffAtt+ABS(poly(i))**igrad 00132 END DO 00133 coeffAtt=coeffAtt**(1.0_dp/REAL(igrad)) 00134 END IF 00135 !! infinity-norm 00136 !coeffAtt=maxval(abs(poly(from_i:to_i-1))) 00137 !! 2-norm (euclidean) 00138 !coeffAtt=sqrt(sum(poly(from_i:to_i-1)**2)) 00139 rAtt=sphRad2(coeffAtt,igrad,alpha,epsilon) ! bounds value 00140 !rAtt=sphRad2(4*pi*coeffAtt,igrad+2,alpha,epsilon) ! bounds surface integral 00141 rAtts(igrad)=rAtt 00142 res=MAX(res,rAtts(igrad)) 00143 gradSize=gradSize+(igrad+2) 00144 from_i=to_i 00145 END DO 00146 residual=0.0_dp 00147 DO igrad=0,grad 00148 ! pade expansion 1,2 of exp(-x) at x=0 00149 t=alpha*(res-rAtts(igrad)) 00150 residual=residual+1.0_dp/(1.0_dp+t+0.5_dp*t**2) 00151 END DO 00152 res=res+LOG(residual)/alpha 00153 DEALLOCATE(rAtts) 00154 END FUNCTION 00155 00156 ! ***************************************************************************** 00159 FUNCTION sphRad2(coeff,l,alpha,epsilon) RESULT(res) 00160 REAL(dp), INTENT(in) :: coeff 00161 INTEGER, INTENT(in) :: l 00162 REAL(dp), INTENT(in) :: alpha, epsilon 00163 REAL(dp) :: res 00164 00165 REAL(dp) :: r1, t, x, xx 00166 00167 IF (ABS(coeff)<small) THEN 00168 res=0.0_dp 00169 RETURN 00170 END IF 00171 xx=MAX(0.0_dp,-LOG(epsilon/coeff)/alpha) 00172 x=MAX(1.0_dp,xx) 00173 t=0.5_dp*REAL(l,dp)/alpha 00174 00175 ! rational functions that approximate the solution of r=x+t*log(r) 00176 r1=MAX(x,(5.8679922317012645_dp - 13.433907351015852_dp*t + t**2 + 2.7189706320638596_dp*x)& 00177 /(2.1224230654286007_dp + 1.5855089862522989_dp*t + 0.05338142797029377_dp*t**2 & 00178 + 1.4761929491905326_dp*x - 0.059575647341275406_dp*t*x + 0.031067088387632627_dp*x**2)& 00179 +(-0.43804767398212585_dp - 0.49900615723405245_dp*t + 0.20629903250691906_dp*t**2& 00180 + 0.9017510689333377_dp*x - 0.07046283868925775_dp*t*x - 0.1743645447869195_dp*x**2)& 00181 /(-3.617044072726269_dp - 0.7662306349036399_dp*t - 0.7614179211580372_dp*t**2& 00182 + 0.17934761179419514_dp*t**3 + 6.788069591421069_dp*x + 0.5260819714537482_dp*t*x& 00183 + 1.4780624656229489_dp*t**2*x - 3.05963073688407_dp*x**2 - 1.5521580811226825_dp*t*x**2& 00184 + x**3) - 0.6840379938652372_dp*MAX(0.0_dp,0.4_dp - t + (1.0_dp - x)/3.0_dp)& 00185 - 0.16927886252742727_dp*MAX(0.0_dp,1.0_dp - t + (1.0_dp - x)/3.0_dp) & 00186 +(-2.1054595688424427_dp + 3.0262780685014503_dp*t + 0.1353476536688748_dp*t**2 + x& 00187 + 0.03244405654399495_dp*t*x + 0.0017806182354747218_dp*x**2)& 00188 /( 0.9552112747891883_dp + 0.016045994433298932_dp*t + 0.0018347038147891919_dp*x& 00189 - 1.8296261802177613e-6_dp*MIN(1000.0_dp,t)**2)) 00190 res=MAX(0.0_dp,xx+t*LOG(r1)) ! one step of the iterative solver, relative error less than 0.004 00191 END FUNCTION 00192 00193 ! ***************************************************************************** 00198 FUNCTION calcBox(h,h_inv,posi,max_r2,& 00199 periodic,gdim,error,guarantee_nearest) RESULT(res) 00200 REAL(dp), DIMENSION(0:2, 0:2), 00201 INTENT(in) :: h, h_inv 00202 REAL(dp), DIMENSION(0:2), INTENT(in) :: posi 00203 REAL(dp), INTENT(in) :: max_r2 00204 INTEGER, DIMENSION(0:2), INTENT(in) :: periodic, gdim 00205 TYPE(cp_error_type), INTENT(inout) :: error 00206 LOGICAL, INTENT(in), OPTIONAL :: guarantee_nearest 00207 INTEGER :: res(2,0:2) 00208 00209 CHARACTER(len=*), PARAMETER :: routineN = 'calcBox', 00210 routineP = moduleN//':'//routineN 00211 00212 INTEGER :: i, imax, imin, j, jmax, jmin, 00213 k, kmax, kmin 00214 INTEGER, DIMENSION(0:2) :: cellShift, fullShift, ndim, 00215 shiftPos 00216 LOGICAL :: g_nearest 00217 REAL(dp) :: cci0, cci1, cci2, ccj0, ccj1, ccj2, cck0, cck1, cck2, 00218 delta_i, delta_j, delta_k, m(0:2,0:2), maxr2, r_0, scaled_h(0:2,0:2), 00219 sqDi, sqDj, sqDk 00220 REAL(dp), DIMENSION(0:2) :: l, normD, resPos, riPos, 00221 rpos, wrPos 00222 00223 g_nearest=.TRUE. 00224 IF (PRESENT(guarantee_nearest)) g_nearest=guarantee_nearest 00225 ndim=gdim 00226 rPos=0.0_dp 00227 DO j=0,2 00228 DO i=0,2 00229 rPos(i)=rPos(i)+h_inv(i,j)*posi(j) 00230 END DO 00231 END DO 00232 cellShift=FLOOR(rPos)*periodic 00233 wrPos=rPos-REAL(cellShift,dp) 00234 riPos=wrPos*ndim 00235 shiftPos=FLOOR(riPos+0.5_dp) 00236 fullShift=shiftPos+ndim*cellShift 00237 resPos=riPos-shiftPos 00238 normD=1.0_dp/REAL(ndim,dp) 00239 00240 IF (g_nearest) THEN 00241 maxr2=0.0_dp 00242 DO j=0,2 00243 DO i=0,2 00244 maxr2=maxr2+(h(i,j)*normD(j))**2 00245 END DO 00246 END DO 00247 maxr2=MAX(max_r2,maxr2) 00248 ELSE 00249 maxr2=max_r2 00250 END IF 00251 00252 scaled_h=0.0_dp 00253 DO j=0,2 00254 DO i=0,2 00255 scaled_h(i,2-j)=scaled_h(i,2-j)+h(i,j)*normD(j) 00256 END DO 00257 END DO 00258 00259 ! build up quadratic form (ellipsoid) 00260 m=0.0_dp 00261 DO j=0,2 00262 DO i=0,2 00263 DO k=0,2 00264 m(i,j)=m(i,j)+scaled_h(k,i)*scaled_h(k,j) 00265 END DO 00266 END DO 00267 END DO 00268 00269 l=0.0_dp 00270 DO j=0,2 00271 DO i=0,2 00272 l(j)=l(j)-2.0*resPos(2-i)*m(i,j) 00273 END DO 00274 END DO 00275 00276 r_0=0.0_dp 00277 DO i=0,2 00278 r_0=r_0-0.5*resPos(2-i)*l(i) 00279 END DO 00280 00281 ! i box bounds 00282 cci2 = (m(2,2) * m(0,0) * m(1,1) - m(2,2) * m(0,1) ** 2 - m(1,1) * m(0,2) ** 2 & 00283 + 2.0_dp * m(0,2) * m(0,1) * m(1,2) - m(0,0) * m(1,2) ** 2) / (m(2,2) * m(1,1) - m(1,2) ** 2) 00284 cci1 = -(-m(2,2) * l(0) * m(1,1) + m(2,2) * m(0,1) * l(1) + l(2) * m(0,2) * m(1,1) & 00285 + l(0) * m(1,2) ** 2 - l(2) * m(0,1) * m(1,2) - m(1,2) * l(1) * m(0,2)) / (m(2,2) * m(1,1) - m(1,2) ** 2) 00286 cci0 = -((-4.0_dp * m(2,2) * r_0 * m(1,1) + m(2,2) * l(1) ** 2 + l(2) ** 2 * m(1,1) & 00287 - 2.0_dp * l(1) * m(1,2) * l(2) + 4.0_dp * r_0 * m(1,2) ** 2) & 00288 / (m(2,2) * m(1,1) - m(1,2) ** 2)) / 4.0_dp-maxr2 00289 delta_i=cci1*cci1-4.0_dp*cci2*cci0 00290 IF (delta_i<0.0_dp) THEN 00291 imin=fullShift(2)+CEILING((-cci1)/(2.0_dp*cci2)) 00292 imax=imin-1 00293 ELSE 00294 sqDi=SQRT(delta_i) 00295 imin=fullShift(2)+CEILING((-cci1-sqDi)/(2.0_dp*cci2)) 00296 imax=fullShift(2)+FLOOR((-cci1+sqDi)/(2.0_dp*cci2)) 00297 END IF 00298 IF (periodic(2)==0) THEN 00299 imin=MAX(0,imin) 00300 imax=MIN(imax,ndim(2)-1) 00301 END IF 00302 res(1,2)=imin 00303 res(2,2)=imax 00304 00305 ! j box bounds 00306 ccj2 = (m(0,0) * m(2,2) * m(1,1) - m(0,0) * m(1,2) ** 2 - m(0,1) ** 2 * m(2,2) & 00307 + 2.0_dp * m(0,1) * m(0,2) * m(1,2) - m(1,1) * m(0,2) ** 2) & 00308 / (m(0,0) * m(2,2) - m(0,2) ** 2) 00309 ccj1 = -(-m(0,0) * l(1) * m(2,2) + m(0,0) * l(2) * m(1,2) + l(0) * m(0,1) * m(2,2) & 00310 - m(0,2) * l(2) * m(0,1) - l(0) * m(0,2) * m(1,2) + l(1) * m(0,2) ** 2) & 00311 / (m(0,0) * m(2,2) - m(0,2) ** 2) 00312 ccj0 = (4.0_dp * m(0,0) * m(2,2) * r_0 - m(0,0) * l(2) ** 2 - m(2,2) * l(0) ** 2 & 00313 + 2.0_dp * m(0,2) * l(2) * l(0) - 4.0_dp * r_0 * m(0,2) ** 2) & 00314 / (m(0,0) * m(2,2) - m(0,2) ** 2) / 4.0_dp-maxr2 00315 delta_j=ccj1*ccj1-4.0_dp*ccj2*ccj0 00316 IF (delta_j<0.0_dp) THEN 00317 jmin=fullShift(1)+CEILING((-ccj1)/(2.0_dp*ccj2)) 00318 jmax=jmin-1 00319 ELSE 00320 sqDj=SQRT(delta_j) 00321 jmin=fullShift(1)+CEILING((-ccj1-sqDj)/(2.0_dp*ccj2)) 00322 jmax=fullShift(1)+FLOOR((-ccj1+sqDj)/(2.0_dp*ccj2)) 00323 END IF 00324 IF (periodic(1)==0) THEN 00325 jmin=MAX(0,jmin) 00326 jmax=MIN(jmax,ndim(1)-1) 00327 END IF 00328 res(1,1)=jmin 00329 res(2,1)=jmax 00330 00331 ! k box bounds 00332 cck2 = (m(0,0) * m(2,2) * m(1,1) - m(0,0) * m(1,2) ** 2 - m(0,1) ** 2 * m(2,2) & 00333 + 2.0_dp * m(0,1) * m(0,2) * m(1,2) - m(1,1) * m(0,2) ** 2) / (m(0,0) * m(1,1) - m(0,1) ** 2) 00334 cck1 = (m(0,0) * l(2) * m(1,1) - m(0,0) * m(1,2) * l(1) - m(0,2) * l(0) * m(1,1) & 00335 + l(0) * m(0,1) * m(1,2) - l(2) * m(0,1) ** 2 + m(0,1) * l(1) * m(0,2)) / (m(0,0) * m(1,1) - m(0,1) ** 2) 00336 cck0 = (4.0_dp * m(0,0) * m(1,1) * r_0 - m(0,0) * l(1) ** 2 - m(1,1) * l(0) ** 2 & 00337 + 2.0_dp * m(0,1) * l(1) * l(0) - 4.0_dp * r_0 * m(0,1) ** 2) & 00338 / (m(0,0) * m(1,1) - m(0,1) ** 2) / 4.0_dp-maxr2 00339 delta_k=cck1*cck1-4.0_dp*cck2*cck0 00340 IF (delta_k<0.0_dp) THEN 00341 kmin=fullShift(0)+CEILING((-cck1)/(2.0_dp*cck2)) 00342 kmax=kmin-1 00343 ELSE 00344 sqDk=SQRT(delta_k) 00345 kmin=fullShift(0)+CEILING((-cck1-sqDk)/(2.0_dp*cck2)) 00346 kmax=fullShift(0)+FLOOR((-cck1+sqDk)/(2.0_dp*cck2)) 00347 END IF 00348 IF (periodic(0)==0) THEN 00349 kmin=MAX(0,kmin) 00350 kmax=MIN(kmax,ndim(0)-1) 00351 END IF 00352 res(1,0)=kmin 00353 res(2,0)=kmax 00354 00355 END FUNCTION 00356 00357 ! ***************************************************************************** 00388 SUBROUTINE collocGauss(h,h_inv,grid,poly,alphai,posi,max_r2,& 00389 periodic,gdim,local_bounds,local_shift,poly_shift,scale,lgrid,error) 00390 REAL(dp), DIMENSION(0:2, 0:2), 00391 INTENT(in) :: h, h_inv 00392 REAL(dp), DIMENSION(0:, 0:, 0:), 00393 INTENT(inout) :: grid 00394 REAL(dp), DIMENSION(:), INTENT(in) :: poly 00395 REAL(dp), INTENT(in) :: alphai 00396 REAL(dp), DIMENSION(0:2), INTENT(in) :: posi 00397 REAL(dp), INTENT(in) :: max_r2 00398 INTEGER, DIMENSION(0:2), INTENT(in) :: periodic 00399 INTEGER, DIMENSION(0:2), INTENT(in), 00400 OPTIONAL :: gdim 00401 INTEGER, DIMENSION(2, 0:2), INTENT(in), 00402 OPTIONAL :: local_bounds 00403 INTEGER, DIMENSION(0:2), INTENT(in), 00404 OPTIONAL :: local_shift 00405 REAL(dp), DIMENSION(0:2), INTENT(in), 00406 OPTIONAL :: poly_shift 00407 REAL(dp), INTENT(in), OPTIONAL :: scale 00408 TYPE(lgrid_type), INTENT(inout), 00409 OPTIONAL :: lgrid 00410 TYPE(cp_error_type), INTENT(inout) :: error 00411 00412 CHARACTER(len=*), PARAMETER :: routineN = 'collocGauss', 00413 routineP = moduleN//':'//routineN 00414 00415 00416 #include "colloc_int_body.f90" 00417 00418 END SUBROUTINE 00419 00420 ! ***************************************************************************** 00426 SUBROUTINE collocGaussFlat(h,h_inv,grid,ngpts,ldim,poly,alphai,posi,max_r2,& 00427 periodic,gdim,local_bounds,local_shift,poly_shift,scale,error) 00428 REAL(dp), DIMENSION(0:2, 0:2), 00429 INTENT(in) :: h, h_inv 00430 INTEGER, INTENT(in) :: ngpts,ldim(0:2) 00431 REAL(dp), DIMENSION(IF_CHECK(ngpts,*)), 00432 INTENT(inout) :: grid 00433 REAL(dp), DIMENSION(:), INTENT(in) :: poly 00434 REAL(dp), INTENT(in) :: alphai 00435 REAL(dp), DIMENSION(0:2), INTENT(in) :: posi 00436 REAL(dp), INTENT(in) :: max_r2 00437 INTEGER, DIMENSION(0:2), INTENT(in) :: periodic 00438 INTEGER, DIMENSION(0:2), INTENT(in), 00439 OPTIONAL :: gdim 00440 INTEGER, DIMENSION(2, 0:2), INTENT(in), 00441 OPTIONAL :: local_bounds 00442 INTEGER, DIMENSION(0:2), INTENT(in), 00443 OPTIONAL :: local_shift 00444 REAL(dp), DIMENSION(0:2), INTENT(in), 00445 OPTIONAL :: poly_shift 00446 REAL(dp), INTENT(in), OPTIONAL :: scale 00447 TYPE(cp_error_type), INTENT(inout) :: error 00448 00449 CHARACTER(len=*), PARAMETER :: routineN = 'collocGaussFlat', 00450 routineP=moduleN//':'//routineN 00451 #define FM_FLAT_GRID 00452 #include "colloc_int_body.f90" 00453 #undef FM_FLAT_GRID 00454 END SUBROUTINE 00455 00456 ! ***************************************************************************** 00466 SUBROUTINE integrateGauss(h,h_inv,grid,poly,alphai,posi,max_r2,& 00467 periodic,npoly,res,gdim,local_bounds,local_shift,poly_shift,scale,error) 00468 INTEGER, INTENT(in) :: npoly 00469 REAL(dp), DIMENSION(npoly), INTENT(out) :: res 00470 REAL(dp), DIMENSION(0:2, 0:2), 00471 INTENT(in) :: h, h_inv 00472 REAL(dp), DIMENSION(0:, 0:, 0:), 00473 INTENT(inout) :: grid 00474 REAL(dp), DIMENSION(:), INTENT(in) :: poly 00475 REAL(dp), INTENT(in) :: alphai 00476 REAL(dp), DIMENSION(0:2), INTENT(in) :: posi 00477 REAL(dp), INTENT(in) :: max_r2 00478 INTEGER, DIMENSION(0:2), INTENT(in) :: periodic 00479 INTEGER, DIMENSION(0:2), INTENT(in), 00480 OPTIONAL :: gdim 00481 INTEGER, DIMENSION(2, 0:2), INTENT(in), 00482 OPTIONAL :: local_bounds 00483 INTEGER, DIMENSION(0:2), INTENT(in), 00484 OPTIONAL :: local_shift 00485 REAL(dp), DIMENSION(0:2), INTENT(in), 00486 OPTIONAL :: poly_shift 00487 REAL(dp), INTENT(in), OPTIONAL :: scale 00488 TYPE(cp_error_type), INTENT(inout) :: error 00489 00490 CHARACTER(len=*), PARAMETER :: routineN = 'integrateGauss', 00491 routineP=moduleN//':'//routineN 00492 00493 #define FMG_INTEGRATE 00494 #include "colloc_int_body.f90" 00495 #undef FMG_INTEGRATE 00496 END SUBROUTINE 00497 00498 ! ***************************************************************************** 00504 SUBROUTINE integrateGaussFull(h,h_inv,grid,poly,alphai,posi,max_r2,& 00505 periodic,gdim,local_bounds,local_shift,poly_shift,scale,error) 00506 REAL(dp), DIMENSION(0:2, 0:2), 00507 INTENT(in) :: h, h_inv 00508 REAL(dp), DIMENSION(0:, 0:, 0:), 00509 INTENT(inout) :: grid 00510 REAL(dp), DIMENSION(:), INTENT(out) :: poly 00511 REAL(dp), INTENT(in) :: alphai 00512 REAL(dp), DIMENSION(0:2), INTENT(in) :: posi 00513 REAL(dp), INTENT(in) :: max_r2 00514 INTEGER, DIMENSION(0:2), INTENT(in) :: periodic 00515 INTEGER, DIMENSION(0:2), INTENT(in), 00516 OPTIONAL :: gdim 00517 INTEGER, DIMENSION(2, 0:2), INTENT(in), 00518 OPTIONAL :: local_bounds 00519 INTEGER, DIMENSION(0:2), INTENT(in), 00520 OPTIONAL :: local_shift 00521 REAL(dp), DIMENSION(0:2), INTENT(in), 00522 OPTIONAL :: poly_shift 00523 REAL(dp), INTENT(in), OPTIONAL :: scale 00524 TYPE(cp_error_type), INTENT(inout) :: error 00525 00526 CHARACTER(len=*), PARAMETER :: routineN = 'integrateGaussFull', 00527 routineP=moduleN//':'//routineN 00528 00529 #define FMG_INTEGRATE_FULL 00530 #include "colloc_int_body.f90" 00531 #undef FMG_INTEGRATE_FULL 00532 END SUBROUTINE 00533 00534 ! ***************************************************************************** 00540 SUBROUTINE integrateGaussFullFlat(h,h_inv,grid,ngpts,ldim,poly,alphai,posi,max_r2,& 00541 periodic,gdim,local_bounds,local_shift,poly_shift,scale,error) 00542 REAL(dp), DIMENSION(0:2, 0:2), 00543 INTENT(in) :: h, h_inv 00544 INTEGER, INTENT(in) :: ngpts,ldim(0:2) 00545 REAL(dp), DIMENSION(IF_CHECK(ngpts,*)), 00546 INTENT(inout) :: grid 00547 REAL(dp), DIMENSION(:), INTENT(out) :: poly 00548 REAL(dp), INTENT(in) :: alphai 00549 REAL(dp), DIMENSION(0:2), INTENT(in) :: posi 00550 REAL(dp), INTENT(in) :: max_r2 00551 INTEGER, DIMENSION(0:2), INTENT(in) :: periodic 00552 INTEGER, DIMENSION(0:2), INTENT(in), 00553 OPTIONAL :: gdim 00554 INTEGER, DIMENSION(2, 0:2), INTENT(in), 00555 OPTIONAL :: local_bounds 00556 INTEGER, DIMENSION(0:2), INTENT(in), 00557 OPTIONAL :: local_shift 00558 REAL(dp), DIMENSION(0:2), INTENT(in), 00559 OPTIONAL :: poly_shift 00560 REAL(dp), INTENT(in), OPTIONAL :: scale 00561 TYPE(cp_error_type), INTENT(inout) :: error 00562 00563 CHARACTER(len=*), PARAMETER :: routineN = 'integrateGaussFullFlat', 00564 routineP=moduleN//':'//routineN 00565 00566 #define FM_FLAT_GRID 00567 #define FMG_INTEGRATE_FULL 00568 #include "colloc_int_body.f90" 00569 #undef FMG_INTEGRATE_FULL 00570 #undef FM_FLAT_GRID 00571 END SUBROUTINE 00572 00573 ! ***************************************************************************** 00576 SUBROUTINE collocGauss_safe(h,h_inv,grid,poly,alphai,posi,max_r2,& 00577 periodic,gdim,local_bounds,local_shift,error) 00578 REAL(dp), DIMENSION(0:2, 0:2), 00579 INTENT(in) :: h, h_inv 00580 REAL(dp), DIMENSION(0:, 0:, 0:), 00581 INTENT(inout) :: grid 00582 REAL(dp), DIMENSION(:), INTENT(in) :: poly 00583 REAL(dp), INTENT(in) :: alphai 00584 REAL(dp), DIMENSION(0:2), INTENT(in) :: posi 00585 REAL(dp), INTENT(in) :: max_r2 00586 INTEGER, DIMENSION(0:2), INTENT(in) :: periodic 00587 INTEGER, DIMENSION(0:2), INTENT(in), 00588 OPTIONAL :: gdim 00589 INTEGER, DIMENSION(2, 0:2), INTENT(in), 00590 OPTIONAL :: local_bounds 00591 INTEGER, DIMENSION(0:2), INTENT(in), 00592 OPTIONAL :: local_shift 00593 TYPE(cp_error_type), INTENT(inout) :: error 00594 00595 CHARACTER(len=*), PARAMETER :: routineN = 'collocGauss_safe', 00596 routineP = moduleN//':'//routineN 00597 00598 INTEGER :: bounds(0:1,0:2), i, idir, ii, 00599 ij, ik, j, k, 00600 l_bounds(0:1,0:2) 00601 INTEGER, DIMENSION(0:2) :: cellShift, l_shift, lb, ndim, 00602 shiftPos, ub 00603 REAL(dp) :: f, max_miss, maxr2, maxr2_0, 00604 maxr2_2, min_miss, p_v(1), 00605 r_2, scaled_h(0:2,0:2) 00606 REAL(dp), DIMENSION(0:2) :: normD, r, resPosReal, ri, 00607 riPos, rj, rPos, wrPos 00608 00609 IF (PRESENT(gdim)) THEN 00610 ndim=gdim 00611 ELSE 00612 ndim(0)=SIZE(grid,1) 00613 ndim(1)=SIZE(grid,2) 00614 ndim(2)=SIZE(grid,3) 00615 END IF 00616 IF (PRESENT(local_bounds)) THEN 00617 l_bounds=local_bounds 00618 ELSE 00619 l_bounds(0,:)=0 00620 l_bounds(1,:)=ndim-1 00621 END IF 00622 IF (PRESENT(local_shift)) THEN 00623 l_shift=local_shift 00624 ELSE 00625 l_shift=0 ! l_bounds(0,:) 00626 END IF 00627 lb=l_bounds(0,:) 00628 ub=l_bounds(1,:)-l_bounds(0,:) 00629 00630 rPos=0 00631 DO j=0,2 00632 DO i=0,2 00633 rPos(i)=rPos(i)+h_inv(i,j)*posi(j) 00634 END DO 00635 END DO 00636 cellShift=FLOOR(rPos)*periodic 00637 wrPos=rPos-cellShift 00638 riPos=wrPos*ndim 00639 shiftPos=FLOOR(riPos+0.5_dp) 00640 normD=1.0_dp/REAL(ndim,dp) 00641 resPosReal=0.0_dp 00642 DO j=0,2 00643 DO i=0,2 00644 resPosReal(i)=resPosReal(i)+h(i,j)*(wrPos(j)-normD(j)*REAL(shiftPos(j),dp)) 00645 END DO 00646 END DO 00647 00648 IF (ALL(poly==0.0_dp)) RETURN ! remove if derivs 00649 00650 maxr2=0.0_dp 00651 DO j=0,2 00652 DO i=0,2 00653 maxr2=maxr2+(h(i,j)*normD(j))**2 00654 END DO 00655 END DO 00656 maxr2_2=maxr2+0.2_dp/alphai 00657 maxr2_0=0.2_dp/alphai 00658 maxr2=MAX(max_r2,maxr2) 00659 maxr2_2=maxr2+maxr2_2+2.0_dp*SQRT(maxr2_2*maxr2) 00660 maxr2_0=maxr2+maxr2_0-2.0_dp*SQRT(maxr2_0*maxr2) 00661 00662 DO j=0,2 00663 DO i=0,2 00664 scaled_h(i,j)=h(i,j)*normD(j) 00665 END DO 00666 END DO 00667 00668 ! sqSize=1+floor(ndim*sqrt(maxr2/minimum.reduce(numpy.linalg.eig(dot(transpose(h),h))[0]))).astype(int) 00669 bounds=calcBox(h,h_inv,posi,max_r2,periodic,gdim,error) 00670 bounds(0,:)=bounds(0,:)-2 00671 bounds(1,:)=bounds(1,:)+2 00672 DO i=0,2 00673 IF (periodic(i)==0) THEN 00674 bounds(0,i)=MAX(l_bounds(0,i),bounds(0,i)) 00675 bounds(1,i)=MIN(l_bounds(1,i),bounds(1,i)) 00676 END IF 00677 END DO 00678 p_v=0 00679 max_miss=-1.0_dp 00680 min_miss=-1.0_dp 00681 DO i=bounds(0,0),bounds(1,0) 00682 ii=MODULO(i-lb(0),ndim(0)) 00683 IF (ii>ub(0)) CYCLE 00684 ii=ii+l_shift(0) 00685 DO idir=0,2 00686 ri(idir)=scaled_h(idir,0)*i-posi(idir) 00687 END DO 00688 DO j=bounds(0,1),bounds(1,1) 00689 ij=MODULO(j-lb(1),ndim(1)) 00690 IF (ij>ub(1)) CYCLE 00691 ij=ij+l_shift(1) 00692 DO idir=0,2 00693 rj(idir)=scaled_h(idir,1)*j+ri(idir) 00694 END DO 00695 DO k=bounds(0,2),bounds(1,2) 00696 ik=MODULO(k-lb(2),ndim(2)) 00697 IF (ik>ub(2)) CYCLE 00698 ik=ik+l_shift(2) 00699 r_2=0.0_dp 00700 DO idir=0,2 00701 r(idir)=scaled_h(idir,2)*k+rj(idir) 00702 r_2=r_2+r(idir)**2 00703 END DO 00704 IF (r_2<maxr2_2) THEN 00705 f=EXP(-alphai*r_2) 00706 CALL poly_eval3(poly,r(0)+posi(0),r(1)+posi(1),r(2)+posi(2),p_v,error=error) 00707 IF (r_2<=maxr2) THEN 00708 IF (r_2>maxr2_0) min_miss=MAX(min_miss,ABS(p_v(1)*f)) 00709 grid(ii,ij,ik)=grid(ii,ij,ik)+p_v(1)*f 00710 ELSE 00711 max_miss=MAX(max_miss,ABS(p_v(1)*f)) 00712 END IF 00713 END IF 00714 END DO 00715 END DO 00716 END DO 00717 PRINT *,'miss2',min_miss,max_miss 00718 END SUBROUTINE 00719 00720 END MODULE gauss_colloc
1.7.3