src/ngmx/nmol.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 3.0
  11  *
  12  * Copyright (c) 1991-2001
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * This program is free software; you can redistribute it and/or
  17  * modify it under the terms of the GNU General Public License
  18  * as published by the Free Software Foundation; either version 2
  19  * of the License, or (at your option) any later version.
  20  *
  21  * If you want to redistribute modifications, please consider that
  22  * scientific software is very special. Version control is crucial -
  23  * bugs must be traceable. We will be happy to consider code for
  24  * inclusion in the official distribution, but derived work must not
  25  * be called official GROMACS. Details are found in the README & COPYING
  26  * files - if they are missing, get the official version at www.gromacs.org.
  27  *
  28  * To help us fund GROMACS development, we humbly ask that you cite
  29  * the papers on the package - you can find them in the top README file.
  30  *
  31  * Do check out http://www.gromacs.org , or mail us at gromacs@gromacs.org .
  32  *
  33  * And Hey:
  34  * Giving Russians Opium May Alter Current Situation
  35  */
  36 static char *SRCID_nmol_c = "$Id$";
  37 #include <math.h>
  38 #include "sysstuff.h"
  39 #include "string.h"
  40 #include "smalloc.h"
  41 #include "macros.h"
  42 #include "xutil.h"
  43 #include "3dview.h"
  44 #include "fatal.h"
  45 #include "buttons.h"
  46 #include "manager.h"
  47 #include "nmol.h"
  48 #include "vec.h"
  49 #include "assert.h"
  50 #include "txtdump.h"
  51 #include "pbc.h"
  52
  53 #define MSIZE 4
  54
  55 static bool MWCallBack(t_x11 *x11,XEvent *event, Window w, void *data)
  56 {
  57   t_molwin *mw;
  58   Window   To;
  59   XEvent   letter;
  60
  61   mw=(t_molwin *)data;
  62   To=mw->wd.Parent;
  63   letter.type=ClientMessage;
  64   letter.xclient.display=x11->disp;
  65   letter.xclient.window=To;
  66   letter.xclient.message_type=0;
  67   letter.xclient.format=32;
  68   switch(event->type) {
  69   case Expose:
  70     /* Do Not draw anything, but signal parent instead, he will
  71      * coordinate drawing.
  72      */
  73     letter.xclient.data.l[0]=IDDRAWMOL;
  74     letter.xclient.data.l[1]=Button1;
  75     XSendEvent(x11->disp,To,True,0,&letter);
  76     break;
  77   case ButtonPress:
  78 #ifdef DEBUG
  79     printf("Molwindow: Buttonpress\n");
  80 #endif
  81     letter.xclient.data.l[0]=IDLABEL;
  82     letter.xclient.data.l[1]=(long)event->xbutton.button;
  83     letter.xclient.data.l[2]=event->xbutton.x;
  84     letter.xclient.data.l[3]=event->xbutton.y;
  85     XSendEvent(x11->disp,To,True,0,&letter);
  86     break;
  87   case ConfigureNotify:
  88     mw->wd.width=event->xconfigure.width;
  89     mw->wd.height=event->xconfigure.height;
  90     break;
  91   default:
  92     break;
  93   }
  94   return FALSE;
  95 }
  96
  97 void set_def (t_molwin *mw,matrix box)
  98 {
  99   mw->bShowHydrogen=TRUE;
 100   mw->bond_type=eBFat;
 101   mw->boxtype=esbRect;
 102   mw->realbox=TRICLINIC(box) ? esbTri : esbRect;
 103 }
 104
 105 t_molwin *init_mw(t_x11 *x11,Window Parent,
 106                   int x,int y,int width,int height,
 107                   unsigned long fg,unsigned long bg,
 108                   matrix box)
 109 {
 110   t_molwin *mw;
 111
 112   snew(mw,1);
 113   set_def(mw,box);
 114
 115   InitWin(&mw->wd,x,y,width,height,1,"Mol Window");
 116
 117   mw->wd.Parent=Parent;
 118   mw->wd.self=XCreateSimpleWindow(x11->disp,Parent,x,y,width,height,1,fg,bg);
 119   x11->RegisterCallback(x11,mw->wd.self,Parent,MWCallBack,mw);
 120   x11->SetInputMask(x11,mw->wd.self,
 121                     ExposureMask | StructureNotifyMask |
 122                     ButtonPressMask);
 123   return mw;
 124 }
 125
 126 void map_mw(t_x11 *x11,t_molwin *mw)
 127 {
 128   XMapWindow(x11->disp,mw->wd.self);
 129 }
 130
 131 bool toggle_hydrogen(t_x11 *x11,t_molwin *mw)
 132 {
 133   mw->bShowHydrogen=!mw->bShowHydrogen;
 134   ExposeWin(x11->disp,mw->wd.self);
 135
 136   return mw->bShowHydrogen;
 137 }
 138
 139 void set_bond_type(t_x11 *x11,t_molwin *mw,int bt)
 140 {
 141   if (bt != mw->bond_type) {
 142     mw->bond_type=bt;
 143     ExposeWin(x11->disp,mw->wd.self);
 144   }
 145 }
 146
 147 void set_box_type (t_x11 *x11,t_molwin *mw,int bt)
 148 {
 149 #ifdef DEBUG
 150   fprintf(stderr,"mw->boxtype = %d, bt = %d\n",mw->boxtype,bt);
 151 #endif
 152   if (bt != mw->boxtype) {
 153     if ((((bt == esbTrunc) || (bt == esbTri)) &&
 154          (mw->realbox == esbTri)) || (bt == esbNone)) {
 155       mw->boxtype = bt;
 156       ExposeWin(x11->disp,mw->wd.self);
 157     }
 158     else
 159       fprintf(stderr,"Can not change rectangular box to triclinic or 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                       bool bBall,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 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[],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   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[],bool bShowHydro,int bond_type,
 351                   bool bPlus)
 352 {
 353   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     fatal_error(0,"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(box,box_center);
 424   if (boxtype == esbTrunc) {
 425     calc_compact_unitcell_vertices(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 == -1)
 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