1*
2* $Id$
3*
4
5*     *******************************************
6*     *						*
7*     *	 	   wvfnc_expand_cell		*
8*     *						*
9*     *******************************************
10
11      logical function wvfnc_expand_cell(rtdb)
12      implicit none
13      integer rtdb
14
15#include "bafdecls.fh"
16#include "btdb.fh"
17
18      logical value
19      integer version
20
21      integer ierr
22
23      integer ne(2),ispin,dne(2)
24
25      character*50 new_wavefunction_filename,filename
26      character*50 old_wavefunction_filename
27      character*255 full_filename,full_filename2
28
29
30      integer ngrid(3)
31      integer dngrid(3)
32      integer cell_expand(3),i,j,k
33      integer cfull(2),dcfull(2)
34      integer nfft3d,n2ft3d
35      integer dnfft3d,dn2ft3d
36      integer ms,n,l,occupation
37
38
39      double precision unita(3,3)
40
41      value = .false.
42      version = 3
43
44*     **** get wavefunction information ****
45      value = btdb_cget(rtdb,'xpndr:old_wavefunction_filename',
46     >                  1,old_wavefunction_filename)
47      value = btdb_cget(rtdb,'xpndr:new_wavefunction_filename',
48     >                  1,new_wavefunction_filename)
49
50      value = btdb_get(rtdb,'nwpw:cell_expand',mt_int,3,cell_expand)
51
52
53      call util_file_name_noprefix(old_wavefunction_filename,
54     >                    .false.,
55     >                    .false.,
56     >                    full_filename)
57
58      l = index(full_filename,' ') - 1
59      call openfile(5,full_filename,l,'r',l)
60      call iread(5,version,1)
61      call iread(5,ngrid,3)
62      call dread(5,unita,9)
63      call iread(5,ispin,1)
64      call iread(5,ne,2)
65      call iread(5,occupation,1)
66
67      dngrid(1) = cell_expand(1)*ngrid(1)
68      dngrid(2) = cell_expand(2)*ngrid(2)
69      dngrid(3) = cell_expand(3)*ngrid(3)
70      dne(1) = cell_expand(1)*cell_expand(2)*cell_expand(3)*ne(1)
71      dne(2) = cell_expand(1)*cell_expand(2)*cell_expand(3)*ne(2)
72      unita(1,1) = unita(1,1)*cell_expand(1)
73      unita(2,1) = unita(2,1)*cell_expand(1)
74      unita(3,1) = unita(3,1)*cell_expand(1)
75      unita(1,2) = unita(1,2)*cell_expand(2)
76      unita(2,2) = unita(2,2)*cell_expand(2)
77      unita(3,2) = unita(3,2)*cell_expand(2)
78      unita(1,3) = unita(1,3)*cell_expand(3)
79      unita(2,3) = unita(2,3)*cell_expand(3)
80      unita(3,3) = unita(3,3)*cell_expand(3)
81
82      l = index(new_wavefunction_filename,' ') - 1
83      filename = new_wavefunction_filename(1:l)//".expander"
84      call util_file_name(filename,
85     >                    .true.,
86     >                    .false.,
87     >                    full_filename)
88      l = index(full_filename,' ') - 1
89      call openfile(6,full_filename,l,'w',l)
90      call iwrite(6,version,1)
91      call iwrite(6,dngrid,3)
92      call dwrite(6,unita,9)
93      call iwrite(6,ispin,1)
94      call iwrite(6,dne,2)
95      call iwrite(6,occupation,1)
96
97
98       nfft3d = ( ngrid(1)/2+1)* ngrid(2)* ngrid(3)
99      dnfft3d = (dngrid(1)/2+1)*dngrid(2)*dngrid(3)
100       n2ft3d = 2* nfft3d
101      dn2ft3d = 2*dnfft3d
102
103      write(*,109) old_wavefunction_filename
104      write(*,110) new_wavefunction_filename
105      write(*,111) ngrid(1), ngrid(2), ngrid(3),
106     >            dngrid(1),dngrid(2),dngrid(3)
107  109 format(' old_filename: ',A)
108  110 format(' new_filename: ',A)
109  111 format(' converting  : ',I3,'x',I3,'x',I3,' --> ',
110     >                     I3,'x',I3,'x',I3)
111
112*     ***** allocate wavefunction memory ****
113      value = BA_alloc_get(mt_dcpl,nfft3d,
114     >                     'cfull',cfull(2),cfull(1))
115
116      value = BA_alloc_get(mt_dcpl,dnfft3d,
117     >                     'dcfull',dcfull(2),dcfull(1))
118
119      do ms=1,ispin
120        do n=1,ne(ms)
121          call dread(5,dcpl_mb(cfull(1)),n2ft3d)
122
123          write(*,'(A,I5,A,I2)') "converting .... psi:", n," spin:",ms
124          do k=0,cell_expand(3)-1
125          do j=0,cell_expand(2)-1
126          do i=0,cell_expand(1)-1
127          write(*,*) "           ....    :", i,j,k
128          call wvfnc_expand_cell_convert(i,j,k,
129     >                               ngrid,dcpl_mb(cfull(1)),
130     >                               dngrid,dcpl_mb(dcfull(1)))
131
132          call dwrite(6,dcpl_mb(dcfull(1)),dn2ft3d)
133          end do
134          end do
135          end do
136
137        end do
138      end do
139      call closefile(5)
140      call closefile(6)
141
142c     *** copy wvfnc_expander to new_wavefunction_filename ****
143      call util_file_name_noprefix(new_wavefunction_filename,
144     >                    .false.,
145     >                    .false.,
146     >                    full_filename2)
147      call util_file_copy(full_filename,full_filename2)
148      call util_file_unlink(full_filename)
149      IERR=0
150      GO TO 9999
151
152 9110 IERR=10
153      GO TO 9999
154 9111 IERR=11
155      GO TO 9999
156
157 9999 value = BA_free_heap(cfull(2))
158      value = BA_free_heap(dcfull(2))
159      !IF(IERR.EQ.0) THEN
160      !  WRITE(6,*) ' JOB HAS BEEN COMPLETED.  CODE=',IERR
161      !ELSE
162      !  WRITE(6,*) ' JOB HAS BEEN TERMINATED DUE TO CODE=',IERR
163      !  value = .false.
164      !ENDIF
165      !call nwpw_message(4)
166
167      wvfnc_expand_cell = value
168      return
169      end
170
171
172*     ***************************************************
173*     *							*
174*     *	 	   wvfnc_expand_cell_convert   		*
175*     *							*
176*     ***************************************************
177
178      subroutine wvfnc_expand_cell_convert(ishift,jshift,kshift,
179     >                              ngrid,psi1,dngrid,psi2)
180      implicit none
181      integer ishift,jshift,kshift
182      integer ngrid(3)
183      complex*16 psi1(*)
184      integer dngrid(3)
185      complex*16 psi2(*)
186
187*     **** local variables ****
188      integer nfft3d,dnfft3d,n2ft3d,dn2ft3d
189      integer inc2,inc3,dinc2,dinc3
190      integer nxh,nyh,nzh
191      integer i,j,k
192      integer i1,j1,k1
193      integer i2,j2,k2
194      integer fi,fj,fk
195      integer indx,dindx
196      real*8  pi,gi,gj,gk
197      complex*16 bw1,bw2,bw3
198      complex*16 cw1,cw2,cw3
199
200       nfft3d = ( ngrid(1)/2+1)* ngrid(2)* ngrid(3)
201      dnfft3d = (dngrid(1)/2+1)*dngrid(2)*dngrid(3)
202       n2ft3d = 2* nfft3d
203      dn2ft3d = 2*dnfft3d
204       inc2 = ( ngrid(1)/2+1)
205      dinc2 = (dngrid(1)/2+1)
206       inc3 = ( ngrid(1)/2+1)* ngrid(2)
207      dinc3 = (dngrid(1)/2+1)*dngrid(2)
208
209      nxh = ngrid(1)/2
210      nyh = ngrid(2)/2
211      nzh = ngrid(3)/2
212      fi  = dngrid(1)/ngrid(1)
213      fj  = dngrid(2)/ngrid(2)
214      fk  = dngrid(3)/ngrid(3)
215      gi  = 1.0d0/dble(fi)
216      gj  = 1.0d0/dble(fj)
217      gk  = 1.0d0/dble(fk)
218
219      pi = 4.0d0*datan(1.0d0)
220      cw1=dcmplx(dcos(0.5d0*ishift*pi),-dsin(0.5d0*ishift*pi))
221      cw2=dcmplx(dcos(0.5d0*jshift*pi),-dsin(0.5d0*jshift*pi))
222      cw3=dcmplx(dcos(0.5d0*kshift*pi),-dsin(0.5d0*kshift*pi))
223c     write(*,*) "cws:",cw1,cw2,cw3
224      cw1=dcmplx(dcos(gi*pi*ishift),-dsin(gi*pi*ishift))
225      cw2=dcmplx(dcos(gj*pi*jshift),-dsin(gj*pi*jshift))
226      cw3=dcmplx(dcos(gk*pi*kshift),-dsin(gk*pi*kshift))
227
228      call dcopy(dn2ft3d,0.0d0,0,psi2,1)
229      do k=-nzh+1,nzh-1
230      do j=-nyh+1,nyh-1
231      do i=0,nxh-1
232
233         i1 = i
234         j1 = j
235         k1 = k
236         !if (i1.ge.0) i2 = fi*i1 + ishift
237         !if (i1.lt.0) i2 = fi*i1 - ishift
238         !if (j1.ge.0) j2 = fj*j1 + jshift
239         !if (j1.lt.0) j2 = fj*j1 - jshift
240         !if (k1.ge.0) k2 = fk*k1 + kshift
241         !if (k1.lt.0) k2 = fk*k1 - kshift
242         i2 = fi*i1 + ishift
243         j2 = fj*j1 + jshift
244         k2 = fk*k1 + kshift
245
246
247
248
249         if (i1 .lt. 0) i1 = i1 + ngrid(1)
250         if (j1 .lt. 0) j1 = j1 + ngrid(2)
251         if (k1 .lt. 0) k1 = k1 + ngrid(3)
252
253         if (i2 .lt. 0) i2 = i2 + dngrid(1)
254         if (j2 .lt. 0) j2 = j2 + dngrid(2)
255         if (k2 .lt. 0) k2 = k2 + dngrid(3)
256
257         indx   = (k1)*inc3  +(j1)*inc2  + i1+1
258         dindx  = (k2)*dinc3 +(j2)*dinc2 + i2+1
259
260
261         psi2(dindx) = psi1(indx)*cw1*cw2*cw3
262
263
264c        i1 = i
265c        j1 = -j
266c        k1 = -k
267c        i2 = fi*i1
268c        j2 = fj*j1
269c        k2 = fk*k1
270
271c        if (i1 .lt. 0) i1 = i1 + ngrid(1)
272c        if (j1 .lt. 0) j1 = j1 + ngrid(2)
273c        if (k1 .lt. 0) k1 = k1 + ngrid(3)
274
275c        if (i2 .lt. 0) i2 = i2 + dngrid(1)
276c        if (j2 .lt. 0) j2 = j2 + dngrid(2)
277c        if (k2 .lt. 0) k2 = k2 + dngrid(3)
278
279c        indx   = (k1)*inc3  +(j1)*inc2  + i1+1
280c        dindx  = (k2)*dinc3 +(j2)*dinc2 + i2+1
281
282c        psi2(dindx) = psi1(indx)
283
284
285      end do
286      end do
287      end do
288
289      return
290      end
291
292
293
294