1*
2* $Id$
3*
4
5      subroutine lattice_min_difference(x,y,z)
6      implicit none
7      real*8 x,y,z
8
9*     **** common block ****
10      real*8 ecut,wcut,omega
11      real*8 ua(3,3),unitg(3,3)
12      common / lattice_block / ua,unitg,ecut,wcut,omega
13
14      real*8 ub(3,3)
15      common / lattice_block2 / ub
16
17*     *** local variables ****
18      real*8 c1,c2,c3
19
20      c1 = x*ub(1,1) + y*ub(2,1) + z*ub(3,1)
21      c2 = x*ub(1,2) + y*ub(2,2) + z*ub(3,2)
22      c3 = x*ub(1,3) + y*ub(2,3) + z*ub(3,3)
23      c1 = c1 - DNINT(c1)
24      c2 = c2 - DNINT(c2)
25      c3 = c3 - DNINT(c3)
26      x = ua(1,1)*c1 + ua(1,2)*c2 + ua(1,3)*c3
27      y = ua(2,1)*c1 + ua(2,2)*c2 + ua(2,3)*c3
28      z = ua(3,1)*c1 + ua(3,2)*c2 + ua(3,3)*c3
29
30      return
31      end
32
33      subroutine lattice_frac_to_r1(n,f1,r1)
34      implicit none
35      integer n
36      real*8 f1(3,*),r1(3,*)
37
38*     **** common block ****
39      real*8 ecut,wcut,omega
40      real*8 ua(3,3),unitg(3,3)
41      common / lattice_block / ua,unitg,ecut,wcut,omega
42
43*     **** local variables ***
44      integer i
45
46      do i=1,n
47         r1(1,i) = f1(1,i)*ua(1,1) + f1(2,i)*ua(2,1) + f1(3,i)*ua(3,1)
48         r1(2,i) = f1(1,i)*ua(1,2) + f1(2,i)*ua(2,2) + f1(3,i)*ua(3,2)
49         r1(3,i) = f1(1,i)*ua(1,3) + f1(2,i)*ua(2,3) + f1(3,i)*ua(3,3)
50      end do
51      return
52      end
53
54      subroutine lattice_r1_to_frac(n,r1,f1)
55      implicit none
56      integer n
57      real*8 r1(3,*),f1(3,*)
58
59*     **** common block ****
60      real*8 ub(3,3)
61      common / lattice_block2 / ub
62
63*     **** local variables ***
64      integer i
65
66      do i=1,n
67         f1(1,i) = r1(1,i)*ub(1,1) + r1(2,i)*ub(2,1) + r1(3,i)*ub(3,1)
68         f1(2,i) = r1(1,i)*ub(1,2) + r1(2,i)*ub(2,2) + r1(3,i)*ub(3,2)
69         f1(3,i) = r1(1,i)*ub(1,3) + r1(2,i)*ub(2,3) + r1(3,i)*ub(3,3)
70      end do
71      return
72      end
73
74
75      subroutine lattice_center_xyz_to_ijk(x,y,z,c1,c2,c3)
76      implicit none
77      real*8 x,y,z
78      integer c1,c2,c3
79
80      integer np1,np2,np3
81      real*8  f1,f2,f3
82
83*     **** common block ****
84      real*8 ub(3,3)
85      common / lattice_block2 / ub
86
87      call D3dB_nx(1,np1)
88      call D3dB_ny(1,np2)
89      call D3dB_nz(1,np3)
90      f1 = x*ub(1,1) + y*ub(2,1) + z*ub(3,1)
91      f2 = x*ub(1,2) + y*ub(2,2) + z*ub(3,2)
92      f3 = x*ub(1,3) + y*ub(2,3) + z*ub(3,3)
93      c1 = (f1-0.5d0+0.5d0/dble(np1))*np1
94      c2 = (f2-0.5d0+0.5d0/dble(np2))*np2
95      c3 = (f3-0.5d0+0.5d0/dble(np3))*np3
96      return
97      end
98
99      subroutine lattice_center0_xyz_to_ijk(x,y,z,c1,c2,c3)
100      implicit none
101      real*8 x,y,z
102      integer c1,c2,c3
103
104      integer np1,np2,np3
105      real*8  f1,f2,f3
106
107*     **** common block ****
108      real*8 ub(3,3)
109      common / lattice_block2 / ub
110
111      call D3dB_nx(1,np1)
112      call D3dB_ny(1,np2)
113      call D3dB_nz(1,np3)
114      f1 = x*ub(1,1) + y*ub(2,1) + z*ub(3,1)
115      f2 = x*ub(1,2) + y*ub(2,2) + z*ub(3,2)
116      f3 = x*ub(1,3) + y*ub(2,3) + z*ub(3,3)
117      c1 = (f1+0.5d0/dble(np1))*np1
118      c2 = (f2+0.5d0/dble(np2))*np2
119      c3 = (f3+0.5d0/dble(np3))*np3
120      return
121      end
122
123
124      subroutine lattice_fragcell(n1,r1)
125      implicit none
126      real*8  rcm(3)
127      integer n1
128      real*8  r1(3,*)
129
130*     **** common block ****
131      real*8 ecut,wcut,omega
132      real*8 ua(3,3),unitg(3,3)
133      common / lattice_block / ua,unitg,ecut,wcut,omega
134
135*     **** common block ****
136      real*8 ub(3,3)
137      common / lattice_block2 / ub
138
139*     **** local variables ****
140      integer i,j
141      real*8 fcm
142
143      do i=2,n1
144         do j=1,3
145            fcm = (r1(1,i)-r1(1,1))*ub(1,j)
146     >          + (r1(2,i)-r1(2,1))*ub(2,j)
147     >          + (r1(3,i)-r1(3,1))*ub(3,j)
148            do while (fcm.gt.0.5)
149               r1(1,i) = r1(1,i) - ua(1,j)
150               r1(2,i) = r1(2,i) - ua(2,j)
151               r1(3,i) = r1(3,i) - ua(3,j)
152               fcm = (r1(1,i)-r1(1,1))*ub(1,j)
153     >             + (r1(2,i)-r1(2,1))*ub(2,j)
154     >             + (r1(3,i)-r1(3,1))*ub(3,j)
155            end do
156            do while (fcm.le.(-0.5))
157               r1(1,i) = r1(1,i) + ua(1,j)
158               r1(2,i) = r1(2,i) + ua(2,j)
159               r1(3,i) = r1(3,i) + ua(3,j)
160               fcm = (r1(1,i)-r1(1,1))*ub(1,j)
161     >             + (r1(2,i)-r1(2,1))*ub(2,j)
162     >             + (r1(3,i)-r1(3,1))*ub(3,j)
163            end do
164         end do
165      end do
166      return
167      end
168
169
170      subroutine lattice_incell1_frag(rcm,n1,r1)
171      implicit none
172      real*8  rcm(3)
173      integer n1
174      real*8  r1(3,*)
175
176*     **** common block ****
177      real*8 ecut,wcut,omega
178      real*8 ua(3,3),unitg(3,3)
179      common / lattice_block / ua,unitg,ecut,wcut,omega
180
181*     **** common block ****
182      real*8 ub(3,3)
183      common / lattice_block2 / ub
184
185*     **** local variables ****
186      integer i
187      real*8 fcm(3)
188
189      call lattice_r1_to_frac(1,rcm,fcm)
190
191      do while (fcm(1).gt.(0.5d0))
192         rcm(1) = rcm(1) - ua(1,1)
193         rcm(2) = rcm(2) - ua(2,1)
194         rcm(3) = rcm(3) - ua(3,1)
195         do i=1,n1
196           r1(1,i) = r1(1,i) - ua(1,1)
197           r1(2,i) = r1(2,i) - ua(2,1)
198           r1(3,i) = r1(3,i) - ua(3,1)
199         end do
200         fcm(1) = rcm(1)*ub(1,1) + rcm(2)*ub(2,1) + rcm(3)*ub(3,1)
201      end do
202      do while (fcm(1).le.(-0.5d0))
203         rcm(1) = rcm(1) + ua(1,1)
204         rcm(2) = rcm(2) + ua(2,1)
205         rcm(3) = rcm(3) + ua(3,1)
206         do i=1,n1
207           r1(1,i) = r1(1,i) + ua(1,1)
208           r1(2,i) = r1(2,i) + ua(2,1)
209           r1(3,i) = r1(3,i) + ua(3,1)
210         end do
211         fcm(1) = rcm(1)*ub(1,1) + rcm(2)*ub(2,1) + rcm(3)*ub(3,1)
212      end do
213
214      do while (fcm(2).gt.(0.5d0))
215         rcm(1) = rcm(1) - ua(1,2)
216         rcm(2) = rcm(2) - ua(2,2)
217         rcm(3) = rcm(3) - ua(3,2)
218         do i=1,n1
219           r1(1,i) = r1(1,i) - ua(1,2)
220           r1(2,i) = r1(2,i) - ua(2,2)
221           r1(3,i) = r1(3,i) - ua(3,2)
222         end do
223         fcm(2) = rcm(1)*ub(1,2) + rcm(2)*ub(2,2) + rcm(3)*ub(3,2)
224      end do
225      do while (fcm(2).le.(-0.5d0))
226         rcm(1) = rcm(1) + ua(1,2)
227         rcm(2) = rcm(2) + ua(2,2)
228         rcm(3) = rcm(3) + ua(3,2)
229         do i=1,n1
230           r1(1,i) = r1(1,i) + ua(1,2)
231           r1(2,i) = r1(2,i) + ua(2,2)
232           r1(3,i) = r1(3,i) + ua(3,2)
233         end do
234         fcm(2) = rcm(1)*ub(1,2) + rcm(2)*ub(2,2) + rcm(3)*ub(3,2)
235      end do
236
237      do while (fcm(3).gt.(0.5d0))
238         rcm(1) = rcm(1) - ua(1,3)
239         rcm(2) = rcm(2) - ua(2,3)
240         rcm(3) = rcm(3) - ua(3,3)
241         do i=1,n1
242           r1(1,i) = r1(1,i) - ua(1,3)
243           r1(2,i) = r1(2,i) - ua(2,3)
244           r1(3,i) = r1(3,i) - ua(3,3)
245         end do
246         fcm(3) = rcm(1)*ub(1,3) + rcm(2)*ub(2,3) + rcm(3)*ub(3,3)
247      end do
248      do while (fcm(3).le.(-0.5d0))
249         rcm(1) = rcm(1) + ua(1,3)
250         rcm(2) = rcm(2) + ua(2,3)
251         rcm(3) = rcm(3) + ua(3,3)
252         do i=1,n1
253           r1(1,i) = r1(1,i) + ua(1,3)
254           r1(2,i) = r1(2,i) + ua(2,3)
255           r1(3,i) = r1(3,i) + ua(3,3)
256         end do
257         fcm(3) = rcm(1)*ub(1,3) + rcm(2)*ub(2,3) + rcm(3)*ub(3,3)
258      end do
259      return
260      end
261
262      subroutine lattice_incell2_frag(rcm,n1,r1,r2)
263      implicit none
264      real*8  rcm(3)
265      integer n1
266      real*8  r1(3,*)
267      real*8  r2(3,*)
268
269*     **** common block ****
270      real*8 ecut,wcut,omega
271      real*8 ua(3,3),unitg(3,3)
272      common / lattice_block / ua,unitg,ecut,wcut,omega
273
274*     **** common block ****
275      real*8 ub(3,3)
276      common / lattice_block2 / ub
277
278*     **** local variables ****
279      integer i
280      real*8 fcm(3)
281
282      call lattice_r1_to_frac(1,rcm,fcm)
283
284      do while (fcm(1).gt.(0.5d0))
285         rcm(1) = rcm(1) - ua(1,1)
286         rcm(2) = rcm(2) - ua(2,1)
287         rcm(3) = rcm(3) - ua(3,1)
288         do i=1,n1
289           r1(1,i) = r1(1,i) - ua(1,1)
290           r1(2,i) = r1(2,i) - ua(2,1)
291           r1(3,i) = r1(3,i) - ua(3,1)
292
293           r2(1,i) = r2(1,i) - ua(1,1)
294           r2(2,i) = r2(2,i) - ua(2,1)
295           r2(3,i) = r2(3,i) - ua(3,1)
296         end do
297         fcm(1) = rcm(1)*ub(1,1) + rcm(2)*ub(2,1) + rcm(3)*ub(3,1)
298      end do
299      do while (fcm(1).le.(-0.5d0))
300         rcm(1) = rcm(1) + ua(1,1)
301         rcm(2) = rcm(2) + ua(2,1)
302         rcm(3) = rcm(3) + ua(3,1)
303         do i=1,n1
304           r1(1,i) = r1(1,i) + ua(1,1)
305           r1(2,i) = r1(2,i) + ua(2,1)
306           r1(3,i) = r1(3,i) + ua(3,1)
307
308           r2(1,i) = r2(1,i) + ua(1,1)
309           r2(2,i) = r2(2,i) + ua(2,1)
310           r2(3,i) = r2(3,i) + ua(3,1)
311         end do
312         fcm(1) = rcm(1)*ub(1,1) + rcm(2)*ub(2,1) + rcm(3)*ub(3,1)
313      end do
314
315      do while (fcm(2).gt.(0.5d0))
316         rcm(1) = rcm(1) - ua(1,2)
317         rcm(2) = rcm(2) - ua(2,2)
318         rcm(3) = rcm(3) - ua(3,2)
319         do i=1,n1
320           r1(1,i) = r1(1,i) - ua(1,2)
321           r1(2,i) = r1(2,i) - ua(2,2)
322           r1(3,i) = r1(3,i) - ua(3,2)
323
324           r2(1,i) = r2(1,i) - ua(1,2)
325           r2(2,i) = r2(2,i) - ua(2,2)
326           r2(3,i) = r2(3,i) - ua(3,2)
327         end do
328         fcm(2) = rcm(1)*ub(1,2) + rcm(2)*ub(2,2) + rcm(3)*ub(3,2)
329      end do
330      do while (fcm(2).le.(-0.5d0))
331         rcm(1) = rcm(1) + ua(1,2)
332         rcm(2) = rcm(2) + ua(2,2)
333         rcm(3) = rcm(3) + ua(3,2)
334         do i=1,n1
335           r1(1,i) = r1(1,i) + ua(1,2)
336           r1(2,i) = r1(2,i) + ua(2,2)
337           r1(3,i) = r1(3,i) + ua(3,2)
338
339           r2(1,i) = r2(1,i) + ua(1,2)
340           r2(2,i) = r2(2,i) + ua(2,2)
341           r2(3,i) = r2(3,i) + ua(3,2)
342         end do
343         fcm(2) = rcm(1)*ub(1,2) + rcm(2)*ub(2,2) + rcm(3)*ub(3,2)
344      end do
345
346      do while (fcm(3).gt.(0.5d0))
347         rcm(1) = rcm(1) - ua(1,3)
348         rcm(2) = rcm(2) - ua(2,3)
349         rcm(3) = rcm(3) - ua(3,3)
350         do i=1,n1
351           r1(1,i) = r1(1,i) - ua(1,3)
352           r1(2,i) = r1(2,i) - ua(2,3)
353           r1(3,i) = r1(3,i) - ua(3,3)
354
355           r2(1,i) = r2(1,i) - ua(1,3)
356           r2(2,i) = r2(2,i) - ua(2,3)
357           r2(3,i) = r2(3,i) - ua(3,3)
358         end do
359         fcm(3) = rcm(1)*ub(1,3) + rcm(2)*ub(2,3) + rcm(3)*ub(3,3)
360      end do
361      do while (fcm(3).le.(-0.5d0))
362         rcm(1) = rcm(1) + ua(1,3)
363         rcm(2) = rcm(2) + ua(2,3)
364         rcm(3) = rcm(3) + ua(3,3)
365         do i=1,n1
366           r1(1,i) = r1(1,i) + ua(1,3)
367           r1(2,i) = r1(2,i) + ua(2,3)
368           r1(3,i) = r1(3,i) + ua(3,3)
369
370           r2(1,i) = r2(1,i) + ua(1,3)
371           r2(2,i) = r2(2,i) + ua(2,3)
372           r2(3,i) = r2(3,i) + ua(3,3)
373         end do
374         fcm(3) = rcm(1)*ub(1,3) + rcm(2)*ub(2,3) + rcm(3)*ub(3,3)
375      end do
376      return
377      end
378
379      subroutine lattice_incell3_frag(rcm,n1,r1,r2,r3)
380      implicit none
381      real*8  rcm(3)
382      integer n1
383      real*8  r1(3,*)
384      real*8  r2(3,*)
385      real*8  r3(3,*)
386
387*     **** common block ****
388      real*8 ecut,wcut,omega
389      real*8 ua(3,3),unitg(3,3)
390      common / lattice_block / ua,unitg,ecut,wcut,omega
391
392*     **** common block ****
393      real*8 ub(3,3)
394      common / lattice_block2 / ub
395
396*     **** local variables ****
397      integer i
398      real*8 fcm(3)
399
400      call lattice_r1_to_frac(1,rcm,fcm)
401
402      do while (fcm(1).gt.(0.5d0))
403         rcm(1) = rcm(1) - ua(1,1)
404         rcm(2) = rcm(2) - ua(2,1)
405         rcm(3) = rcm(3) - ua(3,1)
406         do i=1,n1
407           r1(1,i) = r1(1,i) - ua(1,1)
408           r1(2,i) = r1(2,i) - ua(2,1)
409           r1(3,i) = r1(3,i) - ua(3,1)
410
411           r2(1,i) = r2(1,i) - ua(1,1)
412           r2(2,i) = r2(2,i) - ua(2,1)
413           r2(3,i) = r2(3,i) - ua(3,1)
414
415           r3(1,i) = r3(1,i) - ua(1,1)
416           r3(2,i) = r3(2,i) - ua(2,1)
417           r3(3,i) = r3(3,i) - ua(3,1)
418         end do
419         fcm(1) = rcm(1)*ub(1,1) + rcm(2)*ub(2,1) + rcm(3)*ub(3,1)
420      end do
421      do while (fcm(1).le.(-0.5d0))
422         rcm(1) = rcm(1) + ua(1,1)
423         rcm(2) = rcm(2) + ua(2,1)
424         rcm(3) = rcm(3) + ua(3,1)
425         do i=1,n1
426           r1(1,i) = r1(1,i) + ua(1,1)
427           r1(2,i) = r1(2,i) + ua(2,1)
428           r1(3,i) = r1(3,i) + ua(3,1)
429
430           r2(1,i) = r2(1,i) + ua(1,1)
431           r2(2,i) = r2(2,i) + ua(2,1)
432           r2(3,i) = r2(3,i) + ua(3,1)
433
434           r3(1,i) = r3(1,i) + ua(1,1)
435           r3(2,i) = r3(2,i) + ua(2,1)
436           r3(3,i) = r3(3,i) + ua(3,1)
437         end do
438         fcm(1) = rcm(1)*ub(1,1) + rcm(2)*ub(2,1) + rcm(3)*ub(3,1)
439      end do
440
441      do while (fcm(2).gt.(0.5d0))
442         rcm(1) = rcm(1) - ua(1,2)
443         rcm(2) = rcm(2) - ua(2,2)
444         rcm(3) = rcm(3) - ua(3,2)
445         do i=1,n1
446           r1(1,i) = r1(1,i) - ua(1,2)
447           r1(2,i) = r1(2,i) - ua(2,2)
448           r1(3,i) = r1(3,i) - ua(3,2)
449
450           r2(1,i) = r2(1,i) - ua(1,2)
451           r2(2,i) = r2(2,i) - ua(2,2)
452           r2(3,i) = r2(3,i) - ua(3,2)
453
454           r3(1,i) = r3(1,i) - ua(1,2)
455           r3(2,i) = r3(2,i) - ua(2,2)
456           r3(3,i) = r3(3,i) - ua(3,2)
457         end do
458         fcm(2) = rcm(1)*ub(1,2) + rcm(2)*ub(2,2) + rcm(3)*ub(3,2)
459      end do
460      do while (fcm(2).le.(-0.5d0))
461         rcm(1) = rcm(1) + ua(1,2)
462         rcm(2) = rcm(2) + ua(2,2)
463         rcm(3) = rcm(3) + ua(3,2)
464         do i=1,n1
465           r1(1,i) = r1(1,i) + ua(1,2)
466           r1(2,i) = r1(2,i) + ua(2,2)
467           r1(3,i) = r1(3,i) + ua(3,2)
468
469           r2(1,i) = r2(1,i) + ua(1,2)
470           r2(2,i) = r2(2,i) + ua(2,2)
471           r2(3,i) = r2(3,i) + ua(3,2)
472
473           r3(1,i) = r3(1,i) + ua(1,2)
474           r3(2,i) = r3(2,i) + ua(2,2)
475           r3(3,i) = r3(3,i) + ua(3,2)
476         end do
477         fcm(2) = rcm(1)*ub(1,2) + rcm(2)*ub(2,2) + rcm(3)*ub(3,2)
478      end do
479
480      do while (fcm(3).gt.(0.5d0))
481         rcm(1) = rcm(1) - ua(1,3)
482         rcm(2) = rcm(2) - ua(2,3)
483         rcm(3) = rcm(3) - ua(3,3)
484         do i=1,n1
485           r1(1,i) = r1(1,i) - ua(1,3)
486           r1(2,i) = r1(2,i) - ua(2,3)
487           r1(3,i) = r1(3,i) - ua(3,3)
488
489           r2(1,i) = r2(1,i) - ua(1,3)
490           r2(2,i) = r2(2,i) - ua(2,3)
491           r2(3,i) = r2(3,i) - ua(3,3)
492
493           r3(1,i) = r3(1,i) - ua(1,3)
494           r3(2,i) = r3(2,i) - ua(2,3)
495           r3(3,i) = r3(3,i) - ua(3,3)
496         end do
497         fcm(3) = rcm(1)*ub(1,3) + rcm(2)*ub(2,3) + rcm(3)*ub(3,3)
498      end do
499      do while (fcm(3).le.(-0.5d0))
500         rcm(1) = rcm(1) + ua(1,3)
501         rcm(2) = rcm(2) + ua(2,3)
502         rcm(3) = rcm(3) + ua(3,3)
503         do i=1,n1
504           r1(1,i) = r1(1,i) + ua(1,3)
505           r1(2,i) = r1(2,i) + ua(2,3)
506           r1(3,i) = r1(3,i) + ua(3,3)
507
508           r2(1,i) = r2(1,i) + ua(1,3)
509           r2(2,i) = r2(2,i) + ua(2,3)
510           r2(3,i) = r2(3,i) + ua(3,3)
511
512           r3(1,i) = r3(1,i) + ua(1,3)
513           r3(2,i) = r3(2,i) + ua(2,3)
514           r3(3,i) = r3(3,i) + ua(3,3)
515         end do
516         fcm(3) = rcm(1)*ub(1,3) + rcm(2)*ub(2,3) + rcm(3)*ub(3,3)
517      end do
518      return
519      end
520
521
522      real*8 function lattice_wcut()
523      implicit none
524
525*     **** common block ****
526      real*8 ecut,wcut,omega
527      real*8 unita(3,3),unitg(3,3)
528      common / lattice_block / unita,unitg,ecut,wcut,omega
529
530      lattice_wcut = wcut
531      return
532      end
533
534      real*8 function lattice_ecut()
535      implicit none
536
537*     **** common block ****
538      real*8 ecut,wcut,omega
539      real*8 unita(3,3),unitg(3,3)
540      common / lattice_block / unita,unitg,ecut,wcut,omega
541
542      lattice_ecut = ecut
543      return
544      end
545
546      real*8 function lattice_ggcut()
547      implicit none
548
549*     **** common block ****
550      real*8 ecut,wcut,omega
551      real*8 unita(3,3),unitg(3,3)
552      common / lattice_block / unita,unitg,ecut,wcut,omega
553
554      lattice_ggcut = 2.0d0*ecut
555      return
556      end
557
558      real*8 function lattice_wggcut()
559      implicit none
560
561*     **** common block ****
562      real*8 ecut,wcut,omega
563      real*8 unita(3,3),unitg(3,3)
564      common / lattice_block / unita,unitg,ecut,wcut,omega
565
566      lattice_wggcut = 2.0d0*wcut
567      return
568      end
569
570
571
572      real*8 function lattice_omega()
573      implicit none
574
575*     **** common block ****
576      real*8 ecut,wcut,omega
577      real*8 unita(3,3),unitg(3,3)
578      common / lattice_block / unita,unitg,ecut,wcut,omega
579
580      lattice_omega = omega
581      return
582      end
583
584      real*8 function lattice_unita(i,j)
585      implicit none
586      integer i,j
587
588*     **** common block ****
589      real*8 ecut,wcut,omega
590      real*8 unita(3,3),unitg(3,3)
591      common / lattice_block / unita,unitg,ecut,wcut,omega
592
593      lattice_unita = unita(i,j)
594      return
595      end
596
597
598      real*8 function lattice_unitg(i,j)
599      implicit none
600      integer i,j
601
602*     **** common block ****
603      real*8 ecut,wcut,omega
604      real*8 unita(3,3),unitg(3,3)
605      common / lattice_block / unita,unitg,ecut,wcut,omega
606
607      lattice_unitg = unitg(i,j)
608      return
609      end
610
611*     ********************************************
612*     *                                          *
613*     *          lattice_unita_small             *
614*     *                                          *
615*     ********************************************
616      real*8 function lattice_unita_small(i,j)
617      implicit none
618      integer i,j
619
620*     **** common block ****
621      logical has_small
622      real*8  omega_small
623      real*8  unita_small(3,3),unitg_small(3,3)
624      common /lattice_small_block/ unita_small,unitg_small,
625     >                             omega_small,has_small
626
627      lattice_unita_small = unita_small(i,j)
628      return
629      end
630
631*     ********************************************
632*     *                                          *
633*     *          lattice_unita_frozen_small      *
634*     *                                          *
635*     ********************************************
636      real*8 function lattice_unita_frozen_small(i,j)
637      implicit none
638      integer i,j
639
640*     **** common block ****
641      real*8  omega_frozen_small
642      real*8  unita_frozen_small(3,3),unitg_frozen_small(3,3)
643      common /lattice_small_frozen_block/ unita_frozen_small,
644     >                             unitg_frozen_small,
645     >                             omega_frozen_small
646
647      lattice_unita_frozen_small = unita_frozen_small(i,j)
648      return
649      end
650
651
652*     ********************************************
653*     *                                          *
654*     *          lattice_unitg_small             *
655*     *                                          *
656*     ********************************************
657      real*8 function lattice_unitg_small(i,j)
658      implicit none
659      integer i,j
660
661*     **** common block ****
662      logical has_small
663      real*8  omega_small
664      real*8  unita_small(3,3),unitg_small(3,3)
665      common /lattice_small_block/ unita_small,unitg_small,
666     >                             omega_small,has_small
667
668      lattice_unitg_small = unitg_small(i,j)
669      return
670      end
671
672*     ********************************************
673*     *                                          *
674*     *          lattice_unitg_frozen_small      *
675*     *                                          *
676*     ********************************************
677      real*8 function lattice_unitg_frozen_small(i,j)
678      implicit none
679      integer i,j
680
681*     **** common block ****
682      real*8  omega_frozen_small
683      real*8  unita_frozen_small(3,3),unitg_frozen_small(3,3)
684      common /lattice_small_frozen_block/ unita_frozen_small,
685     >                             unitg_frozen_small,
686     >                             omega_frozen_small
687
688      lattice_unitg_frozen_small = unitg_frozen_small(i,j)
689      return
690      end
691
692
693
694*     ********************************************
695*     *                                          *
696*     *          lattice_omega_small             *
697*     *                                          *
698*     ********************************************
699      real*8 function lattice_omega_small(i,j)
700      implicit none
701      integer i,j
702
703*     **** common block ****
704      logical has_small
705      real*8  omega_small
706      real*8  unita_small(3,3),unitg_small(3,3)
707      common /lattice_small_block/ unita_small,unitg_small,
708     >                             omega_small,has_small
709
710      lattice_omega_small = omega_small
711      return
712      end
713
714*     ********************************************
715*     *                                          *
716*     *          lattice_omega_frozen_small      *
717*     *                                          *
718*     ********************************************
719      real*8 function lattice_omega_frozen_small()
720      implicit none
721      integer i,j
722
723*     **** common block ****
724      real*8  omega_frozen_small
725      real*8  unita_frozen_small(3,3),unitg_frozen_small(3,3)
726      common /lattice_small_frozen_block/ unita_frozen_small,
727     >                             unitg_frozen_small,
728     >                             omega_frozen_small
729
730      lattice_omega_frozen_small = omega_frozen_small
731      return
732      end
733
734*     ********************************************
735*     *                                          *
736*     *          lattice_has_small               *
737*     *                                          *
738*     ********************************************
739      logical function lattice_has_small()
740      implicit none
741
742*     **** common block ****
743      logical has_small
744      real*8  omega_small
745      real*8  unita_small(3,3),unitg_small(3,3)
746      common /lattice_small_block/ unita_small,unitg_small,
747     >                             omega_small,has_small
748
749      lattice_has_small = has_small
750      return
751      end
752
753
754
755
756
757*     ********************************************
758*     *                                          *
759*     *          lattice_unitg_frozen            *
760*     *                                          *
761*     ********************************************
762
763*     frozen lattice structure - used to determine masking arrays
764
765      real*8 function lattice_unitg_frozen(i,j)
766      implicit none
767      integer i,j
768
769*     **** common blocks ****
770      logical frozen
771      real*8 ecut_frozen,wcut_frozen,omega_frozen
772      real*8 unita_frozen(3,3),unitg_frozen(3,3)
773      common / lattice_froze / unita_frozen,unitg_frozen,
774     >                         ecut_frozen,wcut_frozen,
775     >                         omega_frozen,frozen
776
777      lattice_unitg_frozen = unitg_frozen(i,j)
778      return
779      end
780
781*     ********************************************
782*     *                                          *
783*     *          lattice_unita_frozen            *
784*     *                                          *
785*     ********************************************
786
787*     frozen lattice structure - used to determine masking arrays
788
789      real*8 function lattice_unita_frozen(i,j)
790      implicit none
791      integer i,j
792
793*     **** common blocks ****
794      logical frozen
795      real*8 ecut_frozen,wcut_frozen,omega_frozen
796      real*8 unita_frozen(3,3),unitg_frozen(3,3)
797      common / lattice_froze / unita_frozen,unitg_frozen,
798     >                         ecut_frozen,wcut_frozen,
799     >                         omega_frozen,frozen
800
801      lattice_unita_frozen = unita_frozen(i,j)
802      return
803      end
804
805*     ********************************************
806*     *                                          *
807*     *          lattice_omega_frozen            *
808*     *                                          *
809*     ********************************************
810
811*     frozen lattice structure - used to determine masking arrays
812
813      real*8 function lattice_omega_frozen()
814      implicit none
815
816*     **** common blocks ****
817      logical frozen
818      real*8 ecut_frozen,wcut_frozen,omega_frozen
819      real*8 unita_frozen(3,3),unitg_frozen(3,3)
820      common / lattice_froze / unita_frozen,unitg_frozen,
821     >                         ecut_frozen,wcut_frozen,
822     >                         omega_frozen,frozen
823
824      lattice_omega_frozen = omega_frozen
825      return
826      end
827
828*     ********************************************
829*     *                                          *
830*     *          lattice_wcut_frozen             *
831*     *                                          *
832*     ********************************************
833
834*     frozen lattice structure - used to determine masking arrays
835
836      real*8 function lattice_wcut_frozen()
837      implicit none
838
839*     **** common blocks ****
840      logical frozen
841      real*8 ecut_frozen,wcut_frozen,omega_frozen
842      real*8 unita_frozen(3,3),unitg_frozen(3,3)
843      common / lattice_froze / unita_frozen,unitg_frozen,
844     >                         ecut_frozen,wcut_frozen,
845     >                         omega_frozen,frozen
846
847      lattice_wcut_frozen = wcut_frozen
848      return
849      end
850
851*     ********************************************
852*     *                                          *
853*     *          lattice_ecut_frozen             *
854*     *                                          *
855*     ********************************************
856
857*     frozen lattice structure - used to determine masking arrays
858
859      real*8 function lattice_ecut_frozen()
860      implicit none
861
862*     **** common blocks ****
863      logical frozen
864      real*8 ecut_frozen,wcut_frozen,omega_frozen
865      real*8 unita_frozen(3,3),unitg_frozen(3,3)
866      common / lattice_froze / unita_frozen,unitg_frozen,
867     >                         ecut_frozen,wcut_frozen,
868     >                         omega_frozen,frozen
869
870      lattice_ecut_frozen = ecut_frozen
871      return
872      end
873
874*     ********************************************
875*     *                                          *
876*     *          lattice_ggcut_frozen            *
877*     *                                          *
878*     ********************************************
879
880*     frozen lattice structure - used to determine masking arrays
881
882      real*8 function lattice_ggcut_frozen()
883      implicit none
884
885*     **** common blocks ****
886      logical frozen
887      real*8 ecut_frozen,wcut_frozen,omega_frozen
888      real*8 unita_frozen(3,3),unitg_frozen(3,3)
889      common / lattice_froze / unita_frozen,unitg_frozen,
890     >                         ecut_frozen,wcut_frozen,
891     >                         omega_frozen,frozen
892
893      lattice_ggcut_frozen = 2.0d0*ecut_frozen
894      return
895      end
896
897*     ********************************************
898*     *                                          *
899*     *         lattice_wggcut_frozen            *
900*     *                                          *
901*     ********************************************
902
903*     frozen lattice structure - used to determine masking arrays
904
905      real*8 function lattice_wggcut_frozen()
906      implicit none
907
908*     **** common blocks ****
909      logical frozen
910      real*8 ecut_frozen,wcut_frozen,omega_frozen
911      real*8 unita_frozen(3,3),unitg_frozen(3,3)
912      common / lattice_froze / unita_frozen,unitg_frozen,
913     >                         ecut_frozen,wcut_frozen,
914     >                         omega_frozen,frozen
915
916      lattice_wggcut_frozen = 2.0d0*wcut_frozen
917      return
918      end
919
920*     ********************************************
921*     *                                          *
922*     *         lattice_frozen                   *
923*     *                                          *
924*     ********************************************
925
926*     frozen lattice structure - used to determine masking arrays
927
928      logical function lattice_frozen()
929      implicit none
930
931*     **** common blocks ****
932      logical frozen
933      real*8 ecut_frozen,wcut_frozen,omega_frozen
934      real*8 unita_frozen(3,3),unitg_frozen(3,3)
935      common / lattice_froze / unita_frozen,unitg_frozen,
936     >                         ecut_frozen,wcut_frozen,
937     >                         omega_frozen,frozen
938
939      lattice_frozen = frozen
940      return
941      end
942
943
944
945
946*     *********************************************
947*     *                                           *
948*     *               lattice_init                *
949*     *                                           *
950*     *********************************************
951
952      subroutine lattice_init()
953      implicit none
954
955*     **** common block ****
956      real*8 ecut,wcut,omega
957      real*8 unita(3,3),unitg(3,3)
958      common / lattice_block / unita,unitg,ecut,wcut,omega
959
960      real*8 ub(3,3)
961      common / lattice_block2 / ub
962
963      logical frozen
964      real*8 ecut_frozen,wcut_frozen,omega_frozen
965      real*8 unita_frozen(3,3),unitg_frozen(3,3)
966      common / lattice_froze / unita_frozen,unitg_frozen,
967     >                         ecut_frozen,wcut_frozen,
968     >                         omega_frozen,frozen
969
970      logical has_small
971      real*8  omega_small
972      real*8  unita_small(3,3),unitg_small(3,3)
973      common /lattice_small_block/ unita_small,unitg_small,
974     >                             omega_small,has_small
975
976      real*8  omega_frozen_small
977      real*8  unita_frozen_small(3,3),unitg_frozen_small(3,3)
978      common /lattice_small_frozen_block/ unita_frozen_small,
979     >                             unitg_frozen_small,
980     >                             omega_frozen_small
981
982
983*     **** local variables ****
984      integer nx,ny,nz
985      integer nxh,nyh,nzh
986      real*8  gx,gy,gz,gg
987      real*8  gg1,gg2,gg3
988      real*8  ecut0,wcut0
989
990*     **** external functions ****
991      logical  control_frozen,control_has_ngrid_small
992      integer  control_ngrid,control_ngrid_small
993      real*8   control_unita,control_ecut,control_wcut
994      real*8   control_unita_frozen
995      external control_frozen,control_has_ngrid_small
996      external control_ngrid,control_ngrid_small
997      external control_unita,control_ecut,control_wcut
998      external control_unita_frozen
999
1000      ecut0 = control_ecut()
1001      wcut0 = control_wcut()
1002
1003*     **** define lattice ****
1004      unita(1,1) = control_unita(1,1)
1005      unita(2,1) = control_unita(2,1)
1006      unita(3,1) = control_unita(3,1)
1007      unita(1,2) = control_unita(1,2)
1008      unita(2,2) = control_unita(2,2)
1009      unita(3,2) = control_unita(3,2)
1010      unita(1,3) = control_unita(1,3)
1011      unita(2,3) = control_unita(2,3)
1012      unita(3,3) = control_unita(3,3)
1013      call get_cube(unita,unitg,omega)
1014
1015*     **** define frozen lattice - note this is only different from unita if nwpw:frozen_lattice is set ****
1016      frozen = control_frozen()
1017      unita_frozen(1,1) = control_unita_frozen(1,1)
1018      unita_frozen(2,1) = control_unita_frozen(2,1)
1019      unita_frozen(3,1) = control_unita_frozen(3,1)
1020      unita_frozen(1,2) = control_unita_frozen(1,2)
1021      unita_frozen(2,2) = control_unita_frozen(2,2)
1022      unita_frozen(3,2) = control_unita_frozen(3,2)
1023      unita_frozen(1,3) = control_unita_frozen(1,3)
1024      unita_frozen(2,3) = control_unita_frozen(2,3)
1025      unita_frozen(3,3) = control_unita_frozen(3,3)
1026      call get_cube(unita_frozen,unitg_frozen,omega_frozen)
1027
1028*     **** define ub ****
1029      call dcopy(9,unitg,1,ub,1)
1030      call dscal(9,(1.0d0/(8.0d0*datan(1.0d0))),ub,1)
1031
1032
1033
1034*     *** set the ecut variable using the frozen lattice***
1035c     call D3dB_nx(1,nx)
1036c     call D3dB_ny(1,ny)
1037c     call D3dB_nz(1,nz)
1038      nx = control_ngrid(1)
1039      ny = control_ngrid(2)
1040      nz = control_ngrid(3)
1041      nxh = nx/2
1042      nyh = ny/2
1043      nzh = nz/2
1044
1045      gx = unitg_frozen(1,1)*dble(nxh)
1046      gy = unitg_frozen(2,1)*dble(nxh)
1047      gz = unitg_frozen(3,1)*dble(nxh)
1048      gg1 = gx*gx + gy*gy + gz*gz
1049
1050      gx = unitg_frozen(1,2)*dble(nyh)
1051      gy = unitg_frozen(2,2)*dble(nyh)
1052      gz = unitg_frozen(3,2)*dble(nyh)
1053      gg2 = gx*gx + gy*gy + gz*gz
1054
1055      gx = unitg_frozen(1,3)*dble(nzh)
1056      gy = unitg_frozen(2,3)*dble(nzh)
1057      gz = unitg_frozen(3,3)*dble(nzh)
1058      gg3 = gx*gx + gy*gy + gz*gz
1059
1060      gg = gg1
1061      if (gg2.lt.gg) gg=gg2
1062      if (gg3.lt.gg) gg=gg3
1063
1064      ecut = 0.5d0*gg
1065      if (ecut0.lt.ecut) then
1066         ecut = ecut0
1067      end if
1068
1069      wcut = ecut
1070      if (wcut0.lt.wcut) then
1071         wcut = wcut0
1072      end if
1073
1074*     **** set ecut,wcut for frozen lattice ****
1075      wcut_frozen = wcut
1076      ecut_frozen = ecut
1077
1078
1079*     **** small cell ****
1080      has_small = control_has_ngrid_small()
1081      if (has_small) then
1082         gg1 = dble(control_ngrid_small(1))/dble(nx)
1083         gg2 = dble(control_ngrid_small(2))/dble(ny)
1084         gg3 = dble(control_ngrid_small(3))/dble(nz)
1085         unita_small(1,1) = unita(1,1)*gg1
1086         unita_small(2,1) = unita(2,1)*gg1
1087         unita_small(3,1) = unita(3,1)*gg1
1088         unita_small(1,2) = unita(1,2)*gg2
1089         unita_small(2,2) = unita(2,2)*gg2
1090         unita_small(3,2) = unita(3,2)*gg2
1091         unita_small(1,3) = unita(1,3)*gg3
1092         unita_small(2,3) = unita(2,3)*gg3
1093         unita_small(3,3) = unita(3,3)*gg3
1094         call get_cube(unita_small,unitg_small,omega_small)
1095         unita_frozen_small(1,1) = unita_frozen(1,1)*gg1
1096         unita_frozen_small(2,1) = unita_frozen(2,1)*gg1
1097         unita_frozen_small(3,1) = unita_frozen(3,1)*gg1
1098         unita_frozen_small(1,2) = unita_frozen(1,2)*gg2
1099         unita_frozen_small(2,2) = unita_frozen(2,2)*gg2
1100         unita_frozen_small(3,2) = unita_frozen(3,2)*gg2
1101         unita_frozen_small(1,3) = unita_frozen(1,3)*gg3
1102         unita_frozen_small(2,3) = unita_frozen(2,3)*gg3
1103         unita_frozen_small(3,3) = unita_frozen(3,3)*gg3
1104         call get_cube(unita_frozen_small,
1105     >                 unitg_frozen_small,
1106     >                 omega_frozen_small)
1107      end if
1108
1109      return
1110      end
1111
1112
1113*     *******************************
1114*     *                             *
1115*     *         lattice_p_grid      *
1116*     *                             *
1117*     *******************************
1118*
1119*     This routine computes coordinates of grid points in
1120*     the unit cell
1121*
1122*     Uses -
1123*          Parallel_taskid --- processor number
1124*          D3dB_nx --- number of grid points in direction 1
1125*          D3dB_ny --- number of grid points in direction 2
1126*          D3dB_nz --- number of grid points in direction 2
1127*          lattice_unita -- primitive lattice vectors in real space
1128*
1129*     Exit -
1130*          xs  --- coordinates of grid points (Rx,Ry,Rz)
1131*
1132*
1133      subroutine lattice_p_grid(xs)
1134      implicit none
1135      real*8 xs(*)
1136
1137*     **** local variables ****
1138      integer nfft3d,n2ft3d
1139      integer i,j,k,p,taskid
1140      integer idx,k1,k2,k3
1141      integer np1,np2,np3
1142      integer nph1,nph2,nph3
1143      real*8  a(3,3),dk1,dk2,dk3,twopi
1144
1145
1146*     **** constants ****
1147      call Parallel2d_taskid_i(taskid)
1148      call D3dB_nfft3d(1,nfft3d)
1149      n2ft3d = 2*nfft3d
1150      call D3dB_nx(1,np1)
1151      call D3dB_ny(1,np2)
1152      call D3dB_nz(1,np3)
1153      twopi = 8.0d0*datan(1.0d0)
1154
1155      nph1 = np1/2
1156      nph2 = np2/2
1157      nph3 = np3/2
1158
1159*     **** elemental vectors ****
1160      call dcopy(9,0.0d0,0,a,1)
1161      a(1,1) = twopi/dble(np1)
1162      a(2,2) = twopi/dble(np2)
1163      a(3,3) = twopi/dble(np3)
1164
1165      call dcopy(6*n2ft3d,0.0d0,0,xs,1)
1166
1167*     **** grid points in coordination space ****
1168      do k3 = -nph3, nph3-1
1169        do k2 = -nph2, nph2-1
1170          do k1 = -nph1, nph1-1
1171
1172               i = k1 + nph1
1173               j = k2 + nph2
1174               k = k3 + nph3
1175
1176               call D3dB_ijktoindex2p(1,i+1,j+1,k+1,idx,p)
1177               if (p .eq. taskid) then
1178                dk1=dble(k1)
1179                dk2=dble(k2)
1180                dk3=dble(k3)
1181                xs(idx)         =dcos(a(1,1)*dk1+a(1,2)*dk2+a(1,3)*dk3)
1182                xs(idx+n2ft3d)  =dsin(a(1,1)*dk1+a(1,2)*dk2+a(1,3)*dk3)
1183                xs(idx+2*n2ft3d)=dcos(a(2,1)*dk1+a(2,2)*dk2+a(2,3)*dk3)
1184                xs(idx+3*n2ft3d)=dsin(a(2,1)*dk1+a(2,2)*dk2+a(2,3)*dk3)
1185                xs(idx+4*n2ft3d)=dcos(a(3,1)*dk1+a(3,2)*dk2+a(3,3)*dk3)
1186                xs(idx+5*n2ft3d)=dsin(a(3,1)*dk1+a(3,2)*dk2+a(3,3)*dk3)
1187               end if
1188          end do
1189        end do
1190      end do
1191
1192      return
1193      end
1194
1195
1196
1197
1198*     *******************************
1199*     *                             *
1200*     *         lattice_r_grid      *
1201*     *                             *
1202*     *******************************
1203*
1204*     This routine computes coordinates of grid points in
1205*     the unit cell
1206*
1207*     Uses -
1208*          Parallel_taskid --- processor number
1209*          D3dB_nx --- number of grid points in direction 1
1210*          D3dB_ny --- number of grid points in direction 2
1211*          D3dB_nz --- number of grid points in direction 2
1212*          lattice_unita -- primitive lattice vectors in real space
1213*
1214*     Exit -
1215*          r  --- coordinates of grid points (Rx,Ry,Rz)
1216*
1217*
1218      subroutine lattice_r_grid(r)
1219      implicit none
1220      real*8 r(3,*)
1221
1222*     **** local variables ****
1223      integer nfft3d,n2ft3d
1224      integer i,j,k,p,taskid,tid,nthreads
1225      integer index,k1,k2,k3,it
1226      integer np1,np2,np3
1227      integer nph1,nph2,nph3
1228      real*8  a(3,3),dk1,dk2,dk3
1229
1230*     **** external functions ****
1231      real*8   lattice_unita
1232      external lattice_unita
1233      integer  Parallel_threadid,Parallel_nthreads
1234      external Parallel_threadid,Parallel_nthreads
1235
1236
1237*     **** constants ****
1238      call Parallel2d_taskid_i(taskid)
1239      tid      = Parallel_threadid()
1240      nthreads = Parallel_nthreads()
1241      call D3dB_nfft3d(1,nfft3d)
1242      n2ft3d = 2*nfft3d
1243      call D3dB_nx(1,np1)
1244      call D3dB_ny(1,np2)
1245      call D3dB_nz(1,np3)
1246
1247      nph1 = np1/2
1248      nph2 = np2/2
1249      nph3 = np3/2
1250
1251*     **** elemental vectors ****
1252      do i=1,3
1253         a(i,1) = lattice_unita(i,1)/dble(np1)
1254         a(i,2) = lattice_unita(i,2)/dble(np2)
1255         a(i,3) = lattice_unita(i,3)/dble(np3)
1256      end do
1257
1258      !call dcopy(3*n2ft3d,0.0d0,0,r,1)
1259      call Parallel_shared_vector_zero(.true.,2*n2ft3d,r)
1260
1261*     **** grid points in coordination space ****
1262      it = 0
1263      do k3 = -nph3, nph3-1
1264        do k2 = -nph2, nph2-1
1265          do k1 = -nph1, nph1-1
1266
1267               i = k1 + nph1
1268               j = k2 + nph2
1269               k = k3 + nph3
1270
1271               !call D3dB_ktoqp(1,k+1,q,p)
1272               call D3dB_ijktoindex2p(1,i+1,j+1,k+1,index,p)
1273               if ((p.eq.taskid).and.(it.eq.tid)) then
1274c                 index = (q-1)*(np1+2)*np2
1275c    >                  + j    *(np1+2)
1276c    >                  + i+1
1277                dk1=dble(k1)
1278                dk2=dble(k2)
1279                dk3=dble(k3)
1280                r(1,index) = a(1,1)*dk1 + a(1,2)*dk2 + a(1,3)*dk3
1281                r(2,index) = a(2,1)*dk1 + a(2,2)*dk2 + a(2,3)*dk3
1282                r(3,index) = a(3,1)*dk1 + a(3,2)*dk2 + a(3,3)*dk3
1283
1284c*               **** reverse y and z ****
1285c                r(1,index) = a(1,1)*k1 + a(1,2)*k3 + a(1,3)*k2
1286c                r(2,index) = a(2,1)*k1 + a(2,2)*k3 + a(2,3)*k2
1287c                r(3,index) = a(3,1)*k1 + a(3,2)*k3 + a(3,3)*k2
1288
1289               end if
1290               it = mod(it+1,nthreads)
1291          end do
1292        end do
1293      end do
1294
1295      return
1296      end
1297
1298
1299      subroutine lattice_r_grid_sym(r)
1300      implicit none
1301      real*8 r(3,*)
1302
1303*     **** local variables ****
1304      integer nfft3d,n2ft3d
1305      integer i,j,k,p,taskid
1306      integer index,k1,k2,k3
1307      integer np1,np2,np3
1308      integer nph1,nph2,nph3
1309      real*8  a(3,3),dk1,dk2,dk3
1310
1311*     **** external functions ****
1312      real*8   lattice_unita
1313      external lattice_unita
1314
1315
1316*     **** constants ****
1317      call Parallel2d_taskid_i(taskid)
1318      call D3dB_nfft3d(1,nfft3d)
1319      n2ft3d = 2*nfft3d
1320      call D3dB_nx(1,np1)
1321      call D3dB_ny(1,np2)
1322      call D3dB_nz(1,np3)
1323
1324      nph1 = np1/2
1325      nph2 = np2/2
1326      nph3 = np3/2
1327
1328*     **** elemental vectors ****
1329      do i=1,3
1330         a(i,1) = lattice_unita(i,1)/dble(np1)
1331         a(i,2) = lattice_unita(i,2)/dble(np2)
1332         a(i,3) = lattice_unita(i,3)/dble(np3)
1333      end do
1334
1335      call dcopy(3*n2ft3d,0.0d0,0,r,1)
1336
1337*     **** grid points in coordination space ****
1338      do k3 = -nph3+1, nph3-1
1339        do k2 = -nph2+1, nph2-1
1340          do k1 = -nph1+1, nph1-1
1341
1342               i = k1 + nph1
1343               j = k2 + nph2
1344               k = k3 + nph3
1345
1346               !call D3dB_ktoqp(1,k+1,q,p)
1347               call D3dB_ijktoindex2p(1,i+1,j+1,k+1,index,p)
1348               if (p .eq. taskid) then
1349c                 index = (q-1)*(np1+2)*np2
1350c    >                  + j    *(np1+2)
1351c    >                  + i+1
1352                dk1=dble(k1)
1353                dk2=dble(k2)
1354                dk3=dble(k3)
1355                r(1,index) = a(1,1)*dk1 + a(1,2)*dk2 + a(1,3)*dk3
1356                r(2,index) = a(2,1)*dk1 + a(2,2)*dk2 + a(2,3)*dk3
1357                r(3,index) = a(3,1)*dk1 + a(3,2)*dk2 + a(3,3)*dk3
1358
1359c*               **** reverse y and z ****
1360c                r(1,index) = a(1,1)*k1 + a(1,2)*k3 + a(1,3)*k2
1361c                r(2,index) = a(2,1)*k1 + a(2,2)*k3 + a(2,3)*k2
1362c                r(3,index) = a(3,1)*k1 + a(3,2)*k3 + a(3,3)*k2
1363
1364               end if
1365          end do
1366        end do
1367      end do
1368
1369      return
1370      end
1371
1372
1373      subroutine lattice_r_grid_sym0(r)
1374      implicit none
1375      real*8 r(*)
1376
1377*     **** local variables ****
1378      integer nfft3d,n2ft3d
1379      integer i,j,k,p,taskid
1380      integer index,k1,k2,k3
1381      integer np1,np2,np3
1382      integer nph1,nph2,nph3
1383      real*8  a(3,3),dk1,dk2,dk3
1384
1385*     **** external functions ****
1386      real*8   lattice_unita
1387      external lattice_unita
1388
1389
1390*     **** constants ****
1391      call Parallel2d_taskid_i(taskid)
1392      call D3dB_nfft3d(1,nfft3d)
1393      n2ft3d = 2*nfft3d
1394      call D3dB_nx(1,np1)
1395      call D3dB_ny(1,np2)
1396      call D3dB_nz(1,np3)
1397
1398      nph1 = np1/2
1399      nph2 = np2/2
1400      nph3 = np3/2
1401
1402*     **** elemental vectors ****
1403      do i=1,3
1404         a(i,1) = lattice_unita(i,1)/dble(np1)
1405         a(i,2) = lattice_unita(i,2)/dble(np2)
1406         a(i,3) = lattice_unita(i,3)/dble(np3)
1407      end do
1408
1409      call dcopy(3*n2ft3d,0.0d0,0,r,1)
1410
1411*     **** grid points in coordination space ****
1412      do k3 = -nph3+1, nph3-1
1413        do k2 = -nph2+1, nph2-1
1414          do k1 = -nph1+1, nph1-1
1415
1416               i = k1 + nph1
1417               j = k2 + nph2
1418               k = k3 + nph3
1419
1420               !call D3dB_ktoqp(1,k+1,q,p)
1421               call D3dB_ijktoindex2p(1,i+1,j+1,k+1,index,p)
1422               if (p .eq. taskid) then
1423c                 index = (q-1)*(np1+2)*np2
1424c    >                  + j    *(np1+2)
1425c    >                  + i+1
1426                dk1=dble(k1)
1427                dk2=dble(k2)
1428                dk3=dble(k3)
1429                r(index)          = a(1,1)*dk1 + a(1,2)*dk2 + a(1,3)*dk3
1430                r(index+  n2ft3d) = a(2,1)*dk1 + a(2,2)*dk2 + a(2,3)*dk3
1431                r(index+2*n2ft3d) = a(3,1)*dk1 + a(3,2)*dk2 + a(3,3)*dk3
1432               end if
1433          end do
1434        end do
1435      end do
1436
1437      return
1438      end
1439
1440
1441*     *******************************
1442*     *                             *
1443*     *         lattice_i_grid      *
1444*     *                             *
1445*     *******************************
1446*
1447*     This routine computes coordinates of grid points in
1448*     the unit cell
1449*
1450*     Uses -
1451*          Parallel_taskid --- processor number
1452*          D3dB_nx --- number of grid points in direction 1
1453*          D3dB_ny --- number of grid points in direction 2
1454*          D3dB_nz --- number of grid points in direction 2
1455*
1456*     Entry - nb
1457*     Exit -
1458*          ijk --- coordinates of grid points (k1,k2,k3)
1459*
1460*
1461      subroutine lattice_i_grid(nb,ijk)
1462      implicit none
1463      integer nb
1464      integer ijk(4,*)
1465
1466*     **** local variables ****
1467      integer nfft3d,n2ft3d
1468      integer i,j,k,p,taskid
1469      integer index,k1,k2,k3
1470      integer np1,np2,np3
1471      integer nph1,nph2,nph3
1472
1473*     **** constants ****
1474      call Parallel2d_taskid_i(taskid)
1475      call D3dB_nfft3d(nb,nfft3d)
1476
1477      n2ft3d = 2*nfft3d
1478      call D3dB_nx(nb,np1)
1479      call D3dB_ny(nb,np2)
1480      call D3dB_nz(nb,np3)
1481
1482      nph1 = np1/2
1483      nph2 = np2/2
1484      nph3 = np3/2
1485
1486      call icopy(4*n2ft3d,0,0,ijk,1)
1487
1488*     **** grid points  ****
1489      do k3 = -nph3, nph3-1
1490        do k2 = -nph2, nph2-1
1491          do k1 = -nph1, nph1-1
1492
1493               i = k1 + nph1
1494               j = k2 + nph2
1495               k = k3 + nph3
1496
1497               call D3dB_ijktoindex2p(nb,i+1,j+1,k+1,index,p)
1498               if (p .eq. taskid) then
1499                  ijk(1,index) = k1
1500                  ijk(2,index) = k2
1501                  ijk(3,index) = k3
1502                  ijk(4,index) = 1
1503               end if
1504          end do
1505        end do
1506      end do
1507
1508      return
1509      end
1510
1511
1512
1513
1514
1515
1516
1517      subroutine lattice_mask_sym(r)
1518      implicit none
1519      real*8 r(*)
1520
1521*     **** local variables ****
1522      integer nfft3d,n2ft3d
1523      integer i,j,k,p,taskid
1524      integer index,k1,k2,k3
1525      integer np1,np2,np3
1526      integer nph1,nph2,nph3
1527
1528
1529*     **** constants ****
1530      call Parallel2d_taskid_i(taskid)
1531      call D3dB_nfft3d(1,nfft3d)
1532      n2ft3d = 2*nfft3d
1533      call D3dB_nx(1,np1)
1534      call D3dB_ny(1,np2)
1535      call D3dB_nz(1,np3)
1536
1537      nph1 = np1/2
1538      nph2 = np2/2
1539      nph3 = np3/2
1540
1541
1542      call dcopy(n2ft3d,0.0d0,0,r,1)
1543
1544*     **** grid points in coordination space ****
1545      do k3 = -nph3+1, nph3-1
1546        do k2 = -nph2+1, nph2-1
1547          do k1 = -nph1+1, nph1-1
1548
1549               i = k1 + nph1
1550               j = k2 + nph2
1551               k = k3 + nph3
1552
1553               !call D3dB_ktoqp(1,k+1,q,p)
1554               call D3dB_ijktoindex2p(1,i+1,j+1,k+1,index,p)
1555               if (p .eq. taskid) then
1556c                 index = (q-1)*(np1+2)*np2
1557c    >                  + j    *(np1+2)
1558c    >                  + i+1
1559
1560                r(index) =  1.0d0
1561               end if
1562          end do
1563        end do
1564      end do
1565
1566      return
1567      end
1568
1569
1570
1571
1572      subroutine get_cube(unita,unitg,volume)
1573
1574******************************************************************************
1575*                                                                            *
1576*     This routine computes primitive vectors both in coordination           *
1577*     space and in reciporocal space and the volume of primitive cell.       *
1578*                                                                            *
1579*     Inputs:                                                                *
1580*             type --- type of cube (1=SC, 2=FCC, 3=BCC, 4=linear)           *
1581*             unit --- lattice constants                                     *
1582*                                                                            *
1583*     Outputs:                                                               *
1584*             volume --- volume of primitive cell                            *
1585*             unita  --- primitive vectors in coordination space             *
1586*             unitg  --- primitive vectors in reciprocal space               *
1587*                                                                            *
1588*     Library:  DSCAL from BLAS                                              *
1589*                                                                            *
1590*     Last modification:  7/03/93  by R. Kawai                               *
1591*                                                                            *
1592******************************************************************************
1593
1594      implicit none
1595
1596*     ------------------
1597*     argument variables
1598*     ------------------
1599      double precision unita(3,3), unitg(3,3)
1600      double precision volume
1601
1602*     ---------------
1603*     local variables
1604*     ---------------
1605      double precision twopi
1606
1607      twopi = 8.0d0*datan(1.0d0)
1608
1609
1610*     -----------------------------------------
1611*     primitive vectors in the reciprocal space
1612*     -----------------------------------------
1613      unitg(1,1) = unita(2,2)*unita(3,3) - unita(3,2)*unita(2,3)
1614      unitg(2,1) = unita(3,2)*unita(1,3) - unita(1,2)*unita(3,3)
1615      unitg(3,1) = unita(1,2)*unita(2,3) - unita(2,2)*unita(1,3)
1616      unitg(1,2) = unita(2,3)*unita(3,1) - unita(3,3)*unita(2,1)
1617      unitg(2,2) = unita(3,3)*unita(1,1) - unita(1,3)*unita(3,1)
1618      unitg(3,2) = unita(1,3)*unita(2,1) - unita(2,3)*unita(1,1)
1619      unitg(1,3) = unita(2,1)*unita(3,2) - unita(3,1)*unita(2,2)
1620      unitg(2,3) = unita(3,1)*unita(1,2) - unita(1,1)*unita(3,2)
1621      unitg(3,3) = unita(1,1)*unita(2,2) - unita(2,1)*unita(1,2)
1622*     ---------------------
1623*     volume of a unit cell
1624*     ---------------------
1625      volume = unita(1,1)*unitg(1,1)
1626     >       + unita(2,1)*unitg(2,1)
1627     >       + unita(3,1)*unitg(3,1)
1628      call dscal(9,twopi/volume,unitg,1)
1629
1630      volume=dabs(volume)
1631      return
1632      end
1633
1634*     *******************************
1635*     *                             *
1636*     *     lattice_abc_abg         *
1637*     *                             *
1638*     *******************************
1639*
1640*     This routine computes a,b,c,alpha,beta,gamma.
1641*
1642      subroutine lattice_abc_abg(a,b,c,alpha,beta,gamma)
1643      implicit none
1644      real*8 a,b,c
1645      real*8 alpha,beta,gamma
1646
1647*     *** local variables ****
1648      real*8 d2,pi
1649
1650*     **** external functions ****
1651      real*8   lattice_unita
1652      external lattice_unita
1653
1654*     **** determine a,b,c,alpha,beta,gmma ***
1655      pi = 4.0d0*datan(1.0d0)
1656      a = dsqrt(lattice_unita(1,1)**2
1657     >        + lattice_unita(2,1)**2
1658     >        + lattice_unita(3,1)**2)
1659      b = dsqrt(lattice_unita(1,2)**2
1660     >        + lattice_unita(2,2)**2
1661     >        + lattice_unita(3,2)**2)
1662      c = dsqrt(lattice_unita(1,3)**2
1663     >        + lattice_unita(2,3)**2
1664     >        + lattice_unita(3,3)**2)
1665
1666      d2 = (lattice_unita(1,2)-lattice_unita(1,3))**2
1667     >   + (lattice_unita(2,2)-lattice_unita(2,3))**2
1668     >   + (lattice_unita(3,2)-lattice_unita(3,3))**2
1669      alpha = (b*b + c*c - d2)/(2.0d0*b*c)
1670      alpha = dacos(alpha)*180.0d0/pi
1671
1672      d2 = (lattice_unita(1,3)-lattice_unita(1,1))**2
1673     >   + (lattice_unita(2,3)-lattice_unita(2,1))**2
1674     >   + (lattice_unita(3,3)-lattice_unita(3,1))**2
1675      beta = (c*c + a*a - d2)/(2.0d0*c*a)
1676      beta = dacos(beta)*180.0d0/pi
1677
1678      d2 = (lattice_unita(1,1)-lattice_unita(1,2))**2
1679     >   + (lattice_unita(2,1)-lattice_unita(2,2))**2
1680     >   + (lattice_unita(3,1)-lattice_unita(3,2))**2
1681      gamma = (a*a + b*b - d2)/(2.0d0*a*b)
1682      gamma = dacos(gamma)*180.0d0/pi
1683
1684      return
1685      end
1686
1687
1688
1689*     *************************************
1690*     *                                   *
1691*     *     lattice_small_abc_abg         *
1692*     *                                   *
1693*     *************************************
1694*
1695*     This routine computes a,b,c,alpha,beta,gamma.
1696*
1697      subroutine lattice_small_abc_abg(a,b,c,alpha,beta,gamma)
1698      implicit none
1699      real*8 a,b,c
1700      real*8 alpha,beta,gamma
1701
1702*     *** local variables ****
1703      real*8 d2,pi
1704
1705*     **** external functions ****
1706      real*8   lattice_unita_small
1707      external lattice_unita_small
1708
1709*     **** determine a,b,c,alpha,beta,gmma ***
1710      pi = 4.0d0*datan(1.0d0)
1711      a = dsqrt(lattice_unita_small(1,1)**2
1712     >        + lattice_unita_small(2,1)**2
1713     >        + lattice_unita_small(3,1)**2)
1714      b = dsqrt(lattice_unita_small(1,2)**2
1715     >        + lattice_unita_small(2,2)**2
1716     >        + lattice_unita_small(3,2)**2)
1717      c = dsqrt(lattice_unita_small(1,3)**2
1718     >        + lattice_unita_small(2,3)**2
1719     >        + lattice_unita_small(3,3)**2)
1720
1721      d2 = (lattice_unita_small(1,2)-lattice_unita_small(1,3))**2
1722     >   + (lattice_unita_small(2,2)-lattice_unita_small(2,3))**2
1723     >   + (lattice_unita_small(3,2)-lattice_unita_small(3,3))**2
1724      alpha = (b*b + c*c - d2)/(2.0d0*b*c)
1725      alpha = dacos(alpha)*180.0d0/pi
1726
1727      d2 = (lattice_unita_small(1,3)-lattice_unita_small(1,1))**2
1728     >   + (lattice_unita_small(2,3)-lattice_unita_small(2,1))**2
1729     >   + (lattice_unita_small(3,3)-lattice_unita_small(3,1))**2
1730      beta = (c*c + a*a - d2)/(2.0d0*c*a)
1731      beta = dacos(beta)*180.0d0/pi
1732
1733      d2 = (lattice_unita_small(1,1)-lattice_unita_small(1,2))**2
1734     >   + (lattice_unita_small(2,1)-lattice_unita_small(2,2))**2
1735     >   + (lattice_unita_small(3,1)-lattice_unita_small(3,2))**2
1736      gamma = (a*a + b*b - d2)/(2.0d0*a*b)
1737      gamma = dacos(gamma)*180.0d0/pi
1738
1739      return
1740      end
1741
1742
1743
1744
1745