]> hackdaworld.org Git - physik/posic.git/blob - visual_atoms.c
timestep -> tau
[physik/posic.git] / 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 }