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