CP2K: ai_fermi_contact.f90 Source File

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 ! *****************************************************************************
00027 MODULE ai_fermi_contact
00028   USE f77_blas
00029   USE kinds,                           ONLY: dp
00030   USE mathconstants,                   ONLY: pi
00031   USE orbital_pointers,                ONLY: coset,&
00032                                              ncoset
00033 #include "cp_common_uses.h"
00034 
00035   IMPLICIT NONE
00036 
00037   PRIVATE
00038 
00039   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_fermi_contact'
00040 
00041 ! *** Public subroutines ***
00042   PUBLIC :: fermi_contact
00043 
00044 CONTAINS
00045 
00046 ! *****************************************************************************
00053   SUBROUTINE fermi_contact(la_max,la_min,npgfa,rpgfa,zeta,&
00054                            lb_max,lb_min,npgfb,rpgfb,zetb,&
00055                            rac,rbc,dab,fcab,ldfc)
00056     INTEGER, INTENT(IN)                      :: la_max, la_min, npgfa
00057     REAL(KIND=dp), DIMENSION(:), INTENT(IN)  :: rpgfa, zeta
00058     INTEGER, INTENT(IN)                      :: lb_max, lb_min, npgfb
00059     REAL(KIND=dp), DIMENSION(:), INTENT(IN)  :: rpgfb, zetb
00060     REAL(KIND=dp), DIMENSION(3), INTENT(IN)  :: rac, rbc
00061     REAL(KIND=dp)                            :: dab
00062     INTEGER, INTENT(IN)                      :: ldfc
00063     REAL(KIND=dp), DIMENSION(ldfc, *), 
00064       INTENT(INOUT)                          :: fcab
00065 
00066     CHARACTER(len=*), PARAMETER :: routineN = 'fermi_contact', 
00067       routineP = moduleN//':'//routineN
00068 
00069     INTEGER                                  :: ax, ay, az, bx, by, bz, coa, 
00070                                                 cob, i, ipgf, j, jpgf, la, 
00071                                                 lb, ma, mb, na, nb
00072     REAL(KIND=dp)                            :: dac2, dbc2, f0, fax, fay, 
00073                                                 faz, fbx, fby, fbz
00074 
00075 ! *** Calculate some prefactors ***
00076 
00077     dac2 = rac(1)**2 + rac(2)**2 + rac(3)**2
00078     dbc2 = rbc(1)**2 + rbc(2)**2 + rbc(3)**2
00079 
00080     ! *** Loop over all pairs of primitive Gaussian-type functions ***
00081 
00082     na = 0
00083 
00084     DO ipgf=1,npgfa
00085 
00086        nb = 0
00087 
00088        DO jpgf=1,npgfb
00089 
00090           ! *** Screening ***
00091 
00092           IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) THEN
00093              DO j=nb+1,nb+ncoset(lb_max)
00094                 DO i=na+1,na+ncoset(la_max)
00095                    fcab(i,j) = 0.0_dp
00096                 ENDDO
00097              ENDDO
00098              nb = nb + ncoset(lb_max)
00099              CYCLE
00100           ENDIF
00101 
00102           ! *** Calculate some prefactors ***
00103 
00104           f0 = 4.0_dp/3.0_dp*pi*EXP(-zeta(ipgf)*dac2-zetb(jpgf)*dbc2)
00105 
00106           ! *** Calculate the primitive Fermi contact integrals ***
00107 
00108           DO lb = lb_min,lb_max
00109           DO bx = 0,lb
00110              fbx = 1.0_dp
00111              IF(bx.GT.0) fbx = (rbc(1))**bx
00112           DO by = 0,lb-bx
00113              bz = lb - bx - by
00114              cob = coset(bx,by,bz)
00115              mb = nb + cob
00116              fby = 1.0_dp
00117              IF(by.GT.0) fby = (rbc(2))**by
00118              fbz = 1.0_dp
00119              IF(bz.GT.0) fbz = (rbc(3))**bz
00120              DO la = la_min,la_max
00121              DO ax = 0,la
00122                 fax = fbx
00123                 IF(ax.GT.0) fax = fbx * (rac(1))**ax
00124              DO ay = 0,la-ax
00125                 az = la - ax - ay
00126                 coa = coset(ax,ay,az)
00127                 ma = na + coa
00128                 fay = fby
00129                 IF(ay.GT.0) fay = fby * (rac(2))**ay
00130                 faz = fbz
00131                 IF(az.GT.0) faz = fbz * (rac(3))**az
00132 
00133                 fcab(ma,mb) = f0 * fax * fay * faz
00134 
00135              ENDDO
00136              ENDDO
00137              ENDDO !la
00138 
00139           ENDDO
00140           ENDDO
00141           ENDDO !lb
00142 
00143           nb = nb + ncoset(lb_max)
00144 
00145        ENDDO
00146 
00147        na = na + ncoset(la_max)
00148 
00149     ENDDO
00150 
00151   END SUBROUTINE fermi_contact
00152 
00153 END MODULE ai_fermi_contact