1* 2* $Id$ 3* 4 5* ************************************** 6* * * 7* * hghppv1 * 8* * * 9* ************************************** 10 11 logical function hghppv1(oprint_in,version, 12 > psp_filename,formatted_filename, 13 > ngrid,unita,rlocal) 14 implicit none 15 16#include "bafdecls.fh" 17#include "tcgmsg.fh" 18#include "msgtypesf.h" 19#include "errquit.fh" 20#include "util.fh" 21#include "stdio.fh" 22 23 logical oprint_in,filter 24 integer version 25 character*50 psp_filename,formatted_filename 26 integer ngrid(3) 27 double precision unita(3,3) 28 real*8 rlocal 29 character*255 full_filename 30 31 logical value,mprint,hprint,oprint 32 33 integer taskid,MASTER,msglen 34 parameter (MASTER=0) 35 36* **** 1d pseudopotential data **** 37 character*80 comment 38 39 character*2 atom 40 integer Zion 41 double precision zv,amass 42 43 double precision rloc,C1,C2,C3,C4 44 double precision r(0:3),H(3,0:3),K(3,0:3) 45 double precision Gijl(9*4),Gtmp(3,3,0:3) 46 double precision rcore 47 double precision ecut,wcut 48 49 integer lmax,locp 50 51 integer n_prj_indx,l_prj_indx,m_prj_indx,b_prj_indx 52 integer n_prj_hndl,l_prj_hndl,m_prj_hndl,b_prj_hndl 53 54* ***** ngrid data ***** 55 integer Ylm_indx,vl_indx,vnl_indx,G_indx 56 integer Ylm_hndl,vl_hndl,vnl_hndl,G_hndl 57 58* **** other variables **** 59 double precision unitg(3,3) 60 integer nsize,i,j,n,l,m,psp_type,indx 61 integer nfft1,nfft2,nfft3 62 integer nmax,nprj 63 64* **** external functions **** 65 logical control_print,control_kbpp_filter 66 external control_print,control_kbpp_filter 67 real*8 control_ecut,control_wcut 68 external control_ecut,control_wcut 69 70 71 72 73c call Parallel_init() 74 call Parallel_taskid(taskid) 75 hprint = (taskid.eq.MASTER).and.control_print(print_high) 76 mprint = (taskid.eq.MASTER).and.control_print(print_medium) 77 oprint = (oprint_in.or.hprint) 78 79 value = .false. 80 81* ***** read in pseudopotential data **** 82 if (taskid.eq.MASTER) then 83 call util_file_name_noprefix(psp_filename,.false.,.false., 84 > full_filename) 85 l = index(full_filename,' ') - 1 86 open(unit=11,file=full_filename(1:l), 87 > status='old',form='formatted') 88 read(11,*) psp_type 89 read(11,'(A2)') atom 90 read(11,*) Zion 91 read(11,*) lmax 92 read(11,*) rloc,C1,C2,C3,C4 93 94 if (lmax.ge.0) then 95 read(11,*) r(0),H(1,0),H(2,0),H(3,0) 96 do l=1,lmax 97 read(11,*) r(l),H(1,l),H(2,l),H(3,l) 98 read(11,*) K(1,l),K(2,l),K(3,l) 99 end do 100 end if 101 read(11,'(A)') comment 102 close(11) 103 104 105 !**** determine nmax **** 106 nmax = 0 107 do l=0,lmax 108 do i=1,3 109 if ((H(i,l).ne.0.0d0) .and. (i.gt.nmax)) nmax = i 110 end do 111 end do 112 end if 113 114 115 msglen = 1 116 call BRDCST(9+MSGINT,psp_type,mitob(msglen),MASTER) 117 call BRDCST(9+MSGINT,Zion,mitob(msglen),MASTER) 118 call BRDCST(9+MSGINT,lmax,mitob(msglen),MASTER) 119 call BRDCST(9+MSGINT,nmax,mitob(msglen),MASTER) 120 121 call BRDCST(9+MSGDBL,rloc,mdtob(msglen),MASTER) 122 call BRDCST(9+MSGDBL,C1,mdtob(msglen),MASTER) 123 call BRDCST(9+MSGDBL,C2,mdtob(msglen),MASTER) 124 call BRDCST(9+MSGDBL,C3,mdtob(msglen),MASTER) 125 call BRDCST(9+MSGDBL,C4,mdtob(msglen),MASTER) 126 127 msglen = 4 128 call BRDCST(9+MSGDBL,r,mdtob(msglen),MASTER) 129 130 msglen = 12 131 call BRDCST(9+MSGDBL,H,mdtob(msglen),MASTER) 132 call BRDCST(9+MSGDBL,K,mdtob(msglen),MASTER) 133 134 135 136* **** set the maximum angular momentum **** 137 138* **** set the local potential **** 139 locp = lmax+1 140 141* **** set the local potential **** 142 rlocal = 1.0d0 143 144 !**** determine nprj **** 145 nprj= 0 146 do l=0,lmax 147 !write(luout,*) "???H :",l,(H(i,l),i=1,3) 148 do i=1,3 149 if ((H(i,l).ne.0.0d0)) nprj = nprj + (2*l+1) 150 end do 151 end do 152 !write(luout,*) "???nprj:", nprj 153 154 155 156* **** set the projector coeficients **** 157 call dcopy(9*4,0.0d0,0,Gtmp,1) 158 call dcopy(9*4,0.0d0,0,Gijl,1) 159 Gtmp(1,1,0) = H(1,0) 160 Gtmp(2,2,0) = H(2,0) 161 Gtmp(3,3,0) = H(3,0) 162 Gtmp(1,2,0) = -0.5d0*dsqrt(3.0d0/5.0d0) *Gtmp(2,2,0) 163 Gtmp(2,1,0) = Gtmp(1,2,0) 164 Gtmp(1,3,0) = 0.5d0*dsqrt(5.0d0/21.0d0) *Gtmp(3,3,0) 165 Gtmp(3,1,0) = Gtmp(1,3,0) 166 Gtmp(2,3,0) = -0.5d0*dsqrt(100.0d0/63.0d0)*Gtmp(3,3,0) 167 Gtmp(3,2,0) = Gtmp(2,3,0) 168 169 170 Gtmp(1,1,1) = H(1,1) 171 Gtmp(2,2,1) = H(2,1) 172 Gtmp(3,3,1) = H(3,1) 173 Gtmp(1,2,1) = -0.5d0*dsqrt(5.0d0/7.0d0) *Gtmp(2,2,1) 174 Gtmp(2,1,1) = Gtmp(1,2,1) 175 Gtmp(1,3,1) = (1.0d0/6.0d0)*dsqrt(35.0d0/11.0d0)*Gtmp(3,3,1) 176 Gtmp(3,1,1) = Gtmp(1,3,1) 177 Gtmp(2,3,1) = -(14.0d0/6.0d0)*dsqrt(1.0d0/11.0d0)*Gtmp(3,3,1) 178 Gtmp(3,2,1) = Gtmp(2,3,1) 179 180 Gtmp(1,1,2) = H(1,2) 181 Gtmp(2,2,2) = H(2,2) 182 Gtmp(3,3,2) = H(3,2) 183 Gtmp(1,2,2) = -0.5d0*dsqrt(7.0d0/9.0d0) *Gtmp(2,2,2) 184 Gtmp(2,1,2) = Gtmp(1,2,2) 185 Gtmp(1,3,2) = 0.5d0*dsqrt(63.0d0/143.0d0)*Gtmp(3,3,2) 186 Gtmp(3,1,2) = Gtmp(1,3,2) 187 Gtmp(2,3,2) = -0.5d0*18.0d0*dsqrt(1.0d0/143.0d0)*Gtmp(3,3,2) 188 Gtmp(3,2,2) = Gtmp(2,3,2) 189 190 Gtmp(1,1,3) = H(1,3) 191 Gtmp(2,2,3) = H(2,3) 192 Gtmp(3,3,3) = H(3,3) 193 194 195 196 do l=0,lmax 197 do i=1,nmax 198 do j=1,nmax 199 Gijl(i+(j-1)*nmax+l*nmax*nmax) = Gtmp(i,j,l) 200 end do 201 end do 202 end do 203 204 205* **** allocate vl, vnl, G **** 206 nsize = (ngrid(1)/2+1)*ngrid(2)*ngrid(3) 207 value = BA_alloc_get(mt_dbl,nsize, 208 > 'vl',vl_hndl,vl_indx) 209 value = value.and.BA_alloc_get(mt_dbl,nsize*(nprj), 210 > 'vnl',vnl_hndl, vnl_indx) 211 value = value.and.BA_alloc_get(mt_dbl,nsize, 212 > 'Ylm',Ylm_hndl, Ylm_indx) 213 214 value = value.and.BA_alloc_get(mt_dbl,nsize*(3), 215 > 'G',G_hndl, G_indx) 216 value = value.and.BA_alloc_get(mt_int,nprj, 217 > 'n_prj', n_prj_hndl, n_prj_indx) 218 value = value.and.BA_alloc_get(mt_int,nprj, 219 > 'l_prj', l_prj_hndl, l_prj_indx) 220 value = value.and.BA_alloc_get(mt_int,nprj, 221 > 'm_prj', m_prj_hndl, m_prj_indx) 222 value = value.and.BA_alloc_get(mt_int,nprj, 223 > 'b_prj', b_prj_hndl, b_prj_indx) 224 if(.not.value) 225 > call errquit('hghppv1: out of heap memory', 0, MA_ERR) 226 227 228 !**** determine n_prj, l_prj, and m_prj arrays **** 229 indx = 0 230 nfft1 = 1 231 do l=0,lmax 232 do i=1,3 233 if ((H(i,l).ne.0.0d0)) then 234 do m=-l,l 235 int_mb(n_prj_indx+indx) = i 236 int_mb(l_prj_indx+indx) = l 237 int_mb(m_prj_indx+indx) = m 238 int_mb(b_prj_indx+indx) = nfft1 239 indx = indx + 1 240 end do 241 nfft1=nfft1+1 242 end if 243 end do 244 end do 245 246 filter = control_kbpp_filter() 247 ecut = control_ecut() 248 wcut = control_wcut() 249 250* **** preparation of constants **** 251 nfft1=ngrid(1) 252 nfft2=ngrid(2) 253 nfft3=ngrid(3) 254 call setup_kbpp(nfft1,nfft2,nfft3,unita,unitg,dbl_mb(G_indx)) 255 zv = dble(Zion) 256 257 call HGH_local(version,rlocal, 258 > zv,rloc,C1,C2,C3,C4, 259 > nfft1,nfft2,nfft3, 260 > dbl_mb(G_indx), 261 > dbl_mb(vl_indx)) 262 if (filter) 263 > call kbpp_filter(nfft1,nfft2,nfft3, 264 > dbl_mb(G_indx), 265 > ecut, 266 > dbl_mb(vl_indx)) 267 268 do i=1,nprj 269 n=int_mb(n_prj_indx+i-1) 270 l=int_mb(l_prj_indx+i-1) 271 m=int_mb(m_prj_indx+i-1) 272 call Tesseral(l,m, 273 > nfft1,nfft2,nfft3, 274 > dbl_mb(G_indx), 275 > dbl_mb(Ylm_indx)) 276 call HGH_nonlocal(n,l, 277 > r(l), 278 > nfft1,nfft2,nfft3, 279 > dbl_mb(G_indx), 280 > dbl_mb(vnl_indx+(i-1)*nsize)) 281 if ((taskid.eq.MASTER).and.(oprint)) 282 > write(luout,*) "creating projector:",n,l,m 283 do j=1,nsize 284 dbl_mb(vnl_indx+(i-1)*nsize+j-1) 285 > = dbl_mb(vnl_indx+(i-1)*nsize+j-1)*dbl_mb(Ylm_indx+j-1) 286 end do 287 288 if (filter) 289 > call kbpp_filter(nfft1,nfft2,nfft3, 290 > dbl_mb(G_indx), 291 > wcut, 292 > dbl_mb(vnl_indx+(i-1)*nsize)) 293 end do 294 295 296 if ((taskid.eq.MASTER).and.(oprint)) then 297 write(luout,*) " ********************************************" 298 write(luout,*) " * *" 299 write(luout,*) " * HGHPPV1 - Pseudopotential Formatter *" 300 write(luout,*) " * *" 301 write(luout,*) " * version last updated 11/13/03 *" 302 write(luout,*) " * *" 303 write(luout,*) " * developed by Eric J. Bylaska *" 304 write(luout,*) " * *" 305 write(luout,*) " ********************************************" 306 call nwpw_message(1) 307 write(luout,*) 308 write(luout,*) "Pseudpotential Data" 309 write(luout,*) "-------------------" 310 write(luout,*) " atom :",atom 311 write(luout,*) " charge :",Zion 312 write(luout,*) " highest angular component used :",lmax 313 write(luout,*) " highest radial component used :",nmax 314 write(luout,*) " number of non-local projectors :",nprj 315 write(luout,111) " projector cutoffs: ",(r(i), i=0,lmax) 316 write(luout,*) 317 write(luout,111) " local psp cutoff : ",rloc 318 write(luout,111) " local psp coefficients : ",C1,C2,C3,C4 319 if (version.eq.4) 320 > write(luout,*) " aperiodic cutoff radius :",rlocal 321 322 write(luout,*) 323 write(luout,*) "Simulation Cell" 324 write(luout,*) "---------------" 325 if (version.eq.3) write(luout,112) " boundry: periodic" 326 if (version.eq.4) write(luout,112) " boundry: aperiodic" 327 write(luout,113) " ngrid :",ngrid 328 write(luout,114) " unita :",unita(1,1),unita(2,1),unita(3,1) 329 write(luout,114) " ",unita(1,2),unita(2,2),unita(3,2) 330 write(luout,114) " ",unita(1,3),unita(2,3),unita(3,3) 331 write(luout,*) 332 111 format(a,10f10.3) 333 112 format(a) 334 113 format(a,3I4) 335 114 format(a,3F10.3) 336 115 format(a,2E14.6) 337 end if 338 339 340 341 342 if (taskid.eq.MASTER) then 343 call util_file_name_noprefix(formatted_filename, 344 > .false., 345 > .false., 346 > full_filename) 347 l = index(full_filename,' ') - 1 348 if (mprint) then 349 write(luout,*) 350 write(luout,*) "Generated formatted_filename: ",full_filename(1:l) 351 if (filter) write(luout,*) "- filtering pseudopotential -" 352 !write(luout,*) 353 end if 354 call openfile(2,full_filename,l,'w',l) 355 356 call cwrite(2,comment,80) 357 call iwrite(2,psp_type,1) 358 call iwrite(2,version,1) 359 call iwrite(2,ngrid,3) 360 call dwrite(2,unita,9) 361 call cwrite(2,atom,2) 362 call dwrite(2,amass,1) 363 call dwrite(2,zv,1) 364 call iwrite(2,lmax,1) 365 call iwrite(2,locp,1) 366 367 call iwrite(2,nmax,1) 368 call dwrite(2,r,lmax+1) 369 370 371 call iwrite(2,nprj,1) 372 if (nprj.gt.0) then 373 call iwrite(2,int_mb(n_prj_indx),nprj) 374 call iwrite(2,int_mb(l_prj_indx),nprj) 375 call iwrite(2,int_mb(m_prj_indx),nprj) 376 call iwrite(2,int_mb(b_prj_indx),nprj) 377 call dwrite(2,Gijl,(nmax*nmax*(lmax+1))) 378 end if 379 380 if (version.eq.4) call dwrite(2,rlocal,1) 381 rcore = 0.0d0 382 call dwrite(2,rcore,1) 383 384 call dwrite(2,dbl_mb(vl_indx),nsize) 385 if (nprj.gt.0) then 386 call dwrite(2,dbl_mb(vnl_indx),nsize*nprj) 387 end if 388 389 call closefile(2) 390 end if 391 392* **** free heap space **** 393 value = BA_free_heap(vl_hndl) 394 value = value.and.BA_free_heap(vnl_hndl) 395 value = value.and.BA_free_heap(Ylm_hndl) 396 value = value.and.BA_free_heap(G_hndl) 397 value = value.and.BA_free_heap(n_prj_hndl) 398 value = value.and.BA_free_heap(l_prj_hndl) 399 value = value.and.BA_free_heap(m_prj_hndl) 400 value = value.and.BA_free_heap(b_prj_hndl) 401 if(.not.value) 402 > call errquit('hghppv1: deallocatin heap memory', 0, MA_ERR) 403 404 405 if ((taskid.eq.MASTER).and.(oprint)) call nwpw_message(4) 406 hghppv1 = value 407 return 408 409 9999 call errquit('hghppv1:Error reading psp_filename',0, DISK_ERR) 410 411 END 412 413 414 415 416