Chill

Data Structures
struct	chill::YlmAtom
	This contains a complex vector of length $2l+1$. More...
struct	chill::QlmAtom
	This is the local orientational bond order parameter $q_{lm}$, of length $2l+1$. More...

Functions
molSys::PointCloud< molSys::Point< double >, double >	chill::getCorrel (molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int >> nList, bool isSlice=false)
molSys::PointCloud< molSys::Point< double >, double >	chill::getIceTypeNoPrint (molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int >> nList, bool isSlice=false)
	Classifies each atom according to the CHILL algorithm without printing. More...
molSys::PointCloud< molSys::Point< double >, double >	chill::getIceType (molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int >> nList, std::string path, int firstFrame, bool isSlice=false, std::string outputFileName="chill.txt")
	Classifies each atom according to the CHILL algorithm. More...
molSys::PointCloud< molSys::Point< double >, double >	chill::getCorrelPlus (molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int >> nList, bool isSlice=false)
	Gets c_ij and then classifies bond types according to the CHILL+ algorithm. More...
molSys::PointCloud< molSys::Point< double >, double >	chill::getIceTypePlus (molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int >> nList, std::string path, int firstFrame, bool isSlice=false, std::string outputFileName="chillPlus.txt")
	Classifies each atom according to the CHILL+ algorithm. More...
std::vector< double >	chill::getq6 (molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int >> nList, bool isSlice=false)
molSys::PointCloud< molSys::Point< double >, double >	chill::reclassifyWater (molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< double > *q6)
int	chill::printIceType (molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::string path, int firstFrame, bool isSlice=false, std::string outputFileName="superChill.txt")
	Prints out the iceType for a particular frame onto the terminal. More...
bool	chill::isInterfacial (molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int >> nList, int iatom, int num_staggrd, int num_eclipsd)
int	chill::numStaggered (molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int >> nList, int jatom)
	Finds the number of staggered bonds for a given atom of index jatom. More...
std::vector< std::complex< double > >	sph::spheriHarmo (int orderL, std::array< double, 2 > radialCoord)
std::array< double, 2 >	sph::radialCoord (std::array< double, 3 > cartCoord)
std::vector< std::complex< double > >	sph::lookupTableQ3Vec (std::array< double, 2 > angles)
	Lookup table for Q3. More...
std::complex< double >	sph::lookupTableQ3 (int m, std::array< double, 2 > angles)
	Lookup table for Q3 (m=0 to m=6) More...
std::vector< std::complex< double > >	sph::lookupTableQ6Vec (std::array< double, 2 > angles)
	Lookup table for Q6. More...
std::complex< double >	sph::lookupTableQ6 (int m, std::array< double, 2 > angles)
	Lookup table for Q6 (m=0 to m=12) More...

Variables
std::vector< std::complex< double > >	chill::YlmAtom::ylm
std::vector< YlmAtom >	chill::QlmAtom::ptq

Detailed Description

Function Documentation

getCorrel()

molSys::PointCloud< molSys::Point< double >, double > chill::getCorrel	(	molSys::PointCloud< molSys::Point< double >, double > *	yCloud,
		std::vector< std::vector< int >>	nList,
		bool	isSlice = false
	)

Function for getting the bond order correlations $c_{ij}$ (or $a_{ij}$ in some treatments) according to the CHILL algorithm.

Parameters

[in,out]	yCloud	The output molSys::PointCloud
[in]	nList	The row-ordered neighbour list, by ID. The first element of each row is the particle ID, followed by the IDs of the neighbours
[in]	isSlice	This decides whether there is a slice or not

Uses Boost for spherical harmonics, and gets c_ij according to the CHILL algorithm

Definition at line 299 of file bop.cpp.

                                                                {
  //
  int l = 3;      // TODO: Don't hard-code this; change later
  int iatomID;    // Atom ID (key) of iatom
  int iatomIndex; // Index (value) of iatom
  int jatomID;    // Atom ID (key) of jatom
  int jatomIndex; // Index (value) of nearest neighbour
  std::array<double, 3> delta;
  std::array<double, 2> angles;
  chill::QlmAtom QlmTotal; // Qlm for each iatom
  std::vector<std::complex<double>>
      yl; // temp q_lm for each pair of iatom and jatom
  std::complex<double> dot_product = {0, 0};
  std::complex<double> qI = {0, 0};
  std::complex<double> qJ = {0, 0};
  std::complex<double> Inorm = {0, 0};
  std::complex<double> Jnorm = {0, 0};
  std::complex<double> complexDenominator = {0, 0};
  std::complex<double> complexCij = {0, 0};
  molSys::Result temp_cij; // Holds the c_ij value
  double cij_real;
  int nnumNeighbours; // Number of nearest neighbours for iatom
 
  QlmTotal.ptq.resize(yCloud->nop);
 
  // Loop through the neighbour list
  for (int iatom = 0; iatom < nList.size(); iatom++) {
    // The atom index is iatom
    iatomIndex = iatom;
    iatomID =
        nList[iatomIndex][0]; // The first element in nList is the ID of iatom
    nnumNeighbours = nList[iatomIndex].size() - 1;
    // Now loop over the first four neighbours
    for (int j = 1; j <= nnumNeighbours; j++) {
      // Get the ID of jatom saved in the neighbour list
      jatomID = nList[iatomIndex][j];
 
      // Get the index of jatom
      auto it = yCloud->idIndexMap.find(jatomID);
 
      if (it != yCloud->idIndexMap.end()) {
        jatomIndex = it->second;
      } // found jatom
      else {
        std::cerr << "Something is wrong with the ID and index map.\n";
        return *yCloud;
      } // error handling
 
      // Get the relative distance now that the index values are known
      delta = gen::relDist(yCloud, iatomIndex, jatomIndex);
      double r = std::sqrt(std::pow(delta[0], 2.0) + std::pow(delta[1], 2.0) +
                           std::pow(delta[2], 2.0));
      angles[1] = acos(delta[2] / r);        // theta
      angles[0] = atan2(delta[0], delta[1]); // phi
 
      // Now add over all nearest neighbours
      if (j == 1) {
        QlmTotal.ptq[iatomIndex].ylm = sph::spheriHarmo(3, angles);
        // QlmTotal.ptq[iatom].ylm = sph::lookupTableQ3Vec(angles);
        continue;
      }
      // Not for the first jatom
      yl = sph::spheriHarmo(3, angles);
      for (int m = 0; m < 2 * l + 1; m++) {
        QlmTotal.ptq[iatomIndex].ylm[m] += yl[m];
        // QlmTotal.ptq[iatom].ylm[m] += sph::lookupTableQ3(m, angles);
      }
    } // End of loop over 4 nearest neighbours
 
    // Divide by 4
    QlmTotal.ptq[iatomIndex].ylm =
        gen::avgVector(QlmTotal.ptq[iatom].ylm, l, nnumNeighbours);
  } // End of looping over all iatom
 
  // ------------------------------------------------
  // Now that you have all qlm for the particles,
  // find c_ij
  for (int iatom = 0; iatom < yCloud->nop; iatom++) {
    // if(yCloud->pts[iatom].type!=typeO){continue;}
    // if this is a slice and the particle is not in the slice
    // then skip
    if (isSlice) {
      if (yCloud->pts[iatom].inSlice == false) {
        continue;
      }
    }
    // The index is what we are looping through
    iatomIndex = iatom;
    nnumNeighbours = nList[iatomIndex].size() - 1;
    yCloud->pts[iatomIndex].c_ij.reserve(nnumNeighbours);
    // loop over the 4 nearest neighbours
    for (int j = 1; j <= nnumNeighbours; j++) {
      // Init to zero
      dot_product = {0, 0};
      Inorm = {0, 0};
      Jnorm = {0, 0};
      // Get ID of the nearest neighbour
      jatomID = nList[iatomIndex][j];
      // Get the index (value) from the ID (key)
      auto it = yCloud->idIndexMap.find(jatomID);
 
      if (it != yCloud->idIndexMap.end()) {
        jatomIndex = it->second;
      } // found jatom
      else {
        std::cerr << "Something is wrong with the ID and index map.\n";
        return *yCloud;
      } // error handling
      // Spherical harmonics
      for (int m = 0; m < 2 * l + 1; m++) {
        qI = QlmTotal.ptq[iatomIndex].ylm[m];
        qJ = QlmTotal.ptq[jatomIndex].ylm[m];
        dot_product = dot_product + (qI * std::conj(qJ)); // unnormalized
        Inorm = Inorm + (qI * std::conj(qI));
        Jnorm = Jnorm + (qJ * std::conj(qJ));
      } // end loop over m components
      // Get the denominator
      complexDenominator = std::sqrt(Inorm * Jnorm);
      complexCij = dot_product / complexDenominator;
      // Update c_ij and type
      cij_real = complexCij.real();
      temp_cij.c_value = cij_real;
      if (cij_real < -0.8) {
        temp_cij.classifier = molSys::bond_type::staggered;
      } else if (cij_real > -0.2 && cij_real < -0.05) {
        temp_cij.classifier = molSys::bond_type::eclipsed;
      } else {
        temp_cij.classifier = molSys::bond_type::out_of_range;
      }
      yCloud->pts[iatomIndex].c_ij.push_back(temp_cij);
    } // end loop over nearest neighbours
  }
 
  return *yCloud;
}

Code

getCorrelPlus()

molSys::PointCloud< molSys::Point< double >, double > chill::getCorrelPlus	(	molSys::PointCloud< molSys::Point< double >, double > *	yCloud,
		std::vector< std::vector< int >>	nList,
		bool	isSlice = false
	)

Gets c_ij and then classifies bond types according to the CHILL+ algorithm.

Function for getting the bond order correlations $c_{ij}$ (alternatively $a_{ij}$ in certain texts) using the CHILL+ algorithm

Parameters

[in,out]	yCloud	The output molSys::PointCloud
[in]	nList	Row-ordered neighbour list by atom ID
[in]	isSlice	This decides whether there is a slice or not

Definition at line 602 of file bop.cpp.

                                                                    {
  //
  int l = 3;      // TODO: Don't hard-code this; change later
  int iatomID;    // Atom ID (key) of iatom
  int iatomIndex; // Index (value) of iatom
  int jatomID;    // Atom ID (key) of jatom
  int jatomIndex; // Index (value) of nearest neighbour
  std::array<double, 3> delta;
  std::array<double, 2> angles;
  chill::QlmAtom QlmTotal; // Qlm for each iatom
  std::vector<std::complex<double>>
      yl; // temp q_lm for each pair of iatom and jatom
  std::complex<double> dot_product = {0, 0};
  std::complex<double> qI = {0, 0};
  std::complex<double> qJ = {0, 0};
  std::complex<double> Inorm = {0, 0};
  std::complex<double> Jnorm = {0, 0};
  std::complex<double> complexDenominator = {0, 0};
  std::complex<double> complexCij = {0, 0};
  molSys::Result temp_cij; // Holds the c_ij value
  double cij_real;
  int nnumNeighbours; // Number of nearest neighbours for iatom
 
  QlmTotal.ptq.resize(yCloud->nop);
 
  // Loop through the neighbour list
  for (int iatom = 0; iatom < nList.size(); iatom++) {
    // The atom index is iatom
    iatomIndex = iatom;
    iatomID =
        nList[iatomIndex][0]; // The first element in nList is the ID of iatom
    nnumNeighbours = nList[iatomIndex].size() - 1;
    // Now loop over the first four neighbours
    for (int j = 1; j <= nnumNeighbours; j++) {
      // Get the ID of jatom saved in the neighbour list
      jatomID = nList[iatomIndex][j];
 
      // Get the index of jatom
      auto it = yCloud->idIndexMap.find(jatomID);
 
      if (it != yCloud->idIndexMap.end()) {
        jatomIndex = it->second;
      } // found jatom
      else {
        std::cerr << "Something is wrong with the ID and index map.\n";
        return *yCloud;
      } // error handling
 
      // Get the relative distance now that the index values are known
      delta = gen::relDist(yCloud, iatomIndex, jatomIndex);
      double r = std::sqrt(std::pow(delta[0], 2.0) + std::pow(delta[1], 2.0) +
                           std::pow(delta[2], 2.0));
      angles[1] = acos(delta[2] / r);        // theta
      angles[0] = atan2(delta[0], delta[1]); // phi
 
      // Now add over all nearest neighbours
      if (j == 1) {
        QlmTotal.ptq[iatomIndex].ylm = sph::spheriHarmo(3, angles);
        // QlmTotal.ptq[iatom].ylm = sph::lookupTableQ3Vec(angles);
        continue;
      }
      // Not for the first jatom
      yl = sph::spheriHarmo(3, angles);
      for (int m = 0; m < 2 * l + 1; m++) {
        QlmTotal.ptq[iatomIndex].ylm[m] += yl[m];
        // QlmTotal.ptq[iatom].ylm[m] += sph::lookupTableQ3(m, angles);
      }
    } // End of loop over 4 nearest neighbours
 
    // Divide by 4
    QlmTotal.ptq[iatomIndex].ylm =
        gen::avgVector(QlmTotal.ptq[iatom].ylm, l, nnumNeighbours);
  } // End of looping over all iatom
 
  // ------------------------------------------------
  // Now that you have all qlm for the particles,
  // find c_ij
  for (int iatom = 0; iatom < yCloud->nop; iatom++) {
    // if(yCloud->pts[iatom].type!=typeO){continue;}
    // if this is a slice and the particle is not in the slice
    // then skip
    if (isSlice) {
      if (yCloud->pts[iatom].inSlice == false) {
        continue;
      }
    }
    // The index is what we are looping through
    iatomIndex = iatom;
    nnumNeighbours = nList[iatomIndex].size() - 1;
    yCloud->pts[iatomIndex].c_ij.reserve(nnumNeighbours);
    // loop over the 4 nearest neighbours
    for (int j = 1; j <= nnumNeighbours; j++) {
      // Init to zero
      dot_product = {0, 0};
      Inorm = {0, 0};
      Jnorm = {0, 0};
      // Get ID of the nearest neighbour
      jatomID = nList[iatomIndex][j];
      // Get the index (value) from the ID (key)
      auto it = yCloud->idIndexMap.find(jatomID);
 
      if (it != yCloud->idIndexMap.end()) {
        jatomIndex = it->second;
      } // found jatom
      else {
        std::cerr << "Something is wrong with the ID and index map.\n";
        return *yCloud;
      } // error handling
      // Spherical harmonics
      for (int m = 0; m < 2 * l + 1; m++) {
        qI = QlmTotal.ptq[iatomIndex].ylm[m];
        qJ = QlmTotal.ptq[jatomIndex].ylm[m];
        dot_product = dot_product + (qI * std::conj(qJ)); // unnormalized
        Inorm = Inorm + (qI * std::conj(qI));
        Jnorm = Jnorm + (qJ * std::conj(qJ));
      } // end loop over m components
      // Get the denominator
      complexDenominator = std::sqrt(Inorm * Jnorm);
      complexCij = dot_product / complexDenominator;
      // Update c_ij and type
      cij_real = complexCij.real();
      temp_cij.c_value = cij_real;
      if (cij_real <= -0.8) {
        temp_cij.classifier = molSys::bond_type::staggered;
      } else if (cij_real >= -0.35 && cij_real <= 0.25) {
        temp_cij.classifier = molSys::bond_type::eclipsed;
      } else {
        temp_cij.classifier = molSys::bond_type::out_of_range;
      }
      yCloud->pts[iatomIndex].c_ij.push_back(temp_cij);
    } // end loop over nearest neighbours
  }
 
  // ------------------------------------------------
 
  return *yCloud;
}

Code

getIceType()

molSys::PointCloud< molSys::Point< double >, double > chill::getIceType	(	molSys::PointCloud< molSys::Point< double >, double > *	yCloud,
		std::vector< std::vector< int >>	nList,
		std::string	path,
		int	firstFrame,
		bool	isSlice = false,
		std::string	outputFileName = "chill.txt"
	)

Classifies each atom according to the CHILL algorithm.

Function that classifies every particle's molSys::atom_state_type ice type, according to the CHILL algorithm.

Parameters

[in,out]	yCloud	The output molSys::PointCloud
[in]	nList	Row-ordered neighbour list by atom ID
[in]	path	Path to the output directory to which ice types are written out to
[in]	firstFrame	First frame to be analyzed
[in]	isSlice	This decides whether there is a slice or not
[in]	outputFileName	Name of the output file, to which the ice types will be written out.

Definition at line 505 of file bop.cpp.

                                                                          {
  int ih, ic, water, interIce, unknown, total; // No. of particles of each type
  ih = ic = water = unknown = interIce = total = 0;
  int num_staggrd, num_eclipsd, na;
  molSys::bond_type bondType;
  int nnumNeighbours; // Number of nearest neighbours
 
  for (int iatom = 0; iatom < yCloud->nop; iatom++) {
    // if(yCloud->pts[iatom].type!=typeO){continue;}
    // if this is a slice and the particle is not in the slice
    // then skip
    if (isSlice) {
      if (yCloud->pts[iatom].inSlice == false) {
        continue;
      }
    }
    total++; // Update the total number of atoms considered. Change this to
    // check for slices
    num_staggrd = num_eclipsd = na =
        0; // init to zero before loop through neighbours
 
    nnumNeighbours = nList[iatom].size() - 1;
    // Loop through the bond cij and get the number of staggered, eclipsed bonds
    for (int j = 0; j < nnumNeighbours; j++) {
      bondType = yCloud->pts[iatom].c_ij[j].classifier;
      if (bondType == molSys::bond_type::eclipsed) {
        num_eclipsd++;
      } else if (bondType == molSys::bond_type::staggered) {
        num_staggrd++;
      } else {
        na++;
      }
    } // End of loop through neighbours
 
    // Add more tests later
    yCloud->pts[iatom].iceType = molSys::atom_state_type::unclassified; // default
    // Cubic ice
    // if (num_eclipsd==0 && num_staggrd==4){
    //  yCloud->pts[iatom].iceType = molSys::cubic;
    //  ic++;
    // }
    if (num_staggrd >= 4) {
      yCloud->pts[iatom].iceType = molSys::atom_state_type::cubic;
      ic++;
    }
    // Hexagonal
    else if (num_eclipsd == 1 && num_staggrd == 3) {
      yCloud->pts[iatom].iceType = molSys::atom_state_type::hexagonal;
      ih++;
    }
    // Interfacial
    else if (isInterfacial(yCloud, nList, iatom, num_staggrd, num_eclipsd)) {
      yCloud->pts[iatom].iceType = molSys::atom_state_type::interfacial;
      interIce++;
    } else {
      yCloud->pts[iatom].iceType = molSys::atom_state_type::water;
      water++;
    }
 
  } // End of loop through every iatom
 
  // water = total - ic -ih;
 
  // --------------------
  // Create the directories if needed
  sout::makePath(path);
  std::string outputDirName = path + "bop";
  sout::makePath(outputDirName);
  // --------------------
 
  // Print to file
  std::ofstream outputFile;
  outputFile.open(path + "bop/" + outputFileName, std::ios_base::app);
  // --------------------
  // Write out the comment line for the first frame
  if (yCloud->currentFrame == firstFrame) {
    outputFile << "Frame Ic Ih Interfacial Water Total \n";
  }
  // --------------------
  outputFile << yCloud->currentFrame << " " << ic << " " << ih << " "
             << interIce << " " << water << " " << total << "\n";
  outputFile.close();
 
  return *yCloud;
}

Code

getIceTypeNoPrint()

molSys::PointCloud< molSys::Point< double >, double > chill::getIceTypeNoPrint	(	molSys::PointCloud< molSys::Point< double >, double > *	yCloud,
		std::vector< std::vector< int >>	nList,
		bool	isSlice = false
	)

Classifies each atom according to the CHILL algorithm without printing.

Function that classifies every particle's molSys::atom_state_type ice type, according to the CHILL algorithm. Does not print out the information.

Parameters

[in,out]	yCloud	The output molSys::PointCloud
[in]	isSlice	This decides whether there is a slice or not
[in]	nList	Row-ordered neighbour list by atom ID

Definition at line 437 of file bop.cpp.

                                                   {
  int ih, ic, water, interIce, unknown, total; // No. of particles of each type
  ih = ic = water = unknown = interIce = total = 0;
  int num_staggrd, num_eclipsd, na;
  molSys::bond_type bondType;
  int nnumNeighbours; // Number of nearest neighbours
 
  for (int iatom = 0; iatom < yCloud->nop; iatom++) {
    // if(yCloud->pts[iatom].type!=typeO){continue;}
    // if this is a slice and the particle is not in the slice
    // then skip
    if (isSlice) {
      if (yCloud->pts[iatom].inSlice == false) {
        continue;
      }
    }
    total++; // Update the total number of atoms considered. Change this to
    // check for slices
    num_staggrd = num_eclipsd = na =
        0; // init to zero before loop through neighbours
 
    nnumNeighbours = nList[iatom].size() - 1;
    // Loop through the bond cij and get the number of staggered, eclipsed bonds
    for (int j = 0; j < nnumNeighbours; j++) {
      bondType = yCloud->pts[iatom].c_ij[j].classifier;
      if (bondType == molSys::bond_type::eclipsed) {
        num_eclipsd++;
      } else if (bondType == molSys::bond_type::staggered) {
        num_staggrd++;
      } else {
        na++;
      }
    } // End of loop through neighbours
 
    // Add more tests later
    yCloud->pts[iatom].iceType = molSys::atom_state_type::unclassified; // default
    // Cubic ice
    // if (num_eclipsd==0 && num_staggrd==4){
    //  yCloud->pts[iatom].iceType = molSys::cubic;
    //  ic++;
    // }
    if (num_staggrd >= 4) {
      yCloud->pts[iatom].iceType = molSys::atom_state_type::cubic;
      ic++;
    }
    // Hexagonal
    else if (num_eclipsd == 1 && num_staggrd == 3) {
      yCloud->pts[iatom].iceType = molSys::atom_state_type::hexagonal;
      ih++;
    }
    // Interfacial
    else if (isInterfacial(yCloud, nList, iatom, num_staggrd, num_eclipsd)) {
      yCloud->pts[iatom].iceType = molSys::atom_state_type::interfacial;
      interIce++;
    } else {
      yCloud->pts[iatom].iceType = molSys::atom_state_type::water;
      water++;
    }
 
  } // End of loop through every iatom
 
  return *yCloud;
}

Code

getIceTypePlus()

molSys::PointCloud< molSys::Point< double >, double > chill::getIceTypePlus	(	molSys::PointCloud< molSys::Point< double >, double > *	yCloud,
		std::vector< std::vector< int >>	nList,
		std::string	path,
		int	firstFrame,
		bool	isSlice = false,
		std::string	outputFileName = "chillPlus.txt"
	)

Classifies each atom according to the CHILL+ algorithm.

Function that classifies the molSys::atom_state_type ice type of each particle, according to the CHILL+ algorithm.

Parameters

[in,out]	yCloud	The output molSys::PointCloud
[in]	nList	Row-ordered neighbour list by atom ID
[in]	path	Path to the output directory to which ice types are written out to
[in]	firstFrame	The first frame to be analyzed
[in]	isSlice	This decides whether there is a slice or not
[in]	outputFileName	Name of the output file, to which the ice types will be written out. The default file name is "chillPlus.txt"

Definition at line 755 of file bop.cpp.

                                                {
  int ih, ic, interIce, water, unknown, clath, interClath,
      total; // No. of particles of each type
  ih = ic = water = unknown = interIce = total = 0;
  clath = interClath = 0;
  int num_staggrd, num_eclipsd, na;
  molSys::bond_type bondType;
  int nnumNeighbours; // number of nearest neighbours
 
  for (int iatom = 0; iatom < yCloud->nop; iatom++) {
    // if(yCloud->pts[iatom].type!=typeO){continue;}
    // if this is a slice and the particle is not in the slice
    // then skip
    if (isSlice) {
      if (yCloud->pts[iatom].inSlice == false) {
        continue;
      }
    }
    total++; // Update the total number of atoms considered. Change this to a
             // check for slices
    nnumNeighbours =
        nList[iatom].size() - 1; // number of nearest neighbours (total -1)
    num_staggrd = num_eclipsd = na =
        0; // init to zero before loop through neighbours
    // Loop through the bond cij and get the number of staggered, eclipsed bonds
    for (int j = 0; j < nnumNeighbours; j++) {
      bondType = yCloud->pts[iatom].c_ij[j].classifier;
      if (bondType == molSys::bond_type::eclipsed) {
        num_eclipsd++;
      } else if (bondType == molSys::bond_type::staggered) {
        num_staggrd++;
      } else {
        na++;
      }
    } // End of loop through neighbours
 
    // Add more tests later
    yCloud->pts[iatom].iceType = molSys::atom_state_type::unclassified; // default
    if (nnumNeighbours == 4) {
      // Cubic ice
      if (num_eclipsd == 0 && num_staggrd == 4) {
        yCloud->pts[iatom].iceType = molSys::atom_state_type::cubic;
        ic++;
      }
      // Hexagonal
      else if (num_eclipsd == 1 && num_staggrd == 3) {
        yCloud->pts[iatom].iceType = molSys::atom_state_type::hexagonal;
        ih++;
      }
      // Interfacial
      else if (isInterfacial(yCloud, nList, iatom, num_staggrd, num_eclipsd)) {
        yCloud->pts[iatom].iceType = molSys::atom_state_type::interfacial;
        interIce++;
      }
      // Clathrate
      else if (num_eclipsd == 4 && num_staggrd == 0) {
        yCloud->pts[iatom].iceType = molSys::atom_state_type::clathrate;
        clath++;
      }
      // Interfacial clathrate
      else if (num_eclipsd == 3) {
        yCloud->pts[iatom].iceType = molSys::atom_state_type::interClathrate;
        interClath++;
      }
      // Water
      else {
        yCloud->pts[iatom].iceType = molSys::atom_state_type::water;
        water++;
      }
    } else {
      yCloud->pts[iatom].iceType = molSys::atom_state_type::water;
      water++;
    }
 
  } // End of loop through every iatom
 
  // water = total - ic -ih;
 
  // --------------------
  // Create the directories if needed
  sout::makePath(path);
  std::string outputDirName = path + "bop";
  sout::makePath(outputDirName);
  // --------------------
 
  std::ofstream outputFile;
  outputFile.open(path + "bop/" + outputFileName, std::ios_base::app);
  // --------------------
  // Comment line for the first line
  if (yCloud->currentFrame == firstFrame) {
    outputFile << "Frame Ic Ih Interfacial Clath InterClath Water Total\n";
  }
  // --------------------
  outputFile << yCloud->currentFrame << " " << ic << " " << ih << " "
             << interIce << " " << clath << " " << interClath << " " << water
             << " " << total << "\n";
  outputFile.close();
 
  return *yCloud;
}

Code

getq6()

std::vector< double > chill::getq6	(	molSys::PointCloud< molSys::Point< double >, double > *	yCloud,
		std::vector< std::vector< int >>	nList,
		bool	isSlice = false
	)

q6 can distinguish between water and ice. Use this for the largest ice cluster

Function for getting the averaged $q_6$ parameter.

Parameters

[in,out]	yCloud	The output molSys::PointCloud
[in]	nList	Row-ordered neighbour list by atom ID
[in]	isSlice	This decides whether there is a slice or not

Returns

a double vector of the averaged $q_6$ values.

Definition at line 868 of file bop.cpp.

                                                            {
  //
  int l = 6;      // We're using q6 here
  int jatomID;    // Atom ID of the nearest neighbour
  int jatomIndex; // Index of nearest neighbour
  std::array<double, 3> delta;
  std::array<double, 2> angles;
  chill::QlmAtom QlmTotal; // Qlm for each iatom
  std::vector<std::complex<double>>
      yl; // temp q_lm for each pair of iatom and jatom
  std::complex<double> dot_product = {0, 0};
  std::complex<double> qI = {0, 0};
  std::complex<double> qJ = {0, 0};
  std::complex<double> Inorm = {0, 0};
  std::complex<double> Jnorm = {0, 0};
  std::complex<double> complexDenominator = {0, 0};
  std::complex<double> complexCij = {0, 0};
  double cij_real;
  int nnumNeighbours;
  std::vector<double> resultQ; // Vector with averaged q values
  double q_value = 0.0;        // Averaged q value per neighbour pair
 
  QlmTotal.ptq.resize(yCloud->nop);
  resultQ.resize(yCloud->nop);
 
  // Loop through every index in yCloud
  for (int iatom = 0; iatom < yCloud->nop; iatom++) {
    // if(yCloud->pts[iatom].type!=typeO){continue;}
 
    nnumNeighbours = nList[iatom].size() - 1; // One less than the actual length
    // Now loop over the first four neighbours
    for (int j = 1; j <= nnumNeighbours; j++) {
      // Get the atom ID
      jatomID = nList[iatom][j]; // Atom ID (key)
 
      // Get the atom index (value) from the atom ID (key)
      auto it = yCloud->idIndexMap.find(jatomID);
 
      if (it != yCloud->idIndexMap.end()) {
        jatomIndex = it->second;
      } else {
        std::cerr << "Your map must be wrong.\n";
        return resultQ; // return with error
      }
 
      // Get the relative distances
      delta = gen::relDist(yCloud, iatom, jatomIndex);
 
      // angles = sph::radialCoord(delta);
      double r = std::sqrt(std::pow(delta[0], 2.0) + std::pow(delta[1], 2.0) +
                           std::pow(delta[2], 2.0));
      angles[1] = acos(delta[2] / r);        // theta
      angles[0] = atan2(delta[0], delta[1]); // phi
 
      // Now add over all nearest neighbours
      if (j == 1) {
        // QlmTotal.ptq[iatom].ylm = sph::spheriHarmo(l, angles);
        QlmTotal.ptq[iatom].ylm = sph::lookupTableQ6Vec(angles);
        continue;
      }
      // Not for the first jatom
      yl = sph::spheriHarmo(l, angles);
      for (int m = 0; m < 2 * l + 1; m++) {
        QlmTotal.ptq[iatom].ylm[m] += yl[m];
        // QlmTotal.ptq[iatom].ylm[m] += sph::lookupTableQ6(m, angles);
      }
    } // End of loop over 4 nearest neighbours
 
    // Divide by 4
    QlmTotal.ptq[iatom].ylm =
        gen::avgVector(QlmTotal.ptq[iatom].ylm, l, nnumNeighbours);
  } // End of looping over all iatom
 
  // ------------------------------------------------
  // Now that you have all qlm for the particles,
  // find c_ij
  for (int iatom = 0; iatom < yCloud->nop; iatom++) {
    // if(yCloud->pts[iatom].type!=typeO){continue;}
    // if this is a slice and the particle is not in the slice
    // then skip TODO: UNCOMMENT
    // if(isSlice){
    //  if(yCloud->pts[iatom].inSlice==false){continue;}
    // }
 
    nnumNeighbours = nList[iatom].size() - 1; // Number of nearest neighbours
    q_value = 0.0;                            // initialize to zero
    // yCloud->pts[iatom].c_ij.reserve(nnumNeighbours);
    // loop over the 4 nearest neighbours
    for (int j = 1; j <= nnumNeighbours; j++) {
      // Init to zero
      dot_product = {0, 0};
      Inorm = {0, 0};
      Jnorm = {0, 0};
      // Get index of the nearest neighbour!
      jatomID = nList[iatom][j]; // Atom ID (key)
 
      // Get the index jatomIndex
      auto it = yCloud->idIndexMap.find(jatomID);
 
      if (it != yCloud->idIndexMap.end()) {
        jatomIndex = it->second;
      } // end of getting the index of jatom
 
      for (int m = 0; m < 2 * l + 1; m++) {
        qI = QlmTotal.ptq[iatom].ylm[m];
        qJ = QlmTotal.ptq[jatomIndex].ylm[m];
        dot_product = dot_product + (qI * std::conj(qJ)); // unnormalized
        Inorm = Inorm + (qI * std::conj(qI));
        Jnorm = Jnorm + (qJ * std::conj(qJ));
      } // end loop over m components
      // Get the denominator
      complexDenominator = std::sqrt(Inorm * Jnorm);
      complexCij = dot_product / complexDenominator;
      // Update c_ij and type
      cij_real = complexCij.real();
 
      q_value += cij_real;
 
    } // end loop over nearest neighbours
 
    // Average q_value over all nearest neighbours
    q_value /= (double)nnumNeighbours;
 
    resultQ[iatom] = q_value; // Update the vector of averaged q6
  }
 
  // ------------------------------------------------
 
  return resultQ;
}

Code

isInterfacial()

bool chill::isInterfacial	(	molSys::PointCloud< molSys::Point< double >, double > *	yCloud,
		std::vector< std::vector< int >>	nList,
		int	iatom,
		int	num_staggrd,
		int	num_eclipsd
	)

Checks if a given iatom is interfacial ice or not, according to the CHILL algorithm

Function that checks if the particle with the given atom index is interfacial or not.

Parameters

[in]	yCloud	The input molSys::PointCloud
[in]	nList	Row-ordered neighbour list by atom ID
[in]	iatom	The vector index of the current particle
[in]	num_staggrd	The number of staggered bonds that the current particle participates in
[in]	num_eclipsd	The number of eclipsed bonds that the current particle participates in

Returns

a bool; true if the particle is interfacial and otherwise false

Definition at line 1129 of file bop.cpp.

                     {
  int nnumNeighbours =
      nList[iatom].size() - 1; // number of nearest neighbours of iatom
  int neighStaggered =
      0;          // number of staggered bonds in the neighbours of iatom
  int jatomID;    // ID of the nearest neighbour
  int jatomIndex; // Index (value) of nearest neighbour
 
  // INTERFACIAL
  // Condition 1 : only two staggered bonds and at least
  // one neighbor with more than two staggered bonds
  if (num_staggrd == 2) {
    // Loop over the nearest neighbours
    for (int j = 1; j <= nnumNeighbours; j++) {
      // Get index of the nearest neighbour
      jatomID = nList[iatom][j];
      // Get the index (value) from the ID (key)
      auto it = yCloud->idIndexMap.find(jatomID);
 
      if (it != yCloud->idIndexMap.end()) {
        jatomIndex = it->second;
      } else {
        std::cerr << "Something is gravely wrong with your map.\n";
        return false;
      }
      //
      neighStaggered = chill::numStaggered(yCloud, nList, jatomIndex);
      if (neighStaggered > 2) {
        return true;
      }
    } // End loop over nearest neighbours
  }   // end condition 1
  // Condition 2 : three staggered bonds, no eclipsed bond,
  // and at least one neighbor with two staggered bonds
  if (num_staggrd == 3 && num_eclipsd == 0) {
    // Loop over the nearest neighbours
    for (int j = 1; j <= nnumNeighbours; j++) {
      // Get index of the nearest neighbour
      // ID of the nearest neighbour
      jatomID = nList[iatom][j];
      // Get the index (value) from the ID (key)
      auto it = yCloud->idIndexMap.find(jatomID);
 
      if (it != yCloud->idIndexMap.end()) {
        jatomIndex = it->second;
      } else {
        std::cerr << "Something is gravely wrong with your map.\n";
        return false;
      }
      //
      neighStaggered = chill::numStaggered(yCloud, nList, jatomIndex);
      if (neighStaggered == 2) {
        return true;
      }
    }
  }
  // not interfacial
  return false;
}

Code

lookupTableQ3()

std::complex< double > sph::lookupTableQ3	(	int	m,
		std::array< double, 2 >	angles
	)

Lookup table for Q3 (m=0 to m=6)

Look-up hard-coded values for $Q_3$

It is recommended to use the Boost version of this function, sph::spheriHarmo, instead.

Parameters

[in]	m	An int such that $-3<=m<=3$
[in]	angles	The azimuth and polar angles for a particular particle

Returns

a complex vector, of length $7$, calculated using hard-coded values

Definition at line 133 of file bop.cpp.

                                                                       {
  std::complex<double> result(0.0, 0.0);
  const double pi = std::acos(-1);
  const std::complex<double> i(0.0, 1.0);
  double constant;
  double theta = angles[1];
  double phi = angles[0];
 
  if (m == 0) {
    constant = 0.125 * std::sqrt(35 / pi);
    result = constant * std::exp(-3.0 * i * phi) * std::pow(sin(theta), 3.0);
    return result;
  } else if (m == 1) {
    constant = 0.25 * std::sqrt(105 / (2 * pi));
    result = constant * std::exp(-2.0 * i * phi) * std::pow(sin(theta), 2.0) *
             cos(theta);
    return result;
  } else if (m == 2) {
    constant = 0.125 * std::sqrt(21 / pi);
    result = constant * std::exp(-1.0 * i * phi) * sin(theta) *
             (5.0 * std::pow(cos(theta), 2.0) - 1);
    return result;
  } else if (m == 3) {
    constant = 0.25 * std::sqrt(7 / pi);
    result = constant * (5.0 * std::pow(cos(theta), 3.0) - 3.0 * cos(theta));
    return result;
  } else if (m == 4) {
    constant = -0.125 * std::sqrt(21 / pi);
    result = constant * std::exp(i * phi) * sin(theta) *
             (5.0 * std::pow(cos(theta), 2.0) - 1);
    return result;
  } else if (m == 5) {
    constant = 0.25 * std::sqrt(105 / (2 * pi));
    result = constant * std::exp(2.0 * i * phi) * std::pow(sin(theta), 2.0) *
             cos(theta);
    return result;
  } else if (m == 6) {
    constant = -0.125 * std::sqrt(35 / pi);
    result = constant * std::exp(3.0 * i * phi) * std::pow(sin(theta), 3.0);
    return result;
  }
 
  return result;
}

Code

lookupTableQ3Vec()

std::vector< std::complex< double > > sph::lookupTableQ3Vec

(

std::array< double, 2 >

angles

)

Lookup table for Q3.

Calculates $Q_3$ using hard-coded look-up values.

Deprecated:

It is recommended to use the Boost version of this function, sph::spheriHarmo, instead.

Parameters

[in]

angles

The azimuth and polar angles of a particular point

Returns

a complex vector, of length $7$, calculated using spherical harmonics

Definition at line 107 of file bop.cpp.

                                                {
  // For keeping track of the index of the output vector
  std::vector<std::complex<double>> result;
  double theta = angles[1];
  double phi = angles[0];
 
  result.resize(7);
 
  for (int m = 0; m < 7; m++) {
    result[m] = sph::lookupTableQ3(m, angles);
  }
 
  return result;
}

Code

lookupTableQ6()

std::complex< double > sph::lookupTableQ6	(	int	m,
		std::array< double, 2 >	angles
	)

Lookup table for Q6 (m=0 to m=12)

Hard-coded calculations for determining $Q_6$.

It is recommended to use the general Boost version of this function, sph::spheriHarmo, instead.

Parameters

[in]	m	An int such that $-6<=m<=6$
[in]	angles	The azimuth and polar angles for a particular particle

Returns

a complex vector, of length $13$, calculated using hard-coded values

Definition at line 215 of file bop.cpp.

                                                                       {
  std::complex<double> result(0.0, 0.0);
  const double pi = std::acos(-1);
  const std::complex<double> i(0.0, 1.0);
  double constant;
  double theta = angles[1];
  double phi = angles[0];
 
  if (m == 0) {
    constant = 0.015625 * std::sqrt(3003 / pi);
    result = constant * std::exp(-6.0 * i * phi) * std::pow(sin(theta), 6.0);
    return result;
  } else if (m == 1) {
    constant = 0.09375 * std::sqrt(1001 / pi);
    result = constant * std::exp(-5.0 * i * phi) * std::pow(sin(theta), 5.0) *
             cos(theta);
    return result;
  } else if (m == 2) {
    constant = 0.09375 * std::sqrt(91 / (2 * pi));
    result = constant * std::exp(-4.0 * i * phi) * std::pow(sin(theta), 4.0) *
             (11.0 * std::pow(cos(theta), 2.0) - 1);
    return result;
  } else if (m == 3) {
    constant = 0.03125 * std::sqrt(1365 / pi);
    result = constant * std::exp(-3.0 * i * phi) * std::pow(sin(theta), 3.0) *
             (11.0 * std::pow(cos(theta), 3.0) - 3.0 * cos(theta));
    return result;
  } else if (m == 4) {
    constant = 0.015625 * std::sqrt(1365 / pi);
    result = constant * std::exp(-2.0 * i * phi) * std::pow(sin(theta), 2.0) *
             (33.0 * std::pow(cos(theta), 4.0) -
              18.0 * std::pow(cos(theta), 2.0) + 1.0);
    return result;
  } else if (m == 5) {
    constant = 0.0625 * std::sqrt(273 / (2 * pi));
    result = constant * std::exp(-1.0 * i * phi) * sin(theta) *
             (33.0 * std::pow(cos(theta), 5.0) -
              30.0 * std::pow(cos(theta), 3.0) + 5.0 * cos(theta));
    return result;
  } else if (m == 6) {
    constant = 0.03125 * std::sqrt(13 / pi);
    result = constant * (231.0 * std::pow(cos(theta), 6.0) -
                         315.0 * std::pow(cos(theta), 4.0) +
                         105.0 * std::pow(cos(theta), 2.0) - 5.0);
    return result;
  } else if (m == 7) {
    constant = -0.0625 * std::sqrt(273 / (2 * pi));
    result = constant * std::exp(i * phi) * sin(theta) *
             (33.0 * std::pow(cos(theta), 5.0) -
              30.0 * std::pow(cos(theta), 3.0) + 5.0 * cos(theta));
    return result;
  } else if (m == 8) {
    constant = 0.015625 * std::sqrt(1365 / pi);
    result = constant * std::exp(2.0 * i * phi) * std::pow(sin(theta), 2.0) *
             (33.0 * std::pow(cos(theta), 4.0) -
              18.0 * std::pow(cos(theta), 2.0) + 1.0);
    return result;
  } else if (m == 9) {
    constant = -0.03125 * std::sqrt(1365 / pi);
    result = constant * std::exp(3.0 * i * phi) * std::pow(sin(theta), 3.0) *
             (11.0 * std::pow(cos(theta), 3.0) - 3.0 * cos(theta));
    return result;
  } else if (m == 10) {
    constant = 0.09375 * std::sqrt(91 / (2 * pi));
    result = constant * std::exp(4.0 * i * phi) * std::pow(sin(theta), 4.0) *
             (11.0 * std::pow(cos(theta), 2.0) - 1);
    return result;
  } else if (m == 11) {
    constant = -0.09375 * std::sqrt(1001 / pi);
    result = constant * std::exp(5.0 * i * phi) * std::pow(sin(theta), 5.0) *
             cos(theta);
    return result;
  } else if (m == 12) {
    constant = 0.015625 * std::sqrt(3003 / pi);
    result = constant * std::exp(6.0 * i * phi) * std::pow(sin(theta), 6.0);
    return result;
  }
 
  return result;
}

Code

lookupTableQ6Vec()

std::vector< std::complex< double > > sph::lookupTableQ6Vec

(

std::array< double, 2 >

angles

)

Lookup table for Q6.

Calculates $Q_6$ using hard-coded values.

It is recommended to use the Boost version of this function, sph::spheriHarmo, instead.

Parameters

[in]

angles

The azimuth and polar angles for a particular particle

Returns

a complex vector, of length $13$, calculated using hard-coded values

Definition at line 189 of file bop.cpp.

                                                {
  // For keeping track of the index of the output vector
  std::vector<std::complex<double>> result;
  double theta = angles[1];
  double phi = angles[0];
 
  result.resize(13);
 
  for (int m = 0; m < 13; m++) {
    result[m] = sph::lookupTableQ6(m, angles);
  }
 
  return result;
}

Code

numStaggered()

int chill::numStaggered	(	molSys::PointCloud< molSys::Point< double >, double > *	yCloud,
		std::vector< std::vector< int >>	nList,
		int	jatom
	)

Finds the number of staggered bonds for a given atom of index jatom.

Calculates the number of staggered bonds of an atom with the given index.

Parameters

[in]	yCloud	The input molSys::PointCloud
[in]	nList	Row-ordered neighbour list by atom ID
[in]	jatom	The vector index of the current particle

Returns

an int value, holding the number of staggered bonds of the given particle

Definition at line 1201 of file bop.cpp.

                                                {
  int num_staggrd = 0;        // Number of staggered bonds
  molSys::bond_type bondType; // Bond type
  int num_bonds;              // No. of bonds of the jatom
  int nnumNeighbours =
      nList[jatom].size() - 1; // No. of nearest neighbours of index jatom
 
  // Loop over all bonds
  for (int i = 0; i < nnumNeighbours; i++) {
    bondType = yCloud->pts[jatom].c_ij[i].classifier;
    // If the bond is staggered increment the number of staggered bonds
    if (bondType == molSys::bond_type::staggered) {
      num_staggrd++;
    }
  } // end of loop over c_ij
 
  return num_staggrd;
}

Code

printIceType()

int chill::printIceType	(	molSys::PointCloud< molSys::Point< double >, double > *	yCloud,
		std::string	path,
		int	firstFrame,
		bool	isSlice = false,
		std::string	outputFileName = "superChill.txt"
	)

Prints out the iceType for a particular frame onto the terminal.

Prints out the molSys::atom_state_type per-particle ice type, for a particular frame, to a file.

Parameters

[in]	yCloud	The input molSys::PointCloud for the current frame
[in]	path	Path to the output directory to which ice types are written out to
[in]	firstFrame	First frame to be analyzed
[in]	isSlice	Determines whether there is a slice or not
[in]	outputFileName	File name of the output file, to which the per-particle ice types will be written out. The default file name is "superChill.txt"

Definition at line 1060 of file bop.cpp.

                                                            {
  int ih, ic, interIce, water, unknown, clath, interClath,
      total; // No. of particles of each type
  ih = ic = water = unknown = interIce = total = 0;
  clath = interClath = 0;
 
  for (int iatom = 0; iatom < yCloud->nop; iatom++) {
    // if(yCloud->pts[iatom].type != typeO){continue;}
    // check for slice
    if (isSlice) {
      if (yCloud->pts[iatom].inSlice == false) {
        continue;
      }
    }
    total++;
    if (yCloud->pts[iatom].iceType == molSys::atom_state_type::cubic) {
      ic++;
    } else if (yCloud->pts[iatom].iceType == molSys::atom_state_type::hexagonal) {
      ih++;
    } else if (yCloud->pts[iatom].iceType == molSys::atom_state_type::water) {
      water++;
    } else if (yCloud->pts[iatom].iceType == molSys::atom_state_type::interfacial) {
      interIce++;
    } else if (yCloud->pts[iatom].iceType == molSys::atom_state_type::clathrate) {
      clath++;
    } else if (yCloud->pts[iatom].iceType == molSys::atom_state_type::interClathrate) {
      interClath++;
    } else {
      unknown++;
    }
  }
 
  // --------------------
  // Create the directories if needed
  sout::makePath(path);
  std::string outputDirName = path + "bop";
  sout::makePath(outputDirName);
  // --------------------
  // Print to file
  std::ofstream outputFile;
  outputFile.open(path + "bop/" + outputFileName, std::ios_base::app);
  // --------------------
  // Write out the comment line
  if (yCloud->currentFrame == firstFrame) {
    outputFile << "Frame Ic Ih Interfacial Clath InterClath Water Total\n";
  }
  // --------------------
  outputFile << yCloud->currentFrame << " " << ic << " " << ih << " "
             << interIce << " " << clath << " " << interClath << " " << water
             << " " << total << "\n";
  outputFile.close();
 
  return 0;
}

Code

radialCoord()

std::array< double, 2 > sph::radialCoord

(

std::array< double, 3 >

cartCoord

)

Function for the azimuth and polar angles, given the Cartesian coordinates This function uses the Boost libraries.

Parameters

[in]

cartCoord

The Cartesian coordinates of a particular point

Returns

a double array, holding the azimuth and polar angles

Definition at line 81 of file bop.cpp.

                                                                  {
  // The output
  std::array<double, 2> result;
  // Point Definitions
  bg::model::point<long double, 3, bg::cs::cartesian> cartesianPoint;
  bg::model::point<long double, 3, bg::cs::spherical<bg::radian>> azuPoint;
  // Set Value (TODO: Recurse this)
  bg::set<0>(cartesianPoint, cartCoord[0]);
  bg::set<1>(cartesianPoint, cartCoord[1]);
  bg::set<2>(cartesianPoint, cartCoord[2]);
  // Transform
  bg::transform(cartesianPoint, azuPoint);
  result[0] = bg::get<0>(azuPoint);
  result[1] = bg::get<1>(azuPoint);
  return result;
}

Code

reclassifyWater()

molSys::PointCloud< molSys::Point< double >, double > chill::reclassifyWater	(	molSys::PointCloud< molSys::Point< double >, double > *	yCloud,
		std::vector< double > *	q6
	)

'Test' condition for classifying hexagonal ice using averaged q6 and q3 Checks water According to https://!pubs.rsc.org/en/content/articlehtml/2011/cp/c1cp22167a Gets c_ij and then classifies bond types according to the CHILL+ algorithm

Reclassifies atoms which may have been mis-classified as water using the averaged $q_6$ and $q_3$ parameters. This function can be called after both averaged $q_6$ and bond order correlation function $c_{ij}$ have been calculated as described here.

Parameters

[in,out]	yCloud	The output molSys::PointCloud
[in]	q6	Vector containing the previously calculated averaged $q_6$ values (using chill::getq6)

Definition at line 1011 of file bop.cpp.

                           {
  // If averaged q6 > 0.5, then consider it to be ice
  // If averaged q3 < -0.75 then it is ih or ic. If q3 < -0.85 then it is cubic,
  // otherwise it is hexagonal
  double avgQ3 = 0.0;
  int nnumNeighbours;
 
  for (int iatom = 0; iatom < yCloud->nop; iatom++) {
    // Check if it has been classified as water
    if (yCloud->pts[iatom].iceType == molSys::atom_state_type::water) {
      if ((*q6)[iatom] > 0.5) {
        avgQ3 = 0.0; // init to zero
        // Loop through all c_ij
        nnumNeighbours = yCloud->pts[iatom].c_ij.size();
        for (int j = 0; j < nnumNeighbours; j++) {
          avgQ3 += yCloud->pts[iatom].c_ij[j].c_value;
        }
        avgQ3 /= (double)nnumNeighbours;
 
        // If averaged q3 < -0.75, then reclassify
        if (avgQ3 <= -0.75) {
          if (avgQ3 < -0.85) {
            yCloud->pts[iatom].iceType = molSys::atom_state_type::reCubic;
          } // molSys::cubic
          else {
            yCloud->pts[iatom].iceType = molSys::atom_state_type::reHex;
          } // molSys::hexagonal
        }   // end of reclassification
      }     // check for solid atom!
    }       // end of check for water
  }         // End loop through every iatom
 
  return *yCloud;
}

Code

spheriHarmo()

std::vector< std::complex< double > > sph::spheriHarmo	(	int	orderL,
		std::array< double, 2 >	radialCoord
	)

Function for calculating spherical harmonics, that works for a general $l$.

This function uses the Boost libraries.

Parameters

[in]	orderL	The int value of $l$
[in]	radialCoord	Array containing the polar and azimuth angles

Returns

a complex vector, holding the complex spherical harmonics values, of length $2l+1$

Definition at line 57 of file bop.cpp.

                                                            {
  // For keeping track of the index of the output vector
  std::vector<std::complex<double>> result;
  std::complex<double> b; // Boost temp value
  int m;
 
  result.resize(2 * orderL + 1);
 
  for (int k = 0; k < 2 * orderL + 1; k++) {
    double theta = radialCoord[1];
    double phi = radialCoord[0];
    m = k - orderL;
    result[k] = boost::math::spherical_harmonic(orderL, m, theta, phi);
  }
 
  return result;
}

Code

Variable Documentation

ptq

std::vector<YlmAtom> chill::QlmAtom::ptq

Definition at line 170 of file bop.hpp.

ylm

std::vector<std::complex<double> > chill::YlmAtom::ylm

Definition at line 149 of file bop.hpp.