Tek-Tips is the largest IT community on the Internet today!

Members share and learn making Tek-Tips Forums the best source of peer-reviewed technical information on the Internet!

  • Congratulations strongm on being selected by the Tek-Tips community for having the most helpful posts in the forums last week. Way to Go!

NaN output outside but not inside DO loop

Status
Not open for further replies.

Vahid9

Technical User
Jul 28, 2016
23
CA
Dear Fortran Users,

In the following code, I can print all the values for conc_h inside the first DO loop over i before ENDDO and the total conc_h right after the ENDDO.
When I try to print conc_e inside the second DO loop over i, I get REAL values but when I put the print statement right after the ENDDO,
I get NaN. Th two DO loops over i seem identical. I have attached the input file just in case. I may be missing something obvious but I cannot find it. Any hints would be appreciated.

Thank you,
Vahid

PROGRAM conc_fermi
IMPLICIT NONE

INTEGER, PARAMETER :: max_rows = 200000
INTEGER :: row, rowc_min, rowc_max, rowv_min, rowv_max
INTEGER :: i,j,io,total_row
REAL, DIMENSION(max_rows) :: E, DOS, DOSTOT
REAL :: EVmin, EVmax, ECmin, ECmax, E_F, conc, conc_e, conc_h
REAL :: DOS_ECmin, DOS_ECmax, DOS_EVmin, DOS_EVmax, V, dE, kT
CHARACTER (20) :: input_file
REAL, PARAMETER :: temp=300 ! K
REAL, PARAMETER :: kb = 1.3806488e-23 ! J/K
REAL, PARAMETER :: qq = 1.60217657e-19 ! eV
REAL, PARAMETER :: bohr_2_cm = 5.2917721067e-9

kT = kb*temp/qq ! eV
EVmax=6.28
EVmin=-5.83
ECmax=21.47
ECmin=6.8
V=262.8707

PRINT*, "Enter DOS file name"
READ*, input_file

! Convert V from bohr^3 to cm^3
V=V*(bohr_2_cm)**3

OPEN(UNIT=100, FILE=input_file, STATUS='old', ACTION='read')

! Read the DOS file
E:))=0
DOS:))=0
DOSTOT:))=0
DO row=1,max_rows
READ(100,*,IOSTAT=io) E(row), DOS(row), DOSTOT(row)

IF (io .LT. 0) THEN
total_row = row - 1
EXIT
ENDIF

! Define row number and DOS at EVmax
IF ( E(row) .EQ. EVmax ) THEN
rowv_max = row
DOS_EVmax = DOS(rowv_max)
ENDIF

! Define row number and DOS at EVmin
IF ( E(row) .EQ. EVmin ) THEN
rowv_min = row
DOS_EVmin = DOS(rowv_min)
ENDIF

! Define row number and DOS at ECmax
IF ( E(row) .EQ. ECmax ) THEN
rowc_max = row
DOS_ECmax = DOS(rowc_max)
ENDIF

! Define row number and DOS at ECmin
IF ( E(row) .EQ. ECmin ) THEN
rowc_min = row
DOS_ECmin = DOS(rowc_min)
ENDIF

ENDDO
CLOSE(100)

! energy interval in DOS
dE = ABS(E(1)-E(2))

OPEN(1,file='concentration-Fermi.dat')
WRITE(1,*) "Fermi Level (eV) Electron concentration (in 1/cm^3)"

! Vary E_F from 0.3eV below to 0.3eV above ECmin

DO j = -3, 3

E_F=ECmin+(REAL(j)/10)

! Trapezoid Rule
conc_h = 0.0
DO i = rowv_min, rowv_max
conc_h = conc_h + DOS(i) * (1-(1/(1+exp((E(i)-E_F)/kT))))
END DO
PRINT*, conc_h
conc_h=dE*(conc_h-0.5*(DOS_EVmin*(1-(1/(1+exp((EVmin-E_F)/kT)))) &
+ DOS_EVmax*(1-(1/(1+exp((EVmax-E_F)/kT))))))

conc_e = 0.0
DO i = rowc_min, rowc_max
conc_e = conc_e + DOS(i) * (1/(1+exp((E(i)-E_F)/kT)))
! PRINT*, conc_e ! prints all the conc_e
END DO
PRINT*, conc_e ! prints NaN
conc_e = dE*(conc_e-0.5*(DOS_ECmin*(1/(1+exp((ECmin-E_F)/kT))) &
+ DOS_ECmax*(1/(1+exp((ECmax-E_F)/kT)))))
conc=conc_h-conc_e

WRITE(1,*) E_F, " ", conc/V
ENDDO !j

CLOSE(1)
END
 
 http://files.engineering.com/getfile.aspx?folder=e648489e-2948-4c0d-aa81-bde8300ef12f&file=si.dos
I modified the print statements and tried your program
Code:
...
     DO i = rowv_min, rowv_max
       conc_h = conc_h + DOS(i) * (1-(1/(1+exp((E(i)-E_F)/kT)))) 
     END DO
     [highlight]PRINT*, 'conc_h =', conc_h[/highlight]
...
     DO i = rowc_min, rowc_max
       conc_e = conc_e + DOS(i) * (1/(1+exp((E(i)-E_F)/kT)))
       ! PRINT*, conc_e ! prints all the conc_e
     END DO
     [highlight]PRINT*, 'conc_e =', conc_e ! prints NaN[/highlight]
...

and got normal output without an error:

Code:
$ gfortran vahid9.f95 -o vahid9

$ vahid9
 Enter DOS file name
si.dos
 conc_h =   1.76246031E-05
 conc_e =   6.62863044E-07
 conc_h =   3.68310367E-07
 conc_e =   3.17181075E-05
 conc_h =   7.69643638E-09
 conc_e =   1.50968332E-03
 conc_h =   1.60829211E-10
 conc_e =   5.95530570E-02
 conc_h =   3.36078019E-12
 conc_e =  0.701390862
 conc_h =   7.02287620E-14
 conc_e =   2.32466984
 conc_h =   1.46748299E-15
 conc_e =   4.66214800
 
Thanks Mikrom. The NaN seem to be related to compiling with ifort which is the compiler I used. When I compile with gfortran, I get the same result as yours.

Best wishes,
Vahid
 
Status
Not open for further replies.

Part and Inventory Search

Sponsor

Back
Top