1* 2* $Id$ 3* 4 5 6* ************************************** 7* * * 8* * kbpp_band_stress * 9* * * 10* ************************************** 11 12 logical function kbpp_band_stress(oprint_in,version, 13 > psp_filename,formatted_filename, 14 > ngrid,unita,locp,lmax, 15 > nbrillioun,kvectors) 16 implicit none 17#include "errquit.fh" 18#include "bafdecls.fh" 19#include "tcgmsg.fh" 20#include "msgtypesf.h" 21#include "util.fh" 22 23 logical oprint_in 24 integer version 25 character*50 psp_filename,formatted_filename 26 integer ngrid(3) 27 double precision unita(3,3) 28 integer locp,lmax 29 integer nbrillioun 30 real*8 kvectors(3,*) 31 32* **** local variables **** 33 character*255 full_filename 34 35 integer taskid,MASTER,msglen 36 parameter (MASTER=0) 37 38* **** 1d pseudopotential data **** 39 character*2 atom 40 character*80 comment 41 double precision zv,amass 42 integer nb,lmax0,lmmax,lmax1,locp1 43 double precision rc(0:9),rlocal1 44 integer nrho 45 double precision drho 46 integer rho_indx,vp_indx,wp_indx,sc_r_indx,sc_k_indx 47 integer rho_hndl,vp_hndl,wp_hndl,sc_r_hndl,sc_k_hndl 48 integer isemicore 49 logical semicore 50 double precision rcore 51 52 integer f_indx,cs_indx,sn_indx 53 integer f_hndl,cs_hndl,sn_hndl 54 55* ***** ngrid data ***** 56 integer vl_indx,vnl_indx,G_indx 57 integer vl_hndl,vnl_hndl,G_hndl 58 59* **** ray data **** 60 integer nray,G_ray_hndl,tmp_ray_hndl 61 integer rho_sc_k_ray_hndl,dvnl_ray_hndl,dvl_ray_hndl 62 integer G_ray_indx,tmp_ray_indx 63 integer rho_sc_k_ray_indx,dvnl_ray_indx,dvl_ray_indx 64 65* **** other variables **** 66 logical value,mprint,hprint,oprint,filter,kray 67 double precision unitg(3,3) 68 integer nsize,i,l,ierr,psp_type 69 integer nfft1,nfft2,nfft3 70 71* **** external functions **** 72 logical control_print,control_kbpp_ray,control_kbpp_filter 73 external control_print,control_kbpp_ray,control_kbpp_filter 74 integer kbpp_band_calc_nray 75 external kbpp_band_calc_nray 76 77 78c call Parallel_init() 79 call Parallel_taskid(taskid) 80 hprint = (taskid.eq.MASTER).and.control_print(print_high) 81 mprint = (taskid.eq.MASTER).and.control_print(print_medium) 82 oprint = (oprint_in.or.hprint) 83 84 85 value = .false. 86 87* ***** read in pseudopotential data **** 88 89 if (taskid.eq.MASTER) then 90 call util_file_name_noprefix(psp_filename,.false.,.false., 91 > full_filename) 92 l = index(full_filename,' ') - 1 93 open(unit=11,file=full_filename(1:l), 94 > status='old',form='formatted') 95 96 read(11,'(A2)') atom 97 read(11,*) zv,amass,lmax0,lmax1,locp1,rlocal1 98 read(11,*) (rc(i),i=0,lmax0) 99 read(11,*) nrho,drho 100 read(11,'(A)') comment 101 end if 102 103 msglen = 1 104 call BRDCST(9+MSGDBL,zv,mdtob(msglen),MASTER) 105 call BRDCST(9+MSGDBL,amass,mdtob(msglen),MASTER) 106 call BRDCST(9+MSGINT,lmax0,mitob(msglen),MASTER) 107 call BRDCST(9+MSGINT,lmax1,mitob(msglen),MASTER) 108 call BRDCST(9+MSGINT,locp1,mitob(msglen),MASTER) 109 msglen = lmax0+1 110 call BRDCST(9+MSGDBL,rc,mdtob(msglen),MASTER) 111 msglen = 1 112 call BRDCST(9+MSGINT,nrho,mitob(msglen),MASTER) 113 call BRDCST(9+MSGDBL,drho,mdtob(msglen),MASTER) 114 115 116* **** set the maximum angular momentum **** 117 if (lmax.eq.-1) lmax = lmax1 118 if (lmax.gt.lmax0) lmax = lmax0 119 if (lmax.lt.0) lmax = lmax0 120 121* **** set the local potential **** 122 if (locp.eq.-1) locp = locp1 123 if (locp.gt.lmax) locp = lmax 124 if (locp.lt.0) locp = lmax 125 126* **** allocate rho, vp, and wp **** 127 value = BA_alloc_get(mt_dbl,nrho, 128 > 'rho',rho_hndl,rho_indx) 129 value = value.and.BA_alloc_get(mt_dbl,nrho*(lmax+1), 130 > 'vp',vp_hndl, vp_indx) 131 value = value.and.BA_alloc_get(mt_dbl,nrho*(lmax+1), 132 > 'wp', wp_hndl, wp_indx) 133 value = value.and.BA_alloc_get(mt_dbl,2*nrho, 134 > 'sc', sc_r_hndl, sc_r_indx) 135 if (.not.value) 136 > call errquit('kbpp_band_stress:out of heap memory',0,MA_ERR) 137 138 if (taskid.eq.MASTER) then 139 call read_vpwp_band(11,nrho,lmax,dbl_mb(rho_indx), 140 > dbl_mb(vp_indx), 141 > dbl_mb(wp_indx)) 142 call read_semicore_band(11,isemicore,rcore,nrho,dbl_mb(sc_r_indx)) 143 close(11) 144 end if 145 146 msglen = nrho 147 call BRDCST(9+MSGDBL,dbl_mb(rho_indx),mdtob(msglen),MASTER) 148 msglen = nrho*(lmax+1) 149 call BRDCST(9+MSGDBL,dbl_mb(vp_indx),mdtob(msglen),MASTER) 150 call BRDCST(9+MSGDBL,dbl_mb(wp_indx),mdtob(msglen),MASTER) 151 152 msglen = 1 153 call BRDCST(9+MSGINT,isemicore,mitob(msglen),MASTER) 154 semicore = (isemicore.eq.1) 155 if (semicore) then 156 msglen = 2*nrho 157 call BRDCST(9+MSGDBL,dbl_mb(sc_r_indx),mdtob(msglen),MASTER) 158 else 159 rcore = 0.0d0 160 end if 161 162 163* **** more temporary space **** 164 value = BA_alloc_get(mt_dbl,nrho, 165 > 'f',f_hndl,f_indx) 166 value = value.and.BA_alloc_get(mt_dbl,nrho, 167 > 'cs',cs_hndl,cs_indx) 168 value = value.and.BA_alloc_get(mt_dbl,nrho, 169 > 'sn',sn_hndl,sn_indx) 170 if (.not.value) 171 > call errquit('kbpp_band_stress:out of heap memory',0,MA_ERR) 172 173* **** allocate vl,vnl,vnlnrm G **** 174 lmmax = (lmax+1)**2 - (2*locp+1) 175 nsize = ngrid(1)*ngrid(2)*ngrid(3) 176 value = BA_alloc_get(mt_dbl,nsize, 177 > 'vl',vl_hndl,vl_indx) 178 value = value.and.BA_alloc_get(mt_dbl,nsize*(lmmax)*3, 179 > 'vnl',vnl_hndl, vnl_indx) 180 value = value.and.BA_alloc_get(mt_dbl,nsize*(3), 181 > 'G',G_hndl, G_indx) 182 value = value.and.BA_alloc_get(mt_dbl,4*nsize, 183 > 'sc_k',sc_k_hndl,sc_k_indx) 184 if (.not.value) 185 > call errquit('kbpp_band_stress:out of heap memory',0,MA_ERR) 186 187 188* **** preparation of constants **** 189 nfft1=ngrid(1) 190 nfft2=ngrid(2) 191 nfft3=ngrid(3) 192 call setup_kbpp_band(nfft1,nfft2,nfft3,unita,unitg,dbl_mb(G_indx)) 193 filter = control_kbpp_filter() 194 kray = control_kbpp_ray() 195 196 if (kray) then 197 !**** allocate memory for rays **** 198 nray = kbpp_band_calc_nray(nfft1,nfft2,nfft3,dbl_mb(G_indx)) 199 200 value = BA_alloc_get(mt_dbl,nray, 201 > 'G_ray2',G_ray_hndl,G_ray_indx) 202 value = value.and.BA_alloc_get(mt_dbl,2*nray, 203 > 'dvl_ray2',dvl_ray_hndl,dvl_ray_indx) 204 value = value.and.BA_alloc_get(mt_dbl,4*nray*(lmax+1), 205 > 'dvnl_ray2',dvnl_ray_hndl,dvnl_ray_indx) 206 value = value.and.BA_alloc_get(mt_dbl,2*nray, 207 > 'rho_sc_k_ray2',rho_sc_k_ray_hndl,rho_sc_k_ray_indx) 208 value = value.and.BA_alloc_get(mt_dbl,nray, 209 > 'tmp_ray2',tmp_ray_hndl,tmp_ray_indx) 210 if (.not.value) 211 > call errquit('kbpp_band_stress:out of heap memory',0,MA_ERR) 212 213 call kbpp_band_generate_G_ray(nfft1,nfft2,nfft3, 214 > dbl_mb(G_indx), 215 > dbl_mb(G_ray_indx)) 216 217 call integrate_kbpp_band_stress_local_new(version, 218 > nrho,drho,lmax,locp,zv, 219 > dbl_mb(vp_indx), 220 > dbl_mb(wp_indx), 221 > dbl_mb(rho_indx), 222 > dbl_mb(f_indx), 223 > dbl_mb(cs_indx), 224 > dbl_mb(sn_indx), 225 > nfft1,nfft2,nfft3,lmmax, 226 > dbl_mb(G_indx), 227 > dbl_mb(vl_indx), 228 > semicore, 229 > dbl_mb(sc_r_indx), 230 > dbl_mb(sc_k_indx), 231 > nray, 232 > dbl_mb(G_ray_indx), 233 > dbl_mb(dvl_ray_indx), 234 > dbl_mb(dvnl_ray_indx), 235 > dbl_mb(rho_sc_k_ray_indx), 236 > dbl_mb(tmp_ray_indx), 237 > filter, 238 > ierr) 239 240 241 else 242 call integrate_kbpp_band_stress_local(version, 243 > nrho,drho,lmax,locp,zv, 244 > dbl_mb(vp_indx), 245 > dbl_mb(wp_indx), 246 > dbl_mb(rho_indx), 247 > dbl_mb(f_indx), 248 > dbl_mb(cs_indx), 249 > dbl_mb(sn_indx), 250 > nfft1,nfft2,nfft3,lmmax, 251 > dbl_mb(G_indx), 252 > dbl_mb(vl_indx), 253 > semicore, 254 > dbl_mb(sc_r_indx), 255 > dbl_mb(sc_k_indx), 256 > ierr) 257 end if 258 259 260 if ((taskid.eq.MASTER).and.(oprint)) then 261 write(*,*) " ************************************************" 262 write(*,*) " * *" 263 write(*,*) " * KBPP_Band_Stress - Pseudopotential Formatter *" 264 write(*,*) " * *" 265 write(*,*) " * version last updated 4/15/99 *" 266 write(*,*) " * *" 267 write(*,*) " ************************************************" 268 call nwpw_message(1) 269 write(*,*) 270 write(*,*) "Pseudpotential Data" 271 write(*,*) "-------------------" 272 write(*,*) " atom :",atom 273 write(*,*) " charge :",zv 274 write(*,*) " mass no. :",amass 275 write(*,*) " highest angular component :",lmax0 276 write(*,*) " highest angular component used :",lmax 277 write(*,*) " local potential used :",locp 278 279 if (semicore) then 280 write(*,*) 281 write(*,115) " semi-core stress derivative added" 282 end if 283 write(*,*) 284 write(*,*) "Simulation Cell" 285 write(*,*) "---------------" 286 if (version.eq.3) write(*,112) " boundry: periodic" 287 write(*,113) " ngrid :",ngrid 288 write(*,114) " unita :",unita(1,1),unita(2,1),unita(3,1) 289 write(*,114) " ",unita(1,2),unita(2,2),unita(3,2) 290 write(*,114) " ",unita(1,3),unita(2,3),unita(3,3) 291 write(*,*) 292 111 format(a,10f10.3) 293 112 format(a) 294 113 format(a,3I4) 295 114 format(a,3F10.3) 296 115 format(a,2E14.6) 297 end if 298 299 300 301 if (taskid.eq.MASTER) then 302 call util_file_name_noprefix(formatted_filename, 303 > .false., 304 > .false., 305 > full_filename) 306 l = index(full_filename,' ') - 1 307 if (mprint) then 308 write(*,*) 309 write(*,*) "Generated formatted_stress_filename: ", 310 > full_filename(1:l) 311 if (kray)write(*,'(A,I5,A)')" - Spline fitted, nray=",nray," -" 312 if (filter) write(*,'(A)') " - filtered -" 313 end if 314 call openfile(2,full_filename,l,'w',l) 315 call iwrite(2,version,1) 316 call iwrite(2,ngrid,3) 317 call dwrite(2,unita,9) 318 319 call iwrite(2,nbrillioun,1) 320 call dwrite(2,kvectors,3*nbrillioun) 321 322 call dwrite(2,dbl_mb(vl_indx),nsize) 323 end if 324 325 326 do nb=1,nbrillioun 327 328 if ((oprint).and.(taskid.eq.MASTER)) 329 > write(*,*) "generating brillioun #",nb, 330 > kvectors(1,nb), 331 > kvectors(2,nb), 332 > kvectors(3,nb) 333 334 if (kray) then 335 call integrate_kbpp_band_stress_nonlocal_new(version, 336 > kvectors(1,nb), 337 > nrho,drho,lmax,locp,zv, 338 > dbl_mb(vp_indx), 339 > dbl_mb(wp_indx), 340 > dbl_mb(rho_indx), 341 > dbl_mb(f_indx), 342 > dbl_mb(cs_indx), 343 > dbl_mb(sn_indx), 344 > nfft1,nfft2,nfft3,lmmax, 345 > dbl_mb(G_indx), 346 > dbl_mb(vnl_indx), 347 > nray, 348 > dbl_mb(G_ray_indx), 349 > dbl_mb(dvnl_ray_indx), 350 > ierr) 351 352 else 353 call integrate_kbpp_band_stress_nonlocal(version, 354 > kvectors(1,nb), 355 > nrho,drho,lmax,locp,zv, 356 > dbl_mb(vp_indx), 357 > dbl_mb(wp_indx), 358 > dbl_mb(rho_indx), 359 > dbl_mb(f_indx), 360 > dbl_mb(cs_indx), 361 > dbl_mb(sn_indx), 362 > nfft1,nfft2,nfft3,lmmax, 363 > dbl_mb(G_indx), 364 > dbl_mb(vnl_indx), 365 > ierr) 366 367 end if 368 if (taskid.eq.MASTER) 369 > call dwrite(2,dbl_mb(vnl_indx),nsize*lmmax*3) 370 end do 371 372 373 if (taskid.eq.MASTER) then 374 if (semicore) then 375 call dwrite(2,dbl_mb(sc_k_indx),4*nsize) 376 end if 377 call closefile(2) 378 end if 379 380 381* **** free ray heap space **** 382 if (kray) then 383 value = BA_free_heap(tmp_ray_hndl) 384 value = value.and.BA_free_heap(rho_sc_k_ray_hndl) 385 value = value.and.BA_free_heap(dvnl_ray_hndl) 386 value = value.and.BA_free_heap(dvl_ray_hndl) 387 value = value.and.BA_free_heap(G_ray_hndl) 388 if (.not.value) 389 > call errquit('kbpp_band_stress:Error freeing memory',0,MA_ERR) 390 end if 391 392* **** free heap space **** 393 value = BA_free_heap(rho_hndl) 394 value = value.and.BA_free_heap(vp_hndl) 395 value = value.and.BA_free_heap(wp_hndl) 396 value = value.and.BA_free_heap(sc_r_hndl) 397 value = value.and.BA_free_heap(sc_k_hndl) 398 value = value.and.BA_free_heap(f_hndl) 399 value = value.and.BA_free_heap(cs_hndl) 400 value = value.and.BA_free_heap(sn_hndl) 401 402 403 value = value.and.BA_free_heap(vl_hndl) 404 value = value.and.BA_free_heap(vnl_hndl) 405 value = value.and.BA_free_heap(G_hndl) 406 if (.not.value) 407 > call errquit('kbpp_band_stress:Error freeing heap',0,MA_ERR) 408 409c call Parallel_Finalize() 410 411 if ((taskid.eq.MASTER).and.(oprint)) call nwpw_message(4) 412 413 kbpp_band_stress = value 414 return 415 416 9999 kbpp_band_stress = .false. 417 call errquit('kbpp_band_stress:Error reading psp_filename', 418 > 0,DISK_ERR) 419 END 420 421 422 423 424