[求助] 建立一个通道

program main
implicit none

integer,parameter :: nmax=80000
!input
integer :: ncx,ncy,ncz
real(8) :: a,a1,a2,a3
!output
integer :: n
real(8) :: rx(nmax),ry(nmax),rz(nmax)
!local
integer :: i

call files_read(ncx,ncy,ncz,a)
n=8*ncx*ncy*ncz
a1=dble(ncx)*a
a2=dble(ncy)*a
a3=dble(ncz)*a

call set_silicon(nmax,ncx,ncy,ncz,a,rx,ry,rz)
call writepos(nmax,n,a1,a2,a3,rx,ry,rz)

stop
end program main


!读取文件
subroutine files_read(ncx,ncy,ncz,a)
implicit none
!global
integer :: ncx,ncy,ncz
real(8) :: a
! locala1
integer :: i
open(unit=10,file='input.in',status='old',action='read',position='rewind')
read(10,*) ncx
read(10,*) ncy
read(10,*) ncz
read(10,*) a

close(10)
return
end subroutine files_read


!写入文件
subroutine writepos(nmax,n,a1,a2,a3,rx,ry,rz)
implicit none
! global
integer,intent(in) :: nmax,n
real(8),intent(in) :: a1,a2,a3,rx(nmax),ry(nmax),rz(nmax)
! local
integer :: i,j
open(unit=20,file='supercell.dat',status='replace',action='write')
201 format(I5,A5,A5,I5, F8.3, F8.3, F8.3)
open(unit=21,file='guding.dat')
j=0

do i=1,n
if(rx(i)*rx(i)/100+ry(i)*ry(i)/100>=2.25d0) then
j=j+1
write(20,201) 1,'SIL','SI',j,rx(i)/10,ry(i)/10,rz(i)/10
end if
end do

close(20)

return
end subroutine writepos


subroutine set_silicon(nmax,ncx,ncy,ncz,a,rx,ry,rz)
implicit none
! global
integer,intent(in) :: nmax,ncx,ncy,ncz
real(8),intent(in) :: a
real(8),intent(out) :: rx(nmax),ry(nmax),rz(nmax)
! locala1
real(8) :: unitlen,halflenx,halfleny,halflenz
integer :: counter,xcount,ycount,zcount
integer :: nx,ny,nz
integer :: i,j,k,judge

unitlen=a/4.0d0

nx=4*ncx
ny=4*ncy
nz=4*ncz

xcount=0
ycount=0
zcount=0
counter=0

do i=1,nx
        do j=1,nz
                do k=1,ny
                        judge=0
                        select case(xcount)
                                case(0)
                                        if((ycount==0.and.zcount==0).or.(ycount==2.and.zcount==2)) judge=1
                                case(1)
                                        if((ycount==1.and.zcount==1).or.(ycount==3.and.zcount==3)) judge=1
                                case(2)
                                        if((ycount==2.and.zcount==0).or.(ycount==0.and.zcount==2)) judge=1
                                case(3)
                                        if((ycount==1.and.zcount==3).or.(ycount==3.and.zcount==1)) judge=1
                                case default
                                        stop 'setdiamond error !'
                        end select
                        if(judge==1) then
                                counter=counter+1
                                rx(counter)=dble(i-1)*unitlen-dble(ncx)*a/2
                                ry(counter)=dble(k-1)*unitlen-dble(ncy)*a/2
                                rz(counter)=dble(j-1)*unitlen-dble(ncz)*a/2
                        end if
                        ycount=ycount+1
                        if(ycount>3.or.k==ny) ycount=0
                end do
                zcount=zcount+1
                if(zcount>3) zcount=0
        end do
        xcount=xcount+1
        if(xcount>3) xcount=0
end do

if(counter>nmax) then
    stop 'nmax is too small !'
end if

return
end subroutine set_silicon