1* 2* $Id$ 3* 4 subroutine integrate_kbpp_band_orb(version,kvec, 5 > nrho,drho,lmax,locp, 6 > wp,rho,f,cs,sn, 7 > nfft1,nfft2,nfft3,lmmax, 8 > G,borbs, 9 > ierr) 10 implicit none 11 integer version 12 double precision kvec(3) 13 integer nrho 14 double precision drho 15 integer lmax 16 integer locp 17 double precision wp(nrho,0:lmax) 18 double precision rho(nrho) 19 double precision f(nrho) 20 double precision cs(nrho) 21 double precision sn(nrho) 22 23 integer nfft1,nfft2,nfft3,lmmax 24 double precision G(nfft1,nfft2,nfft3,3) 25 double precision borbs(nfft1,nfft2,nfft3,lmmax) 26 27 integer ierr 28 29 integer np,taskid,MASTER 30 parameter (MASTER=0) 31 32#include "bafdecls.fh" 33 34* *** local variables **** 35 integer lcount,task_count,nfft3d 36 integer k1,k2,k3,i,l 37 double precision pi,twopi,forpi 38 double precision p0,p1,p2,p3,p 39 double precision gx,gy,gz,a,q,d 40 41 42* **** external functions **** 43 double precision dsum,simp 44 external dsum,simp 45 logical value 46 47 call Parallel_np(np) 48 call Parallel_taskid(taskid) 49 50 nfft3d = (nfft1)*nfft2*nfft3 51 pi=4.0d0*datan(1.0d0) 52 twopi=2.0d0*pi 53 forpi=4.0d0*pi 54 55 IF(LMMAX.GT.16) THEN 56 WRITE(*,*)"non-local psp not generated: lmmax exceeds 16" 57 IERR=1 58 RETURN 59 ENDIF 60 IF((NRHO/2)*2.EQ.NRHO) THEN 61 WRITE(*,*)"non-local psp not generated: nrho is not odd" 62 IERR=2 63 RETURN 64 ENDIF 65 66 P0=DSQRT(FORPI) 67 P1=DSQRT(3.0d0*FORPI) 68 P2=DSQRT(15.0d0*FORPI) 69 P3=DSQRT(105.0d0*FORPI) 70 71 72*====================== Fourier transformation ====================== 73 call dcopy(lmmax*nfft3d,0.0d0,0,borbs,1) 74 75 task_count = -1 76 DO 700 k3=1,nfft3 77 DO 700 k2=1,nfft2 78 DO 700 k1=1,nfft1 79 task_count = task_count + 1 80 if (mod(task_count,np).ne.taskid) go to 700 81 gx=G(k1,k2,k3,1)+kvec(1) 82 gy=G(k1,k2,k3,2)+kvec(2) 83 gz=G(k1,k2,k3,3)+kvec(3) 84 85 Q=DSQRT(gx**2 + gy**2 + gz**2) 86 87 if (dabs(Q).gt.1.0d-9) then 88 89 gx=gx/Q 90 gy=gy/Q 91 gz=gz/Q 92 DO I=1,NRHO 93 CS(I)=DCOS(Q*RHO(I)) 94 SN(I)=DSIN(Q*RHO(I)) 95 END DO 96 97 lcount = lmmax+1 98 GO TO (500,400,300,200,100), LMAX+1 99 100 101*:::::::::::::::::::::::::::::: g-wave :::::::::::::::::::::::::::::: 102 100 CONTINUE 103*:::::::::::::::::::::::::::::: f-wave :::::::::::::::::::::::::::::: 104 200 CONTINUE 105 if (locp.ne.3) then 106 F(1)=0.0d0 107 do I=2,NRHO 108 A=SN(I)/(Q*RHO(I)) 109 A=15.0d0*(A-CS(I))/(Q*RHO(I))**2 - 6*A + CS(I) 110 F(I)=A*WP(I,3) 111 end do 112 D=P3*SIMP(NRHO,F,DRHO)/Q 113 114* **** fy(3x2-y2) component **** 115 lcount = lcount-1 116 borbs(k1,k2,k3,lcount)=D*GY*(3.0d0*(1.0d0-GZ*GZ)-4.0d0*GY*GY) 117 > /dsqrt(24.0d0) 118 119* **** fxyz component **** 120 lcount = lcount-1 121 borbs(k1,k2,k3,lcount)=D*GX*GY*GZ 122 123* **** fy(5z2-1) component **** 124 lcount = lcount-1 125 borbs(k1,k2,k3,lcount)=D*GY*(5.0d0*GZ*GZ-1.0d0) 126 > /dsqrt(40.0d0) 127 128* **** fz(5z2-3) component **** 129 lcount = lcount-1 130 borbs(k1,k2,k3,lcount)=D*GZ*(5.0d0*GZ*GZ-3.0d0) 131 > /dsqrt(60.0d0) 132 133* **** fx(5z2-1) component **** 134 lcount = lcount-1 135 borbs(k1,k2,k3,lcount)=D*GX*(5.0d0*GZ*GZ-1.0d0) 136 > /dsqrt(40.0d0) 137 138* **** fz(x2-y2) component **** 139 lcount = lcount-1 140 borbs(k1,k2,k3,lcount)=D*GZ*(GX*GX - GY*GY) 141 > /2.0d0 142 143* **** fx(x2-3y2) component **** 144 lcount = lcount-1 145 borbs(k1,k2,k3,lcount)=D*GX*(4.0d0*GX*GX-3.0d0*(1.0d0-GZ*GZ)) 146 > /dsqrt(24.0d0) 147 end if 148 149 150*:::::::::::::::::::::::::::::: d-wave :::::::::::::::::::::::::::::: 151 300 CONTINUE 152 if (locp.ne.2) then 153 F(1)=0.0d0 154 DO I=2,NRHO 155 A=3.0d0*(SN(I)/(Q*RHO(I))-CS(I))/(Q*RHO(I))-SN(I) 156 F(I)=A*WP(I,2) 157 END DO 158 D=P2*SIMP(NRHO,F,DRHO)/Q 159 160* **** dxy component **** 161 lcount = lcount-1 162 borbs(k1,k2,k3,lcount)=D*GX*GY 163 164* **** dyz component **** 165 lcount = lcount-1 166 borbs(k1,k2,k3,lcount)=D*GY*GZ 167 168* **** d3z2-1 component **** 169 lcount = lcount-1 170 borbs(k1,k2,k3,lcount)=D*(3.0d0*GZ*GZ-1.0d0) 171 > /(2.0d0*dsqrt(3.0d0)) 172 173* **** dzx component **** 174 lcount = lcount-1 175 borbs(k1,k2,k3,lcount)=D*GZ*GX 176 177* **** dx2-y2 component **** 178 lcount = lcount-1 179 borbs(k1,k2,k3,lcount)=D*(GX*GX-GY*GY)/(2.0d0) 180 181 end if 182 183*:::::::::::::::::::::::::::::: p-wave :::::::::::::::::::::::::::::: 184 400 CONTINUE 185 if (locp.ne.1) then 186 F(1)=0.0d0 187 DO I=2,NRHO 188 F(I)=(SN(I)/(Q*RHO(I))-CS(I))*WP(I,1) 189 END DO 190 P=P1*SIMP(NRHO,F,DRHO)/Q 191 192* **** py component **** 193 lcount = lcount-1 194 borbs(k1,k2,k3,lcount)=P*GY 195 196* **** pz component **** 197 lcount = lcount-1 198 borbs(k1,k2,k3,lcount)=P*GZ 199 200* **** px component **** 201 lcount = lcount-1 202 borbs(k1,k2,k3,lcount)=P*GX 203 end if 204 205*:::::::::::::::::::::::::::::: s-wave ::::::::::::::::::::::::::::::: 206 500 CONTINUE 207 if (locp.ne.0) then 208 DO I=1,NRHO 209 F(I)=SN(I)*WP(I,0) 210 END DO 211 lcount = lcount-1 212 borbs(k1,k2,k3,lcount)=P0*SIMP(NRHO,F,DRHO)/Q 213 end if 214 215 600 CONTINUE 216 217*::::::::::::::::::::::::::::::: G+k=0 :::::::::::::::::::::::::::::::: 218 else 219 220 do l=1,lmmax 221 borbs(k1,k2,k3,l)=0.0d0 222 end do 223* *** only j0 is non-zero at zero **** 224 if (locp.ne.0) then 225 DO I=1,NRHO 226 F(I)=RHO(I)*WP(I,0) 227 END DO 228 borbs(k1,k2,k3,1)=P0*SIMP(NRHO,F,DRHO) 229 end if 230 231 end if 232 233 234 700 CONTINUE 235 236 call Parallel_Vector_SumAll(lmmax*nfft3d,borbs) 237 238 IERR=0 239 RETURN 240 END 241 242 243 244