POSCAR转换为XYZ格式的Fortran小代码

POSCAR转换为XYZ，且可以扩展原胞大小

program morePOSCAR
implicit none
real :: vector(3,3),lattice,x(3),y(3)
INTEGER::flagname,flagc2d,flagd2c
integer ::numin,species,i,j,k,temp,a,b
character(len=50) ::filenamein,filenameout  
character(len=50) ::name               
character(len=70) ::test,backup
character(len=50) ::coordinate
integer ::nx,ny,nz,numout,atomnow,supernow
real,allocatable ::input(:,:),output(:,:)
integer,allocatable ::numspecies(:),whichspecie(:)
character(2),allocatable::speciename(:)
integer ::fileid,outid

write(*,*) "文件名：1. 手工输入；2. 默认文件“POSCAR"
read(*,*) flagname
if(flagname==1) then
write(*,*) "请输入需要转换为XYZ的POSCAR的文件名称:"
read(*,*) filenamein
filenamein=trim(filenamein)
else
filenamein='POSCAR'
filenamein=trim(filenamein)
end if
write(*,*) "正在打开",filenamein
fileid=3208
outid=3301
open(unit=fileid, file=filenamein,status='old')
read(fileid,*) name
name=adjustl(name)
name=trim(name)
write(*,*) name
read(fileid,*) lattice
write(*,*)  lattice
read(fileid,*) vector(1,1),vector(1,2),vector(1,3)
read(fileid,*) vector(2,1),vector(2,2),vector(2,3)
read(fileid,*) vector(3,1),vector(3,2),vector(3,3)
write(*,*) vector(1,1),vector(1,2),vector(1,3)
write(*,*) vector(2,1),vector(2,2),vector(2,3)
write(*,*) vector(3,1),vector(3,2),vector(3,3)
read(fileid,"(a)") test
backup=test
species=0
27   j=len_trim(test)
     if(j.gt.0) then
     test=adjustl(test)
     i=index(test,' ')
  species=species+1
     test=test(i:j)
     goto 27
     endif
allocate(numspecies(species))
allocate(speciename(species))
allocate(whichspecie(species))
backspace(fileid)
read(fileid,*) (numspecies(i),i=1,species)
write(*,*) (numspecies(i),i=1,species)
numin=0
do i=1,species
numin=numin+numspecies(i)
if(i==1) then
whichspecie(i)=numspecies(i)
else
whichspecie(i)=whichspecie(i-1)+numspecies(i)
end if
end do
allocate(input(numin,3))
read(fileid,*) coordinate
coordinate=adjustl(coordinate)
write(*,*) coordinate
if(index(coordinate,'C')==1 .or. index(coordinate,'c')==1) then
write(*,*) "在POSCAR里发现了",species,"类不同原子，请给出元素名"
do i=1,species
write(*,*) i,"类原子名称："
read(*,*) speciename(i)
end do
do i=1,numin
read(fileid,*) (input(i,j),j=1,3)
end do
write(*,*) "在第一个基矢方向上扩展数目"
read(*,*) nx
write(*,*) "在第二个基矢方向上扩展数目"
read(*,*) ny
write(*,*) "在第三个基矢方向上扩展数目"
read(*,*) nz
numout=numin*nx*ny*nz
allocate(output(numout,3))
supernow=0
do i=1,nx
do j=1,ny
do k=1,nz
 do a=1,numin
 atomnow=numin*supernow+a
    output(atomnow,1)=input(a,1)+lattice*(vector(1,1)*(i-1)+vector(2,1)*(j-1)+vector(3,1)*(k-1))
 output(atomnow,2)=input(a,2)+lattice*(vector(1,2)*(i-1)+vector(2,2)*(j-1)+vector(3,2)*(k-1))
 output(atomnow,3)=input(a,3)+lattice*(vector(1,3)*(i-1)+vector(2,3)*(j-1)+vector(3,3)*(k-1))
 end do
    supernow=supernow+1
end do
end do
end do
filenameout=trim(filenamein)//".xyz"
write(*,*) '生成新文件',filenameout,'记录结果'
open(unit=outid, file=filenameout)
close(outid,status='DELETE')
open(unit=outid,file=filenameout)
write(outid,*) numout
write(outid,*) name
do a=1,numout
 b=mod(a,numin)
 if(b==0) then
 b=numin
 end if
 do i=1,species
  if(i==1) then
   if(b<=whichspecie(i)) then
   k=1
   goto 2007
   end if
  else
   if(b>whichspecie(i-1) .AND. b<=whichspecie(i)) then
   k=i
   goto 2007
   end if
  end if
 end do

2007 write(outid,*) speciename(k),output(a,1),output(a,2),output(a,3)
end do
close(outid)
else
write(*,*) "ERROR!"
write(*,*) "非笛卡儿坐标，无法转换成XYZ格式！"
end if
close(fileid)
end program