Removed unnecessary declarations in compute_born_matrix numdiff part. Commented virial_addon method.

src/EXTRA-COMPUTE/compute_born_matrix.cpp

@@ -301,6 +301,8 @@ void ComputeBornMatrix::init()
     compute_virial = modify->compute[icompute];
 
     // set up reverse index lookup
+    // This table is used for consistency between numdiff and analytical
+    // ordering of the terms.
 
     for (int m = 0; m < nvalues; m++) {
       int a = C_albe[m][0];
@@ -309,47 +311,6 @@ void ComputeBornMatrix::init()
       revalbe[b][a] = m;
     }
 
-    // for (int a = 0; a < NDIR_VIRIAL; a++) {
-    //   for (int b = 0; b < NDIR_VIRIAL; b++) {
-    //     printf("%d ",revalbe[a][b]);
-    //   }
-    //   printf("\n");
-    // }
-
-    // voigt3VtoM notation in normal physics sense,
-    // 3x3 matrix and vector indexing
-    //  i-j:      (1-1), (2-2), (3-3), (2-3), (1-3), (1-2)
-    //  voigt3VtoM: 1      2      3      4      5      6
-
-    voigt3VtoM[0][0]=0;  //  for 1
-    voigt3VtoM[0][1]=0;
-    voigt3VtoM[1][0]=1;  //  for 2
-    voigt3VtoM[1][1]=1;
-    voigt3VtoM[2][0]=2;  //  for 3
-    voigt3VtoM[2][1]=2;
-    voigt3VtoM[3][0]=1;  //  for 4 
-    voigt3VtoM[3][1]=2;
-    voigt3VtoM[4][0]=0;  //  for 5 
-    voigt3VtoM[4][1]=2;
-    voigt3VtoM[5][0]=0;  //  for 6 
-    voigt3VtoM[5][1]=1;
-
-    // to convert to vector indexing:
-    // matrix index to vector index, double -> single index
-    // this is not used at all
-
-    voigt3MtoV[0][0]=0;    voigt3MtoV[0][1]=5;    voigt3MtoV[0][2]=4;
-    voigt3MtoV[1][0]=5;    voigt3MtoV[1][1]=1;    voigt3MtoV[1][2]=3;
-    voigt3MtoV[2][0]=4;    voigt3MtoV[2][1]=3;    voigt3MtoV[2][2]=2;
-
-    // this is just for the virial.
-    // since they use the xx,yy,zz,xy,xz,yz
-    // order not the ordinary voigt
-
-    virialMtoV[0][0]=0;    virialMtoV[0][1]=3;    virialMtoV[0][2]=4;
-    virialMtoV[1][0]=3;    virialMtoV[1][1]=1;    virialMtoV[1][2]=5;
-    virialMtoV[2][0]=4;    virialMtoV[2][1]=5;    virialMtoV[2][2]=2;
-
     // reorder LAMMPS virial vector to Voigt order
 
     virialVtoV[0] = 0;
@@ -359,24 +320,6 @@ void ComputeBornMatrix::init()
     virialVtoV[4] = 4;
     virialVtoV[5] = 3;
 
-    // the following is for 6x6 matrix and vector indexing converter
-    // this is clearly different order form albe[][] and revalbe[]
-    // should not be used
-
-    int indcounter = 0;
-    for(int row = 0; row < NDIR_VIRIAL; row++)
-      for(int col = row; col< NDIR_VIRIAL; col++) {
-        voigt6MtoV[row][col] = voigt6MtoV[col][row] = indcounter;
-        indcounter++;
-      }
-    // printf("Voigt6MtoV:\n");
-    // for (int a = 0; a < NDIR_VIRIAL; a++) {
-    //   for (int b = 0; b < NDIR_VIRIAL; b++) {
-    //     printf("%d ", voigt6MtoV[a][b]);
-    //   }
-    //   printf("\n");
-    // }
-
     // set up 3x3 kronecker deltas
 
     for(int row = 0; row < NXYZ_VIRIAL; row++)
@@ -419,20 +362,12 @@ void ComputeBornMatrix::compute_vector()
 
     MPI_Allreduce(values_local, values_global, nvalues, MPI_DOUBLE, MPI_SUM, world);
 
-    // // convert to pressure units
-    // // As discussed, it might be better to keep it as energy units.
-    // // but this is to be defined
-
-    // double nktv2p = force->nktv2p;
-    // double inv_volume = 1.0 / (domain->xprd * domain->yprd * domain->zprd);
-    // for (int m = 0; m < nvalues; m++) {
-    //   values_global[m] *= (nktv2p * inv_volume);
-    // }
   } else {
 
     // calculate Born matrix using stress finite differences
     compute_numdiff();
 
+    // compute_numdiff output is in pressure units
     // for consistency this is returned in energy units
     double inv_nktv2p = 1.0/force->nktv2p;
     double volume = domain->xprd * domain->yprd * domain->zprd;
@@ -630,18 +565,6 @@ void ComputeBornMatrix::displace_atoms(int nall, int idir, double magnitude)
 {
   double **x = atom->x;
 
-  // A.T.
-  // this works for vector indices  7,  8,  9, 12, 14, 18 and 15, 16, 17
-  // corresponding i,j indices     12, 13, 14, 23, 25, 36 and 26, 34, 35
-  // int k = dirlist[idir][1];
-  // int l = dirlist[idir][0];
-
-  // A.T.
-  // this works for vector indices  7,  8,  9, 12, 14, 18 and 10, 11, 13
-  // corresponding i,j indices     12, 13, 14, 23, 25, 36 and 15, 16, 24
-  // G.C.:
-  // I see no difference with a 0 step simulation between both
-  // methods.
   int k = dirlist[idir][0];
   int l = dirlist[idir][1];
   for (int i = 0; i < nall; i++)
@@ -671,7 +594,6 @@ void ComputeBornMatrix::restore_atoms(int nall, int idir)
 void ComputeBornMatrix::update_virial()
 {
   int eflag = 0;
-  // int vflag = VIRIAL_FDOTR; // Need to generalize this 
   int vflag = 1;
 
   if (force->pair) force->pair->compute(eflag, vflag);
@@ -691,95 +613,49 @@ void ComputeBornMatrix::update_virial()
 /* ----------------------------------------------------------------------
    calculate virial stress addon terms to the Born matrix
    this is based on original code of Dr. Yubao Zhen
-   described here: Comp. Phys. Comm. 183 (2012) 261-265
 ------------------------------------------------------------------------- */
 
 void ComputeBornMatrix::virial_addon()
 {
-  // compute the contribution due to perturbation
-  // here the addon parts are put into born
-  // delta_il sigv_jk +  delta_ik sigv_jl +
-  // delta_jl sigv_ik +  delta_jk sigv_il
-  // Note: in calculation kl is all there from 0 to 6, and ij=(id,jd)
-  //       each time there are six numbers passed for (Dijkl+Djikl)
-  //       and the term I need should be div by 2.
-  //       Job is to arrange the 6 numbers with ij indexing to the 21
-  //       element data structure.
-  // the sigv is the virial stress at current time. It is never touched.
-  // Note the symmetry of (i-j), (k-n), and (ij, kn)
-  // so we only need to evaluate 6x6 matrix with symmetry
-
   int kd, nd, id, jd;
   int m;
 
   double* sigv = compute_virial->vector;
-  double modefactor[6] = {1.0, 1.0, 1.0, 0.5, 0.5, 0.5};
 
-  // Back to the ugly way
+  // This way of doing is not very elegant but is correct.
-  // You can compute these factor by looking at
+  // The complete Cijkl terms are the sum of symmetric terms
-  // every Dijkl terms and adding the proper virials
+  // computed in compute_numdiff and virial stress terms.
-  // Take into account the symmetries. For example:
+  // The viral terms are not symmetric in the tensor computation.
-  // B2323 = s33+D2323; B3232= s22+D3232;
+  // For example:
-  // but D3232=D2323 (computed in compute_numdiff)
+  // C2323 = s33+D2323; C3232 = s22+D3232 etc...
-  // and Cijkl = (Bijkl+Bjikl+Bijlk+Bjilk)/4. = (Bijkl+Bjilk)/2.
+  // However there are two symmetry breaking when reducing
-  // see Yoshimoto eq 15.and eq A3.
+  // the 4-rank tensor to a 2-rank tensor
+  // Cijkl = (Bijkl+Bjikl+Bijlk+Bjilk)/4. = (Bijkl+Bjilk)/2.
+  // and when computing only the 21 independant term.
+  // see Comp. Phys. Comm. 183 (2012) 261–265
+  // and Phys. Rev. B 71, 184108 (2005)
   values_global[0]  += 2.0*sigv[0];
   values_global[1]  += 2.0*sigv[1];
   values_global[2]  += 2.0*sigv[2];
-  // values_global[3]  += 0.5*(sigv[1]+sigv[2]);
-  // values_global[4]  += 0.5*(sigv[0]+sigv[2]);
-  // values_global[5]  += 0.5*(sigv[0]+sigv[1]);
   values_global[3]  += sigv[2];
   values_global[4]  += sigv[2];
   values_global[5]  += sigv[1];
   values_global[6]  += 0.0;
   values_global[7]  += 0.0;
   values_global[8]  += 0.0;
-  //  values_global[9]  += sigv[4];
   values_global[9]  += 2.0*sigv[4];
-  //  values_global[10] += sigv[3];
   values_global[10] += 2.0*sigv[3];
   values_global[11] += 0.0;
-  //  values_global[12] += sigv[5];
   values_global[12] += 2.0*sigv[5];
   values_global[13] += 0.0;
-  //  values_global[14] += sigv[3];
   values_global[14] += 0.0;
-  //  values_global[15] += sigv[5];
   values_global[15] += 0.0;
-  //  values_global[16] += sigv[4];
   values_global[16] += 0.0;
   values_global[17] += 0.0;
   values_global[18] += 0.0;
-  //  values_global[19] += sigv[4];
   values_global[19] += 0.0;
   values_global[20] += sigv[5];
 
-  // This loop is actually bogus.
-  //
-  // for (int idir = 0; idir < NDIR_VIRIAL; idir++) {
-  //   int ijvgt = idir; // this is it.
-  //   double addon;
-
-  //   // extract the two indices composing the voigt representation
-
-  //   id  = voigt3VtoM[ijvgt][0];
-  //   jd  = voigt3VtoM[ijvgt][1];
-
-  //   for (int knvgt=ijvgt; knvgt < NDIR_VIRIAL; knvgt++) {
-  //     kd = voigt3VtoM[knvgt][0];
-  //     nd = voigt3VtoM[knvgt][1];
-  //     addon = kronecker[id][nd]*sigv[virialMtoV[jd][kd]] +
-  //             kronecker[id][kd]*sigv[virialMtoV[jd][nd]];
-  //     if(id != jd)
-  //     addon += kronecker[jd][nd]*sigv[virialMtoV[id][kd]] +
-  //              kronecker[jd][kd]*sigv[virialMtoV[id][nd]];
-
-  //     m = revalbe[ijvgt][knvgt];
-
-  //     values_global[revalbe[ijvgt][knvgt]] += 0.5*modefactor[idir]*addon;
-  //   }
-  // }
 }
 
 /* ----------------------------------------------------------------------

src/EXTRA-COMPUTE/compute_born_matrix.h

@@ -72,11 +72,7 @@ namespace LAMMPS_NS {
      static constexpr int NDIR_VIRIAL = 6;    // dimension of virial and strain vectors
      static constexpr int NXYZ_VIRIAL = 3;    // number of Cartesian coordinates
      int revalbe[NDIR_VIRIAL][NDIR_VIRIAL];
-     int voigt3VtoM[NDIR_VIRIAL][2];
-     int voigt3MtoV[NXYZ_VIRIAL][NXYZ_VIRIAL];
-     int virialMtoV[NXYZ_VIRIAL][NXYZ_VIRIAL];
      int virialVtoV[NDIR_VIRIAL];
-     int voigt6MtoV[NDIR_VIRIAL][NDIR_VIRIAL];
      int kronecker[NXYZ_VIRIAL][NXYZ_VIRIAL];
      double **temp_x;                // original coords
      double **temp_f;                // original forces