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