|
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 ! ***************************************************************************** 00010 MODULE cp_ddapc_methods 00011 USE cell_types, ONLY: cell_type 00012 USE cp_para_types, ONLY: cp_para_env_type 00013 USE erf_fn, ONLY: erfc 00014 USE f77_blas 00015 USE input_constants, ONLY: weight_type_mass,& 00016 weight_type_unit 00017 USE input_section_types, ONLY: section_vals_type,& 00018 section_vals_val_get,& 00019 section_vals_val_set 00020 USE kahan_sum, ONLY: accurate_sum 00021 USE kinds, ONLY: dp 00022 USE mathconstants, ONLY: fourpi,& 00023 oorootpi,& 00024 pi,& 00025 twopi 00026 USE mathlib, ONLY: diamat_all,& 00027 invert_matrix 00028 !NB for paralelizaion of ewald_ddapc_pot() 00029 USE message_passing, ONLY: mp_sum 00030 USE particle_types, ONLY: particle_type 00031 USE pw_spline_utils, ONLY: Eval_Interp_Spl3_pbc 00032 USE pw_types, ONLY: pw_type 00033 USE spherical_harmonics, ONLY: legendre 00034 USE timings, ONLY: timeset,& 00035 timestop 00036 #include "cp_common_uses.h" 00037 00038 IMPLICIT NONE 00039 PRIVATE 00040 LOGICAL, PRIVATE, PARAMETER :: debug_this_module=.FALSE. 00041 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_ddapc_methods' 00042 PUBLIC :: ddapc_eval_gfunc,& 00043 build_b_vector,& 00044 build_der_b_vector,& 00045 build_A_matrix,& 00046 build_der_A_matrix_rows,& 00047 prep_g_dot_rvec_sin_cos,& 00048 cleanup_g_dot_rvec_sin_cos,& 00049 ddapc_eval_AmI,& 00050 ewald_ddapc_pot,& 00051 solvation_ddapc_pot 00052 00053 CONTAINS 00054 00055 ! ***************************************************************************** 00056 SUBROUTINE ddapc_eval_gfunc(gfunc, w, gcut, rho_tot_g, radii, error) 00057 REAL(KIND=dp), DIMENSION(:, :), POINTER :: gfunc 00058 REAL(kind=dp), DIMENSION(:), POINTER :: w 00059 REAL(KIND=dp), INTENT(IN) :: gcut 00060 TYPE(pw_type), POINTER :: rho_tot_g 00061 REAL(kind=dp), DIMENSION(:), POINTER :: radii 00062 TYPE(cp_error_type), INTENT(inout) :: error 00063 00064 CHARACTER(len=*), PARAMETER :: routineN = 'ddapc_eval_gfunc', 00065 routineP = moduleN//':'//routineN 00066 00067 INTEGER :: e_dim, handle, ig, igauss, 00068 s_dim, stat 00069 LOGICAL :: failure 00070 REAL(KIND=dp) :: g2, gcut2, rc, rc2 00071 00072 CALL timeset(routineN,handle) 00073 gcut2 = gcut * gcut 00074 failure = .FALSE. 00075 IF (.NOT.failure) THEN 00076 ! 00077 s_dim = rho_tot_g % pw_grid % first_gne0 00078 e_dim = rho_tot_g % pw_grid % ngpts_cut_local 00079 ALLOCATE(gfunc(s_dim:e_dim,SIZE(radii)), stat=stat) 00080 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00081 ALLOCATE(w(s_dim:e_dim), stat=stat) 00082 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00083 gfunc = 0.0_dp 00084 w = 0.0_dp 00085 DO igauss = 1, SIZE(radii) 00086 rc = radii(igauss) 00087 rc2 = rc*rc 00088 DO ig = s_dim, e_dim 00089 g2 = rho_tot_g % pw_grid % gsq ( ig ) 00090 IF (g2 > gcut2) EXIT 00091 gfunc(ig,igauss) = EXP(-g2*rc2/4.0_dp) 00092 END DO 00093 END DO 00094 DO ig = s_dim, e_dim 00095 g2 = rho_tot_g % pw_grid % gsq ( ig ) 00096 IF (g2 > gcut2) EXIT 00097 w(ig) = fourpi*(g2 -gcut2)**2/(g2*gcut2) 00098 END DO 00099 END IF 00100 CALL timestop(handle) 00101 END SUBROUTINE ddapc_eval_gfunc 00102 00103 ! ***************************************************************************** 00109 SUBROUTINE build_b_vector(bv, gfunc, w, particle_set, radii, rho_tot_g, gcut, error) 00110 REAL(KIND=dp), DIMENSION(:), 00111 INTENT(INOUT) :: bv 00112 REAL(KIND=dp), DIMENSION(:, :), POINTER :: gfunc 00113 REAL(KIND=dp), DIMENSION(:), POINTER :: w 00114 TYPE(particle_type), DIMENSION(:), 00115 POINTER :: particle_set 00116 REAL(KIND=dp), DIMENSION(:), POINTER :: radii 00117 TYPE(pw_type), POINTER :: rho_tot_g 00118 REAL(KIND=dp), INTENT(IN) :: gcut 00119 TYPE(cp_error_type), INTENT(inout) :: error 00120 00121 CHARACTER(len=*), PARAMETER :: routineN = 'build_b_vector', 00122 routineP = moduleN//':'//routineN 00123 00124 COMPLEX(KIND=dp) :: phase 00125 INTEGER :: e_dim, handle, idim, ig, 00126 igauss, igmax, iparticle, 00127 s_dim, stat 00128 LOGICAL :: failure 00129 REAL(KIND=dp) :: arg, g2, gcut2, gvec(3), 00130 rvec(3) 00131 REAL(KIND=dp), DIMENSION(:), POINTER :: my_bv, my_bvw 00132 00133 failure = .FALSE. 00134 CALL timeset(routineN,handle) 00135 IF (.NOT.failure) THEN 00136 NULLIFY(my_bv, my_bvw) 00137 gcut2 = gcut * gcut 00138 s_dim = rho_tot_g % pw_grid % first_gne0 00139 e_dim = rho_tot_g % pw_grid % ngpts_cut_local 00140 igmax = 0 00141 DO ig = s_dim, e_dim 00142 g2 = rho_tot_g % pw_grid % gsq ( ig ) 00143 IF (g2 > gcut2) EXIT 00144 igmax = ig 00145 ENDDO 00146 ALLOCATE(my_bv(s_dim:igmax),stat=stat) 00147 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00148 ALLOCATE(my_bvw(s_dim:igmax),stat=stat) 00149 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00150 ! 00151 DO iparticle = 1, SIZE(particle_set) 00152 rvec = particle_set(iparticle)%r 00153 my_bv = 0.0_dp 00154 DO ig = s_dim, igmax 00155 gvec = rho_tot_g % pw_grid % g (:,ig) 00156 arg = DOT_PRODUCT(gvec,rvec) 00157 phase = CMPLX(COS(arg),-SIN(arg),KIND=dp) 00158 my_bv (ig) = w(ig) * REAL(CONJG(rho_tot_g%cc(ig))*phase,KIND=dp) 00159 END DO 00160 DO igauss = 1, SIZE(radii) 00161 idim = (iparticle-1)*SIZE(radii)+igauss 00162 DO ig = s_dim, igmax 00163 my_bvw(ig) = my_bv(ig)* gfunc(ig,igauss) 00164 END DO 00165 bv (idim) = accurate_sum(my_bvw) 00166 END DO 00167 END DO 00168 DEALLOCATE(my_bvw,stat=stat) 00169 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00170 DEALLOCATE(my_bv,stat=stat) 00171 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00172 END IF 00173 CALL timestop(handle) 00174 END SUBROUTINE build_b_vector 00175 00176 ! ***************************************************************************** 00183 SUBROUTINE build_A_matrix(Am, gfunc, w, particle_set, radii, rho_tot_g, gcut, g_dot_rvec_sin, g_dot_rvec_cos, error) 00184 REAL(KIND=dp), DIMENSION(:, :), 00185 INTENT(INOUT) :: Am 00186 REAL(KIND=dp), DIMENSION(:, :), POINTER :: gfunc 00187 REAL(KIND=dp), DIMENSION(:), POINTER :: w 00188 TYPE(particle_type), DIMENSION(:), 00189 POINTER :: particle_set 00190 REAL(KIND=dp), DIMENSION(:), POINTER :: radii 00191 TYPE(pw_type), POINTER :: rho_tot_g 00192 REAL(KIND=dp), INTENT(IN) :: gcut 00193 REAL(KIND=dp), DIMENSION(:, :), 00194 INTENT(IN) :: g_dot_rvec_sin, g_dot_rvec_cos 00195 TYPE(cp_error_type), INTENT(inout) :: error 00196 00197 CHARACTER(len=*), PARAMETER :: routineN = 'build_A_matrix', 00198 routineP = moduleN//':'//routineN 00199 00200 INTEGER :: e_dim, handle, idim1, idim2, ig, igauss1, igauss2, igmax, 00201 iparticle1, iparticle2, istart_g, s_dim, stat 00202 LOGICAL :: failure 00203 REAL(KIND=dp) :: g2, gcut2, tmp 00204 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: my_Am, my_Amw 00205 REAL(KIND=dp), ALLOCATABLE, 00206 DIMENSION(:, :, :) :: gfunc_sq(:,:,:) 00207 00208 !NB precalculate as many things outside of the innermost loop as possible, in particular w(ig)*gfunc(ig,igauus1)*gfunc(ig,igauss2) 00209 00210 failure = .FALSE. 00211 00212 CALL timeset(routineN,handle) 00213 IF (.NOT. failure) THEN 00214 gcut2 = gcut * gcut 00215 s_dim = rho_tot_g % pw_grid % first_gne0 00216 e_dim = rho_tot_g % pw_grid % ngpts_cut_local 00217 igmax=0 00218 DO ig = s_dim, e_dim 00219 g2 = rho_tot_g % pw_grid % gsq ( ig ) 00220 IF (g2 > gcut2) EXIT 00221 igmax = ig 00222 ENDDO 00223 ALLOCATE(my_Am(s_dim:igmax),stat=stat) 00224 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00225 ALLOCATE(my_Amw(s_dim:igmax),stat=stat) 00226 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00227 ALLOCATE(gfunc_sq(s_dim:igmax,SIZE(radii),SIZE(radii)),stat=stat) 00228 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00229 00230 DO igauss1 = 1, SIZE(radii) 00231 DO igauss2 = 1, SIZE(radii) 00232 gfunc_sq(s_dim:igmax,igauss1,igauss2) = w(s_dim:igmax)*gfunc(s_dim:igmax,igauss1)*gfunc(s_dim:igmax,igauss2) 00233 END DO 00234 END DO 00235 00236 DO iparticle1 = 1, SIZE(particle_set) 00237 DO iparticle2 = iparticle1, SIZE(particle_set) 00238 DO ig = s_dim, igmax 00239 !NB replace explicit dot product and cosine with cos(A+B) formula - much faster 00240 my_Am(ig) = (g_dot_rvec_cos(ig-s_dim+1,iparticle1)*g_dot_rvec_cos(ig-s_dim+1,iparticle2) + & 00241 g_dot_rvec_sin(ig-s_dim+1,iparticle1)*g_dot_rvec_sin(ig-s_dim+1,iparticle2)) 00242 END DO 00243 DO igauss1 = 1, SIZE(radii) 00244 idim1 = (iparticle1-1)*SIZE(radii)+igauss1 00245 istart_g = 1 00246 IF (iparticle2==iparticle1) istart_g = igauss1 00247 DO igauss2 = istart_g, SIZE(radii) 00248 idim2 = (iparticle2-1)*SIZE(radii)+igauss2 00249 my_Amw(s_dim:igmax) = my_Am(s_dim:igmax)*gfunc_sq(s_dim:igmax,igauss1,igauss2) 00250 !NB no loss of accuracy in my test cases 00251 #ifdef ACCURATE_SUMS 00252 tmp = accurate_sum(my_Amw) 00253 #else 00254 tmp = SUM(my_Amw) 00255 #endif 00256 Am (idim2,idim1) = tmp 00257 Am (idim1,idim2) = tmp 00258 END DO 00259 END DO 00260 END DO 00261 END DO 00262 DEALLOCATE(gfunc_sq,stat=stat) 00263 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00264 DEALLOCATE(my_Amw,stat=stat) 00265 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00266 DEALLOCATE(my_Am,stat=stat) 00267 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00268 END IF 00269 CALL timestop(handle) 00270 END SUBROUTINE build_A_matrix 00271 00272 ! ***************************************************************************** 00278 SUBROUTINE build_der_b_vector(dbv, gfunc, w, particle_set, radii, rho_tot_g, gcut, iparticle0, error) 00279 REAL(KIND=dp), DIMENSION(:, :), 00280 INTENT(INOUT) :: dbv 00281 REAL(KIND=dp), DIMENSION(:, :), POINTER :: gfunc 00282 REAL(KIND=dp), DIMENSION(:), POINTER :: w 00283 TYPE(particle_type), DIMENSION(:), 00284 POINTER :: particle_set 00285 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: radii 00286 TYPE(pw_type), POINTER :: rho_tot_g 00287 REAL(KIND=dp), INTENT(IN) :: gcut 00288 INTEGER, INTENT(IN) :: iparticle0 00289 TYPE(cp_error_type), INTENT(inout) :: error 00290 00291 CHARACTER(len=*), PARAMETER :: routineN = 'build_der_b_vector', 00292 routineP = moduleN//':'//routineN 00293 00294 COMPLEX(KIND=dp) :: dphase 00295 INTEGER :: e_dim, handle, idim, ig, 00296 igauss, igmax, iparticle, 00297 s_dim, stat 00298 LOGICAL :: failure 00299 REAL(KIND=dp) :: arg, g2, gcut2 00300 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: my_dbvw 00301 REAL(KIND=dp), ALLOCATABLE, 00302 DIMENSION(:, :) :: my_dbv 00303 REAL(KIND=dp), DIMENSION(3) :: gvec, rvec 00304 00305 failure = .FALSE. 00306 CALL timeset(routineN,handle) 00307 IF (.NOT.failure) THEN 00308 gcut2 = gcut * gcut 00309 s_dim = rho_tot_g % pw_grid % first_gne0 00310 e_dim = rho_tot_g % pw_grid % ngpts_cut_local 00311 igmax = 0 00312 DO ig = s_dim, e_dim 00313 g2 = rho_tot_g % pw_grid % gsq ( ig ) 00314 IF (g2 > gcut2) EXIT 00315 igmax = ig 00316 ENDDO 00317 ALLOCATE(my_dbv(3,s_dim:igmax),stat=stat) 00318 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00319 ALLOCATE(my_dbvw(s_dim:igmax),stat=stat) 00320 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00321 DO iparticle = 1, SIZE(particle_set) 00322 IF (iparticle /= iparticle0) CYCLE 00323 rvec = particle_set(iparticle)%r 00324 DO ig = s_dim, igmax 00325 gvec = rho_tot_g % pw_grid % g (:,ig) 00326 arg = DOT_PRODUCT(gvec,rvec) 00327 dphase = - CMPLX(SIN(arg),COS(arg),KIND=dp) 00328 my_dbv (:,ig) = w(ig) * REAL(CONJG(rho_tot_g%cc(ig))*dphase,KIND=dp) * gvec(:) 00329 END DO 00330 DO igauss = 1, SIZE(radii) 00331 idim = (iparticle-1)*SIZE(radii)+igauss 00332 DO ig = s_dim, igmax 00333 my_dbvw(ig) = my_dbv(1,ig)*gfunc(ig,igauss) 00334 END DO 00335 dbv (idim,1) = accurate_sum(my_dbvw) 00336 DO ig = s_dim, igmax 00337 my_dbvw(ig) = my_dbv(2,ig)*gfunc(ig,igauss) 00338 END DO 00339 dbv (idim,2) = accurate_sum(my_dbvw) 00340 DO ig = s_dim, igmax 00341 my_dbvw(ig) = my_dbv(3,ig)*gfunc(ig,igauss) 00342 END DO 00343 dbv (idim,3) = accurate_sum(my_dbvw) 00344 END DO 00345 END DO 00346 DEALLOCATE(my_dbvw,stat=stat) 00347 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00348 DEALLOCATE(my_dbv,stat=stat) 00349 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00350 END IF 00351 CALL timestop(handle) 00352 END SUBROUTINE build_der_b_vector 00353 00354 ! ***************************************************************************** 00362 SUBROUTINE build_der_A_matrix_rows(dAm, gfunc, w, particle_set, radii,& 00363 rho_tot_g, gcut, iparticle0, nparticles, g_dot_rvec_sin, g_dot_rvec_cos, error) 00364 REAL(KIND=dp), DIMENSION(:, :, :), 00365 INTENT(INOUT) :: dAm 00366 REAL(KIND=dp), DIMENSION(:, :), POINTER :: gfunc 00367 REAL(KIND=dp), DIMENSION(:), POINTER :: w 00368 TYPE(particle_type), DIMENSION(:), 00369 POINTER :: particle_set 00370 REAL(KIND=dp), DIMENSION(:), POINTER :: radii 00371 TYPE(pw_type), POINTER :: rho_tot_g 00372 REAL(KIND=dp), INTENT(IN) :: gcut 00373 INTEGER, INTENT(IN) :: iparticle0, nparticles 00374 REAL(KIND=dp), DIMENSION(:, :), 00375 INTENT(IN) :: g_dot_rvec_sin, g_dot_rvec_cos 00376 TYPE(cp_error_type), INTENT(inout) :: error 00377 00378 CHARACTER(len=*), PARAMETER :: routineN = 'build_der_A_matrix_rows', 00379 routineP = moduleN//':'//routineN 00380 00381 INTEGER :: e_dim, handle, ig, igauss2, 00382 igmax, iparticle1, 00383 iparticle2, s_dim, stat 00384 LOGICAL :: failure 00385 REAL(KIND=dp) :: g2, gcut2 00386 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: arg1_v_sin 00387 00388 !NB calculate derivatives for a block of particles, just the row parts (since derivative matrix is symmetric) 00389 !NB Use DGEMM to speed up calculation, can't do accurate_sum() anymore because dgemm does the sum over g 00390 00391 EXTERNAL DGEMM 00392 REAL(KIND=dp), ALLOCATABLE :: lhs(:,:), rhs(:,:) 00393 INTEGER :: Nr, Np, Ng, icomp, ipp 00394 00395 failure = .FALSE. 00396 CALL timeset(routineN,handle) 00397 IF (.NOT. failure) THEN 00398 gcut2 = gcut * gcut 00399 s_dim = rho_tot_g % pw_grid % first_gne0 00400 e_dim = rho_tot_g % pw_grid % ngpts_cut_local 00401 igmax = 0 00402 DO ig = s_dim, e_dim 00403 g2 = rho_tot_g % pw_grid % gsq ( ig ) 00404 IF (g2 > gcut2) EXIT 00405 igmax = ig 00406 ENDDO 00407 ALLOCATE(arg1_v_sin(s_dim:igmax),stat=stat) 00408 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00409 00410 Nr=SIZE(radii) 00411 Np=SIZE(particle_set) 00412 Ng=igmax-s_dim+1 00413 ALLOCATE(lhs(nparticles*Nr,Ng), stat=stat) 00414 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00415 ALLOCATE(rhs(Ng,Np*Nr), stat=stat) 00416 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00417 00418 ! rhs with first term of sin(g.(rvec1-rvec2)), which used to be called arg1_v_sin() 00419 ! rhs has all parts that depend on iparticle2 00420 DO iparticle2=1,Np 00421 DO igauss2=1,Nr 00422 rhs(1:Ng,(iparticle2-1)*Nr+igauss2) = g_dot_rvec_sin(1:Ng, iparticle2)*gfunc(s_dim:igmax,igauss2) 00423 END DO 00424 END DO 00425 DO icomp=1,3 00426 ! create lhs, which has all parts that depend on iparticle1 00427 DO ipp=1, nparticles 00428 iparticle1 = iparticle0 + ipp - 1 00429 DO ig=s_dim,igmax 00430 lhs((ipp-1)*Nr+1:(ipp-1)*Nr+Nr,ig-s_dim+1) = w(ig) * rho_tot_g%pw_grid%g(icomp,ig) * & 00431 gfunc(ig,1:Nr) * g_dot_rvec_cos(ig-s_dim+1, iparticle1) 00432 END DO 00433 END DO ! ipp 00434 ! do main multiply 00435 CALL DGEMM('N', 'N', nparticles*Nr, Np*Nr, Ng, 1.0D0, lhs(1,1), nparticles*Nr, rhs(1,1), & 00436 Ng, 0.0D0, dAm((iparticle0-1)*Nr+1,1,icomp), Np*Nr) 00437 ! do extra multiplies to compensate for missing factor of 2 00438 DO ipp=1, nparticles 00439 iparticle1 = iparticle0 + ipp - 1 00440 CALL DGEMM('N', 'N', Nr, Nr, Ng, 1.0D0, lhs((ipp-1)*Nr+1,1), nparticles*Nr, rhs(1,(iparticle1-1)*Nr+1), & 00441 Ng, 1.0D0, dAm((iparticle1-1)*Nr+1,(iparticle1-1)*Nr+1,icomp), Np*Nr) 00442 END DO 00443 ! now extra columns to account for factor of 2 in some rhs columns 00444 END DO ! icomp 00445 00446 ! rhs with second term of sin(g.(rvec1-rvec2)), which used to be called arg1_v_sin() 00447 ! rhs has all parts that depend on iparticle2 00448 DO iparticle2=1,Np 00449 DO igauss2=1,Nr 00450 rhs(1:Ng,(iparticle2-1)*Nr+igauss2) = -g_dot_rvec_cos(1:Ng, iparticle2)*gfunc(s_dim:igmax,igauss2) 00451 END DO 00452 END DO 00453 DO icomp=1,3 00454 ! create lhs, which has all parts that depend on iparticle1 00455 DO ipp=1, nparticles 00456 iparticle1 = iparticle0 + ipp - 1 00457 DO ig=s_dim,igmax 00458 lhs((ipp-1)*Nr+1:(ipp-1)*Nr+Nr,ig-s_dim+1) = w(ig) * rho_tot_g%pw_grid%g(icomp,ig) * gfunc(ig,1:Nr) * & 00459 g_dot_rvec_sin(ig-s_dim+1, iparticle1) 00460 END DO 00461 ENDDO 00462 ! do main multiply 00463 CALL DGEMM('N', 'N', nparticles*Nr, Np*Nr, Ng, 1.0D0, lhs(1,1), nparticles*Nr, rhs(1,1), & 00464 Ng, 1.0D0, dAm((iparticle0-1)*Nr+1,1,icomp), Np*Nr) 00465 ! do extra multiples to compensate for missing factor of 2 00466 DO ipp=1, nparticles 00467 iparticle1 = iparticle0 + ipp - 1 00468 CALL DGEMM('N', 'N', Nr, Nr, Ng, 1.0D0, lhs((ipp-1)*Nr+1,1), nparticles*Nr, rhs(1,(iparticle1-1)*Nr+1), & 00469 Ng, 1.0D0, dAm((iparticle1-1)*Nr+1,(iparticle1-1)*Nr+1,icomp), Np*Nr) 00470 END DO 00471 END DO 00472 00473 DEALLOCATE(rhs,stat=stat) 00474 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00475 DEALLOCATE(lhs,stat=stat) 00476 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00477 !!! DEALLOCATE(arg1_v,stat=stat) 00478 !!! CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00479 DEALLOCATE(arg1_v_sin,stat=stat) 00480 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00481 END IF 00482 CALL timestop(handle) 00483 END SUBROUTINE build_der_A_matrix_rows 00484 00485 !NB deallocate g_dot_rvec_* arrays 00486 SUBROUTINE cleanup_g_dot_rvec_sin_cos(g_dot_rvec_sin, g_dot_rvec_cos, error) 00487 REAL(KIND=dp), ALLOCATABLE, 00488 DIMENSION(:, :) :: g_dot_rvec_sin, g_dot_rvec_cos 00489 TYPE(cp_error_type), INTENT(inout) :: error 00490 00491 CHARACTER(len=*), PARAMETER :: routineN = 'cleanup_g_dot_rvec_sin_cos', 00492 routineP = moduleN//':'//routineN 00493 00494 INTEGER :: stat 00495 LOGICAL :: failure 00496 00497 DEALLOCATE(g_dot_rvec_sin,stat=stat) 00498 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00499 DEALLOCATE(g_dot_rvec_cos,stat=stat) 00500 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00501 END SUBROUTINE cleanup_g_dot_rvec_sin_cos 00502 00503 !NB precompute sin(g.r) and cos(g.r) for quicker evaluations of sin(g.(r1-r2)) and cos(g.(r1-r2)) 00504 SUBROUTINE prep_g_dot_rvec_sin_cos(rho_tot_g, particle_set, gcut, g_dot_rvec_sin, g_dot_rvec_cos, error) 00505 TYPE(pw_type), POINTER :: rho_tot_g 00506 TYPE(particle_type), DIMENSION(:), 00507 POINTER :: particle_set 00508 REAL(KIND=dp), INTENT(IN) :: gcut 00509 REAL(KIND=dp), ALLOCATABLE, 00510 DIMENSION(:, :) :: g_dot_rvec_sin, g_dot_rvec_cos 00511 TYPE(cp_error_type), INTENT(inout) :: error 00512 00513 CHARACTER(len=*), PARAMETER :: routineN = 'prep_g_dot_rvec_sin_cos', 00514 routineP = moduleN//':'//routineN 00515 00516 INTEGER :: e_dim, ig, igmax, iparticle, 00517 s_dim, stat 00518 LOGICAL :: failure 00519 REAL(KIND=dp) :: g2, g_dot_rvec, gcut2, rvec(3) 00520 00521 gcut2 = gcut * gcut 00522 s_dim = rho_tot_g % pw_grid % first_gne0 00523 e_dim = rho_tot_g % pw_grid % ngpts_cut_local 00524 igmax = 0 00525 DO ig = s_dim, e_dim 00526 g2 = rho_tot_g % pw_grid % gsq ( ig ) 00527 IF (g2 > gcut2) EXIT 00528 igmax = ig 00529 ENDDO 00530 00531 ALLOCATE(g_dot_rvec_sin(1:igmax-s_dim+1,SIZE(particle_set)),stat=stat) 00532 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00533 ALLOCATE(g_dot_rvec_cos(1:igmax-s_dim+1,SIZE(particle_set)),stat=stat) 00534 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00535 00536 DO iparticle=1, SIZE(particle_set) 00537 rvec = particle_set(iparticle)%r 00538 DO ig=s_dim, igmax 00539 g_dot_rvec = DOT_PRODUCT(rho_tot_g%pw_grid%g(:,ig),rvec) 00540 g_dot_rvec_sin(ig-s_dim+1,iparticle) = SIN(g_dot_rvec) 00541 g_dot_rvec_cos(ig-s_dim+1,iparticle) = COS(g_dot_rvec) 00542 END DO 00543 END DO 00544 00545 END SUBROUTINE prep_g_dot_rvec_sin_cos 00546 00547 ! ***************************************************************************** 00553 SUBROUTINE ddapc_eval_AmI(GAmI, c0, gfunc, w, particle_set, gcut,& 00554 rho_tot_g, radii, iw, Vol, error) 00555 REAL(KIND=dp), DIMENSION(:, :), POINTER :: GAmI 00556 REAL(KIND=dp), INTENT(OUT) :: c0 00557 REAL(KIND=dp), DIMENSION(:, :), POINTER :: gfunc 00558 REAL(KIND=dp), DIMENSION(:), POINTER :: w 00559 TYPE(particle_type), DIMENSION(:), 00560 POINTER :: particle_set 00561 REAL(KIND=dp), INTENT(IN) :: gcut 00562 TYPE(pw_type), POINTER :: rho_tot_g 00563 REAL(KIND=dp), DIMENSION(:), POINTER :: radii 00564 INTEGER, INTENT(IN) :: iw 00565 REAL(KIND=dp), INTENT(IN) :: Vol 00566 TYPE(cp_error_type), INTENT(inout) :: error 00567 00568 CHARACTER(len=*), PARAMETER :: routineN = 'ddapc_eval_AmI', 00569 routineP = moduleN//':'//routineN 00570 00571 INTEGER :: handle, ndim, stat 00572 LOGICAL :: failure 00573 REAL(KIND=dp) :: condition_number, inv_error 00574 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: AmE, cv 00575 REAL(KIND=dp), ALLOCATABLE, 00576 DIMENSION(:, :) :: Am, AmI, Amw, g_dot_rvec_cos, 00577 g_dot_rvec_sin 00578 00579 !NB for precomputation of sin(g.r) and cos(g.r) 00580 00581 failure = .FALSE. 00582 CALL timeset(routineN,handle) 00583 IF (.NOT. failure) THEN 00584 ndim = SIZE(particle_set)*SIZE(radii) 00585 ALLOCATE(Am(ndim, ndim), stat=stat) 00586 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00587 ALLOCATE(AmI(ndim, ndim), stat=stat) 00588 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00589 ALLOCATE(GAmI(ndim, ndim), stat=stat) 00590 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00591 ALLOCATE(cv(ndim), stat=stat) 00592 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00593 Am = 0.0_dp 00594 AmI = 0.0_dp 00595 cv = 1.0_dp/Vol 00596 !NB precompute sin(g.r) and cos(g.r) for faster evaluation of cos(g.(r1-r2)) in build_A_matrix() 00597 CALL prep_g_dot_rvec_sin_cos(rho_tot_g, particle_set, gcut, g_dot_rvec_sin, g_dot_rvec_cos, error) 00598 CALL build_A_matrix(Am, gfunc, w, particle_set, radii, rho_tot_g, gcut, g_dot_rvec_sin, g_dot_rvec_cos, error) 00599 CALL cleanup_g_dot_rvec_sin_cos(g_dot_rvec_sin, g_dot_rvec_cos, error) 00600 Am = Am / (Vol*Vol) 00601 CALL mp_sum(Am,rho_tot_g%pw_grid%para%group) 00602 IF (iw>0) THEN 00603 ! Checking conditions numbers and eigenvalues 00604 ALLOCATE(Amw(ndim, ndim), stat=stat) 00605 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00606 ALLOCATE(AmE(ndim), stat=stat) 00607 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00608 Amw = Am 00609 CALL diamat_all(Amw, AmE,error=error) 00610 condition_number = MAXVAL(ABS(AmE))/MINVAL(ABS(AmE)) 00611 WRITE(iw,'(T3,A)')" Eigenvalues of Matrix A:" 00612 WRITE(iw,'(T3,4E15.8)') AmE 00613 WRITE(iw,'(T3,A,1E15.9)')" Condition number:",condition_number 00614 IF (condition_number > 1.0E12_dp) THEN 00615 WRITE (iw,FMT="(/,T2,A)")& 00616 "WARNING: high condition number => possibly ill-conditioned matrix" 00617 END IF 00618 DEALLOCATE(Amw, stat=stat) 00619 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00620 DEALLOCATE(AmE, stat=stat) 00621 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00622 END IF 00623 CALL invert_matrix(Am, AmI, inv_error, "N",error=error, improve=.FALSE.) 00624 IF (iw>0) THEN 00625 WRITE(iw,'(T3,A,F15.9)')" Error inverting the A matrix: ", inv_error 00626 END IF 00627 c0 = DOT_PRODUCT(cv,MATMUL(AmI,cv)) 00628 DEALLOCATE(Am, stat=stat) 00629 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00630 DEALLOCATE(cv, stat=stat) 00631 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00632 GAmI = AmI 00633 DEALLOCATE(AmI, stat=stat) 00634 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00635 END IF 00636 CALL timestop(handle) 00637 END SUBROUTINE ddapc_eval_AmI 00638 00639 ! ***************************************************************************** 00647 RECURSIVE SUBROUTINE ewald_ddapc_pot(cp_para_env, coeff, factor, cell, multipole_section,& 00648 particle_set, M, radii, error) 00649 TYPE(cp_para_env_type), POINTER :: cp_para_env 00650 TYPE(pw_type), POINTER :: coeff 00651 REAL(KIND=dp), INTENT(IN) :: factor 00652 TYPE(cell_type), POINTER :: cell 00653 TYPE(section_vals_type), POINTER :: multipole_section 00654 TYPE(particle_type), DIMENSION(:), 00655 POINTER :: particle_set 00656 REAL(KIND=dp), DIMENSION(:, :), POINTER :: M 00657 REAL(KIND=dp), DIMENSION(:), POINTER :: radii 00658 TYPE(cp_error_type), INTENT(inout) :: error 00659 00660 CHARACTER(len=*), PARAMETER :: routineN = 'ewald_ddapc_pot', 00661 routineP = moduleN//':'//routineN 00662 00663 INTEGER :: ewmdim, handle, idim, idim1, idim2, idimo, igauss1, igauss2, 00664 ip1, ip2, iparticle1, iparticle2, istart_g, k1, k2, k3, n_rep, ndim, 00665 nmax1, nmax2, nmax3, r1, r2, r3, rmax1, rmax2, rmax3, sfact, stat 00666 LOGICAL :: analyt, failure 00667 REAL(KIND=dp) :: alpha, eps, ew_neut, fac, fac3, fs, g_ewald, galpha, 00668 gsq, gsqi, ij_fac, my_val, r, r2tmp, r_ewald, rc1, rc12, rc2, rc22, 00669 rcut, rcut2, t1, tol, tol1 00670 REAL(KIND=dp), DIMENSION(3) :: fvec, gvec, ra, rvec 00671 REAL(KIND=dp), DIMENSION(:), POINTER :: EwM 00672 00673 failure = .FALSE. 00674 NULLIFY(EwM) 00675 CALL timeset(routineN,handle) 00676 IF (.NOT.failure) THEN 00677 CPPostcondition(.NOT.ASSOCIATED(M),cp_failure_level,routineP,error,failure) 00678 CPPostcondition(ASSOCIATED(radii),cp_failure_level,routineP,error,failure) 00679 CPPostcondition(cell%orthorhombic,cp_failure_level,routineP,error,failure) 00680 sfact = factor 00681 rcut = MIN(cell%hmat(1,1),cell%hmat(2,2),cell%hmat(3,3))/2.0_dp 00682 CALL section_vals_val_get(multipole_section,"RCUT",n_rep_val=n_rep,error=error) 00683 IF (n_rep==1) CALL section_vals_val_get(multipole_section,"RCUT",r_val=rcut,error=error) 00684 CALL section_vals_val_get(multipole_section,"EWALD_PRECISION",r_val=eps,error=error) 00685 CALL section_vals_val_get(multipole_section,"ANALYTICAL_GTERM",l_val=analyt,error=error) 00686 rcut2 = rcut**2 00687 ! 00688 ! Setting-up parameters for Ewald summation 00689 ! 00690 eps = MIN(ABS(eps),0.5_dp) 00691 tol = SQRT(ABS(LOG(eps*rcut))) 00692 alpha = SQRT(ABS(LOG(eps*rcut*tol)))/rcut 00693 galpha = 1.0_dp/(4.0_dp*alpha*alpha) 00694 tol1 = SQRT(-LOG(eps*rcut*(2.0_dp*tol*alpha)**2)) 00695 nmax1 = NINT(0.25_dp + cell%hmat(1,1)*alpha*tol1/pi) 00696 nmax2 = NINT(0.25_dp + cell%hmat(2,2)*alpha*tol1/pi) 00697 nmax3 = NINT(0.25_dp + cell%hmat(3,3)*alpha*tol1/pi) 00698 00699 rmax1 = CEILING(rcut/cell%hmat(1,1)) 00700 rmax2 = CEILING(rcut/cell%hmat(2,2)) 00701 rmax3 = CEILING(rcut/cell%hmat(3,3)) 00702 fac = 1.e0_dp/cell%deth 00703 fac3 = fac*pi 00704 fvec = twopi / (/cell%hmat(1,1),cell%hmat(2,2),cell%hmat(3,3)/) 00705 ew_neut = - fac * pi / alpha ** 2 00706 ! 00707 ewmdim = SIZE(particle_set) * (SIZE(particle_set)+1) / 2 00708 ndim = SIZE(particle_set) * SIZE(radii) 00709 ALLOCATE(EwM(ewmdim), stat=stat) 00710 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00711 ALLOCATE(M(ndim, ndim), stat=stat) 00712 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00713 M = 0.0_dp 00714 ! 00715 idim = 0 00716 !NB zero EwM so mp_sum(EwM) will give right answer 00717 EwM = 0.0_dp 00718 DO iparticle1 = 1, SIZE(particle_set) 00719 ip1 = (iparticle1-1)*SIZE(radii) 00720 DO iparticle2 = 1, iparticle1 00721 ij_fac = 1.0_dp 00722 IF (iparticle1==iparticle2) ij_fac=0.5_dp 00723 00724 ip2 = (iparticle2-1)*SIZE(radii) 00725 idim = idim + 1 00726 !NB parallelization, done here so indexing is right 00727 IF (MOD(iparticle1, cp_para_env%num_pe) /= cp_para_env%mepos) CYCLE 00728 ! 00729 ! Real-Space Contribution 00730 ! 00731 my_val = 0.0_dp 00732 rvec = particle_set(iparticle1)%r - particle_set(iparticle2)%r 00733 r_ewald = 0.0_dp 00734 IF (iparticle1 /= iparticle2) THEN 00735 ra = rvec 00736 r2tmp = DOT_PRODUCT(ra,ra) 00737 IF (r2tmp<=rcut2) THEN 00738 r = SQRT(r2tmp) 00739 t1 = erfc(alpha*r) / r 00740 r_ewald = t1 00741 END IF 00742 END IF 00743 DO r1 = -rmax1, rmax1 00744 DO r2 = -rmax2, rmax2 00745 DO r3 = -rmax3, rmax3 00746 IF ((r1==0).AND.(r2==0).AND.(r3==0)) CYCLE 00747 ra(1) = rvec(1) + cell%hmat(1,1)*r1 00748 ra(2) = rvec(2) + cell%hmat(2,2)*r2 00749 ra(3) = rvec(3) + cell%hmat(3,3)*r3 00750 r2tmp = DOT_PRODUCT(ra,ra) 00751 IF (r2tmp<=rcut2) THEN 00752 r = SQRT(r2tmp) 00753 t1 = erfc(alpha*r) / r 00754 r_ewald = r_ewald + t1*ij_fac 00755 END IF 00756 END DO 00757 END DO 00758 END DO 00759 ! 00760 ! G-space Contribution 00761 ! 00762 IF (analyt) THEN 00763 g_ewald = 0.0_dp 00764 DO k1 = 0, nmax1 00765 DO k2 = -nmax2, nmax2 00766 DO k3 = -nmax3, nmax3 00767 IF (k1 == 0.AND.k2 == 0.AND.k3 == 0) CYCLE 00768 fs = 2.0_dp; IF (k1==0) fs = 1.0_dp 00769 gvec = fvec * (/REAL(k1,KIND=dp),REAL(k2,KIND=dp),REAL(k3,KIND=dp)/) 00770 gsq = DOT_PRODUCT(gvec,gvec) 00771 gsqi = fs/gsq 00772 t1 = fac * gsqi * EXP(-galpha*gsq) 00773 g_ewald = g_ewald + t1 * COS(DOT_PRODUCT(gvec,rvec)) 00774 END DO 00775 END DO 00776 END DO 00777 ELSE 00778 g_ewald = Eval_Interp_Spl3_pbc(rvec, coeff, error) 00779 END IF 00780 ! 00781 ! G-EWALD, R-EWALD 00782 ! 00783 g_ewald = r_ewald + fourpi*g_ewald 00784 ! 00785 ! Self Contribution 00786 ! 00787 IF (iparticle1 == iparticle2) THEN 00788 g_ewald = g_ewald - 2.0_dp * alpha * oorootpi 00789 END IF 00790 ! 00791 IF (iparticle1/=iparticle2) THEN 00792 ra = rvec 00793 r = SQRT(DOT_PRODUCT(ra,ra)) 00794 my_val = sfact / r 00795 END IF 00796 EwM(idim) = my_val - sfact * g_ewald 00797 END DO ! iparticle2 00798 END DO ! iparticle1 00799 !NB sum over parallelized contributions of different nodes 00800 CALL mp_sum(EwM, cp_para_env%group) 00801 idim = 0 00802 DO iparticle2 = 1, SIZE(particle_set) 00803 ip2 = (iparticle2-1)*SIZE(radii) 00804 idimo = (iparticle2-1) 00805 idimo = idimo *(idimo+1) / 2 00806 DO igauss2 = 1, SIZE(radii) 00807 idim2 = ip2 + igauss2 00808 rc2 = radii(igauss2) 00809 rc22 = rc2*rc2 00810 DO iparticle1 = 1, iparticle2 00811 ip1 = (iparticle1-1)*SIZE(radii) 00812 idim = idimo + iparticle1 00813 istart_g = 1 00814 IF (iparticle1==iparticle2) istart_g = igauss2 00815 DO igauss1 = istart_g, SIZE(radii) 00816 idim1 = ip1 + igauss1 00817 rc1 = radii(igauss1) 00818 rc12 = rc1*rc1 00819 M(idim1,idim2) = EwM(idim) - sfact * ew_neut - sfact * fac3 * (rc12+rc22) 00820 M(idim2,idim1) = M(idim1,idim2) 00821 END DO 00822 END DO 00823 END DO ! iparticle2 00824 END DO ! iparticle1 00825 DEALLOCATE(EwM, stat=stat) 00826 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00827 END IF 00828 CALL timestop(handle) 00829 END SUBROUTINE ewald_ddapc_pot 00830 00831 ! ***************************************************************************** 00838 SUBROUTINE solvation_ddapc_pot( solvation_section, particle_set, M, radii, error) 00839 TYPE(section_vals_type), POINTER :: solvation_section 00840 TYPE(particle_type), DIMENSION(:), 00841 POINTER :: particle_set 00842 REAL(KIND=dp), DIMENSION(:, :), POINTER :: M 00843 REAL(KIND=dp), DIMENSION(:), POINTER :: radii 00844 TYPE(cp_error_type), INTENT(inout) :: error 00845 00846 CHARACTER(len=*), PARAMETER :: routineN = 'solvation_ddapc_pot', 00847 routineP = moduleN//':'//routineN 00848 00849 INTEGER :: i, idim, idim1, idim2, igauss1, igauss2, ip1, ip2, iparticle1, 00850 iparticle2, istart_g, j, l, lmax, n_rep1, n_rep2, ndim, output_unit, 00851 stat, weight 00852 INTEGER, DIMENSION(:), POINTER :: list 00853 LOGICAL :: failure, fixed_center 00854 REAL(KIND=dp) :: center(3), eps_in, eps_out, 00855 factor, mass, mycos, r1, r2, 00856 Rs, rvec(3) 00857 REAL(KIND=dp), DIMENSION(:), POINTER :: pos, R0 00858 REAL(KIND=dp), DIMENSION(:, :), POINTER :: cost, LocP 00859 TYPE(cp_logger_type), POINTER :: logger 00860 00861 failure = .FALSE. 00862 fixed_center = .FALSE. 00863 NULLIFY(logger) 00864 logger => cp_error_get_logger(error) 00865 output_unit= cp_logger_get_default_io_unit(logger) 00866 IF (.NOT.failure) THEN 00867 ndim = SIZE(particle_set) * SIZE(radii) 00868 ALLOCATE(M(ndim, ndim), stat=stat) 00869 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00870 M = 0.0_dp 00871 eps_in = 1.0_dp 00872 CALL section_vals_val_get(solvation_section,"EPS_OUT",r_val=eps_out,error=error) 00873 CALL section_vals_val_get(solvation_section,"LMAX",i_val=lmax,error=error) 00874 CALL section_vals_val_get(solvation_section,"SPHERE%RADIUS",r_val=Rs,error=error) 00875 CALL section_vals_val_get(solvation_section,"SPHERE%CENTER%XYZ",n_rep_val=n_rep1,& 00876 error=error) 00877 IF (n_rep1/=0) THEN 00878 CALL section_vals_val_get(solvation_section,"SPHERE%CENTER%XYZ",r_vals=R0,& 00879 error=error) 00880 center = R0 00881 ELSE 00882 CALL section_vals_val_get(solvation_section,"SPHERE%CENTER%ATOM_LIST",& 00883 n_rep_val=n_rep2,error=error) 00884 IF (n_rep2/=0) THEN 00885 CALL section_vals_val_get(solvation_section,"SPHERE%CENTER%ATOM_LIST",i_vals=list,& 00886 error=error) 00887 CALL section_vals_val_get(solvation_section,"SPHERE%CENTER%WEIGHT_TYPE",i_val=weight,& 00888 error=error) 00889 ALLOCATE(R0(SIZE(list)),stat=stat) 00890 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00891 SELECT CASE(weight) 00892 CASE (weight_type_unit) 00893 R0 = 0.0_dp 00894 DO i = 1, SIZE(list) 00895 R0 = R0 + particle_set(list(i))%r 00896 END DO 00897 R0 = R0 / REAL(SIZE(list),KIND=dp) 00898 CASE (weight_type_mass) 00899 R0 = 0.0_dp 00900 mass = 0.0_dp 00901 DO i = 1, SIZE(list) 00902 R0 = R0 + particle_set(list(i))%r * particle_set(list(i))%atomic_kind%mass 00903 mass = mass + particle_set(list(i))%atomic_kind%mass 00904 END DO 00905 R0 = R0 / mass 00906 END SELECT 00907 center = R0 00908 CALL section_vals_val_get(solvation_section,"SPHERE%CENTER%FIXED",l_val=fixed_center,& 00909 error=error) 00910 IF (fixed_center) THEN 00911 CALL section_vals_val_set(solvation_section,"SPHERE%CENTER%XYZ",& 00912 r_vals_ptr=R0,error=error) 00913 ELSE 00914 DEALLOCATE(R0, stat=stat) 00915 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00916 END IF 00917 END IF 00918 END IF 00919 CPPostcondition(n_rep1/=0.OR.n_rep2/=0,cp_failure_level,routineP,error,failure) 00920 ! Potential calculation 00921 ALLOCATE(LocP(0:lmax,SIZE(particle_set)),stat=stat) 00922 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00923 ALLOCATE(pos(SIZE(particle_set)),stat=stat) 00924 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00925 ALLOCATE(cost(SIZE(particle_set),SIZE(particle_set)),stat=stat) 00926 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00927 ! Determining the single atomic contribution to the dielectric dipole 00928 DO i = 1, SIZE(particle_set) 00929 rvec = particle_set(i)%r-center 00930 r2 = DOT_PRODUCT(rvec,rvec) 00931 r1 = SQRT(r2) 00932 IF (r1>=Rs) THEN 00933 IF (output_unit>0) THEN 00934 WRITE(output_unit,'(A,I6,A)')"Atom number :: ",i," is out of the solvation sphere" 00935 WRITE(output_unit,'(2(A,F12.6))')"Distance from the center::",r1," Radius of the sphere::",rs 00936 END IF 00937 CPPostcondition(.FALSE.,cp_failure_level,routineP,error,failure) 00938 END IF 00939 LocP(:,i) = 0.0_dp 00940 IF (r1/=0.0_dp) THEN 00941 DO l = 0, lmax 00942 LocP(l,i) = ( r1**l * REAL(l+1,KIND=dp) * (eps_in-eps_out))/ 00943 (Rs**(2*l+1)*eps_in*(REAL(l,KIND=dp)*eps_in+REAL(l+1,KIND=dp)*eps_out)) 00944 END DO 00945 ELSE 00946 ! limit for r->0 00947 LocP(0,i) = (eps_in-eps_out)/(Rs*eps_in*eps_out) 00948 END IF 00949 pos(i) = r1 00950 END DO 00951 ! Particle-Particle potential energy matrix 00952 cost = 0.0_dp 00953 DO i = 1, SIZE(particle_set) 00954 DO j = 1, i 00955 factor = 0.0_dp 00956 IF (pos(i)*pos(j)/=0.0_dp) THEN 00957 mycos = DOT_PRODUCT(particle_set(i)%r-center,particle_set(j)%r-center)/(pos(i)*pos(j)) 00958 IF (ABS(mycos) > 1.0_dp) mycos = SIGN(1.0_dp, mycos) 00959 DO l = 0, lmax 00960 factor = factor + LocP(l,i) * pos(j)**l * legendre(mycos,l,0) 00961 END DO 00962 ELSE 00963 factor = LocP(0,i) 00964 END IF 00965 cost(i,j) = factor 00966 cost(j,i) = factor 00967 END DO 00968 END DO 00969 ! Computes the full potential energy matrix 00970 idim = 0 00971 DO iparticle2 = 1, SIZE(particle_set) 00972 ip2 = (iparticle2-1)*SIZE(radii) 00973 DO igauss2 = 1, SIZE(radii) 00974 idim2 = ip2 + igauss2 00975 DO iparticle1 = 1, iparticle2 00976 ip1 = (iparticle1-1)*SIZE(radii) 00977 istart_g = 1 00978 IF (iparticle1==iparticle2) istart_g = igauss2 00979 DO igauss1 = istart_g, SIZE(radii) 00980 idim1 = ip1 + igauss1 00981 M(idim1,idim2) = cost(iparticle1,iparticle2) 00982 M(idim2,idim1) = M(idim1,idim2) 00983 END DO 00984 END DO 00985 END DO 00986 END DO 00987 DEALLOCATE(cost,stat=stat) 00988 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00989 DEALLOCATE(pos,stat=stat) 00990 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00991 DEALLOCATE(LocP,stat=stat) 00992 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00993 END IF 00994 END SUBROUTINE solvation_ddapc_pot 00995 00996 END MODULE cp_ddapc_methods
1.7.3