src/contrib/mkice.c

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