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