first commit, Readme.html更新して、git管理することにした。

2023-05-26 20:13:05 +09:00
commit ce85e3907f
26 changed files with 158505 additions and 0 deletions
@@ -0,0 +1,130 @@
+#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 Gr 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);
+}