1      subroutine gridQpqPotential(zora_Qpq,      ! in : type of EFG potential
2     &                            xyz_EFGcoords, ! in : EFG-nuclear coord.
3     &                            nqpts,         ! in : nr. grid points
4     &                            qxyz,          ! in : grid points
5     &                            amat_Qnucl,    ! out: EFG potential
6     &                            closegridpts)
7c    Purpose: Evaluates Q_pq(\vec{r},\vec{r}')
8c    Input:
9c           zora_Qpq = 1  -> Evaluates Qxx
10c                    = 2  -> Evaluates Qyy
11c                    = 3  -> Evaluates Qzz
12c                    = 4  -> Evaluates Qxy
13c                    = 5  -> Evaluates Qxz
14c                    = 6  -> Evaluates Qyz
15c           nqpts    , number of grid points
16c           qxyz     , grid points
17c           xyz_EFGcoords, Quadrupole potential is evaluated
18c                          in this point (\vec{r})
19c    Output:
20c           amat_Qnucl, Quadrupole potential ev. in the grid
21c                       for integration purpose (\vec{r}')
22c    Author: Fredy Aquino
23c    Date  : 10-29-10
24
25      implicit none
26#include "global.fh"
27#include "stdio.fh"
28#include "zora.fh"
29      integer igrid,nqpts
30      integer closegridpts(*)
31      integer zora_Qpq
32      double precision xyz_EFGcoords(3)
33      double precision qxyz(3,nqpts)
34      double precision rx,ry,rz,dist,dist5
35      double precision amat_Qnucl(nqpts)
36
37c     == loop over the grid points ==
38      if (zora_Qpq.eq.1) then        ! Qxx
39       do igrid = 1,nqpts
40c     == distance from the grid points to given xyz_EFGcoords() ==
41        rx = xyz_EFGcoords(1) - qxyz(1,igrid)
42        ry = xyz_EFGcoords(2) - qxyz(2,igrid)
43        rz = xyz_EFGcoords(3) - qxyz(3,igrid)
44        dist = dsqrt(rx*rx + ry*ry + rz*rz)
45        if (dist.gt.zoracutoff_EFG) then ! check-cutoff
46          dist5=dist*dist*dist*dist*dist
47          amat_Qnucl(igrid) = (2.d0*rx*rx-(ry*ry+rz*rz))/dist5
48        else
49          closegridpts(igrid) = igrid
50        end if
51       end do ! end-grid
52      else if (zora_Qpq.eq.2) then    ! Qyy
53       do igrid = 1,nqpts
54c     == distance from the grid points to given xyz_EFGcoords() ==
55        rx = xyz_EFGcoords(1) - qxyz(1,igrid)
56        ry = xyz_EFGcoords(2) - qxyz(2,igrid)
57        rz = xyz_EFGcoords(3) - qxyz(3,igrid)
58        dist = dsqrt(rx*rx + ry*ry + rz*rz)
59        if (dist.gt.zoracutoff_EFG) then ! check-cutoff
60          dist5=dist*dist*dist*dist*dist
61         amat_Qnucl(igrid) = (2.d0*ry*ry-(rx*rx+rz*rz))/dist5
62        else
63          closegridpts(igrid) = igrid
64        end if
65       end do ! end-grid
66      else if (zora_Qpq.eq.3) then    ! Qzz
67       do igrid = 1,nqpts
68c     == distance from the grid points to given xyz_EFGcoords() ==
69        rx = xyz_EFGcoords(1) - qxyz(1,igrid)
70        ry = xyz_EFGcoords(2) - qxyz(2,igrid)
71        rz = xyz_EFGcoords(3) - qxyz(3,igrid)
72        dist = dsqrt(rx*rx + ry*ry + rz*rz)
73        if (dist.gt.zoracutoff_EFG) then ! check-cutoff
74          dist5=dist*dist*dist*dist*dist
75         amat_Qnucl(igrid) = (2.d0*rz*rz-(rx*rx+ry*ry))/dist5
76        else
77          closegridpts(igrid) = igrid
78        end if
79       end do ! end-grid
80      else if (zora_Qpq.eq.4) then    ! Qxy
81       do igrid = 1,nqpts
82c     == distance from the grid points to given xyz_EFGcoords() ==
83        rx = xyz_EFGcoords(1) - qxyz(1,igrid)
84        ry = xyz_EFGcoords(2) - qxyz(2,igrid)
85        rz = xyz_EFGcoords(3) - qxyz(3,igrid)
86        dist = dsqrt(rx*rx + ry*ry + rz*rz)
87        if (dist.gt.zoracutoff_EFG) then ! check-cutoff
88          dist5=dist*dist*dist*dist*dist
89          amat_Qnucl(igrid) = (3.d0*rx*ry)/dist5
90        else
91          closegridpts(igrid) = igrid
92        end if
93       end do ! end-grid
94      else if (zora_Qpq.eq.5) then    ! Qxz
95       do igrid = 1,nqpts
96c     == distance from the grid points to given xyz_EFGcoords() ==
97        rx = xyz_EFGcoords(1) - qxyz(1,igrid)
98        ry = xyz_EFGcoords(2) - qxyz(2,igrid)
99        rz = xyz_EFGcoords(3) - qxyz(3,igrid)
100        dist = dsqrt(rx*rx + ry*ry + rz*rz)
101        if (dist.gt.zoracutoff_EFG) then ! check-cutoff
102          dist5=dist*dist*dist*dist*dist
103          amat_Qnucl(igrid) = (3.d0*rx*rz)/dist5
104        else
105          closegridpts(igrid) = igrid
106        end if
107       end do ! end-grid
108      else if (zora_Qpq.eq.6) then    ! Qyz
109       do igrid = 1,nqpts
110c     == distance from the grid points to given xyz_EFGcoords() ==
111        rx = xyz_EFGcoords(1) - qxyz(1,igrid)
112        ry = xyz_EFGcoords(2) - qxyz(2,igrid)
113        rz = xyz_EFGcoords(3) - qxyz(3,igrid)
114        dist = dsqrt(rx*rx + ry*ry + rz*rz)
115        if (dist.gt.zoracutoff_EFG) then ! check-cutoff
116          dist5=dist*dist*dist*dist*dist
117          amat_Qnucl(igrid) = (3.d0*ry*rz)/dist5
118        else
119          closegridpts(igrid) = igrid
120        end if
121       end do ! end-grid
122      endif
123      return
124      end
125c $Id$
126