src/contrib/mkice.c

   1 /*
   2  * $Id$
   3  *
   4  *       This source code is part of
   5  *
   6  *        G   R   O   M   A   C   S
   7  *
   8  * GROningen MAchine for Chemical Simulations
   9  *
  10  *               VERSION 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * Great Red Oystrich Makes All Chemists Sane
  28  */
  29 static char *SRCID_mkice_c = "$Id$";
  30
  31 #include <stdio.h>
  32 #include <math.h>
  33 #include "typedefs.h"
  34 #include "statutil.h"
  35 #include "copyrite.h"
  36 #include "fatal.h"
  37 #include "pdbio.h"
  38 #include "macros.h"
  39 #include "smalloc.h"
  40 #include "vec.h"
  41 #include "pbc.h"
  42 #include "physics.h"
  43 #include "names.h"
  44 #include "txtdump.h"
  45 #include "trnio.h"
  46 #include "symtab.h"
  47 #include "atomprop.h"
  48 #include "confio.h"
  49
  50 #define TET   109.47
  51 #define DCONS 0.117265878
  52
  53 typedef struct {
  54   int n,aa[4];
  55 } t_bbb;
  56
  57 static char *watname[] = { "OW ", "HW1", "HW2", "DW", "SW" };
  58 static char *diamname[] = { "C1", "H2", "H3", "H4", "H5", "H2", "H3", "H4", "H5" };
  59 static real qspc[]     = { -0.82, 0.41, 0.41 };
  60 static real qyaw[]     = { 1.24588, 0.62134, 0.62134, 0.0, -2.48856 };
  61 static real spc_lj[3][6] = {
  62   { 2.6171e-3, 2.6331e-6, 0, 0, 0, 0 },
  63   { 0,      0,      0, 0, 0, 0 },
  64   { 0,      0,      0, 0, 0, 0 }
  65 };
  66 #define CHH 2e-9
  67 #define CHS 2e-9
  68 #define COS 2e-6
  69 static real yaw_lj[5][10] = {
  70   { 0, 0,   0, 0,   0,   0, 0, 0, 0,      COS },
  71   { 0, 0,   0, CHH, 0, CHH, 0, 0, 0,      CHS },
  72   { 0, 0,   0, CHH, 0, CHH, 0, 0, 0,      CHS },
  73   { 0, 0,   0, 0,   0,   0, 0, 0, 0,      0   },
  74   { 0, COS, 0, CHS, 0, CHS, 0, 0, 2.6e-3, 0   }
  75 };
  76
  77 void unitcell(rvec x[],rvec box,bool bYaw,real odist,real hdist)
  78 {
  79 #define cx  0.81649658
  80 #define cy  0.47140452
  81 #define cy2 0.94280904
  82 #define cz  0.33333333
  83
  84   rvec xx[24] = {
  85     { 0,   0,         0 }, /* O1 */
  86     { 0,   0,         1 }, /* H relative to Oxygen */
  87     { cx, cy,       -cz },
  88     { cx, cy,       -cz }, /* O2 */
  89     { 0, 0,       -1    }, /* H relative to Oxygen */
  90     { cx,-cy,       +cz },
  91     { cx, cy+cy2,     0 }, /* O3 */
  92     { -cx, cy,    -cz   }, /* H relative to Oxygen */
  93     { 0,   -cy2,    -cz },
  94     { 0,  2*cy+cy2, -cz }, /* O4 */
  95     {-cx,-cy,       +cz }, /* H relative to Oxygen */
  96     { 0 , cy2,      +cz },
  97     { 0,   0,         1 }, /* O5 */
  98     {-cx, cy,       +cz }, /* H relative to Oxygen */
  99     { 0 , -cy2,     +cz },
 100     { cx, cy,      1+cz }, /* O6 */
 101     { -cx, -cy,     -cz }, /* H relative to Oxygen */
 102     { 0,   cy2,     -cz },
 103     { cx, cy+cy2,     1 }, /* O7 */
 104     { 0,  0,       -1   }, /* H relative to Oxygen */
 105     { cx, cy,       +cz },
 106     { 0,  2*cy+cy2,1+cz }, /* O8 */
 107     { 0,   0,         1 }, /* H relative to Oxygen */
 108     { cx,   -cy,    -cz }
 109   };
 110   int  i,iin,iout,j,m;
 111   rvec tmp,t2,dip;
 112
 113   clear_rvec(dip);
 114   for(i=0; (i<8); i++) {
 115     iin = 3*i;
 116     if (bYaw)
 117       iout = 5*i;
 118     else
 119       iout = iin;
 120     svmul(odist,xx[iin],x[iout]);
 121     svmul(-0.82,x[iout],t2);
 122     rvec_inc(dip,t2);
 123     for(j=1; (j<=2); j++) {
 124       svmul(hdist,xx[iin+j],tmp);
 125       rvec_add(x[iout],tmp,x[iout+j]);
 126       svmul(0.41,x[iout+j],t2);
 127       rvec_inc(dip,t2);
 128     }
 129     if (bYaw)
 130       for(m=0; (m<DIM); m++)
 131         x[iout+3][m] = x[iout+4][m] =
 132           (1-2*DCONS)*x[iout][m]+DCONS*(x[iout+1][m]+x[iout+2][m]);
 133   }
 134
 135   box[XX] = 2*cx;
 136   box[YY] = 2*(cy2+cy);
 137   box[ZZ] = 2*(1+cz);
 138   for(i=0; (i<DIM); i++)
 139     box[i] *= odist;
 140
 141   printf("Unitcell:  %10.5f  %10.5f  %10.5f\n",box[XX],box[YY],box[ZZ]);
 142   printf("Dipole:    %10.5f  %10.5f  %10.5f (e nm)\n",dip[XX],dip[YY],dip[ZZ]);
 143 }
 144
 145 void unitcell_d(rvec x[],rvec box,real odist)
 146 {
 147   rvec cc[8] = {
 148     { 0,   0,         0 }, /* C1 */
 149     { cx, cy,       -cz }, /* C2 */
 150     { cx, cy+cy2,     0 }, /* C3 */
 151     { 0,  2*cy+cy2, -cz }, /* C4 */
 152     { 0,   0,         1 }, /* C5 */
 153     { cx, cy,      1+cz }, /* C6 */
 154     { cx, cy+cy2,     1 }, /* C7 */
 155     { 0,  2*cy+cy2,1+cz }, /* C8 */
 156   };
 157   rvec hh[4] = {
 158     { 0,   0,         1  }, /* H relative to C */
 159     { cx,  cy,       -cz },
 160     { cx, -cy,       -cz },
 161     {-cy2,  0,       -cz }
 162   };
 163   int  i,iin,iout,j,m;
 164   rvec tmp,t2,dip;
 165
 166   clear_rvec(dip);
 167   for(i=0; (i<8); i++) {
 168     iin  = i;
 169     iout = i;
 170     svmul(odist,cc[iin],x[iout]);
 171   }
 172   box[XX] = 2*cx;
 173   box[YY] = 2*(cy2+cy);
 174   box[ZZ] = 2*(1+cz);
 175   for(i=0; (i<DIM); i++)
 176     box[i] *= odist;
 177
 178   printf("Unitcell:  %10.5f  %10.5f  %10.5f\n",box[XX],box[YY],box[ZZ]);
 179 }
 180
 181 static t_bbb *mk_bonds(int natoms,rvec x[],real odist)
 182 {
 183   real  od2 = odist*odist+1e-5;
 184   t_bbb *bbb;
 185   int   i,j;
 186   rvec  dx;
 187
 188   snew(bbb,natoms);
 189   for(i=0; (i<natoms); i++) {
 190     for(j=i+1; (j<natoms); j++) {
 191       rvec_sub(x[i],x[j],dx);
 192       if (iprod(dx,dx) <= od2) {
 193         bbb[i].aa[bbb[i].n++] = j;
 194         bbb[j].aa[bbb[j].n++] = i;
 195       }
 196     }
 197   }
 198   if (debug)
 199 #define PRB(nn) (bbb[(i)].n >= (nn)) ? bbb[i].aa[nn-1] : -1
 200     for(i=0; (i<natoms); i++)
 201       fprintf(debug,"bonds from %3d:  %d %d %d %d\n",
 202               i,PRB(1),PRB(2),PRB(3),PRB(4));
 203 #undef PRB
 204   return bbb;
 205 }
 206
 207 static void mk_diamond(t_atoms *a,rvec x[],real odist,t_symtab *symtab)
 208 {
 209
 210   int   i,ib,j,k,l,m,nrm=0;
 211   t_bbb *bbb;
 212   bool  *bRemove;
 213   rvec  dx;
 214
 215   do {
 216     nrm = 0;
 217     bbb = mk_bonds(a->nr,x,odist);
 218
 219     for(i=0; (i<a->nr); i++) {
 220       if (bbb[i].n < 2) {
 221         for(k=0; (k<bbb[i].n); k++) {
 222           ib = bbb[i].aa[k];
 223           for(j=0; (j<bbb[ib].n); j++)
 224             if (bbb[ib].aa[j] == i)
 225               break;
 226           if (j == bbb[ib].n)
 227             fatal_error(0,"Bond inconsistency (%d not in list of %d)!\n",i,ib);
 228           for( ; (j<bbb[ib].n-1); j++)
 229             bbb[ib].aa[j] = bbb[ib].aa[j+1];
 230           bbb[ib].n--;
 231           nrm++;
 232         }
 233         bbb[i].n = 0;
 234       }
 235     }
 236
 237     for(i=j=0; (i<a->nr); i++) {
 238       if (bbb[i].n >= 2) {
 239         copy_rvec(x[i],x[j]);
 240         j++;
 241       }
 242     }
 243     fprintf(stderr,"Kicking out %d carbon atoms (out of %d)\n",
 244             a->nr-j,a->nr);
 245     a->nr = j;
 246     sfree(bbb);
 247   } while (nrm > 0);
 248   /* Rename atoms */
 249   bbb = mk_bonds(a->nr,x,odist);
 250   for(i=0; (i<a->nr); i++) {
 251     switch (bbb[i].n) {
 252     case 4:
 253       a->atomname[i] = put_symtab(symtab,"C");
 254       break;
 255     case 3:
 256       a->atomname[i] = put_symtab(symtab,"CH1");
 257       break;
 258     case 2:
 259       a->atomname[i] = put_symtab(symtab,"CH2");
 260       break;
 261     default:
 262       fatal_error(0,"This atom (%d) has %d bonds only",i,bbb[i].n);
 263     }
 264   }
 265   sfree(bbb);
 266 }
 267
 268 static real calc_ener(real c6,real c12,rvec dx,tensor vir)
 269 {
 270   real r,e,f;
 271   int  m,n;
 272
 273   r  = norm(dx);
 274   e  = c12*pow(r,-12.0) - c6*pow(r,-6.0);
 275   f  = 12*c12*pow(r,-14.0) - 6*c6*pow(r,-8.0);
 276   for(m=0; (m<DIM); m++)
 277     for(n=0; (n<DIM); n++)
 278       vir[m][n] -= 0.5*f*dx[m]*dx[n];
 279
 280   return e;
 281 }
 282
 283 void virial(FILE *fp,bool bFull,int nmol,rvec x[],matrix box,real rcut,
 284             bool bYaw,real q[],bool bLJ)
 285 {
 286   int  i,j,im,jm,natmol,ik,jk,m,ninter;
 287   rvec dx,f,ftot,dvir,vir,pres,xcmi,xcmj,*force;
 288   real dx6,dx2,dx1,fscal,c6,c12,vcoul,v12,v6,vctot,v12tot,v6tot;
 289
 290   init_pbc(box,FALSE);
 291   /* Initiate the periodic boundary conditions. Set bTruncOct to
 292    * TRUE when using a truncated octahedron box.
 293    */
 294   fprintf(fp,"%3s   -  %3s: %6s %6s %6s  %6s %8s %8s %8s\n",
 295           "ai","aj","dx","dy","dz","|d|","virx","viry","virz");
 296   clear_rvec(ftot);
 297   clear_rvec(vir);
 298   ninter = 0;
 299   vctot  = 0;
 300   v12tot = 0;
 301   v6tot  = 0;
 302   natmol = bYaw ? 5 : 3;
 303   snew(force,nmol*natmol);
 304
 305   for(i=0; (i<nmol); i++) {
 306     im = natmol*i;
 307     /* Center of geometry */
 308     clear_rvec(xcmi);
 309     for(ik=0; (ik<natmol); ik++)
 310       rvec_inc(xcmi,x[im+ik]);
 311     for(m=0; (m<DIM); m++)
 312       xcmi[m] /= natmol;
 313
 314     for(j=i+1; (j<nmol); j++) {
 315       jm = natmol*j;
 316       /* Center of geometry */
 317       clear_rvec(xcmj);
 318       for(jk=0; (jk<natmol); jk++)
 319         rvec_inc(xcmj,x[jm+jk]);
 320       for(m=0; (m<DIM); m++)
 321         xcmj[m] /= natmol;
 322
 323       /* First check COM-COM distance */
 324       pbc_dx(xcmi,xcmj,dx);
 325       if (norm(dx) < rcut) {
 326         ninter++;
 327         /* Neirest neighbour molecules! */
 328         clear_rvec(dvir);
 329         for(ik=0; (ik<natmol); ik++) {
 330           for(jk=0; (jk<natmol); jk++) {
 331             pbc_dx(x[im+ik],x[jm+jk],dx);
 332             dx2    = iprod(dx,dx);
 333             dx1    = sqrt(dx2);
 334             vcoul  = q[ik]*q[jk]*ONE_4PI_EPS0/dx1;
 335             vctot += vcoul;
 336
 337             if (bLJ) {
 338               if (bYaw) {
 339                 c6  = yaw_lj[ik][2*jk];
 340                 c12 = yaw_lj[ik][2*jk+1];
 341               }
 342               else {
 343                 c6  = spc_lj[ik][2*jk];
 344                 c12 = spc_lj[ik][2*jk+1];
 345               }
 346               dx6    = dx2*dx2*dx2;
 347               v6     = c6/dx6;
 348               v12    = c12/(dx6*dx6);
 349               v6tot -= v6;
 350               v12tot+= v12;
 351               fscal  = (vcoul+12*v12-6*v6)/dx2;
 352             }
 353             else
 354               fscal  = vcoul/dx2;
 355             for(m=0; (m<DIM); m++) {
 356               f[m]     = dx[m]*fscal;
 357               dvir[m] -= 0.5*dx[m]*f[m];
 358             }
 359             rvec_inc(force[ik+im],f);
 360             rvec_dec(force[jk+jm],f);
 361             /*if (bFull)
 362               fprintf(fp,"%3s%4d-%3s%4d: %6.3f %6.3f %6.3f %6.3f"
 363                       " %8.3f %8.3f %8.3f\n",
 364                       watname[ik],im+ik,watname[jk],jm+jk,
 365                       dx[XX],dx[YY],dx[ZZ],norm(dx),
 366                       dvir[XX],dvir[YY],dvir[ZZ]);*/
 367           }
 368         }
 369         if (bFull)
 370           fprintf(fp,"%3s%4d-%3s%4d: "
 371                   " %8.3f %8.3f %8.3f\n",
 372                   "SOL",i,"SOL",j,dvir[XX],dvir[YY],dvir[ZZ]);
 373         rvec_inc(vir,dvir);
 374       }
 375     }
 376   }
 377   fprintf(fp,"There were %d interactions between the %d molecules (%.2f %%)\n",
 378           ninter,nmol,(real)ninter/(0.5*nmol*(nmol-1)));
 379   fprintf(fp,"Vcoul: %10.4e  V12: %10.4e  V6: %10.4e  Vtot: %10.4e (kJ/mol)\n",
 380           vctot/nmol,v12tot/nmol,v6tot/nmol,(vctot+v12tot+v6tot)/nmol);
 381   pr_rvec(fp,0,"vir ",vir,DIM);
 382
 383   for(m=0; (m<DIM); m++)
 384     pres[m] = -2*PRESFAC/(det(box))*vir[m];
 385   pr_rvec(fp,0,"pres",pres,DIM);
 386   pr_rvecs(fp,0,"force",force,natmol*nmol);
 387   sfree(force);
 388 }
 389
 390 static real calc_mass(t_atoms *atoms,bool bGetMass)
 391 {
 392   real tmass;
 393   int i;
 394
 395   tmass = 0;
 396   for(i=0; (i<atoms->nr); i++) {
 397     if (bGetMass)
 398       atoms->atom[i].m = get_mass(*atoms->resname[atoms->atom[i].resnr],
 399                                   *atoms->atomname[i]);
 400     tmass += atoms->atom[i].m;
 401   }
 402
 403   return tmass;
 404 }
 405
 406 static real density(t_atoms *atoms,matrix box)
 407 {
 408   real tm;
 409   real vol;
 410
 411   vol = det(box);
 412   tm  = calc_mass(atoms,TRUE);
 413
 414   return (tm*AMU)/(vol*NANO*NANO*NANO);
 415 }
 416
 417 int main(int argc,char *argv[])
 418 {
 419   static char *desc[] = {
 420     "mkice generates an ice crystal in the Ih crystal form which is the",
 421     "most stable form. The rectangular unitcell contains eight molecules",
 422     "and all oxygens are tetrahedrally coordinated."
 423   };
 424   static int nx=1,ny=1,nz=1;
 425   static bool bYaw=FALSE,bLJ=TRUE,bFull=TRUE,bSeries=FALSE;
 426   static bool bOrdered=TRUE,bDiamond=FALSE;
 427   static real rcut=0.3,odist=0.274,hdist=0.09572;
 428   t_pargs pa[] = {
 429     { "-nx",    FALSE, etINT,  {&nx}, "nx" },
 430     { "-ny",    FALSE, etINT,  {&ny}, "ny" },
 431     { "-nz",    FALSE, etINT,  {&nz}, "nz" },
 432     { "-yaw",   FALSE, etBOOL, {&bYaw},
 433       "Generate dummies and shell positions" },
 434     { "-lj",    FALSE, etBOOL, {&bLJ},
 435       "Use LJ as well as coulomb for virial calculation" },
 436     { "-rcut",  FALSE,etREAL,  {&rcut},"Cut-off for virial calculations" },
 437     { "-full",  FALSE,etBOOL,  {&bFull},"Full virial output" },
 438     { "-odist", FALSE, etREAL, {&odist}, "Distance (nm) between oxygens" },
 439     { "-hdist", FALSE, etREAL, {&hdist}, "Bondlength (nm) for OH bond" },
 440     { "-diamond",FALSE,etBOOL, {&bDiamond}, "Make a diamond instead" },
 441     { "-order", FALSE,etBOOL,  {&bOrdered}, "Make a proton-ordered ice lattice" },
 442     { "-series",FALSE, etBOOL, {&bSeries},
 443       "Do a series of virial calculations with different cut-off (from 0.3 up till the specified one)" }
 444   };
 445   t_filenm fnm[] = {
 446     { efPDB, "-p", "ice", ffWRITE },
 447     { efSTO, "-c", NULL,  ffOPTRD },
 448     { efTRN, "-o", "ice", ffOPTWR }
 449   };
 450 #define NFILE asize(fnm)
 451
 452   FILE      *fp;
 453   int       i,j,k,n,nmax,m,natom,natmol;
 454   t_atoms   *pdba;
 455   t_atoms   atoms;
 456   t_symtab  symtab;
 457   rvec      box,tmp,*xx;
 458   matrix    boxje;
 459
 460   CopyRight(stdout,argv[0]);
 461   parse_common_args(&argc,argv,0,FALSE,NFILE,fnm,asize(pa),pa,asize(desc),
 462                     desc,0,NULL);
 463   if (debug) {
 464     fprintf(debug,"nx  = %3d, ny  = %3d,  nz   = %3d\n",nx,ny,nz);
 465     fprintf(debug,"YAW = %3s, LJ  = %3s,  rcut = %g\n",yesno_names[bYaw],
 466             yesno_names[bLJ],rcut);
 467   }
 468
 469   if (bYaw)
 470     natmol = 5;
 471   else if (bDiamond)
 472     natmol = 1;
 473   else
 474     natmol = 3;
 475   natom = natmol*8;
 476   nmax = natom*nx*ny*nz;
 477   snew(xx,nmax);
 478   snew(pdba,1);
 479   init_t_atoms(pdba,nmax,TRUE);
 480   pdba->nr = nmax;
 481   open_symtab(&symtab);
 482   for(n=0; (n<nmax); n++) {
 483     pdba->pdbinfo[n].type   = epdbATOM;
 484     pdba->pdbinfo[n].atomnr = n;
 485     pdba->atom[n].resnr     = n/natmol;
 486     pdba->atomname[n] = put_symtab(&symtab,
 487                                    bDiamond ? diamname[(n % natmol)] : watname[(n % natmol)]);
 488     if (bDiamond)
 489       pdba->resname[n] = put_symtab(&symtab,"DIA");
 490     else
 491       pdba->resname[n] = put_symtab(&symtab,"SOL");
 492     sprintf(pdba->pdbinfo[n].pdbresnr,"%d",n);
 493     pdba->atom[n].chain = ' ';
 494   }
 495
 496   /* Generate the unit cell */
 497   if (bDiamond)
 498     unitcell_d(xx,box,odist);
 499   else
 500     unitcell(xx,box,bYaw,odist,hdist);
 501   if (debug) {
 502     clear_mat(boxje);
 503     boxje[XX][XX] = box[XX];
 504     boxje[YY][YY] = box[YY];
 505     boxje[ZZ][ZZ] = box[ZZ];
 506   }
 507   n=0;
 508   for(i=0; (i<nx); i++) {
 509     tmp[XX] = i*box[XX];
 510     for(j=0; (j<ny); j++) {
 511       tmp[YY] = j*box[YY];
 512       for(k=0; (k<nz); k++) {
 513         tmp[ZZ] = k*box[ZZ];
 514         for(m=0; (m<natom); m++,n++) {
 515           if ((!bOrdered && ((m % natmol) == 0)) || bOrdered)
 516             rvec_add(xx[n % natom],tmp,xx[n]);
 517           else
 518             ;
 519         }
 520       }
 521     }
 522   }
 523
 524   clear_mat(boxje);
 525   boxje[XX][XX] = box[XX]*nx;
 526   boxje[YY][YY] = box[YY]*ny;
 527   boxje[ZZ][ZZ] = box[ZZ]*nz;
 528
 529   printf("Crystal:   %10.5f  %10.5f  %10.5f\n",
 530          nx*box[XX],ny*box[YY],nz*box[ZZ]);
 531
 532   if (debug && !bDiamond) {
 533     if (bSeries)
 534       for(i=3; (i<=10*rcut); i++) {
 535         fprintf(debug,"This is with rcut = %g\n",i*0.1);
 536         virial(debug,bFull,nmax/natmol,xx,boxje,
 537                0.1*i,bYaw,bYaw ? qyaw : qspc,bLJ);
 538       }
 539     else
 540       virial(debug,bFull,nmax/natmol,xx,boxje,
 541              rcut,bYaw,bYaw ? qyaw : qspc,bLJ);
 542   }
 543
 544   if (bDiamond)
 545     mk_diamond(pdba,xx,odist,&symtab);
 546
 547   fp = ftp2FILE(efPDB,NFILE,fnm,"w");
 548   if (bDiamond)
 549     fprintf(fp,"HEADER    This is a *diamond*\n");
 550   else
 551     fprintf(fp,"HEADER    A beautiful Ice Ih crystal\n");
 552   fprintf(fp,"REMARK    Generated by mkice with the following options:\n"
 553           "REMARK    nx = %d, ny = %d, nz = %d, odist = %g, hdist = %g\n",
 554           nx,ny,nz,odist,hdist);
 555   fprintf(fp,"REMARK    Density of this crystal is %g (g/l)\n",
 556           density(pdba,boxje));
 557   write_pdbfile(fp,bromacs(),pdba,xx,boxje,' ',-1);
 558   fclose(fp);
 559
 560   if (ftp2bSet(efTRN,NFILE,fnm))
 561     write_trn(ftp2fn(efTRN,NFILE,fnm),0,0,0,boxje,nmax,xx,NULL,NULL);
 562
 563   return 0;
 564 }
 565