]> hackdaworld.org Git - physik/posic.git/blob - visual_atoms.c
pse handling again
[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 /* 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 }