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