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