1* 2* $Id$ 3* 4 5* ******************************************* 6* * * 7* * wvfnc_expander * 8* * * 9* ******************************************* 10 11 logical function wvfnc_expander(rtdb) 12 implicit none 13 integer rtdb 14 15#include "bafdecls.fh" 16#include "btdb.fh" 17#include "stdio.fh" 18#include "util.fh" 19 20 logical value 21 integer version 22 23 integer ierr 24 25 integer ne(2),ispin 26 27 character*50 new_wavefunction_filename 28 character*50 old_wavefunction_filename 29 character*255 full_filename,full_filename2 30 31 32 integer ngrid(3) 33 integer dngrid(3) 34 integer cfull(2),dcfull(2),occ1(2) 35 integer nfft3d,n2ft3d 36 integer dnfft3d,dn2ft3d 37 integer ms,n,l,occupation 38 double precision unita(3,3) 39 40 logical control_print 41 external control_print 42 43 value = .false. 44 version = 3 45 46* **** get wavefunction information **** 47 value = btdb_cget(rtdb,'xpndr:old_wavefunction_filename', 48 > 1,old_wavefunction_filename) 49 value = btdb_cget(rtdb,'xpndr:new_wavefunction_filename', 50 > 1,new_wavefunction_filename) 51 52 value = btdb_get(rtdb,'xpndr:ngrid',mt_int,3,dngrid) 53 54 55 call util_file_name_noprefix(old_wavefunction_filename, 56 > .false., 57 > .false., 58 > full_filename) 59 60 l = index(full_filename,' ') - 1 61 call openfile(5,full_filename,l,'r',l) 62 call iread(5,version,1) 63 call iread(5,ngrid,3) 64 call dread(5,unita,9) 65 call iread(5,ispin,1) 66 call iread(5,ne,2) 67 call iread(5,occupation,1) 68 69 call util_file_name('wvfnc_expander', 70 > .true., 71 > .false., 72 > full_filename) 73 l = index(full_filename,' ') - 1 74 call openfile(6,full_filename,l,'w',l) 75 call iwrite(6,version,1) 76 call iwrite(6,dngrid,3) 77 call dwrite(6,unita,9) 78 call iwrite(6,ispin,1) 79 call iwrite(6,ne,2) 80 call iwrite(6,occupation,1) 81 82 83 nfft3d = ( ngrid(1)/2+1)* ngrid(2)* ngrid(3) 84 dnfft3d = (dngrid(1)/2+1)*dngrid(2)*dngrid(3) 85 n2ft3d = 2* nfft3d 86 dn2ft3d = 2*dnfft3d 87 if (control_print(print_medium)) then 88 write(luout,109) old_wavefunction_filename 89 write(luout,110) new_wavefunction_filename 90 write(luout,111) ngrid(1), ngrid(2), ngrid(3), 91 > dngrid(1),dngrid(2),dngrid(3) 92 end if 93 109 format(' old_filename: ',A) 94 110 format(' new_filename: ',A) 95 111 format(' converting : ',I3,'x',I3,'x',I3,' --> ', 96 > I3,'x',I3,'x',I3) 97 98* ***** allocate wavefunction memory **** 99 value = BA_alloc_get(mt_dcpl,nfft3d, 100 > 'cfull',cfull(2),cfull(1)) 101 102 value = BA_alloc_get(mt_dcpl,dnfft3d, 103 > 'dcfull',dcfull(2),dcfull(1)) 104 value = BA_alloc_get(mt_dbl,ne(1)+ne(2), 105 > 'occ1',occ1(2),occ1(1)) 106 107 do ms=1,ispin 108 do n=1,ne(ms) 109 call dread(5,dcpl_mb(cfull(1)),n2ft3d) 110 if (control_print(print_medium)) 111 > write(luout,'(A,I5,A,I2)') 112 > "converting .... psi:", n," spin:",ms 113 call wvfnc_expander_convert(ngrid,dcpl_mb(cfull(1)), 114 > dngrid,dcpl_mb(dcfull(1))) 115 116 call dwrite(6,dcpl_mb(dcfull(1)),dn2ft3d) 117 118 end do 119 end do 120 if (occupation.gt.0) then 121 call dread(5,dbl_mb(occ1(1)),(ne(1)+ne(2))) 122 call dwrite(6,dbl_mb(occ1(1)),(ne(1)+ne(2))) 123 end if 124 call closefile(5) 125 call closefile(6) 126 127c *** copy wvfnc_expander to new_wavefunction_filename **** 128 call util_file_name_noprefix(new_wavefunction_filename, 129 > .false., 130 > .false., 131 > full_filename2) 132 call util_file_copy(full_filename,full_filename2) 133 call util_file_unlink(full_filename) 134 IERR=0 135 GO TO 9999 136 137 9110 IERR=10 138 GO TO 9999 139 9111 IERR=11 140 GO TO 9999 141 142 9999 value = BA_free_heap(cfull(2)) 143 value = BA_free_heap(dcfull(2)) 144 value = BA_free_heap(occ1(2)) 145 !IF(IERR.EQ.0) THEN 146 ! WRITE(6,*) ' JOB HAS BEEN COMPLETED. CODE=',IERR 147 !ELSE 148 ! WRITE(6,*) ' JOB HAS BEEN TERMINATED DUE TO CODE=',IERR 149 ! value = .false. 150 !ENDIF 151 !call nwpw_message(4) 152 153 wvfnc_expander = value 154 return 155 end 156 157 158 159* *************************************************** 160* * * 161* * wvfnc_expander_convert * 162* * * 163* *************************************************** 164 165 subroutine wvfnc_expander_convert(ngrid,psi1,dngrid,psi2) 166 implicit none 167 integer ngrid(3) 168 complex*16 psi1(*) 169 integer dngrid(3) 170 complex*16 psi2(*) 171 172* **** local variables **** 173 integer nfft3d,dnfft3d,n2ft3d,dn2ft3d 174 integer inc2,inc3,dinc2,dinc3 175 integer i,j,k,i2,j2,k2,n1,n2,n3 176 integer jdiff,kdiff,indx,dindx 177 logical jreverse,kreverse 178 179 nfft3d = ( ngrid(1)/2+1)* ngrid(2)* ngrid(3) 180 dnfft3d = (dngrid(1)/2+1)*dngrid(2)*dngrid(3) 181 n2ft3d = 2* nfft3d 182 dn2ft3d = 2*dnfft3d 183 inc2 = ( ngrid(1)/2+1) 184 dinc2 = (dngrid(1)/2+1) 185 inc3 = ( ngrid(1)/2+1)* ngrid(2) 186 dinc3 = (dngrid(1)/2+1)*dngrid(2) 187 188 189 n1 = ngrid(1) 190 n2 = ngrid(2) 191 n3 = ngrid(3) 192 if (n1.gt.dngrid(1)) n1 = dngrid(1) 193 if (n2.gt.dngrid(2)) n2 = dngrid(2) 194 if (n3.gt.dngrid(3)) n3 = dngrid(3) 195 196 jdiff = dngrid(2) - ngrid(2) 197 kdiff = dngrid(3) - ngrid(3) 198 jreverse = (jdiff.lt.0) 199 kreverse = (kdiff.lt.0) 200 if (jreverse) jdiff = -jdiff 201 if (kreverse) kdiff = -kdiff 202 203 call dcopy(dn2ft3d,0.0d0,0,psi2,1) 204 do k=0,n3-1 205 do j=0,n2-1 206 do i=0,(n1/2) 207 indx = i 208 dindx = i 209 210 if (k.lt. (n3/2)) then 211 k2 = k 212 else 213 k2 = kdiff + k 214 end if 215 216 if (j.lt. (n2/2)) then 217 j2 = j 218 else 219 j2 = jdiff + j 220 end if 221 if (jreverse) then 222 indx = indx + j2*inc2 223 dindx = dindx + j *dinc2 224 else 225 indx = indx + j *inc2 226 dindx = dindx + j2*dinc2 227 end if 228 if (kreverse) then 229 indx = indx + k2*inc3 230 dindx = dindx + k *dinc3 231 else 232 indx = indx + k *inc3 233 dindx = dindx + k2*dinc3 234 end if 235 236 psi2(dindx+1) = psi1(indx+1) 237 end do 238 end do 239 end do 240 241 return 242 end 243 244 245