|
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 ! ***************************************************************************** 00012 MODULE spme 00013 00014 USE atomic_kind_types, ONLY: get_atomic_kind 00015 USE atprop_types, ONLY: atprop_type 00016 USE bibliography, ONLY: Essmann1995,& 00017 cite_reference 00018 USE cell_types, ONLY: cell_type,& 00019 real_to_scaled 00020 USE dgs, ONLY: dg_sum_patch,& 00021 dg_sum_patch_force_1d,& 00022 dg_sum_patch_force_3d 00023 USE ewald_environment_types, ONLY: ewald_env_get,& 00024 ewald_environment_type 00025 USE ewald_pw_types, ONLY: ewald_pw_get,& 00026 ewald_pw_type 00027 USE f77_blas 00028 USE kinds, ONLY: dp 00029 USE mathconstants, ONLY: fourpi 00030 USE particle_types, ONLY: particle_type 00031 USE pme_tools, ONLY: get_center,& 00032 set_list 00033 USE pw_grid_types, ONLY: pw_grid_type 00034 USE pw_grids, ONLY: get_pw_grid_info 00035 USE pw_methods, ONLY: pw_integral_a2b,& 00036 pw_transfer 00037 USE pw_poisson_methods, ONLY: pw_poisson_rebuild,& 00038 pw_poisson_solve 00039 USE pw_poisson_types, ONLY: greens_fn_type,& 00040 pw_poisson_type 00041 USE pw_pool_types, ONLY: pw_pool_create_pw,& 00042 pw_pool_give_back_pw,& 00043 pw_pool_type 00044 USE pw_types, ONLY: COMPLEXDATA1D,& 00045 REALDATA3D,& 00046 REALSPACE,& 00047 RECIPROCALSPACE,& 00048 pw_p_type,& 00049 pw_type 00050 USE realspace_grid_types, ONLY: & 00051 pw2rs, realspace_grid_desc_type, realspace_grid_p_type, & 00052 realspace_grid_type, rs2pw, rs_grid_create, rs_grid_release, & 00053 rs_grid_set_box, rs_grid_zero, rs_pw_transfer 00054 USE shell_potential_types, ONLY: shell_kind_type 00055 USE termination, ONLY: stop_program 00056 USE timings, ONLY: timeset,& 00057 timestop 00058 #include "cp_common_uses.h" 00059 00060 IMPLICIT NONE 00061 00062 PRIVATE 00063 PUBLIC :: spme_evaluate, get_patch 00064 00065 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'spme' 00066 00067 CONTAINS 00068 00069 ! ***************************************************************************** 00074 SUBROUTINE spme_evaluate ( ewald_env, ewald_pw, box, particle_set, & 00075 fg_coulomb, vg_coulomb, pv_g, shell_particle_set, core_particle_set,& 00076 fgshell_coulomb, fgcore_coulomb, use_virial, charges, atprop, error ) 00077 00078 TYPE(ewald_environment_type), POINTER :: ewald_env 00079 TYPE(ewald_pw_type), POINTER :: ewald_pw 00080 TYPE(cell_type), POINTER :: box 00081 TYPE(particle_type), DIMENSION(:), 00082 INTENT(IN) :: particle_set 00083 REAL(KIND=dp), DIMENSION(:, :), 00084 INTENT(OUT) :: fg_coulomb 00085 REAL(KIND=dp), INTENT(OUT) :: vg_coulomb 00086 REAL(KIND=dp), DIMENSION(:, :), 00087 INTENT(OUT) :: pv_g 00088 TYPE(particle_type), DIMENSION(:), 00089 OPTIONAL, POINTER :: shell_particle_set, 00090 core_particle_set 00091 REAL(KIND=dp), DIMENSION(:, :), 00092 INTENT(OUT), OPTIONAL :: fgshell_coulomb, 00093 fgcore_coulomb 00094 LOGICAL, INTENT(IN) :: use_virial 00095 REAL(KIND=dp), DIMENSION(:), OPTIONAL, 00096 POINTER :: charges 00097 TYPE(atprop_type), POINTER :: atprop 00098 TYPE(cp_error_type), INTENT(inout) :: error 00099 00100 CHARACTER(len=*), PARAMETER :: routineN = 'spme_evaluate', 00101 routineP = moduleN//':'//routineN 00102 00103 INTEGER :: group, handle, i, ipart, j, 00104 n, ncore, npart, nshell, 00105 o_spline, p1, p1_shell, stat 00106 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: center, core_center, 00107 shell_center 00108 INTEGER, DIMENSION(3) :: npts 00109 LOGICAL :: do_shell, failure 00110 REAL(KIND=dp) :: alpha, dvols, fat1, ffa 00111 REAL(KIND=dp), ALLOCATABLE, 00112 DIMENSION(:, :, :) :: rhos 00113 REAL(KIND=dp), DIMENSION(3) :: fat 00114 REAL(KIND=dp), DIMENSION(3, 3) :: f_stress, h_stress 00115 TYPE(greens_fn_type), POINTER :: green 00116 TYPE(pw_grid_type), POINTER :: grid_spme 00117 TYPE(pw_p_type), DIMENSION(3) :: dphi_g 00118 TYPE(pw_poisson_type), POINTER :: poisson_env 00119 TYPE(pw_pool_type), POINTER :: pw_pool 00120 TYPE(pw_type), POINTER :: phi_g, rhob_g, rhob_r 00121 TYPE(realspace_grid_desc_type), POINTER :: rs_desc 00122 TYPE(realspace_grid_p_type), 00123 DIMENSION(:), POINTER :: drpot 00124 TYPE(realspace_grid_type), POINTER :: rden, rpot 00125 00126 failure = .FALSE. 00127 NULLIFY(drpot, grid_spme, green, poisson_env, phi_g, rhob_g, rhob_r, pw_pool, rden, rpot) 00128 CALL timeset ( routineN, handle ) 00129 CALL cite_reference(Essmann1995) 00130 00131 !-------------- INITIALISATION --------------------- 00132 CALL ewald_env_get ( ewald_env, alpha=alpha, o_spline = o_spline,& 00133 group = group, error=error) 00134 CALL ewald_pw_get ( ewald_pw, pw_big_pool=pw_pool, rs_desc=rs_desc,& 00135 poisson_env=poisson_env) 00136 CALL pw_poisson_rebuild(poisson_env,error=error) 00137 green => poisson_env%green_fft 00138 grid_spme => pw_pool % pw_grid 00139 00140 npart = SIZE ( particle_set ) 00141 00142 CALL get_pw_grid_info(grid_spme,npts=npts,dvol=dvols,error=error) 00143 00144 n = o_spline 00145 ALLOCATE ( rhos(n,n,n), STAT = stat ) 00146 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00147 CALL rs_grid_create (rden, rs_desc, error=error) 00148 CALL rs_grid_set_box ( grid_spme, rs=rden, error=error ) 00149 CALL rs_grid_zero ( rden ) 00150 00151 ALLOCATE ( center(3, npart), STAT = stat ) 00152 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00153 CALL get_center ( particle_set, box, center, npts, n ) 00154 IF(PRESENT(shell_particle_set).AND.(PRESENT(core_particle_set))) THEN 00155 CPPostcondition(ASSOCIATED(shell_particle_set),cp_failure_level,routineP,error,failure) 00156 CPPostcondition(ASSOCIATED(core_particle_set),cp_failure_level,routineP,error,failure) 00157 nshell=SIZE(shell_particle_set) 00158 ncore =SIZE(core_particle_set) 00159 CPPostcondition(nshell==ncore,cp_failure_level,routineP,error,failure) 00160 ALLOCATE ( shell_center ( 3, nshell ), STAT = stat ) 00161 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00162 CALL get_center ( shell_particle_set, box, shell_center, npts, n ) 00163 ALLOCATE ( core_center ( 3, nshell ), STAT = stat ) 00164 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00165 CALL get_center ( core_particle_set, box, core_center, npts, n ) 00166 END IF 00167 00168 !-------------- DENSITY CALCULATION ---------------- 00169 ipart = 0 00170 ! Particles 00171 DO 00172 CALL set_list ( particle_set, npart, center, p1, rden, ipart ) 00173 IF (p1 == 0) EXIT 00174 00175 do_shell = (particle_set(p1)%shell_index/=0) 00176 IF(do_shell) CYCLE 00177 ! calculate function on small boxes 00178 CALL get_patch ( particle_set, box, green, npts, p1, rhos, is_core=.FALSE.,& 00179 is_shell=.FALSE., unit_charge=.FALSE., charges=charges) 00180 00181 ! add boxes to real space grid (big box) 00182 CALL dg_sum_patch ( rden, rhos, center(:,p1) ) 00183 END DO 00184 ! Shell-Model 00185 IF(PRESENT(shell_particle_set).AND.PRESENT(core_particle_set)) THEN 00186 ipart = 0 00187 DO 00188 ! calculate function on small boxes 00189 CALL set_list ( shell_particle_set, nshell, shell_center, p1_shell,& 00190 rden, ipart ) 00191 IF (p1_shell == 0) EXIT 00192 CALL get_patch(shell_particle_set, box, green, npts, p1_shell, rhos,& 00193 is_core=.FALSE., is_shell=.TRUE., unit_charge=.FALSE.) 00194 00195 ! add boxes to real space grid (big box) 00196 CALL dg_sum_patch ( rden, rhos, shell_center(:,p1_shell) ) 00197 END DO 00198 ipart = 0 00199 DO 00200 ! calculate function on small boxes 00201 CALL set_list ( core_particle_set, nshell, core_center, p1_shell,& 00202 rden, ipart ) 00203 IF ( p1_shell == 0 ) EXIT 00204 CALL get_patch ( core_particle_set, box, green, npts, p1_shell, rhos,& 00205 is_core=.TRUE., is_shell=.FALSE., unit_charge=.FALSE.) 00206 00207 ! add boxes to real space grid (big box) 00208 CALL dg_sum_patch ( rden, rhos, core_center(:,p1_shell) ) 00209 END DO 00210 END IF 00211 !----------- END OF DENSITY CALCULATION ------------- 00212 00213 CALL pw_pool_create_pw ( pw_pool, rhob_r, use_data = REALDATA3D, & 00214 in_space = REALSPACE ,error=error) 00215 CALL rs_pw_transfer ( rden, rhob_r, rs2pw, error=error) 00216 ! transform density to G space and add charge function 00217 CALL pw_pool_create_pw ( pw_pool, rhob_g, & 00218 use_data = COMPLEXDATA1D, & 00219 in_space = RECIPROCALSPACE ,error=error) 00220 CALL pw_transfer ( rhob_r, rhob_g, error=error) 00221 ! update charge function 00222 rhob_g % cc = rhob_g % cc * green % p3m_charge % cr 00223 00224 !-------------- ELECTROSTATIC CALCULATION ----------- 00225 ! allocate intermediate arrays 00226 DO i = 1, 3 00227 NULLIFY(dphi_g(i)%pw) 00228 CALL pw_pool_create_pw ( pw_pool, dphi_g ( i )%pw,& 00229 use_data = COMPLEXDATA1D,& 00230 in_space = RECIPROCALSPACE,& 00231 error=error) 00232 END DO 00233 CALL pw_pool_create_pw ( pw_pool, phi_g, & 00234 use_data = COMPLEXDATA1D,& 00235 in_space = RECIPROCALSPACE,& 00236 error=error) 00237 CALL pw_poisson_solve ( poisson_env, rhob_g, vg_coulomb, phi_g, dphi_g, & 00238 h_stress ,error=error) 00239 ! Atomic Energy 00240 IF (atprop%energy .OR. atprop%stress) THEN 00241 CALL rs_grid_create (rpot, rs_desc, error=error) 00242 phi_g%cc = phi_g%cc * green%p3m_charge%cr 00243 CALL pw_transfer ( phi_g, rhob_r, error=error) 00244 CALL rs_pw_transfer (rpot, rhob_r, pw2rs, error=error) 00245 ipart = 0 00246 DO 00247 CALL set_list ( particle_set, npart, center, p1, rden, ipart ) 00248 IF ( p1 == 0 ) EXIT 00249 do_shell = (particle_set(p1)%shell_index /= 0) 00250 IF (do_shell) CYCLE 00251 ! calculate function on small boxes 00252 CALL get_patch ( particle_set, box, green, grid_spme % npts, p1, rhos, is_core=.FALSE.,& 00253 is_shell=.FALSE., unit_charge=.FALSE., charges=charges) 00254 ! integrate box and potential 00255 CALL dg_sum_patch_force_1d ( rpot, rhos, center(:,p1), fat1 ) 00256 IF (atprop%energy) THEN 00257 atprop%atener(p1) = atprop%atener(p1) + 0.5_dp*fat1*dvols 00258 END IF 00259 IF (atprop%stress) THEN 00260 atprop%atstress(1,1,p1) = atprop%atstress(1,1,p1) + 0.5_dp*fat1*dvols 00261 atprop%atstress(2,2,p1) = atprop%atstress(2,2,p1) + 0.5_dp*fat1*dvols 00262 atprop%atstress(3,3,p1) = atprop%atstress(3,3,p1) + 0.5_dp*fat1*dvols 00263 END IF 00264 END DO 00265 ! Core-shell model 00266 IF (PRESENT(shell_particle_set).AND.PRESENT(core_particle_set)) THEN 00267 ipart = 0 00268 DO 00269 CALL set_list(shell_particle_set,nshell,shell_center,p1_shell,rden,ipart) 00270 IF (p1_shell == 0) EXIT 00271 CALL get_patch(shell_particle_set,box,green,npts,p1_shell,rhos,& 00272 is_core=.FALSE.,is_shell=.TRUE.,unit_charge=.FALSE.) 00273 CALL dg_sum_patch_force_1d(rpot,rhos,shell_center(:,p1_shell),fat1) 00274 p1 = shell_particle_set(p1_shell)%atom_index 00275 IF (atprop%energy) THEN 00276 atprop%atener(p1) = atprop%atener(p1) + 0.5_dp*fat1*dvols 00277 END IF 00278 END DO 00279 ipart = 0 00280 DO 00281 CALL set_list(core_particle_set,nshell,core_center,p1_shell,rden,ipart) 00282 IF (p1_shell == 0) EXIT 00283 CALL get_patch(core_particle_set,box,green,npts,p1_shell,rhos,& 00284 is_core=.TRUE.,is_shell=.FALSE.,unit_charge=.FALSE.) 00285 CALL dg_sum_patch_force_1d(rpot,rhos,core_center(:,p1_shell),fat1) 00286 p1 = core_particle_set(p1_shell)%atom_index 00287 IF (atprop%energy) THEN 00288 atprop%atener(p1) = atprop%atener(p1) + 0.5_dp*fat1*dvols 00289 END IF 00290 END DO 00291 END IF 00292 CALL rs_grid_release ( rpot, error=error) 00293 END IF 00294 00295 CALL pw_pool_give_back_pw ( pw_pool, phi_g ,error=error) 00296 CALL pw_pool_give_back_pw ( pw_pool, rhob_g ,error=error) 00297 !---------- END OF ELECTROSTATIC CALCULATION -------- 00298 00299 !------------- STRESS TENSOR CALCULATION ------------ 00300 IF (use_virial) THEN 00301 DO i = 1, 3 00302 DO j = i, 3 00303 f_stress ( i, j ) = pw_integral_a2b ( dphi_g ( i ) % pw, & 00304 dphi_g ( j ) % pw, error=error) 00305 f_stress ( j, i ) = f_stress ( i, j ) 00306 END DO 00307 END DO 00308 ffa = ( 1.0_dp / fourpi ) * ( 0.5_dp / alpha ) ** 2 00309 f_stress = -ffa * f_stress 00310 pv_g = h_stress + f_stress 00311 END IF 00312 !--------END OF STRESS TENSOR CALCULATION ----------- 00313 ! move derivative of potential to real space grid and 00314 ! multiply by charge function in g-space 00315 ALLOCATE ( drpot(1:3), STAT=stat ) 00316 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00317 DO i = 1, 3 00318 CALL rs_grid_create( drpot(i)%rs_grid, rs_desc, error=error ) 00319 CALL rs_grid_set_box ( grid_spme, rs=drpot(i)%rs_grid, error=error ) 00320 dphi_g ( i ) % pw % cc = dphi_g ( i ) % pw % cc * green % p3m_charge % cr 00321 CALL pw_transfer ( dphi_g ( i )%pw, rhob_r, error=error) 00322 CALL pw_pool_give_back_pw ( pw_pool, dphi_g ( i )%pw ,error=error) 00323 CALL rs_pw_transfer ( drpot ( i ) % rs_grid, rhob_r, pw2rs,error=error) 00324 END DO 00325 00326 CALL pw_pool_give_back_pw ( pw_pool, rhob_r ,error=error) 00327 !----------------- FORCE CALCULATION ---------------- 00328 ! initialize the forces 00329 fg_coulomb = 0.0_dp 00330 ! Particles 00331 ipart = 0 00332 DO 00333 CALL set_list ( particle_set, npart, center, p1, rden, ipart ) 00334 IF ( p1 == 0 ) EXIT 00335 00336 do_shell = (particle_set(p1)%shell_index/=0) 00337 IF(do_shell) CYCLE 00338 ! calculate function on small boxes 00339 CALL get_patch ( particle_set, box, green, grid_spme % npts, p1, rhos, is_core=.FALSE.,& 00340 is_shell=.FALSE., unit_charge=.FALSE., charges=charges) 00341 00342 ! add boxes to real space grid (big box) 00343 CALL dg_sum_patch_force_3d ( drpot, rhos, center(:,p1), fat ) 00344 fg_coulomb ( 1, p1 ) = fg_coulomb ( 1, p1 ) - fat ( 1 ) * dvols 00345 fg_coulomb ( 2, p1 ) = fg_coulomb ( 2, p1 ) - fat ( 2 ) * dvols 00346 fg_coulomb ( 3, p1 ) = fg_coulomb ( 3, p1 ) - fat ( 3 ) * dvols 00347 END DO 00348 ! Shell-Model 00349 IF(PRESENT(shell_particle_set).AND.(PRESENT(core_particle_set))) THEN 00350 IF(PRESENT(fgshell_coulomb)) THEN 00351 ipart = 0 00352 fgshell_coulomb = 0.0_dp 00353 DO 00354 ! calculate function on small boxes 00355 CALL set_list ( shell_particle_set, nshell, shell_center, p1_shell,& 00356 rden, ipart ) 00357 IF ( p1_shell == 0 ) EXIT 00358 00359 CALL get_patch ( shell_particle_set, box, green, grid_spme % npts, & 00360 p1_shell, rhos, is_core=.FALSE., is_shell=.TRUE., unit_charge=.FALSE.) 00361 00362 ! add boxes to real space grid (big box) 00363 CALL dg_sum_patch_force_3d ( drpot, rhos, shell_center(:,p1_shell), fat ) 00364 fgshell_coulomb ( 1, p1_shell ) = fgshell_coulomb ( 1, p1_shell ) - fat ( 1 ) * dvols 00365 fgshell_coulomb ( 2, p1_shell ) = fgshell_coulomb ( 2, p1_shell ) - fat ( 2 ) * dvols 00366 fgshell_coulomb ( 3, p1_shell ) = fgshell_coulomb ( 3, p1_shell ) - fat ( 3 ) * dvols 00367 00368 END DO 00369 END IF 00370 IF(PRESENT(fgcore_coulomb)) THEN 00371 ipart = 0 00372 fgcore_coulomb = 0.0_dp 00373 DO 00374 ! calculate function on small boxes 00375 CALL set_list ( core_particle_set, nshell, core_center, p1_shell,& 00376 rden, ipart ) 00377 IF ( p1_shell == 0 ) EXIT 00378 00379 CALL get_patch ( core_particle_set, box, green, grid_spme % npts, & 00380 p1_shell, rhos, is_core=.TRUE., is_shell=.FALSE., unit_charge=.FALSE.) 00381 00382 ! add boxes to real space grid (big box) 00383 CALL dg_sum_patch_force_3d ( drpot, rhos, core_center(:,p1_shell), fat ) 00384 fgcore_coulomb ( 1, p1_shell ) = fgcore_coulomb ( 1, p1_shell ) - fat ( 1 ) * dvols 00385 fgcore_coulomb ( 2, p1_shell ) = fgcore_coulomb ( 2, p1_shell ) - fat ( 2 ) * dvols 00386 fgcore_coulomb ( 3, p1_shell ) = fgcore_coulomb ( 3, p1_shell ) - fat ( 3 ) * dvols 00387 END DO 00388 END IF 00389 00390 END IF 00391 !--------------END OF FORCE CALCULATION ------------- 00392 00393 !------------------CLEANING UP ---------------------- 00394 CALL rs_grid_release (rden, error=error) 00395 IF (ASSOCIATED(drpot)) THEN 00396 DO i = 1, 3 00397 CALL rs_grid_release ( drpot(i)%rs_grid, error=error) 00398 END DO 00399 DEALLOCATE ( drpot, STAT = stat ) 00400 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00401 END IF 00402 00403 DEALLOCATE ( rhos, STAT = stat ) 00404 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00405 DEALLOCATE ( center, STAT = stat ) 00406 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00407 IF(ALLOCATED(shell_center)) THEN 00408 DEALLOCATE ( shell_center, STAT = stat ) 00409 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00410 END IF 00411 IF(ALLOCATED(core_center)) THEN 00412 DEALLOCATE ( core_center, STAT = stat ) 00413 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00414 END IF 00415 CALL timestop(handle) 00416 00417 END SUBROUTINE spme_evaluate 00418 00419 ! ***************************************************************************** 00425 SUBROUTINE get_patch ( part, box, green, npts, p, rhos, is_core, is_shell,& 00426 unit_charge, charges ) 00427 00428 TYPE(particle_type), DIMENSION(:), 00429 INTENT(IN) :: part 00430 TYPE(cell_type), POINTER :: box 00431 TYPE(greens_fn_type), POINTER :: green 00432 INTEGER, DIMENSION(3), INTENT(IN) :: npts 00433 INTEGER, INTENT(IN) :: p 00434 REAL(KIND=dp), DIMENSION(:, :, :), 00435 INTENT(OUT) :: rhos 00436 LOGICAL, INTENT(IN) :: is_core, is_shell, unit_charge 00437 REAL(KIND=dp), DIMENSION(:), OPTIONAL, 00438 POINTER :: charges 00439 00440 CHARACTER(len=*), PARAMETER :: routineN = 'get_patch', 00441 routineP = moduleN//':'//routineN 00442 00443 INTEGER :: nbox 00444 LOGICAL :: use_charge_array 00445 REAL(KIND=dp) :: q 00446 REAL(KIND=dp), DIMENSION(3) :: delta, r 00447 TYPE(shell_kind_type), POINTER :: shell 00448 00449 NULLIFY(shell) 00450 use_charge_array=PRESENT(charges) 00451 IF (use_charge_array) use_charge_array=ASSOCIATED(charges) 00452 IF(is_core .AND. is_shell) THEN 00453 CALL stop_program(routineN,moduleN,__LINE__,& 00454 "Shell-model: cannot be core and shell simultaneously") 00455 END IF 00456 00457 nbox = SIZE ( rhos, 1 ) 00458 r = part(p)% r 00459 q = 1.0_dp 00460 IF (.NOT.unit_charge) THEN 00461 IF (is_core) THEN 00462 CALL get_atomic_kind ( atomic_kind=part(p)%atomic_kind, shell=shell) 00463 q = shell%charge_core 00464 ELSE IF (is_shell) THEN 00465 CALL get_atomic_kind ( atomic_kind=part(p)%atomic_kind, shell=shell) 00466 q = shell%charge_shell 00467 ELSE 00468 CALL get_atomic_kind ( atomic_kind=part(p)%atomic_kind, qeff=q) 00469 END IF 00470 IF (use_charge_array) q=charges(p) 00471 END IF 00472 CALL get_delta ( box, r, npts, delta, nbox ) 00473 CALL spme_get_patch ( rhos, nbox, delta, q, green % p3m_coeff ) 00474 00475 END SUBROUTINE get_patch 00476 00477 ! ***************************************************************************** 00483 SUBROUTINE spme_get_patch ( rhos, n, delta, q, coeff ) 00484 00485 REAL(KIND=dp), DIMENSION(:, :, :), 00486 INTENT(OUT) :: rhos 00487 INTEGER, INTENT(IN) :: n 00488 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: delta 00489 REAL(KIND=dp), INTENT(IN) :: q 00490 REAL(KIND=dp), 00491 DIMENSION(-(n-1):n-1, 0:n-1), 00492 INTENT(IN) :: coeff 00493 00494 CHARACTER(len=*), PARAMETER :: routineN = 'spme_get_patch', 00495 routineP = moduleN//':'//routineN 00496 INTEGER, PARAMETER :: nmax = 12 00497 00498 INTEGER :: i, i1, i2, i3, j, l 00499 REAL(KIND=dp) :: r2, r3 00500 REAL(KIND=dp), DIMENSION(3, -nmax:nmax) :: w_assign 00501 REAL(KIND=dp), DIMENSION(3, 0:nmax-1) :: deltal 00502 REAL(KIND=dp), DIMENSION(3, 1:nmax) :: f_assign 00503 00504 IF ( n > nmax ) THEN 00505 CALL stop_program(routineN,moduleN,__LINE__,"nmax value too small") 00506 END IF 00507 ! calculate the assignment function values and 00508 ! the charges on the grid (small box) 00509 00510 deltal ( 1, 0 ) = 1.0_dp 00511 deltal ( 2, 0 ) = 1.0_dp 00512 deltal ( 3, 0 ) = 1.0_dp 00513 DO l = 1, n-1 00514 deltal ( 1, l ) = deltal ( 1, l-1 ) * delta ( 1 ) 00515 deltal ( 2, l ) = deltal ( 2, l-1 ) * delta ( 2 ) 00516 deltal ( 3, l ) = deltal ( 3, l-1 ) * delta ( 3 ) 00517 END DO 00518 00519 w_assign = 0.0_dp 00520 DO j = -(n-1), n-1, 2 00521 DO l = 0, n-1 00522 w_assign ( 1, j ) = w_assign ( 1, j ) + coeff ( j, l ) * deltal ( 1, l ) 00523 w_assign ( 2, j ) = w_assign ( 2, j ) + coeff ( j, l ) * deltal ( 2, l ) 00524 w_assign ( 3, j ) = w_assign ( 3, j ) + coeff ( j, l ) * deltal ( 3, l ) 00525 END DO 00526 END DO 00527 DO i = 1, n 00528 j = n + 1 - 2 * i 00529 f_assign (1, i ) = w_assign ( 1, j ) 00530 f_assign (2, i ) = w_assign ( 2, j ) 00531 f_assign (3, i ) = w_assign ( 3, j ) 00532 END DO 00533 00534 DO i3 = 1, n 00535 r3 = q * f_assign ( 3, i3 ) 00536 DO i2 = 1, n 00537 r2 = r3 * f_assign ( 2, i2 ) 00538 DO i1 = 1, n 00539 rhos ( i1, i2, i3 ) = r2 * f_assign ( 1, i1 ) 00540 END DO 00541 END DO 00542 END DO 00543 00544 END SUBROUTINE spme_get_patch 00545 00546 ! ***************************************************************************** 00547 SUBROUTINE get_delta ( box, r, npts, delta, n ) 00548 00549 TYPE(cell_type), POINTER :: box 00550 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: r 00551 INTEGER, DIMENSION(3), INTENT(IN) :: npts 00552 REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: delta 00553 INTEGER, INTENT(IN) :: n 00554 00555 CHARACTER(len=*), PARAMETER :: routineN = 'get_delta', 00556 routineP = moduleN//':'//routineN 00557 00558 INTEGER :: mp 00559 INTEGER, DIMENSION(3) :: center 00560 REAL(KIND=dp) :: rmp 00561 REAL(KIND=dp), DIMENSION(3) :: ca, grid_i, s 00562 00563 mp = MAXVAL ( npts(:) ) 00564 rmp = REAL ( mp,KIND=dp) 00565 ! compute the scaled coordinate of atomi 00566 CALL real_to_scaled(s,r,box) 00567 s = s - REAL ( NINT ( s ),KIND=dp) 00568 00569 ! find the continuous ``grid'' point 00570 grid_i ( 1:3 ) = REAL ( npts ( 1:3 ),KIND=dp) * s ( 1:3 ) 00571 00572 ! find the closest grid point 00573 00574 IF ( MOD ( n, 2 ) == 0 ) THEN 00575 center ( : ) = INT ( grid_i ( : ) + rmp ) - mp 00576 ca ( : ) = REAL ( center ( : ) ) + 0.5_dp 00577 ELSE 00578 center ( : ) = NINT ( grid_i ( : ) ) 00579 ca ( : ) = REAL ( center ( : ) ) 00580 END IF 00581 00582 ! find the distance vector 00583 delta ( : ) = grid_i ( : ) - ca ( : ) 00584 00585 END SUBROUTINE get_delta 00586 00587 END MODULE spme 00588
1.7.3