|
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 ! ***************************************************************************** 00018 MODULE semi_empirical_int_gks 00019 00020 USE dg_rho0_types, ONLY: dg_rho0_type 00021 USE dg_types, ONLY: dg_get,& 00022 dg_type 00023 USE erf_fn, ONLY: erfc 00024 USE f77_blas 00025 USE input_constants, ONLY: do_multipole_none 00026 USE kinds, ONLY: dp 00027 USE mathconstants, ONLY: fourpi,& 00028 oorootpi 00029 USE pw_grid_types, ONLY: pw_grid_type 00030 USE pw_pool_types, ONLY: pw_pool_type 00031 USE semi_empirical_int_arrays, ONLY: indexb,& 00032 rij_threshold 00033 USE semi_empirical_mpole_types, ONLY: semi_empirical_mpole_type 00034 USE semi_empirical_types, ONLY: se_int_control_type,& 00035 se_taper_type,& 00036 semi_empirical_type,& 00037 setup_se_int_control_type 00038 #include "cp_common_uses.h" 00039 00040 IMPLICIT NONE 00041 PRIVATE 00042 00043 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'semi_empirical_int_gks' 00044 LOGICAL, PARAMETER, PRIVATE :: debug_this_module=.FALSE. 00045 00046 PUBLIC :: corecore_gks, rotnuc_gks, drotnuc_gks, rotint_gks, drotint_gks 00047 00048 CONTAINS 00049 00050 ! ***************************************************************************** 00053 SUBROUTINE rotnuc_gks (sepi,sepj,rij,e1b,e2a,se_int_control,se_taper,error) 00054 TYPE(semi_empirical_type), POINTER :: sepi, sepj 00055 REAL(dp), DIMENSION(3), INTENT(IN) :: rij 00056 REAL(dp), DIMENSION(45), INTENT(OUT), 00057 OPTIONAL :: e1b, e2a 00058 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 00059 TYPE(se_taper_type), POINTER :: se_taper 00060 TYPE(cp_error_type), INTENT(inout) :: error 00061 00062 CHARACTER(len=*), PARAMETER :: routineN = 'rotnuc_gks', 00063 routineP = moduleN//':'//routineN 00064 00065 INTEGER :: i, mu, nu 00066 REAL(KIND=dp), DIMENSION(3) :: rab 00067 REAL(kind=dp), DIMENSION(45, 45) :: Coul 00068 00069 rab= -rij 00070 00071 IF(se_int_control%do_ewald_gks) THEN 00072 IF ( DOT_PRODUCT(rij,rij) > rij_threshold) THEN 00073 CALL makeCoulE(rab,sepi,sepj,Coul,se_int_control,error) 00074 ELSE 00075 CALL makeCoulE0(sepi,Coul,se_int_control,error) 00076 END IF 00077 ELSE 00078 CALL makeCoul(rab,sepi,sepj,Coul,se_int_control,error) 00079 END IF 00080 00081 i = 0 00082 DO mu = 1, sepi%natorb 00083 DO nu = 1, mu 00084 i = i + 1 00085 e1b(i)= -Coul(i,1)*sepj%zeff 00086 END DO 00087 END DO 00088 00089 i = 0 00090 DO mu = 1, sepj%natorb 00091 DO nu = 1, mu 00092 i = i + 1 00093 e2a(i)= -Coul(1,i)*sepi%zeff 00094 END DO 00095 END DO 00096 00097 END SUBROUTINE rotnuc_gks 00098 00099 ! ***************************************************************************** 00102 SUBROUTINE rotint_gks (sepi,sepj,rij,w,se_int_control,se_taper,error) 00103 TYPE(semi_empirical_type), POINTER :: sepi, sepj 00104 REAL(dp), DIMENSION(3), INTENT(IN) :: rij 00105 REAL(dp), DIMENSION(2025), INTENT(OUT), 00106 OPTIONAL :: w 00107 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 00108 TYPE(se_taper_type), POINTER :: se_taper 00109 TYPE(cp_error_type), INTENT(inout) :: error 00110 00111 CHARACTER(len=*), PARAMETER :: routineN = 'rotint_gks', 00112 routineP = moduleN//':'//routineN 00113 00114 INTEGER :: i, ind1, ind2, lam, mu, nu, 00115 sig 00116 REAL(KIND=dp), DIMENSION(3) :: rab 00117 REAL(kind=dp), DIMENSION(45, 45) :: Coul 00118 00119 rab= -rij 00120 00121 IF(se_int_control%do_ewald_gks) THEN 00122 IF ( DOT_PRODUCT(rij,rij) > rij_threshold) THEN 00123 CALL makeCoulE(rab,sepi,sepj,Coul,se_int_control,error) 00124 ELSE 00125 CALL makeCoulE0(sepi,Coul,se_int_control,error) 00126 END IF 00127 ELSE 00128 CALL makeCoul(rab,sepi,sepj,Coul,se_int_control,error) 00129 END IF 00130 00131 i = 0 00132 ind1 = 0 00133 DO mu = 1, sepi%natorb 00134 DO nu = 1, mu 00135 ind1 = ind1 + 1 00136 ind2 = 0 00137 DO lam = 1, sepj%natorb 00138 DO sig = 1, lam 00139 i = i + 1 00140 ind2 = ind2 + 1 00141 w(i) = Coul(ind1,ind2) 00142 END DO 00143 END DO 00144 END DO 00145 END DO 00146 00147 END SUBROUTINE rotint_gks 00148 00149 ! ***************************************************************************** 00152 SUBROUTINE drotnuc_gks(sepi,sepj,rij,de1b,de2a,se_int_control,se_taper,error) 00153 TYPE(semi_empirical_type), POINTER :: sepi, sepj 00154 REAL(dp), DIMENSION(3), INTENT(IN) :: rij 00155 REAL(dp), DIMENSION(3, 45), 00156 INTENT(OUT), OPTIONAL :: de1b, de2a 00157 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 00158 TYPE(se_taper_type), POINTER :: se_taper 00159 TYPE(cp_error_type), INTENT(inout) :: error 00160 00161 CHARACTER(len=*), PARAMETER :: routineN = 'drotnuc_gks', 00162 routineP = moduleN//':'//routineN 00163 00164 INTEGER :: i, mu, nu 00165 REAL(KIND=dp), DIMENSION(3) :: rab 00166 REAL(kind=dp), DIMENSION(3, 45, 45) :: dCoul 00167 00168 rab= -rij 00169 00170 IF(se_int_control%do_ewald_gks) THEN 00171 CALL makedCoulE(rab,sepi,sepj,dCoul,se_int_control,error) 00172 ELSE 00173 CALL makedCoul(rab,sepi,sepj,dCoul,se_int_control,error) 00174 END IF 00175 00176 i = 0 00177 DO mu = 1, sepi%natorb 00178 DO nu = 1, mu 00179 i = i + 1 00180 de1b(1,i)= dCoul(1,i,1)*sepj%zeff 00181 de1b(2,i)= dCoul(2,i,1)*sepj%zeff 00182 de1b(3,i)= dCoul(3,i,1)*sepj%zeff 00183 END DO 00184 END DO 00185 00186 i = 0 00187 DO mu = 1, sepj%natorb 00188 DO nu = 1, mu 00189 i = i + 1 00190 de2a(1,i)= dCoul(1,1,i)*sepi%zeff 00191 de2a(2,i)= dCoul(2,1,i)*sepi%zeff 00192 de2a(3,i)= dCoul(3,1,i)*sepi%zeff 00193 END DO 00194 END DO 00195 00196 END SUBROUTINE drotnuc_gks 00197 00198 ! ***************************************************************************** 00201 SUBROUTINE drotint_gks(sepi,sepj,rij,dw,se_int_control,se_taper, error) 00202 TYPE(semi_empirical_type), POINTER :: sepi, sepj 00203 REAL(dp), DIMENSION(3), INTENT(IN) :: rij 00204 REAL(dp), DIMENSION(3, 2025), 00205 INTENT(OUT) :: dw 00206 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 00207 TYPE(se_taper_type), POINTER :: se_taper 00208 TYPE(cp_error_type), INTENT(inout) :: error 00209 00210 CHARACTER(len=*), PARAMETER :: routineN = 'drotint_gks', 00211 routineP = moduleN//':'//routineN 00212 00213 INTEGER :: i, ind1, ind2, lam, mu, nu, 00214 sig 00215 REAL(KIND=dp), DIMENSION(3) :: rab 00216 REAL(kind=dp), DIMENSION(3, 45, 45) :: dCoul 00217 00218 rab= -rij 00219 00220 IF(se_int_control%do_ewald_gks) THEN 00221 CALL makedCoulE(rab,sepi,sepj,dCoul,se_int_control,error) 00222 ELSE 00223 CALL makedCoul(rab,sepi,sepj,dCoul,se_int_control,error) 00224 END IF 00225 00226 i = 0 00227 ind1 = 0 00228 DO mu = 1, sepi%natorb 00229 DO nu = 1, mu 00230 ind1 = ind1 + 1 00231 ind2 = 0 00232 DO lam = 1, sepj%natorb 00233 DO sig = 1, lam 00234 i = i + 1 00235 ind2 = ind2 + 1 00236 dw(1,i)= -dCoul(1,ind1,ind2) 00237 dw(2,i)= -dCoul(2,ind1,ind2) 00238 dw(3,i)= -dCoul(3,ind1,ind2) 00239 END DO 00240 END DO 00241 END DO 00242 END DO 00243 00244 END SUBROUTINE drotint_gks 00245 00246 ! ***************************************************************************** 00249 SUBROUTINE makeCoul(RAB,sepi,sepj,Coul,se_int_control,error) 00250 REAL(kind=dp), DIMENSION(3) :: RAB 00251 TYPE(semi_empirical_type), POINTER :: sepi, sepj 00252 REAL(kind=dp), DIMENSION(45, 45), 00253 INTENT(OUT) :: Coul 00254 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 00255 TYPE(cp_error_type), INTENT(inout) :: error 00256 00257 CHARACTER(len=*), PARAMETER :: routineN = 'makeCoul', 00258 routineP = moduleN//':'//routineN 00259 00260 INTEGER :: iA, iB, imA, imB, jA, jB, k1, 00261 k2, k3, k4 00262 LOGICAL :: shortrange 00263 REAL(kind=dp) :: a2, ACOULA, ACOULB, d1, 00264 d1f(3), d2, d2f(3,3), d3, 00265 d3f(3,3,3), d4, d4f(3,3,3,3), 00266 f, rr, w, w0, w1, w2, w3, w4, 00267 w5 00268 REAL(kind=dp), DIMENSION(3) :: v 00269 REAL(kind=dp), DIMENSION(3, 3, 45) :: M2A, M2B 00270 REAL(kind=dp), DIMENSION(3, 45) :: M1A, M1B 00271 REAL(kind=dp), DIMENSION(45) :: M0A, M0B 00272 00273 shortrange = se_int_control%shortrange 00274 CALL get_se_slater_multipole(sepi, M0A, M1A, M2A, ACOULA) 00275 CALL get_se_slater_multipole(sepj, M0B, M1B, M2B, ACOULB) 00276 00277 v(1)=RAB(1) 00278 v(2)=RAB(2) 00279 v(3)=RAB(3) 00280 rr=SQRT(DOT_PRODUCT(v,v)) 00281 00282 a2=0.5_dp*(1.0_dp/ACOULA+1.0_dp/ACOULB) 00283 w0= a2*rr 00284 w= EXP(-w0) 00285 w1= (1.0_dp+0.5_dp*w0) 00286 w2= (w1+0.5_dp*w0+0.5_dp*w0**2) 00287 w3= (w2+w0**3/6.0_dp) 00288 w4= (w3+w0**4/30.0_dp) 00289 w5= (w3+8.0_dp*w0**4/210.0_dp+w0**5/210.0_dp) 00290 00291 IF(shortrange)THEN 00292 f= (-w*w1)/rr 00293 d1= -1.0_dp*(-w*w2)/rr**3 00294 d2= 3.0_dp*(-w*w3)/rr**5 00295 d3= -15.0_dp*(-w*w4)/rr**7 00296 d4= 105.0_dp*(-w*w5)/rr**9 00297 ELSE 00298 f= (1.0_dp-w*w1)/rr 00299 d1= -1.0_dp*(1.0_dp-w*w2)/rr**3 00300 d2= 3.0_dp*(1.0_dp-w*w3)/rr**5 00301 d3= -15.0_dp*(1.0_dp-w*w4)/rr**7 00302 d4= 105.0_dp*(1.0_dp-w*w5)/rr**9 00303 ENDIF 00304 00305 CALL build_d_tensor_gks(d1f, d2f, d3f, d4f, v=v, d1=d1, d2=d2, d3=d3, d4=d4) 00306 00307 imA = 0 00308 DO iA = 1, sepi%natorb 00309 DO jA = 1, iA 00310 imA = imA + 1 00311 00312 imB = 0 00313 DO iB = 1, sepj%natorb 00314 DO jB = 1, iB 00315 imB = imB + 1 00316 00317 w= M0A(imA)*M0B(imB)*f 00318 DO k1=1,3 00319 w=w+( M1A(k1,imA)*M0B(imB)-M0A(imA)*M1B(k1,imB) )*d1f(k1) 00320 ENDDO 00321 DO k2=1,3 00322 DO k1=1,3 00323 w=w+( M2A(k1,k2,imA)*M0B(imB)-M1A(k1,imA)*M1B(k2,imB)+M0A(imA)*M2B(k1,k2,imB) )*d2f(k1,k2) 00324 ENDDO 00325 ENDDO 00326 DO k3=1,3 00327 DO k2=1,3 00328 DO k1=1,3 00329 w=w+( -M2A(k1,k2,imA)*M1B(k3,imB)+M1A(k1,imA)*M2B(k2,k3,imB) )*d3f(k1,k2,k3) 00330 ENDDO 00331 ENDDO 00332 ENDDO 00333 00334 DO k4=1,3 00335 DO k3=1,3 00336 DO k2=1,3 00337 DO k1=1,3 00338 w=w+ M2A(k1,k2,imA)*M2B(k3,k4,imB)*d4f(k1,k2,k3,k4) 00339 ENDDO 00340 ENDDO 00341 ENDDO 00342 ENDDO 00343 00344 Coul(imA,imB)=w 00345 ENDDO 00346 ENDDO 00347 ENDDO 00348 ENDDO 00349 00350 END SUBROUTINE makeCoul 00351 00352 ! ***************************************************************************** 00356 SUBROUTINE makedCoul(RAB,sepi,sepj,dCoul,se_int_control,error) 00357 REAL(kind=dp), DIMENSION(3) :: RAB 00358 TYPE(semi_empirical_type), POINTER :: sepi, sepj 00359 REAL(kind=dp), DIMENSION(3, 45, 45), 00360 INTENT(OUT) :: dCoul 00361 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 00362 TYPE(cp_error_type), INTENT(inout) :: error 00363 00364 CHARACTER(len=*), PARAMETER :: routineN = 'makedCoul', 00365 routineP = moduleN//':'//routineN 00366 00367 INTEGER :: iA, iB, imA, imB, jA, jB, k1, 00368 k2, k3, k4 00369 LOGICAL :: shortrange 00370 REAL(kind=dp) :: a2, ACOULA, ACOULB, d1, d1f(3), d2, d2f(3,3), d3, 00371 d3f(3,3,3), d4, d4f(3,3,3,3), d5, d5f(3,3,3,3,3), f, rr, tmp, w, w0, 00372 w1, w2, w3, w4, w5, w6 00373 REAL(kind=dp), DIMENSION(3) :: v, wv 00374 REAL(kind=dp), DIMENSION(3, 3, 45) :: M2A, M2B 00375 REAL(kind=dp), DIMENSION(3, 45) :: M1A, M1B 00376 REAL(kind=dp), DIMENSION(45) :: M0A, M0B 00377 00378 shortrange = se_int_control%shortrange 00379 CALL get_se_slater_multipole(sepi, M0A, M1A, M2A, ACOULA) 00380 CALL get_se_slater_multipole(sepj, M0B, M1B, M2B, ACOULB) 00381 00382 v(1)=RAB(1) 00383 v(2)=RAB(2) 00384 v(3)=RAB(3) 00385 rr=SQRT(DOT_PRODUCT(v,v)) 00386 00387 a2=0.5_dp*(1.0_dp/ACOULA+1.0_dp/ACOULB) 00388 w0= a2*rr 00389 w= EXP(-w0) 00390 w1= (1.0_dp+0.5_dp*w0) 00391 w2= (w1+0.5_dp*w0+0.5_dp*w0**2) 00392 w3= (w2+w0**3/6.0_dp) 00393 w4= (w3+w0**4/30.0_dp) 00394 w5= (w3+4.0_dp*w0**4/105.0_dp+w0**5/210.0_dp) 00395 w6= (w3+15.0_dp*w0**4/378.0_dp+2.0_dp*w0**5/315.0_dp+w0**6/1890.0_dp) 00396 00397 IF(shortrange)THEN 00398 f= (-w*w1)/rr 00399 d1= -1.0_dp*(-w*w2)/rr**3 00400 d2= 3.0_dp*(-w*w3)/rr**5 00401 d3= -15.0_dp*(-w*w4)/rr**7 00402 d4= 105.0_dp*(-w*w5)/rr**9 00403 d5= -945.0_dp*(-w*w6)/rr**11 00404 ELSE 00405 f= (1.0_dp-w*w1)/rr 00406 d1= -1.0_dp*(1.0_dp-w*w2)/rr**3 00407 d2= 3.0_dp*(1.0_dp-w*w3)/rr**5 00408 d3= -15.0_dp*(1.0_dp-w*w4)/rr**7 00409 d4= 105.0_dp*(1.0_dp-w*w5)/rr**9 00410 d5= -945.0_dp*(1.0_dp-w*w6)/rr**11 00411 ENDIF 00412 00413 CALL build_d_tensor_gks(d1f, d2f, d3f, d4f, d5f, v, d1, d2, d3, d4, d5) 00414 00415 imA = 0 00416 DO iA = 1, sepi%natorb 00417 DO jA = 1, iA 00418 imA = imA + 1 00419 00420 imB = 0 00421 DO iB = 1, sepj%natorb 00422 DO jB = 1, iB 00423 imB = imB + 1 00424 00425 tmp = M0A(imA)*M0B(imB) 00426 wv(1)= tmp*d1f(1) 00427 wv(2)= tmp*d1f(2) 00428 wv(3)= tmp*d1f(3) 00429 DO k1=1,3 00430 tmp = M1A(k1,imA)*M0B(imB)-M0A(imA)*M1B(k1,imB) 00431 wv(1)=wv(1)+tmp*d2f(1,k1) 00432 wv(2)=wv(2)+tmp*d2f(2,k1) 00433 wv(3)=wv(3)+tmp*d2f(3,k1) 00434 ENDDO 00435 DO k2=1,3 00436 DO k1=1,3 00437 tmp = M2A(k1,k2,imA)*M0B(imB)-M1A(k1,imA)*M1B(k2,imB)+M0A(imA)*M2B(k1,k2,imB) 00438 wv(1)=wv(1)+tmp*d3f(1,k1,k2) 00439 wv(2)=wv(2)+tmp*d3f(2,k1,k2) 00440 wv(3)=wv(3)+tmp*d3f(3,k1,k2) 00441 ENDDO 00442 ENDDO 00443 DO k3=1,3 00444 DO k2=1,3 00445 DO k1=1,3 00446 tmp = -M2A(k1,k2,imA)*M1B(k3,imB)+M1A(k1,imA)*M2B(k2,k3,imB) 00447 wv(1)=wv(1)+tmp*d4f(1,k1,k2,k3) 00448 wv(2)=wv(2)+tmp*d4f(2,k1,k2,k3) 00449 wv(3)=wv(3)+tmp*d4f(3,k1,k2,k3) 00450 ENDDO 00451 ENDDO 00452 ENDDO 00453 00454 DO k4=1,3 00455 DO k3=1,3 00456 DO k2=1,3 00457 DO k1=1,3 00458 tmp = M2A(k1,k2,imA)*M2B(k3,k4,imB) 00459 wv(1)=wv(1)+ tmp*d5f(1,k1,k2,k3,k4) 00460 wv(2)=wv(2)+ tmp*d5f(2,k1,k2,k3,k4) 00461 wv(3)=wv(3)+ tmp*d5f(3,k1,k2,k3,k4) 00462 ENDDO 00463 ENDDO 00464 ENDDO 00465 ENDDO 00466 00467 dCoul(1,imA,imB)=wv(1) 00468 dCoul(2,imA,imB)=wv(2) 00469 dCoul(3,imA,imB)=wv(3) 00470 ENDDO 00471 ENDDO 00472 ENDDO 00473 ENDDO 00474 00475 END SUBROUTINE makedCoul 00476 00477 ! ***************************************************************************** 00480 SUBROUTINE corecore_gks (sepi,sepj,rijv,enuc,denuc,se_int_control,se_taper,error) 00481 TYPE(semi_empirical_type), POINTER :: sepi, sepj 00482 REAL(dp), DIMENSION(3), INTENT(IN) :: rijv 00483 REAL(dp), INTENT(OUT), OPTIONAL :: enuc 00484 REAL(dp), DIMENSION(3), INTENT(OUT), 00485 OPTIONAL :: denuc 00486 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 00487 TYPE(se_taper_type), POINTER :: se_taper 00488 TYPE(cp_error_type), INTENT(inout) :: error 00489 00490 CHARACTER(len=*), PARAMETER :: routineN = 'corecore_gks', 00491 routineP = moduleN//':'//routineN 00492 00493 LOGICAL :: failure, l_denuc, l_enuc 00494 REAL(dp) :: alpi, alpj, dscale, rij, 00495 scale, zz 00496 REAL(kind=dp), DIMENSION(3, 45, 45) :: dCoul, dCoulE 00497 REAL(kind=dp), DIMENSION(45, 45) :: Coul, CoulE 00498 TYPE(se_int_control_type) :: se_int_control_off 00499 00500 failure = .FALSE. 00501 00502 rij=DOT_PRODUCT(rijv,rijv) 00503 00504 l_enuc = PRESENT(enuc) 00505 l_denuc= PRESENT(denuc) 00506 IF ((rij > rij_threshold).AND.(l_enuc.OR.l_denuc)) THEN 00507 00508 rij = SQRT(rij) 00509 00510 IF(se_int_control%shortrange) THEN 00511 CALL setup_se_int_control_type(se_int_control_off, shortrange=.FALSE., do_ewald_r3=.FALSE.,& 00512 do_ewald_gks=.FALSE., integral_screening=se_int_control%integral_screening,& 00513 max_multipole=do_multipole_none, pc_coulomb_int=.FALSE.) 00514 CALL makeCoul(rijv,sepi,sepj,Coul,se_int_control_off,error) 00515 IF (l_denuc) CALL makedCoul(rijv,sepi,sepj,dCoul,se_int_control_off,error) 00516 IF (se_int_control%do_ewald_gks) THEN 00517 CALL makeCoulE(rijv,sepi,sepj,CoulE,se_int_control,error) 00518 IF (l_denuc) CALL makedCoulE(rijv,sepi,sepj,dCoulE,se_int_control,error) 00519 ELSE 00520 CALL makeCoul(rijv,sepi,sepj,CoulE,se_int_control,error) 00521 IF (l_denuc) CALL makedCoul(rijv,sepi,sepj,dCoulE,se_int_control,error) 00522 END IF 00523 ELSE 00524 CALL makeCoul(rijv,sepi,sepj,Coul,se_int_control,error) 00525 CoulE = Coul 00526 IF (l_denuc) CALL makedCoul(rijv,sepi,sepj,dCoul,se_int_control,error) 00527 IF (l_denuc) dCoulE = dCoul 00528 END IF 00529 00530 scale = 0.0_dp 00531 dscale= 0.0_dp 00532 zz = sepi%zeff*sepj%zeff 00533 alpi = sepi%alp 00534 alpj = sepj%alp 00535 scale = EXP(-alpi*rij)+EXP(-alpj*rij) 00536 IF (l_enuc) THEN 00537 enuc=zz*CoulE(1,1)+scale*zz*Coul(1,1) 00538 END IF 00539 IF (l_denuc) THEN 00540 dscale= -alpi*EXP(-alpi*rij)-alpj*EXP(-alpj*rij) 00541 denuc(1)= zz*dCoulE(1,1,1)+dscale*(rijv(1)/rij)*zz*Coul(1,1)+scale*zz*dCoul(1,1,1) 00542 denuc(2)= zz*dCoulE(2,1,1)+dscale*(rijv(2)/rij)*zz*Coul(1,1)+scale*zz*dCoul(2,1,1) 00543 denuc(3)= zz*dCoulE(3,1,1)+dscale*(rijv(3)/rij)*zz*Coul(1,1)+scale*zz*dCoul(3,1,1) 00544 END IF 00545 00546 ELSE 00547 00548 IF (se_int_control%do_ewald_gks) THEN 00549 zz = sepi%zeff*sepi%zeff 00550 CALL makeCoulE0(sepi,CoulE,se_int_control,error) 00551 IF (l_enuc) THEN 00552 enuc=zz*CoulE(1,1) 00553 END IF 00554 IF (l_denuc) THEN 00555 denuc=0._dp 00556 END IF 00557 END IF 00558 00559 ENDIF 00560 END SUBROUTINE corecore_gks 00561 00562 ! ***************************************************************************** 00565 SUBROUTINE makeCoulE(RAB,sepi,sepj,Coul,se_int_control,error) 00566 REAL(KIND=dp), DIMENSION(3) :: RAB 00567 TYPE(semi_empirical_type), POINTER :: sepi, sepj 00568 REAL(KIND=dp), DIMENSION(45, 45), 00569 INTENT(OUT) :: Coul 00570 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 00571 TYPE(cp_error_type), INTENT(inout) :: error 00572 00573 CHARACTER(len=*), PARAMETER :: routineN = 'makeCoulE', 00574 routineP = moduleN//':'//routineN 00575 00576 INTEGER :: gpt, imA, imB, k1, k2, k3, 00577 k4, lp, mp, np 00578 INTEGER, DIMENSION(:, :), POINTER :: bds 00579 REAL(KIND=dp) :: a2, ACOULA, ACOULB, alpha, cc, d1, d1f(3), d2, d2f(3,3), 00580 d3, d3f(3,3,3), d4, d4f(3,3,3,3), f, ff, kr, kr2, r1, r2, r3, r5, r7, 00581 r9, rr, ss, w, w0, w1, w2, w3, w4, w5 00582 REAL(KIND=dp), DIMENSION(3) :: kk, v 00583 REAL(KIND=dp), DIMENSION(3, 3, 45) :: M2A, M2B 00584 REAL(KIND=dp), DIMENSION(3, 45) :: M1A, M1B 00585 REAL(KIND=dp), DIMENSION(45) :: M0A, M0B 00586 REAL(KIND=dp), DIMENSION(:, :, :), 00587 POINTER :: rho0 00588 TYPE(dg_rho0_type), POINTER :: dg_rho0 00589 TYPE(dg_type), POINTER :: dg 00590 TYPE(pw_grid_type), POINTER :: pw_grid 00591 TYPE(pw_pool_type), POINTER :: pw_pool 00592 00593 alpha = se_int_control%ewald_gks%alpha 00594 pw_pool => se_int_control%ewald_gks%pw_pool 00595 dg => se_int_control%ewald_gks%dg 00596 CALL dg_get (dg, dg_rho0=dg_rho0) 00597 rho0 => dg_rho0%density%pw%cr3d 00598 pw_grid => pw_pool%pw_grid 00599 bds => pw_grid%bounds 00600 00601 CALL get_se_slater_multipole(sepi, M0A, M1A, M2A, ACOULA) 00602 CALL get_se_slater_multipole(sepj, M0B, M1B, M2B, ACOULB) 00603 00604 v(1)=RAB(1) 00605 v(2)=RAB(2) 00606 v(3)=RAB(3) 00607 rr=SQRT(DOT_PRODUCT(v,v)) 00608 00609 r1=1.0_dp/rr 00610 r2=r1*r1 00611 r3=r2*r1 00612 r5=r3*r2 00613 r7=r5*r2 00614 r9=r7*r2 00615 00616 a2=0.5_dp*(1.0_dp/ACOULA+1.0_dp/ACOULB) 00617 00618 w0= a2*rr 00619 w= EXP(-w0) 00620 w1= (1.0_dp+0.5_dp*w0) 00621 w2= (w1+0.5_dp*w0+0.5_dp*w0**2) 00622 w3= (w2+w0**3/6.0_dp) 00623 w4= (w3+w0**4/30.0_dp) 00624 w5= (w3+8.0_dp*w0**4/210.0_dp+w0**5/210.0_dp) 00625 00626 f= (1.0_dp-w*w1)*r1 00627 d1= -1.0_dp*(1.0_dp-w*w2)*r3 00628 d2= 3.0_dp*(1.0_dp-w*w3)*r5 00629 d3= -15.0_dp*(1.0_dp-w*w4)*r7 00630 d4= 105.0_dp*(1.0_dp-w*w5)*r9 00631 00632 00633 kr= alpha*rr 00634 kr2= kr*kr 00635 w0= 1.0_dp-erfc(kr) 00636 w1= 2.0_dp*oorootpi*EXP(-kr2) 00637 w2= w1*kr 00638 00639 00640 f=f -w0*r1 00641 d1=d1+ (-w2+w0)*r3 00642 d2=d2+ (w2*(3.0_dp+kr2*2.0_dp)-3.0_dp*w0)*r5 00643 d3=d3+ (-w2*(15.0_dp+kr2*(10.0_dp+kr2*4.0_dp))+15.0_dp*w0)*r7 00644 d4=d4+ (w2*(105.0_dp+kr2*(70.0_dp+kr2*(28.0_dp+kr2*8.0_dp)))-105.0_dp*w0)*r9 00645 00646 CALL build_d_tensor_gks(d1f, d2f, d3f, d4f, v=v, d1=d1, d2=d2, d3=d3, d4=d4) 00647 00648 DO imA=1, (sepi%natorb*(sepi%natorb+1))/2 00649 DO imB=1, (sepj%natorb*(sepj%natorb+1))/2 00650 00651 w= M0A(imA)*M0B(imB)*f 00652 DO k1=1,3 00653 w=w+( M1A(k1,imA)*M0B(imB)-M0A(imA)*M1B(k1,imB) )*d1f(k1) 00654 ENDDO 00655 DO k2=1,3 00656 DO k1=1,3 00657 w=w+( M2A(k1,k2,imA)*M0B(imB)-M1A(k1,imA)*M1B(k2,imB)+M0A(imA)*M2B(k1,k2,imB) )*d2f(k1,k2) 00658 ENDDO 00659 ENDDO 00660 DO k3=1,3 00661 DO k2=1,3 00662 DO k1=1,3 00663 w=w+( -M2A(k1,k2,imA)*M1B(k3,imB)+M1A(k1,imA)*M2B(k2,k3,imB) )*d3f(k1,k2,k3) 00664 ENDDO 00665 ENDDO 00666 ENDDO 00667 00668 DO k4=1,3 00669 DO k3=1,3 00670 DO k2=1,3 00671 DO k1=1,3 00672 w=w+ M2A(k1,k2,imA)*M2B(k3,k4,imB)*d4f(k1,k2,k3,k4) 00673 ENDDO 00674 ENDDO 00675 ENDDO 00676 ENDDO 00677 00678 Coul(imA,imB)=w 00679 ENDDO 00680 ENDDO 00681 00682 v(1)=RAB(1) 00683 v(2)=RAB(2) 00684 v(3)=RAB(3) 00685 00686 f= 0.0_dp 00687 d1f = 0.0_dp 00688 d2f = 0.0_dp 00689 d3f = 0.0_dp 00690 d4f = 0.0_dp 00691 00692 DO gpt = 1, pw_grid%ngpts_cut 00693 lp = pw_grid%mapl%pos(pw_grid%g_hat(1,gpt)) 00694 mp = pw_grid%mapm%pos(pw_grid%g_hat(2,gpt)) 00695 np = pw_grid%mapn%pos(pw_grid%g_hat(3,gpt)) 00696 00697 lp = lp + bds(1,1) 00698 mp = mp + bds(1,2) 00699 np = np + bds(1,3) 00700 00701 IF (pw_grid%gsq(gpt) == 0.0_dp) CYCLE 00702 kk(:)=pw_grid%g(:,gpt) 00703 ff = 2.0_dp * fourpi * rho0(lp,mp,np)**2 * pw_grid%vol / pw_grid%gsq(gpt) 00704 00705 kr=DOT_PRODUCT(kk,v) 00706 cc=COS(kr) 00707 ss=SIN(kr) 00708 00709 f=f+cc*ff 00710 DO k1=1,3 00711 d1f(k1)=d1f(k1)-kk(k1)*ss*ff 00712 ENDDO 00713 DO k2=1,3 00714 DO k1=1,3 00715 d2f(k1,k2)=d2f(k1,k2)-kk(k1)*kk(k2)*cc*ff 00716 ENDDO 00717 ENDDO 00718 DO k3=1,3 00719 DO k2=1,3 00720 DO k1=1,3 00721 d3f(k1,k2,k3)=d3f(k1,k2,k3)+kk(k1)*kk(k2)*kk(k3)*ss*ff 00722 ENDDO 00723 ENDDO 00724 ENDDO 00725 DO k4=1,3 00726 DO k3=1,3 00727 DO k2=1,3 00728 DO k1=1,3 00729 d4f(k1,k2,k3,k4)=d4f(k1,k2,k3,k4)+kk(k1)*kk(k2)*kk(k3)*kk(k4)*cc*ff 00730 ENDDO 00731 ENDDO 00732 ENDDO 00733 ENDDO 00734 00735 ENDDO 00736 00737 DO imA=1, (sepi%natorb*(sepi%natorb+1))/2 00738 DO imB=1, (sepj%natorb*(sepj%natorb+1))/2 00739 00740 w= M0A(imA)*M0B(imB)*f 00741 DO k1=1,3 00742 w=w+( M1A(k1,imA)*M0B(imB)-M0A(imA)*M1B(k1,imB) )*d1f(k1) 00743 ENDDO 00744 DO k2=1,3 00745 DO k1=1,3 00746 w=w+( M2A(k1,k2,imA)*M0B(imB)-M1A(k1,imA)*M1B(k2,imB)+M0A(imA)*M2B(k1,k2,imB) )*d2f(k1,k2) 00747 ENDDO 00748 ENDDO 00749 DO k3=1,3 00750 DO k2=1,3 00751 DO k1=1,3 00752 w=w+( -M2A(k1,k2,imA)*M1B(k3,imB)+M1A(k1,imA)*M2B(k2,k3,imB) )*d3f(k1,k2,k3) 00753 ENDDO 00754 ENDDO 00755 ENDDO 00756 00757 DO k4=1,3 00758 DO k3=1,3 00759 DO k2=1,3 00760 DO k1=1,3 00761 w=w+ M2A(k1,k2,imA)*M2B(k3,k4,imB)*d4f(k1,k2,k3,k4) 00762 ENDDO 00763 ENDDO 00764 ENDDO 00765 ENDDO 00766 00767 Coul(imA,imB)=Coul(imA,imB)+w 00768 00769 ENDDO 00770 ENDDO 00771 00772 DO imA=1, (sepi%natorb*(sepi%natorb+1))/2 00773 DO imB=1, (sepj%natorb*(sepj%natorb+1))/2 00774 w= -M0A(imA)*M0B(imB)*0.25_dp*fourpi/(pw_grid%vol*alpha**2) 00775 Coul(imA,imB)=Coul(imA,imB)+w 00776 ENDDO 00777 ENDDO 00778 00779 END SUBROUTINE makeCoulE 00780 00781 ! ***************************************************************************** 00785 SUBROUTINE makedCoulE(RAB,sepi,sepj,dCoul,se_int_control,error) 00786 REAL(KIND=dp), DIMENSION(3) :: RAB 00787 TYPE(semi_empirical_type), POINTER :: sepi, sepj 00788 REAL(KIND=dp), DIMENSION(3, 45, 45), 00789 INTENT(OUT) :: dCoul 00790 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 00791 TYPE(cp_error_type), INTENT(inout) :: error 00792 00793 CHARACTER(len=*), PARAMETER :: routineN = 'makedCoulE', 00794 routineP = moduleN//':'//routineN 00795 00796 INTEGER :: gpt, imA, imB, k1, k2, k3, 00797 k4, k5, lp, mp, np 00798 INTEGER, DIMENSION(:, :), POINTER :: bds 00799 REAL(KIND=dp) :: a2, ACOULA, ACOULB, alpha, cc, d1, d1f(3), d2, d2f(3,3), 00800 d3, d3f(3,3,3), d4, d4f(3,3,3,3), d5, d5f(3,3,3,3,3), f, ff, kr, kr2, 00801 r1, r11, r2, r3, r5, r7, r9, rr, ss, tmp, w, w0, w1, w2, w3, w4, w5, w6 00802 REAL(KIND=dp), DIMENSION(3) :: kk, v, wv 00803 REAL(kind=dp), DIMENSION(3, 3, 45) :: M2A, M2B 00804 REAL(kind=dp), DIMENSION(3, 45) :: M1A, M1B 00805 REAL(kind=dp), DIMENSION(45) :: M0A, M0B 00806 REAL(KIND=dp), DIMENSION(:, :, :), 00807 POINTER :: rho0 00808 TYPE(dg_rho0_type), POINTER :: dg_rho0 00809 TYPE(dg_type), POINTER :: dg 00810 TYPE(pw_grid_type), POINTER :: pw_grid 00811 TYPE(pw_pool_type), POINTER :: pw_pool 00812 00813 alpha = se_int_control%ewald_gks%alpha 00814 pw_pool => se_int_control%ewald_gks%pw_pool 00815 dg => se_int_control%ewald_gks%dg 00816 CALL dg_get (dg, dg_rho0=dg_rho0) 00817 rho0 => dg_rho0%density%pw%cr3d 00818 pw_grid => pw_pool%pw_grid 00819 bds => pw_grid%bounds 00820 00821 CALL get_se_slater_multipole(sepi, M0A, M1A, M2A, ACOULA) 00822 CALL get_se_slater_multipole(sepj, M0B, M1B, M2B, ACOULB) 00823 00824 v(1)=RAB(1) 00825 v(2)=RAB(2) 00826 v(3)=RAB(3) 00827 rr=SQRT(DOT_PRODUCT(v,v)) 00828 00829 a2=0.5_dp*(1.0_dp/ACOULA+1.0_dp/ACOULB) 00830 00831 r1=1.0_dp/rr 00832 r2=r1*r1 00833 r3=r2*r1 00834 r5=r3*r2 00835 r7=r5*r2 00836 r9=r7*r2 00837 r11=r9*r2 00838 00839 w0= a2*rr 00840 w= EXP(-w0) 00841 w1= (1.0_dp+0.5_dp*w0) 00842 w2= (w1+0.5_dp*w0+0.5_dp*w0**2) 00843 w3= (w2+w0**3/6.0_dp) 00844 w4= (w3+w0**4/30.0_dp) 00845 w5= (w3+8.0_dp*w0**4/210.0_dp+w0**5/210.0_dp) 00846 w6= (w3+5.0_dp*w0**4/126.0_dp+2.0_dp*w0**5/315.0_dp+w0**6/1890.0_dp) 00847 00848 f= (1.0_dp-w*w1)*r1 00849 d1= -1.0_dp*(1.0_dp-w*w2)*r3 00850 d2= 3.0_dp*(1.0_dp-w*w3)*r5 00851 d3= -15.0_dp*(1.0_dp-w*w4)*r7 00852 d4= 105.0_dp*(1.0_dp-w*w5)*r9 00853 d5= -945.0_dp*(1.0_dp-w*w6)*r11 00854 00855 kr= alpha*rr 00856 kr2= kr*kr 00857 w0= 1.0_dp-erfc(kr) 00858 w1= 2.0_dp*oorootpi*EXP(-kr2) 00859 w2= w1*kr 00860 00861 f=f -w0*r1 00862 d1=d1+ (-w2+w0)*r3 00863 d2=d2+ (w2*(3.0_dp+kr2*2.0_dp)-3.0_dp*w0)*r5 00864 d3=d3+ (-w2*(15.0_dp+kr2*(10.0_dp+kr2*4.0_dp))+15.0_dp*w0)*r7 00865 d4=d4+ (w2*(105.0_dp+kr2*(70.0_dp+kr2*(28.0_dp+kr2*8.0_dp)))-105.0_dp*w0)*r9 00866 d5=d5+ (-w2*(945.0_dp+kr2*(630.0_dp+kr2*(252.0_dp+kr2*(72.0_dp+kr2*16.0_dp))))+945.0_dp*w0)*r11 00867 00868 CALL build_d_tensor_gks(d1f, d2f, d3f, d4f, d5f, v, d1, d2, d3, d4, d5) 00869 00870 DO imA=1, (sepi%natorb*(sepi%natorb+1))/2 00871 DO imB=1, (sepj%natorb*(sepj%natorb+1))/2 00872 00873 tmp = M0A(imA)*M0B(imB) 00874 wv(1)= tmp*d1f(1) 00875 wv(2)= tmp*d1f(2) 00876 wv(3)= tmp*d1f(3) 00877 00878 DO k1=1,3 00879 tmp = M1A(k1,imA)*M0B(imB)-M0A(imA)*M1B(k1,imB) 00880 wv(1)=wv(1)+tmp*d2f(1,k1) 00881 wv(2)=wv(2)+tmp*d2f(2,k1) 00882 wv(3)=wv(3)+tmp*d2f(3,k1) 00883 ENDDO 00884 DO k2=1,3 00885 DO k1=1,3 00886 tmp = M2A(k1,k2,imA)*M0B(imB)-M1A(k1,imA)*M1B(k2,imB)+M0A(imA)*M2B(k1,k2,imB) 00887 wv(1)=wv(1)+tmp*d3f(1,k1,k2) 00888 wv(2)=wv(2)+tmp*d3f(2,k1,k2) 00889 wv(3)=wv(3)+tmp*d3f(3,k1,k2) 00890 ENDDO 00891 ENDDO 00892 DO k3=1,3 00893 DO k2=1,3 00894 DO k1=1,3 00895 tmp = -M2A(k1,k2,imA)*M1B(k3,imB)+M1A(k1,imA)*M2B(k2,k3,imB) 00896 wv(1)=wv(1)+tmp*d4f(1,k1,k2,k3) 00897 wv(2)=wv(2)+tmp*d4f(2,k1,k2,k3) 00898 wv(3)=wv(3)+tmp*d4f(3,k1,k2,k3) 00899 ENDDO 00900 ENDDO 00901 ENDDO 00902 00903 DO k4=1,3 00904 DO k3=1,3 00905 DO k2=1,3 00906 DO k1=1,3 00907 tmp = M2A(k1,k2,imA)*M2B(k3,k4,imB) 00908 wv(1)=wv(1)+ tmp*d5f(1,k1,k2,k3,k4) 00909 wv(2)=wv(2)+ tmp*d5f(2,k1,k2,k3,k4) 00910 wv(3)=wv(3)+ tmp*d5f(3,k1,k2,k3,k4) 00911 ENDDO 00912 ENDDO 00913 ENDDO 00914 ENDDO 00915 00916 dCoul(1,imA,imB)=wv(1) 00917 dCoul(2,imA,imB)=wv(2) 00918 dCoul(3,imA,imB)=wv(3) 00919 ENDDO 00920 ENDDO 00921 00922 v(1)=RAB(1) 00923 v(2)=RAB(2) 00924 v(3)=RAB(3) 00925 00926 f= 0.0_dp 00927 d1f = 0.0_dp 00928 d2f = 0.0_dp 00929 d3f = 0.0_dp 00930 d4f = 0.0_dp 00931 d5f = 0.0_dp 00932 00933 DO gpt = 1, pw_grid%ngpts_cut 00934 lp = pw_grid%mapl%pos(pw_grid%g_hat(1,gpt)) 00935 mp = pw_grid%mapm%pos(pw_grid%g_hat(2,gpt)) 00936 np = pw_grid%mapn%pos(pw_grid%g_hat(3,gpt)) 00937 00938 lp = lp + bds(1,1) 00939 mp = mp + bds(1,2) 00940 np = np + bds(1,3) 00941 00942 IF (pw_grid%gsq(gpt) == 0.0_dp) CYCLE 00943 kk(:)=pw_grid%g(:,gpt) 00944 ff = 2.0_dp * fourpi * rho0(lp,mp,np)**2 * pw_grid%vol / pw_grid%gsq(gpt) 00945 00946 kr=DOT_PRODUCT(kk,v) 00947 cc=COS(kr) 00948 ss=SIN(kr) 00949 00950 f=f+cc*ff 00951 DO k1=1,3 00952 d1f(k1)=d1f(k1)-kk(k1)*ss*ff 00953 ENDDO 00954 DO k2=1,3 00955 DO k1=1,3 00956 d2f(k1,k2)=d2f(k1,k2)-kk(k1)*kk(k2)*cc*ff 00957 ENDDO 00958 ENDDO 00959 DO k3=1,3 00960 DO k2=1,3 00961 DO k1=1,3 00962 d3f(k1,k2,k3)=d3f(k1,k2,k3)+kk(k1)*kk(k2)*kk(k3)*ss*ff 00963 ENDDO 00964 ENDDO 00965 ENDDO 00966 DO k4=1,3 00967 DO k3=1,3 00968 DO k2=1,3 00969 DO k1=1,3 00970 d4f(k1,k2,k3,k4)=d4f(k1,k2,k3,k4)+kk(k1)*kk(k2)*kk(k3)*kk(k4)*cc*ff 00971 ENDDO 00972 ENDDO 00973 ENDDO 00974 ENDDO 00975 DO k5=1,3 00976 DO k4=1,3 00977 DO k3=1,3 00978 DO k2=1,3 00979 DO k1=1,3 00980 d5f(k1,k2,k3,k4,k5)=d5f(k1,k2,k3,k4,k5)-kk(k1)*kk(k2)*kk(k3)*kk(k4)*kk(k5)*ss*ff 00981 ENDDO 00982 ENDDO 00983 ENDDO 00984 ENDDO 00985 ENDDO 00986 ENDDO 00987 00988 DO imA=1, (sepi%natorb*(sepi%natorb+1))/2 00989 DO imB=1, (sepj%natorb*(sepj%natorb+1))/2 00990 tmp = M0A(imA)*M0B(imB) 00991 wv(1)= tmp*d1f(1) 00992 wv(2)= tmp*d1f(2) 00993 wv(3)= tmp*d1f(3) 00994 DO k1=1,3 00995 tmp = M1A(k1,imA)*M0B(imB)-M0A(imA)*M1B(k1,imB) 00996 wv(1)=wv(1)+tmp*d2f(1,k1) 00997 wv(2)=wv(2)+tmp*d2f(2,k1) 00998 wv(3)=wv(3)+tmp*d2f(3,k1) 00999 ENDDO 01000 DO k2=1,3 01001 DO k1=1,3 01002 tmp = M2A(k1,k2,imA)*M0B(imB)-M1A(k1,imA)*M1B(k2,imB)+M0A(imA)*M2B(k1,k2,imB) 01003 wv(1)=wv(1)+tmp*d3f(1,k1,k2) 01004 wv(2)=wv(2)+tmp*d3f(2,k1,k2) 01005 wv(3)=wv(3)+tmp*d3f(3,k1,k2) 01006 ENDDO 01007 ENDDO 01008 DO k3=1,3 01009 DO k2=1,3 01010 DO k1=1,3 01011 tmp = -M2A(k1,k2,imA)*M1B(k3,imB)+M1A(k1,imA)*M2B(k2,k3,imB) 01012 wv(1)=wv(1)+tmp*d4f(1,k1,k2,k3) 01013 wv(2)=wv(2)+tmp*d4f(2,k1,k2,k3) 01014 wv(3)=wv(3)+tmp*d4f(3,k1,k2,k3) 01015 ENDDO 01016 ENDDO 01017 ENDDO 01018 01019 DO k4=1,3 01020 DO k3=1,3 01021 DO k2=1,3 01022 DO k1=1,3 01023 tmp = M2A(k1,k2,imA)*M2B(k3,k4,imB) 01024 wv(1)=wv(1)+tmp*d5f(1,k1,k2,k3,k4) 01025 wv(2)=wv(2)+tmp*d5f(2,k1,k2,k3,k4) 01026 wv(3)=wv(3)+tmp*d5f(3,k1,k2,k3,k4) 01027 ENDDO 01028 ENDDO 01029 ENDDO 01030 ENDDO 01031 01032 dCoul(1,imA,imB)=dCoul(1,imA,imB)+wv(1) 01033 dCoul(2,imA,imB)=dCoul(2,imA,imB)+wv(2) 01034 dCoul(3,imA,imB)=dCoul(3,imA,imB)+wv(3) 01035 ENDDO 01036 ENDDO 01037 01038 END SUBROUTINE makedCoulE 01039 01040 ! ***************************************************************************** 01044 SUBROUTINE build_d_tensor_gks(d1f, d2f, d3f, d4f, d5f, v, d1, d2, d3, d4, d5) 01045 REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: d1f 01046 REAL(KIND=dp), DIMENSION(3, 3), 01047 INTENT(OUT) :: d2f 01048 REAL(KIND=dp), DIMENSION(3, 3, 3), 01049 INTENT(OUT) :: d3f 01050 REAL(KIND=dp), DIMENSION(3, 3, 3, 3), 01051 INTENT(OUT) :: d4f 01052 REAL(KIND=dp), 01053 DIMENSION(3, 3, 3, 3, 3), 01054 INTENT(OUT), OPTIONAL :: d5f 01055 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: v 01056 REAL(KIND=dp), INTENT(IN) :: d1, d2, d3, d4 01057 REAL(KIND=dp), INTENT(IN), OPTIONAL :: d5 01058 01059 INTEGER :: k1, k2, k3, k4, k5 01060 REAL(KIND=dp) :: w 01061 01062 d1f = 0.0_dp 01063 d2f = 0.0_dp 01064 d3f = 0.0_dp 01065 d4f = 0.0_dp 01066 DO k1=1,3 01067 d1f(k1)=d1f(k1)+v(k1)*d1 01068 ENDDO 01069 DO k1=1,3 01070 DO k2=1,3 01071 d2f(k2,k1)=d2f(k2,k1)+v(k1)*v(k2)*d2 01072 ENDDO 01073 d2f(k1,k1)=d2f(k1,k1)+ d1 01074 ENDDO 01075 DO k1=1,3 01076 DO k2=1,3 01077 DO k3=1,3 01078 d3f(k3,k2,k1)=d3f(k3,k2,k1)+v(k1)*v(k2)*v(k3)*d3 01079 ENDDO 01080 w=v(k1)*d2 01081 d3f(k1,k2,k2)=d3f(k1,k2,k2)+w 01082 d3f(k2,k1,k2)=d3f(k2,k1,k2)+w 01083 d3f(k2,k2,k1)=d3f(k2,k2,k1)+w 01084 ENDDO 01085 ENDDO 01086 DO k1=1,3 01087 DO k2=1,3 01088 DO k3=1,3 01089 DO k4=1,3 01090 d4f(k4,k3,k2,k1)=d4f(k4,k3,k2,k1)+v(k1)*v(k2)*v(k3)*v(k4)*d4 01091 ENDDO 01092 w=v(k1)*v(k2)*d3 01093 d4f(k1,k2,k3,k3)=d4f(k1,k2,k3,k3)+w 01094 d4f(k1,k3,k2,k3)=d4f(k1,k3,k2,k3)+w 01095 d4f(k3,k1,k2,k3)=d4f(k3,k1,k2,k3)+w 01096 d4f(k1,k3,k3,k2)=d4f(k1,k3,k3,k2)+w 01097 d4f(k3,k1,k3,k2)=d4f(k3,k1,k3,k2)+w 01098 d4f(k3,k3,k1,k2)=d4f(k3,k3,k1,k2)+w 01099 ENDDO 01100 d4f(k1,k1,k2,k2)=d4f(k1,k1,k2,k2)+d2 01101 d4f(k1,k2,k1,k2)=d4f(k1,k2,k1,k2)+d2 01102 d4f(k1,k2,k2,k1)=d4f(k1,k2,k2,k1)+d2 01103 ENDDO 01104 ENDDO 01105 IF (PRESENT(d5f).AND.PRESENT(d5)) THEN 01106 d5f = 0.0_dp 01107 01108 DO k1=1,3 01109 DO k2=1,3 01110 DO k3=1,3 01111 DO k4=1,3 01112 DO k5=1,3 01113 d5f(k5,k4,k3,k2,k1)=d5f(k5,k4,k3,k2,k1)+v(k1)*v(k2)*v(k3)*v(k4)*v(k5)*d5 01114 ENDDO 01115 w=v(k1)*v(k2)*v(k3)*d4 01116 d5f(k1,k2,k3,k4,k4)=d5f(k1,k2,k3,k4,k4)+w 01117 d5f(k1,k2,k4,k3,k4)=d5f(k1,k2,k4,k3,k4)+w 01118 d5f(k1,k4,k2,k3,k4)=d5f(k1,k4,k2,k3,k4)+w 01119 d5f(k4,k1,k2,k3,k4)=d5f(k4,k1,k2,k3,k4)+w 01120 d5f(k1,k2,k4,k4,k3)=d5f(k1,k2,k4,k4,k3)+w 01121 d5f(k1,k4,k2,k4,k3)=d5f(k1,k4,k2,k4,k3)+w 01122 d5f(k4,k1,k2,k4,k3)=d5f(k4,k1,k2,k4,k3)+w 01123 d5f(k1,k4,k4,k2,k3)=d5f(k1,k4,k4,k2,k3)+w 01124 d5f(k4,k1,k4,k2,k3)=d5f(k4,k1,k4,k2,k3)+w 01125 d5f(k4,k4,k1,k2,k3)=d5f(k4,k4,k1,k2,k3)+w 01126 ENDDO 01127 w=v(k1)*d3 01128 d5f(k1,k2,k2,k3,k3)=d5f(k1,k2,k2,k3,k3)+w 01129 d5f(k1,k2,k3,k2,k3)=d5f(k1,k2,k3,k2,k3)+w 01130 d5f(k1,k2,k3,k3,k2)=d5f(k1,k2,k3,k3,k2)+w 01131 d5f(k2,k1,k2,k3,k3)=d5f(k2,k1,k2,k3,k3)+w 01132 d5f(k2,k1,k3,k2,k3)=d5f(k2,k1,k3,k2,k3)+w 01133 d5f(k2,k1,k3,k3,k2)=d5f(k2,k1,k3,k3,k2)+w 01134 d5f(k2,k2,k1,k3,k3)=d5f(k2,k2,k1,k3,k3)+w 01135 d5f(k2,k3,k1,k2,k3)=d5f(k2,k3,k1,k2,k3)+w 01136 d5f(k2,k3,k1,k3,k2)=d5f(k2,k3,k1,k3,k2)+w 01137 d5f(k2,k2,k3,k1,k3)=d5f(k2,k2,k3,k1,k3)+w 01138 d5f(k2,k3,k2,k1,k3)=d5f(k2,k3,k2,k1,k3)+w 01139 d5f(k2,k3,k3,k1,k2)=d5f(k2,k3,k3,k1,k2)+w 01140 d5f(k2,k2,k3,k3,k1)=d5f(k2,k2,k3,k3,k1)+w 01141 d5f(k2,k3,k2,k3,k1)=d5f(k2,k3,k2,k3,k1)+w 01142 d5f(k2,k3,k3,k2,k1)=d5f(k2,k3,k3,k2,k1)+w 01143 ENDDO 01144 ENDDO 01145 ENDDO 01146 END IF 01147 END SUBROUTINE build_d_tensor_gks 01148 01149 SUBROUTINE makeCoulE0(sepi,Coul,se_int_control,error) 01150 TYPE(semi_empirical_type), POINTER :: sepi 01151 REAL(KIND=dp), DIMENSION(45, 45), 01152 INTENT(OUT) :: Coul 01153 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 01154 TYPE(cp_error_type), INTENT(inout) :: error 01155 01156 CHARACTER(len=*), PARAMETER :: routineN = 'makeCoulE0', 01157 routineP = moduleN//':'//routineN 01158 01159 INTEGER :: gpt, imA, imB, k1, k2, k3, 01160 k4, lp, mp, np 01161 INTEGER, DIMENSION(:, :), POINTER :: bds 01162 REAL(KIND=dp) :: alpha, d2f(3,3), 01163 d4f(3,3,3,3), f, ff, w 01164 REAL(KIND=dp), DIMENSION(3) :: kk 01165 REAL(KIND=dp), DIMENSION(3, 3, 45) :: M2A 01166 REAL(KIND=dp), DIMENSION(3, 45) :: M1A 01167 REAL(KIND=dp), DIMENSION(45) :: M0A 01168 REAL(KIND=dp), DIMENSION(:, :, :), 01169 POINTER :: rho0 01170 TYPE(dg_rho0_type), POINTER :: dg_rho0 01171 TYPE(dg_type), POINTER :: dg 01172 TYPE(pw_grid_type), POINTER :: pw_grid 01173 TYPE(pw_pool_type), POINTER :: pw_pool 01174 01175 alpha = se_int_control%ewald_gks%alpha 01176 pw_pool => se_int_control%ewald_gks%pw_pool 01177 dg => se_int_control%ewald_gks%dg 01178 CALL dg_get (dg, dg_rho0=dg_rho0) 01179 rho0 => dg_rho0%density%pw%cr3d 01180 pw_grid => pw_pool%pw_grid 01181 bds => pw_grid%bounds 01182 01183 CALL get_se_slater_multipole(sepi, M0A, M1A, M2A) 01184 01185 f= 0.0_dp 01186 d2f = 0.0_dp 01187 d4f = 0.0_dp 01188 01189 DO gpt = 1, pw_grid%ngpts_cut 01190 lp = pw_grid%mapl%pos(pw_grid%g_hat(1,gpt)) 01191 mp = pw_grid%mapm%pos(pw_grid%g_hat(2,gpt)) 01192 np = pw_grid%mapn%pos(pw_grid%g_hat(3,gpt)) 01193 01194 lp = lp + bds(1,1) 01195 mp = mp + bds(1,2) 01196 np = np + bds(1,3) 01197 01198 IF (pw_grid%gsq(gpt) == 0.0_dp) CYCLE 01199 kk(:)=pw_grid%g(:,gpt) 01200 ff = 2.0_dp * fourpi * rho0(lp,mp,np)**2 * pw_grid%vol / pw_grid%gsq(gpt) 01201 01202 f=f+ff 01203 DO k2=1,3 01204 DO k1=1,3 01205 d2f(k1,k2)=d2f(k1,k2)-kk(k1)*kk(k2)*ff 01206 ENDDO 01207 ENDDO 01208 DO k4=1,3 01209 DO k3=1,3 01210 DO k2=1,3 01211 DO k1=1,3 01212 d4f(k1,k2,k3,k4)=d4f(k1,k2,k3,k4)+kk(k1)*kk(k2)*kk(k3)*kk(k4)*ff 01213 ENDDO 01214 ENDDO 01215 ENDDO 01216 ENDDO 01217 01218 ENDDO 01219 01220 DO imA=1, (sepi%natorb*(sepi%natorb+1))/2 01221 DO imB=1, (sepi%natorb*(sepi%natorb+1))/2 01222 01223 w= M0A(imA)*M0A(imB)*f 01224 DO k2=1,3 01225 DO k1=1,3 01226 w=w+( M2A(k1,k2,imA)*M0A(imB)-M1A(k1,imA)*M1A(k2,imB)+M0A(imA)*M2A(k1,k2,imB) )*d2f(k1,k2) 01227 ENDDO 01228 ENDDO 01229 01230 DO k4=1,3 01231 DO k3=1,3 01232 DO k2=1,3 01233 DO k1=1,3 01234 w=w+ M2A(k1,k2,imA)*M2A(k3,k4,imB)*d4f(k1,k2,k3,k4) 01235 ENDDO 01236 ENDDO 01237 ENDDO 01238 ENDDO 01239 01240 Coul(imA,imB)=w 01241 01242 ENDDO 01243 ENDDO 01244 01245 DO imA=1, (sepi%natorb*(sepi%natorb+1))/2 01246 DO imB=1, (sepi%natorb*(sepi%natorb+1))/2 01247 w= -M0A(imA)*M0A(imB)*0.25_dp*fourpi/(pw_grid%vol*alpha**2) 01248 Coul(imA,imB)=Coul(imA,imB)+w 01249 ENDDO 01250 ENDDO 01251 01252 DO imA=1, (sepi%natorb*(sepi%natorb+1))/2 01253 DO imB=1, (sepi%natorb*(sepi%natorb+1))/2 01254 01255 w= M0A(imA)*M0A(imB) 01256 Coul(imA,imB)=Coul(imA,imB)-2.0_dp*alpha*oorootpi*w 01257 01258 w= 0.0_dp 01259 DO k1=1,3 01260 w=w+ M1A(k1,imA)*M1A(k1,imB) 01261 w=w- M0A(imA)*M2A(k1,k1,imB) 01262 w=w- M2A(k1,k1,imA)*M0A(imB) 01263 ENDDO 01264 Coul(imA,imB)=Coul(imA,imB)-4.0_dp*alpha**3*oorootpi*w/3.0_dp 01265 01266 w= 0.0_dp 01267 DO k2=1,3 01268 DO k1=1,3 01269 w=w+ 2.0_dp*M2A(k1,k2,imA)*M2A(k1,k2,imB) 01270 w=w+ M2A(k1,k1,imA)*M2A(k2,k2,imB) 01271 ENDDO 01272 ENDDO 01273 Coul(imA,imB)=Coul(imA,imB)-8.0_dp*alpha**5*oorootpi*w/5.0_dp 01274 ENDDO 01275 ENDDO 01276 END SUBROUTINE makeCoulE0 01277 01278 ! ***************************************************************************** 01281 SUBROUTINE get_se_slater_multipole(sepi, M0, M1, M2, ACOUL) 01282 TYPE(semi_empirical_type), POINTER :: sepi 01283 REAL(kind=dp), DIMENSION(45), 01284 INTENT(OUT) :: M0 01285 REAL(kind=dp), DIMENSION(3, 45), 01286 INTENT(OUT) :: M1 01287 REAL(kind=dp), DIMENSION(3, 3, 45), 01288 INTENT(OUT) :: M2 01289 REAL(KIND=dp), INTENT(OUT), OPTIONAL :: ACOUL 01290 01291 INTEGER :: i, j, jint, size_1c_int 01292 TYPE(semi_empirical_mpole_type), POINTER :: mpole 01293 01294 NULLIFY(mpole) 01295 size_1c_int = SIZE(sepi%w_mpole) 01296 DO jint = 1, size_1c_int 01297 mpole => sepi%w_mpole(jint)%mpole 01298 i=mpole%indi 01299 j=mpole%indj 01300 M0(indexb(i,j)) = -mpole%cs 01301 01302 M1(1,indexb(i,j)) = -mpole%ds(1) 01303 M1(2,indexb(i,j)) = -mpole%ds(2) 01304 M1(3,indexb(i,j)) = -mpole%ds(3) 01305 01306 M2(1,1,indexb(i,j))= -mpole%qq(1,1)/3.0_dp 01307 M2(2,1,indexb(i,j))= -mpole%qq(2,1)/3.0_dp 01308 M2(3,1,indexb(i,j))= -mpole%qq(3,1)/3.0_dp 01309 01310 M2(1,2,indexb(i,j))= -mpole%qq(1,2)/3.0_dp 01311 M2(2,2,indexb(i,j))= -mpole%qq(2,2)/3.0_dp 01312 M2(3,2,indexb(i,j))= -mpole%qq(3,2)/3.0_dp 01313 01314 M2(1,3,indexb(i,j))= -mpole%qq(1,3)/3.0_dp 01315 M2(2,3,indexb(i,j))= -mpole%qq(2,3)/3.0_dp 01316 M2(3,3,indexb(i,j))= -mpole%qq(3,3)/3.0_dp 01317 ENDDO 01318 IF ( PRESENT (ACOUL) ) ACOUL = sepi%acoul 01319 END SUBROUTINE get_se_slater_multipole 01320 01321 END MODULE semi_empirical_int_gks
1.7.3