1* 2* $Id$ 3* 4 subroutine integrate_kbpp_band_nonlocal(version,kvec, 5 > nrho,drho,lmax,locp,zv, 6 > vp,wp,rho,f,cs,sn, 7 > nfft1,nfft2,nfft3,lmmax, 8 > G,vnl, 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 zv 18 double precision vp(nrho,0:lmax) 19 double precision wp(nrho,0:lmax) 20 double precision rho(nrho) 21 double precision f(nrho) 22 double precision cs(nrho) 23 double precision sn(nrho) 24 25 integer nfft1,nfft2,nfft3,lmmax 26 double precision G(nfft1,nfft2,nfft3,3) 27 double precision vnl(nfft1,nfft2,nfft3,lmmax) 28 29 integer ierr 30 31 integer np,taskid,MASTER 32 parameter (MASTER=0) 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#include "bafdecls.fh" 47 48 call Parallel_np(np) 49 call Parallel_taskid(taskid) 50 51 nfft3d = (nfft1)*nfft2*nfft3 52 pi=4.0d0*datan(1.0d0) 53 twopi=2.0d0*pi 54 forpi=4.0d0*pi 55 56 IF(LMMAX.GT.16) THEN 57 WRITE(*,*)"non-local psp not generated: lmmax exceeds 16" 58 IERR=1 59 RETURN 60 ENDIF 61 IF((NRHO/2)*2.EQ.NRHO) THEN 62 WRITE(*,*)"non-local psp not generated: nrho is not odd" 63 IERR=2 64 RETURN 65 ENDIF 66 67 P0=DSQRT(FORPI) 68 P1=DSQRT(3.0d0*FORPI) 69 P2=DSQRT(15.0d0*FORPI) 70 P3=DSQRT(105.0d0*FORPI) 71 72 73*====================== Fourier transformation ====================== 74 call dcopy(lmmax*nfft3d,0.0d0,0,VNL,1) 75 76 task_count = -1 77 DO 700 k3=1,nfft3 78 DO 700 k2=1,nfft2 79 DO 700 k1=1,nfft1 80 task_count = task_count + 1 81 if (mod(task_count,np).ne.taskid) go to 700 82 gx=G(k1,k2,k3,1)+kvec(1) 83 gy=G(k1,k2,k3,2)+kvec(2) 84 gz=G(k1,k2,k3,3)+kvec(3) 85 86 Q=DSQRT(gx**2 + gy**2 + gz**2) 87 88 if (dabs(Q).gt.1.0d-9) then 89 90 gx=gx/Q 91 gy=gy/Q 92 gz=gz/Q 93 DO I=1,NRHO 94 CS(I)=DCOS(Q*RHO(I)) 95 SN(I)=DSIN(Q*RHO(I)) 96 END DO 97 98 lcount = lmmax+1 99 GO TO (500,400,300,200), LMAX+1 100 101 102*:::::::::::::::::::::::::::::: f-wave :::::::::::::::::::::::::::::: 103 200 CONTINUE 104 if (locp.ne.3) then 105 F(1)=0.0d0 106 do I=2,NRHO 107 A=SN(I)/(Q*RHO(I)) 108 A=15.0d0*(A-CS(I))/(Q*RHO(I))**2 - 6*A + CS(I) 109 F(I)=A*WP(I,3)*VP(I,3) 110 end do 111 D=P3*SIMP(NRHO,F,DRHO)/Q 112 lcount = lcount-1 113 VNL(k1,k2,k3,lcount)=D*GX*(4.0d0*GX*GX-3.0d0*(1.0d0-GZ*GZ)) 114 > /dsqrt(24.0d0) 115 lcount = lcount-1 116 VNL(k1,k2,k3,lcount)=D*GY*(3.0d0*(1.0d0-GZ*GZ)-4.0d0*GY*GY) 117 > /dsqrt(24.0d0) 118 lcount = lcount-1 119 VNL(k1,k2,k3,lcount)=D*GZ*(GX*GX - GY*GY) 120 > /2.0d0 121 lcount = lcount-1 122 VNL(k1,k2,k3,lcount)=D*GX*GY*GZ 123 lcount = lcount-1 124 VNL(k1,k2,k3,lcount)=D*GX*(5.0d0*GZ*GZ-1.0d0) 125 > /dsqrt(40.0d0) 126 lcount = lcount-1 127 VNL(k1,k2,k3,lcount)=D*GY*(5.0d0*GZ*GZ-1.0d0) 128 > /dsqrt(40.0d0) 129 lcount = lcount-1 130 VNL(k1,k2,k3,lcount)=D*GZ*(5.0d0*GZ*GZ-3.0d0) 131 > /dsqrt(60.0d0) 132 end if 133 134 135 136*:::::::::::::::::::::::::::::: d-wave :::::::::::::::::::::::::::::: 137 300 CONTINUE 138 if (locp.ne.2) then 139 F(1)=0.0d0 140 DO I=2,NRHO 141 A=3.0d0*(SN(I)/(Q*RHO(I))-CS(I))/(Q*RHO(I))-SN(I) 142 F(I)=A*WP(I,2)*VP(I,2) 143 END DO 144 D=P2*SIMP(NRHO,F,DRHO)/Q 145 lcount = lcount-1 146 VNL(k1,k2,k3,lcount)=D*(3.0d0*GZ*GZ-1.0d0) 147 > /(2.0d0*dsqrt(3.0d0)) 148c VNL(k1,k2,k3,lcount)=D*(2.0d0*GZ*GZ-GX*GX-GY*GY) 149c > /(2.0d0*dsqrt(3.0d0)) 150 lcount = lcount-1 151 VNL(k1,k2,k3,lcount)=D*GX*GY 152 lcount = lcount-1 153 VNL(k1,k2,k3,lcount)=D*GY*GZ 154 lcount = lcount-1 155 VNL(k1,k2,k3,lcount)=D*GZ*GX 156 lcount = lcount-1 157 VNL(k1,k2,k3,lcount)=D*(GX*GX-GY*GY)/(2.0d0) 158 end if 159 160*:::::::::::::::::::::::::::::: p-wave :::::::::::::::::::::::::::::: 161 400 CONTINUE 162 if (locp.ne.1) then 163 F(1)=0.0d0 164 DO I=2,NRHO 165 F(I)=(SN(I)/(Q*RHO(I))-CS(I))*WP(I,1)*VP(I,1) 166 END DO 167 P=P1*SIMP(NRHO,F,DRHO)/Q 168 lcount = lcount-1 169 VNL(k1,k2,k3,lcount)=P*GX 170 lcount = lcount-1 171 VNL(k1,k2,k3,lcount)=P*GY 172 lcount = lcount-1 173 VNL(k1,k2,k3,lcount)=P*GZ 174 end if 175 176*:::::::::::::::::::::::::::::: s-wave ::::::::::::::::::::::::::::::: 177 500 CONTINUE 178 if (locp.ne.0) then 179 DO I=1,NRHO 180 F(I)=SN(I)*WP(I,0)*VP(I,0) 181 END DO 182 lcount = lcount-1 183 VNL(k1,k2,k3,lcount)=P0*SIMP(NRHO,F,DRHO)/Q 184 end if 185 186 600 CONTINUE 187 188*::::::::::::::::::::::::::::::: G+k=0 :::::::::::::::::::::::::::::::: 189 else 190 191 do l=1,lmmax 192 vnl(k1,k2,k3,l)=0.0d0 193 end do 194* *** only j0 is non-zero at zero **** 195 if (locp.ne.0) then 196 DO I=1,NRHO 197 F(I)=RHO(I)*WP(I,0)*VP(I,0) 198 END DO 199 VNL(k1,k2,k3,1)=P0*SIMP(NRHO,F,DRHO) 200 end if 201 202 end if 203 204 205 700 CONTINUE 206 207 call Parallel_Vector_SumAll(lmmax*nfft3d,VNL) 208 209 210 IERR=0 211 RETURN 212 END 213 214 215 216