PROGRESS  master
prg_ewald_mod.F90
Go to the documentation of this file.
1 ! Ewald sum routines for kernel calculation
2 module prg_ewald_mod
3 
4  use bml
5  use prg_timer_mod
6  use prg_parallel_mod
7 
8  implicit none
9 
10  private !Everything is private by default
11 
12  integer, parameter :: dp = kind(1.0d0)
13 
14  public :: ewald_real_space_single
15  public :: ewald_real_space_single_latte
16  public :: ewald_real_space
17  public :: ewald_real_space_latte
18  public :: ewald_real_space_test
19  public :: ewald_real_space_matrix_latte
20  public :: ewald_k_space_latte_single
21  public :: ewald_k_space_latte
22  public :: ewald_k_space_test
23  public :: ewald_k_space_matrix_latte
24 
25 contains
26 
28  subroutine ewald_real_space_single_latte(COULOMBV,I,RXYZ,Box,Nr_elem, &
29  DELTAQ,J,U,Element_Pointer,Nr_atoms,COULACC,HDIM,Max_Nr_Neigh)
30 
31  implicit none
32 
33  integer, parameter :: prec = 8
34  real(prec), parameter :: one = 1.d0, two = 2.d0, zero = 0.d0, six = 6.d0, three = 3.d0
35  real(prec), parameter :: fourtyeight = 48.d0, eleven = 11.d0, sixteen = 16.d0
36  real(prec), parameter :: pi = 3.14159265358979323846264d0
37  real(prec), parameter :: sqrtpi = 1.772453850905516d0
38  integer, intent(in) :: nr_atoms, nr_elem, hdim, max_nr_neigh, i, j, element_pointer(nr_atoms)
39  real(prec), intent(in) :: coulacc, deltaq(nr_atoms)
40  real(prec) :: tfact, relperm, keconst
41  real(prec), intent(in) :: rxyz(3,nr_atoms), box(3,3)
42  real(prec), intent(in) :: u(nr_elem)
43  real(prec) :: coulcut, coulcut2
44  real(prec), intent(out) :: coulombv
45  real(prec) :: ra(3), rb(3), dr, rab(3)
46  real(prec) :: coulvol, sqrtx, calpha, dc(3), magr, magr2, z
47  real(prec) :: calpha2, ti,ti2,ti3,ti4,ti6,ssa,ssb,ssc,ssd,sse
48  real(prec) :: numrep_erfc, ca, force, expti, exptj
49  real(prec) :: tj,tj2,tj3,tj4,tj6,ti2mtj2a,sa,sb,sc,sd,se,sf
50  real(prec) :: ti2mtj2, ti2mti2, tj2mti2
51 
52  integer :: k, ccnt,l,m,n
53 
54  coulvol = box(1,1)*box(2,2)*box(3,3)
55  sqrtx = sqrt(-log(coulacc))
56 
57  ccnt = 0
58  coulcut = 12.d0
59  calpha = sqrtx/coulcut
60  coulcut2 = coulcut*coulcut
61  calpha2 = calpha*calpha
62 
63  relperm = one
64  keconst = 14.3996437701414d0*relperm
65  tfact = 16.0d0/(5.0d0*keconst)
66 
67  coulombv = zero
68 
69  ti = tfact*u(element_pointer(i))
70  ti2 = ti*ti
71  ti3 = ti2*ti
72  ti4 = ti2*ti2
73  ti6 = ti4*ti2
74 
75  ssa = ti
76  ssb = ti3/48.d0
77  ssc = 3.d0*ti2/16.d0
78  ssd = 11.d0*ti/16.d0
79  sse = 1.d0
80 
81  ra(1) = rxyz(1,i)
82  ra(2) = rxyz(2,i)
83  ra(3) = rxyz(3,i)
84 
85  do k = -1,1
86  do m = -1,1
87  do l = -1,1
88 
89  rb(1) = rxyz(1,j)+k*box(1,1)
90  rb(2) = rxyz(2,j)+m*box(2,2)
91  rb(3) = rxyz(3,j)+l*box(3,3)
92  rab = rb-ra ! OBS b - a !!!
93  dr = norm2(rab)
94  magr = dr
95  magr2 = dr*dr
96 
97  if ((dr <= coulcut).and.(dr > 1e-12)) then
98 
99  tj = tfact*u(element_pointer(j))
100  dc = rab/dr
101 
102  z = abs(calpha*magr)
103  numrep_erfc = erfc(z)
104 
105  ca = numrep_erfc/magr
106  coulombv = coulombv + deltaq(j)*ca
107  ccnt = ccnt + 1
108  ca = ca + two*calpha*exp( -calpha2*magr2 )/sqrtpi
109  expti = exp(-ti*magr )
110 
111  if (element_pointer(i).eq.element_pointer(j)) then
112  coulombv = coulombv - deltaq(j)*expti*(ssb*magr2 + ssc*magr + ssd + sse/magr)
113  ccnt = ccnt + 1
114  else
115  tj2 = tj*tj
116  tj3 = tj2*tj
117  tj4 = tj2*tj2
118  tj6 = tj4*tj2
119  exptj = exp( -tj*magr )
120  ti2mtj2 = ti2 - tj2
121  tj2mti2 = -ti2mtj2
122  sa = ti
123  sb = tj4*ti/(two*ti2mtj2*ti2mtj2)
124  sc = (tj6 - three*tj4*ti2)/(ti2mtj2*ti2mtj2*ti2mtj2)
125  sd = tj
126  se = ti4*tj/(two * tj2mti2 * tj2mti2)
127  sf = (ti6 - three*ti4*tj2)/(tj2mti2*tj2mti2*tj2mti2)
128 
129  coulombv = coulombv - (deltaq(j)*(expti*(sb - (sc/magr)) + exptj*(se - (sf/magr))))
130  endif
131  endif
132  enddo
133  enddo
134  enddo
135 
136  coulombv = keconst*coulombv
137 
138  end subroutine ewald_real_space_single_latte
139 
140  subroutine ewald_real_space_single(COULOMBV,FCOUL,I,RX,RY,RZ,LBox, &
141  DELTAQ,J,U,Element_Type,Nr_atoms,COULACC,TIMERATIO,HDIM,Max_Nr_Neigh)
142 
143 
144  implicit none
145 
146  integer, parameter :: prec = 8
147  real(prec), parameter :: one = 1.d0, two = 2.d0, zero = 0.d0, six = 6.d0, three = 3.d0
148  real(prec), parameter :: fourtyeight = 48.d0, eleven = 11.d0, sixteen = 16.d0
149  real(prec), parameter :: pi = 3.14159265358979323846264d0
150  real(prec), parameter :: sqrtpi = 1.772453850905516d0
151  integer, intent(in) :: nr_atoms, hdim, max_nr_neigh, i, j
152  real(prec), intent(in) :: coulacc, timeratio,deltaq(nr_atoms)
153  real(prec) :: tfact, relperm, keconst
154  real(prec), intent(in) :: rx(nr_atoms), ry(nr_atoms), rz(nr_atoms), lbox(3)
155  real(prec), intent(in) :: u(nr_atoms)
156  real(prec) :: coulcut, coulcut2
157  character(10), intent(in) :: element_type(nr_atoms)
158  real(prec), intent(out) :: coulombv, fcoul(3)
159  real(prec) :: ra(3), rb(3), dr, rab(3)
160  real(prec) :: coulvol, sqrtx, calpha, dc(3), magr, magr2, z
161  real(prec) :: calpha2, ti,ti2,ti3,ti4,ti6,ssa,ssb,ssc,ssd,sse
162  real(prec) :: numrep_erfc, ca, force, expti, exptj
163  real(prec) :: tj,tj2,tj3,tj4,tj6,ti2mtj2a,sa,sb,sc,sd,se,sf
164  real(prec) :: ti2mtj2, ti2mti2, tj2mti2
165 
166  integer :: k, ccnt,l,m,n
167 
168  coulvol = lbox(1)*lbox(2)*lbox(3)
169  sqrtx = sqrt(-log(coulacc))
170 
171  ccnt = 0
172  coulcut = 12.d0
173  calpha = sqrtx/coulcut
174  coulcut2 = coulcut*coulcut
175  calpha2 = calpha*calpha
176 
177  relperm = one
178  keconst = 14.3996437701414d0*relperm
179  tfact = 16.0d0/(5.0d0*keconst)
180 
181  fcoul = zero
182  coulombv = zero
183 
184  ti = tfact*u(i)
185  ti2 = ti*ti
186  ti3 = ti2*ti
187  ti4 = ti2*ti2
188  ti6 = ti4*ti2
189 
190  ssa = ti
191  ssb = ti3/48.d0
192  ssc = 3.d0*ti2/16.d0
193  ssd = 11.d0*ti/16.d0
194  sse = 1.d0
195 
196  ra(1) = rx(i)
197  ra(2) = ry(i)
198  ra(3) = rz(i)
199 
200  do k = -1,1
201  do m = -1,1
202  do l = -1,1
203 
204  rb(1) = rx(j)+k*lbox(1)
205  rb(2) = ry(j)+m*lbox(2)
206  rb(3) = rz(j)+l*lbox(3)
207  rab = rb-ra ! OBS b - a !!!
208  dr = norm2(rab)
209  magr = dr
210  magr2 = dr*dr
211 
212  if ((dr <= coulcut).and.(dr > 1e-12)) then
213 
214  tj = tfact*u(j)
215  dc = rab/dr
216 
217  z = abs(calpha*magr)
218  numrep_erfc = erfc(z)
219 
220  ca = numrep_erfc/magr
221  coulombv = coulombv + deltaq(j)*ca
222  ccnt = ccnt + 1
223  ca = ca + two*calpha*exp( -calpha2*magr2 )/sqrtpi
224  force = -keconst*deltaq(i)*deltaq(j)*ca/magr
225  expti = exp(-ti*magr )
226 
227  if (element_type(i).eq.element_type(j)) then
228  coulombv = coulombv - deltaq(j)*expti*(ssb*magr2 + ssc*magr + ssd + sse/magr)
229  ccnt = ccnt + 1
230  force = force + (keconst*deltaq(i)*deltaq(j)*expti)*((sse/magr2 - two*ssb*magr - ssc) &
231  + ssa*(ssb*magr2 + ssc*magr + ssd + sse/magr))
232  else
233  tj2 = tj*tj
234  tj3 = tj2*tj
235  tj4 = tj2*tj2
236  tj6 = tj4*tj2
237  exptj = exp( -tj*magr )
238  ti2mtj2 = ti2 - tj2
239  tj2mti2 = -ti2mtj2
240  sa = ti
241  sb = tj4*ti/(two*ti2mtj2*ti2mtj2)
242  sc = (tj6 - three*tj4*ti2)/(ti2mtj2*ti2mtj2*ti2mtj2)
243  sd = tj
244  se = ti4*tj/(two * tj2mti2 * tj2mti2)
245  sf = (ti6 - three*ti4*tj2)/(tj2mti2*tj2mti2*tj2mti2)
246 
247  coulombv = coulombv - (deltaq(j)*(expti*(sb - (sc/magr)) + exptj*(se - (sf/magr))))
248  force = force + keconst*deltaq(i)*deltaq(j)*((expti*(sa*(sb - (sc/magr)) - (sc/magr2))) &
249  + (exptj*(sd*(se - (sf/magr)) - (sf/magr2))))
250  endif
251 
252  fcoul(1) = fcoul(1) + dc(1)*force
253  fcoul(2) = fcoul(2) + dc(2)*force
254  fcoul(3) = fcoul(3) + dc(3)*force
255  endif
256  enddo
257  enddo
258  enddo
259 
260  coulombv = keconst*coulombv
261 
262  end subroutine ewald_real_space_single
263 
264  subroutine ewald_real_space_matrix_latte(E,RXYZ,Box,U,Element_Pointer,Nr_atoms,COULACC,nebcoul,totnebcoul,HDIM,Max_Nr_Neigh,Nr_Elem)
265 
266  implicit none
267 
268  integer, parameter :: prec = 8
269  real(prec), parameter :: one = 1.d0, two = 2.d0, zero = 0.d0, six = 6.d0, three = 3.d0
270  real(prec), parameter :: fourtyeight = 48.d0, eleven = 11.d0, sixteen = 16.d0
271  real(prec), parameter :: pi = 3.14159265358979323846264d0
272  real(prec), parameter :: sqrtpi = 1.772453850905516d0
273  integer, intent(in) :: nr_atoms, hdim, max_nr_neigh, nr_elem
274  real(prec), intent(in) :: coulacc
275  real(prec), intent(out) :: e(nr_atoms,nr_atoms)
276  real(prec) :: tfact, relperm, keconst
277  real(prec), intent(in) :: rxyz(3,nr_atoms), box(3,3)
278  real(prec), intent(in) :: u(nr_elem)
279  real(prec) :: coulcut, coulcut2
280  integer, intent(in) :: element_pointer(nr_atoms)
281  integer, intent(in) :: totnebcoul(nr_atoms), nebcoul(4,max_nr_neigh,nr_atoms)
282  real(prec) :: ra(3), rb(3), dr, rab(3)
283  real(prec) :: coulvol, sqrtx, calpha, dc(3), magr, magr2, z
284  real(prec) :: calpha2, ti,ti2,ti3,ti4,ti6,ssa,ssb,ssc,ssd,sse
285  real(prec) :: numrep_erfc, ca, force, expti, exptj
286  real(prec) :: tj,tj2,tj3,tj4,tj6,ti2mtj2a,sa,sb,sc,sd,se,sf
287  real(prec) :: ti2mtj2, ti2mti2, tj2mti2
288  integer :: i,j,k, ccnt, newj, pbci,pbcj,pbck
289 
290  coulvol = box(1,1)*box(2,2)*box(3,3)
291  sqrtx = sqrt(-log(coulacc))
292 
293  ccnt = 0
294  coulcut = 12.d0
295  calpha = sqrtx/coulcut
296  coulcut2 = coulcut*coulcut
297  calpha2 = calpha*calpha
298 
299  relperm = one
300  keconst = 14.3996437701414d0*relperm
301  tfact = 16.0d0/(5.0d0*keconst)
302 
303  e = 0.0
304 
305  !$OMP PARALLEL DO DEFAULT(PRIVATE) SHARED(E,U,Element_Pointer,RXYZ,totnebcoul,nebcoul,coulcut,calpha)
306 
307  do i = 1,nr_atoms
308 
309  ti = tfact*u(element_pointer(i))
310  ti2 = ti*ti
311  ti3 = ti2*ti
312  ti4 = ti2*ti2
313  ti6 = ti4*ti2
314 
315  ssa = ti
316  ssb = ti3/48.d0
317  ssc = 3.d0*ti2/16.d0
318  ssd = 11.d0*ti/16.d0
319  sse = 1.d0
320 
321  ra(1) = rxyz(1,i)
322  ra(2) = rxyz(2,i)
323  ra(3) = rxyz(3,i)
324 
325  do newj = 1,totnebcoul(i)
326  j = nebcoul(1, newj, i)
327  pbci = nebcoul(2, newj, i)
328  pbcj = nebcoul(3, newj, i)
329  pbck = nebcoul(4, newj, i)
330  rb(1) = rxyz(1,j) + real(pbci)*box(1,1) + real(pbcj)*box(2,1) + &
331  REAL(pbck)*box(3,1)
332 
333  rb(2) = rxyz(2,j) + real(pbci)*box(1,2) + real(pbcj)*box(2,2) + &
334  REAL(pbck)*box(3,2)
335 
336  rb(3) = rxyz(3,j) + real(pbci)*box(1,3) + real(pbcj)*box(2,3) + &
337  REAL(pbck)*box(3,3)
338  rab = rb-ra ! OBS b - a !!!
339  dr = norm2(rab)
340  magr = dr
341  magr2 = dr*dr
342 
343  if ((dr <= coulcut).and.(dr > 1e-12)) then
344 
345  tj = tfact*u(element_pointer(j))
346  dc = rab/dr
347 
348  z = abs(calpha*magr)
349  numrep_erfc = erfc(z)
350 
351  ca = numrep_erfc/magr
352  e(i,j) = e(i,j) + ca
353  ccnt = ccnt + 1
354  ca = ca + two*calpha*exp( -calpha2*magr2 )/sqrtpi
355  expti = exp(-ti*magr )
356 
357  if (element_pointer(i).eq.element_pointer(j)) then
358  e(i,j) = e(i,j) - expti*(ssb*magr2 + ssc*magr + ssd + sse/magr)
359  ccnt = ccnt + 1
360  else
361  tj2 = tj*tj
362  tj3 = tj2*tj
363  tj4 = tj2*tj2
364  tj6 = tj4*tj2
365  exptj = exp( -tj*magr )
366  ti2mtj2 = ti2 - tj2
367  tj2mti2 = -ti2mtj2
368  sa = ti
369  sb = tj4*ti/(two*ti2mtj2*ti2mtj2)
370  sc = (tj6 - three*tj4*ti2)/(ti2mtj2*ti2mtj2*ti2mtj2)
371  sd = tj
372  se = ti4*tj/(two * tj2mti2 * tj2mti2)
373  sf = (ti6 - three*ti4*tj2)/(tj2mti2*tj2mti2*tj2mti2)
374 
375  e(i,j) = e(i,j) - ((expti*(sb - (sc/magr)) + exptj*(se - (sf/magr))))
376  endif
377 
378  endif
379  enddo
380  enddo
381  !$OMP END PARALLEL DO
382 
383  e = keconst*e
384 
385  end subroutine ewald_real_space_matrix_latte
386 
387  subroutine ewald_real_space_latte(COULOMBV,I,RXYZ,Box, &
388  DELTAQ,U,Element_Pointer,Nr_atoms,COULACC,nebcoul,totnebcoul,HDIM,Max_Nr_Neigh,Nr_Elem)
389 
390  implicit none
391 
392  integer, parameter :: prec = 8
393  real(prec), parameter :: one = 1.d0, two = 2.d0, zero = 0.d0, six = 6.d0, three = 3.d0
394  real(prec), parameter :: fourtyeight = 48.d0, eleven = 11.d0, sixteen = 16.d0
395  real(prec), parameter :: pi = 3.14159265358979323846264d0
396  real(prec), parameter :: sqrtpi = 1.772453850905516d0
397  integer, intent(in) :: nr_atoms, hdim, max_nr_neigh, i, nr_elem
398  real(prec), intent(in) :: coulacc
399  real(prec) :: tfact, relperm, keconst
400  real(prec), intent(in) :: rxyz(3,nr_atoms), box(3,3), deltaq(nr_atoms)
401  real(prec), intent(in) :: u(nr_elem)
402  real(prec) :: coulcut, coulcut2
403  integer, intent(in) :: element_pointer(nr_atoms)
404  integer, intent(in) :: totnebcoul(nr_atoms), nebcoul(4,max_nr_neigh,nr_atoms)
405  real(prec), intent(out) :: coulombv
406  real(prec) :: ra(3), rb(3), dr, rab(3)
407  real(prec) :: coulvol, sqrtx, calpha, dc(3), magr, magr2, z
408  real(prec) :: calpha2, ti,ti2,ti3,ti4,ti6,ssa,ssb,ssc,ssd,sse
409  real(prec) :: numrep_erfc, ca, force, expti, exptj
410  real(prec) :: tj,tj2,tj3,tj4,tj6,ti2mtj2a,sa,sb,sc,sd,se,sf
411  real(prec) :: ti2mtj2, ti2mti2, tj2mti2
412  integer :: j,k, ccnt, newj, pbci,pbcj,pbck
413 
414  coulvol = box(1,1)*box(2,2)*box(3,3)
415  sqrtx = sqrt(-log(coulacc))
416 
417  ccnt = 0
418  coulcut = 12.d0
419  calpha = sqrtx/coulcut
420  coulcut2 = coulcut*coulcut
421  calpha2 = calpha*calpha
422 
423  relperm = one
424  keconst = 14.3996437701414d0*relperm
425  tfact = 16.0d0/(5.0d0*keconst)
426 
427  coulombv = zero
428 
429  ti = tfact*u(element_pointer(i))
430  ti2 = ti*ti
431  ti3 = ti2*ti
432  ti4 = ti2*ti2
433  ti6 = ti4*ti2
434 
435  ssa = ti
436  ssb = ti3/48.d0
437  ssc = 3.d0*ti2/16.d0
438  ssd = 11.d0*ti/16.d0
439  sse = 1.d0
440 
441  ra(1) = rxyz(1,i)
442  ra(2) = rxyz(2,i)
443  ra(3) = rxyz(3,i)
444 
445  do newj = 1,totnebcoul(i)
446  j = nebcoul(1, newj, i)
447  pbci = nebcoul(2, newj, i)
448  pbcj = nebcoul(3, newj, i)
449  pbck = nebcoul(4, newj, i)
450  rb(1) = rxyz(1,j) + real(pbci)*box(1,1) + real(pbcj)*box(2,1) + &
451  REAL(pbck)*box(3,1)
452 
453  rb(2) = rxyz(2,j) + real(pbci)*box(1,2) + real(pbcj)*box(2,2) + &
454  REAL(pbck)*box(3,2)
455 
456  rb(3) = rxyz(3,j) + real(pbci)*box(1,3) + real(pbcj)*box(2,3) + &
457  REAL(pbck)*box(3,3)
458  rab = rb-ra ! OBS b - a !!!
459  dr = norm2(rab)
460  magr = dr
461  magr2 = dr*dr
462 
463  if ((dr <= coulcut).and.(dr > 1e-12)) then
464 
465  tj = tfact*u(element_pointer(j))
466  dc = rab/dr
467 
468  z = abs(calpha*magr)
469  numrep_erfc = erfc(z)
470 
471  ca = numrep_erfc/magr
472  coulombv = coulombv + deltaq(j)*ca
473  ccnt = ccnt + 1
474  ca = ca + two*calpha*exp( -calpha2*magr2 )/sqrtpi
475  expti = exp(-ti*magr )
476 
477  if (element_pointer(i).eq.element_pointer(j)) then
478  coulombv = coulombv - deltaq(j)*expti*(ssb*magr2 + ssc*magr + ssd + sse/magr)
479  ccnt = ccnt + 1
480  else
481  tj2 = tj*tj
482  tj3 = tj2*tj
483  tj4 = tj2*tj2
484  tj6 = tj4*tj2
485  exptj = exp( -tj*magr )
486  ti2mtj2 = ti2 - tj2
487  tj2mti2 = -ti2mtj2
488  sa = ti
489  sb = tj4*ti/(two*ti2mtj2*ti2mtj2)
490  sc = (tj6 - three*tj4*ti2)/(ti2mtj2*ti2mtj2*ti2mtj2)
491  sd = tj
492  se = ti4*tj/(two * tj2mti2 * tj2mti2)
493  sf = (ti6 - three*ti4*tj2)/(tj2mti2*tj2mti2*tj2mti2)
494 
495  coulombv = coulombv - (deltaq(j)*(expti*(sb - (sc/magr)) + exptj*(se - (sf/magr))))
496  endif
497 
498  endif
499  enddo
500  coulombv = keconst*coulombv
501 
502  end subroutine ewald_real_space_latte
503 
504  subroutine ewald_real_space_test(COULOMBV,I,RX,RY,RZ,LBox, &
505  DELTAQ,U,Element_Type,Nr_atoms,COULACC,nnRx,nnRy,nnRz,nrnnlist,nnType,Max_Nr_Neigh)
506 
507  implicit none
508 
509  integer, parameter :: prec = 8
510  real(prec), parameter :: one = 1.d0, two = 2.d0, zero = 0.d0, six = 6.d0, three = 3.d0
511  real(prec), parameter :: fourtyeight = 48.d0, eleven = 11.d0, sixteen = 16.d0
512  real(prec), parameter :: pi = 3.14159265358979323846264d0
513  real(prec), parameter :: sqrtpi = 1.772453850905516d0
514  integer, intent(in) :: nr_atoms, max_nr_neigh, i
515  real(prec), intent(in) :: coulacc
516  real(prec) :: tfact, relperm, keconst
517  real(prec), intent(in) :: rx(nr_atoms), ry(nr_atoms), rz(nr_atoms), lbox(3), deltaq(nr_atoms)
518  real(prec), intent(in) :: u(nr_atoms)
519  real(prec) :: coulcut, coulcut2
520  character(10), intent(in) :: element_type(nr_atoms)
521  integer, intent(in) :: nrnnlist(nr_atoms), nntype(nr_atoms,max_nr_neigh)
522  real(prec), intent(in) :: nnrx(nr_atoms,max_nr_neigh)
523  real(prec), intent(in) :: nnry(nr_atoms,max_nr_neigh), nnrz(nr_atoms,max_nr_neigh)
524  real(prec), intent(out) :: coulombv
525  real(prec) :: ra(3), rb(3), dr, rab(3)
526  real(prec) :: coulvol, sqrtx, calpha, dc(3), magr, magr2, z
527  real(prec) :: calpha2, ti,ti2,ti3,ti4,ti6,ssa,ssb,ssc,ssd,sse
528  real(prec) :: numrep_erfc, ca, force, expti, exptj
529  real(prec) :: tj,tj2,tj3,tj4,tj6,ti2mtj2a,sa,sb,sc,sd,se,sf
530  real(prec) :: ti2mtj2, ti2mti2, tj2mti2
531  integer :: j,k, ccnt, nni
532 
533  coulvol = lbox(1)*lbox(2)*lbox(3)
534  sqrtx = sqrt(-log(coulacc))
535 
536  ccnt = 0
537  coulcut = 12.d0
538  calpha = sqrtx/coulcut
539  coulcut2 = coulcut*coulcut
540  calpha2 = calpha*calpha
541 
542  relperm = one
543  keconst = 14.3996437701414d0*relperm
544  tfact = 16.0d0/(5.0d0*keconst)
545 
546  coulombv = zero
547 
548  ti = tfact*u(i)
549  ti2 = ti*ti
550  ti3 = ti2*ti
551  ti4 = ti2*ti2
552  ti6 = ti4*ti2
553 
554  ssa = ti
555  ssb = ti3/48.d0
556  ssc = 3.d0*ti2/16.d0
557  ssd = 11.d0*ti/16.d0
558  sse = 1.d0
559 
560  ra(1) = rx(i)
561  ra(2) = ry(i)
562  ra(3) = rz(i)
563 
564  ! !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(nnI,J,Rb,Rab,dR,MAGR,MAGR2,TJ,DC,Z,NUMREP_ERFC,CA) &
565  ! !$OMP REDUCTION(+:COULOMBV)
566  do nni = 1,nrnnlist(i)
567  rb(1) = nnrx(i,nni)
568  rb(2) = nnry(i,nni)
569  rb(3) = nnrz(i,nni)
570  j = nntype(i,nni)
571  rab = rb-ra ! OBS b - a !!!
572  dr = norm2(rab)
573  magr = dr
574  magr2 = dr*dr
575 
576  if ((dr <= coulcut).and.(dr > 1e-12)) then
577 
578  tj = tfact*u(j)
579  dc = rab/dr
580 
581  ! Not Using Numerical Recipes ERFC
582  z = abs(calpha*magr)
583  numrep_erfc = erfc(z)
584 
585  ca = numrep_erfc/magr
586  coulombv = coulombv + deltaq(j)*ca
587  ccnt = ccnt + 1
588  !TEST(ccnt) = DELTAQ(J)*CA
589  ca = ca + two*calpha*exp( -calpha2*magr2 )/sqrtpi
590  expti = exp(-ti*magr )
591 
592  if (element_type(i).eq.element_type(j)) then
593  coulombv = coulombv - deltaq(j)*expti*(ssb*magr2 + ssc*magr + ssd + sse/magr)
594  ccnt = ccnt + 1
595  !TEST(ccnt) = - DELTAQ(J)*EXPTI*(SSB*MAGR2 + SSC*MAGR + SSD + SSE/MAGR)
596  else
597  tj2 = tj*tj
598  tj3 = tj2*tj
599  tj4 = tj2*tj2
600  tj6 = tj4*tj2
601  exptj = exp( -tj*magr )
602  ti2mtj2 = ti2 - tj2
603  tj2mti2 = -ti2mtj2
604  sa = ti
605  sb = tj4*ti/(two*ti2mtj2*ti2mtj2)
606  sc = (tj6 - three*tj4*ti2)/(ti2mtj2*ti2mtj2*ti2mtj2)
607  sd = tj
608  se = ti4*tj/(two * tj2mti2 * tj2mti2)
609  sf = (ti6 - three*ti4*tj2)/(tj2mti2*tj2mti2*tj2mti2)
610 
611  coulombv = coulombv - (deltaq(j)*(expti*(sb - (sc/magr)) + exptj*(se - (sf/magr))))
612  endif
613 
614  endif
615  enddo
616  ! !$OMP END PARALLEL DO
617  coulombv = keconst*coulombv
618 
619  end subroutine ewald_real_space_test
620 
621  subroutine ewald_real_space(COULOMBV,FCOUL,I,RX,RY,RZ,LBox, &
622  DELTAQ,U,Element_Type,Nr_atoms,COULACC,TIMERATIO,nnRx,nnRy,nnRz,nrnnlist,nnType,HDIM,Max_Nr_Neigh)
623 
624  implicit none
625 
626  integer, parameter :: prec = 8
627  real(prec), parameter :: one = 1.d0, two = 2.d0, zero = 0.d0, six = 6.d0, three = 3.d0
628  real(prec), parameter :: fourtyeight = 48.d0, eleven = 11.d0, sixteen = 16.d0
629  real(prec), parameter :: pi = 3.14159265358979323846264d0
630  real(prec), parameter :: sqrtpi = 1.772453850905516d0
631  integer, intent(in) :: nr_atoms, hdim, max_nr_neigh, i
632  real(prec), intent(in) :: coulacc, timeratio
633  real(prec) :: tfact, relperm, keconst
634  real(prec), intent(in) :: rx(nr_atoms), ry(nr_atoms), rz(nr_atoms), lbox(3), deltaq(nr_atoms)
635  real(prec), intent(in) :: u(nr_atoms)
636  real(prec) :: coulcut, coulcut2
637  character(10), intent(in) :: element_type(nr_atoms)
638  integer, intent(in) :: nrnnlist(nr_atoms), nntype(nr_atoms,max_nr_neigh)
639  real(prec), intent(in) :: nnrx(nr_atoms,max_nr_neigh)
640  real(prec), intent(in) :: nnry(nr_atoms,max_nr_neigh), nnrz(nr_atoms,max_nr_neigh)
641  real(prec), intent(out) :: coulombv, fcoul(3)
642  real(prec) :: ra(3), rb(3), dr, rab(3)
643  real(prec) :: coulvol, sqrtx, calpha, dc(3), magr, magr2, z
644  real(prec) :: calpha2, ti,ti2,ti3,ti4,ti6,ssa,ssb,ssc,ssd,sse
645  real(prec) :: numrep_erfc, ca, force, expti, exptj
646  real(prec) :: tj,tj2,tj3,tj4,tj6,ti2mtj2a,sa,sb,sc,sd,se,sf
647  real(prec) :: ti2mtj2, ti2mti2, tj2mti2
648  integer :: j,k, ccnt, nni
649 
650  coulvol = lbox(1)*lbox(2)*lbox(3)
651  sqrtx = sqrt(-log(coulacc))
652 
653  ccnt = 0
654  coulcut = 12.d0
655  calpha = sqrtx/coulcut
656  coulcut2 = coulcut*coulcut
657  calpha2 = calpha*calpha
658 
659  relperm = one
660  keconst = 14.3996437701414d0*relperm
661  tfact = 16.0d0/(5.0d0*keconst)
662 
663  fcoul = zero
664  coulombv = zero
665 
666  ti = tfact*u(i)
667  ti2 = ti*ti
668  ti3 = ti2*ti
669  ti4 = ti2*ti2
670  ti6 = ti4*ti2
671 
672  ssa = ti
673  ssb = ti3/48.d0
674  ssc = 3.d0*ti2/16.d0
675  ssd = 11.d0*ti/16.d0
676  sse = 1.d0
677 
678  ra(1) = rx(i)
679  ra(2) = ry(i)
680  ra(3) = rz(i)
681 
682  do nni = 1,nrnnlist(i)
683  rb(1) = nnrx(i,nni)
684  rb(2) = nnry(i,nni)
685  rb(3) = nnrz(i,nni)
686  j = nntype(i,nni)
687  rab = rb-ra ! OBS b - a !!!
688  dr = norm2(rab)
689  magr = dr
690  magr2 = dr*dr
691 
692  if ((dr <= coulcut).and.(dr > 1e-12)) then
693 
694  tj = tfact*u(j)
695  dc = rab/dr
696 
697  ! Not Using Numerical Recipes ERFC
698  z = abs(calpha*magr)
699  numrep_erfc = erfc(z)
700 
701  ca = numrep_erfc/magr
702  coulombv = coulombv + deltaq(j)*ca
703  ccnt = ccnt + 1
704  !TEST(ccnt) = DELTAQ(J)*CA
705  ca = ca + two*calpha*exp( -calpha2*magr2 )/sqrtpi
706  force = -keconst*deltaq(i)*deltaq(j)*ca/magr
707  expti = exp(-ti*magr )
708 
709  if (element_type(i).eq.element_type(j)) then
710  coulombv = coulombv - deltaq(j)*expti*(ssb*magr2 + ssc*magr + ssd + sse/magr)
711  ccnt = ccnt + 1
712  !TEST(ccnt) = - DELTAQ(J)*EXPTI*(SSB*MAGR2 + SSC*MAGR + SSD + SSE/MAGR)
713  force = force + (keconst*deltaq(i)*deltaq(j)*expti)*((sse/magr2 - two*ssb*magr - ssc) &
714  + ssa*(ssb*magr2 + ssc*magr + ssd + sse/magr))
715  else
716  tj2 = tj*tj
717  tj3 = tj2*tj
718  tj4 = tj2*tj2
719  tj6 = tj4*tj2
720  exptj = exp( -tj*magr )
721  ti2mtj2 = ti2 - tj2
722  tj2mti2 = -ti2mtj2
723  sa = ti
724  sb = tj4*ti/(two*ti2mtj2*ti2mtj2)
725  sc = (tj6 - three*tj4*ti2)/(ti2mtj2*ti2mtj2*ti2mtj2)
726  sd = tj
727  se = ti4*tj/(two * tj2mti2 * tj2mti2)
728  sf = (ti6 - three*ti4*tj2)/(tj2mti2*tj2mti2*tj2mti2)
729 
730  coulombv = coulombv - (deltaq(j)*(expti*(sb - (sc/magr)) + exptj*(se - (sf/magr))))
731  force = force + keconst*deltaq(i)*deltaq(j)*((expti*(sa*(sb - (sc/magr)) - (sc/magr2))) &
732  + (exptj*(sd*(se - (sf/magr)) - (sf/magr2))))
733  endif
734 
735  fcoul(1) = fcoul(1) + dc(1)*force
736  fcoul(2) = fcoul(2) + dc(2)*force
737  fcoul(3) = fcoul(3) + dc(3)*force
738  endif
739  enddo
740  coulombv = keconst*coulombv
741 
742  end subroutine ewald_real_space
743 
744  subroutine ewald_k_space_latte(COULOMBV,RXYZ,Box,DELTAQ,Nr_atoms,COULACC,Max_Nr_Neigh)
745 
746  implicit none
747 
748  integer, parameter :: prec = 8
749  real(prec), parameter :: one = 1.d0, two = 2.d0, zero = 0.d0, six = 6.d0, three = 3.d0, four = 4.d0
750  real(prec), parameter :: fourtyeight = 48.d0, eleven = 11.d0, sixteen = 16.d0, eight = 8.d0
751  real(prec), parameter :: pi = 3.14159265358979323846264d0
752  real(prec), parameter :: sqrtpi = 1.772453850905516d0
753  integer, intent(in) :: nr_atoms, max_nr_neigh
754  real(prec), intent(in) :: coulacc
755  real(prec) :: keconst, tfact, relperm
756  real(prec), intent(in) :: rxyz(3,nr_atoms), box(3,3), deltaq(nr_atoms)
757  real(prec) :: coulcut, coulcut2
758  real(prec), intent(out) :: coulombv(nr_atoms)
759  real(prec) :: ra(3), rb(3), dr, rab(3)
760  real(prec) :: coulvol, sqrtx, calpha, dc(3), magr, magr2, z
761  real(prec) :: calpha2, ti,ti2,ti3,ti4,ti6,ssa,ssb,ssc,ssd,sse
762  real(prec) :: corrfact,fourcalpha2, force
763  real(prec) :: recipvecs(3,3),sinlist(nr_atoms),coslist(nr_atoms)
764  real(prec) :: k(3),l11,l12,l13,m21,m22,m23,k2,kcutoff,kcutoff2,prefactor
765  real(prec) :: previr, cossum,sinsum,dot, kepref, cossum2, sinsum2
766 
767  integer :: i,j,l,m,n, ccnt, nni, lmax,mmax,nmax,nmin,mmin
768 
769  coulvol = box(1,1)*box(2,2)*box(3,3)
770  sqrtx = sqrt(-log(coulacc))
771 
772  ccnt = 0
773 
774  coulcut = 12.0d0
775  calpha = sqrtx/coulcut
776 
777  coulcut2 = coulcut*coulcut
778  kcutoff = two*calpha*sqrtx
779  kcutoff2 = kcutoff*kcutoff
780  calpha2 = calpha*calpha
781  fourcalpha2 = four*calpha2
782 
783  recipvecs = zero
784  recipvecs(1,1) = two*pi/box(1,1)
785  recipvecs(2,2) = two*pi/box(2,2)
786  recipvecs(3,3) = two*pi/box(3,3)
787  lmax = floor(kcutoff / sqrt(recipvecs(1,1)*recipvecs(1,1)))
788  mmax = floor(kcutoff / sqrt(recipvecs(2,2)*recipvecs(2,2)))
789  nmax = floor(kcutoff / sqrt(recipvecs(3,3)*recipvecs(3,3)))
790 
791  relperm = 1.d0
792  keconst = 14.3996437701414d0*relperm
793 
794  coulombv = zero
795  sinlist = zero
796  coslist = zero
797 
798  do l = 0,lmax
799 
800  if (l.eq.0) then
801  mmin = 0
802  else
803  mmin = -mmax
804  endif
805 
806  l11 = l*recipvecs(1,1)
807  l12 = l*recipvecs(1,2)
808  l13 = l*recipvecs(1,3)
809 
810  do m = mmin,mmax
811 
812  nmin = -nmax
813  if ((l==0).and.(m==0)) then
814  nmin = 1
815  endif
816 
817  m21 = l11 + m*recipvecs(2,1)
818  m22 = l12 + m*recipvecs(2,2)
819  m23 = l13 + m*recipvecs(2,3)
820 
821  do n = nmin,nmax
822  k(1) = m21 + n*recipvecs(3,1)
823  k(2) = m22 + n*recipvecs(3,2)
824  k(3) = m23 + n*recipvecs(3,3)
825  k2 = k(1)*k(1) + k(2)*k(2) + k(3)*k(3)
826  if (k2.le.kcutoff2) then
827  prefactor = eight*pi*exp(-k2/(4.d0*calpha2))/(coulvol*k2)
828  previr = (2.d0/k2) + (2.d0/(4.d0*calpha2));
829 
830  cossum = 0.d0
831  sinsum = 0.d0
832 
833  ! Doing the sin and cos sums
834  !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(I,DOT) &
835  !$OMP REDUCTION(+:COSSUM) &
836  !$OMP REDUCTION(+:SINSUM)
837  do i = 1,nr_atoms
838  dot = k(1)*rxyz(1,i) + k(2)*rxyz(2,i) + k(3)*rxyz(3,i)
839  ! We re-use these in the next loop...
840  sinlist(i) = sin(dot)
841  coslist(i) = cos(dot)
842  cossum = cossum + deltaq(i)*coslist(i)
843  sinsum = sinsum + deltaq(i)*sinlist(i)
844  enddo
845  !$OMP END PARALLEL DO
846  cossum2 = cossum*cossum
847  sinsum2 = sinsum*sinsum
848 
849  ! Add up energy and force contributions
850 
851  kepref = keconst*prefactor
852  !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(I)
853  do i = 1,nr_atoms
854  coulombv(i) = coulombv(i) + kepref*(coslist(i)*cossum + sinlist(i)*sinsum)
855  enddo
856  !$OMP END PARALLEL DO
857 
858  kepref = kepref*(cossum2 + sinsum2)
859  endif
860  enddo
861  enddo
862  enddo
863 
864  ! Point self energy
865  corrfact = 2.d0*keconst*calpha/sqrtpi;
866  coulombv = coulombv - corrfact*deltaq;
867 
868  end subroutine ewald_k_space_latte
869 
870  subroutine ewald_k_space_matrix_latte(E,RXYZ,Box,Nr_atoms,COULACC,Max_Nr_Neigh,nebcoul,totnebcoul)
871 
872  implicit none
873 
874  integer, parameter :: prec = 8
875  real(prec), parameter :: one = 1.d0, two = 2.d0, zero = 0.d0, six = 6.d0, three = 3.d0, four = 4.d0
876  real(prec), parameter :: fourtyeight = 48.d0, eleven = 11.d0, sixteen = 16.d0, eight = 8.d0
877  real(prec), parameter :: pi = 3.14159265358979323846264d0
878  real(prec), parameter :: sqrtpi = 1.772453850905516d0
879  integer, intent(in) :: nr_atoms, max_nr_neigh,nebcoul(4,max_nr_neigh,nr_atoms),totnebcoul(nr_atoms)
880  real(prec), intent(in) :: coulacc
881  real(prec) :: keconst, tfact, relperm
882  real(prec), intent(in) :: rxyz(3,nr_atoms), box(3,3)
883  real(prec) :: coulcut, coulcut2
884  real(prec), intent(out) :: e(nr_atoms,nr_atoms)
885  real(prec) :: ra(3), rb(3), dr, rab(3)
886  real(prec) :: coulvol, sqrtx, calpha, dc(3), magr, magr2, z
887  real(prec) :: calpha2, ti,ti2,ti3,ti4,ti6,ssa,ssb,ssc,ssd,sse
888  real(prec) :: corrfact,fourcalpha2, force
889  real(prec) :: recipvecs(3,3),sinlist(nr_atoms),coslist(nr_atoms)
890  real(prec) :: k(3),l11,l12,l13,m21,m22,m23,k2,kcutoff,kcutoff2,prefactor
891  real(prec) :: previr, cossum,sinsum,dot, kepref, cossum2, sinsum2
892 
893  integer :: i,j,l,m,n, newj, ccnt, nni, lmax,mmax,nmax,nmin,mmin
894 
895  coulvol = box(1,1)*box(2,2)*box(3,3)
896  sqrtx = sqrt(-log(coulacc))
897 
898  ccnt = 0
899 
900  coulcut = 12.0d0
901  calpha = sqrtx/coulcut
902 
903  coulcut2 = coulcut*coulcut
904  kcutoff = two*calpha*sqrtx
905  kcutoff2 = kcutoff*kcutoff
906  calpha2 = calpha*calpha
907  fourcalpha2 = four*calpha2
908 
909  recipvecs = zero
910  recipvecs(1,1) = two*pi/box(1,1)
911  recipvecs(2,2) = two*pi/box(2,2)
912  recipvecs(3,3) = two*pi/box(3,3)
913  lmax = floor(kcutoff / sqrt(recipvecs(1,1)*recipvecs(1,1)))
914  mmax = floor(kcutoff / sqrt(recipvecs(2,2)*recipvecs(2,2)))
915  nmax = floor(kcutoff / sqrt(recipvecs(3,3)*recipvecs(3,3)))
916 
917  relperm = 1.d0
918  keconst = 14.3996437701414d0*relperm
919 
920  !COULOMBV = ZERO
921  sinlist = zero
922  coslist = zero
923 
924  e = 0.0
925 
926  do l = 0,lmax
927 
928  if (l.eq.0) then
929  mmin = 0
930  else
931  mmin = -mmax
932  endif
933 
934  l11 = l*recipvecs(1,1)
935  l12 = l*recipvecs(1,2)
936  l13 = l*recipvecs(1,3)
937 
938  do m = mmin,mmax
939 
940  nmin = -nmax
941  if ((l==0).and.(m==0)) then
942  nmin = 1
943  endif
944 
945  m21 = l11 + m*recipvecs(2,1)
946  m22 = l12 + m*recipvecs(2,2)
947  m23 = l13 + m*recipvecs(2,3)
948 
949  do n = nmin,nmax
950  k(1) = m21 + n*recipvecs(3,1)
951  k(2) = m22 + n*recipvecs(3,2)
952  k(3) = m23 + n*recipvecs(3,3)
953  k2 = k(1)*k(1) + k(2)*k(2) + k(3)*k(3)
954  if (k2.le.kcutoff2) then
955  prefactor = eight*pi*exp(-k2/(4.d0*calpha2))/(coulvol*k2)
956  previr = (2.d0/k2) + (2.d0/(4.d0*calpha2));
957 
958 
959  ! Doing the sin and cos sums
960  !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(I,DOT)
961  do i = 1,nr_atoms
962  dot = k(1)*rxyz(1,i) + k(2)*rxyz(2,i) + k(3)*rxyz(3,i)
963  ! We re-use these in the next loop...
964  sinlist(i) = sin(dot)
965  coslist(i) = cos(dot)
966  ! COSSUM = COSSUM + DELTAQ(I)*COSLIST(I)
967  ! SINSUM = SINSUM + DELTAQ(I)*SINLIST(I)
968  enddo
969  !$OMP END PARALLEL DO
970 
971  ! Add up energy and force contributions
972 
973  kepref = keconst*prefactor
974  corrfact = 2.d0*keconst*calpha/sqrtpi;
975  !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(I,J,NEWJ)
976  do i = 1,nr_atoms
977  do newj = 1,totnebcoul(i)
978  j = nebcoul(1, newj, i)
979  e(i,j) = e(i,j) + kepref*(coslist(i)*coslist(j)+sinlist(i)*sinlist(j))
980  ! COULOMBV(I) = COULOMBV(I) + KEPREF*(COSLIST(I)*COSSUM + SINLIST(I)*SINSUM)
981  enddo
982  e(i,i) = e(i,i) - corrfact
983  enddo
984  !$OMP END PARALLEL DO
985 
986  !KEPREF = KEPREF*(COSSUM2 + SINSUM2)
987  endif
988  enddo
989  enddo
990  enddo
991 
992  ! Point self energy
993  !CORRFACT = 2.D0*KECONST*CALPHA/SQRTPI;
994  !COULOMBV = COULOMBV - CORRFACT*DELTAQ;
995 
996  end subroutine ewald_k_space_matrix_latte
997 
998  subroutine ewald_k_space_latte_single(COULOMBV,J,RXYZ,Box,DELTAQ,Nr_atoms,COULACC)
999 
1000  implicit none
1001 
1002  integer, parameter :: prec = 8
1003  real(prec), parameter :: one = 1.d0, two = 2.d0, zero = 0.d0, six = 6.d0, three = 3.d0, four = 4.d0
1004  real(prec), parameter :: fourtyeight = 48.d0, eleven = 11.d0, sixteen = 16.d0, eight = 8.d0
1005  real(prec), parameter :: pi = 3.14159265358979323846264d0
1006  real(prec), parameter :: sqrtpi = 1.772453850905516d0
1007  integer, intent(in) :: nr_atoms, j
1008  real(prec), intent(in) :: coulacc
1009  real(prec) :: keconst, tfact, relperm
1010  real(prec), intent(in) :: rxyz(3,nr_atoms), box(3,3), deltaq(nr_atoms)
1011  real(prec) :: coulcut, coulcut2
1012  real(prec), intent(out) :: coulombv(nr_atoms)
1013  real(prec) :: ra(3), rb(3), dr, rab(3)
1014  real(prec) :: coulvol, sqrtx, calpha, dc(3), magr, magr2, z
1015  real(prec) :: calpha2
1016  real(prec) :: corrfact,fourcalpha2
1017  real(prec) :: recipvecs(3,3)
1018  real(prec) :: k(3),l11,l12,l13,m21,m22,m23,k2,kcutoff,kcutoff2,prefactor
1019  real(prec) :: idot, jdot, cosjdot, sinjdot, kepref
1020 
1021  integer :: i,l,m,n, ccnt, nni, lmax,mmax,nmax,nmin,mmin
1022 
1023  coulvol = box(1,1)*box(2,2)*box(3,3)
1024  sqrtx = sqrt(-log(coulacc))
1025 
1026  ccnt = 0
1027 
1028  coulcut = 12.0d0
1029  calpha = sqrtx/coulcut
1030 
1031  coulcut2 = coulcut*coulcut
1032  kcutoff = two*calpha*sqrtx
1033  kcutoff2 = kcutoff*kcutoff
1034  calpha2 = calpha*calpha
1035  fourcalpha2 = four*calpha2
1036 
1037  recipvecs = zero
1038  recipvecs(1,1) = two*pi/box(1,1)
1039  recipvecs(2,2) = two*pi/box(2,2)
1040  recipvecs(3,3) = two*pi/box(3,3)
1041  lmax = floor(kcutoff / sqrt(recipvecs(1,1)*recipvecs(1,1)))
1042  mmax = floor(kcutoff / sqrt(recipvecs(2,2)*recipvecs(2,2)))
1043  nmax = floor(kcutoff / sqrt(recipvecs(3,3)*recipvecs(3,3)))
1044 
1045  relperm = 1.d0
1046  keconst = 14.3996437701414d0*relperm
1047 
1048  coulombv = zero
1049 
1050  !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(I,IDOT,JDOT,COSJDOT,SINJDOT,L,M,N,MMIN,MMAX,NMIN,NMAX,L11,M22,K,K2,PREFACTOR,KEPREF)
1051  do l = 0,lmax
1052 
1053  if (l.eq.0) then
1054  mmin = 0
1055  else
1056  mmin = -mmax
1057  endif
1058 
1059  l11 = l*recipvecs(1,1)
1060 
1061  do m = mmin,mmax
1062 
1063  nmin = -nmax
1064  if ((l==0).and.(m==0)) then
1065  nmin = 1
1066  endif
1067 
1068  m22 = m*recipvecs(2,2)
1069 
1070  do n = nmin,nmax
1071  k(1) = l11
1072  k(2) = m22
1073  k(3) = n*recipvecs(3,3)
1074  k2 = k(1)*k(1) + k(2)*k(2) + k(3)*k(3)
1075 
1076  prefactor = eight*pi*exp(-k2/(4.d0*calpha2))/(coulvol*k2)
1077  kepref = keconst*prefactor
1078  jdot = k(1)*rxyz(1,j) + k(2)*rxyz(2,j) + k(3)*rxyz(3,j)
1079  sinjdot = sin(jdot)
1080  cosjdot = cos(jdot)
1081  do i = 1,nr_atoms
1082  idot = k(1)*rxyz(1,i) + k(2)*rxyz(2,i) + k(3)*rxyz(3,i)
1083  coulombv(i) = coulombv(i) + kepref*deltaq(j)*(cosjdot*cos(idot)+sinjdot*sin(idot))
1084  enddo
1085  enddo
1086  enddo
1087  enddo
1088  !$OMP END PARALLEL DO
1089 
1090  ! Point self energy
1091  corrfact = 2.d0*keconst*calpha/sqrtpi;
1092  coulombv = coulombv - corrfact*deltaq;
1093 
1094  end subroutine ewald_k_space_latte_single
1095 
1096  subroutine ewald_k_space_test(COULOMBV,RX,RY,RZ,LBox,DELTAQ,Nr_atoms,COULACC,Max_Nr_Neigh)
1097  !
1098  implicit none
1099  !
1100  integer, parameter :: prec = 8
1101  real(prec), parameter :: one = 1.d0, two = 2.d0, zero = 0.d0, six = 6.d0, three = 3.d0, four = 4.d0
1102  real(prec), parameter :: fourtyeight = 48.d0, eleven = 11.d0, sixteen = 16.d0, eight = 8.d0
1103  real(prec), parameter :: pi = 3.14159265358979323846264d0
1104  real(prec), parameter :: sqrtpi = 1.772453850905516d0
1105  integer, intent(in) :: nr_atoms, max_nr_neigh
1106  real(prec), intent(in) :: coulacc
1107  real(prec) :: keconst, tfact, relperm
1108  real(prec), intent(in) :: rx(nr_atoms), ry(nr_atoms), rz(nr_atoms), lbox(3), deltaq(nr_atoms)
1109  real(prec) :: coulcut, coulcut2
1110  real(prec), intent(out) :: coulombv(nr_atoms)
1111  real(prec) :: ra(3), rb(3), dr, rab(3)
1112  real(prec) :: coulvol, sqrtx, calpha, dc(3), magr, magr2, z
1113  real(prec) :: calpha2, ti,ti2,ti3,ti4,ti6,ssa,ssb,ssc,ssd,sse
1114  real(prec) :: corrfact,fourcalpha2, force
1115  real(prec) :: recipvecs(3,3),sinlist(nr_atoms),coslist(nr_atoms)
1116  real(prec) :: k(3),l11,l12,l13,m21,m22,m23,k2,kcutoff,kcutoff2,prefactor
1117  real(prec) :: previr, cossum,sinsum,dot, kepref, cossum2, sinsum2
1118 
1119  integer :: i,j,l,m,n, ccnt, nni, lmax,mmax,nmax,nmin,mmin
1120  !
1121  coulvol = lbox(1)*lbox(2)*lbox(3)
1122  sqrtx = sqrt(-log(coulacc))
1123 
1124  ccnt = 0
1125 
1126  coulcut = 12.0d0
1127  calpha = sqrtx/coulcut
1128 
1129  coulcut2 = coulcut*coulcut
1130  kcutoff = two*calpha*sqrtx
1131  kcutoff2 = kcutoff*kcutoff
1132  calpha2 = calpha*calpha
1133  fourcalpha2 = four*calpha2
1134 
1135  recipvecs = zero
1136  recipvecs(1,1) = two*pi/lbox(1)
1137  recipvecs(2,2) = two*pi/lbox(2)
1138  recipvecs(3,3) = two*pi/lbox(3)
1139  lmax = floor(kcutoff / sqrt(recipvecs(1,1)*recipvecs(1,1)))
1140  mmax = floor(kcutoff / sqrt(recipvecs(2,2)*recipvecs(2,2)))
1141  nmax = floor(kcutoff / sqrt(recipvecs(3,3)*recipvecs(3,3)))
1142 
1143  relperm = 1.d0
1144  keconst = 14.3996437701414d0*relperm
1145 
1146  coulombv = zero
1147  sinlist = zero
1148  coslist = zero
1149 
1150  do l = 0,lmax
1151 
1152  if (l.eq.0) then
1153  mmin = 0
1154  else
1155  mmin = -mmax
1156  endif
1157 
1158  l11 = l*recipvecs(1,1)
1159  l12 = l*recipvecs(1,2)
1160  l13 = l*recipvecs(1,3)
1161 
1162  do m = mmin,mmax
1163 
1164  nmin = -nmax
1165  if ((l==0).and.(m==0)) then
1166  nmin = 1
1167  endif
1168 
1169  m21 = l11 + m*recipvecs(2,1)
1170  m22 = l12 + m*recipvecs(2,2)
1171  m23 = l13 + m*recipvecs(2,3)
1172 
1173  do n = nmin,nmax
1174  k(1) = m21 + n*recipvecs(3,1)
1175  k(2) = m22 + n*recipvecs(3,2)
1176  k(3) = m23 + n*recipvecs(3,3)
1177  k2 = k(1)*k(1) + k(2)*k(2) + k(3)*k(3)
1178  if (k2.le.kcutoff2) then
1179  prefactor = eight*pi*exp(-k2/(4.d0*calpha2))/(coulvol*k2)
1180  previr = (2.d0/k2) + (2.d0/(4.d0*calpha2));
1181 
1182  cossum = 0.d0
1183  sinsum = 0.d0
1184 
1185  ! Doing the sin and cos sums
1186  !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(I,DOT) &
1187  !$OMP REDUCTION(+:COSSUM) &
1188  !$OMP REDUCTION(+:SINSUM)
1189  do i = 1,nr_atoms
1190  dot = k(1)*rx(i) + k(2)*ry(i) + k(3)*rz(i)
1191  ! We re-use these in the next loop...
1192  sinlist(i) = sin(dot)
1193  coslist(i) = cos(dot)
1194  cossum = cossum + deltaq(i)*coslist(i)
1195  sinsum = sinsum + deltaq(i)*sinlist(i)
1196  enddo
1197  !$OMP END PARALLEL DO
1198  cossum2 = cossum*cossum
1199  sinsum2 = sinsum*sinsum
1200 
1201  ! Add up energy and force contributions
1202  !
1203  kepref = keconst*prefactor
1204  !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(I)
1205  do i = 1,nr_atoms
1206  coulombv(i) = coulombv(i) + kepref*(coslist(i)*cossum + sinlist(i)*sinsum)
1207  enddo
1208  !$OMP END PARALLEL DO
1209 
1210  kepref = kepref*(cossum2 + sinsum2)
1211  endif
1212  enddo
1213  enddo
1214  enddo
1215 
1216  ! Point self energy
1217  corrfact = 2.d0*keconst*calpha/sqrtpi;
1218  coulombv = coulombv - corrfact*deltaq;
1219 
1220  end subroutine ewald_k_space_test
1221 
1222 end module prg_ewald_mod
prg_ewald_mod
Definition: prg_ewald_mod.F90:2
prg_ewald_mod::dp
integer, parameter dp
Definition: prg_ewald_mod.F90:12
prg_ewald_mod::ewald_k_space_latte
subroutine, public ewald_k_space_latte(COULOMBV, RXYZ, Box, DELTAQ, Nr_atoms, COULACC, Max_Nr_Neigh)
Definition: prg_ewald_mod.F90:745
prg_ewald_mod::ewald_k_space_latte_single
subroutine, public ewald_k_space_latte_single(COULOMBV, J, RXYZ, Box, DELTAQ, Nr_atoms, COULACC)
Definition: prg_ewald_mod.F90:999
prg_ewald_mod::ewald_real_space_single_latte
subroutine, public ewald_real_space_single_latte(COULOMBV, I, RXYZ, Box, Nr_elem, DELTAQ, J, U, Element_Pointer, Nr_atoms, COULACC, HDIM, Max_Nr_Neigh)
Find Coulomb potential on site I from single charge at site J.
Definition: prg_ewald_mod.F90:30
prg_ewald_mod::ewald_k_space_test
subroutine, public ewald_k_space_test(COULOMBV, RX, RY, RZ, LBox, DELTAQ, Nr_atoms, COULACC, Max_Nr_Neigh)
Definition: prg_ewald_mod.F90:1097
prg_ewald_mod::ewald_real_space
subroutine, public ewald_real_space(COULOMBV, FCOUL, I, RX, RY, RZ, LBox, DELTAQ, U, Element_Type, Nr_atoms, COULACC, TIMERATIO, nnRx, nnRy, nnRz, nrnnlist, nnType, HDIM, Max_Nr_Neigh)
Definition: prg_ewald_mod.F90:623
prg_ewald_mod::ewald_real_space_test
subroutine, public ewald_real_space_test(COULOMBV, I, RX, RY, RZ, LBox, DELTAQ, U, Element_Type, Nr_atoms, COULACC, nnRx, nnRy, nnRz, nrnnlist, nnType, Max_Nr_Neigh)
Definition: prg_ewald_mod.F90:506
prg_ewald_mod::ewald_real_space_latte
subroutine, public ewald_real_space_latte(COULOMBV, I, RXYZ, Box, DELTAQ, U, Element_Pointer, Nr_atoms, COULACC, nebcoul, totnebcoul, HDIM, Max_Nr_Neigh, Nr_Elem)
Definition: prg_ewald_mod.F90:389
prg_ewald_mod::ewald_k_space_matrix_latte
subroutine, public ewald_k_space_matrix_latte(E, RXYZ, Box, Nr_atoms, COULACC, Max_Nr_Neigh, nebcoul, totnebcoul)
Definition: prg_ewald_mod.F90:871
prg_ewald_mod::ewald_real_space_matrix_latte
subroutine, public ewald_real_space_matrix_latte(E, RXYZ, Box, U, Element_Pointer, Nr_atoms, COULACC, nebcoul, totnebcoul, HDIM, Max_Nr_Neigh, Nr_Elem)
Definition: prg_ewald_mod.F90:265
prg_ewald_mod::ewald_real_space_single
subroutine, public ewald_real_space_single(COULOMBV, FCOUL, I, RX, RY, RZ, LBox, DELTAQ, J, U, Element_Type, Nr_atoms, COULACC, TIMERATIO, HDIM, Max_Nr_Neigh)
Definition: prg_ewald_mod.F90:142
prg_parallel_mod
The parallel module.
Definition: prg_parallel_mod.F90:6
prg_timer_mod
The timer module.
Definition: prg_timer_mod.F90:23