       Program zscreatgb

c     This program will create a [001] twist grain boundary in initial input 
c     lattice.

      implicit real*8 (a-h,o-z)
      parameter(natmax=100000)  
      dimension itype(natmax),xcrd(natmax),ycrd(natmax),zcrd(natmax)
      common /all/theta,natom,xcrd,ycrd,zcrd,xsize,ysize,zsize,itype,
     1	nnatom,xolds,yolds,conlat,wxsize,wysize,wzsize
 
      open (9, file='lattice.crd')
      open (10, file='clmps.crd')

      print *,'input rotation angle'
      read (*,*) theta
      print *,'input lattice constant'
      read (*,*) conlat
      print *, ' Total number of atoms: '
      read (*,*) natom
      read(9,*) xsize,ysize,zsize
      read(9,3000)(itype(i),xcrd(i),ycrd(i),zcrd(i),i=1,natom)

      call transform

      write(*,*) 'total number of atoms: ',nnatom
      write(10,*) wxsize,wysize,wzsize
      do 2000 j=1,2*natom
	if(itype(j).eq.1)then
           write(10,3000) itype(j),xcrd(j),ycrd(j),zcrd(j)
	end if
 2000 continue

      do 2200 j=1,2*natom
	if(itype(j).eq.2)then
           write(10,3000) itype(j),xcrd(j),ycrd(j),zcrd(j)
	end if
 2200 continue

 
 3000 format(i3,3x,3e25.16)
      end
C-------------------------------------------------------------------------- 
      subroutine transform

      implicit real*8 (a-h,o-z)
      parameter (natmax=100000)
      dimension itype(natmax),xcrd(natmax),ycrd(natmax),zcrd(natmax)
      common /all/theta,natom,xcrd,ycrd,zcrd,xsize,ysize,zsize,itype,
     1	nnatom,xolds,yolds,conlat,wxsize,wysize,wzsize

c     Duplicate the lattice along z-axis
      do i=1,natom
         zcrd(i)=zcrd(i)-zsize/2.0D0-conlat/4.0D0
      enddo
C     this gets a-a at the interface and b-b at the other
      do 3010 i=1,natom
         j=natom+i
         itype(j)=itype(i)
         xcrd(j)=xcrd(i)
         ycrd(j)=ycrd(i)
         zcrd(j)=-zcrd(i)
 3010    continue

c     rotate half of the lattice around z-axis
      do 4017 j=natom+1,2*natom
            xcrd(j) = xcrd(j)*cos(theta) + ycrd(j)*sin(theta)
            ycrd(j) = -xcrd(j)*sin(theta) + ycrd(j)*cos(theta)
 4017   continue

c     shift x,y,z-axis
      do j=1,2*natom
        xcrd(j)=xcrd(j)-0.25*conlat
	ycrd(j)=ycrd(j)-0.25*conlat 
      enddo
	
c     set the new box size after rotation
      wxsize=conlat*6.0
      wysize=conlat*6.0
      wzsize=2.0*zsize

c     cut the surface layers of atoms
      nnatom=2*natom
      do 4019 j=1,2*natom   
        if (abs(xcrd(j)).gt.(0.5*wxsize))then
           itype(j)=0
           nnatom=nnatom-1
        else if (abs(ycrd(j)).gt.(0.5*wysize))then
           itype(j)=0
           nnatom=nnatom-1
	end if
 4019   continue


      return
      end  

   
