-
Notifications
You must be signed in to change notification settings - Fork 7
/
msd.cpp
188 lines (161 loc) · 6.21 KB
/
msd.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
#include "driver.h"
#include "math.h"
#include <vector>
/*------------------------------------------------------------------------------
* Method to compute the Mean square displacement of selected atoms.
*
* Note that in this case the dump file read must contain the image info.
*----------------------------------------------------------------------------*/
void Driver::compute_msd()
{
char str[MAXLINE], header[MAXLINE];
printf("\n"); for (int i = 0; i < 5; ++i) printf("====");
printf(" Mean-squared displacement calculation ");
for (int i = 0; i < 5; ++i) printf("===="); printf("\n");
// selection commands; the selections will be intersected among frames.
char selcmd[MAXLINE];
one = all[istr];
int natom = one->natom;
int *insel = memory->create(insel, natom+1, "insel");
for (int id = 1; id <= natom; ++id) insel[id] = 1;
int nused = 0;
// selection commands for atoms
while (1){
printf("\nPlease input the selection command for atoms, `h` for help [all]: ");
input->read_stdin(str);
if (count_words(str) > 0){
strcpy(selcmd, str);
char *ptr = strtok(str," \n\t\r\f");
if (strcmp(ptr,"h") == 0){ one->SelHelp(); continue; }
} else strcpy(selcmd,"all\n");
if (strcmp(selcmd, "all") == 0) break;
// determine if selection is based on first frame or all frames
int flag_sel = 2;
printf("Should the selection be applied to:\n");
printf(" 1. the first frame only;\n");
printf(" 2. all frames and get the common selection.\n");
printf("Your choice [2]: ");
input->read_stdin(str);
if (count_words(str) > 0){
char *ptr = strtok(str," \n\t\r\f");
flag_sel = MIN(2, MAX(atoi(ptr),1));
}
// check the selection command
for (int img = istr; img <= iend; img += inc){
one = all[img];
if (one->natom != natom){
printf("\nError: number of atoms from frame %d differs from previous ones!\n", img+1);
for (int i = 0; i < 20; ++i) printf("===="); printf("\n");
return;
}
if (one->image == NULL){
printf("\nError: image info not available for frame %d! MSD calculation is not possible.\n", img+1);
for (int i = 0; i < 20; ++i) printf("===="); printf("\n");
return;
}
if (flag_sel == 2 || img == istr){
one->selection(selcmd);
for (int id = 1; id <= natom; ++id) insel[id] &= one->atsel[id];
}
// cartesian coordinate needed for MSD calculations
one->dir2car();
++nused;
}
int nsel = 0;
for (int id = 1; id <= natom; ++id) nsel += insel[id];
if (nsel < 1){
printf("It seems that no atom is selected, do you wish to make another selection? (y/n)[n]: ");
input->read_stdin(str);
if (count_words(str) > 0){
char *ptr = strtok(str," \n\t\r\f");
if (strcmp(ptr,"y")== 0 || strcmp(ptr,"Y")==0) continue;
}
for (int i = 0; i < 20; ++i) printf("===="); printf("\n");
return;
}
break;
}
double **msd = memory->create(msd, nused, 6, "msd");
for (int img = 0; img < nused; ++img)
for (int idim = 0; idim < 6; ++idim) msd[img][idim] = 0.;
for (int ir = istr; ir <= iend; ir += inc){
one = all[ir];
int istep = 0;
for (int in = ir+inc; in <= iend; in += inc){
DumpAtom *now = all[in];
msd[++istep][0] = double(now->tstep - one->tstep);
if (one->triclinic == 0) {
for (int id = 1; id <= natom; ++id){
if (insel[id] == 0) continue;
int xbox = one->image[id][0];
int ybox = one->image[id][1];
int zbox = one->image[id][2];
double x0 = one->atpos[id][0] + xbox*one->lx;
double y0 = one->atpos[id][1] + ybox*one->ly;
double z0 = one->atpos[id][2] + zbox*one->lz;
xbox = now->image[id][0];
ybox = now->image[id][1];
zbox = now->image[id][2];
double dx = now->atpos[id][0] - x0 + xbox*now->lx;
double dy = now->atpos[id][1] - y0 + ybox*now->ly;
double dz = now->atpos[id][2] - z0 + zbox*now->lz;
msd[istep][1] += dx*dx;
msd[istep][2] += dy*dy;
msd[istep][3] += dz*dz;
msd[istep][4] += dx*dx + dy*dy + dz*dz;
msd[istep][5] += 1.;
}
} else {
for (int id = 1; id <= natom; ++id){
if (insel[id] == 0) continue;
int xbox = one->image[id][0];
int ybox = one->image[id][1];
int zbox = one->image[id][2];
double x0 = one->atpos[id][0] + xbox * one->lx + ybox * one->xy + zbox * one->xz;
double y0 = one->atpos[id][1] + ybox * one->ly + zbox * one->yz;
double z0 = one->atpos[id][2] + zbox * one->lz;
xbox = now->image[id][0];
ybox = now->image[id][1];
zbox = now->image[id][2];
double dx = now->atpos[id][0] - x0 + xbox * now->lx + ybox * now->xy + zbox * now->xz;
double dy = now->atpos[id][1] - y0 + ybox * now->ly + zbox * now->yz;
double dz = now->atpos[id][2] - z0 + zbox * now->lz;
msd[istep][1] += dx*dx;
msd[istep][2] += dy*dy;
msd[istep][3] += dz*dz;
msd[istep][4] += dx*dx + dy*dy + dz*dz;
msd[istep][5] += 1.;
}
}
}
}
// average
for (int it = 1; it < nused; ++it){
if (msd[it][5] > 0.)
for (int idim = 1; idim <= 4; ++idim) msd[it][idim] /= msd[it][5];
}
// output the result
printf("Please input the file to output the MSD info [msd.dat]: ");
input->read_stdin(str);
char *ptr = strtok(str, " \n\t\r\f");
if (ptr == NULL){
strcpy(str, "msd.dat");
ptr = strtok(str, " \n\t\r\f");
}
char *fname = new char [strlen(ptr)+1];
strcpy(fname, ptr);
ConfirmOverwrite(fname);
FILE *fp = fopen(fname,"w");
fprintf(fp,"# MSD for: %s", selcmd);
fprintf(fp,"# d-step x y z r\n");
for (int it = 0; it < nused; ++it){
fprintf(fp,"%d %lg %lg %lg %lg\n", int(msd[it][0]), msd[it][1], msd[it][2], msd[it][3], msd[it][4]);
}
fclose(fp);
memory->destroy(msd);
printf("\n%d images were used in the evaluation of MSD, which is written to %s\n", nused, fname);
delete []fname;
for (int i = 0; i < 20; ++i) printf("===="); printf("\n");
return;
}
/*------------------------------------------------------------------------------*/