C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%   Subroutine "readpos" is used to read the undeformed    %%%%
C%%%%   cell vector and atomic coordinates from the structure  %%%%
C%%%%   file with format "VASP"                                %%%%
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        subroutine readpos(inpos,
     &                     title, 
     &                     alat, 
     &                     privect, 
     &                     sname,
     &                     natomi,               
     &                     ntype,
     &                     nn,  
     &                     Ndyn,
     &                     Ncoord,
     &                     pripos,
     &                     fix,
     &                     atom_types)
        
        double precision alat,privect,pripos
        dimension privect(3,3),pripos(1000,3)
        character*10 title,fix
        character*1 pricoord,seldyn
        character*5 sname
        character*10 inpos        
        dimension sname(10),fix(1000,3)
        integer i,j,k,ntype,natomi,nn
        integer Ncoord,Ndyn,atom_types,anum
        dimension natomi(10),atom_types(1000),anum(11)

        logical filefound

C       To determine whether the file exists        

        inquire (file=inpos, exist=filefound)
        if (.NOT.filefound) then
                open(unit=19,position='Append',FILE='RELAS') 
                write(19,'(a17, a, a17)') 'ERROR!! The file ',
     &                        trim(inpos), ' does not exit!!!'
                close(19)
                stop        
        endif

        open(20,file=Trim(inpos))
        
        call nnumber(7,20,ntype)
        close(20)

        open(20,file=Trim(inpos))
        read(20,*) title
        read(20,*) alat
        do i=1,3
           read(20,*) (privect(i,j),j=1,3)
        enddo
        read(20,*) (sname(i),i=1,ntype)
        read(20,*) (natomi(i),i=1,ntype)

        nn=0
        do i =1, ntype
           nn=nn+natomi(i) 
        enddo
        read(20,*) seldyn

C       Ndyn=1: Selective Dynamics switched on
C       Ndyn=0: Selective Dynamics switched off

        if (seldyn.EQ.'s'.OR.seldyn.EQ.'S') THEN
                Ndyn=1
                read(20,*) pricoord
                do i=1, nn
                read(20,*) (pripos(i,j),j=1,3),(fix(i,j),j=1,3) 
                enddo
        else
                Ndyn=0
                pricoord=seldyn
                do i=1,nn
                   read(20,*) pripos(i,1), pripos(i,2), pripos(i,3)
                   do j=1,3
                        fix(i,j)='T'
                   enddo
                enddo
        endif

        close(20)

C       Ncoord=1: Cartesian mode
C       Ncoord=0: Dirte mode

        if (pricoord.EQ.'c'.OR.pricoord.EQ.'C') THEN
                       Ncoord=1
        else
                if (pricoord.EQ.'d'.OR.pricoord.EQ.'D') THEN
                       Ncoord=0
                else
                       open(unit=19,position='Append',FILE='RELAS') 
                       write(19,*) 'ERROR!! NOT DIRECT OR CARTESIAN!!!'
                       close(19)
                       stop
                endif
        endif

C       to define the atom_type

        anum(1)=0
        do i=1,ntype
           anum(i+1)=0
           do j=1,i
              anum(i+1)=anum(i+1)+natomi(j)
           enddo
        enddo

        do i=1, ntype
           do j=anum(i)+1, anum(i+1)
               atom_types(j)=i           
           enddo
        enddo

C       end of defining the atom_type 

        end
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%%                    End of readpos.f                     %%%%
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%        
