1+ #include < iostream>
2+ #include < algorithm>
3+ #include < fstream>
4+ #include < string>
5+ #include < vector>
6+ #include < array>
7+ #include < cmath>
8+ #include < io.h>
9+ #include < direct.h>
10+
11+ int N;
12+ double Lx,Ly,Lz;
13+ std::vector<int > number_type;
14+ const int bins=1000 ;
15+ double rhos[bins],rhof[bins],rhosav[bins]{0 },rhofav[bins]{0 }; // short and long
16+ int stats=0 ;
17+ int save_step=100 ;
18+
19+ // ==============================================================
20+ // Checks if directory exists
21+ // --------------------------------------------------------------
22+ void dexist (std::string dirname){
23+ if (_access (dirname.c_str (),0 )==-1 )
24+ _mkdir (dirname.c_str ());
25+ return ;
26+ }
27+
28+ void conc (std::vector<std::array<double ,3 >> &x)
29+ {
30+ int zi,n_type=0 ;
31+ double V=Lx*Ly*Lz;
32+ double z, ADD=bins/V;
33+ for (int i=0 ;i<bins;i++)
34+ {
35+ rhos[i]=0 ;
36+ rhof[i]=0 ;
37+ }
38+
39+ stats ++;
40+
41+ if (stats%save_step==0 )
42+ {
43+ dexist (" dumps"
44+ std::string fname;
45+ fname=" rho_ave" std::to_string (stats);
46+ std::fstream write;
47+ write.open (fname,std::ios::out);
48+ for (int i=0 ;i<bins;i++){
49+ double z = (i+0.5 )*Lx/bins;
50+ write<<z<<' \t ' ' \t ' ' \t '
51+ }
52+ write.close ();
53+
54+ for (int i=0 ;i<bins;i++)
55+ {
56+ rhosav[i]=0 ;
57+ rhofav[i]=0 ;
58+ }
59+ }
60+
61+ for (int i=0 ;i<x.size ();i++)
62+ {
63+ z=x[i][2 ]-std::floor (x[i][2 ]/Lz)*Lz;
64+ zi=bins*z/Lz;
65+ if (i>number_type[n_type])
66+ rhof[zi]+=ADD;
67+ else
68+ rhos[zi]+=ADD;
69+ }
70+ for (int i=0 ;i<bins;i++){
71+ rhofav[i] += (rhof[i]-rhofav[i])/(1.0 *stats);
72+ rhosav[i] += (rhos[i]-rhosav[i])/(1.0 *stats);
73+ }
74+ }
75+
76+ int get_para (std::fstream &read)
77+ {
78+ double low_x=0 ,low_y=0 ,low_z=0 ,high_x=0 ,high_y=0 ,high_z=0 ;
79+ std::string aline;
80+ if (!std::getline (read,aline))
81+ return 0 ;
82+
83+ read>>N;
84+ for (int i=0 ;i<8 ;i++)
85+ if (!std::getline (read,aline))
86+ return 0 ;
87+
88+ read>>low_x>>high_x;
89+ if (!std::getline (read,aline))
90+ return 0 ;
91+ read>>low_y>>high_y;
92+ if (!std::getline (read,aline))
93+ return 0 ;
94+ read>>low_z>>high_z;
95+ Lx=high_x-low_x;
96+ Ly=high_y-low_y;
97+ Lz=high_z-low_z;
98+
99+ return 1 ;
100+ }
101+ // 设计成这样不需要new和delete
102+ int read_data (std::vector<std::array<double ,3 >> &x,std::fstream &read)
103+ {
104+ int NB,atom_type,old_atom_type=0 ;
105+
106+ std::string aline;
107+
108+ read>>NB;
109+ if (NB!=N)
110+ {
111+ std::cout<<" atom number isn't same!"
112+ return 0 ;
113+ }
114+
115+ if (!std::getline (read,aline)) return 0 ;
116+ if (!std::getline (read,aline)) return 0 ;
117+
118+ for (int i=0 ;i<NB;i++)
119+ {
120+ x.push_back (std::array<double ,3 >{0 });
121+ read>>atom_type>>x[i][0 ]>>x[i][1 ]>>x[i][2 ];
122+ if (atom_type!=old_atom_type)
123+ {
124+ number_type.push_back (i);
125+ old_atom_type=atom_type;
126+ }
127+ }
128+
129+ return 1 ;
130+ }
131+
132+ void cal_concertration (std::string readFileName,std::string dataName)
133+ {
134+ std::fstream read;
135+ read.open (readFileName,std::ios::in);
136+ if (!get_para (read))
137+ {
138+ std::cout<<" error of reading"
139+ exit (-1 );
140+ }
141+
142+ std::vector<std::array<double ,3 >> x;
143+ read.close ();
144+
145+ read.open (dataName,std::ios::in);
146+ while (read_data (x,read))
147+ {
148+
149+ }
150+
151+ read.close ();
152+
153+ }
0 commit comments