1      SUBROUTINE ewald_recp(ewald2)
2
3      implicit none
4
5      include 'p_input.inc'
6      include 'p_array.inc'
7      include 'p_const.inc'
8      include 'cm_atom.inc'
9      include 'cm_ewld.inc'
10      include 'cm_latt.inc'
11
12      integer i,iatm,k,kx,ky,kz,kminx,kminy,kminz
13
14      real*8 rksq,rx,ry,rz,rkx,rky,rkz
15      real*8 kcoeff,factor,force,ewald2
16
17      double complex rhosum
18      double complex eikr(mxatms)
19      double complex eikx(1:mxatms,       0:mxkvect)
20      double complex eiky(1:mxatms,-mxkvect:mxkvect)
21      double complex eikz(1:mxatms,-mxkvect:mxkvect)
22
23      do i=1,natms
24
25       eikx(i,0)=(1.0,0.0)
26       eiky(i,0)=(1.0,0.0)
27       eikz(i,0)=(1.0,0.0)
28       rx=rlatt(1,1)*ccc(i,1)+rlatt(1,2)*ccc(i,2)+rlatt(1,3)*ccc(i,3)
29       ry=rlatt(2,1)*ccc(i,1)+rlatt(2,2)*ccc(i,2)+rlatt(2,3)*ccc(i,3)
30       rz=rlatt(3,1)*ccc(i,1)+rlatt(3,2)*ccc(i,2)+rlatt(3,3)*ccc(i,3)
31       eikx(i,1)=dcmplx(dcos(twopi*rx),dsin(twopi*rx))
32       eiky(i,1)=dcmplx(dcos(twopi*ry),dsin(twopi*ry))
33       eikz(i,1)=dcmplx(dcos(twopi*rz),dsin(twopi*rz))
34       eiky(i,-1)=conjg(eiky(i,1))
35       eikz(i,-1)=conjg(eikz(i,1))
36
37      enddo
38
39      do i=1,natms
40
41       do k=2,kmaxx
42        eikx(i,k)=eikx(i,k-1)*eikx(i,1)
43       enddo
44       do k=2,kmaxy
45        eiky(i,k)=eiky(i,k-1)*eiky(i,1)
46        eiky(i,-k)=conjg(eiky(i,k))
47       enddo
48       do k=2,kmaxz
49        eikz(i,k)=eikz(i,k-1)*eikz(i,1)
50        eikz(i,-k)=conjg(eikz(i,k))
51       enddo
52
53      enddo
54
55      kminx=0
56      kminy=-kmaxy
57      kminz=-kmaxz
58
59      do kx=kminx,kmaxx
60
61       if(kx.eq.0)then
62        factor=1.0
63       else
64        factor=2.0
65       endif
66
67       do ky=kminy,kmaxy
68
69        do kz=kminz,kmaxz
70
71         rkx=real(kx)*rlatt(1,1)+real(ky)*rlatt(1,2)+real(kz)*rlatt(1,3)
72         rky=real(kx)*rlatt(2,1)+real(ky)*rlatt(2,2)+real(kz)*rlatt(2,3)
73         rkz=real(kx)*rlatt(3,1)+real(ky)*rlatt(3,2)+real(kz)*rlatt(3,3)
74         rkx=twopi*rkx
75         rky=twopi*rky
76         rkz=twopi*rkz
77         rksq=rkx*rkx+rky*rky+rkz*rkz
78
79          if(rksq.lt.rksqmax.and.rksq.ne.0.0)then
80
81           rhosum=(0.0,0.0)
82
83           do i=1,natms
84
85            iatm=atmtype(i)
86            eikr(i)=eikx(i,kx)*eiky(i,ky)*eikz(i,kz)
87            rhosum=rhosum+typchge(iatm)*eikr(i)
88
89           enddo
90
91           kcoeff=exp(rksq*ralphsq)/rksq
92           ewald2=ewald2+factor*kcoeff*conjg(rhosum)*rhosum
93
94           do i=1,natms
95
96            iatm=atmtype(i)
97            force=-factor*2.0*twopi*convfct1/vol*kcoeff*
98     $            dimag(rhosum*dconjg(eikr(i)))*typchge(iatm)
99            fff(i,1)=fff(i,1)+convfct2*rkx*force
100            fff(i,2)=fff(i,2)+convfct2*rky*force
101            fff(i,3)=fff(i,3)+convfct2*rkz*force
102
103           enddo
104
105          endif
106
107         enddo
108
109        enddo
110
111       enddo
112
113       ewald2=twopi*ewald2*convfct1/vol
114
115       return
116
117       END
118c $Id$
119