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