int c,m,n,l,res_num,  i, j, abv_i, npore, nsur, n_max, m_max, crgadd, s;
float radius, x_slide, xp, yp, zp, pore_rad, test_rad, sur_por, epsilon, xb, yb, y_max, ypore, xpore, ysur, xsur, z_slide;
char *atom_name, *res_name, chain_name;
main()
{
  
  c=1;
  radius=3.5;                      /* in Angstroms */
  pore_rad=19.0*radius;              /* enter here the radius of the pore=pore_rad, which is odd number times the radius */
s=(pore_rad/3.5 +1)/2;  /* this is to determine crgadd */
printf("HEADER 416.5x419.228x70 A^3 BLM with the pore of radius %.1fA\n",pore_rad);
  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 */
  epsilon=0.01;                    /* make sure x^2+y^2 <= test_rad^2 +epsilon */
  for(l=0;l<=9;++l) 
    {
      
    for(n= -34;n<=34;++n)
      {
        x_slide=0.0;          /* x position depends on y; that is why this is here */
        for(m= -29;m<=29;++m)
          {
           atom_name="bdy";
	   if(l==0 || l==9)   /* name top and bottom layer scg to be able to assign charge */ 
	     atom_name="scg";
           res_name="fil";
           chain_name='O';
           res_num=l+1;      /* res_num goes from 1 to 10 ie points to a layer */
	                     /* there are 40710 atoms total, 4071 in each "slice" */
           if(n%2 != 0)               /* test y position */
             x_slide=radius;          /* add radius to x position */
          
           xp=2*radius*m + x_slide;
           yp=1.7320508*radius*n;
           zp= -radius -l*2*radius;  /* -3.5 -10.5 .. -66.5 */

           /* x_goes[-203,203] and [-199.5,206.5] */
	   /* x_range[-203-3.5=-206.5, 206.5+3.5=210] */
	   /* x_max - x_min = x_slide = radius */
	   
	   /* y_goes[-206.114, 206.114] */
	   /* y_range[-s,s], s=209.614 */

           /* z_goes[-3.5, -66.5] */
	   /* z_range[0,-70] */

/* 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 by removing it ******/
	  
	   y_max=1.7320508*radius*34;  /* 34 is max(n) */
	   xb=2*radius*0; 
	   if(xb==0)
	     goto chargedont; /* add charges, don't remove them */
/* xb is x boundary within which atoms are marked at +-y_max; subject to change  */
    /* # of atoms marked(res_name="rem") at  +-y_max is (xb/7*2+1)*2*2 */ 

	   if( (zp == -3.5 || zp == -66.5) && (yp==y_max || yp== -y_max) && xp<xb+epsilon && xp> -xb-epsilon ) /* abs(xp) < xb + epsilon */
	     res_name="rem";


	   yb=1.7320508*radius*0;
/* yb is y boundary within which atoms are marked at +-x_max; subject to change  */
/* # of atoms marked(res_name="rem") at +-x_max is (yb/3.5/sqrt(3)*2+1)*2*2 */

	   if( (zp == -3.5 || zp == -66.5) && (xp == -203 || xp == 203 || xp == -199.5 || xp == 206.5) && yp< yb+epsilon && yp> -yb-epsilon )
	                                       /* abs(yp)< yb + epsilon */
	     res_name="rem";  
/* only top and bottom layer atoms (scg) are marked */ 

      /****** end of the routine that conserves charge ******/


	 chargedont:
	   
	   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 charge by adding it *****/


  crgadd = (3*s*s-3*s+1)*2-s*6*8; /* see notes */
  if(crgadd <= 0)
    goto quitcrg;
  m_max=29;
  n_max=34;
 charge_conserve:
  xp = -m_max*7.0-radius;  /* atoms are added at the left edge periodically */
  /* total # of atoms that can be added this way is (33+1)/2*2*8=272 */
  for(n=1;n<=33;n+=2){      /* yp is 'odd' */
    for(l=0;l<=1;++l){      /* atoms are added to layers 1 and 10 only */
      z_slide=0.0;
      if (l==1)
	z_slide=63.0;
      yp=1.7320508*radius*n;
      zp= -radius -z_slide; /* -3.5 or -66.5 */
      atom_name="scg";
      res_name="put";
      res_num=1;
      if(l==1)
	res_num=10;
   /* add atoms symmetricly around x; that is why there are two printf's */
      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; crgadd-=1;
      if (c-1 == 40710)
	goto end;
      if (crgadd == 0)
	goto quitcrg;
	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; crgadd-=1;
      if (c-1 == 40710)
      goto end;
      if (crgadd == 0)
	goto quitcrg;
    }
  }


  /* routine that conserves charge when # of atoms that should be added is bigger than 272; it goes up to 272+280=552 atoms */
  xp = (m_max+1)*7.0; /* atoms are added at the right edge periodically */
  /* total # of atoms that can be added this way is (34)/2*2*8+8=280 */
  for(n=0;n<=34;n+=2){      /* yp is 'even' */
    for(l=0;l<=1;++l){      /* atoms are added to layers 1 and 10 only */
      z_slide=0.0;
      if (l==1)
	z_slide=63.0;
      yp=1.7320508*radius*n;
      zp= -radius -z_slide; /* -3.5 or -66.5 */
      atom_name="scg";
      res_name="put";
      res_num=1;
      if(l==1)
	res_num=10;
/* add atoms symmetricly around x; that is why there are two printf's */
      if(n==0)
	goto cdontprint0twice;
      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; crgadd-=1;
      if (c-1 == 40710)
	goto end;
      if (crgadd == 0)
	goto quitcrg;
    cdontprint0twice:
      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; crgadd-=1;
      if (c-1 == 40710)
	goto end;
      if (crgadd == 0)
	goto quitcrg;
      
    }
  }

  /* routine that conserves charge when # of atoms that should be added is bigger than 552; it goes up to 552+480=1032 atoms */
  yp = radius*1.7320508*(-n_max-1); /* atoms are added at the bottom edge */
  x_slide=0.0;          
  if((n_max+1)%2 != 0)
    x_slide=radius; /* add radius to x position when n_max+1 is odd */
  /* total # of atoms that can be added this way is (30)*2*8=480 */
  for(m=0;m<=29;++m){     
    for(l=0;l<=1;++l){      /* atoms are added to layers 1 and 10 only */
      z_slide=0.0;
      if (l==1)
	z_slide=63.0;
      xp=2*radius*m + x_slide;
      zp= -radius -z_slide; /* -3.5 or -66.5 */
      atom_name="scg";
      res_name="put";
      res_num=1;
      if(l==1)
	res_num=10;
/* add atoms symmetricly around y; that is why there are two printf's */
    
      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; crgadd-=1;
      if (c-1 == 40710)
	goto end;
      if (crgadd == 0)
	goto quitcrg;

      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; crgadd-=1;
      if (c-1 == 40710)
	goto end;
      if (crgadd == 0)
	goto quitcrg;
      
    }
  }

  /* routine that conserves charge when # of atoms that should be added is bigger than 1032; it goes up to 1032+280=1312 atoms */
  xp = -(m_max+1)*7.0; /* atoms are added at the left edge periodically */
  /* total # of atoms that can be added this way is (34)/2*2*8+8=280 */
  for(n=0;n<=34;n+=2){      /* yp is 'even' */
    for(l=0;l<=1;++l){      /* atoms are added to layers 1 and 10 only */
      z_slide=0.0;
      if (l==1)
	z_slide=63.0;
      yp=1.7320508*radius*n;
      zp= -radius -z_slide; /* -3.5 or -66.5 */
      atom_name="scg";
      res_name="put";
      res_num=1;
      if(l==1)
	res_num=10;
/* add atoms symmetricly around x; that is why there are two printf's */
      if(n==0)
	goto cnot0twice;
      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; crgadd-=1;
      if (c-1 == 40710)
	goto end;
      if (crgadd == 0)
	goto quitcrg;
    cnot0twice:
      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; crgadd-=1;
      if (c-1 == 40710)
	goto end;
      if (crgadd == 0)
	goto quitcrg;
      
    }
  }


 /* routine that conserves charge when # of atoms that should be added is bigger than 1312; it goes up to 1312+272=1584 atoms */
  xp = (m_max+1)*7.0+radius; /* atoms are added at the right edge periodically */
  /* total # of atoms that can be added this way is (33+1)/2*2*8=272 */
  for(n=1;n<=33;n+=2){      /* yp is 'odd' */
    for(l=0;l<=1;++l){      /* atoms are added to layers 1 and 10 only */
      z_slide=0.0;
      if (l==1)
	z_slide=63.0;
      yp=1.7320508*radius*n;
      zp= -radius -z_slide; /* -3.5 or -66.5 */
      atom_name="scg";
      res_name="put";
      res_num=1;
      if(l==1)
	res_num=10;
/* add atoms symmetricly around x; that is why there are two printf's */
      if(n==0)
	goto cnottwice;
      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; crgadd-=1;
      if (c-1 == 40710)
	goto end;
      if (crgadd == 0)
	goto quitcrg;
    cnottwice:
      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; crgadd-=1;
      if (c-1 == 40710)
	goto end;
      if (crgadd == 0)
	goto quitcrg;
      
    }
  }
 
 
  /* routine that conserves charge when # of atoms that should be added is bigger than 1584; it goes up to 1584+480=2064 atoms */
  yp = radius*1.7320508*(n_max+1); /* atoms are added at the top edge */
  x_slide=0.0;          /* add radius to x position because -35 is odd */
  if((n_max+1)%2 != 0)
    x_slide=radius;
  /* total # of atoms that can be added this way is (30)*2*8=480 */
  for(m=0;m<=29;++m){      /* m goes all the way to mmax */
    for(l=0;l<=1;++l){      /* atoms are added to layers 1 and 10 only */
      z_slide=0.0;
      if (l==1)
	z_slide=63.0;
      xp=2*radius*m + x_slide;
      zp= -radius -z_slide; /* -3.5 or -66.5 */
      atom_name="scg";
      res_name="put";
      res_num=1;
      if(l==1)
	res_num=10;
/* add atoms symmetricly around y; that is why there are two printf's */
    
      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; crgadd-=1;
      if (c-1 == 40710)
	goto end;
      if (crgadd == 0)
	goto quitcrg;

      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; crgadd-=1;
      if (c-1 == 40710)
	goto end;
      if (crgadd == 0)
	goto quitcrg;
      
    }
  }
  m_max+=1;
  n_max+=1;
  goto charge_conserve; 
 quitcrg:
  ;

  
            /****** routine to conserve mass *******/
  m_max=29;
  n_max=34;
 mass_conserve:
  xp = -m_max*7.0-radius;  /* atoms are added at the left edge periodically */
 /* total # of atoms that can be added this way is (33+1)/2*2*8=272 */
  for(n=1;n<=33;n+=2){      /* yp is 'odd' */
    for(l=1;l<=8;++l){      /* atoms are added to layers 2-9 only */
      yp=1.7320508*radius*n;
      zp= -radius -l*2*radius;  /* -10.5 .. -59.5 */
      atom_name="add";
      res_name="mas";
      res_num=l+1;
   /* add atoms symmetricly around x; that is why there are two printf's */
      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-1 == 40710)
	goto end;
	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-1 == 40710)
      goto end;
      
    }
  }


  /* routine that conserves mass when # of atoms that should be added is bigger than 272; it goes up to 272+280=552 atoms */
  xp = (m_max+1)*7.0; /* atoms are added at the right edge periodically */
  /* total # of atoms that can be added this way is (34)/2*2*8+8=280 */
  for(n=0;n<=34;n+=2){      /* yp is 'even' */
    for(l=1;l<=8;++l){      /* atoms are added to layers 2-9 only */
      yp=1.7320508*radius*n;
      zp= -radius -l*2*radius;  /* -10.5 .. -59.5 */
      atom_name="add";
      res_name="mas";
      res_num=l+1;
/* add atoms symmetricly around x; that is why there are two printf's */
      if(n==0)
	goto dontprint0twice;
      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-1 == 40710)
	goto end;
    dontprint0twice:
      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-1 == 40710)
	goto end;
      
    }
  }

  /* routine that conserves mass when # of atoms that should be added is bigger than 552; it goes up to 552+480=1032 atoms */
  yp = radius*1.7320508*(-n_max-1); /* atoms are added at the bottom edge */
  x_slide=0.0;          
  if((n_max+1)%2 != 0)
    x_slide=radius; /* add radius to x position when n_max+1 is odd */
  /* total # of atoms that can be added this way is (30)*2*8=480 */
  for(m=0;m<=29;++m){     
    for(l=1;l<=8;++l){      /* atoms are added to layers 2-9 only */
      xp=2*radius*m + x_slide;
      zp= -radius -l*2*radius;  /* -10.5 .. -59.5 */
      atom_name="add";
      res_name="mas";
      res_num=l+1;
/* add atoms symmetricly around y; that is why there are two printf's */
    
      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-1 == 40710)
	goto end;

      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-1 == 40710)
	goto end;
      
    }
  }

  /* routine that conserves mass when # of atoms that should be added is bigger than 1032; it goes up to 1032+280=1312 atoms */
  xp = -(m_max+1)*7.0; /* atoms are added at the left edge periodically */
  /* total # of atoms that can be added this way is (34)/2*2*8+8=280 */
  for(n=0;n<=34;n+=2){      /* yp is 'even' */
    for(l=1;l<=8;++l){      /* atoms are added to layers 2-9 only */
      yp=1.7320508*radius*n;
      zp= -radius -l*2*radius;  /* -10.5 .. -59.5 */
      atom_name="add";
      res_name="mas";
      res_num=l+1;
/* add atoms symmetricly around x; that is why there are two printf's */
      if(n==0)
	goto not0twice;
      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-1 == 40710)
	goto end;
    not0twice:
      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-1 == 40710)
	goto end;
      
    }
  }


 /* routine that conserves mass when # of atoms that should be added is bigger than 1312; it goes up to 1312+272=1584 atoms */
  xp = (m_max+1)*7.0+radius; /* atoms are added at the right edge periodically */
  /* total # of atoms that can be added this way is (33+1)/2*2*8=272 */
  for(n=1;n<=33;n+=2){      /* yp is 'odd' */
    for(l=1;l<=8;++l){      /* atoms are added to layers 2-9 only */
      yp=1.7320508*radius*n;
      zp= -radius -l*2*radius;  /* -10.5 .. -59.5 */
      atom_name="add";
      res_name="mas";
      res_num=l+1;
/* add atoms symmetricly around x; that is why there are two printf's */
      if(n==0)
	goto nottwice;
      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-1 == 40710)
	goto end;
    nottwice:
      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-1 == 40710)
	goto end;
      
    }
  }
 
 
  /* routine that conserves mass when # of atoms that should be added is bigger than 1584; it goes up to 1584+480=2064 atoms */
  yp = radius*1.7320508*(n_max+1); /* atoms are added at the top edge */
  x_slide=0.0;          /* add radius to x position because -35 is odd */
  if((n_max+1)%2 != 0)
    x_slide=radius;
  /* total # of atoms that can be added this way is (30)*2*8=480 */
  for(m=0;m<=29;++m){      /* m goes all the way to mmax */
    for(l=1;l<=8;++l){      /* atoms are added to layers 2-9 only */
      xp=2*radius*m + x_slide;
      zp= -radius -l*2*radius;  /* -10.5 .. -59.5 */
      atom_name="add";
      res_name="mas";
      res_num=l+1;
/* add atoms symmetricly around y; that is why there are two printf's */
    
      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-1 == 40710)
	goto end;

      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-1 == 40710)
	goto end;
      
    }
  }
  m_max+=1;
  n_max+=1;
  goto mass_conserve; 
 end:
  ;
}