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