src/contrib/mkice.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2006, The GROMACS development team,
   6  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012,2013, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #include <stdio.h>
  43 #include <math.h>
  44 #include "typedefs.h"
  45 #include "statutil.h"
  46 #include "copyrite.h"
  47 #include "gmx_fatal.h"
  48 #include "pdbio.h"
  49 #include "macros.h"
  50 #include "smalloc.h"
  51 #include "vec.h"
  52 #include "pbc.h"
  53 #include "physics.h"
  54 #include "names.h"
  55 #include "txtdump.h"
  56 #include "trnio.h"
  57 #include "symtab.h"
  58 #include "strdb.h"
  59 #include "atomprop.h"
  60 #include "confio.h"
  61
  62 #define TET   109.47
  63 #define DCONS 0.117265878
  64
  65 typedef struct {
  66   int n,aa[4];
  67 } t_bbb;
  68
  69 static char *watname[] = { "OW ", "HW1", "HW2", "DW", "SW" };
  70 static char *diamname[] = { "C1", "H2", "H3", "H4", "H5", "H2", "H3", "H4", "H5" };
  71 static real qspc[]     = { -0.82, 0.41, 0.41 };
  72 static real qyaw[]     = { 1.24588, 0.62134, 0.62134, 0.0, -2.48856 };
  73 static real spc_lj[3][6] = {
  74   { 2.6171e-3, 2.6331e-6, 0, 0, 0, 0 },
  75   { 0,      0,      0, 0, 0, 0 },
  76   { 0,      0,      0, 0, 0, 0 }
  77 };
  78 #define CHH 2e-9
  79 #define CHS 2e-9
  80 #define COS 2e-6
  81 static real yaw_lj[5][10] = {
  82   { 0, 0,   0, 0,   0,   0, 0, 0, 0,      COS },
  83   { 0, 0,   0, CHH, 0, CHH, 0, 0, 0,      CHS },
  84   { 0, 0,   0, CHH, 0, CHH, 0, 0, 0,      CHS },
  85   { 0, 0,   0, 0,   0,   0, 0, 0, 0,      0   },
  86   { 0, COS, 0, CHS, 0, CHS, 0, 0, 2.6e-3, 0   }
  87 };
  88
  89 void unitcell(rvec x[],rvec box,gmx_bool bYaw,real odist,real hdist)
  90 {
  91 #define cx  0.81649658
  92 #define cy  0.47140452
  93 #define cy2 0.94280904
  94 #define cz  0.33333333
  95
  96   rvec xx[24] = {
  97     { 0,   0,         0 }, /* O1 */
  98     { 0,   0,         1 }, /* H relative to Oxygen */
  99     { cx, cy,       -cz },
 100     { cx, cy,       -cz }, /* O2 */
 101     { 0, 0,       -1    }, /* H relative to Oxygen */
 102     { cx,-cy,       +cz },
 103     { cx, cy+cy2,     0 }, /* O3 */
 104     { -cx, cy,    -cz   }, /* H relative to Oxygen */
 105     { 0,   -cy2,    -cz },
 106     { 0,  2*cy+cy2, -cz }, /* O4 */
 107     {-cx,-cy,       +cz }, /* H relative to Oxygen */
 108     { 0 , cy2,      +cz },
 109     { 0,   0,         1 }, /* O5 */
 110     {-cx, cy,       +cz }, /* H relative to Oxygen */
 111     { 0 , -cy2,     +cz },
 112     { cx, cy,      1+cz }, /* O6 */
 113     { -cx, -cy,     -cz }, /* H relative to Oxygen */
 114     { 0,   cy2,     -cz },
 115     { cx, cy+cy2,     1 }, /* O7 */
 116     { 0,  0,       -1   }, /* H relative to Oxygen */
 117     { cx, cy,       +cz },
 118     { 0,  2*cy+cy2,1+cz }, /* O8 */
 119     { 0,   0,         1 }, /* H relative to Oxygen */
 120     { cx,   -cy,    -cz }
 121   };
 122   int  i,iin,iout,j,m;
 123   rvec tmp,t2,dip;
 124
 125   clear_rvec(dip);
 126   for(i=0; (i<8); i++) {
 127     iin = 3*i;
 128     if (bYaw)
 129       iout = 5*i;
 130     else
 131       iout = iin;
 132     svmul(odist,xx[iin],x[iout]);
 133     svmul(-0.82,x[iout],t2);
 134     rvec_inc(dip,t2);
 135     for(j=1; (j<=2); j++) {
 136       svmul(hdist,xx[iin+j],tmp);
 137       rvec_add(x[iout],tmp,x[iout+j]);
 138       svmul(0.41,x[iout+j],t2);
 139       rvec_inc(dip,t2);
 140     }
 141     if (bYaw)
 142       for(m=0; (m<DIM); m++)
 143         x[iout+3][m] = x[iout+4][m] =
 144           (1-2*DCONS)*x[iout][m]+DCONS*(x[iout+1][m]+x[iout+2][m]);
 145   }
 146
 147   box[XX] = 2*cx;
 148   box[YY] = 2*(cy2+cy);
 149   box[ZZ] = 2*(1+cz);
 150   for(i=0; (i<DIM); i++)
 151     box[i] *= odist;
 152
 153   printf("Unitcell:  %10.5f  %10.5f  %10.5f\n",box[XX],box[YY],box[ZZ]);
 154   printf("Dipole:    %10.5f  %10.5f  %10.5f (e nm)\n",dip[XX],dip[YY],dip[ZZ]);
 155 }
 156
 157 void random_h_coords(int natmol,int nmol,rvec x[],rvec box,
 158                      gmx_bool bYaw,real odist,real hdist)
 159 {
 160 #define cx  0.81649658
 161 #define cy  0.47140452
 162 #define cy2 0.94280904
 163 #define cz  0.33333333
 164
 165   rvec xx[24] = {
 166     { 0,   0,         0 }, /* O1 */
 167     { 0,   0,         1 }, /* H relative to Oxygen */
 168     { cx, cy,       -cz },
 169     { cx, cy,       -cz }, /* O2 */
 170     { 0, 0,       -1    }, /* H relative to Oxygen */
 171     { cx,-cy,       +cz },
 172     { cx, cy+cy2,     0 }, /* O3 */
 173     { -cx, cy,    -cz   }, /* H relative to Oxygen */
 174     { 0,   -cy2,    -cz },
 175     { 0,  2*cy+cy2, -cz }, /* O4 */
 176     {-cx,-cy,       +cz }, /* H relative to Oxygen */
 177     { 0 , cy2,      +cz },
 178     { 0,   0,         1 }, /* O5 */
 179     {-cx, cy,       +cz }, /* H relative to Oxygen */
 180     { 0 , -cy2,     +cz },
 181     { cx, cy,      1+cz }, /* O6 */
 182     { -cx, -cy,     -cz }, /* H relative to Oxygen */
 183     { 0,   cy2,     -cz },
 184     { cx, cy+cy2,     1 }, /* O7 */
 185     { 0,  0,       -1   }, /* H relative to Oxygen */
 186     { cx, cy,       +cz },
 187     { 0,  2*cy+cy2,1+cz }, /* O8 */
 188     { 0,   0,         1 }, /* H relative to Oxygen */
 189     { cx,   -cy,    -cz }
 190   };
 191   int  i,iin,iout,j,m;
 192   rvec tmp,t2,dip;
 193
 194   clear_rvec(dip);
 195   for(i=0; (i<nmol); i++) {
 196     iin = natmol*i;
 197     iout = iin;
 198     svmul(odist,x[iin],x[iout]);
 199     svmul(-0.82,x[iout],t2);
 200     rvec_inc(dip,t2);
 201     for(j=1; (j<=2); j++) {
 202       svmul(hdist,xx[3*(i % 8)+j],tmp);
 203       rvec_add(x[iout],tmp,x[iout+j]);
 204       svmul(0.41,x[iout+j],t2);
 205       rvec_inc(dip,t2);
 206     }
 207   }
 208
 209   box[XX] = 2*cx;
 210   box[YY] = 2*(cy2+cy);
 211   box[ZZ] = 2*(1+cz);
 212   for(i=0; (i<DIM); i++)
 213     box[i] *= odist;
 214
 215   printf("Unitcell:  %10.5f  %10.5f  %10.5f\n",box[XX],box[YY],box[ZZ]);
 216   printf("Dipole:    %10.5f  %10.5f  %10.5f (e nm)\n",dip[XX],dip[YY],dip[ZZ]);
 217 }
 218
 219 void unitcell_d(rvec x[],rvec box,real odist)
 220 {
 221   rvec cc[8] = {
 222     { 0,   0,         0 }, /* C1 */
 223     { cx, cy,       -cz }, /* C2 */
 224     { cx, cy+cy2,     0 }, /* C3 */
 225     { 0,  2*cy+cy2, -cz }, /* C4 */
 226     { 0,   0,         1 }, /* C5 */
 227     { cx, cy,      1+cz }, /* C6 */
 228     { cx, cy+cy2,     1 }, /* C7 */
 229     { 0,  2*cy+cy2,1+cz }, /* C8 */
 230   };
 231   rvec hh[4] = {
 232     { 0,   0,         1  }, /* H relative to C */
 233     { cx,  cy,       -cz },
 234     { cx, -cy,       -cz },
 235     {-cy2,  0,       -cz }
 236   };
 237   int  i,iin,iout,j,m;
 238   rvec tmp,t2,dip;
 239
 240   clear_rvec(dip);
 241   for(i=0; (i<8); i++) {
 242     iin  = i;
 243     iout = i;
 244     svmul(odist,cc[iin],x[iout]);
 245   }
 246   box[XX] = 2*cx;
 247   box[YY] = 2*(cy2+cy);
 248   box[ZZ] = 2*(1+cz);
 249   for(i=0; (i<DIM); i++)
 250     box[i] *= odist;
 251
 252   printf("Unitcell:  %10.5f  %10.5f  %10.5f\n",box[XX],box[YY],box[ZZ]);
 253 }
 254
 255 static t_bbb *mk_bonds(int natoms,rvec x[],real odist,
 256                        gmx_bool bPBC,matrix box)
 257 {
 258   real  od2 = odist*odist+1e-5;
 259   t_pbc pbc;
 260   t_bbb *bbb;
 261   int   i,j;
 262   rvec  dx;
 263
 264   if (bPBC)
 265     set_pbc(&pbc,box);
 266   snew(bbb,natoms);
 267   for(i=0; (i<natoms); i++) {
 268     for(j=i+1; (j<natoms); j++) {
 269       if (bPBC)
 270         pbc_dx(&pbc,x[i],x[j],dx);
 271       else
 272         rvec_sub(x[i],x[j],dx);
 273       if (iprod(dx,dx) <= od2) {
 274         bbb[i].aa[bbb[i].n++] = j;
 275         bbb[j].aa[bbb[j].n++] = i;
 276       }
 277     }
 278   }
 279   if (debug)
 280 #define PRB(nn) (bbb[(i)].n >= (nn)) ? bbb[i].aa[nn-1] : -1
 281     for(i=0; (i<natoms); i++)
 282       fprintf(debug,"bonds from %3d:  %d %d %d %d\n",
 283               i,PRB(1),PRB(2),PRB(3),PRB(4));
 284 #undef PRB
 285   return bbb;
 286 }
 287
 288 static void mk_diamond(t_atoms *a,rvec x[],real odist,t_symtab *symtab,
 289                        gmx_bool bPBC,matrix box)
 290 {
 291
 292   int   i,ib,j,k,l,m,nrm=0;
 293   t_bbb *bbb;
 294   gmx_bool  *bRemove;
 295   rvec  dx;
 296
 297   do {
 298     nrm = 0;
 299     bbb = mk_bonds(a->nr,x,odist,bPBC,box);
 300
 301     for(i=0; (i<a->nr); i++) {
 302       if (bbb[i].n < 2) {
 303         for(k=0; (k<bbb[i].n); k++) {
 304           ib = bbb[i].aa[k];
 305           for(j=0; (j<bbb[ib].n); j++)
 306             if (bbb[ib].aa[j] == i)
 307               break;
 308           if (j == bbb[ib].n)
 309             gmx_fatal(FARGS,"Bond inconsistency (%d not in list of %d)!\n",i,ib);
 310           for( ; (j<bbb[ib].n-1); j++)
 311             bbb[ib].aa[j] = bbb[ib].aa[j+1];
 312           bbb[ib].n--;
 313           nrm++;
 314         }
 315         bbb[i].n = 0;
 316       }
 317     }
 318
 319     for(i=j=0; (i<a->nr); i++) {
 320       if (bbb[i].n >= 2) {
 321         copy_rvec(x[i],x[j]);
 322         j++;
 323       }
 324     }
 325     fprintf(stderr,"Kicking out %d carbon atoms (out of %d)\n",
 326             a->nr-j,a->nr);
 327     a->nr = j;
 328     sfree(bbb);
 329   } while (nrm > 0);
 330   /* Rename atoms */
 331   bbb = mk_bonds(a->nr,x,odist,bPBC,box);
 332   for(i=0; (i<a->nr); i++) {
 333     switch (bbb[i].n) {
 334     case 4:
 335       a->atomname[i] = put_symtab(symtab,"C");
 336       break;
 337     case 3:
 338       a->atomname[i] = put_symtab(symtab,"CH1");
 339       break;
 340     case 2:
 341       a->atomname[i] = put_symtab(symtab,"CH2");
 342       break;
 343     default:
 344       gmx_fatal(FARGS,"This atom (%d) has %d bonds only",i,bbb[i].n);
 345     }
 346   }
 347   sfree(bbb);
 348 }
 349
 350 static real calc_ener(real c6,real c12,rvec dx,tensor vir)
 351 {
 352   real r,e,f;
 353   int  m,n;
 354
 355   r  = norm(dx);
 356   e  = c12*pow(r,-12.0) - c6*pow(r,-6.0);
 357   f  = 12*c12*pow(r,-14.0) - 6*c6*pow(r,-8.0);
 358   for(m=0; (m<DIM); m++)
 359     for(n=0; (n<DIM); n++)
 360       vir[m][n] -= 0.5*f*dx[m]*dx[n];
 361
 362   return e;
 363 }
 364
 365 static int read_rel_coords(char *fn,rvec **xx,int natmol)
 366 {
 367   int    i,nline;
 368   double x,y,z;
 369   char   **strings=NULL;
 370
 371   nline = get_file(fn,&strings);
 372   printf("Read %d lines from %s\n",nline,fn);
 373   snew(*xx,nline*natmol);
 374   for(i=0; (i<nline); i++) {
 375     if (sscanf(strings[i],"%lf%lf%lf",&x,&y,&z) != 3)
 376       gmx_fatal(FARGS,"Not enough arguments on line %d of file %s (should be 3)",i,fn);
 377     (*xx)[natmol*i][XX] = x;
 378     (*xx)[natmol*i][YY] = y;
 379     (*xx)[natmol*i][ZZ] = z;
 380   }
 381   return natmol*nline;
 382 }
 383
 384 void virial(FILE *fp,gmx_bool bFull,int nmol,rvec x[],matrix box,real rcut,
 385             gmx_bool bYaw,real q[],gmx_bool bLJ)
 386 {
 387   int  i,j,im,jm,natmol,ik,jk,m,ninter;
 388   rvec dx,f,ftot,dvir,vir,pres,xcmi,xcmj,*force;
 389   real dx6,dx2,dx1,fscal,c6,c12,vcoul,v12,v6,vctot,v12tot,v6tot;
 390   t_pbc pbc;
 391
 392   set_pbc(&pbc,box);
 393   fprintf(fp,"%3s   -  %3s: %6s %6s %6s  %6s %8s %8s %8s\n",
 394           "ai","aj","dx","dy","dz","|d|","virx","viry","virz");
 395   clear_rvec(ftot);
 396   clear_rvec(vir);
 397   ninter = 0;
 398   vctot  = 0;
 399   v12tot = 0;
 400   v6tot  = 0;
 401   natmol = bYaw ? 5 : 3;
 402   snew(force,nmol*natmol);
 403
 404   for(i=0; (i<nmol); i++) {
 405     im = natmol*i;
 406     /* Center of geometry */
 407     clear_rvec(xcmi);
 408     for(ik=0; (ik<natmol); ik++)
 409       rvec_inc(xcmi,x[im+ik]);
 410     for(m=0; (m<DIM); m++)
 411       xcmi[m] /= natmol;
 412
 413     for(j=i+1; (j<nmol); j++) {
 414       jm = natmol*j;
 415       /* Center of geometry */
 416       clear_rvec(xcmj);
 417       for(jk=0; (jk<natmol); jk++)
 418         rvec_inc(xcmj,x[jm+jk]);
 419       for(m=0; (m<DIM); m++)
 420         xcmj[m] /= natmol;
 421
 422       /* First check COM-COM distance */
 423       pbc_dx(&pbc,xcmi,xcmj,dx);
 424       if (norm(dx) < rcut) {
 425         ninter++;
 426         /* Neirest neighbour molecules! */
 427         clear_rvec(dvir);
 428         for(ik=0; (ik<natmol); ik++) {
 429           for(jk=0; (jk<natmol); jk++) {
 430             pbc_dx(&pbc,x[im+ik],x[jm+jk],dx);
 431             dx2    = iprod(dx,dx);
 432             dx1    = sqrt(dx2);
 433             vcoul  = q[ik]*q[jk]*ONE_4PI_EPS0/dx1;
 434             vctot += vcoul;
 435
 436             if (bLJ) {
 437               if (bYaw) {
 438                 c6  = yaw_lj[ik][2*jk];
 439                 c12 = yaw_lj[ik][2*jk+1];
 440               }
 441               else {
 442                 c6  = spc_lj[ik][2*jk];
 443                 c12 = spc_lj[ik][2*jk+1];
 444               }
 445               dx6    = dx2*dx2*dx2;
 446               v6     = c6/dx6;
 447               v12    = c12/(dx6*dx6);
 448               v6tot -= v6;
 449               v12tot+= v12;
 450               fscal  = (vcoul+12*v12-6*v6)/dx2;
 451             }
 452             else
 453               fscal  = vcoul/dx2;
 454             for(m=0; (m<DIM); m++) {
 455               f[m]     = dx[m]*fscal;
 456               dvir[m] -= 0.5*dx[m]*f[m];
 457             }
 458             rvec_inc(force[ik+im],f);
 459             rvec_dec(force[jk+jm],f);
 460             /*if (bFull)
 461               fprintf(fp,"%3s%4d-%3s%4d: %6.3f %6.3f %6.3f %6.3f"
 462                       " %8.3f %8.3f %8.3f\n",
 463                       watname[ik],im+ik,watname[jk],jm+jk,
 464                       dx[XX],dx[YY],dx[ZZ],norm(dx),
 465                       dvir[XX],dvir[YY],dvir[ZZ]);*/
 466           }
 467         }
 468         if (bFull)
 469           fprintf(fp,"%3s%4d-%3s%4d: "
 470                   " %8.3f %8.3f %8.3f\n",
 471                   "SOL",i,"SOL",j,dvir[XX],dvir[YY],dvir[ZZ]);
 472         rvec_inc(vir,dvir);
 473       }
 474     }
 475   }
 476   fprintf(fp,"There were %d interactions between the %d molecules (%.2f %%)\n",
 477           ninter,nmol,(real)ninter/(0.5*nmol*(nmol-1)));
 478   fprintf(fp,"Vcoul: %10.4e  V12: %10.4e  V6: %10.4e  Vtot: %10.4e (kJ/mol)\n",
 479           vctot/nmol,v12tot/nmol,v6tot/nmol,(vctot+v12tot+v6tot)/nmol);
 480   pr_rvec(fp,0,"vir ",vir,DIM,TRUE);
 481
 482   for(m=0; (m<DIM); m++)
 483     pres[m] = -2*PRESFAC/(det(box))*vir[m];
 484   pr_rvec(fp,0,"pres",pres,DIM,TRUE);
 485   pr_rvecs(fp,0,"force",force,natmol*nmol);
 486   sfree(force);
 487 }
 488
 489
 490
 491 int main(int argc,char *argv[])
 492 {
 493   static char *desc[] = {
 494     "[TT]mkice[tt] generates an ice crystal in the Ih crystal form which is the",
 495     "most stable form. The rectangular unitcell contains eight molecules",
 496     "and all oxygens are tetrahedrally coordinated.[PAR]",
 497     "If an input file is given it is interpreted as a series of oxygen",
 498     "coordinates the distance between which can be scaled by the odist flag.",
 499     "The program then adds hydrogens to the oxygens in random orientation",
 500     "but with proper OH distances and HOH angle. This feature allows to",
 501     "build water clusters based on oxygen coordinates only."
 502   };
 503   static int nx=1,ny=1,nz=1;
 504   static gmx_bool bYaw=FALSE,bLJ=TRUE,bFull=TRUE,bSeries=FALSE;
 505   static gmx_bool bOrdered=TRUE,bDiamond=FALSE,bPBC=TRUE;
 506   static real rcut=0.3,odist=0.274,hdist=0.09572;
 507   t_pargs pa[] = {
 508     { "-nx",    FALSE, etINT,  {&nx}, "nx" },
 509     { "-ny",    FALSE, etINT,  {&ny}, "ny" },
 510     { "-nz",    FALSE, etINT,  {&nz}, "nz" },
 511     { "-yaw",   FALSE, etBOOL, {&bYaw},
 512       "Generate virtual sites and shell positions" },
 513     { "-lj",    FALSE, etBOOL, {&bLJ},
 514       "Use LJ as well as coulomb for virial calculation" },
 515     { "-rcut",  FALSE,etREAL,  {&rcut},"Cut-off for virial calculations" },
 516     { "-full",  FALSE,etBOOL,  {&bFull},"Full virial output" },
 517     { "-odist", FALSE, etREAL, {&odist}, "Distance (nm) between oxygens" },
 518     { "-hdist", FALSE, etREAL, {&hdist}, "Bondlength (nm) for OH bond" },
 519     { "-diamond",FALSE,etBOOL, {&bDiamond}, "Make a diamond instead" },
 520     { "-pbc",   FALSE, etBOOL, {&bPBC},  "Make a periodic diamond" },
 521     { "-order", FALSE,etBOOL,  {&bOrdered}, "Make a proton-ordered ice lattice" },
 522     { "-series",FALSE, etBOOL, {&bSeries},
 523       "Do a series of virial calculations with different cut-off (from 0.3 up till the specified one)" }
 524   };
 525   t_filenm fnm[] = {
 526     { efSTO, "-p", "ice", ffWRITE },
 527     { efSTO, "-c", NULL,  ffOPTRD },
 528     { efDAT, "-f", NULL,  ffOPTRD },
 529     { efTRN, "-o", "ice", ffOPTWR }
 530   };
 531 #define NFILE asize(fnm)
 532
 533   FILE      *fp;
 534   char      *fn,quote[256];
 535   int       i,j,k,n,nmax,m,natom,natmol;
 536   t_atoms   *pdba;
 537   t_atoms   atoms;
 538   t_symtab  symtab;
 539   rvec      box,tmp,*xx;
 540   matrix    boxje;
 541
 542   CopyRight(stdout,argv[0]);
 543   parse_common_args(&argc,argv,0,NFILE,fnm,asize(pa),pa,asize(desc),
 544                     desc,0,NULL);
 545   if (debug) {
 546     fprintf(debug,"nx  = %3d, ny  = %3d,  nz   = %3d\n",nx,ny,nz);
 547     fprintf(debug,"YAW = %3s, LJ  = %3s,  rcut = %g\n",yesno_names[bYaw],
 548             yesno_names[bLJ],rcut);
 549   }
 550
 551   if (bYaw)
 552     natmol = 5;
 553   else if (bDiamond)
 554     natmol = 1;
 555   else
 556     natmol = 3;
 557
 558   if (opt2bSet("-f",NFILE,fnm)) {
 559     natom = read_rel_coords(opt2fn("-f",NFILE,fnm),&xx,natmol);
 560     nmax  = natom;
 561   }
 562   else {
 563     natom = natmol*8;
 564     nmax = natom*nx*ny*nz;
 565     snew(xx,nmax);
 566   }
 567   snew(pdba,1);
 568   init_t_atoms(pdba,nmax,TRUE);
 569   pdba->nr = nmax;
 570   open_symtab(&symtab);
 571   for(n=0; (n<nmax); n++) {
 572     pdba->pdbinfo[n].type   = epdbATOM;
 573     pdba->pdbinfo[n].atomnr = 1+n;
 574     pdba->atom[n].resnr     = 1+(n/natmol);
 575     pdba->atomname[n] = put_symtab(&symtab,
 576                                    bDiamond ? diamname[(n % natmol)] : watname[(n % natmol)]);
 577     if (bDiamond)
 578       pdba->resname[n] = put_symtab(&symtab,"DIA");
 579     else
 580       pdba->resname[n] = put_symtab(&symtab,"SOL");
 581     sprintf(pdba->pdbinfo[n].pdbresnr,"%d",n);
 582     pdba->atom[n].chain = ' ';
 583   }
 584
 585   /* Generate the unit cell */
 586   if (bDiamond)
 587     unitcell_d(xx,box,odist);
 588   else if (opt2bSet("-f",NFILE,fnm)) {
 589     random_h_coords(natmol,natom/natmol,xx,box,bYaw,odist,hdist);
 590   }
 591   else
 592     unitcell(xx,box,bYaw,odist,hdist);
 593   if (debug) {
 594     clear_mat(boxje);
 595     boxje[XX][XX] = box[XX];
 596     boxje[YY][YY] = box[YY];
 597     boxje[ZZ][ZZ] = box[ZZ];
 598   }
 599   n=0;
 600   for(i=0; (i<nx); i++) {
 601     tmp[XX] = i*box[XX];
 602     for(j=0; (j<ny); j++) {
 603       tmp[YY] = j*box[YY];
 604       for(k=0; (k<nz); k++) {
 605         tmp[ZZ] = k*box[ZZ];
 606         for(m=0; (m<natom); m++,n++) {
 607           if ((!bOrdered && ((m % natmol) == 0)) || bOrdered)
 608             rvec_add(xx[n % natom],tmp,xx[n]);
 609           else
 610             ;
 611         }
 612       }
 613     }
 614   }
 615
 616   clear_mat(boxje);
 617   boxje[XX][XX] = box[XX]*nx;
 618   boxje[YY][YY] = box[YY]*ny;
 619   boxje[ZZ][ZZ] = box[ZZ]*nz;
 620
 621   printf("Crystal:   %10.5f  %10.5f  %10.5f\n",
 622          nx*box[XX],ny*box[YY],nz*box[ZZ]);
 623
 624   if (debug && !bDiamond) {
 625     if (bSeries)
 626       for(i=3; (i<=10*rcut); i++) {
 627         fprintf(debug,"This is with rcut = %g\n",i*0.1);
 628         virial(debug,bFull,nmax/natmol,xx,boxje,
 629                0.1*i,bYaw,bYaw ? qyaw : qspc,bLJ);
 630       }
 631     else
 632       virial(debug,bFull,nmax/natmol,xx,boxje,
 633              rcut,bYaw,bYaw ? qyaw : qspc,bLJ);
 634   }
 635
 636   if (bDiamond)
 637     mk_diamond(pdba,xx,odist,&symtab,bPBC,boxje);
 638
 639   fn = ftp2fn(efSTO,NFILE,fnm);
 640   if (fn2ftp(fn) == efPDB) {
 641     fp = ffopen(fn,"w");
 642     if (bDiamond)
 643       fprintf(fp,"HEADER    This is a *diamond*\n");
 644     else
 645       fprintf(fp,"HEADER    A beautiful Ice Ih crystal\n");
 646     fprintf(fp,"REMARK    Generated by mkice with the following options:\n"
 647             "REMARK    nx = %d, ny = %d, nz = %d, odist = %g, hdist = %g\n",
 648             nx,ny,nz,odist,hdist);
 649         bromacs(quote,255);
 650     write_pdbfile(fp,quote,pdba,xx,boxje,' ',-1);
 651     ffclose(fp);
 652   }
 653   else {
 654     bromacs(quote,255);
 655     write_sto_conf(fn,quote,pdba,xx,NULL,boxje);
 656   }
 657
 658   if (ftp2bSet(efTRN,NFILE,fnm))
 659     write_trn(ftp2fn(efTRN,NFILE,fnm),0,0,0,boxje,nmax,xx,NULL,NULL);
 660
 661   return 0;
 662 }
 663