1* 2* $Id$ 3* 4 5 6#define TCGMSG 7 8* *********************************** 9* * * 10* * cpsp_stress_init * 11* * * 12* *********************************** 13 14 subroutine cpsp_stress_init() 15 implicit none 16 17#include "bafdecls.fh" 18#include "cpsp_common.fh" 19#include "errquit.fh" 20 21 integer npack1,npack0,nbrillq 22 logical value 23 24* **** external functions ***** 25 integer ion_nkatm,brillioun_nbrillq 26 external ion_nkatm,brillioun_nbrillq 27 28 29 call Cram_npack(0,npack0) 30 call Cram_max_npack(npack1) 31 nbrillq = brillioun_nbrillq() 32 33 value = BA_alloc_get(mt_dbl,(npsp*npack0), 34 > 'dvl',dvl(2),dvl(1)) 35 value = value.and. 36 > BA_alloc_get(mt_int,npsp,'dvnl',dvnl(2),dvnl(1)) 37 38 if (.not. value) 39 > call errquit('cpsp_stress_init:out of heap memory',0,MA_ERR) 40 41 call dcopy(npsp*npack0,0.0d0,0,dbl_mb(dvl(1)),1) 42 return 43 end 44 45 46 47* *********************************** 48* * * 49* * cpsp_stress_end * 50* * * 51* *********************************** 52 53 subroutine cpsp_stress_end() 54 implicit none 55 56#include "errquit.fh" 57#include "bafdecls.fh" 58#include "cpsp_common.fh" 59 60 61 logical value 62 63 value = BA_free_heap(dvl(2)) 64 value = value.and.BA_free_heap(dvnl(2)) 65 66 if (.not.value) 67 > call errquit('cpsp_stress_end:freeing heap',0, MA_ERR) 68 69 return 70 end 71 72 73* *********************************** 74* * * 75* * cpsp_v_local_euv * 76* * * 77* *********************************** 78 79 subroutine cpsp_v_local_euv(dng,euv) 80 implicit none 81 complex*16 dng(*) 82 real*8 euv(3,3) 83 84#include "bafdecls.fh" 85#include "cpsp_common.fh" 86#include "errquit.fh" 87 88* *** local variables *** 89 integer nfft3d,npack0 90 integer ii,ia,u,v,s 91 integer exi(2),tmp1(2),tmp2(2) 92 integer G(2,3),vll(2) 93 logical value 94 real*8 elocal,ftmp(3) 95 real*8 hm(3,3),Bus(3,3) 96 real*8 ss,sum,pi,fourpi 97 98* **** common block used for coulomb.f **** 99 integer vc_indx,vc_hndl 100 common / c_vc_block / vc_indx,vc_hndl 101 102 103* **** external functions **** 104 integer c_G_indx,ion_nion,ion_katm 105 real*8 lattice_omega,lattice_unitg 106 external c_G_indx,ion_nion,ion_katm 107 external lattice_omega,lattice_unitg 108 109 call nwpw_timing_start(5) 110 111 call C3dB_nfft3d(1,nfft3d) 112 call Cram_npack(0,npack0) 113 114 pi = 4.0d0*datan(1.0d0) 115 fourpi = 4.0d0*pi 116 ss = 1.0d0/(2.0d0*pi) 117 118* *** define hm **** 119 do v=1,3 120 do u=1,3 121 hm(u,v) = ss*lattice_unitg(u,v) 122 end do 123 end do 124 125* **** average Kohn-Sham v_local energy **** 126 value = BA_push_get(mt_dcpl,npack0,'vll',vll(2),vll(1)) 127 if (.not. value) 128 > call errquit('cpsp_v_local_euv:out of stack memory',0,MA_ERR) 129 call cpsp_v_local(dcpl_mb(vll(1)),.false.,dng,ftmp) 130 call Cram_cc_dot(0,dng,dcpl_mb(vll(1)),elocal) 131 value = BA_pop_stack(vll(2)) 132 if (.not. value) 133 > call errquit('cpsp_v_local_euv:error popping stack',0,MA_ERR) 134 135 136 value = BA_push_get(mt_dcpl,nfft3d,'exi', exi(2), exi(1)) 137 value = value.and. 138 > BA_push_get(mt_dbl, npack0,'tmp1',tmp1(2),tmp1(1)) 139 value = value.and. 140 > BA_push_get(mt_dbl, npack0,'tmp2',tmp2(2),tmp2(1)) 141 value = value.and. 142 > BA_push_get(mt_dbl, nfft3d,'G1',G(2,1),G(1,1)) 143 value = value.and. 144 > BA_push_get(mt_dbl, nfft3d,'G2',G(2,2),G(1,2)) 145 value = value.and. 146 > BA_push_get(mt_dbl, nfft3d,'G3',G(2,3),G(1,3)) 147 if (.not. value) 148 > call errquit('cpsp_v_local_euv:out of stack memory',0,MA_ERR) 149 150* **** define Gx,Gy and Gz in packed space **** 151 call C3dB_t_Copy(1,dbl_mb(c_G_indx(1)),dbl_mb(G(1,1))) 152 call C3dB_t_Copy(1,dbl_mb(c_G_indx(2)),dbl_mb(G(1,2))) 153 call C3dB_t_Copy(1,dbl_mb(c_G_indx(3)),dbl_mb(G(1,3))) 154 call Cram_r_pack(0,dbl_mb(G(1,1))) 155 call Cram_r_pack(0,dbl_mb(G(1,2))) 156 call Cram_r_pack(0,dbl_mb(G(1,3))) 157 158 159 call dcopy(9,0.0d0,0,Bus,1) 160 do ii=1,ion_nion() 161 ia=ion_katm(ii) 162 163* **** structure factor and local pseudopotential **** 164 call cstrfac_pack(0,ii,dcpl_mb(exi(1))) 165 166* **** tmp2(G) = Real[ dconjg(dng(G))*exi(G) ] **** 167 call Cram_ccr_conjgMul(0,dng, 168 > dcpl_mb(exi(1)), 169 > dbl_mb(tmp2(1))) 170 171* **** tmp2(G) = tmp2(G)*(dvl(G)) 172c call Cram_rr_Mul(0,dbl_mb(tmp2(1)), 173c > dbl_mb(dvl(1)+(ia-1)*npack0), 174c > dbl_mb(tmp2(1))) 175 call Cram_rr_Mul2(0,dbl_mb(dvl(1)+(ia-1)*npack0), 176 > dbl_mb(tmp2(1))) 177 178* **** tmp2(G) = tmp2(G)/G **** 179 ss = 1.0d0/fourpi 180 call Cram_r_SMul(0,ss,dbl_mb(vc_indx),dbl_mb(tmp1(1))) 181c call Cram_rr_Sqrt(0,dbl_mb(tmp1(1)),dbl_mb(tmp1(1))) 182c call Cram_rr_Mul(0,dbl_mb(tmp1(1)), 183c > dbl_mb(tmp2(1)), 184c > dbl_mb(tmp2(1))) 185 call Cram_rr_Sqrt1(0,dbl_mb(tmp1(1))) 186 call Cram_rr_Mul2(0,dbl_mb(tmp1(1)),dbl_mb(tmp2(1))) 187 188* **** Bus = Bus - Sum(G) tmp2(G)*Gu*Gs *** 189 do u=1,3 190 do s=u,3 191 call Cram_rr_Mul(0,dbl_mb(G(1,u)), 192 > dbl_mb(G(1,s)), 193 > dbl_mb(tmp1(1))) 194 call Cram_rr_dot(0,dbl_mb(tmp1(1)),dbl_mb(tmp2(1)),sum) 195 196 Bus(u,s) = Bus(u,s) - sum 197 end do 198 end do 199 200 end do 201 value = BA_pop_stack(G(2,3)) 202 value = value.and.BA_pop_stack(G(2,2)) 203 value = value.and.BA_pop_stack(G(2,1)) 204 value = value.and.BA_pop_stack(tmp2(2)) 205 value = value.and.BA_pop_stack(tmp1(2)) 206 value = value.and.BA_pop_stack(exi(2)) 207 if (.not. value) call errquit('error popping stack memory',0, 208 & MA_ERR) 209 210 do u=1,3 211 do s=u+1,3 212 Bus(s,u) = Bus(u,s) 213 end do 214 end do 215 do v=1,3 216 do u=1,3 217 euv(u,v) = -elocal*hm(u,v) 218 do s=1,3 219 euv(u,v) = euv(u,v) + Bus(u,s)*hm(s,v) 220 end do 221 end do 222 end do 223 224 225 call nwpw_timing_end(5) 226 return 227 end 228 229 230* *********************************** 231* * * 232* * cpsp_v_nonlocal_euv_2 * 233* * * 234* *********************************** 235 236 subroutine cpsp_v_nonlocal_euv_2(ispin,ne,psi1_tag,euv) 237 implicit none 238 integer ispin,ne(2) 239 integer psi1_tag 240 real*8 euv(3,3) 241 242#include "bafdecls.fh" 243#include "cpsp_common.fh" 244#include "errquit.fh" 245 246 247* *** local variables *** 248 integer nfft3d,npack1,npack,shift,shift2,nbrillq,nb 249 integer nproj,l_prj,psi1_shift,occ1_tag,occ1_shift,occ_shift 250 integer i,ii,ia,k,l,n,nn 251 integer s,u,v 252 real*8 omega,Bus(3,3),hm(3,3),kx,ky,kz,Aus(3,3) 253 complex*16 ctmp,cxr 254 real*8 pi,scal,weight 255 integer exi(2),vtmp(2),tmp1(2),zsw1(2),zsw2(2),zsw3(2) 256 integer G(2,3) 257 logical value,sd_function 258 259* **** external functions **** 260 logical is_sORd 261 integer ion_nion,ion_katm,c_G_indx 262 integer Pneb_nbrillq,cpsp_projector_get_ptr 263 integer cpsi_data_get_chnk,cpsi_data_get_next 264 real*8 brillioun_weight,brillioun_k 265 real*8 lattice_omega,lattice_unitg 266 external is_sORd 267 external ion_nion,ion_katm,c_G_indx 268 external Pneb_nbrillq,cpsp_projector_get_ptr 269 external cpsi_data_get_chnk,cpsi_data_get_next 270 external brillioun_weight,brillioun_k 271 external lattice_omega,lattice_unitg 272 273 call nwpw_timing_start(6) 274 275 occ1_tag = cpsi_data_get_next(psi1_tag) 276 277 omega = lattice_omega() 278 279* **** allocate local memory **** 280 nn = ne(1)+ne(2) 281 nbrillq = Pneb_nbrillq() 282 call C3dB_nfft3d(1,nfft3d) 283 call Cram_max_npack(npack1) 284 285 value = BA_push_get(mt_dcpl,nfft3d,'exi', exi(2), exi(1)) 286 value = value.and. 287 > BA_push_get(mt_dcpl,nfft3d,'vtmp',vtmp(2),vtmp(1)) 288 value = value.and. 289 > BA_push_get(mt_dbl, npack1,'tmp1',tmp1(2),tmp1(1)) 290 value = value.and. 291 > BA_push_get(mt_dbl, nfft3d,'Gx',G(2,1),G(1,1)) 292 value = value.and. 293 > BA_push_get(mt_dbl, nfft3d,'Gy',G(2,2),G(1,2)) 294 value = value.and. 295 > BA_push_get(mt_dbl, nfft3d,'Gz',G(2,3),G(1,3)) 296 value = value.and. 297 > BA_push_get(mt_dcpl,nn*nprj_max, 298 > 'zsw1',zsw1(2),zsw1(1)) 299 value = value.and. 300 > BA_push_get(mt_dcpl,nn*nprj_max, 301 > 'zsw2',zsw2(2),zsw2(1)) 302 value = value.and. 303 > BA_push_get(mt_dcpl,9*nn, 304 > 'zsw3',zsw3(2),zsw3(1)) 305 if (.not. value) 306 > call errquit('cpsp_v_nonlocal_euv_2:out of stack',0,MA_ERR) 307 308 309 310* *********************** 311* **** calculate Bus **** 312* *********************** 313 call dcopy(9,0.0d0,0,Bus,1) 314 do ii=1,ion_nion() 315 ia=ion_katm(ii) 316 nproj = int_mb(nprj(1)+ia-1) 317 318 if (nproj.gt.0) then 319 320 do nb=1,nbrillq 321 call dcopy(9,0.0d0,0,Aus,1) 322 323 psi1_shift = cpsi_data_get_chnk(psi1_tag,nb) 324 if (occ1_tag.gt.0) 325 > occ1_shift = cpsi_data_get_chnk(occ1_tag,nb) 326 call Cram_npack(nb,npack) 327 kx = brillioun_k(1,nb) 328 ky = brillioun_k(2,nb) 329 kz = brillioun_k(3,nb) 330 weight = brillioun_weight(nb) 331 call C3dB_t_Copy(1,dbl_mb(c_G_indx(1)),dbl_mb(G(1,1))) 332 call C3dB_t_Copy(1,dbl_mb(c_G_indx(2)),dbl_mb(G(1,2))) 333 call C3dB_t_Copy(1,dbl_mb(c_G_indx(3)),dbl_mb(G(1,3))) 334 call Cram_r_pack(nb,dbl_mb(G(1,1))) 335 call Cram_r_pack(nb,dbl_mb(G(1,2))) 336 call Cram_r_pack(nb,dbl_mb(G(1,3))) 337 call daxpy(npack,1.0d0,kx,0,dbl_mb(G(1,1)),1) 338 call daxpy(npack,1.0d0,ky,0,dbl_mb(G(1,2)),1) 339 call daxpy(npack,1.0d0,kz,0,dbl_mb(G(1,3)),1) 340 341* **** structure factor and local pseudopotential **** 342 call cstrfac_pack(nb,ii,dcpl_mb(exi(1))) 343 call cstrfac_k(ii,nb,cxr) 344 call zscal(npack,cxr,dcpl_mb(exi(1)),1) 345 346 347* ********************************************* 348* **** calculate F^(lm)_I = <psi|vnl(nlm)> **** 349* ********************************************* 350 do l=1,nproj 351 shift = cpsp_projector_get_ptr(int_mb(vnl(1)+ia-1),nb,l) 352 l_prj = int_mb(l_projector(1)+(l-1)+(ia-1)*jmmax_max) 353cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 354ccccccccccc change this for actinides where we might have g's ccccc 355cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 356#ifdef GCC4 357 k = iand(l_prj,1) 358#else 359 k = and(l_prj,1) 360#endif 361 sd_function = (k.eq.0) 362 363 364* *** current function is s or d **** 365 if (sd_function) then 366 call Cram_rc_Mul(nb,dbl_mb(shift), 367 > dcpl_mb(exi(1)), 368 > dcpl_mb(vtmp(1))) 369* *** current function is p or f **** 370 else 371 call Cram_irc_Mul(nb,dbl_mb(shift), 372 > dcpl_mb(exi(1)), 373 > dcpl_mb(vtmp(1))) 374 end if 375 call Cram_cc_inzdot(nb,nn, 376 > dbl_mb(psi1_shift), 377 > dcpl_mb(vtmp(1)), 378 > dcpl_mb(zsw1(1)+(l-1)*nn)) 379 end do 380 call C3dB_Vector_SumAll((2*nn*nproj),dcpl_mb(zsw1(1))) 381 382* **** zsw2 = Gijl*zsw1 ****** 383 call Multiply_Gijl_zsw1(nn, 384 > nproj, 385 > int_mb(nmax(1)+ia-1), 386 > int_mb(lmax(1)+ia-1), 387 > int_mb(n_projector(1)+(ia-1)*jmmax_max), 388 > int_mb(l_projector(1)+(ia-1)*jmmax_max), 389 > int_mb(m_projector(1)+(ia-1)*jmmax_max), 390 > dbl_mb(Gijl(1)+(ia-1)*gij_stride), 391 > dcpl_mb(zsw1(1)), 392 > dcpl_mb(zsw2(1))) 393 394 if (ispin.eq.1) call dscal(2*nn*nproj,2.0d0,dcpl_mb(zsw2(1)),1) 395 do n=1,nn*nproj 396 dcpl_mb(zsw1(1)+n-1) = dconjg(dcpl_mb(zsw2(1)+n-1)) 397 end do 398 399 400* ********************************** 401* **** calculate dF^(lm)_I/dhus **** 402* ********************************** 403 do l=1,nproj 404 l_prj = int_mb(l_projector(1)+(l-1)+(ia-1)*jmmax_max) 405#ifdef GCC4 406 k = iand(l_prj,1) 407#else 408 k = and(l_prj,1) 409#endif 410 sd_function = (k.eq.0) 411 412 do s=1,3 413 shift2 = cpsp_projector_get_ptr( 414 > int_mb(dvnl(1)+ia-1),nb,(3*(l-1)+s)) 415 do u=1,3 416 call Cram_rr_Mul(nb,dbl_mb(shift2), 417 > dbl_mb(G(1,u)), 418 > dbl_mb(tmp1(1))) 419 420* *** current function is s or d **** 421 if (sd_function) then 422 call Cram_rc_Mul(nb,dbl_mb(tmp1(1)), 423 > dcpl_mb(exi(1)), 424 > dcpl_mb(vtmp(1))) 425 426* *** current function is p or f **** 427 else 428 call Cram_irc_Mul(nb,dbl_mb(tmp1(1)), 429 > dcpl_mb(exi(1)), 430 > dcpl_mb(vtmp(1))) 431 end if 432 call Cram_cc_nzdot(nb,nn, 433 > dbl_mb(psi1_shift), 434 > dcpl_mb(vtmp(1)), 435 > dcpl_mb(zsw3(1)+(u-1)*nn 436 > +(s-1)*nn*3)) 437 end do 438 end do 439 440 if (occ1_tag.gt.0) then 441 occ_shift = occ1_shift 442 do i=1,nn 443 do s=1,3 444 do u=1,3 445 ctmp = dcpl_mb(zsw1(1)+(i-1)+(l-1)*nn) 446 > *dcpl_mb(zsw3(1)+(i-1) 447 > +(u-1)*nn 448 > +(s-1)*nn*3) 449 450 Bus(u,s) = Bus(u,s) 451 > - dbl_mb(occ_shift)*weight*2.0d0*dble(ctmp)/(omega) 452 end do 453 end do 454 occ_shift = occ_shift + 1 455 end do 456 else 457 do i=1,nn 458 do s=1,3 459 do u=1,3 460 ctmp = dcpl_mb(zsw1(1)+(i-1)+(l-1)*nn) 461 > *dcpl_mb(zsw3(1)+(i-1) 462 > +(u-1)*nn 463 > +(s-1)*nn*3) 464 465 Bus(u,s) = Bus(u,s) 466 > - weight*2.0d0*dble(ctmp)/(omega) 467 Aus(u,s) = 468 > - weight*2.0d0*dble(ctmp)/(omega) 469 end do 470 end do 471 end do 472 end if 473 474 end do !** l ** 475 end do !** nb ** 476 477 end if 478 end do !** ii ** 479 call K1dB_Vector_SumAll(9,Bus) 480 481 value = BA_pop_stack(zsw3(2)) 482 value = value.and.BA_pop_stack(zsw2(2)) 483 value = value.and.BA_pop_stack(zsw1(2)) 484 value = value.and.BA_pop_stack(G(2,3)) 485 value = value.and.BA_pop_stack(G(2,2)) 486 value = value.and.BA_pop_stack(G(2,1)) 487 value = value.and.BA_pop_stack(tmp1(2)) 488 value = value.and.BA_pop_stack(vtmp(2)) 489 value = value.and.BA_pop_stack(exi(2)) 490 if (.not. value) 491 > call errquit('cpsp_v_nonlocal_euv_2:error popping stack',0, 492 & MA_ERR) 493 494 495* *** define hm **** 496 pi = 4.0d0*datan(1.0d0) 497 scal = 1.0d0/(2.0d0*pi) 498 do v=1,3 499 do u=1,3 500 hm(u,v) = scal*lattice_unitg(u,v) 501 end do 502 end do 503 504* *** calculate euv = Sum(s) hm(s,v)*Bus(u,s) 505 call dcopy(9,0.0d0,0,euv,1) 506 do u=1,3 507 do v=1,3 508 do s=1,3 509 euv(u,v) = euv(u,v) + Bus(u,s)*hm(s,v) 510 end do 511 end do 512 end do 513 514 515 call nwpw_timing_end(6) 516 return 517 end 518 519 520 521* *********************************** 522* * * 523* * cpsp_stress_read * 524* * * 525* *********************************** 526 527 subroutine cpsp_stress_read(fname, 528 > version, 529 > nfft,unita, 530 > npack0,dvl, 531 > npack1,lmmax,lmmax_max,dvnl_tag, 532 > semicore,dncore, 533 > tmp,tmp2, 534 > ierr) 535 implicit none 536 character*50 fname 537 integer version 538 integer nfft(3) 539 real*8 unita(3,3) 540 integer npack0 541 real*8 dvl(*) 542 integer npack1,lmmax,lmmax_max 543 integer dvnl_tag 544 logical semicore 545 real*8 dncore(*) 546 real*8 tmp(*) 547 real*8 tmp2(*) 548 integer ierr 549 550#ifdef MPI 551 include 'mpif.h' 552 integer mpierr 553#endif 554#ifdef TCGMSG 555#include "tcgmsg.fh" 556#include "msgtypesf.h" 557#endif 558 559* *** local variables *** 560 integer MASTER,taskid,taskid_k 561 parameter(MASTER=0) 562 integer i,n,l 563 integer msglen 564 character*255 full_filename 565 566 real*8 kv(3) 567 integer nbrillioun,nb,nbq,pk 568 569 integer brillioun_nbrillioun,brillioun_nbrillq 570 integer cpsp_projector_alloc 571 real*8 brillioun_all_k 572 external brillioun_nbrillioun,brillioun_nbrillq 573 external cpsp_projector_alloc 574 external brillioun_all_k 575 576 call Parallel_taskid(taskid) 577 call Parallel3d_taskid_k(taskid_k) 578 579* **** open fname binary file **** 580 if (taskid.eq.MASTER) then 581 call util_file_name_noprefix(fname,.false., 582 > .false., 583 > full_filename) 584 l = index(full_filename,' ') - 1 585 call openfile(5,full_filename,l,'r',l) 586 call iread(5,version,1) 587 call iread(5,nfft,3) 588 call dread(5,unita,9) 589 call iread(5,nbrillioun,1) 590 ierr = 0 591 if (nbrillioun.eq.brillioun_nbrillioun()) then 592 do nb=1,nbrillioun 593 call dread(5,kv,3) 594 if ((brillioun_all_k(1,nb).ne.kv(1)).or. 595 > (brillioun_all_k(2,nb).ne.kv(2)).or. 596 > (brillioun_all_k(3,nb).ne.kv(3))) ierr = 1 597 end do 598 else 599 ierr = 1 600 call closefile(5) 601 end if 602 end if 603 604 msglen = 1 605 call Parallel_Brdcst_ivalues(MASTER,msglen,ierr) 606 if (ierr.ne.0) then 607 return 608 end if 609 610 611* **** send header data to all processors **** 612 msglen = 1 613 call Parallel_Brdcst_ivalues(MASTER,msglen,version) 614 msglen = 3 615 call Parallel_Brdcst_ivalues(MASTER,msglen,nfft) 616 msglen = 9 617 call Parallel_Brdcst_values(MASTER,msglen,unita) 618 619 620* *** read in dvl 3d block *** 621 call C3dB_r_read(1,5,tmp2,tmp,-1,-1,.true.) 622 call Cram_r_pack(0,tmp2) 623 call Cram_r_Copy(0,tmp2,dvl) 624 625* **** read in dvnl 3d blocks **** 626 if (lmmax.gt.0) then 627 dvnl_tag = cpsp_projector_alloc(brillioun_nbrillq(), 628 > 3*lmmax,npack1) 629 do nb=1,brillioun_nbrillioun() 630 call K1dB_ktoqp(nb,nbq,pk) 631 do n=1,lmmax 632 do i=1,3 633 call C3dB_r_read(1,5,tmp2,tmp,-1,pk,.true.) 634 if (pk.eq.taskid_k) then 635 call Cram_r_pack(nbq,tmp2) 636 call cpsp_projector_add(dvnl_tag,nbq,(3*(n-1)+i),tmp2) 637 end if 638 end do 639 end do 640 end do 641 end if 642 643* **** read in semicore density block **** 644 if (semicore) then 645 !write(*,*) "reading in semicore block" !debug 646 call C3dB_r_Read(1,5,tmp2,tmp,-1,-1,.true.) 647 call Cram_r_pack(0,tmp2) 648 call Cram_r_Copy(0,tmp2,dncore) 649 end if 650 651* *** close fname binary file *** 652 if (taskid.eq.MASTER) then 653c close(11) 654 call closefile(5) 655 end if 656 657 ierr = 0 658 return 659 end 660 661* *********************************** 662* * * 663* * cpsp_stress_readall * 664* * * 665* *********************************** 666 667 subroutine cpsp_stress_readall() 668 implicit none 669 670#include "bafdecls.fh" 671#include "errquit.fh" 672#include "cpsp_common.fh" 673#include "c_semicore_common.fh" 674 675 676 677* **** local variables **** 678 integer ngp(3) 679 real*8 unita(3,3) 680 integer version,nfft3d,npack1,npack0,nbrill 681 integer ia,l,nproj 682 character*12 boundry 683 integer tmp(2),tmp2(2),ierr 684 logical value,found,correct_box 685 character*5 element 686 character*50 fname 687 688* **** parallel i/o variable **** 689 integer MASTER,taskid 690 parameter(MASTER=0) 691 692* **** external functions **** 693 logical nwpw_filefind 694 integer control_ngrid,brillioun_nbrillioun 695 integer psp_lmmax,cpsp_nprj 696 real*8 control_unita 697 character*12 control_boundry 698 character*4 ion_atom 699 external nwpw_filefind 700 external control_ngrid,brillioun_nbrillioun 701 external psp_lmmax 702 external control_unita 703 external control_boundry 704 external ion_atom 705 external cpsp_nprj 706 707 call Parallel_taskid(taskid) 708 709 call C3dB_nfft3d(1,nfft3d) 710 call Cram_npack(0,npack0) 711 call Cram_max_npack(npack1) 712 nbrill = brillioun_nbrillioun() 713 714ccc corrected to complex by pjn 11-7-06 715 value = BA_push_get(mt_dbl,(4*nfft3d),'tmp',tmp(2),tmp(1)) 716 if (.not. value) 717 > call errquit('cpsp_stress_readall:out of stack memory',0,MA_ERR) 718 719 value = BA_push_get(mt_dbl,(4*nfft3d),'tmp2',tmp2(2),tmp2(1)) 720 if (.not. value) 721 > call errquit('cpsp_stress_readall:out of stack memory',0,MA_ERR) 722 723* **** read pseudopotentials **** 724 do ia=1,npsp 725 726* **** define formatted psp name **** 727 element = ' ' 728 element = ion_atom(ia) 729 l = index(element,' ') - 1 730 fname = element(1:l)//'.cpp2' 731ccccccccccccc possible error here!!!!!! CCCCCCCCCCCCCCCCCCc 732cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 733c !nproj = psp_lmmax(ia) 734 nproj=cpsp_nprj(ia) 735 736 737* **** not finished **** 738 found = .false. 739 do while (.not. found) 740 if (nwpw_filefind(fname)) then 741 call cpsp_stress_read(fname, 742 > version, 743 > ngp,unita, 744 > npack0, 745 > dbl_mb(dvl(1) + (ia-1)*npack0), 746 > npack1,nproj,lmmax_max, 747 > int_mb(dvnl(1)+ia-1), 748 > log_mb(semicore(1)+ia), 749 > dbl_mb(ncore(1) + npack0 + (ia-1)*5*npack0), 750 > dbl_mb(tmp(1)),dbl_mb(tmp2(1)), 751 > ierr) 752 753 754* **** set semicore(0) **** 755 if (log_mb(semicore(1)+ia)) log_mb(semicore(1)) = .true. 756 if (ierr.gt.1) go to 9000 757 758* ************************************************************** 759* ***** logic for finding out if psp is correctly formatted **** 760* ************************************************************** 761 correct_box = .true. 762 if ( (ngp(1).ne.control_ngrid(1)) .or. 763 > (ngp(2).ne.control_ngrid(2)) .or. 764 > (ngp(3).ne.control_ngrid(3)) .or. 765 > (unita(1,1).ne.control_unita(1,1)) .or. 766 > (unita(2,1).ne.control_unita(2,1)) .or. 767 > (unita(3,1).ne.control_unita(3,1)) .or. 768 > (unita(1,2).ne.control_unita(1,2)) .or. 769 > (unita(2,2).ne.control_unita(2,2)) .or. 770 > (unita(3,2).ne.control_unita(3,2)) .or. 771 > (unita(1,3).ne.control_unita(1,3)) .or. 772 > (unita(2,3).ne.control_unita(2,3)) .or. 773 > (unita(3,3).ne.control_unita(3,3)) .or. 774 > ((boundry(1:l).eq.'periodic').and.(version.ne.3)).or. 775 > ((boundry(1:l).eq.'aperiodic').and.(version.ne.4))) then 776 correct_box = .false. 777 if (taskid.eq.MASTER) then 778 write(6,*) "pseudopotential is not correctly formatted:", 779 > fname 780 end if 781 if ((ierr.eq.0).and.(nproj.gt.0)) then 782 call cpsp_projector_dealloc(int_mb(dvnl(1)+ia-1)) 783 end if 784 end if 785 if (correct_box) found = .true. 786 if (ierr.eq.1) then 787 found = .false. 788 if (taskid.eq.MASTER) then 789 write(6,*) 790 write(6,*) "pseudopotential is not correctly formatted-", 791 > "bad brillioun zone:",fname 792 end if 793 end if 794 795 end if 796 797* **** generate formatted pseudopotential atom.cpp2 ***** 798 if (.not.found) then 799 call cpsp_stress_formatter_auto(ion_atom(ia)) 800 end if 801 802 end do !*** do while **** 803 804 end do 805 9000 value = BA_pop_stack(tmp2(2)) 806 value = value.and.BA_pop_stack(tmp(2)) 807 if (.not. value) 808 > call errquit('cpsp_stress_readall:error popping stack',0, 809 & MA_ERR) 810 return 811 end 812