1c 2c $Id$ 3c 4 5* **************************************** 6* * * 7* * pspw_analysis * 8* * * 9* **************************************** 10 11 subroutine pspw_analysis(flag,rtdb,ispin,ne,psi,eig) 12 implicit none 13 integer flag 14 integer rtdb 15 integer ispin,ne(2) 16 complex*16 psi(*) 17 real*8 eig(*) 18 19#include "bafdecls.fh" 20#include "tcgmsg.fh" 21#include "msgtypesf.h" 22#include "btdb.fh" 23#include "errquit.fh" 24#include "stdio.fh" 25 26 27* ================================================ 28* This code is a rewrite of an earlier analysis 29* code written by R. Kawai. 30* 31* VERSION MPI-1.00 11/15/96 by Eric Bylaska 32* 33* ================================================ 34 35* **** parallel variables **** 36 integer taskid 37 integer MASTER 38 parameter(MASTER=0) 39 40 41* **** electronic variables **** 42 integer npack1 43 integer n1(2),n2(2) 44 45 46 47* **** local variables **** 48 integer i,k,l,n,ms,l1,l2,j 49 integer ll,spin,ind 50 real*8 ttl1,ttl2,subttl 51 52 integer h_actlist,l_actlist,nactive_atoms,ma_type 53 integer npoints,ii,lmax,norbs,naos 54 integer weight(2),coef(2) 55 real*8 emin,emax,alpha,lmbda,rcut,w 56 character*255 filename 57 58 logical value,mulliken_kawai,fixatoms 59 character*28 DD 60 character*255 id,test 61 integer npsp,nion,nemax 62 integer lorb(2) ! integer lorb(npsp) 63 integer b0(2) ! real*8 b0(0:5,npsp) 64 integer a(2) ! real*8 a(36,nemax,nion) 65 integer total(2) ! real*8 total(nion) 66 integer sum(2) ! real*8 sum(nemax) 67 integer subtl(2) ! real*8 subtl(0:5,3) 68 69 character*4 spn(2) 70 DATA SPN / 'UP ', 'DOWN' / 71 72 73* **** external functions **** 74 logical control_DOS,nwpw_filefind 75 logical control_mulliken_kawai,aorbs_init,aorbs_readall 76 external control_DOS,nwpw_filefind 77 external control_mulliken_kawai,aorbs_init,aorbs_readall 78 character spdf_name 79 external spdf_name 80 character*4 ion_atom_qm 81 external ion_atom_qm 82 integer ion_nion_qm,ion_katm_qm,ion_nkatm_qm 83 external ion_nion_qm,ion_katm_qm,ion_nkatm_qm 84 real*8 ion_rion 85 external ion_rion 86 real*8 lattice_omega,lattice_ecut,lattice_unita 87 external lattice_omega,lattice_ecut,lattice_unita 88 real*8 control_mulliken_rcut,control_mulliken_lmbda 89 external control_mulliken_rcut,control_mulliken_lmbda 90 real*8 aorbs_rcut,aorbs_lmbda 91 external aorbs_rcut,aorbs_lmbda 92 93 character*7 c_index_name 94 external c_index_name 95 96 97 call Parallel_taskid(taskid) 98 call Pack_npack(1,npack1) 99 mulliken_kawai = control_mulliken_kawai() 100 101 npsp = ion_nkatm_qm() 102 nion = ion_nion_qm() 103 nemax = ne(1)+ne(2) 104 105 n1(1) = 1 106 n2(1) = ne(1) 107 n1(2) = ne(1)+1 108 n2(2) = ne(1)+ne(2) 109 110 value = BA_alloc_get(mt_int,npsp,'lorb',lorb(2),lorb(1)) 111 value = value.and. 112 > BA_alloc_get(mt_dbl,6*npsp,'b0',b0(2),b0(1)) 113 value = value.and. 114 > BA_alloc_get(mt_dbl,36*nemax*nion,'a',a(2),a(1)) 115 value = value.and. 116 > BA_alloc_get(mt_dbl,nion,'total',total(2),total(1)) 117 value = value.and. 118 > BA_alloc_get(mt_dbl,nemax,'sum',sum(2),sum(1)) 119 value = value.and. 120 > BA_alloc_get(mt_dbl,6*3,'subtl',subtl(2),subtl(1)) 121 122 call dcopy(36*nemax*nion,0.0d0,0,dbl_mb(a(1)),1) 123 124 125* *********************************************** 126* **** create psp1 files if they don't exist **** 127* *********************************************** 128 do k=1,npsp 129 DD = ' ' 130 DD = ion_atom_qm(k) 131 ind = index(DD,' ') - 1 132 test = DD(1:ind)//'.psp1' 133 id = DD(1:ind)//'.aorb' 134 call control_mullikenparameters(ion_atom_qm(k),rcut,lmbda) 135 if ((.not.nwpw_filefind(test)).or.(.not.nwpw_filefind(id))) 136 > call aorbs_formatter_auto(ion_atom_qm(k),rcut,lmbda) 137 end do 138 139* **************************************** 140* **** read in expansion coefficients **** 141* **************************************** 142 do k=1,npsp 143 id = 'analysis:lorb'//ion_atom_qm(k) 144 if (.not. btdb_get(rtdb,id,mt_int,1,int_mb(lorb(1)+k-1))) then 145 DD = ' ' 146 DD = ion_atom_qm(k) 147 ind = index(DD,' ') - 1 148 test = DD(1:ind)//'.psp1' 149 150 !write(*,*) "test:",test,ind 151 value = btdb_parallel(.false.) 152 if (taskid.eq.MASTER) then 153 ind = index(test,' ') - 1 154 !write(*,*) "test:",test,ind 155 call analysis_expansion_coef(test,-1,rtdb) 156 end if 157 value = btdb_parallel(.true.) 158 call ga_sync() 159 160 if (.not. btdb_get(rtdb,id,mt_int,1,int_mb(lorb(1)+k-1))) 161 > call errquit( 162 > 'analysis: btdb_get lorb failed', 0, RTDB_ERR) 163 end if 164 165 id = 'analysis:expansion'//ion_atom_qm(k) 166 if (.not. btdb_get(rtdb,id,mt_dbl,(int_mb(lorb(1)+k-1)+1), 167 > dbl_mb(b0(1)+(k-1)*6))) then 168 DD = ' ' 169 DD = ion_atom_qm(k) 170 ind = index(DD,' ') -1 171 test = DD(1:ind)//'.psp1' 172 call analysis_expansion_coef(test,-1,rtdb) 173 174 if (.not. btdb_get(rtdb,id,mt_dbl,(int_mb(lorb(1)+k-1)+1), 175 > dbl_mb(b0(1)+(k-1)*6))) 176 > call errquit( 177 > 'analysis: btdb_get failed', 0, RTDB_ERR) 178 end if 179 end do 180 181 if (.not.mulliken_kawai) then 182 value = aorbs_init() 183 value = value.and.aorbs_readall() 184 if (.not.value) go to 1901 185 186 value = .true. 187 do k=1,npsp 188 call control_mullikenparameters(ion_atom_qm(k),rcut,lmbda) 189 if ((dabs(aorbs_rcut(k)-rcut).gt.1.0d-6).or. 190 > (dabs(aorbs_lmbda(k)-lmbda).gt.1.0d-6)) then 191 call aorbs_formatter_auto(ion_atom_qm(k),rcut,lmbda) 192 value = .false. 193 end if 194 end do 195 if (.not.value) value = aorbs_readall() 196 if (.not.value) go to 1901 197 end if 198 199 200 201 if (taskid.eq.MASTER) then 202 call util_date(DD) 203 204 WRITE(luout,*) 205 WRITE(luout,*) 206 WRITE(luout,*) 207 WRITE(luout,*) 208 > '*************************************************************' 209 WRITE(luout,*) 210 > '** **' 211 WRITE(luout,*) 212 > '** PSPW Mulliken analysis **' 213 WRITE(luout,*) 214 > '** **' 215 if (flag.eq.1) 216 > WRITE(luout,*) 217 > '** (Virtual Orbitals) **' 218 WRITE(luout,*) 219 > '** Population analysis algorithm devloped by Ryoichi Kawai **' 220 WRITE(luout,*) 221 > '** **' 222 WRITE(luout,1000) DD 223 1000 FORMAT( 224 > ' ** ',A16,' **') 225 WRITE(luout,*) 226 > '** **' 227 WRITE(luout,*) 228 > '*************************************************************' 229 end if 230 231 232c **** ouput xyz format **** 233 call ion_Print_XYZ(luout) 234 235 236 value = btdb_parallel(.false.) 237 if (taskid.eq.MASTER) then 238 if (mulliken_kawai) then 239 WRITE(luout,1300) 240 WRITE(luout,1305) 'ATOM ','S','P','D','F' 241 DO k=1,npsp 242 WRITE(luout,1306) ion_atom_qm(K), 243 > (dbl_mb(b0(1)+l+(k-1)*6), 244 > l=0,int_mb(lorb(1)+k-1)) 245 END DO 246 call util_flush(luout) 247 1300 FORMAT(//'== Kawai Expansion Coefficients =='/) 248 1305 FORMAT(A5,6X,A,12X,A,12X,A,12X,A) 249 1306 FORMAT(A4,' : ',4E13.5) 250 else 251 write(luout,1307) 252 do k=1,npsp 253 call control_mullikenparameters(ion_atom_qm(k),rcut,lmbda) 254 if (lmbda.gt.0.0d0) then 255 write(luout,1308) ion_atom_qm(k),"damping",rcut,lmbda 256 else 257 write(luout,*) ion_atom_qm(k)," nodamping" 258 end if 259 end do 260 call util_flush(luout) 261 1307 FORMAT(//'== Atomic Orbital Expansion =='/) 262 1308 FORMAT(A5,A10,4x,"rcut=",F8.3,2x,"lmbda=",F8.3) 263 end if 264 end if 265 value = btdb_parallel(.true.) 266 267 268 269 call util_file_name('ORBOUT', 270 > .true., 271 > .false., 272 > id) 273 if (taskid.eq.MASTER) 274 > OPEN(UNIT=65,FILE=id,FORM='FORMATTED') 275 276 call Orb_Analysis(65,flag,ispin,ne,npack1,nemax,psi, 277 > int_mb(lorb(1)), 278 > dbl_mb(b0(1)), 279 > dbl_mb(a(1)), 280 > dbl_mb(sum(1))) 281 282 if (taskid.eq.MASTER) close(unit=65) 283 284 285 286 if (taskid.eq.MASTER) then 287 WRITE(luout,*) 288 WRITE(luout,*) 289 WRITE(luout,*) 290 > '=====================================================' 291 if (flag.eq.0) 292 >WRITE(luout,*) 293 > '| POPULATION ANALYSIS OF FILLED MOLECULAR ORBITALS |' 294 if (flag.eq.1) 295 >WRITE(luout,*) 296 > '| POPULATION ANALYSIS OF VIRTUAL MOLECULAR ORBITALS |' 297 WRITE(luout,*) 298 > '=====================================================' 299 if (.not.mulliken_kawai) then 300 WRITE(luout,1311) 301 else 302 WRITE(luout,1312) 303 end if 304c WRITE(6,1313) 305c WRITE(6,1314) 306c WRITE(6,1315) 307 1311 FORMAT(//'== Using pseudoatomic orbital expansion ==') 308 1312 FORMAT(//'== Using Kawai projected atomic expansion ==') 309c1313 FORMAT('== order of orbitals: s ==') 310c1314 FORMAT('== px, py, pz ==') 311c1315 FORMAT('== dzz, dx2-y2, dxy, dyz, dzx ==') 312 313 DO SPIN=1,ISPIN 314 DO N=N1(SPIN),N2(SPIN) 315 WRITE(6,1500) 316 IF(ISPIN.EQ.2) THEN 317 WRITE(6,1510) N,SPN(SPIN),dbl_mb(sum(1)+N-1),eig(n), 318 > eig(n)*27.2116d0 319 ELSE 320 WRITE(6,1515) N,dbl_mb(sum(1)+N-1),eig(n), 321 > eig(n)*27.2116d0 322 ENDIF 323 !write(6,1519) 324 WRITE(6,1520) 'NO','ATOM','L','POPULATION' 325 326 DO L=0,5 327 dbl_mb(SUBTL(1)+L)=0.0d0 328 END DO 329 DO I=1,nion 330 dbl_mb(TOTAL(1)+I-1)=0.0d0 331 DO L=0,int_mb(lorb(1)+ion_katm_qm(I)-1) 332 L1=L**2+1 333 L2=(L+1)**2 334 SUBTTL=0.0d0 335 DO LL=L1,L2 336 dbl_mb(TOTAL(1)+I-1)=dbl_mb(TOTAL(1)+I-1)+ 337 > dbl_mb(A(1)+LL-1+(N-1)*36+(I-1)*36*nemax)**2 338 SUBTTL=SUBTTL+ 339 > dbl_mb(A(1)+LL-1+(N-1)*36+(I-1)*36*nemax)**2 340 END DO 341 dbl_mb(SUBTL(1)+L)=dbl_mb(SUBTL(1)+L)+SUBTTL 342 if (l.eq.0) write(6,1516) 343 if (l.eq.1) write(6,1517) 344 if (l.eq.2) write(6,1518) 345 if (l.eq.3) write(6,1519) 346 WRITE(luout,1530) I,ion_atom_qm(ion_katm_qm(I)),L,SUBTTL, 347 > (dbl_mb(A(1)+LL-1+(N-1)*36+(I-1)*36*nemax), 348 > LL=L1,L2) 349 350 END DO 351 END DO 352 353 WRITE(luout,1540) 354 WRITE(luout,1550) (I,ion_atom_qm(ion_katm_qm(I)), 355 > dbl_mb(TOTAL(1)+I-1),I=1,nion) 356 WRITE(luout,1555) 357 WRITE(luout,1560) 's','p','d','f' 358 WRITE(luout,1570) (dbl_mb(SUBTL(1)+L),L=0,3) 359 END DO 360 END DO 361 call util_flush(luout) 362 363 1500 FORMAT(//'------------------------------------------------', 364 > '------------------------------'//) 365 1510 FORMAT('*** ORBITAL=',I4,'*** SPIN=',A4,4X,'SUM=',E12.5, 366 > ' E=',E12.5,' (',F8.3,'eV)'/) 367 1515 FORMAT('*** ORBITAL=',I4,'*** SPIN=BOTH',4X,'SUM=',E12.5, 368 > ' E=',E12.5,' (',F8.3,'eV)'/) 369 1516 FORMAT(27x,' s') 370 1517 FORMAT(27x,' px pz py') 371 1518 FORMAT(27x,' dx2-y2 dzx d3z2-1 dyz dxy') 372 1519 FORMAT(27x, 373 > ' fx(x2-3y2) fz(5z2-1) fx(5z2-1) fz(5z2-3) fy(5z2-1) ', 374 > ' fxyz fy(3x2-y2)') 375 376c 1519 FORMAT(30x,' s', 377c > /30x,' px py pz ' 378c > /27x,'d3z2-1 dxy dyz dzx dx2-y2', 379c > /30x,' fy(3x2-y2) fxyz fy(5z2-1) fz(5z2-3) fx(5z2-1)', 380c > ' fz(5z2-1) fx(x2-3y2)') 381 382 1520 FORMAT(A2,2X,A4,2X,A1,2X,A10) 383 1530 FORMAT(I2,3X,A4,3X,I1,8F11.5) 384 1531 FORMAT(8F11.5) 385 1540 FORMAT(//'=== DISTRIBUTION ==='/) 386 1550 FORMAT(4(I6,'(',A4,')',F9.4)) 387 1555 FORMAT(//'== ANGULAR MOMENTUM POPULATIONS ==='/) 388 1560 FORMAT(6X,A1,3(9X,A1)) 389 1570 FORMAT(4F10.4) 390 391 WRITE(6,*) 392 WRITE(6,*) 393 WRITE(6,*) '========================================' 394 WRITE(6,*) '| POPULATION ANALYSIS ON EACH ATOM |' 395 WRITE(6,*) '========================================' 396 WRITE(6,*) 397 WRITE(6,*) 398 WRITE(6,1610) 'NO','ATOM','SPIN','TOTAL','s','p','d','f' 399 DO I=1,nion 400 TTL1=0.0d0 401 TTL2=0.0d0 402 403 DO L=0,int_mb(lorb(1)+ion_katm_qm(I)-1) 404 DO ms=1,ISPIN 405 dbl_mb(SUBTL(1)+L+(ms-1)*6)=0.0d0 406 DO N=N1(ms),N2(ms) 407 L1=L**2+1 408 L2=(L+1)**2 409 DO LL=L1,L2 410 dbl_mb(SUBTL(1)+L+(ms-1)*6)= 411 > dbl_mb(SUBTL(1)+L+(ms-1)*6) + 412 > dbl_mb(A(1)+LL-1+(N-1)*36+(I-1)*36*nemax)**2 413 414 END DO 415 END DO 416 END DO 417 END DO 418 419 TTL1=0.0d0 420 DO L=0,int_mb(lorb(1)+ion_katm_qm(I)-1) 421 TTL1=TTL1+dbl_mb(SUBTL(1)+L) 422 END DO 423 WRITE(6,1620) I,ion_atom_qm(ion_katm_qm(I)),SPN(1),TTL1, 424 > ( dbl_mb(SUBTL(1)+L), 425 > L=0,int_mb(lorb(1)+ion_katm_qm(I)-1) ) 426 TTL1=0.0d0 427 DO L=0,int_mb(lorb(1)+ion_katm_qm(I)-1) 428 TTL1=TTL1+dbl_mb(SUBTL(1)+L+(ispin-1)*6) 429 END DO 430 WRITE(6,1620) I,ion_atom_qm(ion_katm_qm(I)),SPN(2),TTL1, 431 > ( dbl_mb(SUBTL(1)+L+(ispin-1)*6), 432 > L=0,int_mb(lorb(1)+ion_katm_qm(I)-1) ) 433 434 END DO 435 call util_flush(6) 436 437 1610 FORMAT(A2,2X,A4,2X,A4,4X,A5,7X,A,10X,A,10X,A,10X,A) 438 1620 FORMAT(I2,3X,A4,3X,A4,5F11.5) 439 440 DO L=0,3 441 dbl_mb(SUBTL(1)+L) =0.0d0 442 dbl_mb(SUBTL(1)+L+6)=0.0d0 443 END DO 444 DO I=1,nion 445 DO SPIN=1,ISPIN 446 DO N=N1(SPIN),N2(SPIN) 447 DO L=0,int_mb(lorb(1)+ion_katm_qm(I)-1) 448 L1=L**2+1 449 L2=(L+1)**2 450 DO LL=L1,L2 451 dbl_mb(SUBTL(1)+L+(SPIN-1)*6)= 452 > dbl_mb(SUBTL(1)+L+(SPIN-1)*6)+ 453 > dbl_mb(A(1)+LL-1+(N-1)*36+(I-1)*36*nemax)**2 454 END DO 455 END DO 456 END DO 457 END DO 458 END DO 459 460 DO L=0,3 461 dbl_mb(SUBTL(1)+L+2*6)= 462 > (dbl_mb(SUBTL(1)+L)+dbl_mb(SUBTL(1)+L+(ISPIN-1)*6)) 463 > *100.d0/(NE(1)+NE(ISPIN)) 464 dbl_mb(SUBTL(1)+L)=dbl_mb(SUBTL(1)+L)*100.0d0/dble(NE(1)) 465 IF((ISPIN.EQ.2).and.(NE(2).gt.0)) 466 > dbl_mb(SUBTL(1)+L+6)=dbl_mb(SUBTL(1)+L+6)*100.0d0/dble(NE(2)) 467 END DO 468 469 WRITE(6,1700) 470 WRITE(6,1710) ' SPIN ','s','p','d','f' 471 WRITE(6,1720) SPN(1),(dbl_mb(SUBTL(1)+L),L=0,3) 472 WRITE(6,1720) SPN(ISPIN),(dbl_mb(SUBTL(1)+L+(ISPIN-1)*6),L=0,3) 473 WRITE(6,1720) ' TOTAL',(dbl_mb(SUBTL(1)+L+(3-1)*6),L=0,3) 474 call util_flush(6) 475 1700 FORMAT(///'=== TOTAL ANGULAR MOMENTUM POPULATION ==='/) 476 1710 FORMAT(A6,6X,A1,3(11X,A1)) 477 1720 FORMAT(A6,4(F10.2,'% ')) 478 479 end if 480 481 482* *********************************************** 483* **** generate projected DENSITY OF STATES ***** 484* *********************************************** 485 if (control_DOS()) then 486 if (taskid.eq.MASTER) write(*,*) "into mulliken DOS" 487 488 value=BA_push_get(mt_dbl,(nemax),'weight',weight(2),weight(1)) 489 if (.not. value) 490 > call errquit('analysis:out of stack memory',0, MA_ERR) 491 call dcopy(nemax,1.0d0,0,dbl_mb(weight(1)),1) 492 493 494 if (.not.btdb_get(rtdb,'dos:alpha',mt_dbl,1,alpha)) then 495 alpha = 0.05d0/27.2116d0 496 end if 497 498 if (.not.btdb_get(rtdb,'dos:npoints',mt_int,1,npoints)) then 499 npoints = 500 500 end if 501 502 if (.not.btdb_get(rtdb,'dos:emin',mt_dbl,1,emin)) then 503 emin = 99999.0d0 504 do ii=1,ne(1)+ne(2) 505 if (eig(ii).lt.emin) emin = eig(ii) 506 end do 507 emin = emin - 0.1d0 508 end if 509 510 if (.not.btdb_get(rtdb,'dos:emax',mt_dbl,1,emax)) then 511 emax = -99999.0d0 512 do ii=1,ne(1)+ne(2) 513 if (eig(ii).gt.emax) emax = eig(ii) 514 end do 515 emax = emax + 0.1d0 516 end if 517 518* **** explicit number of atoms have been requested **** 519 fixatoms = .false. 520 if (rtdb_ma_get(rtdb, 'nwpw:dos:actlist', ma_type, 521 > nactive_atoms, h_actlist)) then 522 523 if (.not.BA_get_index(h_actlist,l_actlist)) 524 > call errquit( 525 > 'analysis: ma_get_index failed for actlist',911, 526 & MA_ERR) 527 528 fixatoms = .true. 529 end if 530 531 if (taskid.eq.MASTER) then 532 write(6,1800) 533 if (.not.fixatoms) write(6,1801) 534 if (fixatoms) then 535 write(6,1802) 536 write(6,1803) (int_mb(l_actlist+j-1),j=1,nactive_atoms) 537 end if 538 end if 539 1800 FORMAT(///'=== PROJECTED DENSITY OF STATES ==='/) 540 1801 FORMAT(' All atoms were used to determine weights') 541 1802 FORMAT(' The following atoms were used to determine weights:') 542 1803 FORMAT(2x,8I6) 543 544* **** angular momentum decomposition ***** 545 lmax = -1 546 do k=1,npsp 547 if (lmax.le.int_mb(lorb(1)+k-1)) lmax = int_mb(lorb(1)+k-1) 548 end do 549 550 do L=0,lmax 551 552 call dcopy(nemax,0.0d0,0,dbl_mb(weight(1)),1) 553 554 if (.not.fixatoms) then 555 DO I=1,nion 556 DO SPIN=1,ISPIN 557 DO N=N1(SPIN),N2(SPIN) 558 if (L.le.int_mb(lorb(1)+ion_katm_qm(I)-1)) then 559 L1=L**2+1 560 L2=(L+1)**2 561 DO LL=L1,L2 562 dbl_mb(weight(1)+n-1)= 563 > dbl_mb(weight(1)+n-1)+ 564 > dbl_mb(A(1)+LL-1+(N-1)*36+(I-1)*36*nemax)**2 565 END DO 566 end if 567 END DO 568 END DO 569 END DO 570 else 571 DO j=1,nactive_atoms 572 I=int_mb(l_actlist+j-1) 573 DO SPIN=1,ISPIN 574 DO N=N1(SPIN),N2(SPIN) 575 if (L.le.int_mb(lorb(1)+ion_katm_qm(I)-1)) then 576 L1=L**2+1 577 L2=(L+1)**2 578 DO LL=L1,L2 579 dbl_mb(weight(1)+n-1)= 580 > dbl_mb(weight(1)+n-1)+ 581 > dbl_mb(A(1)+LL-1+(N-1)*36+(I-1)*36*nemax)**2 582 END DO 583 end if 584 END DO 585 END DO 586 END DO 587 end if 588 589 if (ispin.eq.1) then 590 if (flag.eq.0) filename = "mulliken_fdos_both_"//spdf_name(l) 591 if (flag.eq.1) filename = "mulliken_vdos_both_"//spdf_name(l) 592 call densityofstates(filename,.false., 593 > eig,dbl_mb(weight(1)),ne(1), 594 > 1.0d0,alpha,npoints,emin,emax) 595 if (flag.eq.0) value = .true. 596 if (flag.eq.1) value = .false. 597 filename = "mulliken_dos_both_"//spdf_name(l) 598 call densityofstates(filename,value, 599 > eig,dbl_mb(weight(1)),ne(1), 600 > 1.0d0,alpha,npoints,emin,emax) 601 end if 602 if (ispin.eq.2) then 603 if (flag.eq.0) filename= "mulliken_fdos_alpha_"//spdf_name(l) 604 if (flag.eq.1) filename= "mulliken_vdos_alpha_"//spdf_name(l) 605 call densityofstates(filename,.false., 606 > eig,dbl_mb(weight(1)),ne(1), 607 > 1.0d0,alpha,npoints,emin,emax) 608 if (flag.eq.0) value = .true. 609 if (flag.eq.1) value = .false. 610 filename= "mulliken_dos_alpha_"//spdf_name(l) 611 call densityofstates(filename,value, 612 > eig,dbl_mb(weight(1)),ne(1), 613 > 1.0d0,alpha,npoints,emin,emax) 614 615 if (flag.eq.0) filename = "mulliken_fdos_beta_"//spdf_name(l) 616 if (flag.eq.1) filename = "mulliken_vdos_beta_"//spdf_name(l) 617 call densityofstates(filename,.false., 618 > eig(1+ne(1)),dbl_mb(weight(1)+ne(1)),ne(2), 619 > -1.0d0,alpha,npoints,emin,emax) 620 if (flag.eq.0) value = .true. 621 if (flag.eq.1) value = .false. 622 filename= "mulliken_dos_beta_"//spdf_name(l) 623 call densityofstates(filename,value, 624 > eig(1+ne(1)),dbl_mb(weight(1)+ne(1)),ne(2), 625 > -1.0d0,alpha,npoints,emin,emax) 626 end if 627 end do 628 629* **** combined angular momentum decomposition ***** 630 call dcopy(nemax,0.0d0,0,dbl_mb(weight(1)),1) 631 if (.not.fixatoms) then 632 DO I=1,nion 633 DO SPIN=1,ISPIN 634 DO N=N1(SPIN),N2(SPIN) 635 DO L=0,int_mb(lorb(1)+ion_katm_qm(I)-1) 636 L1=L**2+1 637 L2=(L+1)**2 638 DO LL=L1,L2 639 dbl_mb(weight(1)+n-1)= 640 > dbl_mb(weight(1)+n-1)+ 641 > dbl_mb(A(1)+LL-1+(N-1)*36+(I-1)*36*nemax)**2 642 END DO 643 END DO 644 END DO 645 END DO 646 END DO 647 else 648 DO j=1,nactive_atoms 649 I=int_mb(l_actlist+j-1) 650 DO SPIN=1,ISPIN 651 DO N=N1(SPIN),N2(SPIN) 652 DO L=0,int_mb(lorb(1)+ion_katm_qm(I)-1) 653 L1=L**2+1 654 L2=(L+1)**2 655 DO LL=L1,L2 656 dbl_mb(weight(1)+n-1)= 657 > dbl_mb(weight(1)+n-1)+ 658 > dbl_mb(A(1)+LL-1+(N-1)*36+(I-1)*36*nemax)**2 659 END DO 660 END DO 661 END DO 662 END DO 663 END DO 664 end if 665 666 if (ispin.eq.1) then 667 if (flag.eq.0) filename = "mulliken_fdos_both_all" 668 if (flag.eq.1) filename = "mulliken_vdos_both_all" 669 call densityofstates(filename,.false., 670 > eig,dbl_mb(weight(1)),ne(1), 671 > 1.0d0,alpha,npoints,emin,emax) 672 673 if (flag.eq.0) value = .true. 674 if (flag.eq.1) value = .false. 675 filename = "mulliken_dos_both_all" 676 call densityofstates(filename,value, 677 > eig,dbl_mb(weight(1)),ne(1), 678 > 1.0d0,alpha,npoints,emin,emax) 679 end if 680 if (ispin.eq.2) then 681 if (flag.eq.0) filename = "mulliken_fdos_alpha_all" 682 if (flag.eq.1) filename = "mulliken_vdos_alpha_all" 683 call densityofstates(filename,.false., 684 > eig,dbl_mb(weight(1)),ne(1), 685 > 1.0d0,alpha,npoints,emin,emax) 686 687 if (flag.eq.0) value = .true. 688 if (flag.eq.1) value = .false. 689 filename = "mulliken_dos_alpha_all" 690 call densityofstates(filename,value, 691 > eig,dbl_mb(weight(1)),ne(1), 692 > 1.0d0,alpha,npoints,emin,emax) 693 694 if (flag.eq.0) filename = "mulliken_fdos_beta_all" 695 if (flag.eq.1) filename = "mulliken_vdos_beta_all" 696 call densityofstates(filename,.false., 697 > eig(1+ne(1)),dbl_mb(weight(1)+ne(1)),ne(2), 698 > -1.0d0,alpha,npoints,emin,emax) 699 700 if (flag.eq.0) value = .true. 701 if (flag.eq.1) value = .false. 702 filename = "mulliken_dos_beta_all" 703 call densityofstates(filename,.false., 704 > eig(1+ne(1)),dbl_mb(weight(1)+ne(1)),ne(2), 705 > -1.0d0,alpha,npoints,emin,emax) 706 end if 707 708* **** ORBITAL DOS ***** 709 value=BA_push_get(mt_dbl,36,'coef',coef(2),coef(1)) 710 if (.not.value) 711 > call errquit('analysis:out of stack memory',1, MA_ERR) 712 713 id = 'nwpw:dos:orb:norb' 714 if (.not.btdb_get(rtdb,id,mt_int,1,norbs)) norbs = 0 715 716 do k=1,norbs 717 id = 'nwpw:dos:orb:coef'//c_index_name(k) 718 if (.not.rtdb_get_info(rtdb,id,ma_type,naos,test)) 719 > call errquit( 720 > 'analysis: ma_get_index failed for nwpw:dos:orb:coef',911, 721 > MA_ERR) 722 if (.not. btdb_get(rtdb,id,mt_dbl,naos,dbl_mb(coef(1)))) 723 > call errquit( 724 > 'analysis: ma_get_index failed for nwpw:dos:orb:coef',912, 725 > MA_ERR) 726 !*** normalize coef *** 727 w = 0.0 728 do i=1,naos 729 w = w + dbl_mb(coef(1)+i-1)**2 730 end do 731 w = 1.0d0/dsqrt(w) 732 do i=1,naos 733 dbl_mb(coef(1)+i-1) = dbl_mb(coef(1)+i-1)*w 734 end do 735 736 call dcopy(nemax,0.0d0,0,dbl_mb(weight(1)),1) 737 if (.not.fixatoms) then 738 do I=1,nion 739 do SPIN=1,ISPIN 740 do N=N1(SPIN),N2(SPIN) 741 call dgemm('N','N',1,1,naos, 742 > 1.0d0, 743 > dbl_mb(coef(1)),1, 744 > dbl_mb(A(1)+(N-1)*36+(I-1)*36*nemax),naos, 745 > 0.0d0, 746 > w,1) 747 dbl_mb(weight(1)+n-1) = dbl_mb(weight(1)+n-1) + w**2 748 end do 749 end do 750 end do 751 else 752 do j=1,nactive_atoms 753 I=int_mb(l_actlist+j-1) 754 do SPIN=1,ISPIN 755 do N=N1(SPIN),N2(SPIN) 756 call dgemm('N','N',1,1,naos, 757 > 1.0d0, 758 > dbl_mb(coef(1)),1, 759 > dbl_mb(A(1)+(N-1)*36+(I-1)*36*nemax),naos, 760 > 0.0d0, 761 > w,1) 762 dbl_mb(weight(1)+n-1) = dbl_mb(weight(1)+n-1) + w**2 763 end do 764 end do 765 end do 766 end if 767 if (ispin.eq.1) then 768 if (flag.eq.0) 769 > filename = "mulliken_fdos_both_orb"//c_index_name(k) 770 if (flag.eq.1) 771 > filename = "mulliken_vdos_both_orb"//c_index_name(k) 772 call densityofstates(filename,.false., 773 > eig,dbl_mb(weight(1)),ne(1), 774 > 1.0d0,alpha,npoints,emin,emax) 775 if (flag.eq.0) value = .true. 776 if (flag.eq.1) value = .false. 777 filename = "mulliken_dos_both_orb"//c_index_name(k) 778 call densityofstates(filename,value, 779 > eig,dbl_mb(weight(1)),ne(1), 780 > 1.0d0,alpha,npoints,emin,emax) 781 end if 782 if (ispin.eq.2) then 783 if (flag.eq.0) 784 > filename = "mulliken_fdos_alpha_orb"//c_index_name(k) 785 if (flag.eq.1) 786 > filename = "mulliken_vdos_alpha_orb"//c_index_name(k) 787 call densityofstates(filename,.false., 788 > eig,dbl_mb(weight(1)),ne(1), 789 > 1.0d0,alpha,npoints,emin,emax) 790 if (flag.eq.0) value = .true. 791 if (flag.eq.1) value = .false. 792 filename= "mulliken_dos_alpha_orb"//c_index_name(k) 793 call densityofstates(filename,value, 794 > eig,dbl_mb(weight(1)),ne(1), 795 > 1.0d0,alpha,npoints,emin,emax) 796 797 if (flag.eq.0) 798 > filename = "mulliken_fdos_beta_orb"//c_index_name(k) 799 if (flag.eq.1) 800 > filename = "mulliken_vdos_beta_orb"//c_index_name(k) 801 call densityofstates(filename,.false., 802 > eig(1+ne(1)),dbl_mb(weight(1)+ne(1)),ne(2), 803 > -1.0d0,alpha,npoints,emin,emax) 804 if (flag.eq.0) value = .true. 805 if (flag.eq.1) value = .false. 806 filename= "mulliken_dos_beta_orb"//c_index_name(k) 807 call densityofstates(filename,value, 808 > eig(1+ne(1)),dbl_mb(weight(1)+ne(1)),ne(2), 809 > -1.0d0,alpha,npoints,emin,emax) 810 end if 811 end do 812 813 814* *** free heap *** 815 if(fixatoms) then 816 if (.not. BA_free_heap(h_actlist)) 817 > call errquit('h_actlist:error freeing heap memory',0, MA_ERR) 818 end if 819 820c* **** atom decomposition ***** 821c do I=1,nion 822c 823c call dcopy(nemax,0.0d0,0,dbl_mb(weight(1)),1) 824c DO SPIN=1,ISPIN 825c DO N=N1(SPIN),N2(SPIN) 826c DO L=0,int_mb(lorb(1)+ion_katm_qm(I)-1) 827c L1=L**2+1 828c L2=(L+1)**2 829c DO LL=L1,L2 830c dbl_mb(weight(1)+n-1)= 831c > dbl_mb(weight(1)+n-1)+ 832c > dbl_mb(A(1)+LL-1+(N-1)*36+(I-1)*36*nemax)**2 833c END DO 834c END DO 835c END DO 836c END DO 837c 838c if (ispin.eq.1) then 839c write(filename,1801) I 840c call densityofstates(filename, 841c > eig,dbl_mb(weight(1)),ne(1), 842c > 1.0d0,alpha,npoints,emin,emax) 843c end if 844c if (ispin.eq.2) then 845c write(filename,1802) I 846c call densityofstates(filename, 847c > eig,dbl_mb(weight(1)),ne(1), 848c > 1.0d0,alpha,npoints,emin,emax) 849c write(filename,1803) I 850c call densityofstates(filename, 851c > eig(1+ne(1)),dbl_mb(weight(1)+ne(1)),ne(2), 852c > -1.0d0,alpha,npoints,emin,emax) 853c end if 854c 1801 FORMAT('DOS_both_atom_',I4.4) 855c 1802 FORMAT('DOS_alpha_atom_',I4.4) 856c 1803 FORMAT('DOS_beta_atom_',I4.4) 857c end do 858 859 value = BA_pop_stack(coef(2)) 860 value = value.and.BA_pop_stack(weight(2)) 861 if (.not. value) 862 > call errquit('analysis: error freeing stack',0, MA_ERR) 863 864 end if !*** control_DOS *** 865 866 867 868 869 870* **** free heap space **** 871 1901 continue 872 if (.not.mulliken_kawai) call aorbs_end() 873 value = BA_free_heap(lorb(2)) 874 value = BA_free_heap(b0(2)) 875 value = BA_free_heap(a(2)) 876 value = BA_free_heap(total(2)) 877 value = BA_free_heap(sum(2)) 878 value = BA_free_heap(subtl(2)) 879 880 return 881 end 882 883 884