search_bonds.c

   1 /*
   2  * code searching for special types of bonds
   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 #include "potentials/albe.h"
  19
  20 int usage(char *prog) {
  21
  22         printf("\nusage:\n");
  23         printf("  %s <file> <type> <bondlen> <+->\n",prog);
  24         printf("    type: a - 0-0 bond\n");
  25         printf("    type: b - 1-1 bond\n");
  26         printf("    type: c - 0-1 / 1-0  bond\n\n");
  27
  28         return -1;
  29 }
  30
  31 int main(int argc,char **argv) {
  32
  33         t_moldyn moldyn;
  34         t_atom *itom,*jtom;
  35         int i,j;
  36         int ret;
  37         t_list n[27];
  38         t_list *this;
  39         t_linkcell *lc;
  40         t_3dvec dist;
  41         double d,bondlen,bondpm;
  42
  43         if(argc!=5) {
  44                 usage(argv[0]);
  45                 return -1;
  46         }
  47
  48         bondlen=atof(argv[3]);
  49         bondpm=atof(argv[4]);
  50
  51         memset(&moldyn,0,sizeof(t_moldyn));
  52
  53         printf("[search bonds] reading save file ...\n");
  54         ret=moldyn_read_save_file(&moldyn,argv[1]);
  55         if(ret) {
  56                 printf("[search bonds] exit!\n");
  57                 return ret;
  58         }
  59
  60         /* link cell init */
  61         moldyn.cutoff=bondlen+bondpm;
  62         link_cell_init(&moldyn,VERBOSE);
  63         lc=&(moldyn.lc);
  64
  65         /* analyzing ... */
  66         for(i=0;i<moldyn.count;i++) {
  67                 itom=&(moldyn.atom[i]);
  68                 link_cell_neighbour_index(&moldyn,
  69                                           (itom->r.x+moldyn.dim.x/2)/lc->x,
  70                                           (itom->r.y+moldyn.dim.y/2)/lc->y,
  71                                           (itom->r.z+moldyn.dim.z/2)/lc->z,
  72                                           n);
  73                 for(j=0;j<27;j++) {
  74                         this=&(n[j]);
  75                         list_reset_f(this);
  76
  77                         if(this->start==NULL)
  78                                 continue;
  79
  80                         do {
  81
  82                                 jtom=this->current->data;
  83
  84                                 if(jtom<itom)
  85                                         continue;
  86
  87                                 switch(argv[2][0]) {
  88                                         case 'a':
  89                                                 if(itom->brand!=0)
  90                                                         continue;
  91                                                 if(jtom->brand!=0)
  92                                                         continue;
  93                                                 break;
  94                                         case 'b':
  95                                                 if(itom->brand!=1)
  96                                                         continue;
  97                                                 if(jtom->brand!=1)
  98                                                         continue;
  99                                                 break;
 100                                         default:
 101                                                 if(itom->brand==jtom->brand)
 102                                                         continue;
 103                                                 break;
 104                                 }
 105
 106                                 v3_sub(&dist,&(itom->r),&(jtom->r));
 107                                 check_per_bound(&moldyn,&dist);
 108                                 d=v3_norm(&dist);
 109
 110                                 if((d<=bondlen+bondpm)&(d>=bondlen-bondpm)) {
 111
 112                                         printf(" # atoms %d/%d %d/%d - %f\n",
 113                                                itom->tag,itom->brand,
 114                                                jtom->tag,jtom->brand,d);
 115
 116                                 }
 117
 118                         } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
 119                 }
 120         }
 121
 122         link_cell_shutdown(&moldyn);
 123
 124         moldyn_free_save_file(&moldyn);
 125
 126         return 0;
 127 }