program grid_slit
! *** ************************************************************
!        CALCULATE 3D GRID of H2-WALL INTERACTION IN SLIT PORE
!        (Explicit positions of atoms in the walls are nedded)
!        (Feymann-Hibbs quantum correction included)
!	 
!                           Version Feb 09, 2011
!                 
!                  
!        You need: 1. Natom (number of atoms in the walls)
!                  2.T
!                  3 input file (box dimensions,gridstep,ctof,xyz)
!
! ****************************************************************
!      VARIABLES
! ****************************************************************
	parameter (Natom=832,NB=16)
!		Natom - number of atoms in the walls
      
	parameter (epsC=49.,sigmaC=3.34)			!  LJ parameters

	parameter (epsCa=100.0,sigmaCa=3.4)
	parameter (epsCb=350.0,sigmaCb=3.4)
	parameter (epsB=550.0,sigmaB=3.54)

	parameter (epsH2=34.2, sigmaH2=2.96)

	parameter (T=77.)
	parameter (Q=1.0059)

	parameter (nperiod =1)

	integer ntype(Natom)

 	real xbox,ybox,zbox,limit,step,ctof
	real x(Natom),y(Natom),z(Natom)

	real r2 
	
	open(1,file='Pore_8x13_H7.dat',status='old')
	open(2,file='Grid_8x13.dat',status='unknown')

* *** ***************************************************
*     READ-IN slit, grid and atoms'coordinates

	read(1,*)  xbox,ybox,zbox,limit

	read(1,*)  step,ctof

	do Ia=1,Natom
	   read(1,*) x(Ia),y(Ia),z(Ia),ntype(Ia)
	enddo  
* *** ***************************************************
!      INITIALISATIONS

	pi=4*atan(1.)
	max=1e6
	ctof2=ctof*ctof
	
      Nx=int(xbox/step)+1
	Ny=int(ybox/step)+1
	Nz=int((zbox-2*limit)/step)+1

* *** ***************************************************
!     HERE STARTS THE LOOP OF GRID CALCULATIONS

	do Ix=0,Nx
	   xg=Ix * step
	   print*, Ix

	   do Iy=0,Ny
	      yg=Iy * step

	      do Iz=0,Nz
	         zg=limit + Iz * step

	         Exyz = 0

			 do Ia=1,Natom
	            rz=z(Ia) - zg

	            do Iper = -1,1
				   dx = Iper * xbox
	               rx = x(Ia) + dx - xg

	            do Iper1 = -1,1
	               dy = Iper1 * ybox
				   ry =	y(Ia) + dy - yg	   			   	 

				   r2 =rx*rx + ry*ry + rz*rz
				   if(r2.le.ctof2)then

	                 Itrans=ntype(Ia)   

						call pot(E,r2,epsC,epsB,epsH2,sigmaC,sigmaB,
     1				 sigmaH2,Itrans,T,Q,epsCa,epsCb,sigmaCa,sigmaCb)

						Exyz = Exyz + E
				   endif
	             enddo
				enddo
			  enddo
		    
          	if(Exyz.gt.max) Exyz = max
			write(2,'(i4,2x,i4,2x,i4,2x,e14.8)') Ix,Iy,Iz,Exyz
 
		  enddo
	   enddo
	enddo

	close(1)
	close(2)

	end
* *** **********************************************
	subroutine pot(E,r2,epsC,epsB,epsH2,sigmaC,sigmaB,
	1   sigmaH2,Itrans,T,Q,epsCa,epsCb,sigmaCa,sigmaCb)
	r6 = r2*r2*r2
  	r12 = r6*r6

	if(Itrans.eq.1) then                       ! atom is carbon
	    eps = sqrt(epsH2 *epsC)
 	    sigma = (sigmaH2 + sigmaC)/2
	  elseif(Itrans.eq.2) then				   ! atom is boron
		eps = sqrt(epsH2 *epsB)
 	    sigma = (sigmaH2 + sigmaB)/2
	  elseif(Itrans.eq.3) then				   ! atom is alpha carbon
		eps = sqrt(epsH2 *epsCa)
 	    sigma = (sigmaH2 + sigmaCa)/2
	  else				                       ! atom is beta carbon
		eps = sqrt(epsH2 *epsCb)
 	    sigma = (sigmaH2 + sigmaCb)/2
	endif

 	sigma6 = sigma*sigma*sigma*sigma*sigma*sigma
	sigma12 = sigma6*sigma6

   	   CLJsf=4*eps*(sigma12/r12 - sigma6/r6)

         derivLJsf = 4*eps*(-12*sigma12/r12 + 6*sigma6/r6)*(1/sqrt(r2))
	   dderivLJsf = 4*eps*(12*13*sigma12/r12 - 6*7*sigma6/r6)/r2	

	   corr1LJsf = dderivLJsf
	   corr2LJsf = 2*derivLJsf/sqrt(r2)

	   E = CLJsf + ((Q/2)/T)*(corr1LJsf + corr2LJsf)   !for sf, Q = Q/2

	end