src/programs/view/nmol.cpp

   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.2.0
  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-2013, 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  * Gyas ROwers Mature At Cryogenic Speed
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include <math.h>
  40 #include "sysstuff.h"
  41 #include <string.h>
  42 #include "smalloc.h"
  43 #include "macros.h"
  44 #include "xutil.h"
  45 #include "3dview.h"
  46 #include "gmx_fatal.h"
  47 #include "buttons.h"
  48 #include "manager.h"
  49 #include "nmol.h"
  50 #include "vec.h"
  51 #include "txtdump.h"
  52 #include "pbc.h"
  53
  54 #define MSIZE 4
  55
  56 static bool MWCallBack(t_x11 *x11, XEvent *event, Window w, void *data)
  57 {
  58     t_molwin *mw;
  59     Window    To;
  60     XEvent    letter;
  61
  62     mw                          = (t_molwin *)data;
  63     To                          = mw->wd.Parent;
  64     letter.type                 = ClientMessage;
  65     letter.xclient.display      = x11->disp;
  66     letter.xclient.window       = To;
  67     letter.xclient.message_type = 0;
  68     letter.xclient.format       = 32;
  69     switch (event->type)
  70     {
  71         case Expose:
  72             /* Do Not draw anything, but signal parent instead, he will
  73              * coordinate drawing.
  74              */
  75             letter.xclient.data.l[0] = IDDRAWMOL;
  76             letter.xclient.data.l[1] = Button1;
  77             XSendEvent(x11->disp, To, True, 0, &letter);
  78             break;
  79         case ButtonPress:
  80 #ifdef DEBUG
  81             printf("Molwindow: Buttonpress\n");
  82 #endif
  83             letter.xclient.data.l[0] = IDLABEL;
  84             letter.xclient.data.l[1] = (long)event->xbutton.button;
  85             letter.xclient.data.l[2] = event->xbutton.x;
  86             letter.xclient.data.l[3] = event->xbutton.y;
  87             XSendEvent(x11->disp, To, True, 0, &letter);
  88             break;
  89         case ConfigureNotify:
  90             mw->wd.width  = event->xconfigure.width;
  91             mw->wd.height = event->xconfigure.height;
  92             break;
  93         default:
  94             break;
  95     }
  96     return false;
  97 }
  98
  99 void set_def (t_molwin *mw, int ePBC, matrix box)
 100 {
 101     mw->bShowHydrogen = true;
 102     mw->bond_type     = eBFat;
 103     mw->ePBC          = ePBC;
 104     mw->boxtype       = esbRect;
 105     mw->realbox       = TRICLINIC(box) ? esbTri : esbRect;
 106 }
 107
 108 t_molwin *init_mw(t_x11 *x11, Window Parent,
 109                   int x, int y, int width, int height,
 110                   unsigned long fg, unsigned long bg,
 111                   int ePBC, matrix box)
 112 {
 113     t_molwin *mw;
 114
 115     snew(mw, 1);
 116     set_def(mw, ePBC, box);
 117
 118     InitWin(&mw->wd, x, y, width, height, 1, "Mol Window");
 119
 120     mw->wd.Parent = Parent;
 121     mw->wd.self   = XCreateSimpleWindow(x11->disp, Parent, x, y, width, height, 1, fg, bg);
 122     x11->RegisterCallback(x11, mw->wd.self, Parent, MWCallBack, mw);
 123     x11->SetInputMask(x11, mw->wd.self,
 124                       ExposureMask | StructureNotifyMask |
 125                       ButtonPressMask);
 126     return mw;
 127 }
 128
 129 void map_mw(t_x11 *x11, t_molwin *mw)
 130 {
 131     XMapWindow(x11->disp, mw->wd.self);
 132 }
 133
 134 bool toggle_hydrogen(t_x11 *x11, t_molwin *mw)
 135 {
 136     mw->bShowHydrogen = !mw->bShowHydrogen;
 137     ExposeWin(x11->disp, mw->wd.self);
 138
 139     return mw->bShowHydrogen;
 140 }
 141
 142 void set_bond_type(t_x11 *x11, t_molwin *mw, int bt)
 143 {
 144     if (bt != mw->bond_type)
 145     {
 146         mw->bond_type = bt;
 147         ExposeWin(x11->disp, mw->wd.self);
 148     }
 149 }
 150
 151 void set_box_type (t_x11 *x11, t_molwin *mw, int bt)
 152 {
 153 #ifdef DEBUG
 154     fprintf(stderr, "mw->boxtype = %d, bt = %d\n", mw->boxtype, bt);
 155 #endif
 156     if (bt != mw->boxtype)
 157     {
 158         if ((bt == esbTrunc && mw->realbox == esbTri) || bt == esbTri || bt == esbNone)
 159         {
 160             mw->boxtype = bt;
 161             ExposeWin(x11->disp, mw->wd.self);
 162         }
 163         else
 164         {
 165             fprintf(stderr, "Can not change rectangular box to truncated octahedron\n");
 166         }
 167     }
 168 }
 169
 170 void done_mw(t_x11 *x11, t_molwin *mw)
 171 {
 172     x11->UnRegisterCallback(x11, mw->wd.self);
 173     sfree(mw);
 174 }
 175
 176 static void draw_atom(Display *disp, Window w, GC gc,
 177                       atom_id ai, iv2 vec2[], unsigned long col[], int size[],
 178                       bool bBall, bool bPlus)
 179 {
 180     int xi, yi;
 181
 182     xi = vec2[ai][XX];
 183     yi = vec2[ai][YY];
 184     XSetForeground(disp, gc, col[ai]);
 185     if (bBall)
 186     {
 187         XFillCircle(disp, w, gc, xi, yi, size[ai]-1);
 188         XSetForeground(disp, gc, BLACK);
 189         XDrawCircle(disp, w, gc, xi, yi, size[ai]);
 190         /*    XSetForeground(disp,gc,WHITE);
 191            XFillCircle(disp,w,gc,xi+4,yi-4,4); */
 192     }
 193     else if (bPlus)
 194     {
 195         XDrawLine(disp, w, gc, xi-MSIZE, yi, xi+MSIZE+1, yi);
 196         XDrawLine(disp, w, gc, xi, yi-MSIZE, xi, yi+MSIZE+1);
 197     }
 198     else
 199     {
 200         XDrawLine(disp, w, gc, xi-1, yi, xi+1, yi);
 201     }
 202
 203 }
 204
 205 /* Global variables */
 206 static rvec gl_fbox, gl_hbox, gl_mhbox;
 207
 208 static void my_init_pbc(matrix box)
 209 {
 210     int i;
 211
 212     for (i = 0; (i < DIM); i++)
 213     {
 214         gl_fbox[i]  =  box[i][i];
 215         gl_hbox[i]  =  gl_fbox[i]*0.5;
 216         gl_mhbox[i] = -gl_hbox[i];
 217     }
 218 }
 219
 220 static bool local_pbc_dx(rvec x1, rvec x2)
 221 {
 222     int  i;
 223     real dx;
 224
 225     for (i = 0; (i < DIM); i++)
 226     {
 227         dx = x1[i]-x2[i];
 228         if (dx > gl_hbox[i])
 229         {
 230             return false;
 231         }
 232         else if (dx <= gl_mhbox[i])
 233         {
 234             return false;
 235         }
 236     }
 237     return true;
 238 }
 239
 240 static void draw_bond(Display *disp, Window w, GC gc,
 241                       atom_id ai, atom_id aj, iv2 vec2[],
 242                       rvec x[], unsigned long col[], int size[], bool bBalls)
 243 {
 244     unsigned long   ic, jc;
 245     int             xi, yi, xj, yj;
 246     int             xm, ym;
 247
 248     if (bBalls)
 249     {
 250         draw_atom(disp, w, gc, ai, vec2, col, size, true, false);
 251         draw_atom(disp, w, gc, aj, vec2, col, size, true, false);
 252     }
 253     else
 254     {
 255         if (local_pbc_dx(x[ai], x[aj]))
 256         {
 257             ic = col[ai];
 258             jc = col[aj];
 259             xi = vec2[ai][XX];
 260             yi = vec2[ai][YY];
 261             xj = vec2[aj][XX];
 262             yj = vec2[aj][YY];
 263
 264             if (ic != jc)
 265             {
 266                 xm = (xi+xj) >> 1;
 267                 ym = (yi+yj) >> 1;
 268
 269                 XSetForeground(disp, gc, ic);
 270                 XDrawLine(disp, w, gc, xi, yi, xm, ym);
 271                 XSetForeground(disp, gc, jc);
 272                 XDrawLine(disp, w, gc, xm, ym, xj, yj);
 273             }
 274             else
 275             {
 276                 XSetForeground(disp, gc, ic);
 277                 XDrawLine(disp, w, gc, xi, yi, xj, yj);
 278             }
 279         }
 280     }
 281 }
 282
 283 int compare_obj(const void *a, const void *b)
 284 {
 285     t_object *oa, *ob;
 286     real      z;
 287
 288     oa = (t_object *)a;
 289     ob = (t_object *)b;
 290
 291     z = oa->z-ob->z;
 292
 293     if (z < 0)
 294     {
 295         return 1;
 296     }
 297     else if (z > 0)
 298     {
 299         return -1;
 300     }
 301     else
 302     {
 303         return 0;
 304     }
 305 }
 306
 307 void create_visibility(t_manager *man)
 308 {
 309     t_object *obj;
 310     int       i;
 311
 312     for (i = 0, obj = man->obj; (i < man->nobj); i++, obj++)
 313     {
 314         if (obj->eV != eVHidden)
 315         {
 316             man->bVis[obj->ai] = true;
 317             switch (obj->eO)
 318             {
 319                 case eOBond:
 320                 case eOHBond:
 321                     man->bVis[obj->aj] = true;
 322                     break;
 323                 default:
 324                     break;
 325             }
 326         }
 327     }
 328 }
 329
 330 void z_fill(t_manager *man, real *zz)
 331 {
 332     t_object *obj;
 333     int       i;
 334
 335     for (i = 0, obj = man->obj; (i < man->nobj); i++, obj++)
 336     {
 337         switch (obj->eO)
 338         {
 339             case eOSingle:
 340                 obj->z = zz[obj->ai];
 341                 break;
 342             case eOBond:
 343             case eOHBond:
 344                 obj->z = (zz[obj->ai] + zz[obj->aj]) * 0.5;
 345                 break;
 346             default:
 347                 break;
 348         }
 349     }
 350 }
 351
 352 int filter_vis(t_manager *man)
 353 {
 354     int          i, nobj, nvis, nhide;
 355     atom_id      ai;
 356     bool         bAdd, *bVis;
 357     t_object    *obj;
 358     t_object    *newobj;
 359
 360     nobj = man->nobj;
 361     snew(newobj, nobj);
 362     obj   = man->obj;
 363     bVis  = man->bVis;
 364     nvis  = 0;
 365     nhide = nobj-1;
 366     for (i = 0; (i < nobj); i++, obj++)
 367     {
 368         ai   = obj->ai;
 369         bAdd = bVis[ai];
 370         if (bAdd)
 371         {
 372             if (obj->eO != eOSingle)
 373             {
 374                 bAdd = bVis[obj->aj];
 375             }
 376         }
 377         if (bAdd)
 378         {
 379             newobj[nvis++] = *obj;
 380         }
 381         else
 382         {
 383             newobj[nhide--] = *obj;
 384         }
 385     }
 386     sfree(man->obj);
 387     man->obj = newobj;
 388
 389     return nvis;
 390 }
 391
 392 void draw_objects(Display *disp, Window w, GC gc, int nobj,
 393                   t_object objs[], iv2 vec2[], rvec x[],
 394                   unsigned long col[], int size[], bool bShowHydro, int bond_type,
 395                   bool bPlus)
 396 {
 397     bool         bBalls;
 398     int          i;
 399     t_object    *obj;
 400
 401     bBalls = false;
 402     switch (bond_type)
 403     {
 404         case eBThin:
 405             XSetLineAttributes(disp, gc, 1, LineSolid, CapNotLast, JoinRound);
 406             break;
 407         case eBFat:
 408             XSetLineAttributes(disp, gc, 3, LineSolid, CapNotLast, JoinRound);
 409             break;
 410         case eBVeryFat:
 411             XSetLineAttributes(disp, gc, 5, LineSolid, CapNotLast, JoinRound);
 412             break;
 413         case eBSpheres:
 414             bBalls = true;
 415             bPlus  = false;
 416             break;
 417         default:
 418             gmx_fatal(FARGS, "Invalid bond_type selected: %d\n", bond_type);
 419             break;
 420     }
 421     for (i = 0; (i < nobj); i++)
 422     {
 423         obj = &(objs[i]);
 424         switch (obj->eO)
 425         {
 426             case eOSingle:
 427                 draw_atom(disp, w, gc, obj->ai, vec2, col, size, bBalls, bPlus);
 428                 break;
 429             case eOBond:
 430                 draw_bond(disp, w, gc, obj->ai, obj->aj, vec2, x, col, size, bBalls);
 431                 break;
 432             case eOHBond:
 433                 if (bShowHydro)
 434                 {
 435                     draw_bond(disp, w, gc, obj->ai, obj->aj, vec2, x, col, size, bBalls);
 436                 }
 437                 break;
 438             default:
 439                 break;
 440         }
 441     }
 442     XSetLineAttributes(disp, gc, 1, LineSolid, CapNotLast, JoinRound);
 443 }
 444
 445 static void v4_to_iv2(vec4 x4, iv2 v2, int x0, int y0, real sx, real sy)
 446 {
 447     real inv_z;
 448
 449     inv_z  = 1.0/x4[ZZ];
 450     v2[XX] = x0+sx*x4[XX]*inv_z;
 451     v2[YY] = y0-sy*x4[YY]*inv_z;
 452 }
 453
 454 static void draw_box(t_x11 *x11, Window w, t_3dview *view, matrix box,
 455                      int x0, int y0, real sx, real sy, int boxtype)
 456 {
 457     rvec        rect_tri[8] =  {
 458         { 0, 0, 0 }, { 1, 0, 0 }, { 1, 1, 0 }, { 0, 1, 0 },
 459         { 0, 0, 1 }, { 1, 0, 1 }, { 1, 1, 1 }, { 0, 1, 1 }
 460     };
 461     int         tr_bonds[12][2] = {
 462         { 0, 1 }, { 1, 2 }, { 2, 3 }, { 3, 0 },
 463         { 4, 5 }, { 5, 6 }, { 6, 7 }, { 7, 4 },
 464         { 0, 4 }, { 1, 5 }, { 2, 6 }, { 3, 7 }
 465     };
 466     static int *edge = NULL;
 467     int         i, j, k, i0, i1;
 468     rvec        corner[NCUCEDGE], box_center;
 469     vec4        x4;
 470     iv2         vec2[NCUCEDGE], tv2;
 471
 472     calc_box_center(view->ecenter, box, box_center);
 473     if (boxtype == esbTrunc)
 474     {
 475         calc_compact_unitcell_vertices(view->ecenter, box, corner);
 476         if (edge == NULL)
 477         {
 478             edge = compact_unitcell_edges();
 479         }
 480
 481         for (i = 0; (i < NCUCEDGE); i++)
 482         {
 483             m4_op(view->proj, corner[i], x4);
 484             v4_to_iv2(x4, vec2[i], x0, y0, sx, sy);
 485         }
 486         XSetForeground(x11->disp, x11->gc, YELLOW);
 487         for (i = 0; i < NCUCEDGE; i++)
 488         {
 489             i0 = edge[2*i];
 490             i1 = edge[2*i+1];
 491             XDrawLine(x11->disp, w, x11->gc,
 492                       vec2[i0][XX], vec2[i0][YY], vec2[i1][XX], vec2[i1][YY]);
 493         }
 494     }
 495     else
 496     {
 497         if (boxtype == esbRect)
 498         {
 499             for (j = 0; (j < DIM); j++)
 500             {
 501                 box_center[j] -= 0.5*box[j][j];
 502             }
 503         }
 504         else
 505         {
 506             for (i = 0; (i < DIM); i++)
 507             {
 508                 for (j = 0; (j < DIM); j++)
 509                 {
 510                     box_center[j] -= 0.5*box[i][j];
 511                 }
 512             }
 513         }
 514         for (i = 0; (i < 8); i++)
 515         {
 516             clear_rvec(corner[i]);
 517             for (j = 0; (j < DIM); j++)
 518             {
 519                 if (boxtype == esbTri)
 520                 {
 521                     for (k = 0; (k < DIM); k++)
 522                     {
 523                         corner[i][k] += rect_tri[i][j]*box[j][k];
 524                     }
 525                 }
 526                 else
 527                 {
 528                     corner[i][j] = rect_tri[i][j]*box[j][j];
 529                 }
 530             }
 531             rvec_inc(corner[i], box_center);
 532             m4_op(view->proj, corner[i], x4);
 533             v4_to_iv2(x4, vec2[i], x0, y0, sx, sy);
 534         }
 535         if (debug)
 536         {
 537             pr_rvecs(debug, 0, "box", box, DIM);
 538             pr_rvecs(debug, 0, "corner", corner, 8);
 539         }
 540         XSetForeground(x11->disp, x11->gc, YELLOW);
 541         for (i = 0; (i < 12); i++)
 542         {
 543             i0 = tr_bonds[i][0];
 544             i1 = tr_bonds[i][1];
 545             XDrawLine(x11->disp, w, x11->gc,
 546                       vec2[i0][XX], vec2[i0][YY], vec2[i1][XX], vec2[i1][YY]);
 547         }
 548     }
 549 }
 550
 551 void set_sizes(t_manager *man, real sx, real sy)
 552 {
 553     int  i;
 554
 555     for (i = 0; (i < man->natom); i++)
 556     {
 557         if (man->bVis[i])
 558         {
 559             man->size[i] = 180*man->vdw[i];
 560         }
 561     }
 562 }
 563
 564 void draw_mol(t_x11 *x11, t_manager *man)
 565 {
 566     static char tstr[2][20];
 567     static int  ntime = 0;
 568     t_windata  *win;
 569     t_3dview   *view;
 570     t_molwin   *mw;
 571     int         i, x0, y0, nvis;
 572     iv2        *vec2;
 573     real        sx, sy;
 574     vec4        x4;
 575
 576     if (!man->status)
 577     {
 578         return;
 579     }
 580
 581     view = man->view;
 582     mw   = man->molw;
 583
 584     win = &(mw->wd);
 585
 586     vec2 = man->ix;
 587     x0   = win->width/2;
 588     y0   = win->height/2;
 589     sx   = win->width/2*view->sc_x;
 590     sy   = win->height/2*view->sc_y;
 591
 592     my_init_pbc(man->box);
 593
 594     /* create_visibility(man); */
 595
 596     for (i = 0; (i < man->natom); i++)
 597     {
 598         if (man->bVis[i])
 599         {
 600             m4_op(view->proj, man->x[i], x4);
 601             man->zz[i] = x4[ZZ];
 602             v4_to_iv2(x4, vec2[i], x0, y0, sx, sy);
 603         }
 604     }
 605     set_sizes(man, sx, sy);
 606
 607     z_fill (man, man->zz);
 608
 609     /* Start drawing */
 610     XClearWindow(x11->disp, win->self);
 611
 612     /* Draw Time */
 613     sprintf(tstr[ntime], "Time: %.3f ps", man->time);
 614     if (strcmp(tstr[ntime], tstr[1-ntime]) != 0)
 615     {
 616         set_vbtime(x11, man->vbox, tstr[ntime]);
 617         ntime = 1-ntime;
 618     }
 619
 620     if (mw->boxtype != esbNone)
 621     {
 622         draw_box(x11, win->self, view, man->box, x0, y0, sx, sy, mw->boxtype);
 623     }
 624
 625     /* Should sort on Z-Coordinates here! */
 626     nvis = filter_vis(man);
 627     if (nvis && man->bSort)
 628     {
 629         qsort(man->obj, nvis, sizeof(man->obj[0]), compare_obj);
 630     }
 631
 632     /* Draw the objects */
 633     draw_objects(x11->disp, win->self, x11->gc,
 634                  nvis, man->obj, man->ix, man->x, man->col, man->size,
 635                  mw->bShowHydrogen, mw->bond_type, man->bPlus);
 636
 637     /* Draw the labels */
 638     XSetForeground(x11->disp, x11->gc, WHITE);
 639     for (i = 0; (i < man->natom); i++)
 640     {
 641         if (man->bLabel[i] && man->bVis[i])
 642         {
 643             XDrawString(x11->disp, win->self, x11->gc, vec2[i][XX]+2, vec2[i][YY]-2,
 644                         man->szLab[i], strlen(man->szLab[i]));
 645         }
 646     }
 647
 648     XSetForeground(x11->disp, x11->gc, x11->fg);
 649 }