visual_atoms.c

   1 /*
   2  * code visualize atoms
   3  *
   4  * author: frank.zirkelbach@physik.uni-augsburg.de
   5  *
   6  */
   7
   8 #define _GNU_SOURCE
   9 #include <stdio.h>
  10 //#include <stdlib.h>
  11 //#include <unistd.h>
  12 //#include <string.h>
  13 //#include <sys/types.h>
  14 //#include <sys/stat.h>
  15 //#include <fcntl.h>
  16
  17 #include "moldyn.h"
  18
  19 /* pse */
  20 #define PSE_NAME
  21 #define PSE_COL
  22 #include "pse.h"
  23 #undef PSE_NAME
  24 #undef PSE_COL
  25
  26 int usage(char *prog) {
  27
  28         printf("\nusage:\n");
  29         printf("  %s <file> <centre atom> <radius> <lc> [marked atom]\n\n",
  30                prog);
  31
  32         return -1;
  33 }
  34
  35 int main(int argc,char **argv) {
  36
  37         t_moldyn moldyn;
  38         t_atom *itom,*jtom;
  39         int j;
  40         int ret;
  41         t_list n[27];
  42         t_list *this;
  43         t_linkcell *lc;
  44         t_3dvec dist;
  45         double d,radius;
  46         double ox,oy,oz;
  47         int ma,ca;
  48         double lac;
  49
  50         if(argc<4) {
  51                 usage(argv[0]);
  52                 return -1;
  53         }
  54
  55         ca=atoi(argv[2]);
  56         radius=atof(argv[3]);
  57         lac=atof(argv[4]);
  58
  59         ma=-1;
  60         if(argc==6)
  61                 ma=atoi(argv[5]);
  62
  63         memset(&moldyn,0,sizeof(t_moldyn));
  64
  65         printf("[visual atoms] reading save file ...\n");
  66         ret=moldyn_read_save_file(&moldyn,argv[1]);
  67         if(ret) {
  68                 printf("[visual atoms] exit!\n");
  69                 return ret;
  70         }
  71
  72         /* link cell init */
  73         moldyn.cutoff=radius;
  74         link_cell_init(&moldyn,VERBOSE);
  75         lc=&(moldyn.lc);
  76
  77         /* serach atoms */
  78         itom=&(moldyn.atom[ca]);
  79         link_cell_neighbour_index(&moldyn,
  80                                   (itom->r.x+moldyn.dim.x/2)/lc->x,
  81                                   (itom->r.y+moldyn.dim.y/2)/lc->y,
  82                                   (itom->r.z+moldyn.dim.z/2)/lc->z,
  83                                   n);
  84
  85
  86         /* prepare offset */
  87         ox=0.0;
  88         if(itom->r.x<0) {
  89                 while((itom->r.x+ox)<(-lac/2.0))
  90                         ox+=lac;
  91         }
  92         else {
  93                 while((itom->r.x+ox)>(lac/2.0))
  94                         ox-=lac;
  95         }
  96
  97         oy=0.0;
  98         if(itom->r.y<0) {
  99                 while((itom->r.y+oy)<(-lac/2.0))
 100                         oy+=lac;
 101         }
 102         else {
 103                 while((itom->r.y+oy)>(lac/2.0))
 104                         oy-=lac;
 105         }
 106
 107         oz=0.0;
 108         if(itom->r.z<0) {
 109                 while((itom->r.z+oz)<(-lac/2.0))
 110                         oz+=lac;
 111         }
 112         else {
 113                 while((itom->r.z+oz)>(lac/2.0))
 114                         oz-=lac;
 115         }
 116
 117
 118         printf("%s %f %f %f %s %f\n",
 119                pse_name[itom->element],itom->r.x+ox,itom->r.y+oy,itom->r.z+oz,
 120                "Green",itom->ekin);
 121
 122         for(j=0;j<27;j++) {
 123                 this=&(n[j]);
 124                 list_reset_f(this);
 125
 126                 if(this->start==NULL)
 127                         continue;
 128
 129                 do {
 130
 131                         jtom=this->current->data;
 132
 133                         if(jtom==itom)
 134                                 continue;
 135
 136                         v3_sub(&dist,&(itom->r),&(jtom->r));
 137                         check_per_bound(&moldyn,&dist);
 138                         d=v3_norm(&dist);
 139
 140                         if(d<=radius) {
 141                                 printf("%s %f %f %f %s %f\n",
 142                                        pse_name[jtom->element],
 143                                        jtom->r.x+ox,jtom->r.y+oy,jtom->r.z+oz,
 144                                        (jtom->tag==ma)?"Red":pse_col[jtom->element],
 145                                        jtom->ekin);
 146                         }
 147
 148                 } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
 149         }
 150
 151         link_cell_shutdown(&moldyn);
 152
 153         moldyn_free_save_file(&moldyn);
 154
 155         return 0;
 156 }