1* 2* $Id$ 3* 4*********************************************************************** 5* band_structure * 6* * 7* This is a developing band structure parallel code for NWCHEM * 8* + tcgmsg message passing library used * 9* + my own slap-decomposed parallel 3d-FFT(real->complex) used * 10* * 11* * 12*********************************************************************** 13 14 logical function band_structure(rtdb,flag) 15 implicit none 16 integer rtdb 17 integer flag 18 19 20#include "global.fh" 21#include "bafdecls.fh" 22#include "btdb.fh" 23#include "stdio.fh" 24#include "util.fh" 25#include "errquit.fh" 26 27 28* **** parallel variables **** 29 integer taskid,taskid_k,np,np_i,np_j,np_k 30 integer MASTER 31 parameter(MASTER=0) 32 33* **** timing variables **** 34 real*8 cpu1,cpu2,cpu3,cpu4 35 real*8 t1,t2,t3,t4,av 36 37* **** lattice variables **** 38 integer ngrid(3),nwave,nfft3d 39 40* ***** energy variables **** 41 integer vall(2),nn,nb,ispin,ne(2),rho(2),neall 42 real*8 E(20),en(2),ein,eke 43 real*8 dipole(3) 44 real*8 stress(3,3) 45 46 integer eigs_dos(2),dosgrid(3),weight_dos(2),pweight_dos(2) 47 integer pweight_lmax 48 49* **** gradient variables **** 50 integer fion(2) 51 52* **** error variables **** 53 logical value,ortho,mulliken 54 integer ierr 55 56* **** local variables **** 57 logical newpsi,grid3d,spin_orbit,lprint,mprint,hprint 58 real*8 gx,gy,gz,cx,cy,cz 59 real*8 EV,pi,e1,e2,f0,f1,f2,f3,f4,f5,f6,ttl1 60 real*8 pathlength,dist,kold(3),emin,emax,lmbda,rcut 61 integer ii,jj,kk,ll,i,k,ia,nion,vers,nbrillioun,icharge,isize,indx 62 integer npoints,nbrillall,if1,if2,nbrillq,l3 63 integer mapping,tmp(2),norbs_dos 64 character*255 full_filename 65 character*50 filename 66 character*72 cube_comment 67 68 69* **** external functions **** 70* **** external functions **** 71 integer cpsp_nprj 72 external cpsp_nprj 73 real*8 lattice_omega,lattice_unita,lattice_ecut,lattice_wcut 74 real*8 lattice_unitg,ion_amass,ion_TotalCharge 75 logical cpsi_spin_orbit,control_spin_orbit,control_print 76 character*4 ion_aname 77 integer control_ispin 78 external lattice_omega,lattice_unita,lattice_ecut,lattice_wcut 79 external lattice_unitg,ion_amass,ion_TotalCharge 80 external ion_aname,control_print 81 external cpsi_spin_orbit,control_spin_orbit,control_ispin 82 83 84 real*8 control_tole,control_tolc,control_tolr,ion_rion 85 external control_tole,control_tolc,control_tolr,ion_rion 86 real*8 control_time_step,control_fake_mass 87 external control_time_step,control_fake_mass 88 logical control_read,control_move,ion_init 89 external control_read,control_move,ion_init 90 integer cpsp_psp_type 91 external cpsp_psp_type 92 integer control_it_in,control_it_out,control_gga,control_version 93 integer control_ngrid,pack_nwave,ion_nion,ion_natm,ion_katm 94 integer ion_nkatm 95 external control_it_in,control_it_out,control_gga,control_version 96 external control_ngrid,pack_nwave,ion_nion,ion_natm,ion_katm 97 external ion_nkatm 98 99 character*12 control_boundry 100 external control_boundry 101 102 logical brillioun_print 103 integer brillioun_nbrillioun,Pneb_nbrillq 104 real*8 brillioun_weight_brdcst 105 real*8 brillioun_ks_brdcst,brillioun_k_brdcst 106 external brillioun_print 107 external brillioun_nbrillioun,Pneb_nbrillq 108 external brillioun_weight_brdcst 109 external brillioun_ks_brdcst,brillioun_k_brdcst 110 integer c_electron_count,linesearch_count 111 external c_electron_count,linesearch_count 112 113 real*8 nwpw_timing 114 external nwpw_timing 115 integer Cram_nwave_all_brdcst,Cram_nwave_brdcst 116 external Cram_nwave_all_brdcst,Cram_nwave_brdcst 117 118 integer ewald_ncut 119 real*8 ewald_rcut,ewald_mandelung,ewald_e 120 external ewald_ncut 121 external ewald_rcut,ewald_mandelung,ewald_e 122 logical cpsp_semicore,psi_filefind,cpsi_initialize,cpsi_finalize 123 external cpsp_semicore,psi_filefind,cpsi_initialize,cpsi_finalize 124 logical cpsi_band_finalize 125 external cpsi_band_finalize 126 real*8 c_cgsd_noit_energy,cpsi_1energy,cpsi_eigenvalue_brdcst 127 external c_cgsd_noit_energy,cpsi_1energy,cpsi_eigenvalue_brdcst 128 integer cpsp_lmax,cpsp_locp,cpsp_lmmax 129 external cpsp_lmax,cpsp_locp,cpsp_lmmax 130 real*8 cpsp_rcore,cpsp_rc,cpsp_ncore,cpsp_zv 131 external cpsp_rcore,cpsp_rc,cpsp_ncore,cpsp_zv 132 character*4 ion_atom 133 external ion_atom 134 integer cpsi_ispin,cpsi_ne,psi_get_version 135 external cpsi_ispin,cpsi_ne,psi_get_version 136 logical pspw_reformat_c_wvfnc 137 external pspw_reformat_c_wvfnc 138 integer control_mapping,control_np_dimensions 139 external control_mapping,control_np_dimensions 140 integer control_num_kvectors_structure,control_excited_ne 141 external control_num_kvectors_structure,control_excited_ne 142 logical band_dplot_iteration_check,control_Mulliken 143 external band_dplot_iteration_check,control_Mulliken 144 integer cpsi_iptr_psi, cpsi_iptr_dn, c_electron_iptr_psir 145 external cpsi_iptr_psi, cpsi_iptr_dn, c_electron_iptr_psir 146 character spdf_name 147 external spdf_name 148 character*7 c_index_name 149 external c_index_name 150 151*****************************| PROLOGUE |**************************** 152 153 value = .true. 154 pi = 4.0d0*datan(1.0d0) 155 156 call nwpw_timing_init() 157 call dcopy(20,0.0d0,0,E,1) 158 159 160* **** get parallel variables **** 161 call Parallel_Init() 162 call Parallel_np(np) 163 call Parallel_taskid(taskid) 164 if (taskid.eq.MASTER) call current_second(cpu1) 165 166 167 value = control_read(5,rtdb) 168 if (.not. value) 169 > call errquit('error reading control',0, DISK_ERR) 170 171 lprint = ((taskid.eq.MASTER).and.(control_print(print_low))) 172 mprint = ((taskid.eq.MASTER).and.(control_print(print_medium))) 173 hprint = ((taskid.eq.MASTER).and.(control_print(print_high))) 174 mulliken = control_Mulliken() 175 176* ***** print out header **** 177 if (mprint) then 178 write(luout,1000) 179 write(luout,1010) 180 write(luout,1020) 181 write(luout,1010) 182 write(luout,1040) 183 write(luout,1010) 184 write(luout,1041) 185 write(luout,1043) 186 write(luout,1010) 187 write(luout,1000) 188 call nwpw_message(1) 189 write(luout,1110) 190 end if 191 192 193 call Parallel3d_Init(control_np_dimensions(2), 194 > control_np_dimensions(3)) 195 call Parallel3d_np_i(np_i) 196 call Parallel3d_np_j(np_j) 197 call Parallel3d_np_k(np_k) 198 call Parallel3d_taskid_k(taskid_k) 199 200 201 ngrid(1) = control_ngrid(1) 202 ngrid(2) = control_ngrid(2) 203 ngrid(3) = control_ngrid(3) 204 nwave = 0 205 mapping = control_mapping() 206 207 208 ierr = 0 209 210* **** initialize C3dB data structure **** 211 call C3dB_Init(1,ngrid(1),ngrid(2),ngrid(3),mapping) 212 call C3dB_nfft3d(1,nfft3d) 213 214 call cpsi_data_init(20) 215 216* **** read ions **** 217 value = ion_init(rtdb) 218 call center_geom(cx,cy,cz) 219 call center_mass(gx,gy,gz) 220 221* **** initialize lattice data structure **** 222 call lattice_init() 223 call c_G_init() 224 225* **** initalize brillioun **** 226 call brillioun_init() 227 call Cram_Init() 228 call C3dB_pfft_init() 229 230 231* **** initialize D3dB data structure and mask for GGA **** 232 if ((control_gga().ge.10).and.(control_gga().lt.100)) THEN 233 call D3dB_Init(1,ngrid(1),ngrid(2),ngrid(3),mapping) 234 call G_init() 235 call mask_init() 236 end if 237 238c* **** read ions **** 239c value = ion_init(rtdb) 240c call center_geom(cx,cy,cz) 241c call center_mass(gx,gy,gz) 242 243* **** allocate psp data structure and read in psedupotentials into it **** 244 call cpsp_init() 245 call cpsp_readall() 246 if (cpsp_semicore(0)) call c_semicore_check() 247 248 249* **** initialize ke,and coulomb data structures **** 250 call cstrfac_init() 251 call cke_init() 252 call c_coulomb_init() 253 254 call ewald_init() 255 256* **** set up phase factors at the current geometry **** 257 call cphafac() 258 call cphafac_k() 259 call ewald_phafac() 260 261* **** read in wavefunctions and initialize psi **** 262 if (.not.psi_filefind()) then 263 call cpsi_new() 264 newpsi = .true. 265 266 else 267 newpsi = .false. 268 269* **** convert from pspw format to band format **** 270 vers = psi_get_version() 271 if ((vers.eq.3).or.(vers.eq.4)) then 272 nbrillioun = brillioun_nbrillioun() 273 newpsi = .true. 274 if (taskid.eq.MASTER) then 275 value= pspw_reformat_c_wvfnc(1) 276 end if 277 end if 278 end if 279 280 call psi_get_ne(ispin,ne) 281 if (ispin.eq.3) then 282 spin_orbit = .true. 283 ispin=2 284 else 285 spin_orbit = .false. 286 end if 287 nbrillioun = brillioun_nbrillioun() 288 call Pneb_init(ispin,ne,nbrillioun,spin_orbit) 289 value = cpsi_initialize(.true.) 290 291* **** electron and geodesic data structures **** 292 call c_electron_init() 293 call c_geodesic_init() 294 call linesearch_init() 295 call band_dplot_iteration_init() 296 297 298 299 300 301* |**************************| 302****************** summary of input data ********************** 303* |**************************| 304 305 if (mprint) then 306 write(luout,1111) np 307 write(luout,1117) np_i,np_j,np_k 308 if (mapping.eq.1) write(luout,1112) 309 if (mapping.eq.2) write(luout,1113) 310 if (mapping.eq.3) write(luout,1118) 311 312 write(luout,1115) 313 write(luout,1121) control_boundry(),control_version() 314 315 call v_bwexc_print(luout,control_gga()) 316 317 write(luout,1140) 318 do ia = 1,ion_nkatm() 319 write(luout,1150) ia,ion_atom(ia), 320 > cpsp_zv(ia),cpsp_lmax(ia) 321 write(luout,2000) cpsp_psp_type(ia) 322 write(luout,1152) cpsp_lmax(ia) 323 write(luout,1153) cpsp_locp(ia) 324 write(luout,1154) cpsp_nprj(ia) 325 if (cpsp_semicore(ia)) 326 > write(luout,1155) cpsp_rcore(ia),cpsp_ncore(ia) 327 write(luout,1151) (cpsp_rc(i,ia),i=0,cpsp_lmax(ia)) 328 end do 329 if (control_spin_orbit()) then 330 icharge=-cpsi_ne(1)+ion_TotalCharge() 331 else 332 icharge = -(cpsi_ne(1)+cpsi_ne(cpsi_ispin())) 333 > + ion_TotalCharge() 334 end if 335 write(luout,1159) icharge 336 337 write(luout,1180) 338 write(luout,1190) (I,ion_aname(I), 339 > (ion_rion(K,I),K=1,3), 340 > ion_amass(I)/1822.89d0, 341 > I=1,ion_nion()) 342 write(luout,1200) cx,cy,cz 343 write(luout,1210) gx,gy,gz 344 write(luout,1220) cpsi_ne(1),cpsi_ne(cpsi_ispin()), 345 > ' (Fourier space)' 346 347 348 write(luout,1230) 349 write(luout,1241) lattice_unita(1,1), 350 > lattice_unita(2,1), 351 > lattice_unita(3,1) 352 write(luout,1242) lattice_unita(1,2), 353 > lattice_unita(2,2), 354 > lattice_unita(3,2) 355 write(luout,1243) lattice_unita(1,3), 356 > lattice_unita(2,3), 357 > lattice_unita(3,3) 358 write(luout,1244) lattice_unitg(1,1), 359 > lattice_unitg(2,1), 360 > lattice_unitg(3,1) 361 write(luout,1245) lattice_unitg(1,2), 362 > lattice_unitg(2,2), 363 > lattice_unitg(3,2) 364 write(luout,1246) lattice_unitg(1,3), 365 > lattice_unitg(2,3), 366 > lattice_unitg(3,3) 367 write(luout,1231) lattice_omega() 368 write(luout,1260) ewald_rcut(),ewald_ncut() 369 write(luout,1261) ewald_mandelung() 370 371 ia = brillioun_nbrillioun() 372 write(luout,1255) 373 write(luout,1256) ia 374 end if 375 376c **** print brillioun zone - extra logic for distributed kpoints **** 377 if (brillioun_print()) then 378 do i=1,brillioun_nbrillioun() 379 f0 = brillioun_weight_brdcst(i) 380 f1 = brillioun_ks_brdcst(1,i) 381 f2 = brillioun_ks_brdcst(2,i) 382 f3 = brillioun_ks_brdcst(3,i) 383 f4 = brillioun_k_brdcst(1,i) 384 f5 = brillioun_k_brdcst(2,i) 385 f6 = brillioun_k_brdcst(3,i) 386 if (mprint) write(luout,1257) f0,f1,f2,f3,f4,f5,f6 387 end do 388 else 389 if (mprint) write(luout,1258) 390 end if 391 392 if1 = Cram_nwave_all_brdcst(0) 393 if2 = Cram_nwave_brdcst(0) 394 if (mprint) then 395 write(luout,1249) 396 write(luout,1250) lattice_ecut(),ngrid(1),ngrid(2),ngrid(3), 397 > if1,if2 398 end if 399 400 if (brillioun_print()) then 401 do i=1,brillioun_nbrillioun() 402 if1 = Cram_nwave_all_brdcst(i) 403 if2 = Cram_nwave_brdcst(i) 404 if (mprint) then 405 write(luout,1251) i,lattice_wcut(),ngrid(1),ngrid(2),ngrid(3), 406 > if1,if2 407 end if 408 end do 409 else 410 if (mprint) write(luout,1252) lattice_wcut() 411 end if 412 413 414 if (mprint) then 415 write(luout,1270) 416 write(luout,1280) control_time_step(),control_fake_mass() 417 write(luout,1290) control_tole(),control_tolc() 418 write(luout,1300) 419 call util_flush(luout) 420 call flush(luout) 421 end if 422 423 424 425 !**** set the size of the band - include virtual orbitals **** 426 ispin = cpsi_ispin() 427 ne(1) = cpsi_ne(1)+control_excited_ne(1) 428 ne(2) = 0 429 if (ispin.gt.1) ne(2) = cpsi_ne(2)+control_excited_ne(2) 430 431 432 if (taskid.eq.MASTER) call current_second(cpu2) 433 434* **** allocate vall **** 435 value = BA_alloc_get(mt_dcpl,2*nfft3d,'vall',vall(2),vall(1)) 436 if (.not. value) 437 > call errquit('band_structure:out of heap memory',0, MA_ERR) 438 439 440 441 EV = c_cgsd_noit_energy() 442 call c_electron_gen_vall() 443 call c_electron_get_vall(dcpl_mb(vall(1))) 444 445 if (taskid.eq.MASTER) then 446 write(luout,1600) EV 447 write(luout,*) 448 write(luout,*) "Self-Consistent Potential Generated" 449 end if 450 451 value=cpsi_finalize(.true.) 452 call c_electron_finalize() 453 call c_geodesic_finalize() 454 call ewald_end() 455 call cstrfac_end() 456 call c_coulomb_end() 457 call cke_end() 458 call cpsp_end() 459 call C3dB_pfft_end() 460 call Cram_end() 461 call c_G_end() 462 call brillioun_end() 463 464 465* **** produce eigenvalue band file(s) **** 466 if (ispin.eq.1) then 467 call util_file_name('restricted_band', 468 > .false., 469 > .false., 470 > full_filename) 471 if (taskid.eq.MASTER) then 472 open(unit=58,file=full_filename,form='formatted') 473 end if 474 else 475 if (cpsi_spin_orbit()) then 476 call util_file_name('spinor_band', 477 > .false., 478 > .false., 479 > full_filename) 480 if (taskid.eq.MASTER) then 481 open(unit=58,file=full_filename,form='formatted') 482 end if 483 else 484 call util_file_name('alpha_band', 485 > .false., 486 > .false., 487 > full_filename) 488 if (taskid.eq.MASTER) then 489 open(unit=58,file=full_filename,form='formatted') 490 end if 491 call util_file_name('beta_band', 492 > .false., 493 > .false., 494 > full_filename) 495 if (taskid.eq.MASTER) then 496 open(unit=59,file=full_filename,form='formatted') 497 end if 498 end if 499 end if 500 501* **** DOS calculation **** 502 if (flag.eq.1) then 503 if (taskid.eq.MASTER) write(luout,*) "DOS of states calculation" 504 call control_dos_grid_structure(dosgrid) 505 506 507 508* **** DOS_dplot calculation **** 509 else if (flag.eq.2) then 510 if (taskid.eq.MASTER) 511 > write(luout,*) "DOS_dplot of states calculation" 512 513 call control_dos_grid_structure(dosgrid) 514 isize = dosgrid(1)*dosgrid(2)*dosgrid(3)*(ne(1)+ne(2)) 515 value = BA_alloc_get(mt_dbl,isize, 516 > 'weight_dos',weight_dos(2),weight_dos(1)) 517 value = value.and.BA_alloc_get(mt_dbl,isize, 518 > 'eigs_dos',eigs_dos(2),eigs_dos(1)) 519 value = value.and.BA_alloc_get(mt_dbl,ispin*nfft3d, 520 > 'rho',rho(2),rho(1)) 521 if (.not. value) 522 > call errquit('band_structure:out of heap memory',0, MA_ERR) 523 524 !**** get eigs_dos from rtdb **** 525 if (.not.btdb_get(rtdb,'dos:eigs_dos',mt_dbl, 526 > isize,dbl_mb(eigs_dos(1)))) then 527 call errquit('band_structure:cannot read eigs_dos from rtdb', 528 > 0,RTDB_ERR) 529 end if 530 531 !**** get dos:ein from rtdb **** 532 if (.not.btdb_get(rtdb,'dos:ein',mt_dbl,1,ein)) then 533 call errquit('band_structure:cannot read dos:ein from rtdb', 534 > 0,RTDB_ERR) 535 end if 536 537 call band_dos_weights_generate(dosgrid(1),dosgrid(2),dosgrid(3), 538 > dbl_mb(eigs_dos(1)),ne(1), 539 > ein, 540 > dbl_mb(weight_dos(1))) 541 ii = dosgrid(1)*dosgrid(2)*dosgrid(3)*ne(1) 542 if (ispin.eq.2) 543 > call band_dos_weights_generate(dosgrid(1),dosgrid(2),dosgrid(3), 544 > dbl_mb(eigs_dos(1)+ii),ne(2), 545 > ein, 546 > dbl_mb(weight_dos(1)+ii)) 547 call dcopy(ispin*nfft3d,0.0d0,0,dbl_mb(rho(1)),1) 548 549 if (taskid.eq.MASTER) then 550 do k=1,dosgrid(1)*dosgrid(2)*dosgrid(3) 551 write(*,*) 552 write(*,*) "brillioun k=",k 553 write(*,*) "---------------------" 554 do ii=1,ne(1) 555 write(*,*) "weight=",k,ii, 556 > dbl_mb(weight_dos(1) 557 > +(ii-1) 558 > +(k-1)*ne(1)) 559 call flush(6) 560 end do 561 end do 562 end if 563 564 565 566* **** band structure calculation **** 567 else 568 if (taskid.eq.MASTER) 569 > write(luout,*) "band structure calculation" 570 end if 571 572 nbrillall = control_num_kvectors_structure() 573 call control_reset_band_structure() 574 kk = 0 575 do k=1,nbrillall,np_k 576 ortho = .true. 577 kk = kk + np_k 578 if (kk.gt.nbrillall) kk = nbrillall 579 580 if (taskid.eq.MASTER) then 581 do ii=k,kk 582 write(luout,1301) ii 583 end do 584 end if 585 586 587* **** initialize lattice data structure **** 588 call lattice_init() 589 call c_G_init() 590 call brillioun_structure_init(k,kk-k+1) 591 call Cram_Init() 592 call C3dB_pfft_init() 593 594* **** allocate psp data structure and read in psedupotentials into it **** 595 call cpsp_init() 596 call cpsp_readall() 597 if (cpsp_semicore(0)) call c_semicore_check() 598 599* **** initialize ke,and coulomb data structures **** 600 call cstrfac_init() 601 call cke_init() 602 call c_coulomb_init() 603 call ewald_init() 604 605* **** set up phase factors at the current geometry **** 606 call cphafac() 607 call cphafac_k() 608 call ewald_phafac() 609 610* **** read in wavefunctions and initialize psi **** 611 if (.not.psi_filefind()) then 612 call cpsi_new_ne(ispin,ne) 613 newpsi = .true. 614 ortho = .false. 615 616 else 617 newpsi = .false. 618 619* **** convert from pspw format to band format **** 620 vers = psi_get_version() 621 if ((vers.eq.3).or.(vers.eq.4)) then 622 nbrillioun = brillioun_nbrillioun() 623 newpsi = .true. 624 if (taskid.eq.MASTER) then 625 value= pspw_reformat_c_wvfnc(1) 626 end if 627 end if 628 end if 629 630 call psi_get_ne(ispin,ne) 631 if (ispin.eq.3) then 632 spin_orbit = .true. 633 ispin=2 634 else 635 spin_orbit = .false. 636 end if 637 nbrillioun = brillioun_nbrillioun() 638 639 call Pneb_init(ispin,ne,nbrillioun,spin_orbit) 640 value = cpsi_initialize(ortho) 641 642c if (flag.eq.2) call cpsi_dospsi_read(k) 643 644 645* **** allocate eigs_dos if first iteration and band structure calculation **** 646 if ((flag.eq.1).and.(k.eq.1)) then 647 isize = dosgrid(1)*dosgrid(2)*dosgrid(3)*(cpsi_ne(1)+cpsi_ne(2)) 648 value = BA_alloc_get(mt_dbl,isize, 649 > 'eigs_dos',eigs_dos(2),eigs_dos(1)) 650 if (mulliken) then 651 if (.not.btdb_get(rtdb,'nwpw:dos:orb:norb', 652 > mt_int,1,norbs_dos)) 653 > norbs_dos = 0 654 value = value.and.BA_alloc_get(mt_dbl,isize*(4+norbs_dos), 655 > 'pweight_dos',pweight_dos(2),pweight_dos(1)) 656 end if 657 if (.not. value) 658 > call errquit('band_structure:out of heap memory',0, MA_ERR) 659 end if 660 661 662* **** electron and geodesic data structures **** 663 call c_electron_init() 664 call c_geodesic_init() 665 call linesearch_init() 666 667* **** diagonalize hamiltonian and rotate psi **** 668 call c_electron_set_vall(dcpl_mb(vall(1))) 669 670* **** initialize with steepest descent *** 671 call c_sdminimize_noscf(0) 672 673* **** diagonalize current result *** 674 call cpsi_1gen_hml() 675 call cpsi_diagonalize_hml() 676 call cpsi_1rotate2() 677 call cpsi_2to1() 678 679 call cpsi_KS_minimize(1,.false.,control_tole(),control_tolc()) 680 681 if (control_print(print_high)) call cpsi_check_indx(k) 682 683 if (band_dplot_iteration_check(k)) then 684 call band_dplot_iteration(k,ispin,ne,1, 685 > dcpl_mb(cpsi_iptr_psi(1)), 686 > dbl_mb(cpsi_iptr_dn(1)), 687 > dcpl_mb(c_electron_iptr_psir())) 688 end if 689 690 691 NN=cpsi_ne(1)-cpsi_ne(2) 692 EV=27.2116d0 693 do nb=1,brillioun_nbrillioun() 694 695 f1 = brillioun_ks_brdcst(1,nb) 696 f2 = brillioun_ks_brdcst(2,nb) 697 f3 = brillioun_ks_brdcst(3,nb) 698 f4 = brillioun_k_brdcst(1,nb) 699 f5 = brillioun_k_brdcst(2,nb) 700 f6 = brillioun_k_brdcst(3,nb) 701 if ((k+nb).eq.2) then 702 pathlength = 0.0d0 703 else 704 dist=dsqrt((f4-kold(1))**2+(f5-kold(2))**2+(f6-kold(3))**2) 705 pathlength = pathlength + dist 706 end if 707 kold(1) = f4 708 kold(2) = f5 709 kold(3) = f6 710 711 if (taskid.eq.MASTER) then 712 write(luout,1508) k+nb-1,pathlength,f1,f2,f3,f4,f5,f6 713 write(luout,1500) 714 end if 715 716 !*** flag==2 *** 717 if (flag.eq.2) then 718 719 !*** spin-orbit **** 720 if (cpsi_spin_orbit()) then 721 do i=0,cpsi_ne(1)-1 722 ii=cpsi_ne(1)-i 723 e1 = cpsi_eigenvalue_brdcst(nb,1,cpsi_ne(1)-i) 724 if (taskid.eq.MASTER) 725 > write(luout,1511) e1,e1*EV, 726 > dbl_mb(weight_dos(1)+ii+(k+nb-1)*cpsi_ne(1)) 727 end do 728 729 !*** not spin-orbit **** 730 else 731 do i=0,NN-1 732 ii = cpsi_ne(1)-i 733 e1 = cpsi_eigenvalue_brdcst(nb,1,cpsi_ne(1)-i) 734 if (taskid.eq.MASTER) 735 > write(luout,1511) e1,e1*EV, 736 > dbl_mb(weight_dos(1)+(ii-1)+(k+nb-1)*ne(1)) 737 end do 738 do i=0,cpsi_ne(2)-1 739 ii = cpsi_ne(1)-i 740 jj = cpsi_ne(2)-i + 741 > dosgrid(1)*dosgrid(2)*dosgrid(3)*ne(1) 742 e1 = cpsi_eigenvalue_brdcst(nb,1,cpsi_ne(1)-i-NN) 743 e2 = cpsi_eigenvalue_brdcst(nb,2,cpsi_ne(2)-i) 744 if (taskid.eq.MASTER) 745 > write(luout,1511) e1,e1*EV, 746 > dbl_mb(weight_dos(1)+(ii-1)+(k+nb-2)*ne(1)), 747 > e2,e2*EV, 748 > dbl_mb(weight_dos(1)+(jj-1)+(k+nb-2)*ne(2)) 749 end do 750 end if 751 752 !*** flag!=2 *** 753 else 754 755 !*** spin-orbit **** 756 if (cpsi_spin_orbit()) then 757 do i=0,cpsi_ne(1)-1 758 e1 = cpsi_eigenvalue_brdcst(nb,1,cpsi_ne(1)-i) 759 if (taskid.eq.MASTER) 760 > write(luout,1510) e1,e1*EV 761 end do 762 763 !*** not spin-orbit **** 764 else 765 do i=0,NN-1 766 e1 = cpsi_eigenvalue_brdcst(nb,1,cpsi_ne(1)-i) 767 if (taskid.eq.MASTER) 768 > write(luout,1510) e1,e1*EV 769 end do 770 do i=0,cpsi_ne(2)-1 771 e1 = cpsi_eigenvalue_brdcst(nb,1,cpsi_ne(1)-i-NN) 772 e2 = cpsi_eigenvalue_brdcst(nb,2,cpsi_ne(2)-i) 773 if (taskid.eq.MASTER) 774 > write(luout,1510) e1,e1*EV,e2,e2*EV 775 end do 776 end if 777 778 end if 779 780 !*** set eigs_dos *** 781 if (flag.eq.1) then 782 do i=1,cpsi_ne(1) 783 indx = eigs_dos(1) 784 > + (k+nb-2) 785 > + (i-1)*dosgrid(1)*dosgrid(2)*dosgrid(3) 786 dbl_mb(indx) = cpsi_eigenvalue_brdcst(nb,1,i) 787 end do 788 do i=1,cpsi_ne(2) 789 indx = eigs_dos(1) 790 > + (k+nb-2) 791 > + (i-1+cpsi_ne(1))*dosgrid(1)*dosgrid(2)*dosgrid(3) 792 dbl_mb(indx) = cpsi_eigenvalue_brdcst(nb,2,i) 793 end do 794 end if 795 796 if (.not.BA_push_get(mt_dbl,cpsi_ne(1),'tmp',tmp(2),tmp(1))) 797 > call errquit('band_structure:push stack',99,MA_ERR) 798 799 do i=1,cpsi_ne(1) 800 dbl_mb(tmp(1)+i-1) = cpsi_eigenvalue_brdcst(nb,1,i) 801 end do 802 if (taskid.eq.MASTER) 803 > write(58,'(1000E14.6)') pathlength, 804 > (dbl_mb(tmp(1)+i-1),i=1,cpsi_ne(1)) 805 806 if ((.not.cpsi_spin_orbit()).and.(ispin.eq.2)) then 807 do i=1,cpsi_ne(2) 808 dbl_mb(tmp(1)+i-1) = cpsi_eigenvalue_brdcst(nb,2,i) 809 end do 810 if (taskid.eq.MASTER) 811 > write(59,'(1000E14.6)') pathlength, 812 > (dbl_mb(tmp(1)+i-1),i=1,cpsi_ne(2)) 813 end if 814 if (.not.BA_pop_stack(tmp(2))) 815 > call errquit('band_structure:pop stack',99,MA_ERR) 816 817 end do !*** nb *** 818 819 !*** set rho *** 820 if ((flag.eq.2).and.((k+taskid_k).le.nbrillall)) then 821 ii = (k+taskid_k-1)*ne(1) 822 call c_electron_gen_weighted_density(1, 823 > dbl_mb(weight_dos(1)+ii), 824 > dbl_mb(rho(1))) 825 ii = (k+taskid_k-1)*ne(2)+dosgrid(1)*dosgrid(2)*dosgrid(3)*ne(1) 826 if (ispin.eq.2) 827 > call c_electron_gen_weighted_density(2, 828 > dbl_mb(weight_dos(1)+ii), 829 > dbl_mb(rho(1)+nfft3d)) 830 end if 831 832* **** set pweight_dos **** 833 if ((flag.eq.1).and.mulliken) then 834 call cpsi_projected_dos_weights(rtdb,dosgrid,k, 835 > dbl_mb(pweight_dos(1)), 836 > pweight_lmax,norbs_dos) 837 end if 838 839 840* **** writeout dospsi -- needed for task band dos_dplot **** 841cccc if (flag.eq.1) call cpsi_dospsi_write(k) 842 843 844* **** finalize and deallocate cpsi **** 845 value = cpsi_finalize(.true.) 846 847* **** deallocate heap memory **** 848 call c_electron_finalize() 849 call c_geodesic_finalize() 850 call ewald_end() 851 call cstrfac_end() 852 call c_coulomb_end() 853 call cke_end() 854 call cpsp_end() 855 call C3dB_pfft_end() 856 call Cram_end() 857 call c_G_end() 858 call brillioun_end() 859 860 end do 861 if (taskid.eq.MASTER) then 862 close(58) 863 if (ispin.eq.2) close(59) 864 end if 865 866 if (taskid.eq.MASTER) call current_second(cpu3) 867 868* |***************************| 869****************** DOS plotting ********************** 870* |***************************| 871 872 value = btdb_parallel(.false.) 873 if ((flag.eq.1).and.(taskid.eq.MASTER)) then 874 875 876 if (.not.btdb_get(rtdb,'dos:npoints',mt_int,1,npoints)) then 877 npoints = 500 878 end if 879 880 if (.not.btdb_get(rtdb,'dos:emin',mt_dbl,1,emin)) then 881 emin = 99999.0d0 882 do ii=1,(ne(1)+ne(2))*dosgrid(1)*dosgrid(2)*dosgrid(3) 883 if (dbl_mb(eigs_dos(1)+ii-1).lt.emin) 884 > emin = dbl_mb(eigs_dos(1)+ii-1) 885 end do 886 emin = emin - 0.1d0 887 end if 888 889 if (.not.btdb_get(rtdb,'dos:emax',mt_dbl,1,emax)) then 890 emax = -99999.0d0 891 do ii=1,(ne(1)+ne(2))*dosgrid(1)*dosgrid(2)*dosgrid(3) 892 if (dbl_mb(eigs_dos(1)+ii-1).gt.emax) 893 > emax = dbl_mb(eigs_dos(1)+ii-1) 894 end do 895 emax = emax + 0.1d0 896 end if 897 898 call util_file_name('dos', 899 > .false., 900 > .false., 901 > full_filename) 902 open(unit=58,file=full_filename,form='formatted') 903 call band_dos_generate(58,dosgrid(1),dosgrid(2),dosgrid(3), 904 > dbl_mb(eigs_dos(1)),ne(1)+ne(2), 905 > (2.0d0*(3-ispin)), 906 > npoints,emin,emax) 907 close(58) 908 909 if (ispin.eq.2) then 910 call util_file_name('dos_alpha', 911 > .false., 912 > .false., 913 > full_filename) 914 open(unit=58,file=full_filename,form='formatted') 915 call band_dos_generate(58,dosgrid(1),dosgrid(2),dosgrid(3), 916 > dbl_mb(eigs_dos(1)),ne(1), 917 > (1.0d0), 918 > npoints,emin,emax) 919 close(58) 920 call util_file_name('dos_beta', 921 > .false., 922 > .false., 923 > full_filename) 924 open(unit=58,file=full_filename,form='formatted') 925 ii = dosgrid(1)*dosgrid(2)*dosgrid(3)*ne(1) 926 call band_dos_generate(58,dosgrid(1),dosgrid(2),dosgrid(3), 927 > dbl_mb(eigs_dos(1)+ii),ne(2), 928 > (-1.0d0), 929 > npoints,emin,emax) 930 close(58) 931 end if 932 933 if (mulliken) then 934 do ll=0,pweight_lmax 935 call util_file_name('dos_both_'//spdf_name(ll), 936 > .false., 937 > .false., 938 > full_filename) 939 open(unit=58,file=full_filename,form='formatted') 940 jj=dosgrid(1)*dosgrid(2)*dosgrid(3)*(ne(1)+ne(2))*ll 941 call band_projected_dos_generate(58, 942 > dosgrid(1),dosgrid(2),dosgrid(3), 943 > dbl_mb(eigs_dos(1)), 944 > dbl_mb(pweight_dos(1)+jj),ne(1)+ne(2), 945 > (1.0d0*(3-ispin)), 946 > npoints,emin,emax) 947 close(58) 948 949 if (ispin.eq.2) then 950 call util_file_name('dos_alpha_'//spdf_name(ll), 951 > .false., 952 > .false., 953 > full_filename) 954 open(unit=58,file=full_filename,form='formatted') 955 jj=dosgrid(1)*dosgrid(2)*dosgrid(3)*(ne(1)+ne(2))*ll 956 call band_projected_dos_generate(58, 957 > dosgrid(1),dosgrid(2),dosgrid(3), 958 > dbl_mb(eigs_dos(1)), 959 > dbl_mb(pweight_dos(1)+jj),ne(1), 960 > (1.0d0), 961 > npoints,emin,emax) 962 close(58) 963 call util_file_name('dos_beta_'//spdf_name(ll), 964 > .false., 965 > .false., 966 > full_filename) 967 open(unit=58,file=full_filename,form='formatted') 968 ii = dosgrid(1)*dosgrid(2)*dosgrid(3)*ne(1) 969 jj=dosgrid(1)*dosgrid(2)*dosgrid(3)*(ne(1)+ne(2))*ll + ii 970 call band_projected_dos_generate(58, 971 > dosgrid(1),dosgrid(2),dosgrid(3), 972 > dbl_mb(eigs_dos(1)+ii), 973 > dbl_mb(pweight_dos(1)+jj),ne(2), 974 > (-1.0d0), 975 > npoints,emin,emax) 976 close(58) 977 end if 978 end do 979 980 !*** ORBITAL DOS *** 981 do ll=1,norbs_dos 982 l3 = ll+pweight_lmax+1 983 call util_file_name('dos_both_orb'//c_index_name(ll), 984 > .false., 985 > .false., 986 > full_filename) 987 open(unit=58,file=full_filename,form='formatted') 988 jj=dosgrid(1)*dosgrid(2)*dosgrid(3)*(ne(1)+ne(2))*l3 989 call band_projected_dos_generate(58, 990 > dosgrid(1),dosgrid(2),dosgrid(3), 991 > dbl_mb(eigs_dos(1)), 992 > dbl_mb(pweight_dos(1)+jj),ne(1)+ne(2), 993 > (1.0d0*(3-ispin)), 994 > npoints,emin,emax) 995 close(58) 996 997 if (ispin.eq.2) then 998 call util_file_name('dos_alpha_orb'//c_index_name(ll), 999 > .false., 1000 > .false., 1001 > full_filename) 1002 open(unit=58,file=full_filename,form='formatted') 1003 jj=dosgrid(1)*dosgrid(2)*dosgrid(3)*(ne(1)+ne(2))*l3 1004 call band_projected_dos_generate(58, 1005 > dosgrid(1),dosgrid(2),dosgrid(3), 1006 > dbl_mb(eigs_dos(1)), 1007 > dbl_mb(pweight_dos(1)+jj),ne(1), 1008 > (1.0d0), 1009 > npoints,emin,emax) 1010 close(58) 1011 call util_file_name('dos_beta_orb'//c_index_name(ll), 1012 > .false., 1013 > .false., 1014 > full_filename) 1015 open(unit=58,file=full_filename,form='formatted') 1016 ii = dosgrid(1)*dosgrid(2)*dosgrid(3)*ne(1) 1017 jj=dosgrid(1)*dosgrid(2)*dosgrid(3)*(ne(1)+ne(2))*l3+ii 1018 call band_projected_dos_generate(58, 1019 > dosgrid(1),dosgrid(2),dosgrid(3), 1020 > dbl_mb(eigs_dos(1)+ii), 1021 > dbl_mb(pweight_dos(1)+jj),ne(2), 1022 > (-1.0d0), 1023 > npoints,emin,emax) 1024 close(58) 1025 end if 1026 1027 end do 1028 1029 1030 end if 1031 1032 1033 !**** put eigs_dos on rtdb for use by task band dos_dplot *** 1034 if (.not.btdb_put(rtdb,'dos:eigs_dos',mt_dbl, 1035 > isize,dbl_mb(eigs_dos(1)))) then 1036 call errquit('band_structure:cannot write eigs_dos to rtdb', 1037 > 0,RTDB_ERR) 1038 end if 1039 end if 1040 value = btdb_parallel(.true.) 1041 1042 1043* |***************************| 1044****************** DOS_dplot plotting ********************** 1045* |***************************| 1046 if (flag.eq.2) then 1047 1048 1049 grid3d = .false. 1050 if (btdb_get(rtdb,'band_dplot:3d_grid:nx',mt_int,1,i)) 1051 > grid3d = .true. 1052 1053 if (.not.btdb_cget(rtdb,'dos:dplot_up',1,filename)) then 1054 filename = 'dos_up.cube ' 1055 end if 1056 indx = index(filename,' ') - 1 1057 write(cube_comment,'(A,F8.3)') "dos up density, e=",ein 1058 write(*,*) ' writing dos up density E=',ein, 1059 > ' to filename: ',filename(1:11) 1060 if (grid3d) then 1061 call band_dplot_gcube_write3d(rtdb,filename, 1062 > -1,cube_comment,dbl_mb(rho(1))) 1063 else 1064 call band_dplot_gcube_write(rtdb,filename, 1065 > -1,cube_comment,dbl_mb(rho(1))) 1066 endif 1067 1068 if (ispin.eq.2) then 1069 if (.not.btdb_cget(rtdb,'dos:dplot_dn',1,filename)) then 1070 filename = 'dos_dn.cube ' 1071 end if 1072 indx = index(filename,' ') - 1 1073 write(cube_comment,'(A,F8.3)') "dos down density, e=",ein 1074 write(*,*) ' writing dos down density E=',ein, 1075 > ' to filename: ',filename(1:11) 1076 if (grid3d) then 1077 call band_dplot_gcube_write3d(rtdb,filename, 1078 > -2,cube_comment,dbl_mb(rho(1)+nfft3d)) 1079 else 1080 call band_dplot_gcube_write(rtdb,filename, 1081 > -2,cube_comment,dbl_mb(rho(1)+nfft3d)) 1082 endif 1083 end if 1084 1085 end if 1086 1087 1088* |***************************| 1089****************** Epilogue ********************** 1090* |***************************| 1091 1092 1093 1094* **** deallocate heap memory **** 1095 1096 value = BA_free_heap(vall(2)) 1097 1098 !*** deallocate eigs_dos and pweight_dos **** 1099 if (flag.eq.1) then 1100 value = value.and.BA_free_heap(eigs_dos(2)) 1101 if (mulliken) then 1102 value = value.and.BA_free_heap(pweight_dos(2)) 1103 end if 1104 end if 1105 1106 if (flag.eq.2) then 1107 value = value.and.BA_free_heap(weight_dos(2)) 1108 value = value.and.BA_free_heap(eigs_dos(2)) 1109 value = value.and.BA_free_heap(rho(2)) 1110 end if 1111 1112 call ion_write(rtdb) 1113 call ion_end() 1114 call cpsi_data_end() 1115 call C3dB_end(1) 1116 IF ((control_gga().ge.10).and.(control_gga().lt.100)) THEN 1117 call mask_end() 1118 call G_end() 1119 call D3dB_end(1) 1120 end if 1121 1122 1123* |***************************| 1124****************** report consumed cputime ********************** 1125* |***************************| 1126 if (lprint) then 1127 CALL current_second(cpu4) 1128 1129 T1=CPU2-CPU1 1130 T2=CPU3-CPU2 1131 T3=CPU4-CPU3 1132 T4=CPU4-CPU1 1133 write(luout,*) 1134 write(luout,*) '-----------------' 1135 write(luout,*) 'cputime in seconds' 1136 write(luout,*) 'prologue : ',T1 1137 write(luout,*) 'main loop : ',T2 1138 write(luout,*) 'epilogue : ',T3 1139 write(luout,*) 'total : ',T4 1140 1141 write(luout,*) 1142 write(luout,*) '-------------------------------' 1143 write(luout,*) 'Time spent doing:' 1144 write(luout,*) ' FFTs : ', 1145 > nwpw_timing(1) 1146 write(luout,*) ' dot products : ', 1147 > nwpw_timing(2) 1148 write(luout,*) ' geodesic : ', 1149 > nwpw_timing(10) 1150 write(luout,*) ' exchange correlation : ', 1151 > nwpw_timing(4) 1152 write(luout,*) ' local pseudopotentials : ', 1153 > nwpw_timing(5) 1154 write(luout,*) ' non-local pseudopotentials : ', 1155 > nwpw_timing(6) 1156 write(luout,*) ' hartree potentials : ', 1157 > nwpw_timing(7) 1158 write(luout,*) ' structure factors : ', 1159 > nwpw_timing(8) 1160 write(luout,*) ' masking and packing : ', 1161 > nwpw_timing(9) 1162 write(luout,*) 1163 CALL nwpw_MESSAGE(4) 1164 end if 1165 1166 call Parallel3d_Finalize() 1167 call Parallel_Finalize() 1168 band_structure = value 1169 return 1170 1171 1172*::::::::::::::::::::::::::: format ::::::::::::::::::::::::::::::::: 1173 1000 FORMAT(10X,'****************************************************') 1174 1010 FORMAT(10X,'* *') 1175 1020 FORMAT(10X,'* NWPW Band Structure Calculation *') 1176 1040 FORMAT(10X,'* version #2.00 1/20/07 *') 1177 1041 FORMAT(10X,'* Developed by Eric J. Bylaska *') 1178 1043 FORMAT(10X,'* and Patrick Nichols *') 1179 1100 FORMAT(//) 1180 1110 FORMAT(10X,'================ input data ========================') 1181 1111 FORMAT(/' number of processors used:',I16) 1182 1112 FORMAT( ' parallel mapping : 1d slab') 1183 1113 FORMAT( ' parallel mapping : 2d hilbert') 1184 1115 FORMAT(/' options:') 1185 1117 FORMAT( ' processor grid :',I4,' x',I4,' x',I4) 1186 1118 FORMAT( ' parallel mapping : 2d hcurve') 1187 1120 FORMAT(5X,' ionic motion = ',A) 1188 1121 FORMAT(5X,' boundary conditions = ',A,'(version', I1,')') 1189 1130 FORMAT(5X,' electron spin = ',A) 1190 1131 FORMAT(5X,' exchange-correlation = ',A) 1191 1140 FORMAT(/' elements involved in the cluster:') 1192 1150 FORMAT(5X,I2,': ',A4,' core charge:',F4.1,' lmax=',I1) 1193 1151 FORMAT(5X,' cutoff =',4F8.3) 1194 1152 FORMAT(12X,' highest angular component : ',i2) 1195 1153 FORMAT(12X,' local potential used : ',i2) 1196 1154 FORMAT(12X,' number of non-local projections: ',i2) 1197 1155 FORMAT(12X,' semicore corrections included : ', 1198 > F6.3,' (radius) ',F6.3,' (charge)') 1199 1156 FORMAT(12X,' aperiodic cutoff radius : ',F6.3) 1200 1159 FORMAT(/' total charge:',I2) 1201 1160 FORMAT(/' atomic composition:') 1202 1170 FORMAT(7(5X,A4,':',I3)) 1203 1180 FORMAT(/' initial position of ions:') 1204 1190 FORMAT(5X, I4, A5 ,' (',3F11.5,' ) - atomic mass= ',F7.3,' ') 1205 1200 FORMAT(5X,' G.C. ',' (',3F11.5,' )') 1206 1210 FORMAT(5X,' C.O.M.',' (',3F11.5,' )') 1207 1220 FORMAT(/' number of electrons: spin up=',I4,' spin down=',I4,A) 1208 1230 FORMAT(/' supercell:') 1209 1231 FORMAT(5x,' volume : ',F10.1) 1210 1241 FORMAT(5x,' lattice: a1=<',3f8.3,' >') 1211 1242 FORMAT(5x,' a2=<',3f8.3,' >') 1212 1243 FORMAT(5x,' a3=<',3f8.3,' >') 1213 1244 FORMAT(5x,' reciprocal: b1=<',3f8.3,' >') 1214 1245 FORMAT(5x,' b2=<',3f8.3,' >') 1215 1246 FORMAT(5x,' b3=<',3f8.3,' >') 1216 1217 1249 FORMAT(/' computational grids:') 1218 1250 FORMAT(5X,' density cutoff=',F7.3,' fft=',I4,'x',I4,'x',I4, 1219 & '( ',I8,' waves ',I8,' per task)') 1220 1251 FORMAT(5X,' wavefnc ',I4,' cutoff=',F7.3, 1221 & ' fft=',I4,'x',I4,'x',I4, 1222 & '( ',I8,' waves ',I8,' per task)') 1223 1252 FORMAT(5x,' wavefnc cutoff=',F7.3, 1224 > ' wavefunction grids not printed - ', 1225 > 'number of k-points is very large') 1226 1255 FORMAT(/' brillouin zone:') 1227 1256 FORMAT(5x,'number of zone points:',I6) 1228 1257 FORMAT(5x,' weight=',f8.3,' ks=<',3f8.3,' >, k=<',3f8.3,'>') 1229 1258 FORMAT(5x,' number of k-points is very large') 1230 1231 1260 FORMAT(5X,' Ewald summation: cut radius=',F8.2,' and',I3) 1232 1261 FORMAT(5X,' madelung=',f14.8) 1233 1234 1270 FORMAT(/' technical parameters:') 1235 1280 FORMAT(5X, ' time step=',F10.2,5X,'fictitious mass=',F10.1) 1236 1290 FORMAT(5X, ' tolerance=',E8.3,' (energy)',E12.3, 1237 & ' (density)') 1238 1300 FORMAT(//) 1239 1301 FORMAT(//'== Optimizing Brillouin Zone Point:',I6,' =='/) 1240 1304 FORMAT(/) 1241 1305 FORMAT(10X,'================ iteration =========================') 1242 1310 FORMAT(I8,E20.10,3E15.5) 1243 1320 FORMAT(' number of electrons: spin up=',F11.5,' down=',F11.5,A) 1244 1330 FORMAT(/' comparison between hamiltonian and lambda matrix') 1245 1340 FORMAT(I3,2I3,' H=',E16.7,', L=',E16.7,', H-L=',E16.7) 1246 1350 FORMAT(/' orthonormality') 1247 1360 FORMAT(I3,2I3,E18.7) 1248 1370 FORMAT(I3) 1249 1380 FORMAT(' ''',a,'''',I4) 1250 1390 FORMAT(I3) 1251 1400 FORMAT(I3,3E18.8/3X,3E18.8) 1252 1410 FORMAT(10X,'============= summary of results =================') 1253 1420 FORMAT( ' final position of ions:') 1254 1430 FORMAT(/' total energy :',E19.10,' (',E15.5,'/ion)') 1255 1440 FORMAT( ' total orbital energy:',E19.10,' (',E15.5,'/electron)') 1256 1450 FORMAT( ' hartree energy :',E19.10,' (',E15.5,'/electron)') 1257 1460 FORMAT( ' exc-corr energy :',E19.10,' (',E15.5,'/electron)') 1258 1470 FORMAT( ' ion-ion energy :',E19.10,' (',E15.5,'/ion)') 1259 1480 FORMAT(/' K.S. kinetic energy :',E19.10,' (',E15.5,'/electron)') 1260 1490 FORMAT( ' K.S. V_l energy :',E19.10,' (',E15.5,'/electron)') 1261 1495 FORMAT( ' K.S. V_nl energy :',E19.10,' (',E15.5,'/electron)') 1262 1496 FORMAT( ' K.S. V_Hart energy :',E19.10,' (',E15.5,'/electron)') 1263 1497 FORMAT( ' K.S. V_xc energy :',E19.10,' (',E15.5,'/electron)') 1264 1498 FORMAT( ' Virial Coefficient :',E19.10) 1265 1500 FORMAT(/' orbital energies:') 1266 1508 FORMAT(/' Brillouin zone point: ',i6, 1267 > /'pathlength=',f10.6, 1268 > /' k =<',3f8.3,'> . <b1,b2,b3> ', 1269 > /' =<',3f8.3,'>') 1270 1510 FORMAT(2(E18.7,' (',F8.3,'eV)')) 1271 1511 FORMAT(2(E18.7,' (',F8.3,'eV) dplot weight=',F8.3)) 1272 1600 FORMAT(/' Total BAND energy :',E19.10) 1273 2000 FORMAT(12X,' pseudpotential type : ',i2) 1274 2005 FORMAT(12x,' number of non local projectors : ',i3) 1275 9010 FORMAT(//' >> job terminated due to code =',I3,' <<') 1276 1277 9000 if (taskid.eq.MASTER) write(6,9010) ierr 1278 call Parallel_Finalize() 1279 1280 band_structure = value 1281 return 1282 END 1283 1284**************** Definitions of Cubes and Tetrahedrons ***************** 1285* * 1286* * 1287* (011)------------(111) ( 3 )------------( 7 ) * 1288* + + + + * 1289* /. /| /. /| * 1290* / . / | / . / | * 1291* / . / | / . / | * 1292* / . / | / . / | * 1293* (001)------------(101) | <====> ( 1 )------------( 5 ) | * 1294* | . | | | . | | * 1295* | (010)..........|.(110) | ( 2 )..........|.( 6 ) * 1296* | . | / | . | / * 1297* | . | / | . | / * 1298* | . | / | . | / * 1299* |. |/ |. |/ * 1300* + + + + * 1301* (000)------------(100) ( 0 )------------( 4 ) * 1302* * 1303* * 1304* Algorithm to find diagaonals * 1305* * 1306* Given a cube vertice d1 * 1307* then d2 = d1^(111) = d1^7 * 1308* * 1309* Where the cOR bit operator "^" is defined as follows: * 1310* 0^0 = 0 * 1311* 1^1 = 0 * 1312* 1^0 = 1 * 1313* 0^1 = 1 * 1314* * 1315* The four possible cube diagonals are * 1316* (000) --- (111) (0, 7) * 1317* (001) --- (110) <====> 2-tuple (1, 6) * 1318* (010) --- (101) rep. (2, 5) * 1319* (011) --- (100) (3, 4) * 1320* * 1321* Given a 2-tuple (d1,d2) that defines the diagonal of the cube, * 1322* six tetrahedrons are defined, e.g. * 1323* * 1324* (111) * 1325* . / . * 1326* . / . * 1327* . / . * 1328* . / . * 1329* . / . * 1330* . (101) . * 1331* . / | . <====> 4-tuple (0, 7, 4, 5) * 1332* . / | . rep. * 1333* . / | . * 1334* . / | . * 1335* . / | . * 1336* . / | . * 1337* . / | . * 1338* (000)------------(100) * 1339* * 1340* * 1341* Algorithm to find the six tetradedrons * 1342* * 1343* Given the diagonals vertices d1 and d2 such that d2=d1^7, the six * 1344* tetradedrons (six 4-tuples) can be found using the following * 1345* algorithm: * 1346* * 1347* shift(0) = (001) = 1 * 1348* shift(1) = (010) = 2 * 1349* shift(2) = (100) = 4 * 1350* tcount = 0 * 1351* For i=0,2 * 1352* For j=0,2 * 1353* c1 = d1^shift(i) * 1354* c2 = c1^shift(j) * 1355* If (c1 != d1) and (c1 != d2) and (c2!=d1) and (c2!=d2) Then * 1356* tetra(tcount) = (d1,d2,c1,c2) * 1357* tcount = tcount + 1 * 1358* End If * 1359* End For * 1360* End For * 1361* * 1362**************** Definitions of Cubes and Tetrahedrons ***************** 1363 1364 1365* ********************************************* 1366* * * 1367* * band_dos_generate * 1368* * * 1369* ********************************************* 1370 1371 subroutine band_dos_generate(unit,idx,idy,idz,eigs,neigs, 1372 > sign,npoints,emin,emax) 1373 implicit none 1374 integer unit 1375 integer idx,idy,idz 1376 real*8 eigs(idx,idy,idz,*) 1377 integer neigs 1378 real*8 sign 1379 integer npoints 1380 real*8 emin,emax 1381 1382* **** local variables **** 1383 integer dosgrid(3) 1384 integer i,j,k,ii,jj,kk,ncubes,ntetra,count 1385 integer ishft,jshft,kshft 1386 integer k1_d(4),k2_d(4),k3_d(4),k1_dd(4),k2_dd(4),k3_dd(4) 1387 integer id,d1(4),d2(4) 1388 integer itetra(4,6) 1389 real*8 VT,VG 1390 real*8 B(3,3),unitg(3,3),e,ecube(8),f,g,de 1391 real*8 k1,k2,k3,kx,ky,kz,kkx,kky,kkz,r,rmax 1392 1393* **** external functions **** 1394 real*8 lattice_unitg,Dstates_Cube,Nstates_Cube 1395 external lattice_unitg,Dstates_Cube,Nstates_Cube 1396 1397 dosgrid(1) = idx 1398 dosgrid(2) = idy 1399 dosgrid(3) = idz 1400 1401c write(unit,*) "dosgrid:",dosgrid 1402c write(unit,*) "neigs: ",neigs 1403c write(unit,*) "sign: ",sign 1404c write(unit,*) "npoints: ",npoints 1405c write(unit,*) "emin: ", emin 1406c write(unit,*) "emax: ", emax 1407 1408 do j=1,3 1409 do i=1,3 1410 B(i,j) = lattice_unitg(i,j) 1411 end do 1412 end do 1413 1414* **** volume of reciprocal unit cell, VG **** 1415 unitg(1,1) = B(2,2)*B(3,3) - B(3,2)*B(2,3) 1416 unitg(2,1) = B(3,2)*B(1,3) - B(1,2)*B(3,3) 1417 unitg(3,1) = B(1,2)*B(2,3) - B(2,2)*B(1,3) 1418 1419 unitg(1,2) = B(2,3)*B(3,1) - B(3,3)*B(2,1) 1420 unitg(2,2) = B(3,3)*B(1,1) - B(1,3)*B(3,1) 1421 unitg(3,2) = B(1,3)*B(2,1) - B(2,3)*B(1,1) 1422 1423 unitg(1,3) = B(2,1)*B(3,2) - B(3,1)*B(2,2) 1424 unitg(2,3) = B(3,1)*B(1,2) - B(1,1)*B(3,2) 1425 unitg(3,3) = B(1,1)*B(2,2) - B(2,1)*B(1,2) 1426 VG = B(1,1)*unitg(1,1) 1427 > + B(2,1)*unitg(2,1) 1428 > + B(3,1)*unitg(3,1) 1429 1430 ncubes = dosgrid(1)*dosgrid(2)*dosgrid(3) 1431 ntetra = ncubes*6 1432 VT = VG/dble(ntetra) 1433 1434c write(unit,*) "VG: ",VG 1435c write(unit,*) "number of cubes:",ncubes 1436c write(unit,*) "number of tetra:",ntetra 1437c write(unit,*) "VT: ",VT 1438c 1439c count = 0 1440c do k=0,dosgrid(3)-1 1441c do j=0,dosgrid(2)-1 1442c do i=0,dosgrid(1)-1 1443c count = count + 1 1444c k1 = (dble(i)/dble(dosgrid(1))) 1445c k2 = (dble(j)/dble(dosgrid(2))) 1446c k3 = (dble(k)/dble(dosgrid(3))) 1447c kx = k1*B(1,1) + k2*B(1,2) + k3*B(1,3) 1448c ky = k1*B(2,1) + k2*B(2,2) + k3*B(2,3) 1449c kz = k1*B(3,1) + k2*B(3,2) + k3*B(3,3) 1450c write(unit,*) i,j,k 1451c write(unit,3508) count,k1,k2,k3,kx,ky,kz 1452c write(unit,*) 1453c end do 1454c end do 1455c end do 1456 1457* ******************************** 1458* **** find shortest diagonal **** 1459* ******************************** 1460 1461* **** (000) ---- (111) **** 1462 k1_d(1) = 0 1463 k2_d(1) = 0 1464 k3_d(1) = 0 1465 k1_dd(1) = 1 1466 k2_dd(1) = 1 1467 k3_dd(1) = 1 1468 d1(1) = 0 1469 d2(1) = 7 1470 1471* **** (001) ---- (110) **** 1472 k1_d(2) = 1 1473 k2_d(2) = 0 1474 k3_d(2) = 0 1475 k1_dd(2) = 0 1476 k2_dd(2) = 1 1477 k3_dd(2) = 1 1478 d1(2) = 1 1479 d2(2) = 6 1480 1481* **** (010) ---- (101) **** 1482 k1_d(3) = 0 1483 k2_d(3) = 1 1484 k3_d(3) = 0 1485 k1_dd(3) = 1 1486 k2_dd(3) = 0 1487 k3_dd(3) = 1 1488 d1(3) = 2 1489 d2(3) = 5 1490 1491* **** (011) ---- (100) **** 1492 k1_d(4) = 1 1493 k2_d(4) = 1 1494 k3_d(4) = 0 1495 k1_dd(4) = 0 1496 k2_dd(4) = 0 1497 k3_dd(4) = 1 1498 d1(4) = 3 1499 d2(4) = 4 1500 1501 id = 1 1502 rmax = 9.99d9 1503 do i=1,4 1504 kx = k1_d(i)*B(1,1) + k2_d(i)*B(1,2) + k3_d(i)*B(1,3) 1505 ky = k1_d(i)*B(2,1) + k2_d(i)*B(2,2) + k3_d(i)*B(2,3) 1506 kz = k1_d(i)*B(3,1) + k2_d(i)*B(3,2) + k3_d(i)*B(3,3) 1507 1508 kkx = k1_dd(i)*B(1,1) + k2_dd(i)*B(1,2) + k3_dd(i)*B(1,3) 1509 kky = k1_dd(i)*B(2,1) + k2_dd(i)*B(2,2) + k3_dd(i)*B(2,3) 1510 kkz = k1_dd(i)*B(3,1) + k2_dd(i)*B(3,2) + k3_dd(i)*B(3,3) 1511 r = (kx-kkx)**2 + (ky-kky)**2 + (kz-kkz)**2 1512 !write(unit,*) "diagonal distance:",i,r,rmax 1513 if (r.lt.rmax) then 1514 rmax = r 1515 id = i 1516 end if 1517 end do 1518 1519 !write(unit,*) "diagonal d1,d2 =",d1(id),d2(id) 1520 1521* **** define six tetradrons - clunky but don't know defn of cOR in fortran **** 1522 if (id.eq.1) then 1523 itetra(1,1) = 0 +1 1524 itetra(2,1) = 7 +1 1525 itetra(3,1) = 1 +1 1526 itetra(4,1) = 3 +1 1527 1528 itetra(1,2) = 0 +1 1529 itetra(2,2) = 7 +1 1530 itetra(3,2) = 1 +1 1531 itetra(4,2) = 5 +1 1532 1533 itetra(1,3) = 0 +1 1534 itetra(2,3) = 7 +1 1535 itetra(3,3) = 2 +1 1536 itetra(4,3) = 3 +1 1537 1538 itetra(1,4) = 0 +1 1539 itetra(2,4) = 7 +1 1540 itetra(3,4) = 2 +1 1541 itetra(4,4) = 6 +1 1542 1543 itetra(1,5) = 0 +1 1544 itetra(2,5) = 7 +1 1545 itetra(3,5) = 4 +1 1546 itetra(4,5) = 5 +1 1547 1548 itetra(1,6) = 0 +1 1549 itetra(2,6) = 7 +1 1550 itetra(3,6) = 4 +1 1551 itetra(4,6) = 6 +1 1552 else if (id.eq.2) then 1553 itetra(1,1) = 1 +1 1554 itetra(2,1) = 6 +1 1555 itetra(3,1) = 0 +1 1556 itetra(4,1) = 2 +1 1557 1558 itetra(1,2) = 1 +1 1559 itetra(2,2) = 6 +1 1560 itetra(3,2) = 0 +1 1561 itetra(4,2) = 4 +1 1562 1563 itetra(1,3) = 1 +1 1564 itetra(2,3) = 6 +1 1565 itetra(3,3) = 3 +1 1566 itetra(4,3) = 2 +1 1567 1568 itetra(1,4) = 1 +1 1569 itetra(2,4) = 6 +1 1570 itetra(3,4) = 3 +1 1571 itetra(4,4) = 7 +1 1572 1573 itetra(1,5) = 1 +1 1574 itetra(2,5) = 6 +1 1575 itetra(3,5) = 5 +1 1576 itetra(4,5) = 4 +1 1577 1578 itetra(1,6) = 1 +1 1579 itetra(2,6) = 6 +1 1580 itetra(3,6) = 5 +1 1581 itetra(4,6) = 7 +1 1582 else if (id.eq.3) then 1583 itetra(1,1) = 2 +1 1584 itetra(2,1) = 5 +1 1585 itetra(3,1) = 3 +1 1586 itetra(4,1) = 1 +1 1587 1588 itetra(1,2) = 2 +1 1589 itetra(2,2) = 5 +1 1590 itetra(3,2) = 3 +1 1591 itetra(4,2) = 7 +1 1592 1593 itetra(1,3) = 2 +1 1594 itetra(2,3) = 5 +1 1595 itetra(3,3) = 0 +1 1596 itetra(4,3) = 1 +1 1597 1598 itetra(1,4) = 2 +1 1599 itetra(2,4) = 5 +1 1600 itetra(3,4) = 0 +1 1601 itetra(4,4) = 4 +1 1602 1603 itetra(1,5) = 2 +1 1604 itetra(2,5) = 5 +1 1605 itetra(3,5) = 6 +1 1606 itetra(4,5) = 7 +1 1607 1608 itetra(1,6) = 2 +1 1609 itetra(2,6) = 5 +1 1610 itetra(3,6) = 6 +1 1611 itetra(4,6) = 4 +1 1612 else if (id.eq.4) then 1613 itetra(1,1) = 3 +1 1614 itetra(2,1) = 4 +1 1615 itetra(3,1) = 2 +1 1616 itetra(4,1) = 0 +1 1617 1618 itetra(1,2) = 3 +1 1619 itetra(2,2) = 4 +1 1620 itetra(3,2) = 2 +1 1621 itetra(4,2) = 6 +1 1622 1623 itetra(1,3) = 3 +1 1624 itetra(2,3) = 4 +1 1625 itetra(3,3) = 1 +1 1626 itetra(4,3) = 0 +1 1627 1628 itetra(1,4) = 3 +1 1629 itetra(2,4) = 4 +1 1630 itetra(3,4) = 1 +1 1631 itetra(4,4) = 5 +1 1632 1633 itetra(1,5) = 3 +1 1634 itetra(2,5) = 4 +1 1635 itetra(3,5) = 7 +1 1636 itetra(4,5) = 6 +1 1637 1638 itetra(1,6) = 3 +1 1639 itetra(2,6) = 4 +1 1640 itetra(3,6) = 7 +1 1641 itetra(4,6) = 5 +1 1642 end if 1643 1644 1645c do i=1,6 1646c write(unit,*) id,"tetra :",i,"(",(itetra(j,i),j=1,4),")" 1647c end do 1648 1649 de = (emax-emin)/dble(npoints-1) 1650 do k=1,npoints 1651 e = emin + (k-1)*de 1652 1653 f = 0.0d0 1654 g = 0.0d0 1655 do kk=1,dosgrid(3) 1656 do jj=1,dosgrid(2) 1657 do ii=1,dosgrid(1) 1658 ishft = ii+1 1659 jshft = jj+1 1660 kshft = kk+1 1661 if (ishft.gt.dosgrid(1)) ishft=1 1662 if (jshft.gt.dosgrid(2)) jshft=1 1663 if (kshft.gt.dosgrid(3)) kshft=1 1664 do i=1,neigs 1665 ecube(1) = eigs(ii, jj, kk, i) ! (000) 1666 ecube(2) = eigs(ishft, jj, kk, i) ! (001) 1667 ecube(3) = eigs(ii, jshft, kk, i) ! (010) 1668 ecube(4) = eigs(ishft, jshft, kk, i) ! (011) 1669 ecube(5) = eigs(ii, jj, kshft, i) ! (100) 1670 ecube(6) = eigs(ishft, jj, kshft, i) ! (101) 1671 ecube(7) = eigs( ii, jshft, kshft, i) ! (110) 1672 ecube(8) = eigs(ishft, jshft, kshft, i) ! (111) 1673 1674 f = f + Dstates_Cube(e,itetra,ecube) 1675 g = g + Nstates_Cube(e,itetra,ecube) 1676 end do 1677 end do 1678 end do 1679 end do 1680 f = f*(VT/VG) 1681 g = g*(VT/VG) 1682 1683 write(unit,1310) e,f*sign,g*sign 1684 end do 1685 1686 1687 return 1688 1310 FORMAT(3E15.5) 1689 3508 FORMAT(/' Brillouin zone point: ',i5, 1690 > /' k =<',3f8.3,'> . <b1,b2,b3> ', 1691 > /' =<',3f8.3,'>') 1692 end 1693 1694 1695 real*8 function Dstates_Cube(e,itetra,ecube) 1696 implicit none 1697 real*8 e 1698 integer itetra(4,6) 1699 real*8 ecube(8) 1700 1701* **** local variables **** 1702 integer i,j,k 1703 real*8 ds,etetra(4),swap 1704 1705 real*8 Dstates_Tetra 1706 external Dstates_Tetra 1707 1708* **** sum over 6 tetrahedrons **** 1709 ds = 0.0d0 1710 do k=1,6 1711 etetra(1) = ecube(itetra(1,k)) 1712 etetra(2) = ecube(itetra(2,k)) 1713 etetra(3) = ecube(itetra(3,k)) 1714 etetra(4) = ecube(itetra(4,k)) 1715 1716* **** bubble sort **** 1717 do j=1,3 1718 do i=j+1,4 1719 if (etetra(j).gt.etetra(i)) then 1720 swap = etetra(i) 1721 etetra(i) = etetra(j) 1722 etetra(j) = swap 1723 end if 1724 end do 1725 end do 1726 ds = ds + Dstates_Tetra(e,etetra) 1727 end do 1728 1729 Dstates_cube = ds 1730 return 1731 end 1732 1733 1734 real*8 function Dstates_Tetra(e,ee) 1735 implicit none 1736 real*8 e 1737 real*8 ee(4) 1738 1739* **** local variables **** 1740 real*8 ds 1741 real*8 e1,e2,e4 1742 real*8 e21,e31,e41,e32,e42,e43 1743 1744 if ((ee(1).le.e).and.(e.lt.ee(2))) then 1745 e1 = e-ee(1) 1746 e21 = ee(2) - ee(1) 1747 e31 = ee(3) - ee(1) 1748 e41 = ee(4) - ee(1) 1749 ds = 3.0d0*e1*e1/(e21*e31*e41) 1750 else if ((ee(2).le.e).and.(e.lt.ee(3))) then 1751 e2 = e-ee(2) 1752 e21 = ee(2) - ee(1) 1753 e31 = ee(3) - ee(1) 1754 e41 = ee(4) - ee(1) 1755 e32 = ee(3) - ee(2) 1756 e42 = ee(4) - ee(2) 1757 ds = (3.0d0*e21+6.0d0*e2-3.0d0*(e31+e42)*e2*e2/(e32*e42)) 1758 > /(e31*e41) 1759 else if ((ee(3).le.e).and.(e.lt.ee(4))) then 1760 e4 = ee(4)-e 1761 e41 = ee(4) - ee(1) 1762 e42 = ee(4) - ee(2) 1763 e43 = ee(4) - ee(3) 1764 ds = 3.0d0*e4*e4/(e41*e42*e43) 1765 else 1766 ds = 0.0d0 1767 end if 1768 1769 1770 Dstates_Tetra = ds 1771 return 1772 end 1773 1774 1775 real*8 function Nstates_Cube(e,itetra,ecube) 1776 implicit none 1777 real*8 e 1778 integer itetra(4,6) 1779 real*8 ecube(8) 1780 1781* **** local variables **** 1782 integer i,j,k 1783 real*8 ds,etetra(4),swap 1784 1785 real*8 Nstates_Tetra 1786 external Nstates_Tetra 1787 1788* **** sum over 6 tetrahedrons **** 1789 ds = 0.0d0 1790 do k=1,6 1791 etetra(1) = ecube(itetra(1,k)) 1792 etetra(2) = ecube(itetra(2,k)) 1793 etetra(3) = ecube(itetra(3,k)) 1794 etetra(4) = ecube(itetra(4,k)) 1795 1796* **** bubble sort **** 1797 do j=1,3 1798 do i=j+1,4 1799 if (etetra(j).gt.etetra(i)) then 1800 swap = etetra(i) 1801 etetra(i) = etetra(j) 1802 etetra(j) = swap 1803 end if 1804 end do 1805 end do 1806 ds = ds + Nstates_Tetra(e,etetra) 1807 end do 1808 1809 Nstates_cube = ds 1810 return 1811 end 1812 1813 1814 real*8 function Nstates_Tetra(e,ee) 1815 implicit none 1816 real*8 e 1817 real*8 ee(4) 1818 1819* **** local variables **** 1820 real*8 ds 1821 real*8 e1,e2,e4 1822 real*8 e21,e31,e41,e32,e42,e43 1823 1824 if ((ee(1).le.e).and.(e.lt.ee(2))) then 1825 e1 = e-ee(1) 1826 e21 = ee(2) - ee(1) 1827 e31 = ee(3) - ee(1) 1828 e41 = ee(4) - ee(1) 1829 ds = e1*e1*e1/(e21*e31*e41) 1830 else if ((ee(2).le.e).and.(e.lt.ee(3))) then 1831 e2 = e-ee(2) 1832 e21 = ee(2) - ee(1) 1833 e31 = ee(3) - ee(1) 1834 e41 = ee(4) - ee(1) 1835 e32 = ee(3) - ee(2) 1836 e42 = ee(4) - ee(2) 1837 ds = (e21*e21 1838 > + 3.0d0*e21*e2 1839 > + 3.0d0*e2*e2 1840 > - (e31+e42)*e2*e2*e2/(e32*e42)) 1841 > /(e31*e41) 1842 else if ((ee(3).le.e).and.(e.lt.ee(4))) then 1843 e4 = ee(4)-e 1844 e41 = ee(4) - ee(1) 1845 e42 = ee(4) - ee(2) 1846 e43 = ee(4) - ee(3) 1847 ds = 1.0d0 - e4*e4*e4/(e41*e42*e43) 1848 else if (e.ge.ee(4)) then 1849 ds = 1.0d0 1850 else 1851 ds = 0.0d0 1852 end if 1853 1854 1855 Nstates_Tetra = ds 1856 return 1857 end 1858 1859 1860 1861* ********************************************* 1862* * * 1863* * band_dos_weights_generate * 1864* * * 1865* ********************************************* 1866 1867 subroutine band_dos_weights_generate(idx,idy,idz, 1868 > eigs,neigs, 1869 > ein,weight) 1870 implicit none 1871 integer idx,idy,idz 1872 real*8 eigs(idx,idy,idz,*) 1873 integer neigs 1874 real*8 ein 1875 real*8 weight(neigs,idx,idy,idz) 1876 1877* **** local variables **** 1878 integer dosgrid(3) 1879 integer i,j,k,ii,jj,kk,ncubes,ntetra,count 1880 integer ishft,jshft,kshft 1881 integer k1_d(4),k2_d(4),k3_d(4),k1_dd(4),k2_dd(4),k3_dd(4) 1882 integer id,d1(4),d2(4) 1883 integer itetra(4,6) 1884 real*8 VT,VG 1885 real*8 B(3,3),unitg(3,3),e,ecube(8),wcube(8) 1886 real*8 k1,k2,k3,kx,ky,kz,kkx,kky,kkz,r,rmax 1887 1888* **** external functions **** 1889 real*8 lattice_unitg 1890 external lattice_unitg 1891 1892 dosgrid(1) = idx 1893 dosgrid(2) = idy 1894 dosgrid(3) = idz 1895 1896 1897 do j=1,3 1898 do i=1,3 1899 B(i,j) = lattice_unitg(i,j) 1900 end do 1901 end do 1902 1903* **** volume of reciprocal unit cell, VG **** 1904 unitg(1,1) = B(2,2)*B(3,3) - B(3,2)*B(2,3) 1905 unitg(2,1) = B(3,2)*B(1,3) - B(1,2)*B(3,3) 1906 unitg(3,1) = B(1,2)*B(2,3) - B(2,2)*B(1,3) 1907 1908 unitg(1,2) = B(2,3)*B(3,1) - B(3,3)*B(2,1) 1909 unitg(2,2) = B(3,3)*B(1,1) - B(1,3)*B(3,1) 1910 unitg(3,2) = B(1,3)*B(2,1) - B(2,3)*B(1,1) 1911 1912 unitg(1,3) = B(2,1)*B(3,2) - B(3,1)*B(2,2) 1913 unitg(2,3) = B(3,1)*B(1,2) - B(1,1)*B(3,2) 1914 unitg(3,3) = B(1,1)*B(2,2) - B(2,1)*B(1,2) 1915 VG = B(1,1)*unitg(1,1) 1916 > + B(2,1)*unitg(2,1) 1917 > + B(3,1)*unitg(3,1) 1918 1919 ncubes = dosgrid(1)*dosgrid(2)*dosgrid(3) 1920 ntetra = ncubes*6 1921 VT = VG/dble(ntetra) 1922 1923 1924* ******************************** 1925* **** find shortest diagonal **** 1926* ******************************** 1927 1928* **** (000) ---- (111) **** 1929 k1_d(1) = 0 1930 k2_d(1) = 0 1931 k3_d(1) = 0 1932 k1_dd(1) = 1 1933 k2_dd(1) = 1 1934 k3_dd(1) = 1 1935 d1(1) = 0 1936 d2(1) = 7 1937 1938* **** (001) ---- (110) **** 1939 k1_d(2) = 1 1940 k2_d(2) = 0 1941 k3_d(2) = 0 1942 k1_dd(2) = 0 1943 k2_dd(2) = 1 1944 k3_dd(2) = 1 1945 d1(2) = 1 1946 d2(2) = 6 1947 1948* **** (010) ---- (101) **** 1949 k1_d(3) = 0 1950 k2_d(3) = 1 1951 k3_d(3) = 0 1952 k1_dd(3) = 1 1953 k2_dd(3) = 0 1954 k3_dd(3) = 1 1955 d1(3) = 2 1956 d2(3) = 5 1957 1958* **** (011) ---- (100) **** 1959 k1_d(4) = 1 1960 k2_d(4) = 1 1961 k3_d(4) = 0 1962 k1_dd(4) = 0 1963 k2_dd(4) = 0 1964 k3_dd(4) = 1 1965 d1(4) = 3 1966 d2(4) = 4 1967 1968 id = 1 1969 rmax = 9.99d9 1970 do i=1,4 1971 kx = k1_d(i)*B(1,1) + k2_d(i)*B(1,2) + k3_d(i)*B(1,3) 1972 ky = k1_d(i)*B(2,1) + k2_d(i)*B(2,2) + k3_d(i)*B(2,3) 1973 kz = k1_d(i)*B(3,1) + k2_d(i)*B(3,2) + k3_d(i)*B(3,3) 1974 1975 kkx = k1_dd(i)*B(1,1) + k2_dd(i)*B(1,2) + k3_dd(i)*B(1,3) 1976 kky = k1_dd(i)*B(2,1) + k2_dd(i)*B(2,2) + k3_dd(i)*B(2,3) 1977 kkz = k1_dd(i)*B(3,1) + k2_dd(i)*B(3,2) + k3_dd(i)*B(3,3) 1978 r = (kx-kkx)**2 + (ky-kky)**2 + (kz-kkz)**2 1979 !write(unit,*) "diagonal distance:",i,r,rmax 1980 if (r.lt.rmax) then 1981 rmax = r 1982 id = i 1983 end if 1984 end do 1985 1986 1987* **** define six tetradrons - clunky but don't know defn of cOR in fortran **** 1988 if (id.eq.1) then 1989 itetra(1,1) = 0 +1 1990 itetra(2,1) = 7 +1 1991 itetra(3,1) = 1 +1 1992 itetra(4,1) = 3 +1 1993 1994 itetra(1,2) = 0 +1 1995 itetra(2,2) = 7 +1 1996 itetra(3,2) = 1 +1 1997 itetra(4,2) = 5 +1 1998 1999 itetra(1,3) = 0 +1 2000 itetra(2,3) = 7 +1 2001 itetra(3,3) = 2 +1 2002 itetra(4,3) = 3 +1 2003 2004 itetra(1,4) = 0 +1 2005 itetra(2,4) = 7 +1 2006 itetra(3,4) = 2 +1 2007 itetra(4,4) = 6 +1 2008 2009 itetra(1,5) = 0 +1 2010 itetra(2,5) = 7 +1 2011 itetra(3,5) = 4 +1 2012 itetra(4,5) = 5 +1 2013 2014 itetra(1,6) = 0 +1 2015 itetra(2,6) = 7 +1 2016 itetra(3,6) = 4 +1 2017 itetra(4,6) = 6 +1 2018 else if (id.eq.2) then 2019 itetra(1,1) = 1 +1 2020 itetra(2,1) = 6 +1 2021 itetra(3,1) = 0 +1 2022 itetra(4,1) = 2 +1 2023 2024 itetra(1,2) = 1 +1 2025 itetra(2,2) = 6 +1 2026 itetra(3,2) = 0 +1 2027 itetra(4,2) = 4 +1 2028 2029 itetra(1,3) = 1 +1 2030 itetra(2,3) = 6 +1 2031 itetra(3,3) = 3 +1 2032 itetra(4,3) = 2 +1 2033 2034 itetra(1,4) = 1 +1 2035 itetra(2,4) = 6 +1 2036 itetra(3,4) = 3 +1 2037 itetra(4,4) = 7 +1 2038 2039 itetra(1,5) = 1 +1 2040 itetra(2,5) = 6 +1 2041 itetra(3,5) = 5 +1 2042 itetra(4,5) = 4 +1 2043 2044 itetra(1,6) = 1 +1 2045 itetra(2,6) = 6 +1 2046 itetra(3,6) = 5 +1 2047 itetra(4,6) = 7 +1 2048 else if (id.eq.3) then 2049 itetra(1,1) = 2 +1 2050 itetra(2,1) = 5 +1 2051 itetra(3,1) = 3 +1 2052 itetra(4,1) = 1 +1 2053 2054 itetra(1,2) = 2 +1 2055 itetra(2,2) = 5 +1 2056 itetra(3,2) = 3 +1 2057 itetra(4,2) = 7 +1 2058 2059 itetra(1,3) = 2 +1 2060 itetra(2,3) = 5 +1 2061 itetra(3,3) = 0 +1 2062 itetra(4,3) = 1 +1 2063 2064 itetra(1,4) = 2 +1 2065 itetra(2,4) = 5 +1 2066 itetra(3,4) = 0 +1 2067 itetra(4,4) = 4 +1 2068 2069 itetra(1,5) = 2 +1 2070 itetra(2,5) = 5 +1 2071 itetra(3,5) = 6 +1 2072 itetra(4,5) = 7 +1 2073 2074 itetra(1,6) = 2 +1 2075 itetra(2,6) = 5 +1 2076 itetra(3,6) = 6 +1 2077 itetra(4,6) = 4 +1 2078 else if (id.eq.4) then 2079 itetra(1,1) = 3 +1 2080 itetra(2,1) = 4 +1 2081 itetra(3,1) = 2 +1 2082 itetra(4,1) = 0 +1 2083 2084 itetra(1,2) = 3 +1 2085 itetra(2,2) = 4 +1 2086 itetra(3,2) = 2 +1 2087 itetra(4,2) = 6 +1 2088 2089 itetra(1,3) = 3 +1 2090 itetra(2,3) = 4 +1 2091 itetra(3,3) = 1 +1 2092 itetra(4,3) = 0 +1 2093 2094 itetra(1,4) = 3 +1 2095 itetra(2,4) = 4 +1 2096 itetra(3,4) = 1 +1 2097 itetra(4,4) = 5 +1 2098 2099 itetra(1,5) = 3 +1 2100 itetra(2,5) = 4 +1 2101 itetra(3,5) = 7 +1 2102 itetra(4,5) = 6 +1 2103 2104 itetra(1,6) = 3 +1 2105 itetra(2,6) = 4 +1 2106 itetra(3,6) = 7 +1 2107 itetra(4,6) = 5 +1 2108 end if 2109 2110 2111 2112 2113 e = ein 2114 2115 call dcopy(dosgrid(1)*dosgrid(2)*dosgrid(3)*neigs, 2116 > 0.0d0,0,weight,1) 2117 do kk=1,dosgrid(3) 2118 do jj=1,dosgrid(2) 2119 do ii=1,dosgrid(1) 2120 ishft = ii+1 2121 jshft = jj+1 2122 kshft = kk+1 2123 if (ishft.gt.dosgrid(1)) ishft=1 2124 if (jshft.gt.dosgrid(2)) jshft=1 2125 if (kshft.gt.dosgrid(3)) kshft=1 2126 do i=1,neigs 2127 ecube(1) = eigs(ii, jj, kk, i) ! (000) 2128 ecube(2) = eigs(ishft, jj, kk, i) ! (001) 2129 ecube(3) = eigs(ii, jshft, kk, i) ! (010) 2130 ecube(4) = eigs(ishft, jshft, kk, i) ! (011) 2131 ecube(5) = eigs(ii, jj, kshft, i) ! (100) 2132 ecube(6) = eigs(ishft, jj, kshft, i) ! (101) 2133 ecube(7) = eigs( ii, jshft, kshft, i) ! (110) 2134 ecube(8) = eigs(ishft, jshft, kshft, i) ! (111) 2135 2136 call Dstates_weight_Cube(e,itetra,ecube,wcube) 2137 weight(i,ii, jj, kk) = weight(i,ii, jj, kk) 2138 > + wcube(1) 2139 weight(i,ishft, jj, kk) = weight(i,ishft, jj, kk) 2140 > + wcube(2) 2141 weight(i,ii, jshft, kk) = weight(i,ii, jshft, kk) 2142 > + wcube(3) 2143 weight(i,ishft,jshft, kk) = weight(i,ishft,jshft, kk) 2144 > + wcube(4) 2145 weight(i,ii, jj,kshft) = weight(i,ii, jj,kshft) 2146 > + wcube(5) 2147 weight(i,ishft, jj,kshft) = weight(i,ishft, jj,kshft) 2148 > + wcube(6) 2149 weight(i, ii,jshft,kshft) = weight(i, ii,jshft,kshft) 2150 > + wcube(7) 2151 weight(i,ishft,jshft,kshft) = weight(i,ishft,jshft,kshft) 2152 > + wcube(8) 2153 end do 2154 end do 2155 end do 2156 end do 2157 call dscal(dosgrid(1)*dosgrid(2)*dosgrid(3)*neigs, 2158 > (VT/VG),weight,1) 2159 2160 2161 return 2162 end 2163 2164 2165 subroutine Dstates_weight_Cube(e,itetra,ecube,wcube) 2166 implicit none 2167 real*8 e 2168 integer itetra(4,6) 2169 real*8 ecube(8) 2170 real*8 wcube(8) 2171 2172* **** local variables **** 2173 integer i,j,k 2174 real*8 ds,etetra(4),swap 2175 2176* **** external functions **** 2177 real*8 Dstates_Tetra 2178 external Dstates_Tetra 2179 2180 2181* **** sum over 6 tetrahedrons **** 2182 call dcopy(8,0.0d0,0,wcube,1) 2183 do k=1,6 2184 etetra(1) = ecube(itetra(1,k)) 2185 etetra(2) = ecube(itetra(2,k)) 2186 etetra(3) = ecube(itetra(3,k)) 2187 etetra(4) = ecube(itetra(4,k)) 2188 2189* **** bubble sort **** 2190 do j=1,3 2191 do i=j+1,4 2192 if (etetra(j).gt.etetra(i)) then 2193 swap = etetra(i) 2194 etetra(i) = etetra(j) 2195 etetra(j) = swap 2196 end if 2197 end do 2198 end do 2199 2200 ds = Dstates_Tetra(e,etetra) 2201 2202 wcube(itetra(1,k)) = wcube(itetra(1,k)) + 0.25d0*ds 2203 wcube(itetra(2,k)) = wcube(itetra(2,k)) + 0.25d0*ds 2204 wcube(itetra(3,k)) = wcube(itetra(3,k)) + 0.25d0*ds 2205 wcube(itetra(4,k)) = wcube(itetra(4,k)) + 0.25d0*ds 2206 end do 2207 2208 return 2209 end 2210 2211 2212 2213* ********************************************* 2214* * * 2215* * band_projected_dos_generate * 2216* * * 2217* ********************************************* 2218 2219 subroutine band_projected_dos_generate( 2220 > unit,idx,idy,idz,eigs,pweight,neigs, 2221 > sign,npoints,emin,emax) 2222 implicit none 2223 integer unit 2224 integer idx,idy,idz 2225 real*8 eigs(idx,idy,idz,*) 2226 real*8 pweight(idx,idy,idz,*) 2227 integer neigs 2228 real*8 sign 2229 integer npoints 2230 real*8 emin,emax 2231 2232* **** local variables **** 2233 integer dosgrid(3) 2234 integer i,j,k,ii,jj,kk,ncubes,ntetra,count 2235 integer ishft,jshft,kshft 2236 integer k1_d(4),k2_d(4),k3_d(4),k1_dd(4),k2_dd(4),k3_dd(4) 2237 integer id,d1(4),d2(4) 2238 integer itetra(4,6) 2239 real*8 VT,VG 2240 real*8 B(3,3),unitg(3,3),e,ecube(8),pcube(8),f,g,de,ff 2241 real*8 k1,k2,k3,kx,ky,kz,kkx,kky,kkz,r,rmax 2242 2243* **** external functions **** 2244 real*8 lattice_unitg 2245 real*8 Dstates_Cube_projected,Nstates_Cube_projected 2246 external lattice_unitg 2247 external Dstates_Cube_projected,Nstates_Cube_projected 2248 2249 dosgrid(1) = idx 2250 dosgrid(2) = idy 2251 dosgrid(3) = idz 2252 2253 do j=1,3 2254 do i=1,3 2255 B(i,j) = lattice_unitg(i,j) 2256 end do 2257 end do 2258 2259* **** volume of reciprocal unit cell, VG **** 2260 unitg(1,1) = B(2,2)*B(3,3) - B(3,2)*B(2,3) 2261 unitg(2,1) = B(3,2)*B(1,3) - B(1,2)*B(3,3) 2262 unitg(3,1) = B(1,2)*B(2,3) - B(2,2)*B(1,3) 2263 2264 unitg(1,2) = B(2,3)*B(3,1) - B(3,3)*B(2,1) 2265 unitg(2,2) = B(3,3)*B(1,1) - B(1,3)*B(3,1) 2266 unitg(3,2) = B(1,3)*B(2,1) - B(2,3)*B(1,1) 2267 2268 unitg(1,3) = B(2,1)*B(3,2) - B(3,1)*B(2,2) 2269 unitg(2,3) = B(3,1)*B(1,2) - B(1,1)*B(3,2) 2270 unitg(3,3) = B(1,1)*B(2,2) - B(2,1)*B(1,2) 2271 VG = B(1,1)*unitg(1,1) 2272 > + B(2,1)*unitg(2,1) 2273 > + B(3,1)*unitg(3,1) 2274 2275 ncubes = dosgrid(1)*dosgrid(2)*dosgrid(3) 2276 ntetra = ncubes*6 2277 VT = VG/dble(ntetra) 2278 2279 2280* ******************************** 2281* **** find shortest diagonal **** 2282* ******************************** 2283 2284* **** (000) ---- (111) **** 2285 k1_d(1) = 0 2286 k2_d(1) = 0 2287 k3_d(1) = 0 2288 k1_dd(1) = 1 2289 k2_dd(1) = 1 2290 k3_dd(1) = 1 2291 d1(1) = 0 2292 d2(1) = 7 2293 2294* **** (001) ---- (110) **** 2295 k1_d(2) = 1 2296 k2_d(2) = 0 2297 k3_d(2) = 0 2298 k1_dd(2) = 0 2299 k2_dd(2) = 1 2300 k3_dd(2) = 1 2301 d1(2) = 1 2302 d2(2) = 6 2303 2304* **** (010) ---- (101) **** 2305 k1_d(3) = 0 2306 k2_d(3) = 1 2307 k3_d(3) = 0 2308 k1_dd(3) = 1 2309 k2_dd(3) = 0 2310 k3_dd(3) = 1 2311 d1(3) = 2 2312 d2(3) = 5 2313 2314* **** (011) ---- (100) **** 2315 k1_d(4) = 1 2316 k2_d(4) = 1 2317 k3_d(4) = 0 2318 k1_dd(4) = 0 2319 k2_dd(4) = 0 2320 k3_dd(4) = 1 2321 d1(4) = 3 2322 d2(4) = 4 2323 2324 id = 1 2325 rmax = 9.99d9 2326 do i=1,4 2327 kx = k1_d(i)*B(1,1) + k2_d(i)*B(1,2) + k3_d(i)*B(1,3) 2328 ky = k1_d(i)*B(2,1) + k2_d(i)*B(2,2) + k3_d(i)*B(2,3) 2329 kz = k1_d(i)*B(3,1) + k2_d(i)*B(3,2) + k3_d(i)*B(3,3) 2330 2331 kkx = k1_dd(i)*B(1,1) + k2_dd(i)*B(1,2) + k3_dd(i)*B(1,3) 2332 kky = k1_dd(i)*B(2,1) + k2_dd(i)*B(2,2) + k3_dd(i)*B(2,3) 2333 kkz = k1_dd(i)*B(3,1) + k2_dd(i)*B(3,2) + k3_dd(i)*B(3,3) 2334 r = (kx-kkx)**2 + (ky-kky)**2 + (kz-kkz)**2 2335 !write(unit,*) "diagonal distance:",i,r,rmax 2336 if (r.lt.rmax) then 2337 rmax = r 2338 id = i 2339 end if 2340 end do 2341 2342 !write(unit,*) "diagonal d1,d2 =",d1(id),d2(id) 2343 2344* **** define six tetradrons - clunky but don't know defn of cOR in fortran **** 2345 if (id.eq.1) then 2346 itetra(1,1) = 0 +1 2347 itetra(2,1) = 7 +1 2348 itetra(3,1) = 1 +1 2349 itetra(4,1) = 3 +1 2350 2351 itetra(1,2) = 0 +1 2352 itetra(2,2) = 7 +1 2353 itetra(3,2) = 1 +1 2354 itetra(4,2) = 5 +1 2355 2356 itetra(1,3) = 0 +1 2357 itetra(2,3) = 7 +1 2358 itetra(3,3) = 2 +1 2359 itetra(4,3) = 3 +1 2360 2361 itetra(1,4) = 0 +1 2362 itetra(2,4) = 7 +1 2363 itetra(3,4) = 2 +1 2364 itetra(4,4) = 6 +1 2365 2366 itetra(1,5) = 0 +1 2367 itetra(2,5) = 7 +1 2368 itetra(3,5) = 4 +1 2369 itetra(4,5) = 5 +1 2370 2371 itetra(1,6) = 0 +1 2372 itetra(2,6) = 7 +1 2373 itetra(3,6) = 4 +1 2374 itetra(4,6) = 6 +1 2375 else if (id.eq.2) then 2376 itetra(1,1) = 1 +1 2377 itetra(2,1) = 6 +1 2378 itetra(3,1) = 0 +1 2379 itetra(4,1) = 2 +1 2380 2381 itetra(1,2) = 1 +1 2382 itetra(2,2) = 6 +1 2383 itetra(3,2) = 0 +1 2384 itetra(4,2) = 4 +1 2385 2386 itetra(1,3) = 1 +1 2387 itetra(2,3) = 6 +1 2388 itetra(3,3) = 3 +1 2389 itetra(4,3) = 2 +1 2390 2391 itetra(1,4) = 1 +1 2392 itetra(2,4) = 6 +1 2393 itetra(3,4) = 3 +1 2394 itetra(4,4) = 7 +1 2395 2396 itetra(1,5) = 1 +1 2397 itetra(2,5) = 6 +1 2398 itetra(3,5) = 5 +1 2399 itetra(4,5) = 4 +1 2400 2401 itetra(1,6) = 1 +1 2402 itetra(2,6) = 6 +1 2403 itetra(3,6) = 5 +1 2404 itetra(4,6) = 7 +1 2405 else if (id.eq.3) then 2406 itetra(1,1) = 2 +1 2407 itetra(2,1) = 5 +1 2408 itetra(3,1) = 3 +1 2409 itetra(4,1) = 1 +1 2410 2411 itetra(1,2) = 2 +1 2412 itetra(2,2) = 5 +1 2413 itetra(3,2) = 3 +1 2414 itetra(4,2) = 7 +1 2415 2416 itetra(1,3) = 2 +1 2417 itetra(2,3) = 5 +1 2418 itetra(3,3) = 0 +1 2419 itetra(4,3) = 1 +1 2420 2421 itetra(1,4) = 2 +1 2422 itetra(2,4) = 5 +1 2423 itetra(3,4) = 0 +1 2424 itetra(4,4) = 4 +1 2425 2426 itetra(1,5) = 2 +1 2427 itetra(2,5) = 5 +1 2428 itetra(3,5) = 6 +1 2429 itetra(4,5) = 7 +1 2430 2431 itetra(1,6) = 2 +1 2432 itetra(2,6) = 5 +1 2433 itetra(3,6) = 6 +1 2434 itetra(4,6) = 4 +1 2435 else if (id.eq.4) then 2436 itetra(1,1) = 3 +1 2437 itetra(2,1) = 4 +1 2438 itetra(3,1) = 2 +1 2439 itetra(4,1) = 0 +1 2440 2441 itetra(1,2) = 3 +1 2442 itetra(2,2) = 4 +1 2443 itetra(3,2) = 2 +1 2444 itetra(4,2) = 6 +1 2445 2446 itetra(1,3) = 3 +1 2447 itetra(2,3) = 4 +1 2448 itetra(3,3) = 1 +1 2449 itetra(4,3) = 0 +1 2450 2451 itetra(1,4) = 3 +1 2452 itetra(2,4) = 4 +1 2453 itetra(3,4) = 1 +1 2454 itetra(4,4) = 5 +1 2455 2456 itetra(1,5) = 3 +1 2457 itetra(2,5) = 4 +1 2458 itetra(3,5) = 7 +1 2459 itetra(4,5) = 6 +1 2460 2461 itetra(1,6) = 3 +1 2462 itetra(2,6) = 4 +1 2463 itetra(3,6) = 7 +1 2464 itetra(4,6) = 5 +1 2465 end if 2466 2467 2468c do i=1,6 2469c write(unit,*) id,"tetra :",i,"(",(itetra(j,i),j=1,4),")" 2470c end do 2471 2472 ff = 0.0d0 2473 de = (emax-emin)/dble(npoints-1) 2474 do k=1,npoints 2475 e = emin + (k-1)*de 2476 2477 f = 0.0d0 2478 g = 0.0d0 2479 do kk=1,dosgrid(3) 2480 do jj=1,dosgrid(2) 2481 do ii=1,dosgrid(1) 2482 ishft = ii+1 2483 jshft = jj+1 2484 kshft = kk+1 2485 if (ishft.gt.dosgrid(1)) ishft=1 2486 if (jshft.gt.dosgrid(2)) jshft=1 2487 if (kshft.gt.dosgrid(3)) kshft=1 2488 do i=1,neigs 2489 ecube(1) = eigs(ii, jj, kk, i) ! (000) 2490 ecube(2) = eigs(ishft, jj, kk, i) ! (001) 2491 ecube(3) = eigs(ii, jshft, kk, i) ! (010) 2492 ecube(4) = eigs(ishft, jshft, kk, i) ! (011) 2493 ecube(5) = eigs(ii, jj, kshft, i) ! (100) 2494 ecube(6) = eigs(ishft, jj, kshft, i) ! (101) 2495 ecube(7) = eigs( ii, jshft, kshft, i) ! (110) 2496 ecube(8) = eigs(ishft, jshft, kshft, i) ! (111) 2497 2498 pcube(1) = pweight(ii, jj, kk, i) ! (000) 2499 pcube(2) = pweight(ishft, jj, kk, i) ! (001) 2500 pcube(3) = pweight(ii, jshft, kk, i) ! (010) 2501 pcube(4) = pweight(ishft, jshft, kk, i) ! (011) 2502 pcube(5) = pweight(ii, jj, kshft, i) ! (100) 2503 pcube(6) = pweight(ishft, jj, kshft, i) ! (101) 2504 pcube(7) = pweight( ii, jshft, kshft, i) ! (110) 2505 pcube(8) = pweight(ishft, jshft, kshft, i) ! (111) 2506 2507 2508 f = f + Dstates_Cube_projected(e,itetra,ecube,pcube) 2509 !g = g + Nstates_Cube_projected(e,itetra,ecube,pcube) 2510 end do 2511 end do 2512 end do 2513 end do 2514 f = f*(VT/VG) 2515 !g = g*(VT/VG) 2516 if ((k.eq.1).or.(k.eq.npoints)) then 2517 ff = ff + 0.5d0*f*de 2518 else 2519 ff = ff + f*de 2520 end if 2521 2522 !write(unit,1310) e,f*sign,g 2523 write(unit,1310) e,f*sign,ff*sign 2524 end do 2525 2526 2527 return 2528 1310 FORMAT(3E15.5) 2529 3508 FORMAT(/' Brillouin zone point: ',i5, 2530 > /' k =<',3f8.3,'> . <b1,b2,b3> ', 2531 > /' =<',3f8.3,'>') 2532 end 2533 2534 2535 real*8 function Dstates_Cube_projected(e,itetra,ecube,pcube) 2536 implicit none 2537 real*8 e 2538 integer itetra(4,6) 2539 real*8 ecube(8),pcube(8) 2540 2541* **** local variables **** 2542 integer i,j,k 2543 real*8 ds,etetra(4),ptetra(4),swap 2544 2545 real*8 Dstates_Tetra_projected 2546 external Dstates_Tetra_projected 2547 2548* **** sum over 6 tetrahedrons **** 2549 ds = 0.0d0 2550 do k=1,6 2551 etetra(1) = ecube(itetra(1,k)) 2552 etetra(2) = ecube(itetra(2,k)) 2553 etetra(3) = ecube(itetra(3,k)) 2554 etetra(4) = ecube(itetra(4,k)) 2555 2556 ptetra(1) = pcube(itetra(1,k)) 2557 ptetra(2) = pcube(itetra(2,k)) 2558 ptetra(3) = pcube(itetra(3,k)) 2559 ptetra(4) = pcube(itetra(4,k)) 2560 2561* **** bubble sort **** 2562 do j=1,3 2563 do i=j+1,4 2564 if (etetra(j).gt.etetra(i)) then 2565 swap = etetra(i) 2566 etetra(i) = etetra(j) 2567 etetra(j) = swap 2568 swap = ptetra(i) 2569 ptetra(i) = ptetra(j) 2570 ptetra(j) = swap 2571 end if 2572 end do 2573 end do 2574 ds = ds + Dstates_Tetra_projected(e,etetra,ptetra) 2575 end do 2576 2577 Dstates_cube_projected = ds 2578 return 2579 end 2580 2581 2582 real*8 function Dstates_Tetra_projected(e,ee,pp) 2583 implicit none 2584 real*8 e 2585 real*8 ee(4),pp(4) 2586 2587* **** local variables **** 2588 real*8 ds 2589 real*8 e1,e2,e4 2590 real*8 e21,e31,e41,e32,e42,e43 2591 real*8 points(3,3) 2592 2593* **** external functions **** 2594 real*8 Dstate_triangle 2595 external Dstate_triangle 2596 2597 if ((ee(1).le.e).and.(e.lt.ee(2))) then 2598c e1 = e-ee(1) 2599c e21 = ee(2) - ee(1) 2600c e31 = ee(3) - ee(1) 2601c e41 = ee(4) - ee(1) 2602c ds = 3.0d0*e1*e1/(e21*e31*e41) 2603 points(1,1) = (e-ee(1))/(ee(2)-ee(1)) 2604 points(2,1) = 0.0d0 2605 points(3,1) = 0.0d0 2606 points(1,2) = 0.0d0 2607 points(2,2) = (e-ee(1))/(ee(3)-ee(1)) 2608 points(3,2) = 0.0d0 2609 points(1,3) = 0.0d0 2610 points(2,3) = 0.0d0 2611 points(3,3) = (e-ee(1))/(ee(4)-ee(1)) 2612 ds = Dstate_triangle(points,ee,pp) 2613 2614 else if ((ee(2).le.e).and.(e.lt.ee(3))) then 2615c e2 = e-ee(2) 2616c e21 = ee(2) - ee(1) 2617c e31 = ee(3) - ee(1) 2618c e41 = ee(4) - ee(1) 2619c e32 = ee(3) - ee(2) 2620c e42 = ee(4) - ee(2) 2621c ds = (3.0d0*e21+6.0d0*e2-3.0d0*(e31+e42)*e2*e2/(e32*e42)) 2622c > /(e31*e41) 2623 points(1,1) = 0.0d0 2624 points(2,1) = (e-ee(1))/(ee(3)-ee(1)) 2625 points(3,1) = 0.0d0 2626 points(1,2) = 1.0d0 - (e-ee(2))/(ee(3)-ee(2)) 2627 points(2,2) = (e-ee(2))/(ee(3)-ee(2)) 2628 points(3,2) = 0.0d0 2629 points(1,3) = 0.0d0 2630 points(2,3) = 0.0d0 2631 points(3,3) = (e-ee(1))/(ee(4)-ee(1)) 2632 ds = Dstate_triangle(points,ee,pp) 2633 points(1,1) = 1.0d0 - (e-ee(2))/(ee(4)-ee(2)) 2634 points(2,1) = 0.0d0 2635 points(3,1) = (e-ee(2))/(ee(4)-ee(2)) 2636 points(1,2) = 1.0d0 - (e-ee(2))/(ee(3)-ee(2)) 2637 points(2,2) = (e-ee(2))/(ee(3)-ee(2)) 2638 points(3,2) = 0.0d0 2639 points(1,3) = 0.0d0 2640 points(2,3) = 0.0d0 2641 points(3,3) = (e-ee(1))/(ee(4)-ee(1)) 2642 ds = ds + Dstate_triangle(points,ee,pp) 2643 2644 else if ((ee(3).le.e).and.(e.lt.ee(4))) then 2645c e4 = ee(4)-e 2646c e41 = ee(4) - ee(1) 2647c e42 = ee(4) - ee(2) 2648c e43 = ee(4) - ee(3) 2649c ds = 3.0d0*e4*e4/(e41*e42*e43) 2650 points(1,1) = 1.0d0 - (e-ee(2))/(ee(4)-ee(2)) 2651 points(2,1) = 0.0d0 2652 points(3,1) = (e-ee(2))/(ee(4)-ee(2)) 2653 points(1,2) = 0.0d0 2654 points(2,2) = 1.0d0 - (e-ee(3))/(ee(4)-ee(3)) 2655 points(3,2) = (e-ee(3))/(ee(4)-ee(3)) 2656 points(1,3) = 0.0d0 2657 points(2,3) = 0.0d0 2658 points(3,3) = (e-ee(1))/(ee(4)-ee(1)) 2659 ds = Dstate_triangle(points,ee,pp) 2660 2661 else 2662 ds = 0.0d0 2663 end if 2664 2665 2666 Dstates_Tetra_projected = ds 2667 return 2668 end 2669 2670 2671 2672* ****************************************** 2673* * * 2674* * Dstate_triangle * 2675* * * 2676* ****************************************** 2677 real*8 function Dstate_triangle(points,ee,ff) 2678 implicit none 2679 real*8 points(3,3) 2680 real*8 ee(4),ff(4) 2681 2682 real*8 p13x,p13y,p13z,p23x,p23y,p23z 2683 real*8 f0,f1,f2,e10,e20,e30,nde,n2 2684 real*8 n(3) 2685 2686 p13x = points(1,1) - points(1,3) 2687 p13y = points(2,1) - points(2,3) 2688 p13z = points(3,1) - points(3,3) 2689 2690 p23x = points(1,2) - points(1,3) 2691 p23y = points(2,2) - points(2,3) 2692 p23z = points(3,2) - points(3,3) 2693 2694 n(1) = p13y*p23z - p13z*p23y 2695 n(2) = p13z*p23x - p13x*p23z 2696 n(3) = p13x*p23y - p13y*p23x 2697 n2 = n(1)*n(1) + n(2)*n(2) + n(3)*n(3) 2698 2699 e10 = ee(2)-ee(1) 2700 e20 = ee(3)-ee(1) 2701 e30 = ee(4)-ee(1) 2702 nde = dabs(e10*n(1) + e20*n(2) + e30*n(3)) 2703 2704 f0 = ff(1) + (ff(2)-ff(1))*points(1,3) 2705 > + (ff(3)-ff(1))*points(2,3) 2706 > + (ff(4)-ff(1))*points(3,3) 2707 2708 f1 = (ff(2)-ff(1))*(points(1,1)-points(1,3)) 2709 > + (ff(3)-ff(1))*(points(2,1)-points(2,3)) 2710 > + (ff(4)-ff(1))*(points(3,1)-points(3,3)) 2711 2712 f2 = (ff(2)-ff(1))*(points(1,2)-points(1,3)) 2713 > + (ff(3)-ff(1))*(points(2,2)-points(2,3)) 2714 > + (ff(4)-ff(1))*(points(3,2)-points(3,3)) 2715 2716 Dstate_triangle = 6.0d0*(n2/nde)*(f0/2.0d0 + f1/6.0d0 + f2/6.0d0) 2717 return 2718 end 2719 2720 2721 2722 2723 2724 real*8 function Nstates_Cube_projected(e,itetra,ecube,pcube) 2725 implicit none 2726 real*8 e 2727 integer itetra(4,6) 2728 real*8 ecube(8),pcube(8) 2729 2730* **** local variables **** 2731 integer i,j,k 2732 real*8 ds,etetra(4),ptetra(4),swap 2733 2734 real*8 Nstates_Tetra_projected 2735 external Nstates_Tetra_projected 2736 2737* **** sum over 6 tetrahedrons **** 2738 ds = 0.0d0 2739 do k=1,6 2740 etetra(1) = ecube(itetra(1,k)) 2741 etetra(2) = ecube(itetra(2,k)) 2742 etetra(3) = ecube(itetra(3,k)) 2743 etetra(4) = ecube(itetra(4,k)) 2744 ptetra(1) = pcube(itetra(1,k)) 2745 ptetra(2) = pcube(itetra(2,k)) 2746 ptetra(3) = pcube(itetra(3,k)) 2747 ptetra(4) = pcube(itetra(4,k)) 2748 2749* **** bubble sort **** 2750 do j=1,3 2751 do i=j+1,4 2752 if (etetra(j).gt.etetra(i)) then 2753 swap = etetra(i) 2754 etetra(i) = etetra(j) 2755 etetra(j) = swap 2756 swap = ptetra(i) 2757 ptetra(i) = ptetra(j) 2758 ptetra(j) = swap 2759 end if 2760 end do 2761 end do 2762 ds = ds + Nstates_Tetra_projected(e,etetra,ptetra) 2763 end do 2764 2765 Nstates_cube_projected = ds 2766 return 2767 end 2768 2769 2770 real*8 function Nstates_Tetra_projected(e,ee,pp) 2771 implicit none 2772 real*8 e 2773 real*8 ee(4),pp(4) 2774 2775* **** local variables **** 2776 real*8 ds 2777 real*8 e1,e2,e4 2778 real*8 e21,e31,e41,e32,e42,e43 2779 2780 if ((ee(1).le.e).and.(e.lt.ee(2))) then 2781 e1 = e-ee(1) 2782 e21 = ee(2) - ee(1) 2783 e31 = ee(3) - ee(1) 2784 e41 = ee(4) - ee(1) 2785 ds = e1*e1*e1/(e21*e31*e41) 2786 else if ((ee(2).le.e).and.(e.lt.ee(3))) then 2787 e2 = e-ee(2) 2788 e21 = ee(2) - ee(1) 2789 e31 = ee(3) - ee(1) 2790 e41 = ee(4) - ee(1) 2791 e32 = ee(3) - ee(2) 2792 e42 = ee(4) - ee(2) 2793 ds = (e21*e21 2794 > + 3.0d0*e21*e2 2795 > + 3.0d0*e2*e2 2796 > - (e31+e42)*e2*e2*e2/(e32*e42)) 2797 > /(e31*e41) 2798 else if ((ee(3).le.e).and.(e.lt.ee(4))) then 2799 e4 = ee(4)-e 2800 e41 = ee(4) - ee(1) 2801 e42 = ee(4) - ee(2) 2802 e43 = ee(4) - ee(3) 2803 ds = 1.0d0 - e4*e4*e4/(e41*e42*e43) 2804 else if (e.ge.ee(4)) then 2805 ds = 1.0d0 2806 else 2807 ds = 0.0d0 2808 end if 2809 2810 2811 Nstates_Tetra_projected = ds 2812 return 2813 end 2814 2815 2816