program main

!-------------------------------------------------------------------------------------------

! This program can be used to convert myfile.cell file (obtained from Material

! Studio software by "runing" CASTEP modules ) to myfile.fdf

! To run this program ,you have to compile and run it in current folder which

! contains myfile.cell file iiii.

! $ gfortran -o cell2fdf cell2fdf.f90   #maybe you can use other compilers ,such as f95 ..

! $ ./cell2fdf

! $ myfile          # "myfile" is the prefix of your cell file ,input it by the keyboard !!! very important

! Then you will find myfile.fdf file

! In order to run siesta program ,you just have to cp *.psf file to your work folder and set

! calculation parameters in para.fdf

! If U find debugs in this program , you may contact me by haidi@mail.ustc.edu.cn

! Thanks for your use

!-----------------------------------------------------------------------------------------------

implicit none

integer,parameter::N=20  !the max Species N.O.

character(len=90)::str,vector(3)

character(len=30)::fdf,filename,suffix,prefile

character(len=2),allocatable::ElemSymb(:)

character(len=90),allocatable::frac(:)

integer::fileid,i,j,aux,total,SpeciNo(N),k,m,point,debug

integer,allocatable::ElemZ(:)

integer,external:: GetElemZ

character,external::num2char

!=================Debug=============

!debug=1 for debug

debug=0

!===================================

read(5,err=20,end=20,fmt='(a)') str

20     continue

filename=trim(adjustl(str))//'.cell'

fileid=10

open(unit=fileid,file=filename)

!=================debug===============

if(debug==1) print*,'Input string: '//trim(str)

!=====================================

i=0

do while(.true.)

   read(unit=fileid,fmt='(a)')str

   i=i+1

   if(str(1:24)=="%ENDBLOCK POSITIONS_FRAC")then

   close(fileid)

   goto 30

   end if

end do

30 continue

total=i

aux=total-8

allocate(frac(aux))

allocate(ElemSymb(aux))

allocate(ElemZ(aux))

open(unit=fileid,file=filename)

do i=1,total-1

     read(unit=fileid,fmt='(a)')str

     if(i>1 .and. i<5) then

        vector(i-1)=str

     elseif(i>=8)then

        frac(i-7)=str

     end if

end do

close(unit=fileid)

i=index(filename,'.cell')

prefile=filename(1:(i-1))

fdf=trim(prefile)//".fdf"

 do i=1,aux

   str=frac(i)

   ElemSymb(i)=str(2:3)

   ElemZ(i)=GetElemZ(ElemSymb(i))

   frac(i)=str(4:66)

 end do

!========================debug=========================

if(debug==1)then

 print*,"total N.O.of Atom",aux

 print*,'ElemSymb',ElemSymb

 print*,'ElemZ',ElemZ

 print*,'frac',frac

end if

!======================================================

!core for calculate NO of species and SpeciNo

m=1

k=0

do i=1,aux

  if(i==aux)then

     k=k+1

     SpeciNo(m)=k

  else

     j=i+1

       if(ElemZ(i)==ElemZ(j))then

           k=k+1

       else

           SpeciNo(m)=k+1

           k=0

           m=m+1

       endif

  endif

enddo

!========================Debug==========================

if(debug==1)then

print*,'N.O. of species:  '//num2char(m)

print*,'SpeciNo:',speciNo

endif

!=======================================================

call writetofdf(frac,SpeciNo,m,ElemZ,ElemSymb,aux,prefile,fdf,vector)

end

subroutine writetofdf(frac,SpeciNo,m,ElemZ,ElemSymb,aux,filename,fdf,vector)

implicit none

integer aux

character(len=90) frac(aux),vector(3)

character(len=2) ElemSymb(aux)

integer SpeciNo(m),ElemZ(aux),m,i,count,j

character(len=30) fdf,filename

character(len=90)str

integer,external::GetElemZ

character(len=10),external ::num2char

open(unit=10,file=fdf)

write(10,'(a)')'SystemName    '//filename

write(10,'(a)')'SystemLabel   '//filename

write(10,'(a)')'NumberOfSpecies'//'  '//num2char(m)

write(10,'(a)')'NumberOfAtoms'//'  '//num2char(aux)

write(10,'(a)')'%block ChemicalSpeciesLabel'

count=0

if(aux==1)then

 write(10,'(30a)')trim(num2char(1)//'  '//num2char(GetElemZ(ElemSymb(1)))//' '//ElemSymb(1))

 print*,ElemSymb(1)

else

do i=1,m

  count=count+SpeciNo(i)

    str=num2char(i)//'  '//num2char(GetElemz(ElemSymb(count)))//'  '//ElemSymb(count)

   write(10,'(a)')trim(str)

   print*,ElemSymb(count)

 enddo

end if

write(10,'(a)')'%endblock ChemicalSpeciesLabel'

write(10,'(a)')'LatticeConstant    1.0  Ang'

write(10,'(a)')'%block LatticeVectors'

do i=1,3

write(10,'(a)')vector(i)

end do

write(10,'(a)')'%endblock LatticeVectors'

write(10,'(a)')""

write(10,'(a)')'AtomicCoordinatesFormat     Fractional'

write(10,'(a)')'%block AtomicCoordinatesAndAtomicSpecies'

count=1

if(aux==1)then

  write(10,'(a)')trim( frac(1))//"    "//num2char(1)

else

do i=1,m

!  count=count+SpeciNo(i)

   do j=count,SpeciNo(i)+count-1

      write(10,'(a)')trim( frac(j))//"    "//num2char(i)

   end do

  count=count+SpeciNo(i)

end do

end if

write(10,'(a)')'%endblock AtomicCoordinatesAndAtomicSpecies'

write(10,'(a)')'%include para.fdf'

!call copyfile('later.fdf',fdf)

close(10)

end subroutine

subroutine copyfile(file1,file2)

implicit none

character(len=90) line

character(len=30)file1,file2

open(11,file=trim(file1),status='old')

do while(.true.)

 read(11,err=40,end=40,fmt='(a)')line

 write(10,'(a)')line

end do

40     continue

return

end subroutine

character(len=*) function num2char(num)

implicit none

integer num

character(len=10) str

write(str,'(i4)')num

num2char=trim(adjustl(str))

return

end function

integer function GetElemZ(str)

implicit none

integer i

character(len=2) str

if(str=="H" .or. str=="H " .or. str== " H" )then

   i=1

elseif(str=="He")then

   i=2

elseif(str=="Li")then

   i=3

elseif(str=="Be")then

   i=4

elseif(str=="B" .or. str=="B " .or. str== " B" )then

   i=5

elseif(str=="C" .or. str=="C " .or. str== " C" )then

   i=6

elseif(str=="N" .or. str=="N " .or. str== " N" )then

   i=7

elseif(str=="O" .or. str=="O " .or. str== " O" )then

   i=8

elseif(str=="F" .or. str=="F " .or. str== " F" )then

   i=9

elseif(str=="Na")then

   i=11

elseif(str=="Mg")then

   i=12

elseif(str=="Al")then

   i=13

elseif(str=="Si" )then

   i=14

elseif(str=="P" .or. str=="P " .or. str== " P" )then

   i=15

elseif(str=="S" .or. str=="S " .or. str== " S" )then

   i=16

elseif(str=="Cl")then

   i=17

elseif(str=="K" .or. str=="K " .or. str== " K" )then

   i=19

elseif(str=="Ca")then

   i=20

elseif(str=="Sc")then

   i=21

elseif(str=="Ti")then

   i=22

elseif(str=="V" .or. str=="V " .or. str== " V" )then

   i=23

elseif(str=="Cr" )then

   i=24

elseif(str=="Mn")then

   i=25

elseif(str=="Fe")then

   i=26

elseif(str=="Co")then

   i=27

elseif(str=="Ni")then

   i=28

elseif(str=="Cu")then

   i=29

elseif(str=="Zn")then

   i=30

elseif(str=="Ga")then

   i=31

elseif(str=="Ge")then

   i=32

elseif(str=="As")then

   i=33

elseif(str=="Se")then

   i=34

elseif(str=="Br")then

   i=35

elseif(str=="Rb")then

   i=37

elseif(str=="Sr")then

   i=38

elseif(str=="Y" .or. str=="Y " .or. str== " Y" )then

   i=39

elseif(str=="Zr")then

   i=40

elseif(str=="Nb")then

   i=41

elseif(str=="Mo")then

   i=42

elseif(str=="Tc")then

   i=43

elseif(str=="Ru")then

   i=44

elseif(str=="Rh")then

   i=45

elseif(str=="Pd")then

   i=46

elseif(str=="Ag")then

   i=47

elseif(str=="Cd")then

   i=48

elseif(str=="In")then

   i=49

elseif(str=="Sn")then

   i=50

elseif(str=="Sb")then

   i=51

elseif(str=="Te")then

   i=52

elseif(str=="I" .or. str==' I' .or. str=='I ')then

   i=53

elseif(str=="Cs")then

   i=55

elseif(str=="Ba")then

   i=56

elseif(str=="Bi")then

   i=83

else

   i=0

   stop "element not found !"

endif

GetElemZ=i

 if(GetElemZ==0) stop

return

end function