src/ngmx/nmol.c

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