|
CP2K 2.4 (Revision 12889)
|
00001 !-----------------------------------------------------------------------------! 00002 ! CP2K: A general program to perform molecular dynamics simulations ! 00003 ! Copyright (C) 2000 - 2013 CP2K developers group ! 00004 !-----------------------------------------------------------------------------! 00005 00006 ! ***************************************************************************** 00012 MODULE qs_util 00013 00014 USE f77_blas 00015 USE kinds, ONLY: dp,& 00016 dp_size 00017 USE mathconstants, ONLY: & 00018 dfac, fac, fourpi, oorootpi, rootpi, sqrt105, sqrt15, sqrt2, sqrt21, & 00019 sqrt3, sqrt35, sqrt5, sqrt7, sqrthalf 00020 USE orbital_pointers, ONLY: indco,& 00021 nco,& 00022 ncoset,& 00023 nso,& 00024 nsoset 00025 USE orbital_transformation_matrices, ONLY: orbtramat 00026 USE termination, ONLY: stop_memory,& 00027 stop_program 00028 #include "cp_common_uses.h" 00029 00030 IMPLICIT NONE 00031 00032 PRIVATE 00033 00034 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_util' 00035 00036 !MK sqrt* constants moved to mathconstants 00037 00038 REAL(KIND=dp), PARAMETER :: s_root1o4pi=0.5_dp*oorootpi 00039 REAL(KIND=dp), PARAMETER :: root4pi=2.0_dp*rootpi 00040 REAL(KIND=dp), PARAMETER :: s_root3o4pi=sqrt3*s_root1o4pi 00041 REAL(KIND=dp), PARAMETER :: root4pio3=root4pi/sqrt3 00042 REAL(KIND=dp), PARAMETER :: root4pio5=root4pi/sqrt5 00043 REAL(KIND=dp), PARAMETER :: s_root15o4pi=sqrt15*s_root1o4pi 00044 REAL(KIND=dp), PARAMETER :: root4pio15=root4pi/sqrt15 00045 REAL(KIND=dp), PARAMETER :: s_root105o4pi=sqrt105*s_root1o4pi 00046 REAL(KIND=dp), PARAMETER :: root4pio105=root4pi/sqrt105 00047 REAL(KIND=dp), PARAMETER :: s_root1o16pi=0.25_dp*oorootpi 00048 REAL(KIND=dp), PARAMETER :: root16pi=4.0_dp*rootpi 00049 REAL(KIND=dp), PARAMETER :: s_root5o16pi=sqrt5*s_root1o16pi 00050 REAL(KIND=dp), PARAMETER :: root16pio5=root16pi/sqrt5 00051 REAL(KIND=dp), PARAMETER :: s_2root5o16pi=2.0_dp*s_root5o16pi 00052 REAL(KIND=dp), PARAMETER :: root16pio5o2=root16pio5*0.5_dp 00053 REAL(KIND=dp), PARAMETER :: s_3root5o16pi=3.0_dp*s_root5o16pi 00054 REAL(KIND=dp), PARAMETER :: root16pio5o3=root16pio5/3.0_dp 00055 REAL(KIND=dp), PARAMETER :: s_18root5o16pi=18.0_dp*s_root5o16pi 00056 REAL(KIND=dp), PARAMETER :: root16pio5o18=root16pio5/18.0_dp 00057 REAL(KIND=dp), PARAMETER :: s_2root7o16pi=2.0_dp*sqrt7*s_root1o16pi 00058 REAL(KIND=dp), PARAMETER :: root16pio7o2=root16pi/sqrt7*0.5_dp 00059 REAL(KIND=dp), PARAMETER :: s_3root7o16pi=3.0_dp*sqrt7*s_root1o16pi 00060 REAL(KIND=dp), PARAMETER :: root16pio7o3=root16pi/sqrt7/3.0_dp 00061 REAL(KIND=dp), PARAMETER :: s_root15o16pi=sqrt15*s_root1o16pi 00062 REAL(KIND=dp), PARAMETER :: root16pio15=root16pi/sqrt15 00063 REAL(KIND=dp), PARAMETER :: s_3root35o16pi=sqrt5*s_3root7o16pi 00064 REAL(KIND=dp), PARAMETER :: root16pio35o3=root16pio7o3/sqrt5 00065 REAL(KIND=dp), PARAMETER :: s_root105o16pi=0.5_dp*s_root105o4pi 00066 REAL(KIND=dp), PARAMETER :: root16pio105=root4pio105*2.0_dp 00067 REAL(KIND=dp), PARAMETER :: s_root1o32pi=0.25_dp*sqrthalf*oorootpi 00068 REAL(KIND=dp), PARAMETER :: root32pi=root16pi*sqrt2 00069 REAL(KIND=dp), PARAMETER :: s_3root5o32pi=3.0_dp*sqrt5*s_root1o32pi 00070 REAL(KIND=dp), PARAMETER :: root32pio5o3=root32pi/sqrt5/3.0_dp 00071 REAL(KIND=dp), PARAMETER :: s_9root5o32pi=9.0_dp*sqrt5*s_root1o32pi 00072 REAL(KIND=dp), PARAMETER :: root32pio5o9=root32pi/sqrt5/9.0_dp 00073 REAL(KIND=dp), PARAMETER :: s_12root5o32pi=12.0_dp*sqrt5*s_root1o32pi 00074 REAL(KIND=dp), PARAMETER :: root32pio5o12=root32pi/sqrt5/12.0_dp 00075 REAL(KIND=dp), PARAMETER :: s_root21o32pi=sqrt21*s_root1o32pi 00076 REAL(KIND=dp), PARAMETER :: root32pio21=root32pi/sqrt21 00077 REAL(KIND=dp), PARAMETER :: s_4root21o32pi=4.0_dp*s_root21o32pi 00078 REAL(KIND=dp), PARAMETER :: root32pio21o4=root32pio21/4.0_dp 00079 REAL(KIND=dp), PARAMETER :: s_root35o32pi=sqrt35*s_root1o32pi 00080 REAL(KIND=dp), PARAMETER :: root32pio35=root32pi/sqrt35 00081 REAL(KIND=dp), PARAMETER :: s_3root35o32pi=3.0_dp*s_root35o32pi 00082 REAL(KIND=dp), PARAMETER :: s_9root35o32pi=9.0_dp*s_root35o32pi 00083 REAL(KIND=dp), PARAMETER :: s_18root35o32pi=18.0_dp*s_root35o32pi 00084 REAL(KIND=dp), PARAMETER :: s_root1o64pi=0.125_dp*oorootpi 00085 REAL(KIND=dp), PARAMETER :: s_3root5o64pi=3.0_dp*sqrt5*s_root1o64pi 00086 REAL(KIND=dp), PARAMETER :: s_18root5o64pi=18.0_dp*sqrt5*s_root1o64pi 00087 REAL(KIND=dp), PARAMETER :: s_root1o256pi=0.0625_dp*oorootpi 00088 REAL(KIND=dp), PARAMETER :: s_3root1o256pi=3.0_dp*s_root1o256pi 00089 REAL(KIND=dp), PARAMETER :: s_9root1o256pi=9.0_dp*s_root1o256pi 00090 REAL(KIND=dp), PARAMETER :: s_18root1o256pi=18.0_dp*s_root1o256pi 00091 REAL(KIND=dp), PARAMETER :: s_24root1o256pi=24.0_dp*s_root1o256pi 00092 REAL(KIND=dp), PARAMETER :: s_72root1o256pi=72.0_dp*s_root1o256pi 00093 REAL(KIND=dp), PARAMETER :: s_3root35o256pi=3.0_dp*sqrt35*s_root1o256pi 00094 REAL(KIND=dp), PARAMETER :: s_18root35o256pi=18.0_dp*sqrt35*s_root1o256pi 00095 00096 ! *** Public subroutines *** 00097 00098 PUBLIC :: cart2sph_mat,& 00099 exp_radius,& 00100 gauss_exponent,& 00101 gaussint_sph,& 00102 trace_r_AxB, & 00103 trace_r_AxB_new, & 00104 transform_c2s,& 00105 transform_c2s_new,& 00106 transform_s2c 00107 00108 CONTAINS 00109 00110 ! ***************************************************************************** 00124 FUNCTION gauss_exponent(l,radius,threshold,prefactor) RESULT(exponent) 00125 INTEGER, INTENT(IN) :: l 00126 REAL(KIND=dp), INTENT(IN) :: radius, threshold, prefactor 00127 REAL(KIND=dp) :: exponent 00128 00129 exponent = 0.0_dp 00130 00131 IF (radius < 1.0E-6_dp) RETURN 00132 IF (threshold < 1.0E-12_dp) RETURN 00133 00134 exponent = LOG(ABS(prefactor)*radius**l/threshold)/radius**2 00135 00136 END FUNCTION gauss_exponent 00137 00138 ! ***************************************************************************** 00160 FUNCTION exp_radius(l,alpha,threshold,prefactor,epsin) RESULT(radius) 00161 INTEGER, INTENT(IN) :: l 00162 REAL(KIND=dp), INTENT(IN) :: alpha, threshold, prefactor 00163 REAL(KIND=dp), INTENT(IN), OPTIONAL :: epsin 00164 REAL(KIND=dp) :: radius 00165 00166 CHARACTER(len=*), PARAMETER :: routineN = 'exp_radius', 00167 routineP = moduleN//':'//routineN 00168 INTEGER, PARAMETER :: maxiter = 5000 00169 00170 INTEGER :: iter 00171 REAL(KIND=dp) :: a, ar2, d, epsiter, g, r, 00172 rhigh, rlow, rmid, t 00173 00174 IF (PRESENT(epsin)) THEN 00175 epsiter=epsin 00176 ELSE 00177 epsiter=EPSILON(epsiter)*(1.0E-12_dp / 2.22044604925031E-16_dp) 00178 ENDIF 00179 00180 ! Initialize function value 00181 00182 radius = 0.0_dp 00183 00184 ! Load and check parameter values 00185 00186 IF (l < 0) THEN 00187 CALL stop_program(routineN,moduleN,__LINE__,& 00188 "The angular momentum quantum number is negative") 00189 END IF 00190 00191 IF (alpha == 0.0_dp) THEN 00192 CALL stop_program(routineN,moduleN,__LINE__,& 00193 "The Gaussian function exponent is zero") 00194 ELSE 00195 a = ABS(alpha) 00196 END IF 00197 00198 IF (threshold == 0.0_dp) THEN 00199 CALL stop_program(routineN,moduleN,__LINE__,& 00200 "The requested threshold is zero") 00201 ELSE 00202 t = ABS(threshold) 00203 END IF 00204 00205 IF (prefactor == 0.0_dp) THEN 00206 RETURN 00207 ELSE 00208 d = ABS(prefactor) 00209 END IF 00210 00211 ! Calculate the Maximum g(r) 00212 r = SQRT(0.5_dp*REAL(l,dp)/a) 00213 ar2 = a*r*r 00214 00215 IF (l == 0) THEN 00216 g = d 00217 ELSE 00218 g = d*r**l*EXP(-ar2) 00219 END IF 00220 00221 IF (t > g) THEN 00222 RETURN 00223 END IF 00224 00225 rlow = r 00226 rhigh = 2.0_dp*rlow+1.0_dp 00227 iter=0 00228 DO 00229 iter=iter+1 00230 IF (iter.gt.maxiter) THEN 00231 CALL stop_program(routineN,moduleN,__LINE__,& 00232 "Maximum number of iterations exceeded") 00233 END IF 00234 g = d*rhigh**l*EXP(-a*rhigh**2) 00235 IF (g < t) EXIT 00236 rlow = rhigh 00237 rhigh = 2.0_dp*rlow+1.0_dp 00238 ENDDO 00239 00240 DO iter=1,maxiter 00241 rmid= (rlow+rhigh)*0.5_dp 00242 ar2 = a*rmid*rmid 00243 g = d*rmid**l*EXP(-ar2) 00244 IF (g.lt.t) THEN 00245 rhigh=rmid 00246 ELSE 00247 rlow=rmid 00248 ENDIF 00249 IF (ABS(rhigh-rlow).lt.epsiter) THEN 00250 radius=rhigh 00251 RETURN 00252 ENDIF 00253 ENDDO 00254 CALL stop_program(routineN,moduleN,__LINE__,& 00255 "Maximum number of iterations exceeded") 00256 00257 END FUNCTION exp_radius 00258 00259 ! ***************************************************************************** 00260 FUNCTION gaussint_sph(alpha,l) 00261 00262 ! calculates the radial integral over a spherical Gaussian 00263 ! of the form 00264 ! r**(2+l) * exp(-alpha * r**2) 00265 ! 00266 00267 REAL(dp), INTENT(IN) :: alpha 00268 INTEGER, INTENT(IN) :: l 00269 REAL(dp) :: gaussint_sph 00270 00271 IF ((l/2)*2==l) THEN 00272 !even l: 00273 gaussint_sph=ROOTPI * 0.5_dp**(l/2+2) * dfac(l+1 )& 00274 /SQRT(alpha)**(l+3) 00275 ELSE 00276 !odd l: 00277 gaussint_sph=0.5_dp * fac((l+1 )/2) /SQRT(alpha)**(l+3) 00278 ENDIF 00279 00280 END FUNCTION gaussint_sph 00281 00282 ! ***************************************************************************** 00283 FUNCTION trace_r_AxB(A,lda,B,ldb,m,n) 00284 00285 INTEGER, INTENT(in) :: lda 00286 REAL(dp), INTENT(in) :: A(lda,*) 00287 INTEGER, INTENT(in) :: ldb 00288 REAL(dp), INTENT(in) :: B(ldb,*) 00289 INTEGER, INTENT(in) :: m, n 00290 REAL(dp) :: trace_r_AxB 00291 00292 INTEGER :: i1, i2, imod, mminus3 00293 REAL(dp) :: t1, t2, t3, t4 00294 00295 t1=0._dp 00296 t2=0._dp 00297 t3=0._dp 00298 t4=0._dp 00299 imod=MODULO(m,4) 00300 SELECT CASE (imod) 00301 CASE (0) 00302 DO i2=1,n 00303 DO i1=1,m,4 00304 t1=t1+A(i1,i2)*B(i1,i2) 00305 t2=t2+A(i1+1,i2)*B(i1+1,i2) 00306 t3=t3+A(i1+2,i2)*B(i1+2,i2) 00307 t4=t4+A(i1+3,i2)*B(i1+3,i2) 00308 ENDDO 00309 ENDDO 00310 CASE (1) 00311 mminus3=m-3 00312 DO i2=1,n 00313 DO i1=1,mminus3,4 00314 t1=t1+A(i1,i2)*B(i1,i2) 00315 t2=t2+A(i1+1,i2)*B(i1+1,i2) 00316 t3=t3+A(i1+2,i2)*B(i1+2,i2) 00317 t4=t4+A(i1+3,i2)*B(i1+3,i2) 00318 ENDDO 00319 t1=t1+A(m,i2)*B(m,i2) 00320 ENDDO 00321 CASE (2) 00322 mminus3=m-3 00323 DO i2=1,n 00324 DO i1=1,mminus3,4 00325 t1=t1+A(i1,i2)*B(i1,i2) 00326 t2=t2+A(i1+1,i2)*B(i1+1,i2) 00327 t3=t3+A(i1+2,i2)*B(i1+2,i2) 00328 t4=t4+A(i1+3,i2)*B(i1+3,i2) 00329 ENDDO 00330 t1=t1+A(m-1,i2)*B(m-1,i2) 00331 t2=t2+A(m,i2)*B(m,i2) 00332 ENDDO 00333 CASE (3) 00334 mminus3=m-3 00335 DO i2=1,n 00336 DO i1=1,mminus3,4 00337 t1=t1+A(i1,i2)*B(i1,i2) 00338 t2=t2+A(i1+1,i2)*B(i1+1,i2) 00339 t3=t3+A(i1+2,i2)*B(i1+2,i2) 00340 t4=t4+A(i1+3,i2)*B(i1+3,i2) 00341 ENDDO 00342 t1=t1+A(m-2,i2)*B(m-2,i2) 00343 t2=t2+A(m-1,i2)*B(m-1,i2) 00344 t3=t3+A(m,i2)*B(m,i2) 00345 ENDDO 00346 END SELECT 00347 trace_r_AxB=t1+t2+t3+t4 00348 00349 END FUNCTION trace_r_AxB 00350 00351 ! ***************************************************************************** 00352 FUNCTION trace_r_AxB_new(A,lda,ia,ja,B,ldb,ib,jb,m,n) 00353 00354 INTEGER, INTENT(in) :: lda 00355 REAL(dp), INTENT(in) :: A(lda,*) 00356 INTEGER, INTENT(in) :: ia, ja, ldb 00357 REAL(dp), INTENT(in) :: B(ldb,*) 00358 INTEGER, INTENT(in) :: ib, jb, m, n 00359 REAL(dp) :: trace_r_AxB_new 00360 00361 INTEGER :: i1, i2, imod, mminus3 00362 REAL(dp) :: t1, t2, t3, t4 00363 00364 t1=0._dp 00365 t2=0._dp 00366 t3=0._dp 00367 t4=0._dp 00368 imod=MODULO(m,4) 00369 SELECT CASE (imod) 00370 CASE (0) 00371 DO i2=1,n 00372 DO i1=1,m,4 00373 t1=t1+A(i1+ia,i2+ja)*B(i1+ib,i2+jb) 00374 t2=t2+A(i1+ia+1,i2+ja)*B(i1+ib+1,i2+jb) 00375 t3=t3+A(i1+ia+2,i2+ja)*B(i1+ib+2,i2+jb) 00376 t4=t4+A(i1+ia+3,i2+ja)*B(i1+ib+3,i2+jb) 00377 ENDDO 00378 ENDDO 00379 CASE (1) 00380 mminus3=m-3 00381 DO i2=1,n 00382 DO i1=1,mminus3,4 00383 t1=t1+A(i1+ia,i2+ja)*B(i1+ib,i2+jb) 00384 t2=t2+A(i1+ia+1,i2+ja)*B(i1+ib+1,i2+jb) 00385 t3=t3+A(i1+ia+2,i2+ja)*B(i1+ib+2,i2+jb) 00386 t4=t4+A(i1+ia+3,i2+ja)*B(i1+ib+3,i2+jb) 00387 ENDDO 00388 t1=t1+A(m+ia,i2+ja)*B(m+ib,i2+jb) 00389 ENDDO 00390 CASE (2) 00391 mminus3=m-3 00392 DO i2=1,n 00393 DO i1=1,mminus3,4 00394 t1=t1+A(i1+ia,i2+ja)*B(i1+ib,i2+jb) 00395 t2=t2+A(i1+ia+1,i2+ja)*B(i1+ib+1,i2+jb) 00396 t3=t3+A(i1+ia+2,i2+ja)*B(i1+ib+2,i2+jb) 00397 t4=t4+A(i1+ia+3,i2+ja)*B(i1+ib+3,i2+jb) 00398 ENDDO 00399 t1=t1+A(m-1+ia,i2+ja)*B(m-1+ib,i2+jb) 00400 t2=t2+A(m+ia,i2+ja)*B(m+ib,i2+jb) 00401 ENDDO 00402 CASE (3) 00403 mminus3=m-3 00404 DO i2=1,n 00405 DO i1=1,mminus3,4 00406 t1=t1+A(i1+ia,i2+ja)*B(i1+ib,i2+jb) 00407 t2=t2+A(i1+ia+1,i2+ja)*B(i1+ib+1,i2+jb) 00408 t3=t3+A(i1+ia+2,i2+ja)*B(i1+ib+2,i2+jb) 00409 t4=t4+A(i1+ia+3,i2+ja)*B(i1+ib+3,i2+jb) 00410 ENDDO 00411 t1=t1+A(m-2+ia,i2+ja)*B(m-2+ib,i2+jb) 00412 t2=t2+A(m-1+ia,i2+ja)*B(m-1+ib,i2+jb) 00413 t3=t3+A(m+ia,i2+ja)*B(m+ib,i2+jb) 00414 ENDDO 00415 END SELECT 00416 trace_r_AxB_new=t1+t2+t3+t4 00417 00418 END FUNCTION trace_r_AxB_new 00419 00420 ! ***************************************************************************** 00421 SUBROUTINE transform_c2s(CPC_co,CPC_so,maxl,lm1,lm2) 00422 00423 REAL(dp), DIMENSION(:, :), INTENT(IN) :: CPC_co 00424 REAL(dp), DIMENSION(:, :), INTENT(OUT) :: CPC_so 00425 INTEGER :: maxl, lm1, lm2 00426 00427 CHARACTER(len=*), PARAMETER :: routineN = 'transform_c2s', 00428 routineP = moduleN//':'//routineN 00429 00430 INTEGER :: ic1, ic2, is1, is2, istat, l, 00431 lx, ly, lz 00432 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: work 00433 00434 ALLOCATE (work(ncoset(maxl),nsoset(maxl)),STAT=istat) 00435 IF (istat /= 0) THEN 00436 CALL stop_memory(routineN,moduleN,__LINE__,& 00437 "work",dp_size*ncoset(maxl)*nsoset(maxl)) 00438 END IF 00439 CPC_so = 0.0_dp 00440 work = 0.0_dp 00441 00442 ! DO l = 0,maxl 00443 ! ncgf = nco(l) 00444 ! nsgf = nso(l) 00445 ! ico = ncoset(l-1) + 1 00446 ! iso = nsoset(l-1) + 1 00447 ! 00448 ! CALL dgemm("N","T",ncoset(maxl),nsgf,ncgf,& 00449 ! 1.0_dp,CPC_co(1,ico),ncoset(maxl),& 00450 ! orbtramat(l)%c2s(1,1),nsgf,& 00451 ! 0.0_dp,work(1,iso),ncoset(maxl)) 00452 ! ENDDO 00453 ! 00454 ! DO l = 0,maxl 00455 ! ncgf = nco(l) 00456 ! nsgf = nso(l) 00457 ! ico = ncoset(l-1) + 1 00458 ! iso = nsoset(l-1) + 1 00459 ! 00460 ! CALL dgemm("N","N",nsgf,nsoset(maxl),ncgf,& 00461 ! 1.0_dp,orbtramat(l)%c2s(1,1),nsgf,& 00462 ! work(ico,1),ncoset(maxl),& 00463 ! 0.0_dp,CPC_so(iso,1),nsoset(maxl)) 00464 ! ENDDO 00465 00466 ! do l = 0,maxl 00467 ! do is1 = 1,nso(l) 00468 ! write(*,'(A,2I3,10f10.5)') 'or ', l, is1, orbtramat(l)%c2s(is1,1:nco(l)) 00469 ! enddo 00470 ! enddo 00471 ! 00472 ! stop 00473 00474 CPC_so = 0.0_dp 00475 work = 0.0_dp 00476 DO ic1 = 1,ncoset(lm1) 00477 DO l = 0,lm2 00478 DO is2 = 1,nso(l) 00479 DO ic2 = 1,nco(l) 00480 lx = indco(1,ic2+ncoset(l-1)) 00481 ly = indco(2,ic2+ncoset(l-1)) 00482 lz = indco(3,ic2+ncoset(l-1)) 00483 work(ic1,is2+nsoset(l-1)) = & 00484 work(ic1,is2+nsoset(l-1)) + & 00485 CPC_co(ic1,ic2+ncoset(l-1))*& 00486 orbtramat(l)%c2s(is2,ic2)*& 00487 SQRT(fourpi/dfac(2*l+1)*& 00488 dfac(2*lx-1)*dfac(2*ly-1)*dfac(2*lz-1)) 00489 ! write(*,*) 'dfac 1', dfac(2*lx-1)*dfac(2*ly-1)*dfac(2*lz-1) 00490 ENDDO 00491 ENDDO 00492 ENDDO 00493 ENDDO 00494 00495 DO is2 = 1,nsoset(lm2) 00496 DO l = 0,lm1 00497 DO is1 = 1,nso(l) 00498 DO ic1 = 1,nco(l) 00499 lx = indco(1,ic1+ncoset(l-1)) 00500 ly = indco(2,ic1+ncoset(l-1)) 00501 lz = indco(3,ic1+ncoset(l-1)) 00502 CPC_so(is1+nsoset(l-1),is2) = & 00503 CPC_so(is1+nsoset(l-1),is2) + & 00504 work(ic1+ncoset(l-1),is2)* & 00505 orbtramat(l)%c2s(is1,ic1)* & 00506 SQRT(fourpi/dfac(2*l+1)*& 00507 dfac(2*lx-1)*dfac(2*ly-1)*dfac(2*lz-1)) 00508 ! write(*,*) 'dfac 2', dfac(2*lx-1)*dfac(2*ly-1)*dfac(2*lz-1) 00509 ENDDO 00510 ENDDO 00511 ENDDO 00512 ENDDO 00513 00514 DEALLOCATE (work,STAT=istat) 00515 IF (istat /= 0) THEN 00516 CALL stop_memory(routineN,moduleN,__LINE__,"work") 00517 END IF 00518 00519 END SUBROUTINE transform_c2s 00520 00521 ! ***************************************************************************** 00522 SUBROUTINE transform_c2s_new(CPC_co,CPC_so,maxl) 00523 00524 REAL(dp), DIMENSION(:, :), INTENT(IN) :: CPC_co 00525 REAL(dp), DIMENSION(:, :), INTENT(OUT) :: CPC_so 00526 INTEGER :: maxl 00527 00528 CHARACTER(len=*), PARAMETER :: routineN = 'transform_c2s_new', 00529 routineP = moduleN//':'//routineN 00530 00531 INTEGER :: iso, istat, ldc, lds 00532 REAL(dp), ALLOCATABLE, 00533 DIMENSION(:, :, :) :: work 00534 00535 ALLOCATE (work(ncoset(maxl),ncoset(maxl),1),STAT=istat) 00536 IF (istat /= 0) THEN 00537 CALL stop_memory(routineN,moduleN,__LINE__,& 00538 "work",dp_size*ncoset(maxl)*nsoset(maxl)) 00539 END IF 00540 ldc = ncoset(maxl) 00541 lds = nsoset(maxl) 00542 CPC_so = 0.0_dp 00543 work = 0.0_dp 00544 00545 work(1:ldc,1:ldc,1) = CPC_co(1:ldc,1:ldc) 00546 00547 CALL cart2sph_mat(work,ldc,ldc,1,maxl,maxl) 00548 00549 DO iso = 1,nsoset(maxl) 00550 CPC_so(1:nsoset(maxl),iso) = work(1:nsoset(maxl),iso,1) 00551 END DO 00552 00553 DEALLOCATE (work,STAT=istat) 00554 IF (istat /= 0) THEN 00555 CALL stop_memory(routineN,moduleN,__LINE__,"work") 00556 END IF 00557 00558 END SUBROUTINE transform_c2s_new 00559 00560 ! ***************************************************************************** 00561 SUBROUTINE transform_s2c(matso,matco,maxl,lm1,lm2) 00562 00563 REAL(dp), DIMENSION(:, :), INTENT(IN) :: matso 00564 REAL(dp), DIMENSION(:, :), INTENT(OUT) :: matco 00565 INTEGER, INTENT(IN) :: maxl, lm1, lm2 00566 00567 CHARACTER(len=*), PARAMETER :: routineN = 'transform_s2c', 00568 routineP = moduleN//':'//routineN 00569 00570 INTEGER :: ic1, ic2, ico, is1, is2, iso, 00571 istat, l, lx, ly, lz, nc1, 00572 nc2, ns1, ns2 00573 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: matsc 00574 00575 ALLOCATE (matsc(nsoset(maxl),ncoset(maxl)),STAT=istat) 00576 IF (istat /= 0) THEN 00577 CALL stop_memory(routineN,moduleN,__LINE__,& 00578 "work",nsoset(maxl)*ncoset(maxl)) 00579 END IF 00580 matco = 0.0_dp 00581 matsc = 0.0_dp 00582 00583 ns1 = nsoset(lm1) 00584 ns2 = nsoset(lm2) 00585 nc1 = ncoset(lm1) 00586 nc2 = ncoset(lm2) 00587 00588 ! A = matso (ns1 x ns2) 00589 ! C = Tt x A x T = matco (nc1 x nc2) 00590 ! B = A x T = matsc (ns1 x nc2) 00591 00592 ! Calculate B 00593 DO is1 = 1,ns1 00594 DO l = 0,lm2 00595 DO ico = 1,nco(l) 00596 ic2 = ncoset(l-1) + ico 00597 lx = indco(1,ico+ncoset(l-1)) 00598 ly = indco(2,ico+ncoset(l-1)) 00599 lz = indco(3,ico+ncoset(l-1)) 00600 DO iso = 1,nso(l) 00601 is2 = nsoset(l-1) + iso 00602 matsc(is1,ic2) = matsc(is1,ic2) + & 00603 matso(is1,is2) * orbtramat(l)%s2c(iso,ico) *& 00604 SQRT((fourpi)/dfac(2*l+1)*& 00605 dfac(2*lx-1)*dfac(2*ly-1)*dfac(2*lz-1)) 00606 ENDDO ! iso 00607 ENDDO ! ico 00608 ENDDO ! l 00609 ENDDO ! is1 00610 00611 ! Calculate C 00612 DO ic2 = 1,nc2 00613 DO l = 0,lm1 00614 DO ico = 1,nco(l) 00615 ic1 = ncoset(l-1) + ico 00616 lx = indco(1,ico+ncoset(l-1)) 00617 ly = indco(2,ico+ncoset(l-1)) 00618 lz = indco(3,ico+ncoset(l-1)) 00619 DO iso = 1,nso(l) 00620 is1=nsoset(l-1) + iso 00621 matco(ic1,ic2) = matco(ic1,ic2) + & 00622 matsc(is1,ic2) * orbtramat(l)%s2c(iso,ico) *& 00623 SQRT((fourpi)/dfac(2*l+1)*& 00624 dfac(2*lx-1)*dfac(2*ly-1)*dfac(2*lz-1)) 00625 ENDDO ! iso 00626 ENDDO ! ico 00627 ENDDO ! l 00628 ENDDO ! ic2 00629 00630 DEALLOCATE (matsc,STAT=istat) 00631 IF (istat /= 0) THEN 00632 CALL stop_memory(routineN,moduleN,__LINE__,"matsc") 00633 END IF 00634 00635 END SUBROUTINE transform_s2c 00636 !------------------------------------------------------------------------------! 00637 ! ***************************************************************************** 00638 SUBROUTINE sph2cart_mat(mat,ld_mat,sd_mat,n,lmax1,lmax2) 00639 INTEGER :: ld_mat, sd_mat, n 00640 REAL(dp) :: mat(ld_mat,sd_mat,n) 00641 INTEGER :: lmax1, lmax2 00642 00643 CHARACTER(len=*), PARAMETER :: routineN = 'sph2cart_mat', 00644 routineP = moduleN//':'//routineN 00645 00646 INTEGER :: idx_l1, idx_l2 00647 REAL(dp) :: mat_aux(ncoset(lmax1),nsoset(lmax2),n) 00648 00649 IF (lmax1>3.or.lmax2>3) THEN 00650 CALL stop_program(routineN,moduleN,__LINE__,& 00651 "l>3 not implemented") 00652 END IF 00653 DO idx_l2=1,nsoset(lmax2) 00654 00655 mat_aux(1,idx_l2,:)& 00656 =s_root1o4pi*mat(1,idx_l2,:) 00657 00658 IF (lmax1==0) CYCLE 00659 00660 mat_aux(2,idx_l2,:)& 00661 =-s_root3o4pi*mat(4,idx_l2,:) 00662 00663 mat_aux(3,idx_l2,:)& 00664 =-s_root3o4pi*mat(2,idx_l2,:) 00665 00666 mat_aux(4,idx_l2,:)& 00667 =s_root3o4pi*mat(3,idx_l2,:) 00668 00669 IF (lmax1==1) CYCLE 00670 00671 mat_aux(5,idx_l2,:)& 00672 =-s_root5o16pi*mat(7,idx_l2,:)& 00673 +s_root15o16pi*mat(9,idx_l2,:) 00674 00675 mat_aux(6,idx_l2,:)& 00676 =s_root15o4pi*mat(5,idx_l2,:) 00677 00678 mat_aux(7,idx_l2,:)& 00679 =-s_root5o16pi*mat(7,idx_l2,:)& 00680 -s_root15o16pi*mat(9,idx_l2,:) 00681 00682 mat_aux(8,idx_l2,:)& 00683 =-s_root15o4pi*mat(8,idx_l2,:) 00684 00685 mat_aux(9,idx_l2,:)& 00686 =-s_root15o4pi*mat(6,idx_l2,:) 00687 00688 mat_aux(10,idx_l2,:)& 00689 =s_2root5o16pi*mat(7,idx_l2,:) 00690 00691 IF (lmax1==2) CYCLE 00692 00693 mat_aux(11,idx_l2,:)& 00694 =-s_root21o32pi*mat(14,idx_l2,:)& 00695 +s_root35o32pi*mat(16,idx_l2,:) 00696 00697 mat_aux(12,idx_l2,:)& 00698 =-s_3root35o32pi*mat(10,idx_l2,:)& 00699 +s_root21o32pi*mat(12,idx_l2,:) 00700 00701 mat_aux(13,idx_l2,:)& 00702 =s_root21o32pi*mat(14,idx_l2,:)& 00703 +s_3root35o32pi*mat(16,idx_l2,:) 00704 00705 mat_aux(14,idx_l2,:)& 00706 =s_root35o32pi*mat(10,idx_l2,:)& 00707 +s_root21o32pi*mat(12,idx_l2,:) 00708 00709 mat_aux(15,idx_l2,:)& 00710 =-s_3root7o16pi*mat(13,idx_l2,:)& 00711 +s_root105o16pi*mat(15,idx_l2,:) 00712 00713 mat_aux(16,idx_l2,:)& 00714 =s_root105o4pi*mat(11,idx_l2,:) 00715 00716 mat_aux(17,idx_l2,:)& 00717 =-s_3root7o16pi*mat(13,idx_l2,:)& 00718 -s_root105o16pi*mat(15,idx_l2,:) 00719 00720 mat_aux(18,idx_l2,:)& 00721 =-s_4root21o32pi*mat(14,idx_l2,:) 00722 00723 mat_aux(19,idx_l2,:)& 00724 =-s_4root21o32pi*mat(12,idx_l2,:) 00725 00726 mat_aux(20,idx_l2,:)& 00727 =s_2root7o16pi*mat(13,idx_l2,:) 00728 00729 ENDDO 00730 00731 DO idx_l1=1,ncoset(lmax1) 00732 00733 mat(idx_l1,1,:)& 00734 =s_root1o4pi*mat_aux(idx_l1,1,:) 00735 00736 IF (lmax2==0) CYCLE 00737 00738 mat(idx_l1,2,:)& 00739 =-s_root3o4pi*mat_aux(idx_l1,4,:) 00740 00741 mat(idx_l1,3,:)& 00742 =-s_root3o4pi*mat_aux(idx_l1,2,:) 00743 00744 mat(idx_l1,4,:)& 00745 =s_root3o4pi*mat_aux(idx_l1,3,:) 00746 00747 IF (lmax2==1) CYCLE 00748 00749 mat(idx_l1,5,:)& 00750 =-s_root5o16pi*mat_aux(idx_l1,7,:)& 00751 +s_root15o16pi*mat_aux(idx_l1,9,:) 00752 00753 mat(idx_l1,6,:)& 00754 =s_root15o4pi*mat_aux(idx_l1,5,:) 00755 00756 mat(idx_l1,7,:)& 00757 =-s_root5o16pi*mat_aux(idx_l1,7,:)& 00758 -s_root15o16pi*mat_aux(idx_l1,9,:) 00759 00760 mat(idx_l1,8,:)& 00761 =-s_root15o4pi*mat_aux(idx_l1,8,:) 00762 00763 mat(idx_l1,9,:)& 00764 =-s_root15o4pi*mat_aux(idx_l1,6,:) 00765 00766 mat(idx_l1,10,:)& 00767 =s_2root5o16pi*mat_aux(idx_l1,7,:) 00768 00769 IF (lmax2==2) CYCLE 00770 00771 mat(idx_l1,11,:)& 00772 =-s_root21o32pi*mat_aux(idx_l1,14,:)& 00773 +s_root35o32pi*mat_aux(idx_l1,16,:) 00774 00775 mat(idx_l1,12,:)& 00776 =-s_3root35o32pi*mat_aux(idx_l1,10,:)& 00777 +s_root21o32pi*mat_aux(idx_l1,12,:) 00778 00779 mat(idx_l1,13,:)& 00780 =s_root21o32pi*mat_aux(idx_l1,14,:)& 00781 +s_3root35o32pi*mat_aux(idx_l1,16,:) 00782 00783 mat(idx_l1,14,:)& 00784 =s_root35o32pi*mat_aux(idx_l1,10,:)& 00785 +s_root21o32pi*mat_aux(idx_l1,12,:) 00786 00787 mat(idx_l1,15,:)& 00788 =-s_3root7o16pi*mat_aux(idx_l1,13,:)& 00789 +s_root105o16pi*mat_aux(idx_l1,15,:) 00790 00791 mat(idx_l1,16,:)& 00792 =s_root105o4pi*mat_aux(idx_l1,11,:) 00793 00794 mat(idx_l1,17,:)& 00795 =-s_3root7o16pi*mat_aux(idx_l1,13,:)& 00796 -s_root105o16pi*mat_aux(idx_l1,15,:) 00797 00798 mat(idx_l1,18,:)& 00799 =-s_4root21o32pi*mat_aux(idx_l1,14,:) 00800 00801 mat(idx_l1,19,:)& 00802 =-s_4root21o32pi*mat_aux(idx_l1,12,:) 00803 00804 mat(idx_l1,20,:)& 00805 =s_2root7o16pi*mat_aux(idx_l1,13,:) 00806 00807 ENDDO 00808 00809 END SUBROUTINE sph2cart_mat 00810 00811 ! ***************************************************************************** 00812 SUBROUTINE cart2sph_mat(mat,ld_mat,sd_mat,n,lmax1,lmax2) 00813 00814 !in: 00815 INTEGER :: ld_mat, sd_mat, n 00816 REAL(dp) :: mat(ld_mat,sd_mat,n) 00817 INTEGER :: lmax1, lmax2 00818 00819 CHARACTER(len=*), PARAMETER :: routineN = 'cart2sph_mat', 00820 routineP = moduleN//':'//routineN 00821 00822 INTEGER :: idx_l2, idx_lm1 00823 REAL(dp) :: mat_aux(nsoset(lmax1),ncoset(lmax2),n) 00824 00825 IF ((lmax1 > 3).OR.(lmax2 > 3)) THEN 00826 CALL stop_program(routineN,moduleN,__LINE__,& 00827 "l > 3 not implemented") 00828 END IF 00829 00830 DO idx_l2=1,ncoset(lmax2) 00831 00832 mat_aux(1,idx_l2,:)& 00833 =root4pi*mat(1,idx_l2,:) 00834 00835 IF (lmax1==0) CYCLE 00836 00837 mat_aux(2,idx_l2,:)& 00838 =-root4pio3*mat(3,idx_l2,:) 00839 00840 mat_aux(3,idx_l2,:)& 00841 =root4pio3*mat(4,idx_l2,:) 00842 00843 mat_aux(4,idx_l2,:)& 00844 =-root4pio3*mat(2,idx_l2,:) 00845 00846 IF (lmax1==1) CYCLE 00847 00848 mat_aux(5,idx_l2,:)& 00849 =root4pio15*mat(6,idx_l2,:) 00850 00851 mat_aux(6,idx_l2,:)& 00852 =-root4pio15*mat(9,idx_l2,:) 00853 00854 mat_aux(7,idx_l2,:)& 00855 =-0.5_dp*root4pio5*mat(5,idx_l2,:)& 00856 -0.5_dp*root4pio5*mat(7,idx_l2,:)& 00857 +root4pio5*mat(10,idx_l2,:) 00858 00859 mat_aux(8,idx_l2,:)& 00860 =-root4pio15*mat(8,idx_l2,:) 00861 00862 mat_aux(9,idx_l2,:)& 00863 =0.5_dp*root4pio15*mat(5,idx_l2,:)& 00864 -0.5_dp*root4pio15*mat(7,idx_l2,:) 00865 00866 IF (lmax1==2) CYCLE 00867 00868 mat_aux(10,idx_l2,:)& 00869 =-s_3root35o32pi*mat(12,idx_l2,:)& 00870 +s_root35o32pi*mat(14,idx_l2,:) 00871 00872 mat_aux(11,idx_l2,:)& 00873 =s_root105o4pi*mat(16,idx_l2,:) 00874 00875 mat_aux(12,idx_l2,:)& 00876 =s_root21o32pi*mat(12,idx_l2,:)& 00877 +s_root21o32pi*mat(14,idx_l2,:)& 00878 -s_4root21o32pi*mat(19,idx_l2,:) 00879 00880 mat_aux(13,idx_l2,:)& 00881 =-s_3root7o16pi*mat(15,idx_l2,:)& 00882 -s_3root7o16pi*mat(17,idx_l2,:)& 00883 +s_2root7o16pi*mat(20,idx_l2,:) 00884 00885 mat_aux(14,idx_l2,:)& 00886 =s_root21o32pi*mat(11,idx_l2,:)& 00887 +s_root21o32pi*mat(13,idx_l2,:)& 00888 -s_4root21o32pi*mat(18,idx_l2,:) 00889 00890 mat_aux(15,idx_l2,:)& 00891 =s_root105o16pi*mat(15,idx_l2,:)& 00892 -s_root105o16pi*mat(17,idx_l2,:) 00893 00894 mat_aux(16,idx_l2,:)& 00895 =-s_root35o32pi*mat(11,idx_l2,:)& 00896 +s_3root35o32pi*mat(13,idx_l2,:) 00897 00898 ENDDO 00899 00900 DO idx_lm1=1,nsoset(lmax1) 00901 00902 mat(idx_lm1,1,:)& 00903 =root4pi*mat_aux(idx_lm1,1,:) 00904 00905 IF (lmax2==0) CYCLE 00906 00907 mat(idx_lm1,2,:)& 00908 =-root4pio3*mat_aux(idx_lm1,3,:) 00909 00910 mat(idx_lm1,3,:)& 00911 =root4pio3*mat_aux(idx_lm1,4,:) 00912 00913 mat(idx_lm1,4,:)& 00914 =-root4pio3*mat_aux(idx_lm1,2,:) 00915 00916 IF (lmax2==1) CYCLE 00917 00918 mat(idx_lm1,5,:)& 00919 =root4pio15*mat_aux(idx_lm1,6,:) 00920 00921 mat(idx_lm1,6,:)& 00922 =-root4pio15*mat_aux(idx_lm1,9,:) 00923 00924 mat(idx_lm1,7,:)& 00925 =-0.5_dp*root4pio5*mat_aux(idx_lm1,5,:)& 00926 -0.5_dp*root4pio5*mat_aux(idx_lm1,7,:)& 00927 +root4pio5*mat_aux(idx_lm1,10,:) 00928 00929 mat(idx_lm1,8,:)& 00930 =-root4pio15*mat_aux(idx_lm1,8,:) 00931 00932 mat(idx_lm1,9,:)& 00933 =0.5_dp*root4pio15*mat_aux(idx_lm1,5,:)& 00934 -0.5_dp*root4pio15*mat_aux(idx_lm1,7,:) 00935 00936 IF (lmax2==2) CYCLE 00937 00938 mat(idx_lm1,10,:)& 00939 =-s_3root35o32pi*mat_aux(idx_lm1,12,:)& 00940 +s_root35o32pi*mat_aux(idx_lm1,14,:) 00941 00942 mat(idx_lm1,11,:)& 00943 =s_root105o4pi*mat_aux(idx_lm1,16,:) 00944 00945 mat(idx_lm1,12,:)& 00946 =s_root21o32pi*mat_aux(idx_lm1,12,:)& 00947 +s_root21o32pi*mat_aux(idx_lm1,14,:)& 00948 -s_4root21o32pi*mat_aux(idx_lm1,19,:) 00949 00950 mat(idx_lm1,13,:)& 00951 =-s_3root7o16pi*mat_aux(idx_lm1,15,:)& 00952 -s_3root7o16pi*mat_aux(idx_lm1,17,:)& 00953 +s_2root7o16pi*mat_aux(idx_lm1,20,:) 00954 00955 mat(idx_lm1,14,:)& 00956 =s_root21o32pi*mat_aux(idx_lm1,11,:)& 00957 +s_root21o32pi*mat_aux(idx_lm1,13,:)& 00958 -s_4root21o32pi*mat_aux(idx_lm1,18,:) 00959 00960 mat(idx_lm1,15,:)& 00961 =s_root105o16pi*mat_aux(idx_lm1,15,:)& 00962 -s_root105o16pi*mat_aux(idx_lm1,17,:) 00963 00964 mat(idx_lm1,16,:)& 00965 =-s_root35o32pi*mat_aux(idx_lm1,11,:)& 00966 +s_3root35o32pi*mat_aux(idx_lm1,13,:) 00967 00968 ENDDO 00969 00970 END SUBROUTINE cart2sph_mat 00971 00972 END MODULE qs_util
1.7.3