1!!----------------- Plot TDOS/PDOS/OPDOS 2!For .out or plain text file, only one type of spin MOs will be loaded and processed, and then we will not consider spin type 3!For .fch/.molden/.gms, etc., user can switch spin type anytime 4subroutine plotdos 5use defvar 6use util 7use function 8implicit real*8 (a-h,o-z) 9integer,parameter :: nfragmax=10 10integer,parameter :: num2Dpoints=200 !The number of points constituting the X-axis of 2D LDOS 11real*8 :: curvexpos(num1Dpoints),TDOScurve(num1Dpoints),OPDOScurve(num1Dpoints),PDOScurve(num1Dpoints,nfragmax),LDOScurve(num1Dpoints) 12real*8 :: LDOSxpos(num2Dpoints) 13!All ?DOSliney share DOSlinex(:) as X axis 14real*8,allocatable :: DOSlinex(:),TDOSliney(:),PDOSliney(:,:),OPDOSliney(:),LDOSliney(:) 15real*8,allocatable :: compfrag(:,:) !i,k element is the MPA composition of fragment k in MO i 16real*8,allocatable :: OPfrag12(:) !Overlap population between fragment 1 and 2 17real*8,allocatable :: LDOScomp(:) !Composition at a point of each orbital 18real*8,allocatable :: LDOSptscomp(:,:) !Composition of each MO, ipt in a given line 19real*8,allocatable :: LDOS2Dmap(:,:) !LDOS curve, ipt in a given line 20integer :: nfragDOS(nfragmax),iDOSwidth=1 21integer,allocatable :: fragDOS(:,:) !The index of basis functions in fragments. nfragDOS is the number of basis functions in them (0=undefined) 22real*8,pointer :: tmpmat(:,:) 23real*8 HOMOlevx(2),HOMOlevy(2) 24real*8,allocatable :: MOene_dos(:),MOocc_dos(:) !Using the ene/occ in this to plot DOS, the values are scaled when changing between a.u. and eV. The original MOene is remain unchanged 25integer,allocatable :: selorbarr(:) 26character clegend*960 !(10+2) lines * 80 character per line 27character unitstr*5,c200tmp*200,c80tmp*80 28character :: TDOSstring*80="TDOS",OPDOSstring*80="OPDOS" 29character(len=80), dimension(nfragmax) :: PDOSstring = [character(len=80):: "PDOS frag.1","PDOS frag.2","PDOS frag.3","PDOS frag.4", & 30 "PDOS frag.5","PDOS frag.6","PDOS frag.7","PDOS frag.8","PDOS frag.9","PDOS frag.10"] 31integer :: ishowPDOSline(nfragmax),ishowPDOScurve(nfragmax) 32integer :: iclrPDOS(nfragmax)=(/ 1,3,10,14,12,9,13,11,6,7 /) 33 34if (allocated(FWHM)) deallocate(FWHM) !Global array 35if (allocated(str)) deallocate(str) !Global array 36defFWHM=0.05D0 !Default FWHM 37ipopmethod=1 !The method calculated OPDOS, =1 Mulliken, =3 SCPA, stout-politzer is too bad so don't consider it 38ibroadfunc=2 !Default is Gaussian function 39scalecurve=0.1D0 !Multiply curves with this value 40enelow=-0.8D0 !Energy range, a.u. 41enehigh=0.2D0 42stepx=0.1D0 43stepy=2 44gauweigh=0.5D0 !The weight of Gaussian in Pseudo-Voigt function 45nfragDOS=0 46ishowTDOScurve=1 47ishowTDOSline=1 48ishowPDOSline=0 49ishowPDOScurve=0 50ishowOPDOScurve=0 51ishowOPDOSline=0 52ishowlegend=1 53ishowHOMOlev=1 54iunitx=1 55unitstr=" a.u." 56ispin=0 !restricted wavefunction 57if (wfntype==1.or.wfntype==4) ispin=1 !Unrestricted wavefunction, output alpha part by default 58iusersetYleft=0 !If user has set lower and upper range of Y axis by himself 59Yrightscafac=0.5D0 !Scale factor relative to left Y-axis of OPDOS (right Y-axis) 60yxratio=1D0 61 62ireadgautype=1 63if (ifiletype==0) then 64 !Read energy level information from text file, the first number in first row define how many energy levels 65 !in there, the second number in first row if equals to 1, means below data are only energies, if equals to 2, 66 !means both strength and FWHM also present. 67 open(10,file=filename,status="old") 68 call loclabel(10,"Gaussian",igauout) 69 rewind(10) 70 if (igauout==1) then 71 write(*,*) "This is Gaussian output file" 72 if (ireadgautype==1) then !Read energy level from Gaussian output 73 call loclabel(10,"NBsUse=") 74 read(10,*) c200tmp,nbasis 75 nmo=nbasis 76 allocate(MOene(nmo),MOocc(nmo),str(nmo),FWHM(nmo)) 77 call loclabel(10,"Orbital energies and kinetic energies",ifound) !First assume this is close-shell 78 if (ifound==1) then 79 read(10,*) 80 read(10,"(a)") c200tmp 81 do i=1,nbasis 82 read(10,"(a21)",advance="no") c200tmp 83 read(10,*) MOene(i) 84 MOocc(i)=0 85 if (index(c200tmp,'O')/=0) MOocc(i)=2 86 end do 87 call loclabel(10,"Orbital energies and kinetic energies (beta)",ifound) 88 if (ifound==1) then 89 where(MOocc==2) MOocc=1 90 write(*,*) "Read which type of orbitals? 1:alpha 2:beta" 91 read(*,*) inp 92 if (inp==2) then !Read beta energies, overlay read alpha counterpart 93 read(10,*) 94 read(10,*) 95 do i=1,nbasis 96 read(10,"(a21)",advance="no") c200tmp 97 read(10,*) MOene(i) 98 MOocc(i)=0 99 if (index(c200tmp,'O')/=0) MOocc(i)=1 100 end do 101 end if 102 end if 103 write(*,*) "Read orbital energy from the file" 104 else 105 write(*,"(a)") " Error: Cannot find orbital energies from this file, don't forget using pop=full keyword" 106 write(*,*) 107 return 108 end if 109 end if 110 str=1D0 111 FWHM=defFWHM 112 else !Plain text file 113 read(10,*) nmo,inp 114 allocate(MOene(nmo),MOocc(nmo),str(nmo),FWHM(nmo)) 115 if (inp==1) then 116 do imo=1,nmo 117 read(10,*) MOene(imo),MOocc(imo) 118 end do 119 str=1D0 120 FWHM=defFWHM 121 else if (inp==2) then 122 do imo=1,nmo 123 read(10,*) MOene(imo),MOocc(imo),str(imo),FWHM(imo) 124 end do 125 end if 126 end if 127 close(10) 128 allocate(DOSlinex(3*nmo),TDOSliney(3*nmo)) 129else if (allocated(CObasa)) then 130 allocate(str(nbasis),FWHM(nbasis)) !I assume the number of MOs taken into account is equal to nbasis 131 !Allocate all arrays that may be used, don't consider if they will actually be used, because memory consuming is very little 132 allocate(DOSlinex(3*nbasis),TDOSliney(3*nbasis),PDOSliney(3*nbasis,nfragmax),OPDOSliney(3*nbasis),LDOSliney(3*nbasis)) 133 allocate(compfrag(nbasis,nfragmax),OPfrag12(nbasis)) 134 allocate(fragDOS(nbasis,nfragmax+1)) !The last slot is used to exchange fragment 135 allocate(LDOScomp(nbasis)) 136 str=1D0 137 FWHM=defFWHM 138else 139 write(*,*) "Error: Your input file does not contain basis function information!" 140 write(*,*) "Press ENTER to return" 141 read(*,*) 142 return 143end if 144 145!======Set from where to where are active energy levels 146if (ispin==0) imoend=nmo !Text file or restricted .fch 147if (ispin/=0) imoend=nbasis !For unrestricted fch or Gaussian output file, from 1 to imoend is the length of one type of spin orbitals 148 149if (allocated(MOene_dos)) deallocate(MOene_dos,MOocc_dos) !MOene_dos is the working horse 150allocate(MOene_dos(nmo),MOocc_dos(nmo)) 151MOene_dos=MOene 152MOocc_dos=MOocc 153 154do while(.true.) !!!!! main loop 155idoPDOS=0 156if (any(nfragDOS>0)) idoPDOS=1 157idoOPDOS=0 158if (all(nfragDOS(1:2)>0)) idoOPDOS=1 159 160!Unknow text file doesn't contains wavefunction info, couldn't define fragment 161write(*,*) " ================ Plot density-of-states ===============" 162write(*,*) "-10 Return to main menu" 163write(*,*) "-5 Customize energy levels, occupations, strengths and FWHMs for specific MOs" 164write(*,*) "-4 Show all orbital information" 165write(*,*) "-3 Export energy levels, occupations, strengths and FWHMs to plain text file" 166if (allocated(CObasa)) write(*,*) "-1 Define fragments" 167if (idoOPDOS==1) then 168 write(*,*) "0 Draw TDOS+PDOS+OPDOS graph!" 169else if (idoPDOS==1) then 170 write(*,*) "0 Draw TDOS+PDOS graph!" 171else 172 write(*,*) "0 Draw TDOS graph!" !Reading text file can only draw spinless TDOS, because they impossible to define fragment 173end if 174if (ibroadfunc==1) write(*,*) "1 Select broadening function, current: Lorentzian" 175if (ibroadfunc==2) write(*,*) "1 Select broadening function, current: Gaussian" 176if (ibroadfunc==3) write(*,*) "1 Select broadening function, current: Pseudo-Voigt" 177write(*,"(a,f14.5,a,f14.5,a)") " 2 Set energy range, current:",enelow," to",enehigh,unitstr 178if (maxval(FWHM)==minval(FWHM)) then 179 write(*,"(a,f10.5,a)") " 3 Set full width at half maximum (FWHM), current:",FWHM(1),unitstr 180else 181 write(*,"(a,f10.5)") " 3 Set full width at half maximum (FWHM), current: Orbital dependent" 182end if 183write(*,"(a,f10.5)") " 4 Set scale ratio for DOS curve, current:",scalecurve 184if (ibroadfunc==3) write(*,"(a,f10.5)") " 5 Set Gaussian-weighting coefficient, current:",gauweigh 185if (ispin==1) write(*,*) "6 Switch spin, current: Alpha" 186if (ispin==2) write(*,*) "6 Switch spin, current: Beta" 187if (allocated(CObasa)) then 188 if (ipopmethod==1) write(*,*) "7 Switch method for calculating PDOS, current: Mulliken" 189 if (ipopmethod==3) write(*,*) "7 Switch method for calculating PDOS, current: SCPA" 190end if 191write(*,*) "8 Switch unit between a.u. and eV, current:"//unitstr 192write(*,*) "10 Draw local DOS for a point" 193write(*,*) "11 Draw local DOS along a line" 194 195read(*,*) isel 196 197if (isel==-10) then 198 exit 199else if (isel==-5) then 200 do while(.true.) 201 write(*,*) 202 write(*,*) "0 Return" 203 write(*,*) "1 Set orbital energies for specific orbitals" 204 write(*,*) "2 Set occupation numbers for specific orbitals" 205 write(*,*) "3 Set strengths for specific orbitals" 206 write(*,*) "4 Set FWHMs for specific orbitals" 207 read(*,*) isel 208 209 if (isel==0) then 210 exit 211 else 212 write(*,"(a)") " Input orbital indices. e.g. 1,3-6,8,10-11 means the orbital 1,3,4,5,6,8,10,11 will be selected" 213 read(*,"(a)") c200tmp 214 call str2arr(c200tmp,nselorb) 215 if (allocated(selorbarr)) deallocate(selorbarr) 216 allocate(selorbarr(nselorb)) 217 call str2arr(c200tmp,nselorb,selorbarr) 218 if (isel==1) then 219 write(*,*) "Set their energies to which value? e.g. -0.13" 220 write(*,*) "Note: The value should be given in"//unitstr 221 read(*,*) enetmp 222 if (ispin==2) selorbarr=selorbarr+nbasis 223 do imoidx=1,nselorb 224 imo=selorbarr(imoidx) 225 MOene_dos(imo)=enetmp 226 end do 227 else if (isel==2) then 228 write(*,*) "Set their occupation numbers to which value? e.g. 2.0" 229 read(*,*) occtmp 230 if (ispin==2) selorbarr=selorbarr+nbasis 231 do imoidx=1,nselorb 232 imo=selorbarr(imoidx) 233 MOocc_dos(imo)=occtmp 234 end do 235 else if (isel==3) then 236 write(*,*) "Set their strength to which value? e.g. 1.0" 237 read(*,*) strtmp 238 do imoidx=1,nselorb 239 imo=selorbarr(imoidx) 240 str(imo)=strtmp 241 end do 242 else if (isel==4) then 243 write(*,*) "Set their FWHM to which value? e.g. 0.05" 244 read(*,*) FWHMtmp 245 do imoidx=1,nselorb 246 imo=selorbarr(imoidx) 247 FWHM(imo)=FWHMtmp 248 end do 249 end if 250 end if 251 end do 252else if (isel==-4) then 253 iFermi=0 254 if (ispin==1) write(*,*) "Below orbitals are Alpha type" 255 if (ispin==2) write(*,*) "Below orbitals are Beta type" 256 do imo=1,imoend 257 irealmo=imo 258 if (ispin==2) irealmo=imo+nbasis 259 write(*,"('#',i5,' Energy(',a,'):',f14.6,' Occ:',f8.5,' Str:',f8.5,' FWHM:',f8.5)") imo,trim(unitstr),MOene_dos(irealmo),MOocc_dos(irealmo),str(imo),FWHM(imo) 260 if (imo>1) then 261 if (MOocc_dos(irealmo)==0D0.and.MOocc_dos(irealmo-1)/=0D0) iFermi=irealmo-1 262 end if 263 end do 264 if (iFermi/=0) write(*,"(' Fermi energy level:',f12.6,a)") MOene_dos(iFermi),unitstr 265 write(*,*) 266else if (isel==-3) then 267 open(10,file="orginfo.txt",status="replace") 268 write(10,"(2i6)") imoend,2 269 do imo=1,imoend 270 irealmo=imo 271 if (ispin==2) irealmo=imo+nbasis 272 write(10,"(f14.6,3f12.6)") MOene_dos(irealmo),MOocc_dos(irealmo),str(imo),FWHM(imo) 273 end do 274 close(10) 275 write(*,"(a)") " The energy levels, occupation numbers, strengths, FWHMs have been exported to orginfo.txt in current directory, & 276 you can modify it and then load it into Multiwfn again." 277 write(*,*) "Note: The energy unit of energy levels and FWHMs in this file is in"//unitstr 278 write(*,*) 279else if (isel==-1) then 280 write(*,*) " ----------------- Define fragments -----------------" 281 write(*,"(a)") " Note: Up to 10 fragments can be defined for plotting PDOS, but OPDOS will only be plotted for fragments 1 and 2" 282 do while(.true.) 283 do ifrag=1,nfragmax 284 if (nfragDOS(ifrag)==0) then 285 write(*,"(' Fragment',i5,', has not been defined')") ifrag 286 else 287 write(*,"(' Fragment',i5,', the number of basis functions:',i6)") ifrag,nfragDOS(ifrag) 288 end if 289 end do 290 write(*,*) "Input fragment index to define it, e.g. 2" 291 write(*,*) "Input a negative index can unset the fragment, e.g. -2" 292 write(*,*) "Input two indices can exchange the two fragments, e.g. 1,4" 293 write(*,*) "Input ""e"" can export current fragment setting to DOSfrag.txt in current folder" 294 write(*,*) "Input ""i"" can import fragment setting from DOSfrag.txt in current folder" 295 write(*,*) "To return to the last menu, input 0" 296 read(*,"(a)") c80tmp 297 if (index(c80tmp(1:len_trim(c80tmp)),' ')/=0.or.c80tmp==" ") then 298 write(*,*) "Input error!" 299 write(*,*) 300 else if (c80tmp=='e') then 301 open(10,file="DOSfrag.txt",status="replace") 302 do ifrag=1,nfragmax 303 write(10,*) 304 write(10,"(' #Fragment:',i4,' nbasis:',i8)") ifrag,nfragDOS(ifrag) 305 write(10,"(8i8)") fragDOS(1:nfragDOS(ifrag),ifrag) 306 end do 307 close(10) 308 write(*,*) "Export finished!" 309 write(*,*) 310 else if (c80tmp=='i') then 311 open(10,file="DOSfrag.txt",status="old") 312 do ifrag=1,nfragmax 313 read(10,*) 314 read(10,*) c80tmp,inouse,c80tmp,nfragDOS(ifrag) 315 read(10,"(8i8)") fragDOS(1:nfragDOS(ifrag),ifrag) 316 end do 317 close(10) 318 write(*,*) "Import finished!" 319 write(*,*) 320 else if (index(c80tmp,',')==0) then !Set specific fragment 321 read(c80tmp,*) ifragsel 322 if (ifragsel==0) then 323 exit 324 else if (ifragsel>nfragmax) then 325 write(*,*) "ERROR: The index exceeded upper limit!" 326 else if (ifragsel<0) then 327 nfragDOS(abs(ifragsel))=0 328 else !deffrag routine is only able to deal with global array frag1/2, so we use frag1 as intermediate array 329 allocate(frag1(nfragDOS(ifragsel))) 330 frag1(:)=fragDOS(1:nfragDOS(ifragsel),ifragsel) 331 call deffrag(1) 332 if (allocated(frag1)) then 333 nfragDOS(ifragsel)=size(frag1) 334 fragDOS(1:nfragDOS(ifragsel),ifragsel)=frag1(:) 335 deallocate(frag1) 336 else 337 nfragDOS(ifragsel)=0 338 end if 339 end if 340 else !Exchange fragments 341 read(c80tmp,*) ifragsel,jfragsel 342 ntmp=nfragDOS(jfragsel) 343 nfragDOS(jfragsel)=nfragDOS(ifragsel) 344 nfragDOS(ifragsel)=ntmp 345 fragDOS(:,nfragmax+1)=fragDOS(:,jfragsel) 346 fragDOS(:,jfragsel)=fragDOS(:,ifragsel) 347 fragDOS(:,ifragsel)=fragDOS(:,nfragmax+1) 348 end if 349 end do 350else if (isel==1) then 351 write(*,*) "1 Lorentzian" 352 write(*,*) "2 Gaussian" 353 write(*,*) "3 Pseudo-Voigt" 354 read(*,*) ibroadfunc 355else if (isel==2) then 356 if (iunitx==1) then 357 write(*,*) "Input lower, upper limits and stepsize between legends (in a.u.)" 358 write(*,*) "e.g. -1.5,0.2,0.3" 359 else if (iunitx==2) then 360 write(*,*) "Input lower, upper limits and stepsize between legends (in eV)" 361 write(*,*) "e.g. -20,5,2" 362 end if 363 read(*,*) enelow,enehigh,stepx 364else if (isel==3) then 365 write(*,*) "Input a value" 366 read(*,*) FWHMtmp 367 if (FWHMtmp<0D0) write(*,*) "Error: The value should larger than zero, input again" 368 FWHM=FWHMtmp 369else if (isel==4) then 370 write(*,*) "Input a value" 371 read(*,*) scalecurve 372 if (scalecurve<0D0) write(*,*) "Error: The value should larger than zero, input again" 373else if (isel==5) then 374 write(*,*) "Input a value" 375 read(*,*) gauweigh 376 if (gauweigh<0D0) write(*,*) "Error: The value should larger than zero, input again" 377else if (isel==6) then 378 if (ispin==1) then 379 ispin=2 380 else 381 ispin=1 382 end if 383else if (isel==7) then 384 if (ipopmethod==1) then 385 ipopmethod=3 386 else 387 ipopmethod=1 388 end if 389else if (isel==8) then 390 if (iunitx==1) then !a.u.->eV 391 iunitx=2 392 MOene_dos=MOene_dos*au2eV 393 FWHM=FWHM*au2eV 394 enelow=enelow*au2eV 395 enehigh=enehigh*au2eV 396 unitstr=" eV" 397 !After change the unit, in principle, the curve (and hence Y-range) will be automatically reduced by 27.2114.& 398 !str should also be reduced by 27.2114 so that the discrete line can be properly shown in the graph range & 399 !To compensate the reduce of str, scalecurve thus be augmented by corresponding factor 400 str=str/au2eV 401 scalecurve=scalecurve*au2eV 402 stepx=stepx*au2eV 403 stepy=stepy/au2eV 404 else !eV->a.u. 405 iunitx=1 406 MOene_dos=MOene_dos/au2eV 407 FWHM=FWHM/au2eV 408 enelow=enelow/au2eV 409 enehigh=enehigh/au2eV 410 unitstr=" a.u." 411 str=str*au2eV 412 scalecurve=scalecurve/au2eV 413 stepx=stepx/au2eV 414 stepy=stepy*au2eV 415 end if 416 417 418else if (isel==0.or.isel==10) then 419 420 if (isel==10) then 421 write(*,*) "Input coordinate (in Bohr), e.g. 1.0,1.5,0.2" 422 read(*,*) x,y,z 423 end if 424 425 !======Generate fragment composition and overlap population======= 426 tmpmat=>cobasa 427 if (ispin==2) tmpmat=>cobasb 428 429 !Reset display setting 430 ishowPDOScurve=0 431 ishowPDOSline=0 432 do ifrag=1,nfragmax 433 if (nfragDOS(ifrag)>0) then 434 ishowPDOScurve(ifrag)=1 435 ishowPDOSline(ifrag)=1 436 end if 437 end do 438 ishowOPDOScurve=0 439 ishowOPDOSline=0 440 if (idoOPDOS==1) then 441 ishowOPDOScurve=1 442 ishowOPDOSline=1 443 end if 444 ishowLDOScurve=0 445 ishowLDOSline=0 446 if (isel==10) then 447 ishowLDOScurve=1 448 ishowLDOSline=1 449 end if 450 451 if (idoPDOS==1.or.isel==10) then 452 write(*,*) "Calculating orbital composition, please wait..." 453 do istart=1,nbasis 454 enetmp=MOene_dos(istart) 455 if (ispin==2) enetmp=MOene_dos(istart+nbasis) 456 if (enetmp>enelow-3*FWHM(istart)) exit 457 end do 458 do iend=nbasis,1,-1 459 enetmp=MOene_dos(iend) 460 if (ispin==2) enetmp=MOene_dos(iend+nbasis) 461 if (enetmp<enehigh+3*FWHM(iend)) exit 462 end do 463 write(*,"(' Note: MOs from',i7,' to',i7,' are taken into account')") istart,iend 464 compfrag=0 465 LDOScomp=0 466 OPfrag12=0 467 if (idoPDOS==1) then 468nthreads=getNThreads() 469!$OMP PARALLEL DO SHARED(compfrag,OPfrag12) PRIVATE(ifrag,imo,allsqr,i,j,ibas,jbas) schedule(dynamic) NUM_THREADS(nthreads) 470 do imo=istart,iend !Cycle each orbital, don't use nmo, because for unrestricted wavefunction nmo=2*nbasis 471 if (ipopmethod==3) allsqr=sum(tmpmat(:,imo)**2) 472 do ifrag=1,nfragmax 473 if (nfragDOS(ifrag)==0) cycle 474 do i=1,nfragDOS(ifrag) !Cycle each basis in the fragment 475 ibas=fragDOS(i,ifrag) 476 if (ipopmethod==3) then !SCPA 477 compfrag(imo,ifrag)=compfrag(imo,ifrag)+tmpmat(ibas,imo)**2/allsqr 478 else !Mulliken 479 do jbas=1,nbasis !Cycle all basis, included inner&external cross term and local term (when ibas==jbas) 480 compfrag(imo,ifrag)=compfrag(imo,ifrag)+tmpmat(ibas,imo)*tmpmat(jbas,imo)*Sbas(ibas,jbas) 481 end do 482 end if 483 end do 484 end do 485 !Calculate Overlap population between frag 1&2 486 if (idoOPDOS==1) then 487 do i=1,nfragDOS(1) 488 ibas=fragDOS(i,1) 489 do j=1,nfragDOS(2) 490 jbas=fragDOS(j,2) 491 OPfrag12(imo)=OPfrag12(imo)+2*tmpmat(ibas,imo)*tmpmat(jbas,imo)*Sbas(ibas,jbas) 492 end do 493 end do 494 end if 495 end do 496!$OMP END PARALLEL DO 497 end if 498 if (isel==10) then !Calculate LDOS for a point 499 do imo=istart,iend !Cycle each orbital 500 LDOScomp(imo)=fmo(x,y,z,imo)**2 501 end do 502 end if 503 end if 504 505 506 !======Set X position of curves========== 507 enestep=(enehigh-enelow)/(num1Dpoints-1) 508 do i=1,num1Dpoints 509 curvexpos(i)=enelow+(i-1)*enestep 510 end do 511 512 !======Generate energy levels line======= 513 TDOSliney=0 514 LDOSliney=0 515 OPDOSliney=0 516 do imo=1,imoend 517 inow=3*(imo-1) 518 irealmo=imo 519 if (ispin==2) irealmo=imo+nbasis 520 DOSlinex(inow+1:inow+3)=MOene_dos(irealmo) 521 if (isel==0) then 522 TDOSliney(inow+1)=0D0 523 TDOSliney(inow+2)=str(imo) 524 TDOSliney(inow+3)=0D0 525 do ifrag=1,nfragmax 526 if (nfragDOS(ifrag)>0) then 527 PDOSliney(inow+1,ifrag)=0D0 528 PDOSliney(inow+2,ifrag)=str(imo)*compfrag(imo,ifrag) 529 PDOSliney(inow+3,ifrag)=0D0 530 end if 531 end do 532 if (idoOPDOS==1) then 533 OPDOSliney(inow+1)=0D0 534 OPDOSliney(inow+2)=str(imo)*OPfrag12(imo) 535 OPDOSliney(inow+3)=0D0 536 end if 537 else if (isel==10) then 538 LDOSliney(inow+1)=0D0 539 LDOSliney(inow+2)=str(imo)*LDOScomp(imo) 540 LDOSliney(inow+3)=0D0 541 end if 542 end do 543 544 !======Generate DOS curve======= 545 TDOScurve=0D0 546 PDOScurve=0D0 547 OPDOScurve=0D0 548 LDOScurve=0D0 549 if (ibroadfunc==1.or.ibroadfunc==3) then !Lorentzian function, see http://mathworld.wolfram.com/LorentzianFunction.html 550 do imo=1,imoend !Cycle each orbital 551 irealmo=imo 552 if (ispin==2) irealmo=imo+nbasis 553 preterm=str(imo)*0.5D0/pi*FWHM(imo) 554 do ipoint=1,num1Dpoints !Broaden imo as curve 555 tmp=preterm/( (curvexpos(ipoint)-MOene_dos(irealmo))**2+0.25D0*FWHM(imo)**2 ) 556 if (isel==0) then 557 TDOScurve(ipoint)=TDOScurve(ipoint)+tmp 558 do ifrag=1,nfragmax 559 if (nfragDOS(ifrag)>0) PDOScurve(ipoint,ifrag)=PDOScurve(ipoint,ifrag)+tmp*compfrag(imo,ifrag) 560 end do 561 if (idoOPDOS==1) OPDOScurve(ipoint)=OPDOScurve(ipoint)+tmp*OPfrag12(imo) 562 else if (isel==10) then 563 LDOScurve(ipoint)=LDOScurve(ipoint)+tmp*LDOScomp(imo) 564 end if 565 end do 566 end do 567 end if 568 if (ibroadfunc==2.or.ibroadfunc==3) then 569 if (ibroadfunc==3) TDOScurve=(1-gauweigh)*TDOScurve 570 do imo=1,imoend !Cycle each orbital 571 irealmo=imo 572 if (ispin==2) irealmo=imo+nbasis 573 gauss_c=FWHM(imo)/2D0/sqrt(2*dlog(2D0)) 574 gauss_a=str(imo)/(gauss_c*sqrt(2D0*pi)) 575 do ipoint=1,num1Dpoints !Broaden imo as curve 576 !Gaussian function, see http://en.wikipedia.org/wiki/Gaussian_function 577 tmp=gauss_a*dexp( -(curvexpos(ipoint)-MOene_dos(irealmo))**2/(2*gauss_c**2) ) 578 if (ibroadfunc==3) tmp=gauweigh*tmp !Combine Lorentizan and Gaussian function 579 if (isel==0) then 580 TDOScurve(ipoint)=TDOScurve(ipoint)+tmp 581 do ifrag=1,nfragmax 582 if (nfragDOS(ifrag)>0) PDOScurve(ipoint,ifrag)=PDOScurve(ipoint,ifrag)+tmp*compfrag(imo,ifrag) 583 end do 584 if (idoOPDOS==1) OPDOScurve(ipoint)=OPDOScurve(ipoint)+tmp*OPfrag12(imo) 585 else if (isel==10) then 586 LDOScurve(ipoint)=LDOScurve(ipoint)+tmp*LDOScomp(imo) 587 end if 588 end do 589 end do 590 end if 591 TDOScurve=TDOScurve*scalecurve 592 PDOScurve=PDOScurve*scalecurve 593 OPDOScurve=OPDOScurve*scalecurve 594 LDOScurve=LDOScurve*scalecurve 595 596 idraw=1 597 isavepic=0 598 if (iusersetYleft==0) then !Y axis range was not set by user, we use recommended value 599 if (isel==0) then 600 yupperleft=1.1D0*maxval(TDOScurve) 601 ylowerleft=0 602 if (idoPDOS==1) ylowerleft=minval(PDOScurve(:,:)) !PDOS may be negative 603 if (idoOPDOS==1) ylowerleft=-yupperleft/2 !OPDOS may be large negative value 604 if (ylowerleft>0) ylowerleft=0D0 !Don't allow lower plotting limit >0 605 stepy=nint((yupperleft-ylowerleft)*10)/100D0 606 else if (isel==10) then 607 ylowerleft=0 608 yupperleft=1.1D0*maxval(LDOScurve) 609 stepy=(yupperleft-ylowerleft)/10 610 end if 611 end if 612 ylowerright=ylowerleft*Yrightscafac !Lower and upper limit for OPDOS 613 yupperright=yupperleft*Yrightscafac 614 615 do while(.true.) 616 if (idraw==1.and.isilent==0) then 617 !======Draw DOS now======= 618 if (idoOPDOS==1) then 619 else 620 end if 621 numleg=0 622 ileg=1 623 624 if (isel==0) then 625 !Set legends 626 if (ishowTDOScurve==1.or.ishowTDOSline==1) numleg=numleg+1 627 do ifrag=1,nfragmax 628 if (nfragDOS(ifrag)==0) cycle 629 if (ishowPDOScurve(ifrag)==1.or.ishowPDOSline(ifrag)==1) numleg=numleg+1 630 end do 631 if (ishowOPDOScurve==1.or.ishowOPDOSline==1) numleg=numleg+1 632 633 634 !Draw a vertical dashed line to highlight HOMO level 635 if (ishowHOMOlev==1) then 636 do imo=1,imoend 637 irealmo=imo 638 if (ispin==2) irealmo=imo+nbasis 639 if (imo>1) then 640 if (MOocc_dos(irealmo)==0D0.and.MOocc_dos(irealmo-1)/=0D0) iFermi=irealmo-1 641 end if 642 end do 643 HOMOlevx=MOene_dos(iFermi) 644 HOMOlevy(1)=ylowerleft 645 HOMOlevy(2)=yupperleft 646 end if 647 648 !Draw TDOS 649 if (ishowTDOScurve==1.or.ishowTDOSline==1) then 650 ileg=ileg+1 651 end if 652 653 !Draw PDOS of each defined fragment 654 do ifrag=1,nfragmax 655 if (nfragDOS(ifrag)>0) then 656 if (ishowPDOScurve(ifrag)==1.or.ishowPDOSline(ifrag)==1) then 657 ileg=ileg+1 658 end if 659 end if 660 end do 661 662 !Draw OPDOS 663 if (ishowOPDOScurve==1.or.ishowOPDOSline==1) then 664 end if 665 666 !Draw LDOS 667 else if (isel==10) then 668 end if 669 670 if (isavepic==1) write(*,*) "Graphic file has been saved to current folder with ""DISLIN"" prefix" 671 end if 672 idraw=0 673 674 !======Submenu======= 675 write(*,*) 676 write(*,*) " ======== Post-process menu ========" 677 if (isel==0) then !T/P/OPDOS 678 write(*,*) "0 Return" 679 write(*,*) "1 Show graph again" 680 write(*,*) "2 Save the picture to current folder" 681 write(*,*) "3 Export curve data to plain text file in current folder" 682 if (idoPDOS==1) then 683 write(*,"(' 4 Set Y-axis range and step for TDOS+PDOS, current:',f8.2,' to',f8.2,',',f6.2)") ylowerleft,yupperleft,stepy 684 else 685 write(*,"(' 4 Set Y-axis range and step for TDOS, current:',f8.2,' to',f8.2,',',f6.2)") ylowerleft,yupperleft,stepy 686 end if 687 if (ishowTDOScurve==1) write(*,*) "5 Toggle showing TDOS curve, current: Yes" 688 if (ishowTDOScurve==0) write(*,*) "5 Toggle showing TDOS curve, current: No" 689 if (ishowTDOSline==1) write(*,*) "6 Toggle showing TDOS line, current: Yes" 690 if (ishowTDOSline==0) write(*,*) "6 Toggle showing TDOS line, current: No" 691 if (idoPDOS==1) then 692 write(*,*) "7 Toggle showing PDOS curves" 693 write(*,*) "8 Toggle showing PDOS lines" 694 end if 695 if (idoOPDOS==1) then 696 if (ishowOPDOScurve==1) write(*,*) "9 Toggle showing OPDOS curve, current: Yes" 697 if (ishowOPDOScurve==0) write(*,*) "9 Toggle showing OPDOS curve, current: No" 698 if (ishowOPDOSline==1) write(*,*) "10 Toggle showing OPDOS line, current: Yes" 699 if (ishowOPDOSline==0) write(*,*) "10 Toggle showing OPDOS line, current: No" 700 end if 701 if (idoPDOS==1) write(*,*) "11 Set color for PDOS curves and lines" 702 if (ishowlegend==1) write(*,*) "13 Toggle showing legend, current: Yes" 703 if (ishowlegend==0) write(*,*) "13 Toggle showing legend, current: No" 704 if (idoOPDOS==1) write(*,"(a,f10.5)") " 14 Set scale factor of Y-axis for OPDOS, current:",Yrightscafac 705 write(*,*) "15 Toggle showing vertical dashed line to highlight HOMO level" 706 write(*,*) "16 Set the texts in the legends" 707 write(*,"(a,i3)") " 17 Set width of lines and curves, current:",iDOSwidth 708 read(*,*) isel2 709 710 if (isel2==0) then 711 exit 712 else if (isel2==1) then 713 idraw=1 714 isavepic=0 715 else if (isel2==2) then 716 idraw=1 717 isavepic=1 718 else if (isel2==3) then 719 open(10,file="DOS_line.txt",status="replace") 720 if (idoOPDOS==1) then 721 do i=1,3*imoend 722 write(10,"(f10.5,12f12.6)") DOSlinex(i),TDOSliney(i),OPDOSliney(i),PDOSliney(i,:) 723 end do 724 else if (idoPDOS==1) then 725 do i=1,3*imoend 726 write(10,"(f10.5,11f12.6)") DOSlinex(i),TDOSliney(i),PDOSliney(i,:) 727 end do 728 else 729 do i=1,3*imoend 730 write(10,"(f10.5,f12.6)") DOSlinex(i),TDOSliney(i) 731 end do 732 end if 733 close(10) 734 open(10,file="DOS_curve.txt",status="replace") 735 if (idoOPDOS==1) then 736 do i=1,num1Dpoints 737 write(10,"(f10.5,12f12.6)") curvexpos(i),TDOScurve(i),OPDOScurve(i),PDOScurve(i,:) 738 end do 739 else if (idoPDOS==1) then 740 do i=1,num1Dpoints 741 write(10,"(f10.5,11f12.6)") curvexpos(i),TDOScurve(i),PDOScurve(i,:) 742 end do 743 else 744 do i=1,num1Dpoints 745 write(10,"(f10.5,f12.6)") curvexpos(i),TDOScurve(i) 746 end do 747 end if 748 close(10) 749 write(*,*) "Curve data have been written to DOS_curve.txt in current folder" 750 write(*,*) "Discrete line data have been written to DOS_line.txt in current folder" 751 write(*,*) "Column 1: Energy ("//trim(unitstr)//")" 752 write(*,*) "Column 2: TDOS" 753 if (idoOPDOS==1) then 754 write(*,*) "Column 3: OPDOS" 755 write(*,*) "Column 4-13: PDOS of fragment 1-10" 756 else if (idoPDOS==1) then 757 write(*,*) "Column 3-12: PDOS of fragment 1-10" 758 end if 759 else if (isel2==4) then 760 write(*,*) "Input lower and upper limit as well as stepsize, e.g. 0.0,2.4,0.3" 761 read(*,*) ylowerleft,yupperleft,stepy 762 iusersetYleft=1 763 ylowerright=ylowerleft*Yrightscafac !Lower and upper limit for OPDOS. Set it here aims for immediately make effect 764 yupperright=yupperleft*Yrightscafac 765 else if (isel2==5) then 766 if (ishowTDOScurve==0) then 767 ishowTDOScurve=1 768 else 769 ishowTDOScurve=0 770 end if 771 else if (isel2==6) then 772 if (ishowTDOSline==0) then 773 ishowTDOSline=1 774 else 775 ishowTDOSline=0 776 end if 777 else if (isel2==7) then 778 do while(.true.) 779 write(*,*) 780 write(*,*) "0 Return" 781 do ifrag=1,nfragmax 782 if (nfragDOS(ifrag)>0) then 783 if (ishowPDOScurve(ifrag)==1) write(*,"(i2,' Toggle showing PDOS curve of fragment',i3,', current: Yes')") ifrag,ifrag 784 if (ishowPDOScurve(ifrag)==0) write(*,"(i2,' Toggle showing PDOS curve of fragment',i3,', current: No')") ifrag,ifrag 785 end if 786 end do 787 read(*,*) iselfrag 788 if (iselfrag==0) exit 789 if (ishowPDOScurve(iselfrag)==0) then 790 ishowPDOScurve(iselfrag)=1 791 else 792 ishowPDOScurve(iselfrag)=0 793 end if 794 end do 795 else if (isel2==8) then 796 do while(.true.) 797 write(*,*) 798 write(*,*) "0 Return" 799 do ifrag=1,nfragmax 800 if (nfragDOS(ifrag)>0) then 801 if (ishowPDOSline(ifrag)==1) write(*,"(i2,' Toggle showing PDOS line of fragment',i3,', current: Yes')") ifrag,ifrag 802 if (ishowPDOSline(ifrag)==0) write(*,"(i2,' Toggle showing PDOS line of fragment',i3,', current: No')") ifrag,ifrag 803 end if 804 end do 805 read(*,*) iselfrag 806 if (iselfrag==0) exit 807 if (ishowPDOSline(iselfrag)==0) then 808 ishowPDOSline(iselfrag)=1 809 else 810 ishowPDOSline(iselfrag)=0 811 end if 812 end do 813 else if (isel2==9) then 814 if (ishowOPDOScurve==0) then 815 ishowOPDOScurve=1 816 else 817 ishowOPDOScurve=0 818 end if 819 else if (isel2==10) then 820 if (ishowOPDOSline==0) then 821 ishowOPDOSline=1 822 else 823 ishowOPDOSline=0 824 end if 825 else if (isel2==11) then 826 do while(.true.) 827 write(*,*) 828 write(*,*) "0 Return" 829 do ifrag=1,nfragmax 830 if (nfragDOS(ifrag)>0) write(*,"(' Set color for fragment',i3,', current: ',a)") ifrag,colorname(iclrPDOS(ifrag)) 831 end do 832 read(*,*) iselfrag 833 if (iselfrag<0.or.iselfrag>nfragmax) then 834 write(*,*) "ERROR: Index exceeded valid range!" 835 else if (iselfrag==0) then 836 exit 837 else 838 write(*,*) "Use which color?" 839 write(*,*) "1/2/3/4/5 = Red/Green/Blue/White/Black" 840 write(*,*) "6/7/8/9/10 = Gray/Cyan/Yellow/Orange/Magenta" 841 write(*,*) "11/12/13/14 = Crimson/Dark green/Purple/Brown" 842 read(*,*) iclrPDOS(iselfrag) 843 end if 844 end do 845 else if (isel2==13) then 846 if (ishowlegend==0) then 847 ishowlegend=1 848 else 849 ishowlegend=0 850 end if 851 else if (isel2==14) then 852 write(*,*) "Input scale factor, e.g. 0.15" 853 read(*,*) Yrightscafac 854 ylowerright=ylowerleft*Yrightscafac !Lower and upper limit for OPDOS. Set it here aims for immediately make effect 855 yupperright=yupperleft*Yrightscafac 856 else if (isel2==15) then 857 if (ishowHOMOlev==0) then 858 ishowHOMOlev=1 859 else 860 ishowHOMOlev=0 861 end if 862 else if (isel2==16) then 863 do while(.true.) 864 write(*,*) "Select a term to change its legend, e.g. 3" 865 write(*,"(' -2 TDOS, current text: ',a)") trim(TDOSstring) 866 if (idoOPDOS==1) write(*,"(' -1 OPDOS, current text: ',a)") trim(OPDOSstring) 867 write(*,*) " 0 Return" 868 do ifrag=1,nfragmax 869 if (nfragDOS(ifrag)>0) write(*,"(i3,' PDOS of frag',i2,', current text: ',a)") ifrag,ifrag,trim(PDOSstring(ifrag)) 870 end do 871 read(*,*) iseltmp 872 if (iseltmp==0) exit 873 write(*,*) "Input text for the legend" 874 read(*,"(a)") c80tmp 875 if (iseltmp==-2) then 876 TDOSstring=c80tmp 877 else if (iseltmp==-1) then 878 OPDOSstring=c80tmp 879 else 880 PDOSstring(iseltmp)=c80tmp 881 end if 882 end do 883 else if (isel2==17) then 884 write(*,*) "Input width of lines and curves, e.g. 5" 885 read(*,*) iDOSwidth 886 end if 887 888 else if (isel==10) then !LDOS in 1D 889 write(*,*) "0 Return" 890 write(*,*) "1 Show graph again" 891 write(*,*) "2 Save graphical file of the DOS map in current folder" 892 write(*,*) "3 Export curve data to plain text file in current folder" 893 if (ishowLDOSline==1) write(*,*) "6 Toggle showing LDOS line, current: Yes" 894 if (ishowLDOSline==0) write(*,*) "6 Toggle showing LDOS line, current: No" 895 read(*,*) isel2 896 897 if (isel2==0) then 898 exit 899 else if (isel2==1) then 900 idraw=1 901 isavepic=0 902 else if (isel2==2) then 903 idraw=1 904 isavepic=1 905 else if (isel2==3) then 906 open(10,file="LDOS_line.txt",status="replace") 907 do i=1,3*imoend 908 write(10,"(f10.5,1PE15.6)") DOSlinex(i),LDOSliney(i) 909 end do 910 close(10) 911 open(10,file="LDOS_curve.txt",status="replace") 912 do i=1,num1Dpoints 913 write(10,"(f10.5,1PE15.6)") curvexpos(i),LDOScurve(i) 914 end do 915 close(10) 916 write(*,*) "Curve data have been written to LDOS_curve.txt in current folder" 917 write(*,*) "Discrete line data have been written to LDOS_line.txt in current folder" 918 write(*,*) "Column 1: Energy ("//trim(unitstr)//")" 919 write(*,*) "Column 2: LDOS" 920 else if (isel2==6) then 921 if (ishowLDOSline==0) then 922 ishowLDOSline=1 923 else 924 ishowLDOSline=0 925 end if 926 end if 927 end if 928 end do 929 930 931!Plot local DOS along a line (2D color-filled map) 932else if (isel==11) then 933 if (allocated(LDOS2Dmap)) deallocate(LDOS2Dmap,LDOSptscomp) 934 write(*,*) "Input the starting point (in Bohr), e.g. 1.0,0,0.2" 935 read(*,*) orgx,orgy,orgz 936 write(*,*) "Input the end point (in Bohr), e.g. 2.0,0,0.4" 937 read(*,*) endx,endy,endz 938 write(*,*) "Input the number of points along the line" 939 read(*,*) numLDOSpt 940 allocate(LDOS2Dmap(num2Dpoints,numLDOSpt),LDOSptscomp(nbasis,numLDOSpt)) 941 write(*,*) "Calculating orbital composition, please wait..." 942 do istart=1,nbasis 943 enetmp=MOene_dos(istart) 944 if (ispin==2) enetmp=MOene_dos(istart+nbasis) 945 if (enetmp>enelow-3*FWHM(istart)) exit 946 end do 947 do iend=nbasis,1,-1 948 enetmp=MOene_dos(iend) 949 if (ispin==2) enetmp=MOene_dos(iend+nbasis) 950 if (enetmp<enehigh+3*FWHM(iend)) exit 951 end do 952! istart=1 953! iend=nbasis 954 write(*,"(' Note: MOs from',i7,' to',i7,' are taken into account')") istart,iend 955 xlen=endx-orgx 956 dx=xlen/(numLDOSpt-1) 957 ylen=endy-orgy 958 dy=ylen/(numLDOSpt-1) 959 zlen=endz-orgz 960 dz=zlen/(numLDOSpt-1) 961 totlen=dsqrt(xlen**2+ylen**2+zlen**2) 962 dlen=dsqrt(dx**2+dy**2+dz**2) 963 LDOSptscomp=0D0 964 do ipt=1,numLDOSpt 965 x=orgx+dx*(ipt-1) 966 y=orgy+dy*(ipt-1) 967 z=orgz+dz*(ipt-1) 968 do imo=istart,iend 969 LDOSptscomp(imo,ipt)=fmo(x,y,z,imo)**2 970 end do 971 end do 972 973 enestep=(enehigh-enelow)/(num2Dpoints-1) 974 do i=1,num2Dpoints 975 LDOSxpos(i)=enelow+(i-1)*enestep 976 end do 977 978 LDOS2Dmap=0D0 979 if (ibroadfunc==1.or.ibroadfunc==3) then !Lorentzian function, see http://mathworld.wolfram.com/LorentzianFunction.html 980 do ilinept=1,numLDOSpt !Cycle each point in the line 981 do imo=1,imoend !Cycle each orbital 982 irealmo=imo 983 if (ispin==2) irealmo=imo+nbasis 984 preterm=str(imo)*0.5D0/pi*FWHM(imo) 985 do imappt=1,num2Dpoints !Broaden imo as curve 986 tmp=preterm/( (LDOSxpos(imappt)-MOene_dos(irealmo))**2+0.25D0*FWHM(imo)**2 ) 987 LDOS2Dmap(imappt,ilinept)=LDOS2Dmap(imappt,ilinept)+tmp*LDOSptscomp(imo,ilinept) 988 end do 989 end do 990 end do 991 end if 992 if (ibroadfunc==2.or.ibroadfunc==3) then 993 if (ibroadfunc==3) LDOS2Dmap=(1-gauweigh)*LDOS2Dmap 994 do ilinept=1,numLDOSpt !Cycle each point in the line 995 do imo=1,imoend !Cycle each orbital 996 irealmo=imo 997 if (ispin==2) irealmo=imo+nbasis 998 gauss_c=FWHM(imo)/2D0/sqrt(2*dlog(2D0)) 999 gauss_a=str(imo)/(gauss_c*sqrt(2D0*pi)) 1000 do imappt=1,num2Dpoints !Broaden imo as curve 1001 !Gaussian function, see http://en.wikipedia.org/wiki/Gaussian_function 1002 tmp=gauss_a*dexp( -(LDOSxpos(imappt)-MOene_dos(irealmo))**2/(2*gauss_c**2) ) 1003 if (ibroadfunc==3) tmp=gauweigh*tmp !Combine Lorentizan and Gaussian function 1004 LDOS2Dmap(imappt,ilinept)=LDOS2Dmap(imappt,ilinept)+tmp*LDOSptscomp(imo,ilinept) 1005 end do 1006 end do 1007 end do 1008 end if 1009 LDOS2Dmap=LDOS2Dmap*scalecurve 1010 write(*,*) 1011 write(*,"(' Energy range: from',f12.5,' to',f12.5,1x,a)") enelow,enehigh,trim(unitstr) 1012 write(*,"(' Starting point:',3f12.6,' Bohr')") orgx,orgy,orgz 1013 write(*,"(' End point: ',3f12.6,' Bohr')") endx,endy,endz 1014 write(*,"(' Stepsize:',f12.6,' Bohr, total length:',f12.6,' Bohr')") dlen,totlen 1015 1016 idraw=1 1017 isavepic=0 1018 do while(.true.) 1019 if (isilent==0.and.idraw==1) then 1020 lengthx=2300 1021 if (isavepic==0) then 1022 else if (isavepic==1) then 1023 end if 1024 end if 1025 1026 idraw=0 1027 isavepic=0 1028 write(*,*) 1029 write(*,*) " ======== Post-process menu ========" 1030 write(*,*) "0 Return" 1031 write(*,*) "1 Show graph again" 1032 write(*,*) "2 Save the picture to current folder" 1033 write(*,*) "3 Export curve data to plain text file in current folder" 1034 write(*,"(a,f8.3)") " 4 Set ratio of Y-axis relative to X-axis, current:",yxratio 1035 read(*,*) isel2 1036 1037 if (isel2==0) then 1038 exit 1039 else if (isel2==1) then 1040 idraw=1 1041 isavepic=0 1042 else if (isel2==2) then 1043 idraw=1 1044 isavepic=1 1045 else if (isel2==3) then 1046 open(10,file="LDOS.txt",status="replace") 1047 do imappt=1,num2Dpoints 1048 do ipt=1,numLDOSpt 1049 write(10,"(f8.3,f10.4,1PE16.8)") LDOSxpos(imappt),dlen*(ipt-1),LDOS2Dmap(imappt,ipt) 1050 end do 1051 end do 1052 close(10) 1053 write(*,*) "LDOS data have been written to LDOS.txt in current folder" 1054 write(*,*) "Column 1: Energy ("//trim(unitstr)//")" 1055 write(*,*) "Column 2: Coordinate in the line (Bohr)" 1056 write(*,*) "Column 3: LDOS value" 1057 else if (isel2==4) then 1058 write(*,*) "Input the ratio, e.g. 1.4" 1059 read(*,*) yxratio 1060 end if 1061 end do 1062end if 1063 1064 1065 1066 1067end do !End of main loop 1068end subroutine 1069