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