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