1* 2* $Id$ 3* 4 5* ************************************** 6* * * 7* * pawppv1 * 8* * * 9* ************************************** 10 11 logical function pawppv1(oprint_in,version, 12 > psp_filename,formatted_filename, 13 > ngrid,unita,locp,lmax,rlocal) 14 implicit none 15 logical oprint_in 16 integer version 17 character*50 psp_filename,formatted_filename 18 integer ngrid(3) 19 double precision unita(3,3) 20 integer locp,lmax 21 real*8 rlocal 22 23#include "errquit.fh" 24#include "bafdecls.fh" 25#include "util.fh" 26#include "stdio.fh" 27 28* ***** local variables **** 29 integer taskid,MASTER,msglen 30 parameter (MASTER=0) 31 32 33* **** 1d pseudopotential data **** 34 integer psp_type 35 integer n1dgrid,nbasis,icut,nmax 36 real*8 zv,r1,rmax,core_kin_energy,core_ion_energy 37 real*8 zion,sigma,rcore,amass,rc(25) 38 real*8 log_amesh,amesh,ecorez,ecorecore 39 character*2 atom 40 character*80 comment 41 42 integer bprj(2),mprj(2),lprj(2),nprj(2),prj_ps0(2),prj_ps(2) 43 integer v_ps(2),core_ps(2),core_ae(2) 44 integer core_ps_prime(2),core_ae_prime(2) 45 integer dphi_ps(2),phi_ps(2),dphi_ae(2),phi_ae(2) 46 integer eig(2),lps(2),nps(2),nae(2) 47 integer rgrid(2) 48 integer nmaxl(2) 49 50* **** matrix data **** 51 integer Gijl(2) 52 integer hartree_matrix(2),comp_charge_matrix(2),comp_pot_matrix(2) 53 54* ***** ngrid data ***** 55 integer nproj,nsize,nfft1,nfft2,nfft3 56 integer vl(2),vlpaw(2),vnl(2),G_indx,G_hndl 57 integer f1(2),f2(2),f3(2),f4(2),cs(2),sn(2) 58 59* **** ray data **** 60 logical filter 61 integer nray,G_ray_hndl,tmp_ray_hndl 62 integer vnl_ray_hndl,vl_ray_hndl,vlpaw_ray_hndl 63 integer G_ray_indx,tmp_ray_indx 64 integer vnl_ray_indx,vl_ray_indx,vlpaw_ray_indx 65 66 67* **** other variables **** 68 logical hprint,mprint,oprint,value 69 integer idum,l,ii,i,j,ierr 70 character*255 full_filename 71 real*8 zcore,fourpi,unitg(3,3) 72 73* **** external functions **** 74 logical control_print,control_kbpp_filter 75 external control_print,control_kbpp_filter 76 real*8 log_integrate_def,log_coulomb0_energy,log_coulomb_energy 77 external log_integrate_def,log_coulomb0_energy,log_coulomb_energy 78 integer kbpp_calc_nray 79 external kbpp_calc_nray 80 81 82 call Parallel_taskid(taskid) 83 hprint = (taskid.eq.MASTER).and.control_print(print_high) 84 mprint = (taskid.eq.MASTER).and.control_print(print_medium) 85 oprint = (oprint_in.or.hprint) 86 87 88 if (taskid.eq.MASTER) then 89 call util_file_name_noprefix(psp_filename,.false.,.false., 90 > full_filename) 91 l = index(full_filename,' ') - 1 92 open(unit=11,file=full_filename(1:l), 93 > status='old',form='formatted') 94 read(11,*) psp_type 95 read(11,'(A)') atom 96 read(11,*) zv 97 read(11,*) r1 98 read(11,*) rmax 99 read(11,*) n1dgrid 100 read(11,*) nbasis 101 read(11,*) (rc(i),i=1,nbasis) 102 read(11,*) icut 103 read(11,'(A)') comment 104 read(11,*) core_kin_energy 105 end if 106 107 msglen = 1 108 call Parallel_Brdcst_ivalues(MASTER,msglen,psp_type) 109 call Parallel_Brdcst_values(MASTER,msglen,zv) 110 call Parallel_Brdcst_values(MASTER,msglen,r1) 111 call Parallel_Brdcst_values(MASTER,msglen,rmax) 112 call Parallel_Brdcst_ivalues(MASTER,msglen,n1dgrid) 113 call Parallel_Brdcst_ivalues(MASTER,msglen,nbasis) 114 call Parallel_Brdcst_ivalues(MASTER,msglen,icut) 115 call Parallel_Brdcst_values(MASTER,msglen,core_kin_energy) 116 msglen = nbasis 117 call Parallel_Brdcst_values(MASTER,msglen,rc) 118 119 120* **** define rgrid **** 121 log_amesh = dlog(rmax/r1)/dble(n1dgrid-1) 122 amesh = dexp(log_amesh) 123 value = BA_alloc_get(mt_dbl,n1dgrid,'rgrid',rgrid(2),rgrid(1)) 124 if (.not.value) call errquit('pawppv1:out of heap',0,MA_ERR) 125 dbl_mb(rgrid(1)) = r1 126 do i=1,n1dgrid-1 127 dbl_mb(rgrid(1)+i) = dbl_mb(rgrid(1)+i-1)*amesh 128 end do 129 130* **** allocate rest of grid data **** 131 value = BA_alloc_get(mt_int,nbasis,'nae',nae(2),nae(1)) 132 value = value.and.BA_alloc_get(mt_int,nbasis,'nps',nps(2),nps(1)) 133 value = value.and.BA_alloc_get(mt_int,nbasis,'lps',lps(2),lps(1)) 134 value = value.and.BA_alloc_get(mt_dbl,nbasis,'eig',eig(2),eig(1)) 135 value = value.and.BA_alloc_get(mt_dbl,nbasis*n1dgrid,'phi_ae', 136 > phi_ae(2),phi_ae(1)) 137 value = value.and.BA_alloc_get(mt_dbl,nbasis*n1dgrid,'dphi_ae', 138 > dphi_ae(2),dphi_ae(1)) 139 value = value.and.BA_alloc_get(mt_dbl,nbasis*n1dgrid,'phi_ps', 140 > phi_ps(2),phi_ps(1)) 141 value = value.and.BA_alloc_get(mt_dbl,nbasis*n1dgrid,'dphi_ps', 142 > dphi_ps(2),dphi_ps(1)) 143 value = value.and.BA_alloc_get(mt_dbl,n1dgrid,'core_ae', 144 > core_ae(2),core_ae(1)) 145 value = value.and.BA_alloc_get(mt_dbl,n1dgrid,'core_ps', 146 > core_ps(2),core_ps(1)) 147 value = value.and.BA_alloc_get(mt_dbl,n1dgrid,'core_ae_prime', 148 > core_ae_prime(2),core_ae_prime(1)) 149 value = value.and.BA_alloc_get(mt_dbl,n1dgrid,'core_ps_prime', 150 > core_ps_prime(2),core_ps_prime(1)) 151 value = value.and.BA_alloc_get(mt_dbl,n1dgrid,'v_ps', 152 > v_ps(2),v_ps(1)) 153 value = value.and.BA_alloc_get(mt_dbl,nbasis*n1dgrid,'prj_ps', 154 > prj_ps(2),prj_ps(1)) 155 value = value.and.BA_alloc_get(mt_dbl,nbasis*n1dgrid,'prj_ps0', 156 > prj_ps0(2),prj_ps0(1)) 157 if (.not.value) call errquit('pawppv1:out of heap',1,MA_ERR) 158 159 160 if (taskid.eq.MASTER) then 161 read(11,*) (int_mb(nae(1)+j), 162 > dbl_mb(eig(1)+j), 163 > int_mb(nps(1)+j), 164 > int_mb(lps(1)+j),j=0,nbasis-1) 165 read(11,*) ((dbl_mb(phi_ae(1)+i+j*n1dgrid), 166 > i=0,n1dgrid-1), 167 > j=0,nbasis-1) 168 read(11,*) ((dbl_mb(dphi_ae(1)+i+j*n1dgrid), 169 > i=0,n1dgrid-1), 170 > j=0,nbasis-1) 171 read(11,*) ((dbl_mb(phi_ps(1)+i+j*n1dgrid), 172 > i=0,n1dgrid-1), 173 > j=0,nbasis-1) 174 read(11,*) ((dbl_mb(dphi_ps(1)+i+j*n1dgrid), 175 > i=0,n1dgrid-1), 176 > j=0,nbasis-1) 177 read(11,*) ((dbl_mb(prj_ps(1)+i+j*n1dgrid), 178 > i=0,n1dgrid-1), 179 > j=0,nbasis-1) 180 read(11,*) (dbl_mb(core_ae(1)+i), 181 > i=0,n1dgrid-1) 182 read(11,*) (dbl_mb(core_ps(1)+i), 183 > i=0,n1dgrid-1) 184 read(11,*) (dbl_mb(v_ps(1)+i), 185 > i=0,n1dgrid-1) 186 read(11,*) sigma 187 read(11,*) zion 188 read(11,*) ((dbl_mb(prj_ps0(1)+i+j*n1dgrid), 189 > i=0,n1dgrid-1), 190 > j=0,nbasis-1) 191 close(11) 192 end if 193 194 195 msglen = nbasis 196 call Parallel_Brdcst_ivalues(MASTER,msglen,int_mb(nae(1))) 197 call Parallel_Brdcst_ivalues(MASTER,msglen,int_mb(nps(1))) 198 call Parallel_Brdcst_ivalues(MASTER,msglen,int_mb(lps(1))) 199 call Parallel_Brdcst_values(MASTER,msglen,dbl_mb(eig(1))) 200 201 msglen = nbasis*n1dgrid 202 call Parallel_Brdcst_values(MASTER,msglen,dbl_mb(phi_ae(1))) 203 call Parallel_Brdcst_values(MASTER,msglen,dbl_mb(dphi_ae(1))) 204 call Parallel_Brdcst_values(MASTER,msglen,dbl_mb(phi_ps(1))) 205 call Parallel_Brdcst_values(MASTER,msglen,dbl_mb(dphi_ps(1))) 206 call Parallel_Brdcst_values(MASTER,msglen,dbl_mb(prj_ps(1))) 207 call Parallel_Brdcst_values(MASTER,msglen,dbl_mb(prj_ps0(1))) 208 209 msglen = n1dgrid 210 call Parallel_Brdcst_values(MASTER,msglen,dbl_mb(core_ae(1))) 211 call Parallel_Brdcst_values(MASTER,msglen,dbl_mb(core_ps(1))) 212 call Parallel_Brdcst_values(MASTER,msglen,dbl_mb(v_ps(1))) 213 214 msglen = 1 215 call Parallel_Brdcst_values(MASTER,msglen,sigma) 216 call Parallel_Brdcst_values(MASTER,msglen,zion) 217 218 219* **** define nproj and lmax **** 220 locp = -1 221 lmax = -1 222 nproj = 0 223 do ii=1,nbasis 224 l = int_mb(lps(1)+ii-1) 225 nproj = nproj + 2*l+1 226 if (l.gt.lmax) lmax = l 227 end do 228 229* **** define nmax **** 230 if(.not. BA_push_get(mt_int,lmax+1,'nmaxl',nmaxl(2),nmaxl(1))) 231 > call errquit('pawppv1:out of stack',2,MA_ERR) 232 call icopy(lmax+1,0,0,int_mb(nmaxl(1)),1) 233 do ii=1,nbasis 234 l = int_mb(lps(1)+ii-1) 235 int_mb(nmaxl(1)+l) = int_mb(nmaxl(1)+l) + 1 236 end do 237 nmax = 0 238 do l=0,lmax 239 if (int_mb(nmaxl(1)+l).gt.nmax) nmax = int_mb(nmaxl(1)+l) 240 end do 241 if(.not.BA_pop_stack(nmaxl(2))) 242 > call errquit('pawppv1:error popping stack',0,MA_ERR) 243 244 245* **** allocate Gijl,Sijl,Tijl,vcore,Vpseuo **** 246 l = nmax*nmax*(lmax+1) 247 value= BA_alloc_get(mt_dbl,5*l,'Gijl',Gijl(2),Gijl(1)) 248 !value=value.and.BA_alloc_get(mt_dbl,l,'Tijl',Tijl(2),Tijl(1)) 249 l = nbasis*nbasis*nbasis*nbasis*(2*lmax+1) 250 value=value.and. 251 > BA_alloc_get(mt_dbl,l,'hartree_matrix', 252 > hartree_matrix(2),hartree_matrix(1)) 253 l = nbasis*nbasis*(2*lmax+1) 254 value=value.and. 255 > BA_alloc_get(mt_dbl,l,'comp_charge_matrix', 256 > comp_charge_matrix(2),comp_charge_matrix(1)) 257 l = nbasis*nbasis*(2*lmax+1) 258 value=value.and. 259 > BA_alloc_get(mt_dbl,l,'comp_pot_matrix', 260 > comp_pot_matrix(2),comp_pot_matrix(1)) 261 262* **** allocate nprj,lprj,mprj,bprj **** 263 value=value.and.BA_alloc_get(mt_int,nproj,'nprj',nprj(2),nprj(1)) 264 value=value.and.BA_alloc_get(mt_int,nproj,'lprj',lprj(2),lprj(1)) 265 value=value.and.BA_alloc_get(mt_int,nproj,'mprj',mprj(2),mprj(1)) 266 value=value.and.BA_alloc_get(mt_int,nproj,'bprj',bprj(2),bprj(1)) 267 if (.not.value) call errquit('pawppv1:out of heap',2,MA_ERR) 268 269 270* **** more temporary space **** 271 value = BA_alloc_get(mt_dbl,n1dgrid,'f1',f1(2),f1(1)) 272 value = value.and.BA_alloc_get(mt_dbl,n1dgrid,'f2',f2(2),f2(1)) 273 value = value.and.BA_alloc_get(mt_dbl,n1dgrid,'f3',f3(2),f3(1)) 274 value = value.and.BA_alloc_get(mt_dbl,n1dgrid,'f4',f4(2),f4(1)) 275 value = value.and.BA_alloc_get(mt_dbl,n1dgrid,'cs',cs(2),cs(1)) 276 value = value.and.BA_alloc_get(mt_dbl,n1dgrid,'sn',sn(2),sn(1)) 277 if (.not.value)call errquit('pawppv1:out of heap',0,MA_ERR) 278 279* **** allocate vl,vnl,vnlnrm G **** 280 nsize = (ngrid(1)/2+1)*ngrid(2)*ngrid(3) 281 value = BA_alloc_get(mt_dbl,nsize,'vl',vl(2),vl(1)) 282 value = value.and.BA_alloc_get(mt_dbl,nsize, 283 > 'vlpaw',vlpaw(2),vlpaw(1)) 284 value = value.and.BA_alloc_get(mt_dbl,nsize*(nproj), 285 > 'vnl',vnl(2), vnl(1)) 286 value = value.and.BA_alloc_get(mt_dbl,3*nsize,'G',G_hndl,G_indx) 287 if (.not.value)call errquit('pawppv1:out of heap',0,MA_ERR) 288 289* **** preparation of constants **** 290 nfft1=ngrid(1) 291 nfft2=ngrid(2) 292 nfft3=ngrid(3) 293 call setup_kbpp(nfft1,nfft2,nfft3,unita,unitg,dbl_mb(G_indx)) 294 filter = control_kbpp_filter() 295 296 !**** allocate memory for rays **** 297 nray = kbpp_calc_nray(nfft1,nfft2,nfft3,unita) 298 299 value = BA_alloc_get(mt_dbl,nray, 300 > 'G_ray',G_ray_hndl,G_ray_indx) 301 value = value.and.BA_alloc_get(mt_dbl,2*nray, 302 > 'vl_ray',vl_ray_hndl,vl_ray_indx) 303 value = value.and.BA_alloc_get(mt_dbl,2*nray, 304 > 'vlpaw_ray',vlpaw_ray_hndl,vlpaw_ray_indx) 305 value = value.and.BA_alloc_get(mt_dbl,2*nray*(nbasis), 306 > 'vnl_ray',vnl_ray_hndl,vnl_ray_indx) 307 value = value.and.BA_alloc_get(mt_dbl,nray, 308 > 'tmp_ray',tmp_ray_hndl,tmp_ray_indx) 309 if (.not.value) 310 > call errquit('pawppv1:out of heap memory',0,MA_ERR) 311 312 call kbpp_generate_G_ray(nfft1,nfft2,nfft3, 313 > unita,dbl_mb(G_ray_indx)) 314 315 zcore = zion - zv 316 call integrate_pawppv1(version,rlocal, 317 > n1dgrid,log_amesh,nbasis,nmax,lmax,icut, 318 > zv,zcore,sigma, 319 > int_mb(nps(1)),int_mb(lps(1)), 320 > dbl_mb(v_ps(1)), 321 > dbl_mb(prj_ps(1)), 322 > dbl_mb(phi_ae(1)), 323 > dbl_mb(dphi_ae(1)), 324 > dbl_mb(phi_ps(1)), 325 > dbl_mb(dphi_ps(1)), 326 > dbl_mb(core_ae(1)), 327 > dbl_mb(core_ps(1)), 328 > dbl_mb(rgrid(1)), 329 > dbl_mb(f1(1)), 330 > dbl_mb(f2(1)), 331 > dbl_mb(f3(1)), 332 > dbl_mb(f4(1)), 333 > dbl_mb(cs(1)), 334 > dbl_mb(sn(1)), 335 > nfft1,nfft2,nfft3,nproj, 336 > dbl_mb(G_indx), 337 > dbl_mb(vl(1)), 338 > dbl_mb(vlpaw(1)), 339 > dbl_mb(vnl(1)), 340 > int_mb(nprj(1)), 341 > int_mb(lprj(1)), 342 > int_mb(mprj(1)), 343 > int_mb(bprj(1)), 344 > dbl_mb(Gijl(1)), 345 > dbl_mb(hartree_matrix(1)), 346 > dbl_mb(comp_charge_matrix(1)), 347 > dbl_mb(comp_pot_matrix(1)), 348 > nray, 349 > dbl_mb(G_ray_indx), 350 > dbl_mb(vl_ray_indx), 351 > dbl_mb(vlpaw_ray_indx), 352 > dbl_mb(vnl_ray_indx), 353 > dbl_mb(tmp_ray_indx), 354 > filter, 355 > ierr) 356 value = BA_free_heap(tmp_ray_hndl) 357 value = value.and.BA_free_heap(vl_ray_hndl) 358 value = value.and.BA_free_heap(vlpaw_ray_hndl) 359 value = value.and.BA_free_heap(vnl_ray_hndl) 360 value = value.and.BA_free_heap(G_ray_hndl) 361 if (.not.value) 362 > call errquit('pawppv1:Error freeing memory',0,MA_ERR) 363 364 365* **** calculate radial derivatives of core densities **** 366 call pawppv1_derivative_ngrid( 367 > n1dgrid, 368 > log_amesh, 369 > dbl_mb(rgrid(1)), 370 > dbl_mb(core_ae(1)), 371 > dbl_mb(core_ae_prime(1))) 372 call pawppv1_derivative_ngrid( 373 > n1dgrid, 374 > log_amesh, 375 > dbl_mb(rgrid(1)), 376 > dbl_mb(core_ps(1)), 377 > dbl_mb(core_ps_prime(1))) 378 379* *** integrate core density *** 380c fourpi = 16.0d0*datan(1.0d0) 381c zcore = fourpi*log_integrate_def(0, 382c > dbl_mb(core_ae(1)), 383c > 2,dbl_mb(rgrid(1)),log_amesh,n1dgrid) 384c write(luout,*) "Zcore=",zcore,rmax,dbl_mb(rgrid(1)+n1dgrid-1) 385c zcore = fourpi*log_integrate_def(0, 386c > dbl_mb(core_ps(1)), 387c > 2,dbl_mb(rgrid(1)),log_amesh,n1dgrid) 388c write(luout,*) "Zcoreps=",zcore 389 390c* **** compute the core_ion_energy = ecorez + ecorecore **** 391c ecorez = -zion*fourpi 392c > *log_integrate_def(0,dbl_mb(core_ae(1)), 393c > 1,dbl_mb(rgrid(1)), 394c > log_amesh,n1dgrid) 395 if (dabs(zcore).gt.1.0d-9) then 396 ecorez = log_coulomb0_energy(dbl_mb(core_ae(1)),zion-zv, 397 > dbl_mb(rgrid(1)),n1dgrid,log_amesh, 398 > zion) 399 ecorecore = log_coulomb_energy(dbl_mb(core_ae(1)),zion-zv, 400 > dbl_mb(rgrid(1)),n1dgrid,log_amesh) 401 core_ion_energy = ecorez + ecorecore 402 else 403 core_ion_energy = 0.0d0 404 end if 405c write(*,*) "zv,zion,zcore=",zv,zion,zcore 406c write(*,*) "core_ion_energy=",core_ion_energy,ecorez,ecorecore 407 408 409 410 if (taskid.eq.MASTER) then 411 call util_file_name_noprefix(formatted_filename, 412 > .false., 413 > .false., 414 > full_filename) 415 l = index(full_filename,' ') - 1 416 if (mprint) then 417 write(luout,*) 418 write(luout,*) "Generated formatted_filename: ",full_filename(1:l) 419c if (filter) write(luout,*) "- filtering pseudopotential -" 420 end if 421 call openfile(2,full_filename,l,'w',l) 422 423 call cwrite(2,comment,80) 424 call iwrite(2,psp_type,1) 425 call iwrite(2,version,1) 426 call iwrite(2,ngrid,3) 427 call dwrite(2,unita,9) 428 call cwrite(2,atom,2) 429 amass = 0.0d0 430 call dwrite(2,amass,1) 431 call dwrite(2,zv,1) 432 call iwrite(2,lmax,1) 433 !call iwrite(2,locp,1) 434 call iwrite(2,nbasis,1) 435 436 call iwrite(2,nmax,1) 437 call dwrite(2,rc,lmax+1) 438 439 call iwrite(2,nproj,1) 440 if (nproj.gt.0) then 441 call iwrite(2,int_mb(nprj(1)),nproj) 442 call iwrite(2,int_mb(lprj(1)),nproj) 443 call iwrite(2,int_mb(mprj(1)),nproj) 444 call iwrite(2,int_mb(bprj(1)),nproj) 445 call dwrite(2,dbl_mb(Gijl(1)),5*nmax*nmax*(lmax+1)) 446 end if 447 448 if (version.eq.4) call dwrite(2,rlocal,1) 449 rcore = 0.0d0 450 call dwrite(2,rcore,1) 451 452 453* **** other PAW matrices **** 454 l = nbasis*nbasis*nbasis*nbasis*(2*lmax+1) 455 call dwrite(2,dbl_mb(hartree_matrix(1)),l) 456 l = nbasis*nbasis*(2*lmax+1) 457 call dwrite(2,dbl_mb(comp_charge_matrix(1)),l) 458 call dwrite(2,dbl_mb(comp_pot_matrix(1)),l) 459 460* **** miscelaneous PAW energies **** 461 call dwrite(2,core_kin_energy,1) 462 call dwrite(2,core_ion_energy,1) 463 464* **** write 1d-wavefunctions **** 465 call iwrite(2,n1dgrid,1) 466 call iwrite(2,icut,1) 467 call dwrite(2,log_amesh,1) 468 call dwrite(2,r1,1) 469 call dwrite(2,rmax,1) 470 call dwrite(2,sigma,1) 471 call dwrite(2,zion,1) 472 call dwrite(2,dbl_mb(eig(1)),nbasis) 473 call iwrite(2,int_mb(nae(1)),nbasis) 474 call iwrite(2,int_mb(nps(1)),nbasis) 475 call iwrite(2,int_mb(lps(1)),nbasis) 476 477 call dwrite(2,dbl_mb(rgrid(1)),n1dgrid) 478 call dwrite(2,dbl_mb(phi_ae(1)),n1dgrid*nbasis) 479 call dwrite(2,dbl_mb(dphi_ae(1)),n1dgrid*nbasis) 480 call dwrite(2,dbl_mb(phi_ps(1)),n1dgrid*nbasis) 481 call dwrite(2,dbl_mb(dphi_ps(1)),n1dgrid*nbasis) 482 call dwrite(2,dbl_mb(core_ae(1)),n1dgrid) 483 call dwrite(2,dbl_mb(core_ps(1)),n1dgrid) 484 call dwrite(2,dbl_mb(core_ae_prime(1)),n1dgrid) 485 call dwrite(2,dbl_mb(core_ps_prime(1)),n1dgrid) 486 487 call dwrite(2,dbl_mb(vl(1)),nsize) 488 call dwrite(2,dbl_mb(vlpaw(1)),nsize) 489 if (nproj.gt.0) then 490 call dwrite(2,dbl_mb(vnl(1)),nsize*nproj) 491 end if 492 493 call closefile(2) 494 end if 495 496 497 value = .true. 498 value = value.and.BA_free_heap(G_hndl) 499 value = value.and.BA_free_heap(vnl(2)) 500 value = value.and.BA_free_heap(vlpaw(2)) 501 value = value.and.BA_free_heap(vl(2)) 502 value = value.and.BA_free_heap(sn(2)) 503 value = value.and.BA_free_heap(cs(2)) 504 value = value.and.BA_free_heap(f1(2)) 505 value = value.and.BA_free_heap(f2(2)) 506 value = value.and.BA_free_heap(f3(2)) 507 value = value.and.BA_free_heap(f4(2)) 508 509 value = value.and.BA_free_heap(Gijl(2)) 510 value = value.and.BA_free_heap(hartree_matrix(2)) 511 value = value.and.BA_free_heap(comp_charge_matrix(2)) 512 value = value.and.BA_free_heap(comp_pot_matrix(2)) 513 514 value = value.and.BA_free_heap(bprj(2)) 515 value = value.and.BA_free_heap(mprj(2)) 516 value = value.and.BA_free_heap(lprj(2)) 517 value = value.and.BA_free_heap(nprj(2)) 518 519 value = value.and.BA_free_heap(prj_ps0(2)) 520 value = value.and.BA_free_heap(prj_ps(2)) 521 value = value.and.BA_free_heap(v_ps(2)) 522 value = value.and.BA_free_heap(core_ps_prime(2)) 523 value = value.and.BA_free_heap(core_ae_prime(2)) 524 value = value.and.BA_free_heap(core_ps(2)) 525 value = value.and.BA_free_heap(core_ae(2)) 526 value = value.and.BA_free_heap(dphi_ps(2)) 527 value = value.and.BA_free_heap(phi_ps(2)) 528 value = value.and.BA_free_heap(dphi_ae(2)) 529 value = value.and.BA_free_heap(phi_ae(2)) 530 531 value = value.and.BA_free_heap(eig(2)) 532 value = value.and.BA_free_heap(lps(2)) 533 value = value.and.BA_free_heap(nps(2)) 534 value = value.and.BA_free_heap(nae(2)) 535 536 value = value.and.BA_free_heap(rgrid(2)) 537 if (.not.value) call errquit('pawppv1:freeing heap',5,MA_ERR) 538 539 pawppv1 = value 540 return 541 542 9999 call errquit('pawppv1:Error reading psp_filename',0,DISK_ERR) 543 pawppv1 = value 544 return 545 546 END 547 548 549 550 551c ************************************************* 552c * * 553c * pawppv1_derivative_ngrid * 554c * * 555c ************************************************* 556c 557c This routine computes the seven point derivative of f. 558c where f and df are stored on a logarithmic grid. The 559c dimensions of f and df are, f(1:ng), and df(1:ng) 560 561 subroutine pawppv1_derivative_ngrid(ng,log_amesh,r,f,df) 562 implicit none 563 integer ng 564 double precision log_amesh 565 double precision r(ng) 566 double precision f(ng) 567 double precision df(ng) 568 569 double precision one_over_60 570 parameter (one_over_60 = 1.0d0/60.0d0) 571 572 integer i,n1,n2,m1,m2 573 double precision aa 574 575 aa = one_over_60/log_amesh 576 n1 = 1 577 n2 = ng 578 m1 = n1 579 m2 = n2 580 581 582 if (n1.le.3) then 583 if ((n1.eq.1).and.(n1.ge.m1).and.(n1.le.m2)) then 584 df(1) = aa*(-147.0d0*f(1) 585 > + 360.0d0*f(2) 586 > - 450.0d0*f(3) 587 > + 400.0d0*f(4) 588 > - 225.0d0*f(5) 589 > + 72.0d0*f(6) 590 > - 10.0d0*f(7))/r(1) 591 n1 = n1+1 592 end if 593 if ((n1.eq.2).and.(n1.ge.m1).and.(n1.le.m2)) then 594 df(2) = aa*( -10.0d0*f(1) 595 > - 77.0d0*f(2) 596 > + 150.0d0*f(3) 597 > - 100.0d0*f(4) 598 > + 50.0d0*f(5) 599 > - 15.0d0*f(6) 600 > + 2.0d0*f(7))/r(2) 601 n1 = n1+1 602 end if 603 if ((n1.eq.3.and.(n1.ge.m1).and.(n1.le.m2))) then 604 df(3) = aa*( +2.0d0*f(1) 605 > - 24.0d0*f(2) 606 > - 35.0d0*f(3) 607 > + 80.0d0*f(4) 608 > - 30.0d0*f(5) 609 > + 8.0d0*f(6) 610 > - 1.0d0*f(7))/r(3) 611 n1 = n1+1 612 end if 613 end if 614 615 if (n2.ge.(ng-2)) then 616 if ((n2.eq.ng).and.(n2.ge.m1).and.(n2.le.m2)) then 617 df(ng) = aa*( +147.0d0*f(ng) 618 > - 360.0d0*f(ng-1) 619 > + 450.0d0*f(ng-2) 620 > - 400.0d0*f(ng-3) 621 > + 225.0d0*f(ng-4) 622 > - 72.0d0*f(ng-5) 623 > + 10.0d0*f(ng-6))/r(ng) 624 n2 = n2-1 625 end if 626 if ((n2.eq.(ng-1).and.(n2.ge.m1).and.(n2.le.m2))) then 627 df(ng-1) = aa*( +10.0d0*f(ng) 628 > + 77.0d0*f(ng-1) 629 > - 150.0d0*f(ng-2) 630 > + 100.0d0*f(ng-3) 631 > - 50.0d0*f(ng-4) 632 > + 15.0d0*f(ng-5) 633 > - 2.0d0*f(ng-6))/r(ng-1) 634 n2 = n2-1 635 end if 636 if ((n2.eq.(ng-2).and.(n2.ge.m1).and.(n2.le.m2))) then 637 df(ng-2) = aa*( -2.0d0*f(ng) 638 > + 24.0d0*f(ng-1) 639 > + 35.0d0*f(ng-2) 640 > - 80.0d0*f(ng-3) 641 > + 30.0d0*f(ng-4) 642 > - 8.0d0*f(ng-5) 643 > + 1.0d0*f(ng-6))/r(ng-2) 644 n2 = n2-1 645 end if 646 end if 647 648 do i=n1,n2 649 df(i) = aa*( -1.0d0*f(i-3) 650 > + 9.0d0*f(i-2) 651 > - 45.0d0*f(i-1) 652 > + 45.0d0*f(i+1) 653 > - 9.0d0*f(i+2) 654 > + 1.0d0*f(i+3))/r(i) 655 end do 656 657 return 658 end 659 660 661 662 663 664