Error: Unclassifiable statement in if conditional

Question

This is the program I wrote to calculate the distance between carbon and hydrogen atom coordinates as a practice. But some errors occurred and I don't know how to fix them.

    program C_H_bdlength
    implicit none
    integer                 :: ia            ! integer atoms number               
    integer                 :: na            ! number of atoms                    
    real*8                  :: x  ,  y , z   ! x y z coordinates of all atoms     
    real*8                  :: x1 , y1 , z1  ! x y z coordinates of C atoms       
    real*8                  :: x2 , y2 , z2  ! x y z coordinates of H atoms       
    character(len=2)        :: ele           ! element name                       
    real*8                  :: dis           ! distance between C and H           

    open(100,file='cnt.ini',action='read',status='old')
    na = 0
    do
     read(100,*,end=101)
     na = na + 1
    end do
    101 continue

    open(200, file='cnt.ini',action='wtire',status='replace')
    rewind(100)

    write(200,*)

    do ia = 1 , na
     read(100,*) ele, x, y, z
       if ele == 'C' then x1 = x, y1 = y, z1 = z
          do ia = ia , na
             read(100,*) ele, x, y, z
               if ele == 'H' then x2 = x, y2 = y, z2 = z
               dis = sqrt((x1 - x2)**2 + (y1 - y2)**2 + (z1 - z2)**2)

               write(200,*) dis
               end if
          end do
       end if
    end do
    
    end program C_H_bdlength

Errors:

  In file C_H_bdlength.f90:26

       if ele == 'C' then x1 = x, y1 = y, z1 = z
      1
    Error: Unclassifiable statement at (1)
 In file C_H_bdlength.f90:27

          do ia = ia , na
                        1
 In file C_H_bdlength.f90:24

    do ia = 1 , na
               2
Error: Variable 'ia' at (1) cannot be redefined inside loop beginning at (2)
 In file C_H_bdlength.f90:29

               if ele == 'H' then x2 = x, y2 = y, z2 = z
              1
Error: Unclassifiable statement at (1)
 In file C_H_bdlength.f90:33

               end if
                 1
Error: Expecting END DO statement at (1)
 In file C_H_bdlength.f90:35

       end if
         1
Error: Expecting END DO statement at (1)

Answer 1

You have a nested loop, and are using the same variable (ia) to control both the inner and outer loop. Use a variable with a different name in the inner loop.

The if <condition> then construct cannot have additional operations following the then on the same line.

Answer 2

In addition to @Peter's answer, another (not syntactic) problem is that your DO loops do not take into account all the C-H pairs. For example, if "cnt.ini" contains the following data

O  1.0  1.0  1.0
H  8.0  1.0  6.0  # ia=2
C  2.0  3.0  1.0  # ia=3
N  3.0  3.0  3.0
H  4.0  1.0  6.0  # ia=5
O  5.0  5.0  5.0
C  1.0  7.0  8.0  # ia=7

it only takes into account ia = 3 and 5 pair. Further, because you are writing output data onto the input file (cnt.ini) while reading the input from the latter, data transfer may be broken. So I suggest first reading in the input data into arrays, calculate C-H distances, and then write the result onto a different file. For example, the code may look like this.

integer :: ia, ja, na
real*8  :: dist
real*8, allocatable :: pos(:,:)
character(len=2), allocatable :: ele(:)

open(100,file='cnt.ini',status='old')

!! First, determine the number of atoms in the file.
na = 0
do
    read( 100, *, end=101 )
    na = na + 1
end do
101 continue

!! Then, allocate necessary arrays for atomic positions and element names.
allocate( pos( 3, na ), ele( na ) )

!! Now read the actual data in the file.
rewind( 100 )
do ia = 1, na
    read( 100, * ) ele( ia ), pos( 1:3, ia )

    ! Alternative way using implied DO-loop.
    ! read( 100, * ) ele( ia ), ( pos( k, ia ), k = 1, 3 )
enddo
close( 100 )

!! Now calculate C-H distances.
open(200, file='dist.dat')

!! Loop over unique atom pairs.
do ia = 1 ,     na - 1
do ja = ia + 1, na

    if ( ( ele(ia) == 'C' .and. ele(ja) == 'H' ) .or. &
         ( ele(ia) == 'H' .and. ele(ja) == 'C' ) ) then

        dist = norm2( pos( :, ia ) - pos( :, ja ) )
        ! dist = sqrt( sum( (pos(:,ia) - pos(:,ja))**2 ) )   !! if norm2() not available

        write(200,*) ia, ja, dist
    endif

enddo
enddo