1
2      subroutine smd_ewald_recip_generic(
3     >                          na,
4     >                          nk,
5     >                          eikr,
6     >                          eikx,
7     >                          eiky,
8     >                          eikz,
9     >                          ralphsq,
10     >                          rksqmax,
11     >                          vol,
12     >                          rlatt,
13     >                          kmax,
14     >                          ccc,
15     >                          q,
16     >                          fff,
17     >                          ewald2)
18
19      implicit none
20#include "smd_const_data.fh"
21
22      integer na
23      integer nk
24      double precision ralphsq
25      double precision rksqmax
26      double precision vol
27      double complex eikr(na)
28      double complex eikx(1:na,0:nk)
29      double complex eiky(1:na,-nk:nk)
30      double complex eikz(1:na,-nk:nk)
31      double precision rlatt(3,3)
32      double precision ccc(na,3)
33      double precision fff(na,3)
34      double precision q(na)
35      integer kmax(3)
36
37
38
39      integer i,k,kx,ky,kz,kminx,kminy,kminz
40
41      real*8 rksq,rx,ry,rz,rkx,rky,rkz
42      real*8 kcoeff,factor,force,ewald2
43
44      double complex rhosum
45
46      ewald2 = 0.0d0
47      do i=1,na
48
49       eikx(i,0)=(1.0,0.0)
50       eiky(i,0)=(1.0,0.0)
51       eikz(i,0)=(1.0,0.0)
52       rx=rlatt(1,1)*ccc(i,1)+rlatt(1,2)*ccc(i,2)+rlatt(1,3)*ccc(i,3)
53       ry=rlatt(2,1)*ccc(i,1)+rlatt(2,2)*ccc(i,2)+rlatt(2,3)*ccc(i,3)
54       rz=rlatt(3,1)*ccc(i,1)+rlatt(3,2)*ccc(i,2)+rlatt(3,3)*ccc(i,3)
55       eikx(i,1)=dcmplx(dcos(twopi*rx),dsin(twopi*rx))
56       eiky(i,1)=dcmplx(dcos(twopi*ry),dsin(twopi*ry))
57       eikz(i,1)=dcmplx(dcos(twopi*rz),dsin(twopi*rz))
58       eiky(i,-1)=conjg(eiky(i,1))
59       eikz(i,-1)=conjg(eikz(i,1))
60
61      enddo
62
63      do i=1,na
64
65       do k=2,kmax(1)
66        eikx(i,k)=eikx(i,k-1)*eikx(i,1)
67       enddo
68       do k=2,kmax(2)
69        eiky(i,k)=eiky(i,k-1)*eiky(i,1)
70        eiky(i,-k)=conjg(eiky(i,k))
71       enddo
72       do k=2,kmax(3)
73        eikz(i,k)=eikz(i,k-1)*eikz(i,1)
74        eikz(i,-k)=conjg(eikz(i,k))
75       enddo
76
77      enddo
78
79      kminx=0
80      kminy=-kmax(2)
81      kminz=-kmax(3)
82
83      do kx=kminx,kmax(1)
84
85       if(kx.eq.0)then
86        factor=1.0
87       else
88        factor=2.0
89       endif
90
91       do ky=kminy,kmax(2)
92
93        do kz=kminz,kmax(3)
94
95         rkx=real(kx)*rlatt(1,1)+real(ky)*rlatt(1,2)+real(kz)*rlatt(1,3)
96         rky=real(kx)*rlatt(2,1)+real(ky)*rlatt(2,2)+real(kz)*rlatt(2,3)
97         rkz=real(kx)*rlatt(3,1)+real(ky)*rlatt(3,2)+real(kz)*rlatt(3,3)
98         rkx=twopi*rkx
99         rky=twopi*rky
100         rkz=twopi*rkz
101         rksq=rkx*rkx+rky*rky+rkz*rkz
102
103          if(rksq.lt.rksqmax.and.rksq.ne.0.0)then
104
105           rhosum=(0.0,0.0)
106
107           do i=1,na
108
109            eikr(i)=eikx(i,kx)*eiky(i,ky)*eikz(i,kz)
110            rhosum=rhosum+q(i)*eikr(i)
111
112           enddo
113
114           kcoeff=exp(rksq*ralphsq)/rksq
115           ewald2=ewald2+factor*kcoeff*conjg(rhosum)*rhosum
116           do i=1,na
117
118            force=-factor*2.0*twopi*convfct1/vol*kcoeff*
119     $            dimag(rhosum*dconjg(eikr(i)))*q(i)
120            fff(i,1)=fff(i,1)+convfct2*rkx*force
121            fff(i,2)=fff(i,2)+convfct2*rky*force
122            fff(i,3)=fff(i,3)+convfct2*rkz*force
123
124           enddo
125
126          endif
127
128         enddo
129
130        enddo
131
132       enddo
133
134       ewald2=twopi*ewald2*convfct1/vol
135
136       return
137
138       END
139
140
141      subroutine smd_ewald_excl_generic(na,
142     >                                  nl,
143     >                                  alpha,
144     >                                  rcutsq,
145     >                                  latt,
146     >                                  rlatt,
147     >                                  q,
148     >                                  ccc,
149     >                                  fff,
150     >                                  epoint,
151     >                                  elist,
152     >                                  e)
153
154      implicit none
155
156#include "smd_const_data.fh"
157
158
159      integer na
160      integer nl
161      double precision alpha
162      double precision rcutsq
163      double precision rlatt(3,3),latt(3,3)
164      integer epoint(na)
165      double precision q(na)
166      integer elist(nl)
167      double precision ccc(na,3),fff(na,3)
168      double precision e
169c
170      integer i,j,k,jnab
171      integer jbeg,jend
172
173      double precision  dr,ar,rsq
174      double precision  erfxc,force
175
176      double precision x,y,z
177
178      double precision  ssx,ssy,ssz,xss,yss,zss
179
180      e=0
181
182      do i=1,na-1
183
184       jbeg=epoint(i)
185       jend=epoint(i+1)-1
186
187c       write(*,*) "i,jbeg,jend",i,jbeg,jend
188
189      do jnab=jbeg,jend
190
191       j=elist(jnab)
192
193       x=ccc(i,1)-ccc(j,1)
194       y=ccc(i,2)-ccc(j,2)
195       z=ccc(i,3)-ccc(j,3)
196c
197c      reboxing here
198c      ------------
199       ssx=(rlatt(1,1)*x+rlatt(1,2)*y+rlatt(1,3)*z)
200       ssy=(rlatt(2,1)*x+rlatt(2,2)*y+rlatt(2,3)*z)
201       ssz=(rlatt(3,1)*x+rlatt(3,2)*y+rlatt(3,3)*z)
202
203       xss=ssx-nint(ssx)
204       yss=ssy-nint(ssy)
205       zss=ssz-nint(ssz)
206
207       x=(latt(1,1)*xss+latt(1,2)*yss+latt(1,3)*zss)
208       y=(latt(2,1)*xss+latt(2,2)*yss+latt(2,3)*zss)
209       z=(latt(3,1)*xss+latt(3,2)*yss+latt(3,3)*zss)
210c      done reboxing
211
212       rsq=x*x+y*y+z*z
213       if(rsq.lt.rcutsq)then
214
215        dr=sqrt(rsq)
216        ar=alpha*dr
217
218        e=e-convfct1*q(i)*q(j)
219     $       *(1-erfxc(ar))/dr
220
221        force=-convfct1*q(i)*q(j)*
222     $       ((1-erfxc(ar))-2*ar/sqrpi*exp(-ar*ar))
223     $       /(dr*rsq)
224
225        fff(i,1)=fff(i,1)+convfct2*force*x
226        fff(i,2)=fff(i,2)+convfct2*force*y
227        fff(i,3)=fff(i,3)+convfct2*force*z
228
229        fff(j,1)=fff(j,1)-convfct2*force*x
230        fff(j,2)=fff(j,2)-convfct2*force*y
231        fff(j,3)=fff(j,3)-convfct2*force*z
232
233       endif
234
235      end do
236      end do
237
238      return
239
240      END
241
242      subroutine smd_ewald_real_generic(na,
243     >                                  nl,
244     >                                  alpha,
245     >                                  rcutsq,
246     >                                  q,
247     >                                  ccc,
248     >                                  fff,
249     >                                  point,
250     >                                  list,
251     >                                  e)
252
253      implicit none
254
255#include "smd_const_data.fh"
256
257
258      integer na
259      integer nl
260      double precision alpha
261      double precision rcutsq
262      integer point(na)
263      double precision q(na)
264      integer list(nl)
265      double precision ccc(nl,3)
266      double precision fff(na,3)
267      double precision e
268c
269      integer i,j,k,jnab
270      integer jbeg,jend
271      integer nlist
272
273      double precision  dr,ar,rsq
274      double precision  erfxc,force
275
276      double precision x,y,z
277
278      e=0
279
280      nlist = 0
281      do i=1,na-1
282
283       jbeg=point(i)
284       jend=point(i+1)-1
285
286
287      do jnab=jbeg,jend
288
289       j=list(jnab)
290
291       nlist = nlist + 1
292       x=ccc(nlist,1)
293       y=ccc(nlist,2)
294       z=ccc(nlist,3)
295
296       rsq=x*x+y*y+z*z
297
298       if(rsq.lt.rcutsq)then
299
300        dr=sqrt(rsq)
301        ar=alpha*dr
302
303        e=e+convfct1*q(i)*q(j)
304     $       *erfxc(ar)/dr
305c        write(*,*) q(i),q(j),ar,dr
306
307        force=convfct1*q(i)*q(j)
308     $       *(erfxc(ar)+2*ar/sqrpi*exp(-ar*ar))/(dr*rsq)
309
310        fff(i,1)=fff(i,1)+convfct2*force*x
311        fff(i,2)=fff(i,2)+convfct2*force*y
312        fff(i,3)=fff(i,3)+convfct2*force*z
313
314        fff(j,1)=fff(j,1)-convfct2*force*x
315        fff(j,2)=fff(j,2)-convfct2*force*y
316        fff(j,3)=fff(j,3)-convfct2*force*z
317
318       endif
319
320      end do
321      end do
322
323      return
324
325      END
326
327c $Id$
328