|
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 ! ***************************************************************************** 00015 MODULE hfx_libint_interface 00016 00017 USE cell_types, ONLY: cell_type,& 00018 real_to_scaled 00019 USE f77_blas 00020 USE gamma, ONLY: fgamma=>fgamma_0 00021 USE hfx_contraction_methods, ONLY: contract 00022 USE hfx_libint_wrapper, ONLY: get_derivs,& 00023 get_eris 00024 USE hfx_libint_wrapper_types, ONLY: lib_deriv,& 00025 lib_int,& 00026 prim_data,& 00027 prim_data_f_size 00028 USE hfx_types, ONLY: hfx_pgf_product_list,& 00029 hfx_potential_type 00030 USE input_constants, ONLY: & 00031 do_hfx_potential_coulomb, do_hfx_potential_gaussian, & 00032 do_hfx_potential_id, do_hfx_potential_long, do_hfx_potential_mix_cl, & 00033 do_hfx_potential_mix_cl_trunc, do_hfx_potential_mix_lg, & 00034 do_hfx_potential_short, do_hfx_potential_truncated 00035 USE kinds, ONLY: dp,& 00036 int_8 00037 USE mathconstants 00038 USE orbital_pointers 00039 USE t_c_g0, ONLY: t_c_g0_n 00040 #include "cp_common_uses.h" 00041 00042 IMPLICIT NONE 00043 PRIVATE 00044 PUBLIC evaluate_eri,& 00045 evaluate_deriv_eri,& 00046 evaluate_eri_screen 00047 00048 INTEGER, DIMENSION(12), PARAMETER :: full_perm1 = (/1,2,3,4,5,6,7,8,9,10,11,12/) 00049 INTEGER, DIMENSION(12), PARAMETER :: full_perm2 = (/4,5,6,1,2,3,7,8,9,10,11,12/) 00050 INTEGER, DIMENSION(12), PARAMETER :: full_perm3 = (/1,2,3,4,5,6,10,11,12,7,8,9/) 00051 INTEGER, DIMENSION(12), PARAMETER :: full_perm4 = (/4,5,6,1,2,3,10,11,12,7,8,9/) 00052 INTEGER, DIMENSION(12), PARAMETER :: full_perm5 = (/7,8,9,10,11,12,1,2,3,4,5,6/) 00053 INTEGER, DIMENSION(12), PARAMETER :: full_perm6 = (/7,8,9,10,11,12,4,5,6,1,2,3/) 00054 INTEGER, DIMENSION(12), PARAMETER :: full_perm7 = (/10,11,12,7,8,9,1,2,3,4,5,6/) 00055 INTEGER, DIMENSION(12), PARAMETER :: full_perm8 = (/10,11,12,7,8,9,4,5,6,1,2,3/) 00056 00057 00058 !*** 00059 00060 CONTAINS 00061 00062 00063 ! ***************************************************************************** 00069 SUBROUTINE build_quartet_data_screen(A,B,C,D,Zeta_A, Zeta_B, Zeta_C, Zeta_D, m_max, & 00070 potential_parameter, prim, R11,R22) 00071 00072 REAL(KIND=dp) :: A(3), B(3), C(3), D(3) 00073 REAL(KIND=dp), INTENT(IN) :: Zeta_A, Zeta_B, Zeta_C, Zeta_D 00074 INTEGER, INTENT(IN) :: m_max 00075 TYPE(hfx_potential_type) :: potential_parameter 00076 TYPE(prim_data) :: prim 00077 REAL(dp) :: R11, R22 00078 00079 INTEGER :: i 00080 LOGICAL :: use_gamma 00081 REAL(KIND=dp) :: AB(3), AB2, CD(3), CD2, Eta, EtaInv, factor, omega2, 00082 omega_corr, omega_corr2, P(3), PQ(3), PQ2, Q(3), R, R1, R2, Rho, 00083 RhoInv, S1234, T, tmp, W(3), Zeta, ZetaInv, ZetapEtaInv 00084 REAL(KIND=dp), 00085 DIMENSION(prim_data_f_size) :: Fm 00086 00087 Zeta = Zeta_A + Zeta_B 00088 ZetaInv = 1.0_dp/Zeta 00089 Eta = Zeta_C + Zeta_D 00090 EtaInv = 1.0_dp/Eta 00091 ZetapEtaInv = Zeta+Eta 00092 ZetapEtaInv = 1.0_dp/ZetapEtaInv 00093 Rho = Zeta*Eta*ZetapEtaInv 00094 RhoInv = 1.0_dp/Rho 00095 00096 DO i=1,3 00097 P(i) = (Zeta_A*A(i) + Zeta_B*B(i))*ZetaInv 00098 Q(i) = (Zeta_C*C(i) + Zeta_D*D(i))*EtaInv 00099 AB(i) = A(i)-B(i) 00100 CD(i) = C(i)-D(i) 00101 PQ(i) = P(i)-Q(i) 00102 W(i) = (Zeta*P(i) + Eta*Q(i))*ZetapEtaInv 00103 END DO 00104 00105 AB2 = DOT_PRODUCT(AB,AB) 00106 CD2 = DOT_PRODUCT(CD,CD) 00107 PQ2 = DOT_PRODUCT(PQ,PQ) 00108 00109 S1234= EXP((-Zeta_A*Zeta_B*ZetaInv*AB2)+(-Zeta_C*Zeta_D*EtaInv*CD2)) 00110 T = Rho*PQ2 00111 00112 00113 SELECT CASE(potential_parameter%potential_type) 00114 CASE(do_hfx_potential_truncated) 00115 R = potential_parameter%cutoff_radius * SQRT(rho) 00116 R1 = R11 00117 R2 = R22 00118 IF( PQ2 > (R1+R2+potential_parameter%cutoff_radius)**2 ) THEN 00119 RETURN 00120 END IF 00121 CALL t_c_g0_n(prim%F(1),use_gamma,R,T,m_max) 00122 IF( use_gamma ) CALL fgamma(m_max, T, prim%F(1)) 00123 factor = 2.0_dp * Pi * RhoInv 00124 CASE(do_hfx_potential_coulomb) 00125 CALL fgamma(m_max,T,prim%F(1)) 00126 factor = 2.0_dp*Pi*RhoInv 00127 CASE(do_hfx_potential_short) 00128 R = potential_parameter%cutoff_radius * SQRT(rho) 00129 R1 = R11 00130 R2 = R22 00131 IF( PQ2 > (R1+R2+potential_parameter%cutoff_radius)**2 ) THEN 00132 RETURN 00133 END IF 00134 CALL fgamma(m_max,T,prim%F(1)) 00135 omega2 = potential_parameter%omega**2 00136 omega_corr2 = omega2/(omega2+Rho) 00137 omega_corr = SQRT(omega_corr2) 00138 T = T*omega_corr2 00139 CALL fgamma(m_max,T,Fm) 00140 tmp = - omega_corr 00141 DO i=1,m_max+1 00142 prim%F(i)=prim%F(i) + Fm(i)*tmp 00143 tmp = tmp * omega_corr2 00144 END DO 00145 factor = 2.0_dp*Pi*RhoInv 00146 CASE(do_hfx_potential_long) 00147 omega2 = potential_parameter%omega**2 00148 omega_corr2 = omega2/(omega2+Rho) 00149 omega_corr = SQRT(omega_corr2) 00150 T = T*omega_corr2 00151 CALL fgamma(m_max,T,prim%F(1)) 00152 tmp = omega_corr 00153 DO i=1,m_max+1 00154 prim%F(i)= prim%F(i)*tmp 00155 tmp = tmp * omega_corr2 00156 END DO 00157 factor = 2.0_dp*Pi*RhoInv 00158 CASE(do_hfx_potential_mix_cl) 00159 CALL fgamma(m_max,T,prim%F(1)) 00160 omega2 = potential_parameter%omega**2 00161 omega_corr2 = omega2/(omega2+Rho) 00162 omega_corr = SQRT(omega_corr2) 00163 T = T*omega_corr2 00164 CALL fgamma(m_max,T,Fm) 00165 tmp = omega_corr 00166 DO i=1,m_max+1 00167 prim%F(i)= prim%F(i)*potential_parameter%scale_coulomb + Fm(i)*tmp*potential_parameter%scale_longrange 00168 tmp = tmp * omega_corr2 00169 END DO 00170 factor = 2.0_dp*Pi*RhoInv 00171 CASE(do_hfx_potential_mix_cl_trunc) 00172 00173 ! truncated 00174 R = potential_parameter%cutoff_radius * SQRT(rho) 00175 R1 = R11 00176 R2 = R22 00177 IF( PQ2 > (R1+R2+potential_parameter%cutoff_radius)**2 ) THEN 00178 RETURN 00179 END IF 00180 CALL t_c_g0_n(prim%F(1),use_gamma,R,T,m_max) 00181 IF( use_gamma ) CALL fgamma(m_max, T, prim%F(1)) 00182 00183 ! Coulomb 00184 CALL fgamma(m_max,T,Fm) 00185 00186 DO i=1,m_max+1 00187 prim%F(i)=prim%F(i)*(potential_parameter%scale_coulomb+potential_parameter%scale_longrange)- & 00188 Fm(i)*potential_parameter%scale_longrange 00189 ENDDO 00190 00191 ! longrange 00192 omega2 = potential_parameter%omega**2 00193 omega_corr2 = omega2/(omega2+Rho) 00194 omega_corr = SQRT(omega_corr2) 00195 T = T*omega_corr2 00196 CALL fgamma(m_max,T,Fm) 00197 tmp = omega_corr 00198 DO i=1,m_max+1 00199 prim%F(i)= prim%F(i) + Fm(i)*tmp*potential_parameter%scale_longrange 00200 tmp = tmp * omega_corr2 00201 END DO 00202 factor = 2.0_dp*Pi*RhoInv 00203 00204 CASE(do_hfx_potential_gaussian) 00205 omega2 = potential_parameter%omega**2 00206 T = -omega2*T/(Rho+omega2) 00207 tmp = 1.0_dp 00208 DO i=1,m_max+1 00209 prim%F(i) = EXP(T) * tmp 00210 tmp = tmp * omega2/(Rho+omega2) 00211 END DO 00212 factor = (Pi/(Rho+omega2))**(1.5_dp) 00213 CASE(do_hfx_potential_mix_lg) 00214 omega2 = potential_parameter%omega**2 00215 omega_corr2 = omega2/(omega2+Rho) 00216 omega_corr = SQRT(omega_corr2) 00217 T = T*omega_corr2 00218 CALL fgamma(m_max,T,Fm) 00219 tmp = omega_corr * 2.0_dp*Pi*RhoInv * potential_parameter%scale_longrange 00220 DO i=1,m_max+1 00221 Fm(i)= Fm(i)*tmp 00222 tmp = tmp * omega_corr2 00223 END DO 00224 T = Rho*PQ2 00225 T = -omega2*T/(Rho+omega2) 00226 tmp = (Pi/(Rho+omega2))**(1.5_dp) * potential_parameter%scale_gaussian 00227 DO i=1,m_max+1 00228 prim%F(i) = EXP(T) * tmp + Fm(i) 00229 tmp = tmp * omega2/(Rho+omega2) 00230 END DO 00231 factor = 1.0_dp 00232 CASE(do_hfx_potential_id) 00233 prim%F(1) = (Pi*RhoInv)**(1.5_dp) 00234 DO i=2,m_max+1 00235 prim%F(i)= 0.0_dp 00236 END DO 00237 factor = 1.0_dp 00238 END SELECT 00239 00240 tmp = (Pi*ZetapEtaInv)**3 00241 factor = factor*S1234*SQRT(tmp) 00242 00243 DO i=1,m_max+1 00244 prim%F(i)=prim%F(i)*factor 00245 ENDDO 00246 prim%U(1:3,1) = P-A 00247 prim%U(1:3,3) = Q-C 00248 prim%U(1:3,5) = W-P 00249 prim%U(1:3,6) = W-Q 00250 prim%oo2z = 0.5_dp*ZetaInv 00251 prim%oo2n = 0.5_dp*EtaInv 00252 prim%oo2zn = 0.5_dp*ZetapEtaInv 00253 prim%poz = Rho*ZetaInv 00254 prim%pon = Rho*EtaInv 00255 prim%oo2p = 0.5_dp*RhoInv 00256 END SUBROUTINE build_quartet_data_screen 00257 00258 ! ***************************************************************************** 00264 SUBROUTINE build_deriv_data(prim, A, B, C, D,& 00265 Zeta_A, Zeta_B, Zeta_C, Zeta_D,& 00266 ZetaInv,EtaInv,ZetapEtaInv,Rho,RhoInv,& 00267 P,Q,W) 00268 00269 TYPE(prim_data) :: prim 00270 REAL(KIND=dp), INTENT(IN) :: A(3), B(3), C(3), D(3) 00271 REAL(dp), INTENT(IN) :: Zeta_A, Zeta_B, Zeta_C, 00272 Zeta_D, ZetaInv, EtaInv, 00273 ZetapEtaInv, Rho, RhoInv, 00274 P(3), Q(3), W(3) 00275 00276 prim%U(1:3,1) = P-A 00277 prim%U(1:3,2) = P-B 00278 prim%U(1:3,3) = Q-C 00279 prim%U(1:3,4) = Q-D 00280 prim%U(1:3,5) = W-P 00281 prim%U(1:3,6) = W-Q 00282 prim%twozeta_a = 2.0_dp*Zeta_A 00283 prim%twozeta_b = 2.0_dp*Zeta_B 00284 prim%twozeta_c = 2.0_dp*Zeta_C 00285 prim%twozeta_d = 2.0_dp*Zeta_D 00286 prim%oo2z = 0.5_dp*ZetaInv 00287 prim%oo2n = 0.5_dp*EtaInv 00288 prim%oo2zn = 0.5_dp*ZetapEtaInv 00289 prim%poz = Rho*ZetaInv 00290 prim%pon = Rho*EtaInv 00291 prim%oo2p = 0.5_dp*RhoInv 00292 00293 END SUBROUTINE build_deriv_data 00294 00295 ! ***************************************************************************** 00302 SUBROUTINE evaluate_deriv_eri(deriv, nproducts, pgf_product_list, & 00303 n_a,n_b,n_c,n_d,& 00304 ncoa,ncob,ncoc,ncod,& 00305 nsgfa, nsgfb, nsgfc, nsgfd,& 00306 primitives, max_contraction, tmp_max_all, & 00307 eps_schwarz, neris, & 00308 Zeta_A, Zeta_B, Zeta_C, Zeta_D, ZetaInv,EtaInv, & 00309 s_offset_a, s_offset_b, s_offset_c, s_offset_d,& 00310 nl_a, nl_b, nl_c, nl_d, nsoa, nsob, nsoc, nsod,& 00311 sphi_a, sphi_b, sphi_c, sphi_d,& 00312 work,work2,work_forces,buffer1,buffer2,primitives_tmp,& 00313 use_virial, work_virial, work2_virial, primitives_tmp_virial,& 00314 primitives_virial, cell, tmp_max_all_virial) 00315 00316 TYPE(lib_deriv) :: deriv 00317 INTEGER, INTENT(IN) :: nproducts 00318 TYPE(hfx_pgf_product_list), 00319 DIMENSION(nproducts) :: pgf_product_list 00320 INTEGER, INTENT(IN) :: n_a, n_b, n_c, n_d, ncoa, 00321 ncob, ncoc, ncod, nsgfa, 00322 nsgfb, nsgfc, nsgfd 00323 REAL(dp), DIMENSION(nsgfa, nsgfb, nsgfc,& nsgfd, 12) :: primitives 00324 REAL(dp) :: max_contraction, tmp_max_all, 00325 eps_schwarz 00326 INTEGER(int_8) :: neris 00327 REAL(dp), INTENT(IN) :: Zeta_A, Zeta_B, Zeta_C, 00328 Zeta_D, ZetaInv, EtaInv 00329 INTEGER :: s_offset_a, s_offset_b, 00330 s_offset_c, s_offset_d, nl_a, 00331 nl_b, nl_c, nl_d, nsoa, nsob, 00332 nsoc, nsod 00333 REAL(dp), DIMENSION(ncoa, nsoa*nl_a) :: sphi_a 00334 REAL(dp), DIMENSION(ncob, nsob*nl_b) :: sphi_b 00335 REAL(dp), DIMENSION(ncoc, nsoc*nl_c) :: sphi_c 00336 REAL(dp), DIMENSION(ncod, nsod*nl_d) :: sphi_d 00337 REAL(dp) :: work(ncoa*ncob*ncoc*ncod,12), work2(ncoa,ncob,ncoc,ncod,12), 00338 work_forces(ncoa*ncob*ncoc*ncod,12) 00339 REAL(dp), DIMENSION(ncoa*ncob*ncoc*ncod) :: buffer1, buffer2 00340 REAL(dp), DIMENSION(nsoa*nl_a, nsob*nl_b& , nsoc*nl_c, nsod*nl_d) :: primitives_tmp 00341 LOGICAL, INTENT(IN) :: use_virial 00342 REAL(dp) :: work_virial(ncoa*ncob*ncoc*ncod,12,3), 00343 work2_virial(ncoa,ncob,ncoc,ncod,12,3) 00344 REAL(dp), DIMENSION(nsoa*nl_a, nsob*nl_b& , nsoc*nl_c, nsod*nl_d) :: primitives_tmp_virial 00345 REAL(dp), DIMENSION(nsgfa, nsgfb, nsgfc,& nsgfd, 12, 3) :: primitives_virial 00346 TYPE(cell_type), POINTER :: cell 00347 REAL(dp) :: tmp_max_all_virial 00348 00349 INTEGER :: a_mysize(1), i, j, k, l, m, 00350 m_max, mysize, n, p1, p2, p3, 00351 p4, perm_case 00352 REAL(dp) :: A(3), AB(3), B(3), C(3), CD(3), D(3), P(3), Q(3), Rho, 00353 RhoInv, scoord(12), tmp_max, tmp_max_virial, W(3), ZetapEtaInv 00354 TYPE(prim_data), TARGET :: prim 00355 00356 m_max = n_a+n_b+n_c+n_d 00357 m_max = m_max + 1 00358 mysize = ncoa*ncob*ncoc*ncod 00359 a_mysize = mysize 00360 00361 work = 0.0_dp 00362 work2 = 0.0_dp 00363 00364 IF( use_virial ) THEN 00365 work_virial = 0.0_dp 00366 work2_virial = 0.0_dp 00367 END IF 00368 00369 00370 perm_case = 1 00371 IF(n_a<n_b) THEN 00372 perm_case = perm_case + 1 00373 END IF 00374 IF(n_c<n_d) THEN 00375 perm_case = perm_case + 2 00376 END IF 00377 IF(n_a+n_b > n_c+n_d) THEN 00378 perm_case = perm_case + 4 00379 END IF 00380 00381 SELECT CASE(perm_case) 00382 CASE(1) 00383 DO i = 1,nproducts 00384 A = pgf_product_list(i)%ra 00385 B = pgf_product_list(i)%rb 00386 C = pgf_product_list(i)%rc 00387 D = pgf_product_list(i)%rd 00388 Rho = pgf_product_list(i)%Rho 00389 RhoInv = pgf_product_list(i)%RhoInv 00390 P = pgf_product_list(i)%P 00391 Q = pgf_product_list(i)%Q 00392 W = pgf_product_list(i)%W 00393 AB = pgf_product_list(i)%AB 00394 CD = pgf_product_list(i)%CD 00395 ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv 00396 prim%F(1:m_max+1) = pgf_product_list(i)%Fm(1:m_max+1) 00397 00398 CALL build_deriv_data(prim, A, B, C, D,& 00399 Zeta_A, Zeta_B, Zeta_C, Zeta_D, & 00400 ZetaInv, EtaInv, ZetapEtaInv, Rho, RhoInv,& 00401 P, Q, W) 00402 deriv%AB=AB!A-B 00403 deriv%CD=CD!C-D 00404 CALL get_derivs(n_d, n_c, n_b, n_a, deriv, prim, work_forces, a_mysize) 00405 DO k=4,6 00406 DO j=1,mysize 00407 work_forces(j,k) = - 1.0_dp* (work_forces(j,k-3) + & 00408 work_forces(j,k+3) + & 00409 work_forces(j,k+6)) 00410 END DO 00411 END DO 00412 DO k=1,12 00413 DO j=1,mysize 00414 work(j,k) = work(j,k) + work_forces(j,k) 00415 END DO 00416 END DO 00417 neris = neris + 12 * mysize 00418 IF( use_virial ) THEN 00419 CALL real_to_scaled(scoord(1:3),A,cell) 00420 CALL real_to_scaled(scoord(4:6),B,cell) 00421 CALL real_to_scaled(scoord(7:9),C,cell) 00422 CALL real_to_scaled(scoord(10:12),D,cell) 00423 DO k=1,12 00424 DO j=1,mysize 00425 DO m=1,3 00426 work_virial(j,k,m) = work_virial(j,k,m) + work_forces(j,k)*scoord(INT((k-1)/3)*3+m) 00427 END DO 00428 END DO 00429 END DO 00430 END IF 00431 END DO 00432 00433 DO n = 1,12 00434 tmp_max = 0.0_dp 00435 DO i = 1,mysize 00436 tmp_max = MAX(tmp_max, ABS(work(i,n))) 00437 END DO 00438 tmp_max = tmp_max * max_contraction 00439 tmp_max_all = MAX(tmp_max_all, tmp_max) 00440 00441 DO i = 1,ncoa 00442 p1 = (i-1) *ncob 00443 DO j = 1,ncob 00444 p2 = (p1 + j-1)*ncoc 00445 DO k = 1,ncoc 00446 p3 = (p2 + k-1)*ncod 00447 DO l = 1,ncod 00448 p4 = p3 + l 00449 work2(i,j,k,l,full_perm1(n)) = work(p4,n) 00450 END DO 00451 END DO 00452 END DO 00453 END DO 00454 END DO 00455 IF( use_virial ) THEN 00456 DO n=1,12 00457 tmp_max_virial = 0.0_dp 00458 DO i = 1,mysize 00459 tmp_max_virial = MAX(tmp_max_virial, & 00460 ABS(work_virial(i,n,1)), ABS(work_virial(i,n,2)), ABS(work_virial(i,n,3))) 00461 END DO 00462 tmp_max_virial = tmp_max_virial * max_contraction 00463 tmp_max_all_virial = MAX(tmp_max_all_virial,tmp_max_virial) 00464 00465 DO i = 1,ncoa 00466 p1 = (i-1) *ncob 00467 DO j = 1,ncob 00468 p2 = (p1 + j-1)*ncoc 00469 DO k = 1,ncoc 00470 p3 = (p2 + k-1)*ncod 00471 DO l = 1,ncod 00472 p4 = p3 + l 00473 work2_virial(i,j,k,l,full_perm1(n),1:3) = work_virial(p4,n,1:3) 00474 END DO 00475 END DO 00476 END DO 00477 END DO 00478 END DO 00479 END IF 00480 CASE(2) 00481 DO i = 1,nproducts 00482 A = pgf_product_list(i)%ra 00483 B = pgf_product_list(i)%rb 00484 C = pgf_product_list(i)%rc 00485 D = pgf_product_list(i)%rd 00486 Rho = pgf_product_list(i)%Rho 00487 RhoInv = pgf_product_list(i)%RhoInv 00488 P = pgf_product_list(i)%P 00489 Q = pgf_product_list(i)%Q 00490 W = pgf_product_list(i)%W 00491 AB = pgf_product_list(i)%AB 00492 CD = pgf_product_list(i)%CD 00493 ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv 00494 prim%F(1:m_max+1) = pgf_product_list(i)%Fm(1:m_max+1) 00495 00496 CALL build_deriv_data(prim, B, A, C, D,& 00497 Zeta_B, Zeta_A, Zeta_C, Zeta_D,& 00498 ZetaInv, EtaInv, ZetapEtaInv, Rho, RhoInv,& 00499 P, Q, W) 00500 deriv%AB=-AB!B-A 00501 deriv%CD=CD!C-D 00502 CALL get_derivs(n_d, n_c, n_a, n_b, deriv, prim, work_forces, a_mysize) 00503 DO k=4,6 00504 DO j=1,mysize 00505 work_forces(j,k) = - 1.0_dp* (work_forces(j,k-3) + & 00506 work_forces(j,k+3) + & 00507 work_forces(j,k+6) ) 00508 ENDDO 00509 END DO 00510 DO k=1,12 00511 DO j=1,mysize 00512 work(j,k) = work(j,k) + work_forces(j,k) 00513 END DO 00514 END DO 00515 neris = neris + 12 * mysize 00516 IF( use_virial ) THEN 00517 CALL real_to_scaled(scoord(1:3),B,cell) 00518 CALL real_to_scaled(scoord(4:6),A,cell) 00519 CALL real_to_scaled(scoord(7:9),C,cell) 00520 CALL real_to_scaled(scoord(10:12),D,cell) 00521 DO k=1,12 00522 DO j=1,mysize 00523 DO m=1,3 00524 work_virial(j,k,m) = work_virial(j,k,m) + work_forces(j,k)*scoord(INT((k-1)/3)*3+m) 00525 END DO 00526 END DO 00527 END DO 00528 END IF 00529 00530 END DO 00531 DO n = 1,12 00532 tmp_max = 0.0_dp 00533 DO i = 1,mysize 00534 tmp_max = MAX(tmp_max, ABS(work(i,n))) 00535 END DO 00536 tmp_max = tmp_max * max_contraction 00537 tmp_max_all = MAX(tmp_max_all, tmp_max) 00538 00539 DO j = 1,ncob 00540 p1 = (j-1)*ncoa 00541 DO i = 1,ncoa 00542 p2 = (p1 + i-1)*ncoc 00543 DO k = 1,ncoc 00544 p3 = (p2 + k-1)*ncod 00545 DO l = 1,ncod 00546 p4 = p3 + l 00547 work2(i,j,k,l,full_perm2(n)) = work(p4,n) 00548 END DO 00549 END DO 00550 END DO 00551 END DO 00552 END DO 00553 IF( use_virial ) THEN 00554 DO n = 1,12 00555 tmp_max_virial = 0.0_dp 00556 DO i = 1,mysize 00557 tmp_max_virial = MAX(tmp_max_virial, & 00558 ABS(work_virial(i,n,1)), ABS(work_virial(i,n,2)), ABS(work_virial(i,n,3))) 00559 END DO 00560 tmp_max_virial = tmp_max_virial * max_contraction 00561 tmp_max_all_virial = MAX(tmp_max_all_virial,tmp_max_virial) 00562 00563 DO j = 1,ncob 00564 p1 = (j-1)*ncoa 00565 DO i = 1,ncoa 00566 p2 = (p1 + i-1)*ncoc 00567 DO k = 1,ncoc 00568 p3 = (p2 + k-1)*ncod 00569 DO l = 1,ncod 00570 p4 = p3 + l 00571 work2_virial(i,j,k,l,full_perm2(n),1:3) = work_virial(p4,n,1:3) 00572 END DO 00573 END DO 00574 END DO 00575 END DO 00576 END DO 00577 END IF 00578 CASE(3) 00579 DO i = 1,nproducts 00580 A = pgf_product_list(i)%ra 00581 B = pgf_product_list(i)%rb 00582 C = pgf_product_list(i)%rc 00583 D = pgf_product_list(i)%rd 00584 Rho = pgf_product_list(i)%Rho 00585 RhoInv = pgf_product_list(i)%RhoInv 00586 P = pgf_product_list(i)%P 00587 Q = pgf_product_list(i)%Q 00588 W = pgf_product_list(i)%W 00589 AB = pgf_product_list(i)%AB 00590 CD = pgf_product_list(i)%CD 00591 ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv 00592 prim%F(1:m_max+1) = pgf_product_list(i)%Fm(1:m_max+1) 00593 00594 CALL build_deriv_data(prim, A, B, D, C,& 00595 Zeta_A, Zeta_B, Zeta_D, Zeta_C,& 00596 ZetaInv, EtaInv, ZetapEtaInv, Rho, RhoInv,& 00597 P, Q, W) 00598 deriv%AB=AB!A-B 00599 deriv%CD=-CD!D-C 00600 CALL get_derivs(n_c, n_d, n_b, n_a, deriv, prim, work_forces, a_mysize) 00601 DO k=4,6 00602 DO j=1,mysize 00603 work_forces(j,k) = - 1.0_dp* (work_forces(j,k-3) + & 00604 work_forces(j,k+3) + & 00605 work_forces(j,k+6) ) 00606 END DO 00607 END DO 00608 DO k=1,12 00609 DO j=1,mysize 00610 work(j,k) = work(j,k) + work_forces(j,k) 00611 END DO 00612 END DO 00613 neris = neris + 12 * mysize 00614 IF( use_virial ) THEN 00615 CALL real_to_scaled(scoord(1:3),A,cell) 00616 CALL real_to_scaled(scoord(4:6),B,cell) 00617 CALL real_to_scaled(scoord(7:9),D,cell) 00618 CALL real_to_scaled(scoord(10:12),C,cell) 00619 DO k=1,12 00620 DO j=1,mysize 00621 DO m=1,3 00622 work_virial(j,k,m) = work_virial(j,k,m) + work_forces(j,k)*scoord(INT((k-1)/3)*3+m) 00623 END DO 00624 END DO 00625 END DO 00626 END IF 00627 00628 END DO 00629 DO n = 1,12 00630 tmp_max = 0.0_dp 00631 DO i = 1,mysize 00632 tmp_max = MAX(tmp_max, ABS(work(i,n))) 00633 END DO 00634 tmp_max = tmp_max * max_contraction 00635 tmp_max_all = MAX(tmp_max_all, tmp_max) 00636 00637 DO i = 1,ncoa 00638 p1 = (i-1)*ncob 00639 DO j = 1,ncob 00640 p2 = (p1 + j-1)*ncod 00641 DO l = 1,ncod 00642 p3 = (p2 + l-1) * ncoc 00643 DO k = 1,ncoc 00644 p4 = p3+k 00645 work2(i,j,k,l,full_perm3(n)) = work(p4,n) 00646 END DO 00647 END DO 00648 END DO 00649 END DO 00650 END DO 00651 IF( use_virial ) THEN 00652 DO n = 1,12 00653 tmp_max_virial = 0.0_dp 00654 DO i = 1,mysize 00655 tmp_max_virial = MAX(tmp_max_virial, & 00656 ABS(work_virial(i,n,1)), ABS(work_virial(i,n,2)), ABS(work_virial(i,n,3))) 00657 END DO 00658 tmp_max_virial = tmp_max_virial * max_contraction 00659 tmp_max_all_virial = MAX(tmp_max_all_virial,tmp_max_virial) 00660 00661 DO i = 1,ncoa 00662 p1 = (i-1)*ncob 00663 DO j = 1,ncob 00664 p2 = (p1 + j-1)*ncod 00665 DO l = 1,ncod 00666 p3 = (p2 + l-1) * ncoc 00667 DO k = 1,ncoc 00668 p4 = p3+k 00669 work2_virial(i,j,k,l,full_perm3(n),1:3) = work_virial(p4,n,1:3) 00670 END DO 00671 END DO 00672 END DO 00673 END DO 00674 END DO 00675 END IF 00676 CASE(4) 00677 DO i = 1,nproducts 00678 A = pgf_product_list(i)%ra 00679 B = pgf_product_list(i)%rb 00680 C = pgf_product_list(i)%rc 00681 D = pgf_product_list(i)%rd 00682 Rho = pgf_product_list(i)%Rho 00683 RhoInv = pgf_product_list(i)%RhoInv 00684 P = pgf_product_list(i)%P 00685 Q = pgf_product_list(i)%Q 00686 W = pgf_product_list(i)%W 00687 AB = pgf_product_list(i)%AB 00688 CD = pgf_product_list(i)%CD 00689 ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv 00690 prim%F(1:m_max+1) = pgf_product_list(i)%Fm(1:m_max+1) 00691 00692 CALL build_deriv_data(prim, B, A, D, C,& 00693 Zeta_B, Zeta_A, Zeta_D, Zeta_C,& 00694 ZetaInv, EtaInv, ZetapEtaInv, Rho, RhoInv,& 00695 P, Q, W) 00696 deriv%AB=-AB!B-A 00697 deriv%CD=-CD!D-C 00698 CALL get_derivs(n_c, n_d, n_a, n_b, deriv, prim, work_forces, a_mysize) 00699 DO k=4,6 00700 DO j=1,mysize 00701 work_forces(j,k) = - 1.0_dp* (work_forces(j,k-3) + & 00702 work_forces(j,k+3) + & 00703 work_forces(j,k+6) ) 00704 END DO 00705 END DO 00706 DO k=1,12 00707 DO j=1,mysize 00708 work(j,k) = work(j,k) + work_forces(j,k) 00709 END DO 00710 END DO 00711 neris = neris + 12 * mysize 00712 IF( use_virial ) THEN 00713 CALL real_to_scaled(scoord(1:3),B,cell) 00714 CALL real_to_scaled(scoord(4:6),A,cell) 00715 CALL real_to_scaled(scoord(7:9),D,cell) 00716 CALL real_to_scaled(scoord(10:12),C,cell) 00717 DO k=1,12 00718 DO j=1,mysize 00719 DO m=1,3 00720 work_virial(j,k,m) = work_virial(j,k,m) + work_forces(j,k)*scoord(INT((k-1)/3)*3+m) 00721 END DO 00722 END DO 00723 END DO 00724 END IF 00725 00726 END DO 00727 DO n = 1,12 00728 tmp_max = 0.0_dp 00729 DO i = 1,mysize 00730 tmp_max = MAX(tmp_max, ABS(work(i,n))) 00731 END DO 00732 tmp_max = tmp_max * max_contraction 00733 tmp_max_all = MAX(tmp_max_all, tmp_max) 00734 00735 DO j = 1,ncob 00736 p1 = (j-1)*ncoa 00737 DO i = 1,ncoa 00738 p2 = (p1 + i-1)*ncod 00739 DO l = 1,ncod 00740 p3 = (p2 + l-1)*ncoc 00741 DO k = 1,ncoc 00742 p4 = p3 + k 00743 work2(i,j,k,l,full_perm4(n)) = work(p4,n) 00744 END DO 00745 END DO 00746 END DO 00747 END DO 00748 END DO 00749 IF( use_virial ) THEN 00750 DO n = 1,12 00751 tmp_max_virial = 0.0_dp 00752 DO i = 1,mysize 00753 tmp_max_virial = MAX(tmp_max_virial, & 00754 ABS(work_virial(i,n,1)), ABS(work_virial(i,n,2)), ABS(work_virial(i,n,3))) 00755 END DO 00756 tmp_max_virial = tmp_max_virial * max_contraction 00757 tmp_max_all_virial = MAX(tmp_max_all_virial,tmp_max_virial) 00758 00759 DO j = 1,ncob 00760 p1 = (j-1)*ncoa 00761 DO i = 1,ncoa 00762 p2 = (p1 + i-1)*ncod 00763 DO l = 1,ncod 00764 p3 = (p2 + l-1)*ncoc 00765 DO k = 1,ncoc 00766 p4 = p3 + k 00767 work2_virial(i,j,k,l,full_perm4(n),1:3) = work_virial(p4,n,1:3) 00768 END DO 00769 END DO 00770 END DO 00771 END DO 00772 END DO 00773 END IF 00774 CASE(5) 00775 DO i = 1,nproducts 00776 A = pgf_product_list(i)%ra 00777 B = pgf_product_list(i)%rb 00778 C = pgf_product_list(i)%rc 00779 D = pgf_product_list(i)%rd 00780 Rho = pgf_product_list(i)%Rho 00781 RhoInv = pgf_product_list(i)%RhoInv 00782 P = pgf_product_list(i)%P 00783 Q = pgf_product_list(i)%Q 00784 W = pgf_product_list(i)%W 00785 AB = pgf_product_list(i)%AB 00786 CD = pgf_product_list(i)%CD 00787 ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv 00788 prim%F(1:m_max+1) = pgf_product_list(i)%Fm(1:m_max+1) 00789 00790 CALL build_deriv_data(prim, C, D, A, B,& 00791 Zeta_C ,Zeta_D, Zeta_A, Zeta_B,& 00792 EtaInv, ZetaInv, ZetapEtaInv, Rho, RhoInv,& 00793 Q, P, W) 00794 deriv%AB=CD!C-D 00795 deriv%CD=AB!A-B 00796 CALL get_derivs(n_b, n_a, n_d, n_c, deriv, prim, work_forces, a_mysize) 00797 DO k=4,6 00798 DO j=1,mysize 00799 work_forces(j,k) = - 1.0_dp* (work_forces(j,k-3) + & 00800 work_forces(j,k+3) + & 00801 work_forces(j,k+6) ) 00802 END DO 00803 END DO 00804 DO k=1,12 00805 DO j=1,mysize 00806 work(j,k) = work(j,k) + work_forces(j,k) 00807 END DO 00808 END DO 00809 neris = neris + 12 * mysize 00810 IF( use_virial ) THEN 00811 CALL real_to_scaled(scoord(1:3),C,cell) 00812 CALL real_to_scaled(scoord(4:6),D,cell) 00813 CALL real_to_scaled(scoord(7:9),A,cell) 00814 CALL real_to_scaled(scoord(10:12),B,cell) 00815 DO k=1,12 00816 DO j=1,mysize 00817 DO m=1,3 00818 work_virial(j,k,m) = work_virial(j,k,m) + work_forces(j,k)*scoord(INT((k-1)/3)*3+m) 00819 END DO 00820 END DO 00821 END DO 00822 END IF 00823 00824 END DO 00825 DO n = 1,12 00826 tmp_max = 0.0_dp 00827 DO i = 1,mysize 00828 tmp_max = MAX(tmp_max, ABS(work(i,n))) 00829 END DO 00830 tmp_max = tmp_max * max_contraction 00831 tmp_max_all = MAX(tmp_max_all, tmp_max) 00832 00833 DO k = 1,ncoc 00834 p1 = (k-1)*ncod 00835 DO l = 1,ncod 00836 p2 = (p1 + l-1)*ncoa 00837 DO i = 1,ncoa 00838 p3 = (p2 + i-1)*ncob 00839 DO j = 1,ncob 00840 p4 = p3+j 00841 work2(i,j,k,l,full_perm5(n)) = work(p4,n) 00842 END DO 00843 END DO 00844 END DO 00845 END DO 00846 END DO 00847 IF( use_virial ) THEN 00848 DO n = 1,12 00849 tmp_max_virial = 0.0_dp 00850 DO i = 1,mysize 00851 tmp_max_virial = MAX(tmp_max_virial, & 00852 ABS(work_virial(i,n,1)), ABS(work_virial(i,n,2)), ABS(work_virial(i,n,3))) 00853 END DO 00854 tmp_max_virial = tmp_max_virial * max_contraction 00855 tmp_max_all_virial = MAX(tmp_max_all_virial,tmp_max_virial) 00856 00857 DO k = 1,ncoc 00858 p1 = (k-1)*ncod 00859 DO l = 1,ncod 00860 p2 = (p1 + l-1)*ncoa 00861 DO i = 1,ncoa 00862 p3 = (p2 + i-1)*ncob 00863 DO j = 1,ncob 00864 p4 = p3+j 00865 work2_virial(i,j,k,l,full_perm5(n),1:3) = work_virial(p4,n,1:3) 00866 END DO 00867 END DO 00868 END DO 00869 END DO 00870 END DO 00871 END IF 00872 CASE(6) 00873 DO i = 1,nproducts 00874 A = pgf_product_list(i)%ra 00875 B = pgf_product_list(i)%rb 00876 C = pgf_product_list(i)%rc 00877 D = pgf_product_list(i)%rd 00878 Rho = pgf_product_list(i)%Rho 00879 RhoInv = pgf_product_list(i)%RhoInv 00880 P = pgf_product_list(i)%P 00881 Q = pgf_product_list(i)%Q 00882 W = pgf_product_list(i)%W 00883 AB = pgf_product_list(i)%AB 00884 CD = pgf_product_list(i)%CD 00885 ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv 00886 prim%F(1:m_max+1) = pgf_product_list(i)%Fm(1:m_max+1) 00887 00888 CALL build_deriv_data(prim, C, D, B, A, & 00889 Zeta_C, Zeta_D, Zeta_B, Zeta_A,& 00890 EtaInv, ZetaInv, ZetapEtaInv, Rho, RhoInv,& 00891 Q, P, W) 00892 deriv%AB=CD!C-D 00893 deriv%CD=-AB!B-A 00894 CALL get_derivs(n_a, n_b, n_d, n_c, deriv, prim, work_forces, a_mysize) 00895 DO k=4,6 00896 DO j=1,mysize 00897 work_forces(j,k) = - 1.0_dp* (work_forces(j,k-3) + & 00898 work_forces(j,k+3) + & 00899 work_forces(j,k+6) ) 00900 END DO 00901 END DO 00902 DO k=1,12 00903 DO j=1,mysize 00904 work(j,k) = work(j,k) + work_forces(j,k) 00905 END DO 00906 END DO 00907 neris = neris + 12 * mysize 00908 IF( use_virial ) THEN 00909 CALL real_to_scaled(scoord(1:3),C,cell) 00910 CALL real_to_scaled(scoord(4:6),D,cell) 00911 CALL real_to_scaled(scoord(7:9),B,cell) 00912 CALL real_to_scaled(scoord(10:12),A,cell) 00913 DO k=1,12 00914 DO j=1,mysize 00915 DO m=1,3 00916 work_virial(j,k,m) = work_virial(j,k,m) + work_forces(j,k)*scoord(INT((k-1)/3)*3+m) 00917 END DO 00918 END DO 00919 END DO 00920 END IF 00921 00922 END DO 00923 DO n = 1,12 00924 tmp_max = 0.0_dp 00925 DO i = 1,mysize 00926 tmp_max = MAX(tmp_max, ABS(work(i,n))) 00927 END DO 00928 tmp_max = tmp_max * max_contraction 00929 tmp_max_all = MAX(tmp_max_all, tmp_max) 00930 00931 DO k = 1,ncoc 00932 p1 = (k-1)*ncod 00933 DO l = 1,ncod 00934 p2 = (p1 + l-1)*ncob 00935 DO j = 1,ncob 00936 p3 = (p2 + j-1)*ncoa 00937 DO i = 1,ncoa 00938 p4 = p3+i 00939 work2(i,j,k,l,full_perm6(n)) = work(p4,n) 00940 END DO 00941 END DO 00942 END DO 00943 END DO 00944 END DO 00945 IF( use_virial ) THEN 00946 DO n = 1,12 00947 tmp_max_virial = 0.0_dp 00948 DO i = 1,mysize 00949 tmp_max_virial = MAX(tmp_max_virial, & 00950 ABS(work_virial(i,n,1)), ABS(work_virial(i,n,2)), ABS(work_virial(i,n,3))) 00951 END DO 00952 tmp_max_virial = tmp_max_virial * max_contraction 00953 tmp_max_all_virial = MAX(tmp_max_all_virial,tmp_max_virial) 00954 00955 DO k = 1,ncoc 00956 p1 = (k-1)*ncod 00957 DO l = 1,ncod 00958 p2 = (p1 + l-1)*ncob 00959 DO j = 1,ncob 00960 p3 = (p2 + j-1)*ncoa 00961 DO i = 1,ncoa 00962 p4 = p3+i 00963 work2_virial(i,j,k,l,full_perm6(n),1:3) = work_virial(p4,n,1:3) 00964 END DO 00965 END DO 00966 END DO 00967 END DO 00968 END DO 00969 END IF 00970 CASE(7) 00971 DO i = 1,nproducts 00972 A = pgf_product_list(i)%ra 00973 B = pgf_product_list(i)%rb 00974 C = pgf_product_list(i)%rc 00975 D = pgf_product_list(i)%rd 00976 Rho = pgf_product_list(i)%Rho 00977 RhoInv = pgf_product_list(i)%RhoInv 00978 P = pgf_product_list(i)%P 00979 Q = pgf_product_list(i)%Q 00980 W = pgf_product_list(i)%W 00981 AB = pgf_product_list(i)%AB 00982 CD = pgf_product_list(i)%CD 00983 ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv 00984 prim%F(1:m_max+1) = pgf_product_list(i)%Fm(1:m_max+1) 00985 00986 CALL build_deriv_data(prim, D, C, A, B,& 00987 Zeta_D, Zeta_C, Zeta_A, Zeta_B,& 00988 EtaInv, ZetaInv, ZetapEtaInv, Rho, RhoInv,& 00989 Q, P, W) 00990 deriv%AB=-CD!D-C 00991 deriv%CD=AB!A-B 00992 CALL get_derivs(n_b, n_a, n_c, n_d, deriv, prim, work_forces, a_mysize) 00993 DO k=4,6 00994 DO j=1,mysize 00995 work_forces(j,k) = - 1.0_dp* (work_forces(j,k-3) + & 00996 work_forces(j,k+3) + & 00997 work_forces(j,k+6) ) 00998 END DO 00999 END DO 01000 DO k=1,12 01001 DO j=1,mysize 01002 work(j,k) = work(j,k) + work_forces(j,k) 01003 END DO 01004 END DO 01005 neris = neris + 12 * mysize 01006 IF( use_virial ) THEN 01007 CALL real_to_scaled(scoord(1:3),D,cell) 01008 CALL real_to_scaled(scoord(4:6),C,cell) 01009 CALL real_to_scaled(scoord(7:9),A,cell) 01010 CALL real_to_scaled(scoord(10:12),B,cell) 01011 DO k=1,12 01012 DO j=1,mysize 01013 DO m=1,3 01014 work_virial(j,k,m) = work_virial(j,k,m) + work_forces(j,k)*scoord(INT((k-1)/3)*3+m) 01015 END DO 01016 END DO 01017 END DO 01018 END IF 01019 01020 END DO 01021 DO n = 1,12 01022 tmp_max = 0.0_dp 01023 DO i = 1,mysize 01024 tmp_max = MAX(tmp_max, ABS(work(i,n))) 01025 END DO 01026 tmp_max = tmp_max * max_contraction 01027 tmp_max_all = MAX(tmp_max_all, tmp_max) 01028 01029 DO l = 1,ncod 01030 p1 = (l-1)*ncoc 01031 DO k = 1,ncoc 01032 p2 = (p1 + k-1) * ncoa 01033 DO i = 1,ncoa 01034 p3 = (p2 + i-1) *ncob 01035 DO j = 1,ncob 01036 p4 = p3+j 01037 work2(i,j,k,l,full_perm7(n)) = work(p4,n) 01038 END DO 01039 END DO 01040 END DO 01041 END DO 01042 END DO 01043 IF( use_virial ) THEN 01044 DO n = 1,12 01045 tmp_max_virial = 0.0_dp 01046 DO i = 1,mysize 01047 tmp_max_virial = MAX(tmp_max_virial, & 01048 ABS(work_virial(i,n,1)), ABS(work_virial(i,n,2)), ABS(work_virial(i,n,3))) 01049 END DO 01050 tmp_max_virial = tmp_max_virial * max_contraction 01051 tmp_max_all_virial = MAX(tmp_max_all_virial,tmp_max_virial) 01052 01053 DO l = 1,ncod 01054 p1 = (l-1)*ncoc 01055 DO k = 1,ncoc 01056 p2 = (p1 + k-1) * ncoa 01057 DO i = 1,ncoa 01058 p3 = (p2 + i-1) *ncob 01059 DO j = 1,ncob 01060 p4 = p3+j 01061 work2_virial(i,j,k,l,full_perm7(n),1:3) = work_virial(p4,n,1:3) 01062 END DO 01063 END DO 01064 END DO 01065 END DO 01066 END DO 01067 END IF 01068 CASE(8) 01069 DO i = 1,nproducts 01070 A = pgf_product_list(i)%ra 01071 B = pgf_product_list(i)%rb 01072 C = pgf_product_list(i)%rc 01073 D = pgf_product_list(i)%rd 01074 Rho = pgf_product_list(i)%Rho 01075 RhoInv = pgf_product_list(i)%RhoInv 01076 P = pgf_product_list(i)%P 01077 Q = pgf_product_list(i)%Q 01078 W = pgf_product_list(i)%W 01079 AB = pgf_product_list(i)%AB 01080 CD = pgf_product_list(i)%CD 01081 ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv 01082 prim%F(1:m_max+1) = pgf_product_list(i)%Fm(1:m_max+1) 01083 01084 CALL build_deriv_data(prim, D, C, B, A,& 01085 Zeta_D, Zeta_C, Zeta_B, Zeta_A,& 01086 EtaInv, ZetaInv, ZetapEtaInv, Rho, RhoInv,& 01087 Q, P, W) 01088 deriv%AB=-CD!D-C 01089 deriv%CD=-AB!B-A 01090 CALL get_derivs(n_a, n_b, n_c, n_d, deriv, prim, work_forces, a_mysize) 01091 DO k=4,6 01092 DO j=1,mysize 01093 work_forces(j,k) = - 1.0_dp* (work_forces(j,k-3) + & 01094 work_forces(j,k+3) + & 01095 work_forces(j,k+6) ) 01096 END DO 01097 END DO 01098 DO k=1,12 01099 DO j=1,mysize 01100 work(j,k) = work(j,k) + work_forces(j,k) 01101 END DO 01102 END DO 01103 neris = neris + 12 * mysize 01104 IF( use_virial ) THEN 01105 CALL real_to_scaled(scoord(1:3),D,cell) 01106 CALL real_to_scaled(scoord(4:6),C,cell) 01107 CALL real_to_scaled(scoord(7:9),B,cell) 01108 CALL real_to_scaled(scoord(10:12),A,cell) 01109 DO k=1,12 01110 DO j=1,mysize 01111 DO m=1,3 01112 work_virial(j,k,m) = work_virial(j,k,m) + work_forces(j,k)*scoord(INT((k-1)/3)*3+m) 01113 END DO 01114 END DO 01115 END DO 01116 END IF 01117 01118 END DO 01119 DO n = 1,12 01120 tmp_max = 0.0_dp 01121 DO i = 1,mysize 01122 tmp_max = MAX(tmp_max, ABS(work(i,n))) 01123 END DO 01124 tmp_max = tmp_max * max_contraction 01125 tmp_max_all = MAX(tmp_max_all, tmp_max) 01126 01127 DO l = 1,ncod 01128 p1 = (l-1)*ncoc 01129 DO k = 1,ncoc 01130 p2 = (p1 + k-1) * ncob 01131 DO j = 1,ncob 01132 p3 = (p2 + j-1) * ncoa 01133 DO i = 1,ncoa 01134 p4 = p3 + i 01135 work2(i,j,k,l,full_perm8(n)) = work(p4,n) 01136 END DO 01137 END DO 01138 END DO 01139 END DO 01140 END DO 01141 IF( use_virial ) THEN 01142 DO n = 1,12 01143 tmp_max_virial = 0.0_dp 01144 DO i = 1,mysize 01145 tmp_max_virial = MAX(tmp_max_virial, & 01146 ABS(work_virial(i,n,1)), ABS(work_virial(i,n,2)), ABS(work_virial(i,n,3))) 01147 END DO 01148 tmp_max_virial = tmp_max_virial * max_contraction 01149 tmp_max_all_virial = MAX(tmp_max_all_virial,tmp_max_virial) 01150 01151 DO l = 1,ncod 01152 p1 = (l-1)*ncoc 01153 DO k = 1,ncoc 01154 p2 = (p1 + k-1) * ncob 01155 DO j = 1,ncob 01156 p3 = (p2 + j-1) * ncoa 01157 DO i = 1,ncoa 01158 p4 = p3 + i 01159 work2_virial(i,j,k,l,full_perm8(n),1:3) = work_virial(p4,n,1:3) 01160 END DO 01161 END DO 01162 END DO 01163 END DO 01164 END DO 01165 END IF 01166 END SELECT 01167 01168 IF( .NOT. use_virial) THEN 01169 IF( tmp_max_all < eps_schwarz ) RETURN 01170 END IF 01171 01172 IF( tmp_max_all >= eps_schwarz ) THEN 01173 DO i = 1,12 01174 primitives_tmp(:,:,:,:) = 0.0_dp 01175 CALL contract(ncoa, ncob, ncoc, ncod, nsoa, nsob, nsoc, nsod,& 01176 n_a, n_b, n_c, n_d, nl_a, nl_b, nl_c, nl_d, work2(1,1,1,1,i),& 01177 sphi_a, & 01178 sphi_b, & 01179 sphi_c, & 01180 sphi_d, & 01181 primitives_tmp(1,1,1,1), & 01182 buffer1, buffer2) 01183 01184 primitives(s_offset_a+1:s_offset_a+nsoa*nl_a,& 01185 s_offset_b+1:s_offset_b+nsob*nl_b,& 01186 s_offset_c+1:s_offset_c+nsoc*nl_c,& 01187 s_offset_d+1:s_offset_d+nsod*nl_d,i) = & 01188 primitives(s_offset_a+1:s_offset_a+nsoa*nl_a,& 01189 s_offset_b+1:s_offset_b+nsob*nl_b,& 01190 s_offset_c+1:s_offset_c+nsoc*nl_c,& 01191 s_offset_d+1:s_offset_d+nsod*nl_d,i) + primitives_tmp(:,:,:,:) 01192 END DO 01193 END IF 01194 01195 IF( use_virial .AND. tmp_max_all_virial >= eps_schwarz ) THEN 01196 DO i = 1,12 01197 DO m=1,3 01198 primitives_tmp_virial(:,:,:,:) = 0.0_dp 01199 CALL contract(ncoa, ncob, ncoc, ncod, nsoa, nsob, nsoc, nsod,& 01200 n_a, n_b, n_c, n_d, nl_a, nl_b, nl_c, nl_d, work2_virial(1,1,1,1,i,m),& 01201 sphi_a, & 01202 sphi_b, & 01203 sphi_c, & 01204 sphi_d, & 01205 primitives_tmp_virial(1,1,1,1), & 01206 buffer1, buffer2) 01207 01208 primitives_virial(s_offset_a+1:s_offset_a+nsoa*nl_a,& 01209 s_offset_b+1:s_offset_b+nsob*nl_b,& 01210 s_offset_c+1:s_offset_c+nsoc*nl_c,& 01211 s_offset_d+1:s_offset_d+nsod*nl_d,i,m) = & 01212 primitives_virial(s_offset_a+1:s_offset_a+nsoa*nl_a,& 01213 s_offset_b+1:s_offset_b+nsob*nl_b,& 01214 s_offset_c+1:s_offset_c+nsoc*nl_c,& 01215 s_offset_d+1:s_offset_d+nsod*nl_d,i,m) + primitives_tmp_virial(:,:,:,:) 01216 END DO 01217 END DO 01218 END IF 01219 01220 END SUBROUTINE evaluate_deriv_eri 01221 01222 ! ***************************************************************************** 01229 SUBROUTINE evaluate_eri_screen(lib,A,B,C,D,Zeta_A,Zeta_B,Zeta_C,Zeta_D,& 01230 n_a,n_b,n_c,n_d,& 01231 max_val, potential_parameter, R1, R2,& 01232 p_work) 01233 01234 TYPE(lib_int) :: lib 01235 REAL(dp), INTENT(IN) :: A(3), B(3), C(3), D(3), 01236 Zeta_A, Zeta_B, Zeta_C, Zeta_D 01237 INTEGER, INTENT(IN) :: n_a, n_b, n_c, n_d 01238 REAL(dp), INTENT(INOUT) :: max_val 01239 TYPE(hfx_potential_type) :: potential_parameter 01240 REAL(dp) :: R1, R2 01241 REAL(dp), DIMENSION(:), POINTER :: p_work 01242 01243 INTEGER :: a_mysize(1), i, m_max, 01244 mysize, perm_case 01245 TYPE(prim_data), TARGET :: prim 01246 01247 m_max = n_a+n_b+n_c+n_d 01248 mysize = nco(n_a)*nco(n_b)*nco(n_c)*nco(n_d) 01249 a_mysize = mysize 01250 01251 IF(m_max/=0) THEN 01252 perm_case = 1 01253 IF(n_a<n_b) THEN 01254 perm_case = perm_case + 1 01255 END IF 01256 IF(n_c<n_d) THEN 01257 perm_case = perm_case + 2 01258 END IF 01259 IF(n_a+n_b > n_c+n_d) THEN 01260 perm_case = perm_case + 4 01261 END IF 01262 01263 SELECT CASE(perm_case) 01264 CASE(1) 01265 CALL build_quartet_data_screen(A,B,C,D,Zeta_A, Zeta_B, Zeta_C, Zeta_D, m_max,& 01266 potential_parameter, prim, R1, R2) 01267 lib%AB=A-B 01268 lib%CD=C-D 01269 CALL get_eris(n_d, n_c, n_b, n_a, lib, prim, p_work, a_mysize) 01270 DO i=1,mysize 01271 max_val = MAX(max_val, ABS(p_work(i))) 01272 END DO 01273 CASE(2) 01274 CALL build_quartet_data_screen(B,A,C,D,Zeta_B, Zeta_A, Zeta_C, Zeta_D, m_max,& 01275 potential_parameter, prim, R1,R2) 01276 lib%AB=B-A 01277 lib%CD=C-D 01278 CALL get_eris(n_d, n_c, n_a, n_b, lib, prim, p_work, a_mysize) 01279 DO i=1,mysize 01280 max_val = MAX(max_val, ABS(p_work(i))) 01281 END DO 01282 CASE(3) 01283 CALL build_quartet_data_screen(A,B,D,C,Zeta_A, Zeta_B, Zeta_D, Zeta_C, m_max,& 01284 potential_parameter, prim, R1,R2) 01285 lib%AB=A-B 01286 lib%CD=D-C 01287 CALL get_eris(n_c, n_d, n_b, n_a, lib, prim, p_work, a_mysize) 01288 DO i=1,mysize 01289 max_val = MAX(max_val, ABS(p_work(i))) 01290 END DO 01291 CASE(4) 01292 CALL build_quartet_data_screen(B,A,D,C,Zeta_B, Zeta_A, Zeta_D, Zeta_C, m_max,& 01293 potential_parameter, prim, R1,R2) 01294 lib%AB=B-A 01295 lib%CD=D-C 01296 CALL get_eris(n_c, n_d, n_a, n_b, lib, prim, p_work, a_mysize) 01297 DO i=1,mysize 01298 max_val = MAX(max_val, ABS(p_work(i))) 01299 END DO 01300 CASE(5) 01301 CALL build_quartet_data_screen(C,D,A,B,Zeta_C, Zeta_D, Zeta_A, Zeta_B, m_max,& 01302 potential_parameter, prim, R1,R2) 01303 lib%AB=C-D 01304 lib%CD=A-B 01305 CALL get_eris(n_b, n_a, n_d, n_c, lib, prim, p_work, a_mysize) 01306 DO i=1,mysize 01307 max_val = MAX(max_val, ABS(p_work(i))) 01308 END DO 01309 CASE(6) 01310 CALL build_quartet_data_screen(C,D,B,A,Zeta_C, Zeta_D, Zeta_B, Zeta_A, m_max,& 01311 potential_parameter, prim, R1,R2) 01312 lib%AB=C-D 01313 lib%CD=B-A 01314 CALL get_eris(n_a, n_b, n_d, n_c, lib, prim, p_work, a_mysize) 01315 DO i=1,mysize 01316 max_val = MAX(max_val, ABS(p_work(i))) 01317 END DO 01318 CASE(7) 01319 CALL build_quartet_data_screen(D,C,A,B,Zeta_D, Zeta_C, Zeta_A, Zeta_B, m_max,& 01320 potential_parameter, prim, R1,R2) 01321 lib%AB=D-C 01322 lib%CD=A-B 01323 CALL get_eris(n_b, n_a, n_c, n_d, lib, prim, p_work, a_mysize) 01324 DO i=1,mysize 01325 max_val = MAX(max_val, ABS(p_work(i))) 01326 END DO 01327 CASE(8) 01328 CALL build_quartet_data_screen(D,C,B,A,Zeta_D, Zeta_C, Zeta_B, Zeta_A, m_max,& 01329 potential_parameter, prim, R1,R2) 01330 lib%AB=D-C 01331 lib%CD=B-A 01332 CALL get_eris(n_a, n_b, n_c, n_d, lib, prim, p_work, a_mysize) 01333 DO i=1,mysize 01334 max_val = MAX(max_val, ABS(p_work(i))) 01335 END DO 01336 END SELECT 01337 ELSE 01338 CALL build_quartet_data_screen(A,B,C,D,Zeta_A, Zeta_B, Zeta_C, Zeta_D, m_max,& 01339 potential_parameter, prim, R1,R2) 01340 max_val = ABS(prim%F(1)) 01341 END IF 01342 01343 END SUBROUTINE evaluate_eri_screen 01344 01345 01346 ! ***************************************************************************** 01353 01354 SUBROUTINE evaluate_eri(lib, nproducts, pgf_product_list, & 01355 n_a,n_b,n_c,n_d,& 01356 ncoa,ncob,ncoc,ncod,& 01357 nsgfa, nsgfb, nsgfc, nsgfd,& 01358 primitives, max_contraction, tmp_max, & 01359 eps_schwarz, neris, & 01360 ZetaInv,EtaInv, & 01361 s_offset_a, s_offset_b, s_offset_c, s_offset_d,& 01362 nl_a, nl_b, nl_c, nl_d, nsoa, nsob, nsoc, nsod,& 01363 sphi_a, sphi_b, sphi_c, sphi_d, work,work2,buffer1,buffer2, & 01364 primitives_tmp, p_work) 01365 01366 TYPE(lib_int) :: lib 01367 INTEGER, INTENT(IN) :: nproducts 01368 TYPE(hfx_pgf_product_list), 01369 DIMENSION(nproducts) :: pgf_product_list 01370 INTEGER, INTENT(IN) :: n_a, n_b, n_c, n_d, ncoa, 01371 ncob, ncoc, ncod, nsgfa, 01372 nsgfb, nsgfc, nsgfd 01373 REAL(dp), 01374 DIMENSION(nsgfa, nsgfb, nsgfc, nsgfd) :: primitives 01375 REAL(dp) :: max_contraction, tmp_max, 01376 eps_schwarz 01377 INTEGER(int_8) :: neris 01378 REAL(dp), INTENT(IN) :: ZetaInv, EtaInv 01379 INTEGER :: s_offset_a, s_offset_b, 01380 s_offset_c, s_offset_d, nl_a, 01381 nl_b, nl_c, nl_d, nsoa, nsob, 01382 nsoc, nsod 01383 REAL(dp), DIMENSION(ncoa, nsoa*nl_a), 01384 INTENT(IN) :: sphi_a 01385 REAL(dp), DIMENSION(ncob, nsob*nl_b), 01386 INTENT(IN) :: sphi_b 01387 REAL(dp), DIMENSION(ncoc, nsoc*nl_c), 01388 INTENT(IN) :: sphi_c 01389 REAL(dp), DIMENSION(ncod, nsod*nl_d), 01390 INTENT(IN) :: sphi_d 01391 REAL(dp) :: work(ncoa*ncob*ncoc*ncod), 01392 work2(ncoa,ncob,ncoc,ncod) 01393 REAL(dp), DIMENSION(ncoa*ncob*ncoc*ncod) :: buffer1, buffer2 01394 REAL(dp), DIMENSION(nsoa*nl_a, nsob*nl_b& , nsoc*nl_c, nsod*nl_d) :: primitives_tmp 01395 REAL(dp), DIMENSION(:), POINTER :: p_work 01396 01397 INTEGER :: a_mysize(1), i, j, k, l, 01398 m_max, mysize, p1, p2, p3, 01399 p4, perm_case 01400 REAL(dp) :: A(3), AB(3), B(3), C(3), 01401 CD(3), D(3), P(3), Q(3), Rho, 01402 RhoInv, W(3), ZetapEtaInv 01403 TYPE(prim_data), TARGET :: prim 01404 01405 m_max = n_a+n_b+n_c+n_d 01406 mysize = ncoa*ncob*ncoc*ncod 01407 a_mysize = mysize 01408 01409 work = 0.0_dp 01410 IF(m_max/=0) THEN 01411 perm_case = 1 01412 IF(n_a<n_b) THEN 01413 perm_case = perm_case + 1 01414 END IF 01415 IF(n_c<n_d) THEN 01416 perm_case = perm_case + 2 01417 END IF 01418 IF(n_a+n_b > n_c+n_d) THEN 01419 perm_case = perm_case + 4 01420 END IF 01421 SELECT CASE(perm_case) 01422 CASE(1) 01423 DO i = 1,nproducts 01424 A = pgf_product_list(i)%ra 01425 B = pgf_product_list(i)%rb 01426 C = pgf_product_list(i)%rc 01427 D = pgf_product_list(i)%rd 01428 Rho = pgf_product_list(i)%Rho 01429 RhoInv = pgf_product_list(i)%RhoInv 01430 P = pgf_product_list(i)%P 01431 Q = pgf_product_list(i)%Q 01432 W = pgf_product_list(i)%W 01433 AB = pgf_product_list(i)%AB 01434 CD = pgf_product_list(i)%CD 01435 ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv 01436 prim%F(1:m_max+1) = pgf_product_list(i)%Fm(1:m_max+1) 01437 01438 CALL build_quartet_data(prim, A, C,& 01439 ZetaInv, EtaInv, ZetapEtaInv, Rho, RhoInv,& 01440 P, Q, W) 01441 lib%AB=AB!A-B 01442 lib%CD=CD!C-D 01443 CALL get_eris(n_d, n_c, n_b, n_a, lib, prim, p_work, a_mysize) 01444 work(1:mysize) = work(1:mysize) + p_work(1:mysize) 01445 neris = neris + mysize 01446 END DO 01447 DO i=1,mysize 01448 tmp_max = MAX(tmp_max, ABS(work(i))) 01449 END DO 01450 tmp_max = tmp_max*max_contraction 01451 IF(tmp_max<eps_schwarz) THEN 01452 RETURN 01453 END IF 01454 01455 DO i = 1,ncoa 01456 p1 = (i-1) *ncob 01457 DO j = 1,ncob 01458 p2 = (p1 + j-1)*ncoc 01459 DO k = 1,ncoc 01460 p3 = (p2 + k-1)*ncod 01461 DO l = 1,ncod 01462 p4 = p3 + l 01463 work2(i,j,k,l) = work(p4) 01464 END DO 01465 END DO 01466 END DO 01467 END DO 01468 CASE(2) 01469 DO i = 1,nproducts 01470 A = pgf_product_list(i)%ra 01471 B = pgf_product_list(i)%rb 01472 C = pgf_product_list(i)%rc 01473 D = pgf_product_list(i)%rd 01474 Rho = pgf_product_list(i)%Rho 01475 RhoInv = pgf_product_list(i)%RhoInv 01476 P = pgf_product_list(i)%P 01477 Q = pgf_product_list(i)%Q 01478 W = pgf_product_list(i)%W 01479 AB = pgf_product_list(i)%AB 01480 CD = pgf_product_list(i)%CD 01481 ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv 01482 prim%F(1:m_max+1) = pgf_product_list(i)%Fm(1:m_max+1) 01483 01484 CALL build_quartet_data(prim, B, C,& 01485 ZetaInv, EtaInv, ZetapEtaInv, Rho, RhoInv,& 01486 P,Q,W) 01487 lib%AB=-AB!B-A 01488 lib%CD=CD!C-D 01489 CALL get_eris(n_d, n_c, n_a, n_b, lib, prim, p_work, a_mysize) 01490 work(1:mysize) = work(1:mysize) + p_work(1:mysize) 01491 neris = neris + mysize 01492 END DO 01493 DO i=1,mysize 01494 tmp_max = MAX(tmp_max, ABS(work(i))) 01495 END DO 01496 tmp_max = tmp_max*max_contraction 01497 IF(tmp_max<eps_schwarz) THEN 01498 RETURN 01499 END IF 01500 01501 DO j = 1,ncob 01502 p1 = (j-1)*ncoa 01503 DO i = 1,ncoa 01504 p2 = (p1 + i-1)*ncoc 01505 DO k = 1,ncoc 01506 p3 = (p2 + k-1)*ncod 01507 DO l = 1,ncod 01508 p4 = p3 + l 01509 work2(i,j,k,l) = work(p4) 01510 END DO 01511 END DO 01512 END DO 01513 END DO 01514 CASE(3) 01515 DO i = 1,nproducts 01516 A = pgf_product_list(i)%ra 01517 B = pgf_product_list(i)%rb 01518 C = pgf_product_list(i)%rc 01519 D = pgf_product_list(i)%rd 01520 Rho = pgf_product_list(i)%Rho 01521 RhoInv = pgf_product_list(i)%RhoInv 01522 P = pgf_product_list(i)%P 01523 Q = pgf_product_list(i)%Q 01524 W = pgf_product_list(i)%W 01525 AB = pgf_product_list(i)%AB 01526 CD = pgf_product_list(i)%CD 01527 ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv 01528 prim%F(1:m_max+1) = pgf_product_list(i)%Fm(1:m_max+1) 01529 01530 CALL build_quartet_data(prim, A, D,& 01531 ZetaInv, EtaInv, ZetapEtaInv, Rho, RhoInv,& 01532 P,Q,W) 01533 lib%AB=AB!A-B 01534 lib%CD=-CD!D-C 01535 CALL get_eris(n_c, n_d, n_b, n_a, lib, prim, p_work, a_mysize) 01536 work(1:mysize) = work(1:mysize) + p_work(1:mysize) 01537 neris = neris + mysize 01538 END DO 01539 DO i=1,mysize 01540 tmp_max = MAX(tmp_max, ABS(work(i))) 01541 END DO 01542 tmp_max = tmp_max*max_contraction 01543 IF(tmp_max<eps_schwarz) THEN 01544 RETURN 01545 END IF 01546 01547 DO i = 1,ncoa 01548 p1 = (i-1)*ncob 01549 DO j = 1,ncob 01550 p2 = (p1 + j-1)*ncod 01551 DO l = 1,ncod 01552 p3 = (p2 + l-1) * ncoc 01553 DO k = 1,ncoc 01554 p4 = p3+k 01555 work2(i,j,k,l) = work(p4) 01556 END DO 01557 END DO 01558 END DO 01559 END DO 01560 CASE(4) 01561 DO i = 1,nproducts 01562 A = pgf_product_list(i)%ra 01563 B = pgf_product_list(i)%rb 01564 C = pgf_product_list(i)%rc 01565 D = pgf_product_list(i)%rd 01566 Rho = pgf_product_list(i)%Rho 01567 RhoInv = pgf_product_list(i)%RhoInv 01568 P = pgf_product_list(i)%P 01569 Q = pgf_product_list(i)%Q 01570 W = pgf_product_list(i)%W 01571 AB = pgf_product_list(i)%AB 01572 CD = pgf_product_list(i)%CD 01573 ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv 01574 prim%F(1:m_max+1) = pgf_product_list(i)%Fm(1:m_max+1) 01575 01576 CALL build_quartet_data(prim, B, D,& 01577 ZetaInv, EtaInv, ZetapEtaInv, Rho, RhoInv,& 01578 P, Q, W) 01579 lib%AB=-AB!B-A 01580 lib%CD=-CD!D-C 01581 CALL get_eris(n_c, n_d, n_a, n_b, lib, prim, p_work, a_mysize) 01582 work(1:mysize) = work(1:mysize) + p_work(1:mysize) 01583 neris = neris + mysize 01584 END DO 01585 DO i=1,mysize 01586 tmp_max = MAX(tmp_max, ABS(work(i))) 01587 END DO 01588 tmp_max = tmp_max*max_contraction 01589 IF(tmp_max<eps_schwarz) THEN 01590 RETURN 01591 END IF 01592 01593 DO j = 1,ncob 01594 p1 = (j-1)*ncoa 01595 DO i = 1,ncoa 01596 p2 = (p1 + i-1)*ncod 01597 DO l = 1,ncod 01598 p3 = (p2 + l-1)*ncoc 01599 DO k = 1,ncoc 01600 p4 = p3 + k 01601 work2(i,j,k,l) = work(p4) 01602 END DO 01603 END DO 01604 END DO 01605 END DO 01606 CASE(5) 01607 DO i = 1,nproducts 01608 A = pgf_product_list(i)%ra 01609 B = pgf_product_list(i)%rb 01610 C = pgf_product_list(i)%rc 01611 D = pgf_product_list(i)%rd 01612 Rho = pgf_product_list(i)%Rho 01613 RhoInv = pgf_product_list(i)%RhoInv 01614 P = pgf_product_list(i)%P 01615 Q = pgf_product_list(i)%Q 01616 W = pgf_product_list(i)%W 01617 AB = pgf_product_list(i)%AB 01618 CD = pgf_product_list(i)%CD 01619 ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv 01620 prim%F(1:m_max+1) = pgf_product_list(i)%Fm(1:m_max+1) 01621 01622 CALL build_quartet_data(prim, C, A,& 01623 EtaInv, ZetaInv, ZetapEtaInv, Rho, RhoInv,& 01624 Q, P, W) 01625 lib%AB=CD!C-D 01626 lib%CD=AB!A-B 01627 CALL get_eris(n_b, n_a, n_d, n_c, lib, prim, p_work, a_mysize) 01628 work(1:mysize) = work(1:mysize) + p_work(1:mysize) 01629 neris = neris + mysize 01630 END DO 01631 DO i=1,mysize 01632 tmp_max = MAX(tmp_max, ABS(work(i))) 01633 END DO 01634 tmp_max = tmp_max*max_contraction 01635 IF(tmp_max<eps_schwarz) THEN 01636 RETURN 01637 END IF 01638 01639 DO k = 1,ncoc 01640 p1 = (k-1)*ncod 01641 DO l = 1,ncod 01642 p2 = (p1 + l-1)*ncoa 01643 DO i = 1,ncoa 01644 p3 = (p2 + i-1)*ncob 01645 DO j = 1,ncob 01646 p4 = p3+j 01647 work2(i,j,k,l) = work(p4) 01648 END DO 01649 END DO 01650 END DO 01651 END DO 01652 CASE(6) 01653 DO i = 1,nproducts 01654 A = pgf_product_list(i)%ra 01655 B = pgf_product_list(i)%rb 01656 C = pgf_product_list(i)%rc 01657 D = pgf_product_list(i)%rd 01658 Rho = pgf_product_list(i)%Rho 01659 RhoInv = pgf_product_list(i)%RhoInv 01660 P = pgf_product_list(i)%P 01661 Q = pgf_product_list(i)%Q 01662 W = pgf_product_list(i)%W 01663 AB = pgf_product_list(i)%AB 01664 CD = pgf_product_list(i)%CD 01665 ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv 01666 prim%F(1:m_max+1) = pgf_product_list(i)%Fm(1:m_max+1) 01667 01668 CALL build_quartet_data(prim, C, B, & 01669 EtaInv, ZetaInv, ZetapEtaInv, Rho, RhoInv,& 01670 Q, P, W) 01671 lib%AB=CD!C-D 01672 lib%CD=-AB!B-A 01673 CALL get_eris(n_a, n_b, n_d, n_c, lib, prim, p_work, a_mysize) 01674 work(1:mysize) = work(1:mysize) + p_work(1:mysize) 01675 neris = neris + mysize 01676 END DO 01677 DO i=1,mysize 01678 tmp_max = MAX(tmp_max, ABS(work(i))) 01679 END DO 01680 tmp_max = tmp_max*max_contraction 01681 IF(tmp_max<eps_schwarz) THEN 01682 RETURN 01683 END IF 01684 01685 DO k = 1,ncoc 01686 p1 = (k-1)*ncod 01687 DO l = 1,ncod 01688 p2 = (p1 + l-1)*ncob 01689 DO j = 1,ncob 01690 p3 = (p2 + j-1)*ncoa 01691 DO i = 1,ncoa 01692 p4 = p3+i 01693 work2(i,j,k,l) = work(p4) 01694 END DO 01695 END DO 01696 END DO 01697 END DO 01698 CASE(7) 01699 DO i = 1,nproducts 01700 A = pgf_product_list(i)%ra 01701 B = pgf_product_list(i)%rb 01702 C = pgf_product_list(i)%rc 01703 D = pgf_product_list(i)%rd 01704 Rho = pgf_product_list(i)%Rho 01705 RhoInv = pgf_product_list(i)%RhoInv 01706 P = pgf_product_list(i)%P 01707 Q = pgf_product_list(i)%Q 01708 W = pgf_product_list(i)%W 01709 AB = pgf_product_list(i)%AB 01710 CD = pgf_product_list(i)%CD 01711 ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv 01712 prim%F(1:m_max+1) = pgf_product_list(i)%Fm(1:m_max+1) 01713 01714 CALL build_quartet_data(prim, D, A,& 01715 EtaInv, ZetaInv, ZetapEtaInv, Rho, RhoInv,& 01716 Q, P, W) 01717 lib%AB=-CD!D-C 01718 lib%CD=AB!A-B 01719 CALL get_eris(n_b, n_a, n_c, n_d, lib, prim, p_work, a_mysize) 01720 work(1:mysize) = work(1:mysize) + p_work(1:mysize) 01721 neris = neris + mysize 01722 END DO 01723 DO i=1,mysize 01724 tmp_max = MAX(tmp_max, ABS(work(i))) 01725 END DO 01726 tmp_max = tmp_max*max_contraction 01727 IF(tmp_max<eps_schwarz) THEN 01728 RETURN 01729 END IF 01730 01731 DO l = 1,ncod 01732 p1 = (l-1)*ncoc 01733 DO k = 1,ncoc 01734 p2 = (p1 + k-1) * ncoa 01735 DO i = 1,ncoa 01736 p3 = (p2 + i-1) *ncob 01737 DO j = 1,ncob 01738 p4 = p3+j 01739 work2(i,j,k,l) = work(p4) 01740 END DO 01741 END DO 01742 END DO 01743 END DO 01744 CASE(8) 01745 DO i = 1,nproducts 01746 A = pgf_product_list(i)%ra 01747 B = pgf_product_list(i)%rb 01748 C = pgf_product_list(i)%rc 01749 D = pgf_product_list(i)%rd 01750 Rho = pgf_product_list(i)%Rho 01751 RhoInv = pgf_product_list(i)%RhoInv 01752 P = pgf_product_list(i)%P 01753 Q = pgf_product_list(i)%Q 01754 W = pgf_product_list(i)%W 01755 AB = pgf_product_list(i)%AB 01756 CD = pgf_product_list(i)%CD 01757 ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv 01758 prim%F(1:m_max+1) = pgf_product_list(i)%Fm(1:m_max+1) 01759 01760 CALL build_quartet_data(prim, D, B,& 01761 EtaInv, ZetaInv, ZetapEtaInv, Rho, RhoInv,& 01762 Q, P, W) 01763 lib%AB=-CD!D-C 01764 lib%CD=-AB!B-A 01765 CALL get_eris(n_a, n_b, n_c, n_d, lib, prim, p_work, a_mysize) 01766 work(1:mysize) = work(1:mysize) + p_work(1:mysize) 01767 neris = neris + mysize 01768 END DO 01769 DO i=1,mysize 01770 tmp_max = MAX(tmp_max, ABS(work(i))) 01771 END DO 01772 tmp_max = tmp_max*max_contraction 01773 IF(tmp_max<eps_schwarz) THEN 01774 RETURN 01775 END IF 01776 01777 DO l = 1,ncod 01778 p1 = (l-1)*ncoc 01779 DO k = 1,ncoc 01780 p2 = (p1 + k-1) * ncob 01781 DO j = 1,ncob 01782 p3 = (p2 + j-1) * ncoa 01783 DO i = 1,ncoa 01784 p4 = p3 + i 01785 work2(i,j,k,l) = work(p4) 01786 END DO 01787 END DO 01788 END DO 01789 END DO 01790 END SELECT 01791 ELSE 01792 DO i = 1,nproducts 01793 A = pgf_product_list(i)%ra 01794 B = pgf_product_list(i)%rb 01795 C = pgf_product_list(i)%rc 01796 D = pgf_product_list(i)%rd 01797 Rho = pgf_product_list(i)%Rho 01798 RhoInv = pgf_product_list(i)%RhoInv 01799 P = pgf_product_list(i)%P 01800 Q = pgf_product_list(i)%Q 01801 W = pgf_product_list(i)%W 01802 ZetapEtaInv = pgf_product_list(i)%ZetapEtaInv 01803 prim%F(1:m_max+1) = pgf_product_list(i)%Fm(1:m_max+1) 01804 01805 CALL build_quartet_data(prim, A, C,& 01806 ZetaInv, EtaInv, ZetapEtaInv, Rho, RhoInv,& 01807 P, Q, W) 01808 work(1) = work(1) + prim%F(1) 01809 neris = neris + mysize 01810 END DO 01811 work2(1,1,1,1) = work(1) 01812 tmp_max = max_contraction*ABS(work(1)) 01813 IF( tmp_max < eps_schwarz ) RETURN 01814 END IF 01815 01816 IF( tmp_max < eps_schwarz ) RETURN 01817 primitives_tmp = 0.0_dp 01818 01819 CALL contract(ncoa, ncob, ncoc, ncod, nsoa, nsob, nsoc, nsod,& 01820 n_a, n_b, n_c, n_d, nl_a, nl_b, nl_c, nl_d, work2,& 01821 sphi_a, & 01822 sphi_b, & 01823 sphi_c, & 01824 sphi_d, & 01825 primitives_tmp, & 01826 buffer1, buffer2) 01827 01828 primitives(s_offset_a+1:s_offset_a+nsoa*nl_a,& 01829 s_offset_b+1:s_offset_b+nsob*nl_b,& 01830 s_offset_c+1:s_offset_c+nsoc*nl_c,& 01831 s_offset_d+1:s_offset_d+nsod*nl_d) = & 01832 primitives(s_offset_a+1:s_offset_a+nsoa*nl_a,& 01833 s_offset_b+1:s_offset_b+nsob*nl_b,& 01834 s_offset_c+1:s_offset_c+nsoc*nl_c,& 01835 s_offset_d+1:s_offset_d+nsod*nl_d) + primitives_tmp(:,:,:,:) 01836 01837 END SUBROUTINE evaluate_eri 01838 01839 ! ***************************************************************************** 01845 01846 SUBROUTINE build_quartet_data(prim, A, C,& 01847 ZetaInv,EtaInv,ZetapEtaInv,Rho,RhoInv,& 01848 P,Q,W) 01849 01850 TYPE(prim_data) :: prim 01851 REAL(KIND=dp), INTENT(IN) :: A(3), C(3) 01852 REAL(dp), INTENT(IN) :: ZetaInv, EtaInv, ZetapEtaInv, 01853 Rho, RhoInv, P(3), Q(3), W(3) 01854 01855 prim%U(1:3,1) = P-A 01856 prim%U(1:3,3) = Q-C 01857 prim%U(1:3,5) = W-P 01858 prim%U(1:3,6) = W-Q 01859 prim%oo2z = 0.5_dp*ZetaInv 01860 prim%oo2n = 0.5_dp*EtaInv 01861 prim%oo2zn = 0.5_dp*ZetapEtaInv 01862 prim%poz = Rho*ZetaInv 01863 prim%pon = Rho*EtaInv 01864 prim%oo2p = 0.5_dp*RhoInv 01865 END SUBROUTINE build_quartet_data 01866 01867 END MODULE hfx_libint_interface 01868
1.7.3