src/local/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 "confio.h"
  48
  49 #define HDIST 0.1
  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)
  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 int main(int argc,char *argv[])
 391 {
 392   static char *desc[] = {
 393     "mkice generates an ice crystal in the Ih crystal form which is the",
 394     "most stable form. The rectangular unitcell contains eight molecules",
 395     "and all oxygens are tetrahedrally coordinated."
 396   };
 397   static int nx=1,ny=1,nz=1;
 398   static bool bYaw=FALSE,bLJ=TRUE,bFull=TRUE,bSeries=FALSE,bDiamond=FALSE;
 399   static real rcut=0.3,odist=0.274;
 400   t_pargs pa[] = {
 401     { "-nx",    FALSE, etINT,  {&nx}, "nx" },
 402     { "-ny",    FALSE, etINT,  {&ny}, "ny" },
 403     { "-nz",    FALSE, etINT,  {&nz}, "nz" },
 404     { "-yaw",   FALSE, etBOOL, {&bYaw},
 405       "Generate dummies and shell positions" },
 406     { "-lj",    FALSE, etBOOL, {&bLJ},
 407       "Use LJ as well as coulomb for virial calculation" },
 408     { "-rcut",  FALSE,etREAL,  {&rcut},"Cut-off for virial calculations" },
 409     { "-full",  FALSE,etBOOL,  {&bFull},"Full virial output" },
 410     { "-odist", FALSE, etREAL, {&odist}, "Distance between oxygens" },
 411     { "-diamond",FALSE,etBOOL, {&bDiamond}, "Make a diamond instead" },
 412     { "-series",FALSE, etBOOL, {&bSeries},
 413       "Do a series of virial calculations with different cut-off (from 0.3 up till the specified one)" }
 414   };
 415   t_filenm fnm[] = {
 416     { efPDB, "-p", "ice", ffWRITE },
 417     { efSTO, "-c", NULL,  ffOPTRD },
 418     { efTRN, "-o", "ice", ffOPTWR }
 419   };
 420 #define NFILE asize(fnm)
 421
 422   FILE      *fp;
 423   int       i,j,k,n,nmax,m,natom,natmol;
 424   t_atoms   *pdba;
 425   t_atoms   atoms;
 426   t_symtab  symtab;
 427   rvec      box,tmp,*xx;
 428   matrix    boxje;
 429
 430   CopyRight(stdout,argv[0]);
 431   parse_common_args(&argc,argv,0,FALSE,NFILE,fnm,asize(pa),pa,asize(desc),
 432                     desc,0,NULL);
 433   if (debug) {
 434     fprintf(debug,"nx  = %3d, ny  = %3d,  nz   = %3d\n",nx,ny,nz);
 435     fprintf(debug,"YAW = %3s, LJ  = %3s,  rcut = %g\n",yesno_names[bYaw],
 436             yesno_names[bLJ],rcut);
 437   }
 438
 439   if (bYaw)
 440     natmol = 5;
 441   else if (bDiamond)
 442     natmol = 1;
 443   else
 444     natmol = 3;
 445   natom = natmol*8;
 446   nmax = natom*nx*ny*nz;
 447   snew(xx,nmax);
 448   snew(pdba,1);
 449   init_t_atoms(pdba,nmax,TRUE);
 450   pdba->nr = nmax;
 451   open_symtab(&symtab);
 452   for(n=0; (n<nmax); n++) {
 453     pdba->pdbinfo[n].type   = epdbATOM;
 454     pdba->pdbinfo[n].atomnr = n;
 455     pdba->atom[n].resnr     = n/natmol;
 456     pdba->atomname[n] = put_symtab(&symtab,
 457                                    bDiamond ? diamname[(n % natmol)] : watname[(n % natmol)]);
 458     if (bDiamond)
 459       pdba->resname[n] = put_symtab(&symtab,"DIA");
 460     else
 461       pdba->resname[n] = put_symtab(&symtab,"SOL");
 462     sprintf(pdba->pdbinfo[n].pdbresnr,"%d",n);
 463     pdba->atom[n].chain = ' ';
 464   }
 465
 466   /* Generate the unit cell */
 467   if (bDiamond)
 468     unitcell_d(xx,box,odist);
 469   else
 470     unitcell(xx,box,bYaw,odist);
 471   if (debug) {
 472     clear_mat(boxje);
 473     boxje[XX][XX] = box[XX];
 474     boxje[YY][YY] = box[YY];
 475     boxje[ZZ][ZZ] = box[ZZ];
 476   }
 477   n=0;
 478   for(i=0; (i<nx); i++) {
 479     tmp[XX] = i*box[XX];
 480     for(j=0; (j<ny); j++) {
 481       tmp[YY] = j*box[YY];
 482       for(k=0; (k<nz); k++) {
 483         tmp[ZZ] = k*box[ZZ];
 484         for(m=0; (m<natom); m++,n++) {
 485           rvec_add(xx[n % natom],tmp,xx[n]);
 486         }
 487       }
 488     }
 489   }
 490
 491   clear_mat(boxje);
 492   boxje[XX][XX] = box[XX]*nx;
 493   boxje[YY][YY] = box[YY]*ny;
 494   boxje[ZZ][ZZ] = box[ZZ]*nz;
 495
 496   printf("Crystal:   %10.5f  %10.5f  %10.5f\n",
 497          nx*box[XX],ny*box[YY],nz*box[ZZ]);
 498
 499   if (debug && !bDiamond) {
 500     if (bSeries)
 501       for(i=3; (i<=10*rcut); i++) {
 502         fprintf(debug,"This is with rcut = %g\n",i*0.1);
 503         virial(debug,bFull,nmax/natmol,xx,boxje,
 504                0.1*i,bYaw,bYaw ? qyaw : qspc,bLJ);
 505       }
 506     else
 507       virial(debug,bFull,nmax/natmol,xx,boxje,
 508              rcut,bYaw,bYaw ? qyaw : qspc,bLJ);
 509   }
 510
 511   if (bDiamond)
 512     mk_diamond(pdba,xx,odist,&symtab);
 513
 514   fp = ftp2FILE(efPDB,NFILE,fnm,"w");
 515   if (bDiamond)
 516     fprintf(fp,"HEADER    This is a *diamond*\n");
 517   else
 518     fprintf(fp,"HEADER    A beautifil Ice Ih crystal\n");
 519   fprintf(fp,"REMARK    Generated by mkice with the following options:\n"
 520           "REMARK    nx = %d, ny = %d, nz = %d, odist = %g\n",nx,ny,nz,odist);
 521   write_pdbfile(fp,bromacs(),pdba,xx,boxje,' ',FALSE);
 522   fclose(fp);
 523
 524   if (ftp2bSet(efTRN,NFILE,fnm))
 525     write_trn(ftp2fn(efTRN,NFILE,fnm),0,0,0,boxje,nmax,xx,NULL,NULL);
 526
 527   return 0;
 528 }
 529