Commit 8e1b3116 authored by Aidan Thompson's avatar Aidan Thompson
Browse files

Added Clebsch-Gordan coefficients

parent aa9b805c

src/compute_orientorder_atom.cpp

+330 −31
Original line number Diff line number Diff line
@@ -56,6 +56,7 @@ ComputeOrientOrderAtom::ComputeOrientOrderAtom(LAMMPS *lmp, int narg, char **arg 


  nnn = 12;

  cutsq = 0.0;

  wlflag = 0;

  qlcompflag = 0;



  // specify which orders to request
@@ -103,6 +104,13 @@ ComputeOrientOrderAtom::ComputeOrientOrderAtom(LAMMPS *lmp, int narg, char **arg 
        if (qlist[il] > qmax) qmax = qlist[il];

      }

      iarg += nqlist;

    } else if (strcmp(arg[iarg],"wlparams") == 0) {

      if (iarg+2 > narg)

        error->all(FLERR,"Illegal compute orientorder/atom command");

      if (strcmp(arg[iarg+1],"yes") == 0) wlflag = 1;

      else if (strcmp(arg[iarg+1],"no") == 0) wlflag = 0;

      else error->all(FLERR,"Illegal compute orientorder/atom command");

      iarg += 2;

    } else if (strcmp(arg[iarg],"components") == 0) {

      qlcompflag = 1;

      if (iarg+2 > narg)
@@ -130,8 +138,9 @@ ComputeOrientOrderAtom::ComputeOrientOrderAtom(LAMMPS *lmp, int narg, char **arg 
    } else error->all(FLERR,"Illegal compute orientorder/atom command");

  }



  if (qlcompflag) ncol = nqlist + 2*(2*qlcomp+1);

  else ncol = nqlist;

  ncol = nqlist;

  if (wlflag) ncol += nqlist;

  if (qlcompflag) ncol += nqlist + 2*(2*qlcomp+1);



  peratom_flag = 1;

  size_peratom_cols = ncol;
@@ -151,7 +160,7 @@ ComputeOrientOrderAtom::~ComputeOrientOrderAtom() 
  memory->destroy(qlist);

  memory->destroy(qnm_r);

  memory->destroy(qnm_i);



  memory->destroy(cglist);

}



/* ---------------------------------------------------------------------- */
@@ -183,6 +192,8 @@ void ComputeOrientOrderAtom::init() 
    if (strcmp(modify->compute[i]->style,"orientorder/atom") == 0) count++;

  if (count > 1 && comm->me == 0)

    error->warning(FLERR,"More than one compute orientorder/atom");



  if (wlflag) init_clebsch_gordan();

}



/* ---------------------------------------------------------------------- */
@@ -274,8 +285,8 @@ void ComputeOrientOrderAtom::compute_peratom() 
      // if not nnn neighbors, order parameter = 0;



      if ((ncount == 0) || (ncount < nnn)) {

        for (int il = 0; il < nqlist; il++)

          qn[il] = 0.0;

        for (int jj = 0; jj < ncol; jj++)

          qn[jj] = 0.0;

        continue;

      }
@@ -461,52 +472,71 @@ void ComputeOrientOrderAtom::calc_boop(double **rlist, 
    }

  }



  // convert sums to averages



  double facn = 1.0 / ncount;

  for (int il = 0; il < nqlist; il++) {

    int l = qlist[il];

    for(int m = 0; m < 2*l+1; m++) {

      qnm_r[il][m] *= facn;

      qnm_i[il][m] *= facn;

    }

  }



  // calculate Q_l



  double fac = sqrt(MY_4PI) / ncount;

  double facpi = sqrt(MY_4PI);

  double normfac = 0.0;

  int jcount = 0;

  int jj = 0;

  for (int il = 0; il < nqlist; il++) {

    int l = qlist[il];

    double qm_sum = 0.0;

    for(int m = 0; m < 2*l+1; m++)

      qm_sum += qnm_r[il][m]*qnm_r[il][m] + qnm_i[il][m]*qnm_i[il][m];



    qn[jcount++] = fac * sqrt(qm_sum / (2*l+1));

    for(int m = 0; m < 2*l+1; m++)

      printf("Q_l = %d %g %g\n",l, qnm_r[il][m], qnm_i[il][m]);



    qn[jj++] = facpi * sqrt(qm_sum / (2*l+1) );

    if (qlcompflag && iqlcomp == il) normfac = 1.0/sqrt(qm_sum);



  }



  // TODO: 

  // 1. Need to allocate extra memory in qn[] for this option

  // 2. Need to add keyword option

  // 1. [done]Need to allocate extra memory in qnarray[] for this option

  // 2. [done]Need to add keyword option

  // 3. Need to caclulate Clebsch-Gordan/Wigner 3j coefficients

  // 4. Compate to bcc values in /Users/athomps/netapp/codes/MatMiner/matminer/matminer/featurizers/boop



//   // calculate W_l



//   if (wlflag) {

//     for (int il = 0; il < nqlist; il++) {

//       int l = qlist[il];

//       double wlsum = 0.0;

//       for(int m1 = 0; m1 < 2*l+1; m1++) {

//         for(int m2 = MAX(0,l-m1); m2 < MIN(2*l+1,3*l-m1; m2++)) {

//           int m = m1 + m2 - l;

//           qm1qm2_r = qnm_r[il][m1]*qnm_r[il][m2] - qnm_i[il][m1]*qnm_i[il][m2];

//           qm1qm2_i = qnm_r[il][m1]*qnm_i[il][m2] + qnm_i[il][m1]*qnm_r[il][m2];

//           wlsum += (qm1qm2_r*qnm_r[il][m] + qm1qm2_i*qnm_i[il][m])*cg; 

//         }

//       }

//       qn[jcount++] = wlsum;

//     }

//   }

  //     (Can try getting them from boop.py first)

  // 5. Compate to bcc values in /Users/athomps/netapp/codes/MatMiner/matminer/matminer/featurizers/boop.py



  // calculate W_l



  if (wlflag) {

    int idxcg_count = 0;

    for (int il = 0; il < nqlist; il++) {

      int l = qlist[il];

      double wlsum = 0.0;

      for(int m1 = 0; m1 < 2*l+1; m1++) {

        for(int m2 = MAX(0,l-m1); m2 < MIN(2*l+1,3*l-m1+1); m2++) {

          int m = m1 + m2 - l;

          double qm1qm2_r = qnm_r[il][m1]*qnm_r[il][m2] - qnm_i[il][m1]*qnm_i[il][m2];

          double qm1qm2_i = qnm_r[il][m1]*qnm_i[il][m2] + qnm_i[il][m1]*qnm_r[il][m2];

          wlsum += (qm1qm2_r*qnm_r[il][m] + qm1qm2_i*qnm_i[il][m])*cglist[idxcg_count];

          idxcg_count++;

        }

      }

      // This matches boop.py output, with cg==1

      printf("W_l = %d %g\n",l,wlsum);

      qn[jj++] = wlsum/sqrt(2*l+1);

    }

  }



  // output of the complex vector



  if (qlcompflag) {

    for(int m = 0; m < 2*qlcomp+1; m++) {

      qn[jcount++] = qnm_r[iqlcomp][m] * normfac;

      qn[jcount++] = qnm_i[iqlcomp][m] * normfac;

      qn[jj++] = qnm_r[iqlcomp][m] * normfac;

      qn[jj++] = qnm_i[iqlcomp][m] * normfac;

    }

  }
@@ -567,3 +597,272 @@ double ComputeOrientOrderAtom::associated_legendre(int l, int m, double x) 


  return p;

}



/* ----------------------------------------------------------------------

   assign Clebsch-Gordan coefficients with l1=l2=l

   using the quasi-binomial formula VMK 8.2.1(3)

------------------------------------------------------------------------- */



void ComputeOrientOrderAtom::init_clebsch_gordan()

{

  double sum,dcg,sfaccg;

  int m, aa2, bb2, cc2, j1, j2, j;

  int ifac, idxcg_count;



  idxcg_count = 0;

  for (int il = 0; il < nqlist; il++) {

    int l = qlist[il];

    for(int m1 = 0; m1 < 2*l+1; m1++) {

      for(int m2 = MAX(0,l-m1); m2 < MIN(2*l+1,3*l-m1+1); m2++) {

        idxcg_count++;

      }

    }

  }

  idxcg_max = idxcg_count;

  memory->create(cglist, idxcg_max, "computeorientorderatom:cglist");



  idxcg_count = 0;

  for (int il = 0; il < nqlist; il++) {

    int l = qlist[il];

    j1 = j2 = j = l;

    for(int m1 = 0; m1 < 2*l+1; m1++) {

        aa2 = m1 - j1;

        for(int m2 = MAX(0,l-m1); m2 < MIN(2*l+1,3*l-m1+1); m2++) {

          bb2 = m2 - j2;

          m = (aa2 + bb2 + j) / 2;

          

          sum = 0.0;

          for (int z = MAX(0, MAX(-(j - j2 + aa2)

                                  , -(j - j1 - bb2) ));

               z <= MIN((j1 + j2 - j),

                        MIN((j1 - aa2), (j2 + bb2)));

               z++) {

            ifac = z % 2 ? -1 : 1;

            sum += ifac /

              (factorial(z) *

               factorial((j1 + j2 - j) - z) *

               factorial((j1 - aa2) - z) *

               factorial((j2 + bb2) - z) *

               factorial((j - j2 + aa2) + z) *

               factorial((j - j1 - bb2) + z));

          }

          

          cc2 = m - j;

          dcg = deltacg(j1, j2, j);

          sfaccg = sqrt(factorial((j1 + aa2)) *

                        factorial((j1 - aa2)) *

                        factorial((j2 + bb2)) *

                        factorial((j2 - bb2)) *

                        factorial((j  + cc2)) *

                        factorial((j  - cc2)) *

                        (j + 1));

          

          cglist[idxcg_count] = sum * dcg * sfaccg;

          idxcg_count++;

        }

      }

  }

}



/* ----------------------------------------------------------------------

   factorial n, wrapper for precomputed table

------------------------------------------------------------------------- */



double ComputeOrientOrderAtom::factorial(int n)

{

  if (n < 0 || n > nmaxfactorial) {

    char str[128];

    sprintf(str, "Invalid argument to factorial %d", n);

    error->all(FLERR, str);

  }



  return nfac_table[n];

}



/* ----------------------------------------------------------------------

   factorial n table, size SNA::nmaxfactorial+1

------------------------------------------------------------------------- */



const double ComputeOrientOrderAtom::nfac_table[] = {

  1,

  1,

  2,

  6,

  24,

  120,

  720,

  5040,

  40320,

  362880,

  3628800,

  39916800,

  479001600,

  6227020800,

  87178291200,

  1307674368000,

  20922789888000,

  355687428096000,

  6.402373705728e+15,

  1.21645100408832e+17,

  2.43290200817664e+18,

  5.10909421717094e+19,

  1.12400072777761e+21,

  2.5852016738885e+22,

  6.20448401733239e+23,

  1.5511210043331e+25,

  4.03291461126606e+26,

  1.08888694504184e+28,

  3.04888344611714e+29,

  8.8417619937397e+30,

  2.65252859812191e+32,

  8.22283865417792e+33,

  2.63130836933694e+35,

  8.68331761881189e+36,

  2.95232799039604e+38,

  1.03331479663861e+40,

  3.71993326789901e+41,

  1.37637530912263e+43,

  5.23022617466601e+44,

  2.03978820811974e+46,

  8.15915283247898e+47,

  3.34525266131638e+49,

  1.40500611775288e+51,

  6.04152630633738e+52,

  2.65827157478845e+54,

  1.1962222086548e+56,

  5.50262215981209e+57,

  2.58623241511168e+59,

  1.24139155925361e+61,

  6.08281864034268e+62,

  3.04140932017134e+64,

  1.55111875328738e+66,

  8.06581751709439e+67,

  4.27488328406003e+69,

  2.30843697339241e+71,

  1.26964033536583e+73,

  7.10998587804863e+74,

  4.05269195048772e+76,

  2.35056133128288e+78,

  1.3868311854569e+80,

  8.32098711274139e+81,

  5.07580213877225e+83,

  3.14699732603879e+85,

  1.98260831540444e+87,

  1.26886932185884e+89,

  8.24765059208247e+90,

  5.44344939077443e+92,

  3.64711109181887e+94,

  2.48003554243683e+96,

  1.71122452428141e+98,

  1.19785716699699e+100,

  8.50478588567862e+101,

  6.12344583768861e+103,

  4.47011546151268e+105,

  3.30788544151939e+107,

  2.48091408113954e+109,

  1.88549470166605e+111,

  1.45183092028286e+113,

  1.13242811782063e+115,

  8.94618213078297e+116,

  7.15694570462638e+118,

  5.79712602074737e+120,

  4.75364333701284e+122,

  3.94552396972066e+124,

  3.31424013456535e+126,

  2.81710411438055e+128,

  2.42270953836727e+130,

  2.10775729837953e+132,

  1.85482642257398e+134,

  1.65079551609085e+136,

  1.48571596448176e+138,

  1.3520015276784e+140,

  1.24384140546413e+142,

  1.15677250708164e+144,

  1.08736615665674e+146,

  1.03299784882391e+148,

  9.91677934870949e+149,

  9.61927596824821e+151,

  9.42689044888324e+153,

  9.33262154439441e+155,

  9.33262154439441e+157,

  9.42594775983835e+159,

  9.61446671503512e+161,

  9.90290071648618e+163,

  1.02990167451456e+166,

  1.08139675824029e+168,

  1.14628056373471e+170,

  1.22652020319614e+172,

  1.32464181945183e+174,

  1.44385958320249e+176,

  1.58824554152274e+178,

  1.76295255109024e+180,

  1.97450685722107e+182,

  2.23119274865981e+184,

  2.54355973347219e+186,

  2.92509369349301e+188,

  3.3931086844519e+190,

  3.96993716080872e+192,

  4.68452584975429e+194,

  5.5745857612076e+196,

  6.68950291344912e+198,

  8.09429852527344e+200,

  9.8750442008336e+202,

  1.21463043670253e+205,

  1.50614174151114e+207,

  1.88267717688893e+209,

  2.37217324288005e+211,

  3.01266001845766e+213,

  3.8562048236258e+215,

  4.97450422247729e+217,

  6.46685548922047e+219,

  8.47158069087882e+221,

  1.118248651196e+224,

  1.48727070609069e+226,

  1.99294274616152e+228,

  2.69047270731805e+230,

  3.65904288195255e+232,

  5.01288874827499e+234,

  6.91778647261949e+236,

  9.61572319694109e+238,

  1.34620124757175e+241,

  1.89814375907617e+243,

  2.69536413788816e+245,

  3.85437071718007e+247,

  5.5502938327393e+249,

  8.04792605747199e+251,

  1.17499720439091e+254,

  1.72724589045464e+256,

  2.55632391787286e+258,

  3.80892263763057e+260,

  5.71338395644585e+262,

  8.62720977423323e+264,

  1.31133588568345e+267,

  2.00634390509568e+269,

  3.08976961384735e+271,

  4.78914290146339e+273,

  7.47106292628289e+275,

  1.17295687942641e+278,

  1.85327186949373e+280,

  2.94670227249504e+282,

  4.71472363599206e+284,

  7.59070505394721e+286,

  1.22969421873945e+289,

  2.0044015765453e+291,

  3.28721858553429e+293,

  5.42391066613159e+295,

  9.00369170577843e+297,

  1.503616514865e+300, // nmaxfactorial = 167

};



/* ----------------------------------------------------------------------

   the function delta given by VMK Eq. 8.2(1)

------------------------------------------------------------------------- */



double ComputeOrientOrderAtom::deltacg(int j1, int j2, int j)

{

  double sfaccg = factorial((j1 + j2 + j) + 1);

  return sqrt(factorial((j1 + j2 - j)) *

              factorial((j1 - j2 + j)) *

              factorial((-j1 + j2 + j)) / sfaccg);

}

src/compute_orientorder_atom.h

+9 −1
Original line number Diff line number Diff line
@@ -33,7 +33,7 @@ class ComputeOrientOrderAtom : public Compute { 
  void compute_peratom();

  double memory_usage();

  double cutsq;

  int iqlcomp, qlcomp, qlcompflag;

  int iqlcomp, qlcomp, qlcompflag, wlflag;

  int *qlist;

  int nqlist;
@@ -55,6 +55,14 @@ class ComputeOrientOrderAtom : public Compute { 


  double polar_prefactor(int, int, double);

  double associated_legendre(int, int, double);



  static const int nmaxfactorial = 167;

  static const double nfac_table[];

  double factorial(int);

  void init_clebsch_gordan();

  double deltacg(int, int, int);

  double *cglist;

  int idxcg_max;

};



}

CRA Git | Maintained and supported by SUSTech CRA and CCSE