Materials cell文件转换为siesta fdf格式

materials

发布日期: 2021-07-01

文章字数: 1.4k

阅读时长: 8 分

阅读次数:

问题描述

Materials Studio 软件支持的结构格式为.cell，如何将cell转换为 siesta软件的fdf输入文件？

程序

cell2fdf.f90

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 !!!.
! $ gfortran -o cell2sieta cell2fdf.f90   #maybe you can use other compilers ,such as f95 ,ifort ...
! $ ./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 
! 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