大窪 貴洋
131 lines
4.2 KiB
C
131 lines
4.2 KiB
C
#include "dump2analysis.h"
|
|
/* DOCUMENT */
|
|
/* -m angle */
|
|
/* lammpstrjから3原子の角度分布を計算する */
|
|
/* A原子(-a, -x, -s)とB原子(-c, -z, -u)で端の原子を指定する。 */
|
|
/* B原子(-b, -y, -t)で中心原子を指定する。 */
|
|
/* --rcut_abは角度を計算するA-B原子距離のカッとオフ。 default:2.3 */
|
|
/* --rcut_bcは角度を計算するB-C原子距離のカッとオフ。 default:2.3 */
|
|
|
|
void ErrorAngle(){
|
|
printf("Required arguments atoms A, atoms B, and atom C\n");
|
|
printf("(-a, -x, -s) atom A\n");
|
|
printf("(-b, -y, -t) atom B\n");
|
|
printf("(-c, -z, -u) atom B\n");
|
|
printf("------------------------------------------------------------\n");
|
|
printf("Optional arguments for Gr analysis:\n");
|
|
printf("--rcut_ab cut-off for A-B bonds [2.3]\n");
|
|
printf("--rcut_bc cut-off for B-C bonds [2.3]\n");
|
|
printf("------------------------------------------------------------\n");
|
|
printf("Example:\n");
|
|
printf("dump2analysis -m angle -x Si -y O -z Si -i hoge.lammpstrj -o hoge.angle\n");
|
|
exit(0);
|
|
}
|
|
|
|
|
|
void OutputAngle(PARAM *p, int pairs, double *d, double ave, double std){
|
|
FILE *fp;
|
|
int i, j, bunkatu=1800; /* 180degを1800分割, delta_angle = 0.1; */
|
|
double delta_angle, min_angle;
|
|
int *a;
|
|
double ang;
|
|
|
|
a = malloc(sizeof(int) * bunkatu);
|
|
delta_angle = 180.0/(double)bunkatu;
|
|
min_angle = delta_angle * 0.5;
|
|
|
|
for (i=0; i<bunkatu; i++) a[i] = 0;
|
|
for (i=0; i<pairs; i++) {
|
|
j = (int)floor((d[i] - min_angle)/delta_angle);
|
|
a[j]++;
|
|
}
|
|
|
|
fp = fopen(p->outfile, "w");
|
|
fprintf(fp, "# Rcut a-b: %lf\n", RCUT_AB);
|
|
fprintf(fp, "# Rcut b-c: %lf\n", RCUT_BC);
|
|
fprintf(fp, "# Total pairs: %d\n", pairs);
|
|
fprintf(fp, "# angle: %lf +/- %lf\n", ave, std);
|
|
for (i=0; i<bunkatu; i++){
|
|
ang = min_angle + ((double)i+0.5)*delta_angle;
|
|
fprintf(fp, "%lf %d\n", ang, a[i]);
|
|
}
|
|
printf("%s was created.\n", p->outfile);
|
|
|
|
}
|
|
|
|
void EstimateAngle(PARAM *param){
|
|
FILE *f;
|
|
HEAD head;
|
|
ATOMS a, b, c; /* bが中心元素 */
|
|
int i, j, k;
|
|
int s, pairs=0;
|
|
double x, y, z;
|
|
double xa, ya, za;
|
|
double xb, yb, zb;
|
|
double ab, norm_a, norm_b, angle;
|
|
double ra, rb, *d = NULL;
|
|
double ave = 0.0, std = 0.0;
|
|
|
|
if (ArgCheckInputOutput(param) != 0) ErrorAngle();
|
|
SetAtomsSteps(param);
|
|
if (CheckArgAtomSelect(param) != 7) ErrorAngle();
|
|
|
|
Allocate(&a, param->atoms);
|
|
Allocate(&b, param->atoms);
|
|
Allocate(&c, param->atoms);
|
|
f = fopen(param->infile, "r");
|
|
|
|
/* 距離の計算 */
|
|
for (s=0; s<param->steps; s++){
|
|
if (s%10 == 0) {
|
|
printf("Calculating Angle step: %d\r", s);
|
|
fflush(stdout);
|
|
}
|
|
/* 原子の選択 */
|
|
GetData(f, param, &head, &a, &b, &c);
|
|
for(i=0; i<b.atoms; i++){
|
|
for(j=0; j<a.atoms; j++){
|
|
for(k=0; k<c.atoms; k++){
|
|
if ((b.id[i] == a.id[j]) ||
|
|
(b.id[i] == c.id[k]) ||
|
|
(a.id[j] == c.id[k]))
|
|
continue;
|
|
/* a - b bond */
|
|
x = (a.x[j] - b.x[i]) - round(a.x[j] - b.x[i]);
|
|
y = (a.y[j] - b.y[i]) - round(a.y[j] - b.y[i]);
|
|
z = (a.z[j] - b.z[i]) - round(a.z[j] - b.z[i]);
|
|
xa = x*head.A[0][0] + y*head.A[1][0] + z*head.A[2][0];
|
|
ya = x*head.A[0][1] + y*head.A[1][1] + z*head.A[2][1];
|
|
za = x*head.A[0][2] + y*head.A[1][2] + z*head.A[2][2];
|
|
ra = sqrt(pow(xa, 2) + pow(ya, 2) + pow(za, 2));
|
|
|
|
x = (c.x[k] - b.x[i]) - round(c.x[k] - b.x[i]);
|
|
y = (c.y[k] - b.y[i]) - round(c.y[k] - b.y[i]);
|
|
z = (c.z[k] - b.z[i]) - round(c.z[k] - b.z[i]);
|
|
xb = x*head.A[0][0] + y*head.A[1][0] + z*head.A[2][0];
|
|
yb = x*head.A[0][1] + y*head.A[1][1] + z*head.A[2][1];
|
|
zb = x*head.A[0][2] + y*head.A[1][2] + z*head.A[2][2];
|
|
rb = sqrt(pow(xb, 2) + pow(yb, 2) + pow(zb, 2));
|
|
/* 角度の計算 */
|
|
if (ra < RCUT_AB && rb < RCUT_BC){
|
|
norm_a = sqrt(pow(xa, 2) + pow(ya, 2) + pow(za, 2));
|
|
norm_b = sqrt(pow(xb, 2) + pow(yb, 2) + pow(zb, 2));
|
|
ab = xa*xb + ya*yb + za*zb;
|
|
angle = acos(ab/(norm_a*norm_b)) * 180/M_PI;
|
|
d = realloc(d, sizeof(double) * (pairs+1));
|
|
d[pairs] = angle;
|
|
pairs++;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
/* 平均値の計算 */
|
|
for (i=0; i<pairs; i++) ave = ave + d[i];
|
|
ave = ave/pairs;
|
|
/* 標準偏差の計算 */
|
|
for (i=0; i<pairs; i++) std = std + (d[i] - ave)*(d[i] - ave);
|
|
std = sqrt(std/pairs);
|
|
OutputAngle(param, pairs, d, ave, std);
|
|
}
|