int c, m, n, abv_n, l, res_num, npore, abv_i, j, i, nlay, nsur, nbound, mas_tot;

float radius, xp, yp, zp, pore_rad, lay_rad, test_rad, sur_por, xsur, ysur, xpore, ypore, ybound, xbound;

char *atom_name, *res_name, chain_name;

main()
{
  c=1;                      /* counts atoms */
  radius=3.5;                /*in Angstroms*/
  lay_rad=52.5;            /*enter the radius of the layer here!*/
  pore_rad=24.5;            /*enter the radius of the pore here!*/
  test_rad = pore_rad - radius;      /* identify centers of atmos that make a pore */
  sur_por=test_rad + 2.0*radius;    /* identify atoms surronding the pore */
  
  
  for(l=0;l<=9;++l)
    {
      nlay=(lay_rad/radius-1)/2;
      for(n=nlay;n>= -nlay;--n){
        yp=radius*1.7320508*n;
        abv_n=n;
        if(n<0)
          abv_n= -n;
        for(m=2*nlay-abv_n;m>= -2*nlay+abv_n;m=m-2){
	  
          xp=m*radius;
          zp= -radius -l*2*radius + 35.0;  /* -31.5 ... 31.5 */
	  
	  
	  atom_name="bdy";
	  res_name="fil";
	  chain_name='O';	 
	  if(l==0 || l==9)       /* name top and bottom layer scg -- */
	    atom_name="scg";     /*-- in order to be able to assign charge */
	  res_num=l+1;  /* res_num goes from 1 to 10 i.e. points to a layer */


	  /* new routine that removes atoms that make a pore and marks atoms surrounding the pore 8/22/92 */

	  npore=(pore_rad/radius-1)/2;
	  for(i=npore;i>= -npore;--i){
	    ypore=radius*1.7320508*i;
	    abv_i=i;
	    if(i<0)
	      abv_i= -i;
	    for(j=2*npore-abv_i;j>= -2*npore+abv_i;--j){
	      xpore=j*radius;
	      if(xp==xpore && yp==ypore)
		goto skip;  /* don't print atoms that make up a pore of radius (test_rad + radius) */
	    }
	  }
	  nsur=npore+1;
	  for(i=nsur;i>= -nsur;--i){
	    ysur=radius*1.7320508*i;
	    abv_i=i;
	    if(i<0)
	      abv_i= -i;
	    for(j=2*nsur-abv_i;j>= -2*nsur+abv_i;--j){
	      xsur=j*radius;
	      if(xp==xsur && yp==ysur)
		{
		  atom_name="por";
		  res_name="sur";
	 	}
	    }
	  }



      /****** routine to conserve charge  6/8/93  ******/

  /* removing charges from the +-y_edge; change j=j-4 to j=j-2 if needed */
	  for(i=nlay;i >= nlay-1; --i){        
	    ybound=radius*1.7320508*i;
	    abv_i=i;
	    if(i<0)
	      abv_i= -i;
	    for(j=2*nlay-abv_i;j>= -2*nlay+abv_i;j=j-2){
	      xbound=j*radius;
	      if( xp==xbound  && (yp==ybound || yp== -ybound)   && (l==0 || l==9) )
		{
		  atom_name="scg";
		  res_name="rem";
		}
	    }
	  }
  /* removing charges from the other 4 edges; change i=i+2 to i=i+1 if needed*/
	  for(i=0; i < nlay; i=i+1){ 
	    ybound=radius*1.7320508*i;	    
	    j=2*nlay-i;
	    xbound=j*radius;
	    if( (xp==xbound || xp== - xbound) && (yp==ybound || yp== -ybound)   && (l==0 || l==9) )
	      {
		atom_name="scg";
		res_name="rem";
	      }
	  }
	  
      


	  printf("ATOM%7d %-5s%3s%2c%4d%11.2f%9.3f%8.3f\n",   c, atom_name, res_name, chain_name, res_num, xp, yp, zp);
	  c+=1;
	skip:
	  ;   
	}
      }
    }


  /***** routine to conserve mass 8/6/93 *****/


  /* adding mass on the top of the +-y_edge; change j=j-4 to j=j-2 if needed */
  mas_tot=1690; /* total # of atoms for 52.5A hexagonal layer */
  i=nlay+1;        
  yp=radius*1.7320508*i;
  abv_i=i;
  if(i<0)
    abv_i= -i;
  for(l= -4;l<=3;++l){                 /* change l here */
    for(j=0;j<=2*i-abv_i-2;j=j+2){
      xp=j*radius;
      zp= radius +l*2*radius;  /** -24.5 ... 24.5; mas canot go in either of the top layers **/
      atom_name="add";
      res_name="mas";
      printf("ATOM%7d %-5s%3s%2c%4d%11.2f%9.3f%8.3f\n",   c, atom_name, res_name, chain_name, res_num, xp, yp, zp); 
      c+=1;
      printf("ATOM%7d %-5s%3s%2c%4d%11.2f%9.3f%8.3f\n",   c, atom_name, res_name, chain_name, res_num, xp, -yp, zp); 
      c+=1;
      if(xp==0)
	goto not;
      printf("ATOM%7d %-5s%3s%2c%4d%11.2f%9.3f%8.3f\n",   c, atom_name, res_name, chain_name, res_num, -xp, yp, zp); 
      c+=1;
      printf("ATOM%7d %-5s%3s%2c%4d%11.2f%9.3f%8.3f\n",   c, atom_name, res_name, chain_name, res_num, -xp, -yp, zp); 
      c+=1;
      if(c > mas_tot)
	goto end;
    not:
      ;
    }
  }
  /* adding mas at the other 4 edges; change i=i+2 to i=i+1 if needed*/
  for(l= -4;l<= 3;++l){                 /* change l here */  
    zp= radius +l*2*radius;  /** -24.5 ... 24.5; mas canot go in either of the top layers **/ 
    for(i=0; i <= nlay+1; i=i+1){ 
      yp=radius*1.7320508*i;	    
      j=2*(nlay+1)-i;
      xp=j*radius;
      atom_name="add";
      res_name="mas";
      printf("ATOM%7d %-5s%3s%2c%4d%11.2f%9.3f%8.3f\n",   c, atom_name, res_name, chain_name, res_num, xp, yp, zp); 
      c+=1;
      printf("ATOM%7d %-5s%3s%2c%4d%11.2f%9.3f%8.3f\n",   c, atom_name, res_name, chain_name, res_num, -xp, yp, zp); 
      c+=1;
      if(yp==0)
	goto not1;
      printf("ATOM%7d %-5s%3s%2c%4d%11.2f%9.3f%8.3f\n",   c, atom_name, res_name, chain_name, res_num, xp, -yp, zp); 
      c+=1;
      printf("ATOM%7d %-5s%3s%2c%4d%11.2f%9.3f%8.3f\n",   c, atom_name, res_name, chain_name, res_num, -xp, -yp, zp); 
      c+=1;
      if(c > mas_tot)
	goto end;   
    not1:
      ;
    }
  }



 end:
  ;
}