Program Graphite_Generation
* *** **************************************************************
*                Generates surface and slit pores of graphite 
*             You need: 1.Nx, Ny - number of unit cells to generate
*                       2.H - pore width
*                       3.edge - size of the edge if adsorption in
*                         finite (lateral) size patch will be simulated
*
*            				Version January 30, 2011
* *** **************************************************************
*                     PARAMETERS

      parameter(a=1.42, b=2.46)
*          a = C-C bond
*          b = a*sqrt(3)
      parameter(Nx=8,Ny=13)
*          Nx,Ny = number of unit cells in (x,y) direction
      parameter(H=7.)
*          H = pore width
      parameter(edge=15.)
*          edge = distance from the patch edge and box edge
	parameter(nSubstrate_Gr=1)
	parameter(nSubstrate_Pore=1)
	parameter(nSubstrate_PoreCentered=1)
*          flags for surface, pore or finite pore generation

* *** **************************************************************
*                    VARIABLES
	integer atype
*          atype = 1 	--> atom is C 
*          atype = 2 	--> atom is C alpha
*          atype = 3  --> atom is C betha
*          atype = 4  --> atom is B
*          atype = 9  --> atom is H
      real x1,x2,x3,x4, y1,y2,y3,y4,x5,y5,z,zH

* *** **************************************************************
*                    OPEN FILES

      open(1, file ='Graphite_8x13.dat', status = 'unknown')
	open(2, file ='Pore_8x13_H7.dat', status = 'unknown') 
	open(3, file ='Pore_centered_8x13_H7.dat', status = 'unknown')  

* *** **************************************************************
	if(nSubstrate_Gr.eq.1)then
* *** **************************************************************
*                   GRAPHITE SURFACE GENERATION
*
*      (0,0) - center of hexagon
*                      |     C----C	  |
*          Building    |    /      \    |
*          block:      |---C	      C---|	a(graphite)= 3*a 
*                      |    \      /    |	b(graphite)= 2*b 
*                      |     C----C	  |
*                                

      atype = 1                  !all atoms are C
	z = 0.0                    !graphite is flat!!!

      do j=1,Ny
	   do i=1,Nx

     	      x1 = a + (i-1)*3*a	     !first line of atoms
	      y1 = (j-1)*b
	      x2 = 2*a +(i-1)*3*a
            y2 = y1

	      x3 = a/2 + (i-1)*3*a	 !second line of atoms
	      y3 = b/2 + (j-1)*b
	      x4 = 2.5*a +(i-1)*3*a
            y4 = y3

	      write(1,1000) x1,y1,z,atype
	      write(1,1000) x2,y2,z,atype
	      write(1,1000) x3,y3,z,atype
	      write(1,1000) x4,y4,z,atype
 1000	format(2x,3(2x,f12.6),2x,i3)

	   enddo
     	enddo
 	endif

* *** ***********************************************************************
      if(nSubstrate_Pore.eq.1)then
* *** ***********************************************************************
*                   INFINITE SLIT PORE GENERATION (PORE WIDTH = H)                                

      atype = 1                  !all atoms are C
	z = 0.0                    !graphite is flat!!!
	zH = z + H

      do j=1,Ny
	   do i=1,Nx

     	      x1 = a + (i-1)*3*a	     !first line of atoms
	      y1 = (j-1)*b
	      x2 = 2*a +(i-1)*3*a
            y2 = y1

	      x3 = a/2 + (i-1)*3*a	 !second line of atoms
	      y3 = b/2 + (j-1)*b
	      x4 = 2.5*a +(i-1)*3*a
            y4 = y3

	      write(2,1002) x1,y1,z,atype
	      write(2,1002) x2,y2,z,atype
	      write(2,1002) x3,y3,z,atype
	      write(2,1000) x4,y4,z,atype

		  write(2,1002) x1,y1,zH,atype
	      write(2,1002) x2,y2,zH,atype
	      write(2,1002) x3,y3,zH,atype
	      write(2,1002) x4,y4,zH,atype

 1002	format(2x,3(2x,f12.6),2x,i3)

	   enddo
     	enddo
 	endif

* *** ***********************************************************************
      if(nSubstrate_PoreCentered.eq.1)then
* *** ***********************************************************************
*                   FINITE SLIT PORE GENERATION (PORE WIDTH = H)                                

      atype = 1                  !all atoms are C
	z = 0.0                    !graphite is flat!!!
	zH = z + H

      do j=1,Ny
	   do i=1,Nx

     	      x1 = a + (i-1)*3*a + edge		     !first line of atoms
	      y1 = (j-1)*b	+ edge
	      x2 = 2*a +(i-1)*3*a + edge
            y2 = y1

	      x3 = a/2 + (i-1)*3*a + edge		 !second line of atoms
	      y3 = b/2 + (j-1)*b + edge
	      x4 = 2.5*a +(i-1)*3*a	 + edge
            y4 = y3

	      write(3,1003) x1,y1,z,atype
	      write(3,1003) x2,y2,z,atype
	      write(3,1003) x3,y3,z,atype
	      write(3,1003) x4,y4,z,atype

		  write(3,1003) x1,y1,zH,atype
	      write(3,1003) x2,y2,zH,atype
	      write(3,1003) x3,y3,zH,atype
	      write(3,1003) x4,y4,zH,atype

 1003	format(2x,3(2x,f12.6),2x,i3)

	   enddo
     	enddo

	    y5 = Ny*b	+ edge		   !last line of atoms
	    y6= y5
	do i=1,Nx
		x5 = a + (i-1)*3*a + edge	
		x6 = 2*a +(i-1)*3*a + edge
 	   
	    write(3,1003) x5,y5,z,atype
          write(3,1003) x5,y5,zH,atype
		write(3,1003) x6,y6,z,atype
          write(3,1003) x6,y6,zH,atype
	enddo

 	endif
* *** **********************************************************************
      
	end