问题描述
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
fotran注释行编写了使用方法。另外，部分元素没有加入，需要时可以手动加入。
