Ring

Collaboration diagram for Ring:

Modules
	Prism3

Namespaces
namespace	tum3
namespace	ring
	Topological network criteria functions.

Enumerations
enum class	ring::strucType {   ring::strucType::unclassified , ring::strucType::DDC , ring::strucType::HCbasal , ring::strucType::HCprismatic ,   ring::strucType::bothBasal , ring::strucType::bothPrismatic , ring::strucType::Prism , ring::strucType::deformedPrism ,   ring::strucType::mixedPrismRing }

Functions
int	tum3::topoUnitMatchingBulk (std::string path, std::vector< std::vector< int > > rings, std::vector< std::vector< int > > nList, molSys::PointCloud< molSys::Point< double >, double > *yCloud, int firstFrame, bool printClusters, bool onlyTetrahedral)
Eigen::MatrixXd	tum3::buildRefHC (std::string fileName)
	Build a reference Hexagonal cage, reading in from a template XYZ file.
Eigen::MatrixXd	tum3::buildRefDDC (std::string fileName)
	Build a reference Double-Diamond cage, reading in from a template XYZ file.
int	tum3::shapeMatchHC (molSys::PointCloud< molSys::Point< double >, double > *yCloud, const Eigen::MatrixXd &refPoints, cage::Cage cageUnit, std::vector< std::vector< int > > rings, std::vector< std::vector< int > > nList, std::vector< double > *quat, double *rmsd)
	Shape-matching for a target HC.
int	tum3::shapeMatchDDC (molSys::PointCloud< molSys::Point< double >, double > *yCloud, const Eigen::MatrixXd &refPoints, std::vector< cage::Cage > cageList, int cageIndex, std::vector< std::vector< int > > rings, std::vector< double > *quat, double *rmsd)
	Shape-matching for a target DDC.
int	tum3::updateRMSDatom (std::vector< std::vector< int > > rings, cage::Cage cageUnit, double rmsd, std::vector< double > *rmsdPerAtom, std::vector< int > *noOfCommonAtoms, std::vector< cage::iceType > atomTypes)
int	tum3::averageRMSDatom (std::vector< double > *rmsdPerAtom, std::vector< int > *noOfCommonAtoms)
	Average the RMSD per atom.
std::vector< cage::Cage >	tum3::topoBulkCriteria (std::string path, std::vector< std::vector< int > > rings, std::vector< std::vector< int > > nList, molSys::PointCloud< molSys::Point< double >, double > *yCloud, int firstFrame, int *numHC, int *numDDC, std::vector< ring::strucType > *ringType)
int	tum3::clusterCages (molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::string path, std::vector< std::vector< int > > rings, std::vector< cage::Cage > cageList, int numHC, int numDDC)
std::vector< int >	tum3::atomsFromCages (std::vector< std::vector< int > > rings, std::vector< cage::Cage > cageList, std::vector< int > clusterCages)
	Gets the atoms in the cages of a given cluster.
std::vector< std::vector< int > >	ring::getSingleRingSize (std::vector< std::vector< int > > rings, int ringSize)
	Returns a vector of vectors of rings of a single size.
bool	ring::hasCommonElements (std::vector< int > ring1, std::vector< int > ring2)
bool	ring::compareRings (std::vector< int > ring1, std::vector< int > ring2)
bool	ring::findTripletInRing (std::vector< int > ring, std::vector< int > triplet)
	Searches a particular ring for a triplet.
bool	ring::commonElementsInThreeRings (std::vector< int > ring1, std::vector< int > ring2, std::vector< int > ring3)
	Common elements in 3 rings.
std::vector< int >	ring::findsCommonElements (std::vector< int > ring1, std::vector< int > ring2)
	Returns the common elements of two rings.
int	ring::clearRingList (std::vector< std::vector< int > > &rings)
	Erases memory for a vector of vectors for a list of rings.
int	ring::assignPolygonType (std::vector< std::vector< int > > rings, std::vector< int > *atomTypes, std::vector< int > nRings)
molSys::PointCloud< molSys::Point< double >, double >	gen::getPointCloudOneAtomType (molSys::PointCloud< molSys::Point< double >, double > *yCloud, molSys::PointCloud< molSys::Point< double >, double > *outCloud, int atomTypeI, bool isSlice=false, std::array< double, 3 > coordLow=std::array< double, 3 >{0, 0, 0}, std::array< double, 3 > coordHigh=std::array< double, 3 >{0, 0, 0})
void	gen::moleculesInSingleSlice (molSys::PointCloud< molSys::Point< double >, double > *yCloud, bool clearPreviousSliceSelection=true, std::array< double, 3 > coordLow=std::array< double, 3 >{0, 0, 0}, std::array< double, 3 > coordHigh=std::array< double, 3 >{0, 0, 0})
void	gen::atomsInSingleSlice (molSys::PointCloud< molSys::Point< double >, double > *yCloud, bool clearPreviousSliceSelection=true, std::array< double, 3 > coordLow=std::array< double, 3 >{0, 0, 0}, std::array< double, 3 > coordHigh=std::array< double, 3 >{0, 0, 0})
void	gen::setAtomsWithSameMolID (molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::unordered_multimap< int, int > molIDAtomIDmap, int molID, bool inSliceValue=true)
void	ring::getEdgeMoleculesInRings (std::vector< std::vector< int > > rings, molSys::PointCloud< molSys::Point< double >, double > *oCloud, molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::array< double, 3 > coordLow, std::array< double, 3 > coordHigh, bool identicalCloud=false)
void	ring::printSliceGetEdgeMoleculesInRings (std::string path, std::vector< std::vector< int > > rings, molSys::PointCloud< molSys::Point< double >, double > *oCloud, molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::array< double, 3 > coordLow, std::array< double, 3 > coordHigh, bool identicalCloud=false)
std::vector< int >	ring::findPrisms (std::vector< std::vector< int > > rings, std::vector< strucType > *ringType, int *nPerfectPrisms, int *nImperfectPrisms, std::vector< std::vector< int > > nList, molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< double > *rmsdPerAtom, bool doShapeMatching=false)
bool	ring::basalPrismConditions (std::vector< std::vector< int > > nList, std::vector< int > *basal1, std::vector< int > *basal2)
bool	ring::relaxedPrismConditions (std::vector< std::vector< int > > nList, std::vector< int > *basal1, std::vector< int > *basal2)
bool	ring::discardExtraTetragonBlocks (std::vector< int > *basal1, std::vector< int > *basal2, molSys::PointCloud< molSys::Point< double >, double > *yCloud)
std::vector< std::vector< int > >	ring::keepAxialRingsOnly (std::vector< std::vector< int > > rings, molSys::PointCloud< molSys::Point< double >, double > *yCloud)
	Saves only axial rings out of all possible rings.
int	ring::prismAnalysis (std::string path, std::vector< std::vector< int > > rings, std::vector< std::vector< int > > nList, molSys::PointCloud< molSys::Point< double >, double > *yCloud, int maxDepth, int *atomID, int firstFrame, int currentFrame, bool doShapeMatching=false)
int	ring::assignPrismType (std::vector< std::vector< int > > rings, std::vector< int > listPrism, int ringSize, std::vector< ring::strucType > ringType, std::vector< int > *atomTypes, std::vector< ring::strucType > *atomState)
int	ring::deformedPrismTypes (std::vector< ring::strucType > atomState, std::vector< int > *atomTypes, int maxDepth)
	Get the atom type values for deformed prisms.
int	ring::rmAxialTranslations (molSys::PointCloud< molSys::Point< double >, double > *yCloud, int *atomID, int firstFrame, int currentFrame)
	Shift the entire ice nanotube and remove axial translations.
int	ring::polygonRingAnalysis (std::string path, std::vector< std::vector< int > > rings, std::vector< std::vector< int > > nList, molSys::PointCloud< molSys::Point< double >, double > *yCloud, int maxDepth, double sheetArea, int firstFrame)

Detailed Description

Enumeration Type Documentation

strucType

enum class ring::strucType

strong

Enumerator
unclassified	The ring is unclassified, which may be either water or a deformed type which cannot be classified by the criteria.
DDC	The ring belongs to a double-diamond cage (DDC).
HCbasal	The ring belongs only to a hexagonal cage (HC). Specifically, the ring is purely a basal ring of an HC.
HCprismatic	The ring belongs only to a hexagonal cage (HC); specifically the ring is purely a prismatic ring of an HC. It is not shared by a DDC.
bothBasal	The ring belongs to both a DDC and HC. It is a 'mixed' ring. The ring is also one of the basal rings of the HC of which it is part. A mixed ring must be a peripheral ring of the DDC of which it is part by definition.
bothPrismatic	The ring belongs to both a DDC and HC and is, thus, a 'mixed' ring. The ring is also one of the prismatic rings of the HC of which it is part. A mixed ring must be a peripheral ring of the DDC of which it is part by definition (can never be an equatorial ring of a DDC and also be part of an HC).
Prism	The ring belongs to a prism block, classified according to the prism identification scheme.
deformedPrism
mixedPrismRing

Definition at line 116 of file ring.hpp.

Code

                     {
  unclassified,
  DDC,
  HCbasal,
  HCprismatic,
  bothBasal,
  bothPrismatic,
  Prism,
  deformedPrism,
  mixedPrismRing
};

Function Documentation

assignPolygonType()

int ring::assignPolygonType	(	std::vector< std::vector< int > >	rings,
		std::vector< int > *	atomTypes,
		std::vector< int >	nRings
	)

Assign an atomType (equal to the number of nodes in the ring) given n-membered rings.

Assign an atomType (equal to the number of nodes in the ring) given the rings vector.

Parameters

[in]	rings	The vector of vectors containing the primitive rings, of a particular ring size.
[in,out]	atomTypes	A vector which contains a type for each atom, depending on it's type as classified by the prism identification scheme.
[in]	nRings	Number of rings.

Definition at line 40 of file ring.cpp.

Code

                                                   {
  // Every value in listPrism corresponds to an index in rings.
  // Every ring contains atom indices, corresponding to the indices (not atom
  // IDs) in rings
  int iring;        // Index of current ring
  int iatom;        // Index of current atom
  int ringSize;     // Ring size of the current ring
  int prevRingSize; // Ring size previously assigned to a point
 
  // Dummy value corresponds to a value of 1.
  // If an atom is shared by more than one ring type, it is assigned the
  // value 2.
 
  // Loop through every ring in rings
  for (int iring = 0; iring < rings.size(); iring++) {
    ringSize = rings[iring].size();
    // Loop through every element in iring
    for (int j = 0; j < ringSize; j++) {
      iatom = rings[iring][j]; // Atom index
      // Update the atom type
      if ((*atomTypes)[iatom] == 1) {
        (*atomTypes)[iatom] = ringSize;
      } // The atom is unclassified
      else {
        // Only update the ring type if the number is higher
        prevRingSize = (*atomTypes)[iatom]; // Previously assigned ring size
        if (ringSize > prevRingSize) {
          (*atomTypes)[iatom] = ringSize;
        } // end of assigning the new ring size
      }   // only update if the number is higher
    }     // end of loop through every atom in iring
  }       // end of loop through every ring
 
  return 0;
} // end of function

assignPrismType()

int ring::assignPrismType	(	std::vector< std::vector< int > >	rings,
		std::vector< int >	listPrism,
		int	ringSize,
		std::vector< ring::strucType >	ringType,
		std::vector< int > *	atomTypes,
		std::vector< ring::strucType > *	atomState
	)

Assign an atomType (equal to the number of nodes in the ring) given a vector with a list of indices of rings comprising the prisms

Assign an atomType (equal to the number of nodes in the ring) given a vector with a list of indices of rings comprising the prisms. Note that the ring indices in the listPrism vector correspond to the indices in the input rings vector of vectors.

Parameters

[in]	rings	The vector of vectors containing the primitive rings, of a particular ring size.
[in]	listPrism	The list of prism blocks found.
[in]	ringSize	The current ring size or number of nodes in each ring.
[in,out]	atomTypes	A vector which contains a type for each atom, depending on it's type as classified by the prism identification scheme.

Definition at line 696 of file topo_one_dim.cpp.

Code

                                                               {
  // Every value in listPrism corresponds to an index in rings.
  // Every ring contains atom indices, corresponding to the indices (not atom
  // IDs) in rings
  int iring;                    // Index of current ring
  int iatom;                    // Index of current atom
  ring::strucType currentState; // Current state of the ring being filled
 
  // Dummy value corresponds to a value of 1.
  // Each value is initialized to the value of 1.
 
  // Loop through every ring in rings
  for (int i = 0; i < listPrism.size(); i++) {
    iring = listPrism[i];
    // Get the state of all atoms in the ring
    currentState = ringType[iring];
    // Loop through every element in iring
    for (int j = 0; j < ringSize; j++) {
      iatom = rings[iring][j]; // Atom index
      // Update the atom type
      (*atomTypes)[iatom] = ringSize;
      //
      // Update the state of the atom
      if ((*atomState)[iatom] == ring::strucType::unclassified) {
        (*atomState)[iatom] = currentState;
      } // Update the unclassified atom
      else {
        if ((*atomState)[iatom] != currentState) {
          // For mixed, there is a preference
          if (currentState == ring::strucType::mixedPrismRing) {
            (*atomState)[iatom] = currentState;
          } // fill
          else if ((*atomState)[iatom] == ring::strucType::deformedPrism &&
                   currentState == ring::strucType::Prism) {
            (*atomState)[iatom] = ring::strucType::mixedPrismRing;
          } else if ((*atomState)[iatom] == ring::strucType::Prism &&
                     currentState == ring::strucType::deformedPrism) {
            (*atomState)[iatom] = ring::strucType::mixedPrismRing;
          }
        } //
      }   // already filled?
 
    } // end of loop through every atom in iring
  }   // end of loop through every ring
 
  return 0;
} // end of function

atomsFromCages()

std::vector< int > tum3::atomsFromCages	(	std::vector< std::vector< int > >	rings,
		std::vector< cage::Cage >	cageList,
		std::vector< int >	clusterCages
	)

Gets the atoms in the cages of a given cluster.

Returns a vector containing the atom indices of all the atoms in a cluster of cages, according to the input cageList vector of Cages

Definition at line 801 of file bulkTUM.cpp.

Code

                                                                   {
  //
  std::vector<int> atoms; // Contains the atom indices (not IDs) of atoms
  int ringSize = rings[0].size(); // Number of nodes in each ring
  int iring;                      // Index of ring in rings vector of vectors
  int icage;                      // Index of the current cage in cageList
  int iatom;                      // Atom index
  //
  for (int i = 0; i < clusterCages.size(); i++) {
    //
    icage = clusterCages[i];
    // Loop through every ring in the cage
    for (int j = 0; j < cageList[icage].rings.size(); j++) {
      iring = cageList[icage].rings[j]; // Current ring index
      // Loop through every atom in iring
      for (int k = 0; k < ringSize; k++) {
        iatom = rings[iring][k];
        atoms.push_back(iatom); // Add the atom
      }                         // Loop through every atom in iring
    }                           // loop through every ring in the current cage
  }                             // end of loop through all cages
 
  // Remove the duplicates in the atoms vector
  // Duplicate IDs must be removed
  //
  sort(atoms.begin(), atoms.end());
  auto ip = std::unique(atoms.begin(), atoms.end());
  // Resize peripheral rings to remove undefined terms
  atoms.resize(std::distance(atoms.begin(), ip));
  //
  // Return the vector of atoms
  return atoms;
}

atomsInSingleSlice()

void gen::atomsInSingleSlice	(	molSys::PointCloud< molSys::Point< double >, double > *	yCloud,
		bool	clearPreviousSliceSelection = true,
		std::array< double, 3 >	coordLow = std::array<double, 3>{0, 0, 0},
		std::array< double, 3 >	coordHigh = std::array<double, 3>{0, 0, 0}
	)

Given a pointCloud set the inSlice bool for every atom, if the atoms are inside the specified (single) region. Does not handle atoms in molecules straddling the boundary

Function that loops through a given input PointCloud and sets the inSlice bool for every Point according to whether the atom
is in the specified (single) slice or not. Not inclusive of atoms in molecules

Parameters

[in]	yCloud	The given input PointCloud
[in]	clearPreviousSliceSelection	sets all inSlice bool values to false before adding Points to the slice
[in]	coordLow	Contains the lower limits of the slice, if a slice is to be created
[in]	coordHigh	Contains the upper limits of the slice, if a slice is to be created

Definition at line 101 of file selection.cpp.

Code

                                   {
  //
  bool pointInSlice = true; // If the current point is in the slice, this is true (otherwise this is false)
 
  // --------------------
  // Set inSlice to false for every atom first
  if (clearPreviousSliceSelection)
  {
    for (int iatom = 0; iatom < yCloud->nop; iatom++) {
      yCloud->pts[iatom].inSlice = false;
    }
  } // end of init
  // --------------------
 
  // Loop through every iatom and check if the atom is in the slice or not 
  for (int iatom = 0; iatom < yCloud->nop; iatom++) {
    // Check if the current atom is in the slice 
    pointInSlice = sinp::atomInSlice(yCloud->pts[iatom].x, yCloud->pts[iatom].y, yCloud->pts[iatom].z,
                                               coordLow, coordHigh);
 
    yCloud->pts[iatom].inSlice = pointInSlice; // iatom is inside the slice
    //
  }
 
  return;
}

averageRMSDatom()

int tum3::averageRMSDatom	(	std::vector< double > *	rmsdPerAtom,
		std::vector< int > *	noOfCommonAtoms
	)

Average the RMSD per atom.

Average the RMSD per atom, by the number of common elements

Definition at line 529 of file bulkTUM.cpp.

Code

                                                           {
  //
  int nop = (*rmsdPerAtom).size(); // Number of particles
 
  for (int iatom = 0; iatom < nop; iatom++) {
    //
    if ((*noOfCommonAtoms)[iatom] == 0) {
      (*rmsdPerAtom)[iatom] = -1; // Dummy atom
    } else {
      (*rmsdPerAtom)[iatom] /= (*noOfCommonAtoms)[iatom];
    }
  } // end of averaging
 
  return 0;
} // end of function

basalPrismConditions()

bool ring::basalPrismConditions	(	std::vector< std::vector< int > >	nList,
		std::vector< int > *	basal1,
		std::vector< int > *	basal2
	)

Tests whether two rings are basal rings (true) or not (false) for a prism (strict criterion)

A function that checks to see if two basal rings are basal rings of a prism block or not, using the neighbour list information. The neighbour list nList is a row-ordered vector of vectors, containing atom indices (not atom IDs!). The first element of each subvector in nList is the atom index of the particle for which the other elements are the nearest neighbours.

Parameters

[in]	nList	Row-ordered neighbour list by atom index.
[in]	basal1	The vector for one of the basal rings.
[in]	basal2	The vector for the other basal ring.

Returns

A value that is true if the basal rings constitute a prism block, and false if they do not make up a prism block.

Definition at line 351 of file topo_one_dim.cpp.

Code

                                                        {
  int l1 = (*basal1)[0]; // first element of basal1 ring
  int ringSize =
      (*basal1).size(); // Size of the ring; each ring contains n elements
  int m_k;              // Atom ID of element in basal2
  bool l1_neighbour;    // m_k is a neighbour of l1(true) or not (false)
 
  // isNeighbour is initialized to false for all basal2 elements; indication if
  // basal2 elements are neighbours of basal1
  std::vector<bool> isNeighbour(ringSize, false);
  int kIndex;  // m_k index
  int lAtomID; // atomID of the current element of basal1
  int kAtomID; // atomID of the current element of basal2
 
  // ---------------------------------------------
  // COMPARISON OF basal2 ELEMENTS WITH l1
  for (int k = 0; k < ringSize; k++) {
    l1_neighbour = false;
    m_k = (*basal2)[k];
    // =================================
    // Checking to seee if m_k is be a neighbour of l1
    // Find m_k inside l1 neighbour list
    auto it = std::find(nList[l1].begin() + 1, nList[l1].end(), m_k);
 
    // If the element has been found, for l1
    if (it != nList[l1].end()) {
      l1_neighbour = true;
      kIndex = k;
      break;
    }
  } // l1 is a neighbour of m_k
 
  // If there is no nearest neighbour, then the two rings are not part of the
  // prism
  if (!l1_neighbour) {
    return false;
  }
 
  // ---------------------------------------------
  // NEIGHBOURS of basal1 in basal2
  isNeighbour[kIndex] = true;
 
  // All elements of basal1 must be neighbours of basal2
  for (int i = 1; i < ringSize; i++) {
    lAtomID = (*basal1)[i]; // element of basal1 ring
    for (int k = 0; k < ringSize; k++) {
      // Skip if already a neighbour
      if (isNeighbour[k]) {
        continue;
      }
      // Get the comparison basal2 element
      kAtomID = (*basal2)[k]; // element of basal2 ring;
 
      // Checking to see if kAtomID is a neighbour of lAtomID
      // Find kAtomID inside lAtomID neighbour list
      auto it1 =
          std::find(nList[lAtomID].begin() + 1, nList[lAtomID].end(), kAtomID);
 
      // If the element has been found, for l1
      if (it1 != nList[lAtomID].end()) {
        isNeighbour[k] = true;
      }
    } // Loop through basal2
  }   // end of check for neighbours of basal1
 
  // ---------------------------------------------
 
  // They should all be neighbours
  for (int k = 0; k < ringSize; k++) {
    // Check to see if any element is false
    if (!isNeighbour[k]) {
      return false;
    }
  }
 
  // Everything works out!
  return true;
}

buildRefDDC()

Eigen::MatrixXd tum3::buildRefDDC

(

std::string

fileName

)

Build a reference Double-Diamond cage, reading in from a template XYZ file.

Build a reference Double-Diamond cage, reading it in from a templates directory

Definition at line 435 of file bulkTUM.cpp.

Code

                                                  {
  //
  Eigen::MatrixXd refPnts(14, 3); // Reference point set (Eigen matrix)
  // Get the reference HC point set
  molSys::PointCloud<molSys::Point<double>, double>
      setCloud; // PointCloud for holding the reference point values
  // Variables for rings
  std::vector<std::vector<int>> nList; // Neighbour list
  std::vector<std::vector<int>> rings; // Rings
  std::vector<ring::strucType>
      ringType; // This vector will have a value for each ring inside
  std::vector<int> listDDC,
      listHC; // Contains atom indices of atoms making up DDCs and HCs
  std::vector<int> ddcOrder; // Atom indices of particles in the DDC
  // Make a list of all the DDCs and HCs
  std::vector<cage::Cage> cageList;
  int iring, jring;
  //
  // read in the XYZ file into the pointCloud setCloud
  //
  setCloud = sinp::readXYZ(fileName);
 
  // box lengths
  for (int i = 0; i < 3; i++) {
    setCloud.box[i] = 50;
  } // end of setting box lengths
  //
 
  nList = nneigh::neighListO(3.5, &setCloud, 1);
  // Neighbour list by index
  nList = nneigh::neighbourListByIndex(&setCloud, nList);
  // Find the vector of vector of rings
  rings = primitive::ringNetwork(nList, 6);
  // init the ringType vector
  ringType.resize(rings.size());
  // Find the DDCs
  listDDC = ring::findDDC(rings, &ringType, listHC, &cageList);
  // Save the order of the DDC in a vector
  ddcOrder = pntToPnt::relOrderDDC(0, rings, cageList);
  // Get the reference point set
  refPnts = pntToPnt::changeDiaCageOrder(&setCloud, ddcOrder, 0);
  //
  return refPnts;
}

buildRefHC()

Eigen::MatrixXd tum3::buildRefHC

(

std::string

fileName

)

Build a reference Hexagonal cage, reading in from a template XYZ file.

Build a reference Hexagonal cage, reading it in from a templates directory

Definition at line 377 of file bulkTUM.cpp.

Code

                                                 {
  //
  Eigen::MatrixXd refPnts(12, 3); // Reference point set (Eigen matrix)
  // Get the reference HC point set
 // molSys::PointCloud<molSys::Point<double>, double>
  //    setCloud; // PointCloud for holding the reference point values
  // Variables for rings
  std::vector<std::vector<int>> nList; // Neighbour list
  std::vector<std::vector<int>> rings; // Rings
  std::vector<ring::strucType>
      ringType;            // This vector will have a value for each ring inside
  std::vector<int> listHC; // Contains atom indices of atoms making up HCs
  // Make a list of all the DDCs and HCs
  std::vector<cage::Cage> cageList;
  int iring, jring;
  //
  // read in the XYZ file into the pointCloud setCloud
  //
  auto setCloud = sinp::readXYZ(fileName);
  // box lengths
  for (int i = 0; i < 3; i++) {
    setCloud.box[i] = 50;
  } // end of setting box lengths
  //
 
  nList = nneigh::neighListO(3.5, &setCloud, 1);
  // Neighbour list by index
  nList = nneigh::neighbourListByIndex(&setCloud, nList);
  // Find the vector of vector of rings
  rings = primitive::ringNetwork(nList, 6);
  // init the ringType vector
  ringType.resize(rings.size());
  // Find the HCs
  listHC = ring::findHC(rings, &ringType, nList, &cageList);
  // Get the basal rings from cageList
  iring = cageList[0].rings[0];
  jring = cageList[0].rings[1];
  //
  std::vector<int> matchedBasal1,
      matchedBasal2; // Re-ordered basal rings 1 and 2
  // Reordered basal rings
  // Getting the target Eigen vectors
  // Get the re-ordered matched basal rings, ordered with respect to each
  // other
 
  pntToPnt::relOrderHC(&setCloud, rings[iring], rings[jring], nList,
                       &matchedBasal1, &matchedBasal2);
  // Get the reference point set
  refPnts =
      pntToPnt::changeHexCageOrder(&setCloud, matchedBasal1, matchedBasal2, 0);
  //
  return refPnts;
}

clearRingList()

int ring::clearRingList

(

std::vector< std::vector< int > > &

rings

)

Erases memory for a vector of vectors for a list of rings.

Deletes the memory of a vector of vectors (specifically for rings).

Parameters

[in,out]

rings

The vector of vectors to be cleared.

Definition at line 22 of file ring.cpp.

Code

                                                        {
  //
  std::vector<std::vector<int>> tempEmpty;
 
  rings.swap(tempEmpty);
 
  return 0;
}

clusterCages()

int tum3::clusterCages	(	molSys::PointCloud< molSys::Point< double >, double > *	yCloud,
		std::string	path,
		std::vector< std::vector< int > >	rings,
		std::vector< cage::Cage >	cageList,
		int	numHC,
		int	numDDC
	)

Clustering Clusters cages using the Stillinger algorithm and prints out individual XYZ files of clusters.

Groups cages into clusters and prints out each cluster into an XYZ file

Definition at line 644 of file bulkTUM.cpp.

Code

                           {
  //
  std::vector<int> linkedList; // Linked list for the clustered cages
  int nCages = numHC + numDDC; // Number of cages in total
  int j, temp;                 // Variables used inside the loops
  bool isNeighbour; // true if cages are adjacent to each other, false otherwise
  std::vector<bool>
      visited; // To make sure you don't go through the same cages again.
  int singleDDCs, singleHCs; // Number of single DDCs and HCs
  // Units are in Angstrom^3
  double volDDC = 72.7;          // The calculated alpha volume of a single DDC
  double volHC = 40.63;          // The calculated alpha volume of a single HC
  int nextElement;               // value in list
  int index;                     // starting index value
  int currentIndex;              // Current cage index
  int iClusterNumber;            // Number in the current cluster
  std::vector<int> nClusters;    // Number of cages in each cluster
  std::vector<int> clusterCages; // Indices of each cage in the cluster
  std::vector<int> atoms; // Vector containing the atom indices of all the
                          // atoms in the cages
  cage::cageType type;    // Type of the cages in the cluster
  // -----------------------------------------------------------
  // INITIALIZATION
  linkedList.resize(nCages); // init to dummy value
  // Initial values of the list.
  for (int icage = 0; icage < nCages; icage++) {
    // Assign the index as the ID
    linkedList[icage] = icage; // Index
  }                            // init of cluster IDs
  // -----------------------------------------------------------
  // GETTING THE LINKED LIST
  //
  // HCs first
  for (int i = 0; i < nCages - 1; i++) {
    // If iatom is already in a cluster, skip it
    if (linkedList[i] != i) {
      continue;
    } // Already part of a cluster
    //
    j = i; // Init of j
    // Execute the next part of the loop while j is not equal to i
    do {
      //
      // Go through the rest of the atoms (KLOOP)
      for (int k = i + 1; k < nCages; k++) {
        // Skip if already part of a cluster
        if (linkedList[k] != k) {
          continue;
        } // Already part of a cluster
        //
        // k and j must be of the same type
        if (cageList[k].type != cageList[j].type) {
          continue;
        } // skip if k and j are not of the same type
        // Check to see if k is a nearest neighbour of j
        isNeighbour =
            ring::hasCommonElements(cageList[k].rings, cageList[j].rings);
        if (isNeighbour) {
          // Swap!
          temp = linkedList[j];
          linkedList[j] = linkedList[k];
          linkedList[k] = temp;
        } // j and k are neighbouring cages
      }   // end of loop through k (KLOOP)
      //
      j = linkedList[j];
    } while (j != i); // end of control for j!=i
  } // end of getting the linked list (looping through every i)
  // -----------------------------------------------------------
  // WRITE-OUTS
  //
  // init
  visited.resize(nCages);
  singleDDCs = 0;
  singleHCs = 0;
 
  for (int i = 0; i < nCages; i++) {
    //
    if (visited[i]) {
      continue;
    } // Already counted
    // Now that the cage has been visited, set it to true
    visited[i] = true; // Visited
    // If only one cage is in the cluster
    if (linkedList[i] == i) {
      // Add to the number of single DDCs
      if (cageList[i].type == cage::cageType::DoubleDiaC) {
        singleDDCs++;
      } // add to the single DDCs
      else {
        singleHCs++;
      } // single HCs
      continue;
    } // cluster has only one cage
    //
    // init
    clusterCages.resize(0);
    //
    type = cageList[i].type;
    currentIndex = i;
    nextElement = linkedList[currentIndex];
    index = i;
    iClusterNumber = 1;        // at least one value
    clusterCages.push_back(i); // Add the first cage
    // Get the other cages in the cluster
    while (nextElement != index) {
      iClusterNumber++;
      currentIndex = nextElement;
      visited[currentIndex] = true;
      clusterCages.push_back(currentIndex); // Add the first cage
      nextElement = linkedList[currentIndex];
    } // get number
    // Update the number of molecules in the cluster
    nClusters.push_back(iClusterNumber);
    // --------------
    // PRINT OUT
    atoms = tum3::atomsFromCages(rings, cageList, clusterCages);
    int clusterID = nClusters.size();
    // Write out the cluster
    sout::writeXYZcluster(path, yCloud, atoms, clusterID, type);
    // Print out the cages into an XYZ file
    // --------------
  } // end of loop through
  // -----------------------------------------------------------
  // Write out the stuff for single cages
  // ----------------
  std::ofstream outputFile;
  // Make the output directory if it doesn't exist
  sout::makePath(path);
  std::string outputDirName = path + "bulkTopo";
  sout::makePath(outputDirName);
  outputDirName = path + "bulkTopo/clusterXYZ/";
  sout::makePath(outputDirName);
  outputDirName = path + "bulkTopo/clusterXYZ/frame-" +
                  std::to_string(yCloud->currentFrame);
  sout::makePath(outputDirName);
  // ----------------
  // Write output to file inside the output directory
  outputFile.open(path + "bulkTopo/clusterXYZ/frame-" +
                  std::to_string(yCloud->currentFrame) + "/info.dat");
  outputFile << "# volDDC volHC in Angstrom^3\n";
  outputFile << singleDDCs * volDDC << " " << singleHCs * volHC << "\n";
  outputFile << "There are " << singleDDCs << " single DDCs\n";
  outputFile << "There are " << singleHCs << " single HCs\n";
  outputFile.close(); // Close the file
  // -----------------------------------------------------------
  return 0;
} // end of the function

commonElementsInThreeRings()

bool ring::commonElementsInThreeRings	(	std::vector< int >	ring1,
		std::vector< int >	ring2,
		std::vector< int >	ring3
	)

Common elements in 3 rings.

Function that finds the common elements in three input rings

Parameters

[in]	ring1	The first ring.
[in]	ring2	The second ring.
[in]	ring3	The third ring.

Returns

A value which is true if the three rings have at least one common element, and false if the three rings have no elements in common.

Definition at line 114 of file ring.cpp.

Code

                                                            {
  std::vector<int>
      common1; // Vector containing the common elements of the first two rings
  std::vector<int>
      common2; // Vector containing the common elements of the three rings
 
  // Common elements among the first two rings
  common1 = ring::findsCommonElements(ring1, ring2);
  if (common1.size() == 0) {
    return false;
  } // no common elements in ring1 and ring2
 
  // Common elements among all three rings
  common2 = ring::findsCommonElements(common1, ring3);
 
  // If no common elements were found:
  if (common2.size() == 0) {
    return false;
  } // no common elements between ring1, ring2, and ring3
 
  return true; // Common elements found!
}

compareRings()

bool ring::compareRings	(	std::vector< int >	ring1,
		std::vector< int >	ring2
	)

Compares two disordered vectors and checks to see if they contain the same elements

Checks to see if two vectors (ring1, ring2) have the same elements (which may or may not be disordered). The order or sequence of elements is not important, so the rings are sorted. Returns true if the rings have the same elements

Parameters

[in]	ring1	The first ring.
[in]	ring2	The second ring.

Returns

A bool; true if the rings contain the same elements (not necessarily in the same sequence) and false if they do not have the same elements.

Definition at line 258 of file ring.cpp.

Code

                                                                  {
  // Sort the rings first
  sort(ring1.begin(), ring1.end()); // Sort ring1 by ID
  sort(ring2.begin(), ring2.end()); // Sort ring2 by ID
  bool result;
 
  (ring1 == ring2) ? result = true : result = false;
 
  return result;
}

deformedPrismTypes()

int ring::deformedPrismTypes	(	std::vector< ring::strucType >	atomState,
		std::vector< int > *	atomTypes,
		int	maxDepth
	)

Get the atom type values for deformed prisms.

Assign an atomType value for atoms belonging to deformed prisms.

Parameters

[in]	rings	The vector of vectors containing the primitive rings, of a particular ring size.
[in]	listPrism	The list of prism blocks found.
[in]	ringSize	The current ring size or number of nodes in each ring.
[in,out]	atomTypes	A vector which contains a type for each atom, depending on it's type as classified by the prism identification scheme.

Definition at line 759 of file topo_one_dim.cpp.

Code

                                                                      {
  //
  int nop = atomState.size(); // Number of particles
 
  // Loop through all particles
  for (int iatom = 0; iatom < nop; iatom++) {
    // Check the atom state
    // Deformed
    if (atomState[iatom] == ring::strucType::deformedPrism) {
      (*atomTypes)[iatom] += maxDepth - 2;
    } // type for a deformed prism atom
    else if (atomState[iatom] == ring::strucType::mixedPrismRing) {
      (*atomTypes)[iatom] = 2;
    } // type for a mixed prism ring
  }   // end of reassignation
 
  //
  return 0;
}

discardExtraTetragonBlocks()

bool ring::discardExtraTetragonBlocks	(	std::vector< int > *	basal1,
		std::vector< int > *	basal2,
		molSys::PointCloud< molSys::Point< double >, double > *	yCloud
	)

Checks whether two 4-membered rings are parallel in one dimension or not to prevent overcounting

A function that discards basal tetragonal rings which are not oriented perpendicular to the axial direction. This is will only work for the XY, YZ or XZ planes. If true is returned, then the rings are axial. This function is only needed for tetragonal prism blocks because in the other cases, the basal rings and lateral planes have a different number of nodes.

Parameters

[in]	basal1	The first basal ring.
[in]	basal2	The other candidate basal ring.
[in]	yCloud	The input PointCloud.

Returns

A bool value which is true if the candidate basal rings are axial, and is false if the planes are lateral planes.

Definition at line 487 of file topo_one_dim.cpp.

Code

                                                           {
  int ringSize =
      (*basal1).size(); // Size of the ring; each ring contains n elements
  int iatomIndex,
      jatomIndex;  // Indices of the elements in basal1 and basal2 respectively
  double r_i, r_j; // Coordinates in the axial dimension of iatom and jatom of
                   // basal1 and basal2 respectively
  int axialDim;    // 0 for x, 1 for y and 2 for z dimensions respectively
  // Variables for getting the projected area
  bool axialBasal1, axialBasal2; // bools for checking if basal1 and basal2 are
                                 // axial (true) respectively
  double areaXY, areaXZ,
      areaYZ; // Projected area on the XY, XZ and YZ planes respectively
  // ----------------------------------------
  // Find the axial dimension for a quasi-one-dimensional ice nanotube
  // The axial dimension will have the largest box length
  // Index -> axial dimension
  // 0 -> x dim
  // 1 -> y dim
  // 2 -> z dim
  axialDim = std::max_element(yCloud->box.begin(), yCloud->box.end()) -
             yCloud->box.begin();
  // ----------------------------------------
  // Calculate projected area onto the XY, YZ and XZ planes for basal1
  axialBasal1 = false; // Init to false
  axialBasal2 = false; // Init
 
  // Init the projected area
  areaXY = 0.0;
  areaXZ = 0.0;
  areaYZ = 0.0;
 
  jatomIndex = (*basal1)[0];
 
  // All points except the first pair
  for (int k = 1; k < ringSize; k++) {
    iatomIndex = (*basal1)[k]; // Current vertex
 
    // Add to the polygon area
    // ------
    // XY plane
    areaXY += (yCloud->pts[jatomIndex].x + yCloud->pts[iatomIndex].x) *
              (yCloud->pts[jatomIndex].y - yCloud->pts[iatomIndex].y);
    // ------
    // XZ plane
    areaXZ += (yCloud->pts[jatomIndex].x + yCloud->pts[iatomIndex].x) *
              (yCloud->pts[jatomIndex].z - yCloud->pts[iatomIndex].z);
    // ------
    // YZ plane
    areaYZ += (yCloud->pts[jatomIndex].y + yCloud->pts[iatomIndex].y) *
              (yCloud->pts[jatomIndex].z - yCloud->pts[iatomIndex].z);
    // ------
    jatomIndex = iatomIndex;
  }
 
  // Closure point
  iatomIndex = (*basal1)[0];
  // ------
  // XY plane
  areaXY += (yCloud->pts[jatomIndex].x + yCloud->pts[iatomIndex].x) *
            (yCloud->pts[jatomIndex].y - yCloud->pts[iatomIndex].y);
  // ------
  // XZ plane
  areaXZ += (yCloud->pts[jatomIndex].x + yCloud->pts[iatomIndex].x) *
            (yCloud->pts[jatomIndex].z - yCloud->pts[iatomIndex].z);
  // ------
  // YZ plane
  areaYZ += (yCloud->pts[jatomIndex].y + yCloud->pts[iatomIndex].y) *
            (yCloud->pts[jatomIndex].z - yCloud->pts[iatomIndex].z);
  // ------
  // The actual projected area is half of this
  areaXY *= 0.5;
  areaXZ *= 0.5;
  areaYZ *= 0.5;
  // Get the absolute value
  areaXY = fabs(areaXY);
  areaXZ = fabs(areaXZ);
  areaYZ = fabs(areaYZ);
 
  // If the axial dimension is x, y, or z:
  // then the maximum basal area should be in the YZ, XZ and XY dimensions
  // respectively
  // x dim
  if (axialDim == 0) {
    if (areaYZ > areaXY && areaYZ > areaXZ) {
      axialBasal1 = true;
    } // end of check for axial ring for basal1
  }   // x dim
  // y dim
  else if (axialDim == 1) {
    if (areaXZ > areaXY && areaXZ > areaYZ) {
      axialBasal1 = true;
    } // end of check for axial ring for basal1
  }   // x dim
  // z dim
  else if (axialDim == 2) {
    if (areaXY > areaXZ && areaXY > areaYZ) {
      axialBasal1 = true;
    } // end of check for axial ring for basal1
  }   // x dim
  else {
    std::cerr << "Could not find the axial dimension.\n";
    return false;
  }
  // ----------------------------------------
  // Calculate projected area onto the XY, YZ and XZ planes for basal2
 
  // Init the projected area
  areaXY = 0.0;
  areaXZ = 0.0;
  areaYZ = 0.0;
 
  jatomIndex = (*basal2)[0];
 
  // All points except the first pair
  for (int k = 1; k < ringSize; k++) {
    iatomIndex = (*basal2)[k]; // Current vertex
 
    // Add to the polygon area
    // ------
    // XY plane
    areaXY += (yCloud->pts[jatomIndex].x + yCloud->pts[iatomIndex].x) *
              (yCloud->pts[jatomIndex].y - yCloud->pts[iatomIndex].y);
    // ------
    // XZ plane
    areaXZ += (yCloud->pts[jatomIndex].x + yCloud->pts[iatomIndex].x) *
              (yCloud->pts[jatomIndex].z - yCloud->pts[iatomIndex].z);
    // ------
    // YZ plane
    areaYZ += (yCloud->pts[jatomIndex].y + yCloud->pts[iatomIndex].y) *
              (yCloud->pts[jatomIndex].z - yCloud->pts[iatomIndex].z);
    // ------
    jatomIndex = iatomIndex;
  }
 
  // Closure point
  iatomIndex = (*basal2)[0];
  // ------
  // XY plane
  areaXY += (yCloud->pts[jatomIndex].x + yCloud->pts[iatomIndex].x) *
            (yCloud->pts[jatomIndex].y - yCloud->pts[iatomIndex].y);
  // ------
  // XZ plane
  areaXZ += (yCloud->pts[jatomIndex].x + yCloud->pts[iatomIndex].x) *
            (yCloud->pts[jatomIndex].z - yCloud->pts[iatomIndex].z);
  // ------
  // YZ plane
  areaYZ += (yCloud->pts[jatomIndex].y + yCloud->pts[iatomIndex].y) *
            (yCloud->pts[jatomIndex].z - yCloud->pts[iatomIndex].z);
  // ------
  // The actual projected area is half of this
  areaXY *= 0.5;
  areaXZ *= 0.5;
  areaYZ *= 0.5;
  // Get the absolute value
  areaXY = fabs(areaXY);
  areaXZ = fabs(areaXZ);
  areaYZ = fabs(areaYZ);
 
  // Check if xy projected area is the greatest
  // If the axial dimension is x, y, or z:
  // then the maximum basal area should be in the YZ, XZ and XY dimensions
  // respectively
  // x dim
  if (axialDim == 0) {
    if (areaYZ > areaXY && areaYZ > areaXZ) {
      axialBasal2 = true;
    } // end of check for axial ring for basal1
  }   // x dim
  // y dim
  else if (axialDim == 1) {
    if (areaXZ > areaXY && areaXZ > areaYZ) {
      axialBasal2 = true;
    } // end of check for axial ring for basal1
  }   // x dim
  // z dim
  else if (axialDim == 2) {
    if (areaXY > areaXZ && areaXY > areaYZ) {
      axialBasal2 = true;
    } // end of check for axial ring for basal1
  }   // x dim
  else {
    std::cerr << "Could not find the axial dimension.\n";
    return false;
  }
  // ----------------------------------------
 
  // Now check if basal1 and basal2 are axial or not
  if (axialBasal1 == true && axialBasal2 == true) {
    return true;
  } else {
    return false;
  } // Check for basal1 and basal2
}

findPrisms()

std::vector< int > ring::findPrisms	(	std::vector< std::vector< int > >	rings,
		std::vector< strucType > *	ringType,
		int *	nPerfectPrisms,
		int *	nImperfectPrisms,
		std::vector< std::vector< int > >	nList,
		molSys::PointCloud< molSys::Point< double >, double > *	yCloud,
		std::vector< double > *	rmsdPerAtom,
		bool	doShapeMatching = false
	)

Find out which rings are prisms. Returns a vector containing all the ring IDs which are prisms

Determines which rings are n-sided prisms. This function returns a vector which contains the ring indices of all the rings which are prisms. The ring indices correspond to the index of the rings inside the vector of vector rings, starting from 0. Prism rings can be found using a three-step procedure, in which first two basal rings are found. Prismatic rings are simply rings which share every face made by upper and lower triplets of the basal rings The neighbour list is also required as an input, which is a vector of vectors, containing atom IDs. The first element of the neighbour list is the atom index of the atom for which the other elements are nearest neighbours.\

Parameters

[in]	rings	The input vector of vectors containing the primitive rings of a single ring size (number of nodes).
[in]	ringType	A vector containing a ring::strucType value (a classification type) for each ring.
[in]	nPrisms	The number of prism blocks identified for the number of nodes.
[in]	nList	The row-ordered neighbour list (by atom index).
[in]	yCloud	The input PointCloud.

Returns

A vector containing the ring indices of all the rings which have been classified as prisms. The indices are with respect to the input rings vector of vectors.

Definition at line 193 of file topo_one_dim.cpp.

Code

                                                                       {
  std::vector<int> listPrism;
  int totalRingNum = rings.size(); // Total number of rings
  std::vector<int> basal1;         // First basal ring
  std::vector<int> basal2;         // Second basal ring
  bool cond1, cond2; // Conditions for rings to be basal (true) or not (false)
  bool relaxedCond;  // Condition so that at least one bond exists between the
                     // two basal rings
  bool isAxialPair;  // Basal rings should be parallel in one dimension to
                     // prevent overcounting
  int ringSize = rings[0].size(); // Number of nodes in each ring
  *nImperfectPrisms = 0;          // Number of undeformed prisms
  *nPerfectPrisms = 0;            // Number of undeformed prisms
  // Matrix for the reference ring for a given ringSize.
  Eigen::MatrixXd refPointSet(ringSize, 3);
 
  // Get the reference ring point set for a given ring size.
  // Get the axial dimension
  int axialDim = std::max_element(yCloud->box.begin(), yCloud->box.end()) -
                 yCloud->box.begin();
  refPointSet = pntToPnt::getPointSetRefRing(ringSize, axialDim);
  //
 
  // Two loops through all the rings are required to find pairs of basal rings
  for (int iring = 0; iring < totalRingNum - 1; iring++) {
    cond1 = false;
    cond2 = false;
    basal1 = rings[iring]; // Assign iring to basal1
    // Loop through the other rings to get a pair
    for (int jring = iring + 1; jring < totalRingNum; jring++) {
      basal2 = rings[jring]; // Assign jring to basal2
      // ------------
      // Put extra check for axial basal rings if shapeMatching is being done
      if (doShapeMatching == true || ringSize == 4) {
        isAxialPair = false; // init
        isAxialPair =
            ring::discardExtraTetragonBlocks(&basal1, &basal2, yCloud);
        if (isAxialPair == false) {
          continue;
        }
      } // end of check for tetragonal prism blocks
      // ------------
      // Step one: Check to see if basal1 and basal2 have common
      // elements or not. If they don't, then they cannot be basal rings
      cond1 = ring::hasCommonElements(basal1, basal2);
      if (cond1 == true) {
        continue;
      }
      // -----------
      // Step two and three: One of the elements of basal2 must be the nearest
      // neighbour of the first (index0; l1) If m_k is the nearest neighbour of
      // l1, m_{k+1} ... m_{k+(n-1)} must be neighbours of l_i+1 etc or l_i-1
      cond2 = ring::basalPrismConditions(nList, &basal1, &basal2);
      // If cond2 is false, the strict criteria for prisms has not been met
      if (cond2 == false) {
        // Skip if shape-matching is not desired
        if (!doShapeMatching) {
          continue;
        } // shape-matching not desired
        //
        // If shape-matching is to be done:
        // Check for the reduced criteria fulfilment
        relaxedCond = ring::relaxedPrismConditions(nList, &basal1, &basal2);
        // Skip if relaxed criteria are not met
        if (relaxedCond == false) {
          continue;
        } // end of skipping if the prisms do not fulfil relaxed criteria
 
        // Do shape matching here
        bool isDeformedPrism = match::matchPrism(
            yCloud, nList, refPointSet, &basal1, &basal2, rmsdPerAtom, false);
 
        // Success! The rings are basal rings of a deformed prism!
        if (isDeformedPrism) {
          //
          // // Update the number of prism blocks
          *nImperfectPrisms += 1;
          // Update iring
          if ((*ringType)[iring] == ring::strucType::unclassified) {
            (*ringType)[iring] = ring::strucType::deformedPrism;
            listPrism.push_back(iring);
          } else if ((*ringType)[iring] == ring::strucType::Prism) {
            (*ringType)[iring] = ring::strucType::mixedPrismRing;
          } // if it is deformed
          // Update jring
          if ((*ringType)[jring] == ring::strucType::unclassified) {
            (*ringType)[jring] = ring::strucType::deformedPrism;
            listPrism.push_back(jring);
          } else if ((*ringType)[jring] == ring::strucType::Prism) {
            (*ringType)[jring] = ring::strucType::mixedPrismRing;
          } // if it is deformed
        }   // end of update of ring types
 
        // // Write outs
        // // Now write out axial basal rings
        // sout::writeBasalRingsPrism(&basal1, &basal2, nDeformedPrisms, nList,
        //                            yCloud, true);
      } // end of reduced criteria
      // Strict criteria
      else {
        // Update the number of prism blocks
        *nPerfectPrisms += 1;
        // Update iring
        if ((*ringType)[iring] == ring::strucType::unclassified) {
          (*ringType)[iring] = ring::strucType::Prism;
          listPrism.push_back(iring);
        } else if ((*ringType)[iring] == ring::strucType::deformedPrism) {
          (*ringType)[iring] = ring::strucType::mixedPrismRing;
        } // if it is deformed
        // Update jring
        if ((*ringType)[jring] == ring::strucType::unclassified) {
          (*ringType)[jring] = ring::strucType::Prism;
          listPrism.push_back(jring);
        } else if ((*ringType)[jring] == ring::strucType::deformedPrism) {
          (*ringType)[jring] = ring::strucType::mixedPrismRing;
        } // if it is deformed
        //
        // Shape-matching to get the RMSD (if shape-matching is desired)
        if (doShapeMatching) {
          bool isKnownPrism = match::matchPrism(
              yCloud, nList, refPointSet, &basal1, &basal2, rmsdPerAtom, true);
        } // end of shape-matching to get rmsd
        //
        // // Now write out axial basal rings for convex hull calculations
        // sout::writePrisms(&basal1, &basal2, *nPrisms, yCloud);
        // // Write out prisms for shape-matching
        // sout::writeBasalRingsPrism(&basal1, &basal2, *nPrisms, nList, yCloud,
        //                            false);
        // -----------
      } // end of strict criteria
 
    } // end of loop through rest of the rings to get the second basal ring
  }   // end of loop through all rings for first basal ring
 
  sort(listPrism.begin(), listPrism.end());
  auto ip = std::unique(listPrism.begin(), listPrism.end());
  // Resize peripheral rings to remove undefined terms
  listPrism.resize(std::distance(listPrism.begin(), ip));
 
  return listPrism;
}

findsCommonElements()

std::vector< int > ring::findsCommonElements	(	std::vector< int >	ring1,
		std::vector< int >	ring2
	)

Returns the common elements of two rings.

Function that finds and returns a vector containing the elements shared by two input rings (each ring is a vector).

Parameters

[in]	ring1	The first ring.
[in]	ring2	The second ring.

Returns

A vector containing the common elements between the input rings.

Definition at line 85 of file ring.cpp.

Code

                                                                 {
  //
  std::vector<int> common;
  int iatom; // Index to search for
 
  for (int i = 0; i < ring1.size(); i++) {
    iatom = ring1[i];
    // Search for iatom in ring2
 
    auto it = std::find(ring2.begin(), ring2.end(), iatom);
 
    if (it != ring2.end()) {
      common.push_back(iatom);
    } // iatom was found!
  }   // end of loop through every element of ring1
 
  return common;
}

findTripletInRing()

bool ring::findTripletInRing	(	std::vector< int >	ring,
		std::vector< int >	triplet
	)

Searches a particular ring for a triplet.

Function that finds out if a given triplet is is present (in the same order or in reversed order) in a given ring.

Parameters

[in]	ring	The input ring containing the indices of atoms.
[in]	triplet	Vector containing the triplet, for whose presence the input ring vector will be checked.

Returns

A bool value which is true if the triplet is present in the ring, and is false if the triplet is not in the ring.

Definition at line 148 of file ring.cpp.

Code

                                                                        {
  //
  int ringSize = ring.size();   // should be 6
  std::vector<int> ringTriplet; // triplet from the ring to be searched
  int kIndex;                   // Used for making the triplet
 
  // Loop through every possible triplet in the ring to be searched
  for (int i = 0; i < ringSize; i++) {
    ringTriplet.clear(); // init
    // Get the first element of the ring triplet
    ringTriplet.push_back(ring[i]);
    //
    // Get the next two elements
    for (int k = 1; k < 3; k++) {
      kIndex = i + k;
      // Wrap-around
      if (kIndex >= ringSize) {
        kIndex -= ringSize;
      } // end of wrap-around
      ringTriplet.push_back(ring[kIndex]);
    } // next two elements of ringTriplet
    //
    // Obtained ringTriplet!
    // Check equality
    if (triplet == ringTriplet) {
      return true;
    } // triplet matches!
    //
    // Check the reversed triplet too
    std::reverse(ringTriplet.begin(), ringTriplet.end());
 
    // Check for equality
    if (triplet == ringTriplet) {
      return true;
    } // reversed triplet matches
  }   // first element of ringTriplet
 
  return false;
}

getEdgeMoleculesInRings()

void ring::getEdgeMoleculesInRings	(	std::vector< std::vector< int > >	rings,
		molSys::PointCloud< molSys::Point< double >, double > *	oCloud,
		molSys::PointCloud< molSys::Point< double >, double > *	yCloud,
		std::array< double, 3 >	coordLow,
		std::array< double, 3 >	coordHigh,
		bool	identicalCloud = false
	)

Select edge molecules and atoms which are part of rings, such that rings formed with even one atom in the slice will be included in the selection Modifies the inSlice bool of a given PointCloud (this may be the same) as the given oxygen atom PointCloud which was used to construct the neighbour list used to construct the rings vector of vectors. We assume that the PointCloud structs have the inSlice bool values set according to the presence of the atom in the slice (this can be done using the gen::moleculesInSingleSlice function.

Function that loops through the PointCloud used to construct the neighbour list (used to generate primitive rings) and sets the inSlice bool values of edge atoms which belong to rings that are formed by atoms in the slice. The output PointCloud may not be the same as the PointCloud used to construct the nList, and the inSlice bool value can be set for both.

Parameters

[in]	rings	Vector of vectors of the primitive rings (by index) according to oCloud
[in]	oCloud	PointCloud of O atoms, used to construct the rings vector of vectors
[in]	yCloud	The output PointCloud (may contain more than just the O atoms)
[in]	identicalCloud	bool value; if this is true then oCloud and yCloud are the same
[in]	coordLow	Contains the lower limits of the slice (needed to find all the molecules of oCloud in the slice)
[in]	coordHigh	Contains the upper limits of the slice

Definition at line 265 of file selection.cpp.

Code

                                                                                      {
  //
  // A vector of bool values, such that every ring has a value of true (in the slice) or false (not in the slice) 
  std::vector<bool> ringInSlice(rings.size(), false); // all set to false initially.
  int jatomIndex, jatomID; // Index and ID in oCloud 
  int jatomIndex1; // Index in yCloud  
  std::unordered_multimap<int, int>
      molIDAtomIDmap; // Unordered multimap with molecule IDs of the atoms as the keys and the
                  // atom IDs as the values. More than one atom can have the same molecule ID
 
  // --------------------
  // Sets all O atoms within the slice to an inSlice bool value of true. If a single atom of a molecule is in the
  // slice, then all atoms in the molecule will also be inside the slice 
  gen::moleculesInSingleSlice(oCloud, true, coordLow, coordHigh);
  // --------------------
  // Get the unordered map of the atom IDs (keys) and the molecular IDs
  // (values)
  molIDAtomIDmap = molSys::createMolIDAtomIDMultiMap(yCloud);
  // --------------------
 
  // Loop through every iatom present in the slice in oCloud
  // Check to see if iatom is in any ring. If it is present in a ring, set that ring to true 
  // Change the inSlice bool of every iatom in the ring to true if false. Set also the inSlice bool
  // of the output cloud if it is not identical. 
 
  // Loop through every iatom present in the slice in oCloud to determine which rings are in the slice
  // The indices of oCloud and rings should match (or this will produce unexpected results)
  for (int iatom = 0; iatom < oCloud->nop; iatom++) {
    // Skip if iatom is not in the slice
    if (!oCloud->pts[iatom].inSlice)
    {
      continue;
    } // skip for iatom not in slice
    //
    // For iatom in the slice, 
    // loop through all rings to find all the rings it is a part of 
    for (int iring = 0; iring < rings.size(); iring++)
    {
      // Skip if iring is in the slice already
      if (ringInSlice[iring])
      {
        continue;
      } // skip for iring in slice 
      // Check and see if iatom is in iring 
      if(std::find(rings[iring].begin(), rings[iring].end(), iatom)!=rings[iring].end()){
        // Found iatom; ring is part of the slice 
        ringInSlice[iring] = true; // update the vector of bool values
        // --------------------------
        // Change the inSlice bool of every iatom in oCloud 
        // (and optionally, yCloud) in the ring to true
        // Loop through the elements of the ring 
        for (int j = 0; j < rings[iring].size(); j++)
        {
          jatomIndex = rings[iring][j]; // Index of the atom in oCloud 
          // Set this to true in oCloud 
          oCloud->pts[jatomIndex].inSlice = true; // part of slice 
          jatomID = oCloud->pts[jatomIndex].atomID; // Atom ID 
          // Now if oCloud and yCloud are not the same, use the
          // atom ID to set the inSlice bool value in yCloud 
          if (!identicalCloud)
          {
            // Find the index corresponding to the same atom in yCloud 
            auto gotJ = yCloud->idIndexMap.find(jatomID);
            jatomIndex1 = gotJ->second;
            // throw if not found ?
            // Set the jatom inSlice bool to true
            yCloud->pts[jatomIndex1].inSlice = true; // jatomIndex is inside the slice 
            // set the inSlice value of all atoms in yCloud with the current molecule ID 
            gen::setAtomsWithSameMolID(yCloud, molIDAtomIDmap, 
              yCloud->pts[jatomIndex1].molID, true);
          } // end of setting values in yCloud 
        } // end of loop through the elements of the current ring 
        // --------------------------
      } // found iatom in the ring
      //
    } // end of loop through all rings searching for iatom 
    //
  } // end of loop through atoms in oCloud in the slice 
 
 
  // Optionally add support for molecule slice update if yCloud and oCloud are not identical?
  // Some other bool?  
  // IC(ringInSlice);
 
  return;
}

getPointCloudOneAtomType()

molSys::PointCloud< molSys::Point< double >, double > gen::getPointCloudOneAtomType	(	molSys::PointCloud< molSys::Point< double >, double > *	yCloud,
		molSys::PointCloud< molSys::Point< double >, double > *	outCloud,
		int	atomTypeI,
		bool	isSlice = false,
		std::array< double, 3 >	coordLow = std::array<double, 3>{0, 0, 0},
		std::array< double, 3 >	coordHigh = std::array<double, 3>{0, 0, 0}
	)

Given a pointCloud containing certain atom types, this returns a pointCloud containing atoms of only the desired type

Function that loops through a given input pointCloud and returns a new pointCloud only containing atoms of a given atom type ID.
This is registered as a Lua function, and is exposed to the user directly.

Parameters

[in]	yCloud	The given input PointCloud
[out]	outCloud	The output PointCloud
[in]	atomTypeI	The type ID of the atoms to save in the output PointCloud
[in]	isSlice	This decides whether a slice will be used or not
[in]	coordLow	Contains the lower limits of the slice, if a slice is to be created
[in]	coordHigh	Contains the upper limits of the slice, if a slice is to be created

Definition at line 37 of file selection.cpp.

Code

                                   {
  //
  int natoms = 0;             // Number of atoms of the desired type 
  bool pointInSlice = true; // If the current point is in the slice, this is true (otherwise this is false)
 
  // --------
  // Before filling up the PointCloud, if the vectors are filled
  // empty them
  *outCloud = molSys::clearPointCloud(outCloud);
  // --------
 
  // Loop through every iatom and check the type
  for (int iatom = 0; iatom < yCloud->nop; iatom++) {
    // Skip if the atom is not of the desired type 
    if (yCloud->pts[iatom].type != atomTypeI) {
      continue;
    } // end of checking for the type 
    // ----------------
    // only if a slice has been requested
    if (isSlice) {
      pointInSlice = sinp::atomInSlice(yCloud->pts[iatom].x, yCloud->pts[iatom].y, yCloud->pts[iatom].z,
                                               coordLow, coordHigh);
      // Skip if the atom is not part of the slice
      if (!pointInSlice) {
        continue;
      } // skipped for atoms not in the slice 
    } // end of slice handling 
    //
    // Actually add the atoms to the outCloud
    natoms++; // Update the number of atoms in outCloud
    outCloud->pts.push_back(yCloud->pts[iatom]); // Update the pts vector
    outCloud->idIndexMap[yCloud->pts[iatom].atomID] =
                outCloud->pts.size() - 1; // array index
  }
 
  // Update the number of particles in the PointCloud
  outCloud->nop = outCloud->pts.size();
 
  // Box and box lengths 
  outCloud->box = yCloud->box;
  outCloud->boxLow = yCloud->boxLow;
 
  // Update the frame number
  outCloud->currentFrame = yCloud->currentFrame; 
 
  return *outCloud;
}

getSingleRingSize()

std::vector< std::vector< int > > ring::getSingleRingSize	(	std::vector< std::vector< int > >	rings,
		int	ringSize
	)

Returns a vector of vectors of rings of a single size.

Function that gets rings of a single ring size (i.e. a particular number of nodes) from all primitive rings, and returns a vector of vectors containing the rings of the specified size.

Parameters

[in]	rings	The vector of vectors containing the primitive rings of all sizes.
[in]	ringSize	The desired ring size or number of nodes in each ring to be saved.

Returns

A vector of vectors containing primitive rings of one ring size or length.

Definition at line 200 of file ring.cpp.

Code

                                                                     {
  //
  std::vector<std::vector<int>> ringSingleSize; // rings of one size
 
  // rings contains primitive rings of all sizes
  // Only save rings of a given size (ringSize) to the new
  // vector of vectors, ringSingleSize
  for (int iring = 0; iring < rings.size(); iring++) {
    // Check the size of the current ring
    // If it is the correct size, save it in ringSingleSize
    if (rings[iring].size() == ringSize) {
      ringSingleSize.push_back(rings[iring]);
    } // End of check of the size of iring
  }   // end of loop through all rings in rings
 
  return ringSingleSize;
}

hasCommonElements()

bool ring::hasCommonElements	(	std::vector< int >	ring1,
		std::vector< int >	ring2
	)

Check to see if two vectors have common elements or not True, if common elements are present and false if there are no common elements

For two vectors, checks to see if there are common elements (true) or not (false).

Parameters

[in]	ring1	The vector of the first ring.
[in]	ring2	The vector of the second ring.

Returns

A bool which is true if the input vectors have at least one common element, and false if there are no common elements.

Definition at line 226 of file ring.cpp.

Code

                                                                       {
  std::vector<int> commonElements; // Vector containing common elements
 
  // Sort the vectors before finding common elements
  sort(ring1.begin(), ring1.end());
  sort(ring2.begin(), ring2.end());
 
  // Find intersection of sorted vectors
  auto it =
      std::set_intersection(ring1.begin(), ring1.end(), ring2.begin(),
                            ring2.end(), std::back_inserter(commonElements));
 
  // If there are no elements in common, then return false
  if (commonElements.size() == 0) {
    return false;
  }
  // If there are common elements, return true
  else {
    return true;
  }
}

keepAxialRingsOnly()

std::vector< std::vector< int > > ring::keepAxialRingsOnly	(	std::vector< std::vector< int > >	rings,
		molSys::PointCloud< molSys::Point< double >, double > *	yCloud
	)

Saves only axial rings out of all possible rings.

moleculesInSingleSlice()

void gen::moleculesInSingleSlice	(	molSys::PointCloud< molSys::Point< double >, double > *	yCloud,
		bool	clearPreviousSliceSelection = true,
		std::array< double, 3 >	coordLow = std::array<double, 3>{0, 0, 0},
		std::array< double, 3 >	coordHigh = std::array<double, 3>{0, 0, 0}
	)

Given a pointCloud set the inSlice bool for every atom, if the molecules are inside the specified (single) region. If even one atom of a molecule is inside the region, then all atoms of that molecule will be inside the region (irrespective of type)

Function that loops through a given input PointCloud and sets the inSlice bool for every Point according to whether the molecule
is in the specified (single) slice or not. If even one atom of a molecule is inside the region, then all atoms belonging to that molecule should be inside the slice as well (therefore, inSlice would be set to true) NOTE: THIS DOES NOT WORK. ERROR

Parameters

[in]	yCloud	The given input PointCloud
[in]	clearPreviousSliceSelection	sets all inSlice bool values to false before adding Points to the slice
[in]	coordLow	Contains the lower limits of the slice, if a slice is to be created
[in]	coordHigh	Contains the upper limits of the slice, if a slice is to be created

Definition at line 147 of file selection.cpp.

Code

                                   {
  //
  bool pointInSlice = true; // If the current point is in the slice, this is true (otherwise this is false)
  std::unordered_multimap<int, int>
      molIDAtomIDmap; // Unordered multimap with molecule IDs of the atoms as the keys and the
                  // atom IDs as the values. More than one atom can have the same molecule ID
  int iatomMolID; // molID of the current iatom
  int jatomID;    // atom ID of the current jatom
  int jatomIndex; // index of jatom
 
  // --------------------
  // Get the unordered map of the atom IDs (keys) and the molecular IDs
  // (values)
  molIDAtomIDmap = molSys::createMolIDAtomIDMultiMap(yCloud);
  // --------------------
  // Set inSlice to false for every atom first
  if (clearPreviousSliceSelection)
  {
    for (int iatom = 0; iatom < yCloud->nop; iatom++) {
      yCloud->pts[iatom].inSlice = false;
    }
  } // end of init
  // --------------------
 
  // Loop through every iatom and check if the atom is in the slice or not 
  for (int iatom = 0; iatom < yCloud->nop; iatom++) {
    // Check if the current atom is in the slice 
    pointInSlice = sinp::atomInSlice(yCloud->pts[iatom].x, yCloud->pts[iatom].y, yCloud->pts[iatom].z,
                                               coordLow, coordHigh);
    //
    // Set the inSlice bool for the particular Point 
    if (pointInSlice)
    {
      yCloud->pts[iatom].inSlice = true; // iatom is inside the slice
      // Find mol ID and atoms with that particular molecular ID 
      iatomMolID = yCloud->pts[iatom].molID; // molecule ID 
      // -----------
      // Find all atoms with iatomMolID and set the inSlice bool
      // to true
      auto range = molIDAtomIDmap.equal_range(iatomMolID); 
      // Loop through all atoms with iatomMolID
      for (auto it = range.first; it != range.second; it++)
      {
        // it->second gives the value (in this case, the atom ID)
        jatomID = it->second; // Atom ID with molecule ID equal to iatomMolID
        auto gotJ = yCloud->idIndexMap.find(jatomID);
        jatomIndex = gotJ->second;
        // Set the jatom inSlice bool to true
        yCloud->pts[jatomIndex].inSlice = true; // jatomIndex is inside the slice 
      }
      // -----------
    } // the atom is in the slice 
    else{
      yCloud->pts[iatom].inSlice = false; // iatom is not in the slice  
    } // atom is not in the slice
  }
 
  return;
}

polygonRingAnalysis()

int ring::polygonRingAnalysis	(	std::string	path,
		std::vector< std::vector< int > >	rings,
		std::vector< std::vector< int > >	nList,
		molSys::PointCloud< molSys::Point< double >, double > *	yCloud,
		int	maxDepth,
		double	sheetArea,
		int	firstFrame
	)

Find out which rings are prisms, looping through all ring sizes upto the maxDepth The input ringsAllSizes array has rings of every size.

Function that loops through the primitive rings (which is a vector of vectors) of all sizes, upto maxDepth (the largest ring size). The function returns a vector which contains the ring indices of all the rings which constitute prism blocks. The ring IDs correspond to the index of the rings inside the vector of vector rings, starting from 0. This is registered as a Lua function, and is exposed to the user directly. The function calls the following functions internally:

• ring::clearRingList (Clears the vector of vectors for rings of a single type, to prevent excessive memory being blocked).
• ring::getSingleRingSize (Fill a vector of vectors for rings of a particular ring size).
• ring::findPrisms (Now that rings of a particular size have been obtained, prism blocks are found and saved).
• topoparam::normHeightPercent (Gets the height% for the prism blocks).
• ring::assignPrismType (Assigns a prism type to each atom type).
• sout::writePrismNum (Write out the prism information for the current frame).
• sout::writeLAMMPSdataAllPrisms (Writes out a LAMMPS data file for the current frame, which can be visualized in OVITO).

  Parameters

  [in]	path	The string to the output directory, in which files will be written out.
  [in]	rings	Row-ordered vector of vectors for rings of a single type.
  [in]	nList	Row-ordered neighbour list by index.
  [in]	yCloud	The input PointCloud.
  [in]	maxDepth	The maximum possible size of the primitive rings.
  [in]	sheetArea	Area calculated using the two significant dimensions of the quasi-two-dimensional sheet.
  [in]	firstFrame	The first frame to be analyzed

Definition at line 37 of file topo_two_dim.cpp.

Code

                                      {
  //
  std::vector<std::vector<int>>
      ringsOneType;           // Vector of vectors of rings of a single size
  int nRings;                 // Number of rings of the current type
  std::vector<int> nRingList; // Vector of the values of the number of prisms
                              // for a particular frame
  std::vector<double> coverageAreaXY; // Coverage area percent on the XY plane
                                      // for a particular n and frame
  std::vector<double> coverageAreaXZ; // Coverage area percent on the XZ plane
                                      // for a particular n and frame
  std::vector<double> coverageAreaYZ; // Coverage area percent on the YZ plane
                                      // for a particular n and frame
  std::vector<double> singleAreas; // Coverage area on the XY, XZ and YZ planes
                                   // for a single ring
  std::vector<int>
      atomTypes; // contains int values for each prism type considered
  // -------------------------------------------------------------------------------
  // Init
  nRingList.resize(
      maxDepth -
      2); // Has a value for every value of ringSize from 3, upto maxDepth
  coverageAreaXY.resize(maxDepth - 2);
  coverageAreaXZ.resize(maxDepth - 2);
  coverageAreaYZ.resize(maxDepth - 2);
  // The atomTypes vector is the same size as the pointCloud atoms
  atomTypes.resize(yCloud->nop, 1); // The dummy or unclassified value is 1
  // -------------------------------------------------------------------------------
  // Run this loop for rings of sizes upto maxDepth
  // The smallest possible ring is of size 3
  for (int ringSize = 3; ringSize <= maxDepth; ringSize++) {
    // Clear ringsOneType
    ring::clearRingList(ringsOneType);
    // Get rings of the current ring size
    ringsOneType = ring::getSingleRingSize(rings, ringSize);
    //
    // Continue if there are zero rings of ringSize
    if (ringsOneType.size() == 0) {
      nRingList[ringSize - 3] = 0; // Update the number of prisms
      // Update the coverage area%
      coverageAreaXY[ringSize - 3] = 0.0;
      coverageAreaXZ[ringSize - 3] = 0.0;
      coverageAreaYZ[ringSize - 3] = 0.0;
      continue;
    } // skip if there are no rings
    //
    // -------------
    // Number of rings with n nodes
    nRings = ringsOneType.size();
    // -------------
    // Now that you have rings of a certain size:
    nRingList[ringSize - 3] = nRings; // Update the number of n-membered rings
    // Update the coverage area% for the n-membered ring phase.
    singleAreas = topoparam::calcCoverageArea(yCloud, ringsOneType, sheetArea);
    // XY plane
    coverageAreaXY[ringSize - 3] = singleAreas[0];
    // XZ plane
    coverageAreaXZ[ringSize - 3] = singleAreas[1];
    // YZ plane
    coverageAreaYZ[ringSize - 3] = singleAreas[2];
    // -------------
  } // end of loop through every possible ringSize
 
  // Get the atom types for all the ring types
  ring::assignPolygonType(rings, &atomTypes, nRingList);
 
  // Write out the ring information
  sout::writeRingNum(path, yCloud->currentFrame, nRingList, coverageAreaXY,
                     coverageAreaXZ, coverageAreaYZ, maxDepth, firstFrame);
  // Write out the lammps data file for the particular frame
  sout::writeLAMMPSdataAllRings(yCloud, nList, atomTypes, maxDepth, path);
 
  return 0;
}

printSliceGetEdgeMoleculesInRings()

void ring::printSliceGetEdgeMoleculesInRings	(	std::string	path,
		std::vector< std::vector< int > >	rings,
		molSys::PointCloud< molSys::Point< double >, double > *	oCloud,
		molSys::PointCloud< molSys::Point< double >, double > *	yCloud,
		std::array< double, 3 >	coordLow,
		std::array< double, 3 >	coordHigh,
		bool	identicalCloud = false
	)

Master function for selecting edge molecules and atoms which are part of rings, such that rings formed with even one atom in the slice will be included in the selection Modifies the inSlice bool of a given PointCloud (this may be the same) as the given oxygen atom PointCloud which was used to construct the neighbour list used to construct the rings vector of vectors (calls ring::getEdgeMoleculesInRings) Prints out molecule IDs individually of molecules in the slice, and also prints out a LAMMPS data file of just the molecules and atoms in the slice

Function that loops through the PointCloud used to construct the neighbour list (used to generate primitive rings) and sets the inSlice bool values of edge atoms which belong to rings that are formed by atoms in the slice. The selected molecule IDs are printed to a separate file, and the molecules and atoms in the slice will be output to a LAMMPS data file (with the original box volume) The output PointCloud may not be the same as the PointCloud used to construct the nList, and the inSlice bool value can be set for both.

Parameters

[in]	rings	Vector of vectors of the primitive rings (by index) according to oCloud
[in]	oCloud	PointCloud of O atoms, used to construct the rings vector of vectors
[in]	yCloud	The output PointCloud (may contain more than just the O atoms)
[in]	identicalCloud	bool value; if this is true then oCloud and yCloud are the same
[in]	coordLow	Contains the lower limits of the slice (needed to find all the molecules of oCloud in the slice)
[in]	coordHigh	Contains the upper limits of the slice

Definition at line 375 of file selection.cpp.

Code

                                                                                      {
  //
 
  //Given the full yCloud PointCloud, set the inSlice bool for every atom,
  // if the molecules are inside the specified (single) region. 
  // gen::atomsInSingleSlice(yCloud, true, coordLow, coordHigh);
  gen::moleculesInSingleSlice(yCloud, true, coordLow, coordHigh);
 
  // Make sure that molecules which participate in the rings inside the slice are also 
  // in the selection 
  ring::getEdgeMoleculesInRings(rings, oCloud, yCloud, coordLow, coordHigh, identicalCloud);
 
  // Print out the molecule IDs of all the atoms in the slice
  sout::writeMoleculeIDsInSlice(path, yCloud);
  // Print out a command that could be used for an expression select command in OVITO
  sout::writeMoleculeIDsExpressionSelectOVITO(path, yCloud);
 
  // Print out the dump of all atoms and molecules, with an inSlice value printed in a separate column
  // H atoms not included in the slice (TODO: fix)
  sout::writeLAMMPSdumpSlice(yCloud, path); 
  
  IC(rings[0]);
 
  return;
}

prismAnalysis()

int ring::prismAnalysis	(	std::string	path,
		std::vector< std::vector< int > >	rings,
		std::vector< std::vector< int > >	nList,
		molSys::PointCloud< molSys::Point< double >, double > *	yCloud,
		int	maxDepth,
		int *	atomID,
		int	firstFrame,
		int	currentFrame,
		bool	doShapeMatching = false
	)

Find out which rings are prisms, looping through all ring sizes upto the maxDepth The input ringsAllSizes array has rings of every size.

Function that loops through the primitive rings (which is a vector of vectors) of all sizes, upto maxDepth (the largest ring size). The function returns a vector which contains the ring indices of all the rings which constitute prism blocks. The ring IDs correspond to the index of the rings inside the vector of vector rings, starting from 0. This is registered as a Lua function, and is exposed to the user directly. The function calls the following functions internally:

• ring::clearRingList (Clears the vector of vectors for rings of a single type, to prevent excessive memory being blocked).
• ring::getSingleRingSize (Fill a vector of vectors for rings of a particular ring size).
• ring::findPrisms (Now that rings of a particular size have been obtained, prism blocks are found and saved).
• topoparam::normHeightPercent (Gets the height% for the prism blocks).
• ring::assignPrismType (Assigns a prism type to each atom type).
• sout::writePrismNum (Write out the prism information for the current frame).
• sout::writeLAMMPSdataAllPrisms (Writes out a LAMMPS data file for the current frame, which can be visualized in OVITO).

  Parameters

  [in]	path	The string to the output directory, in which files will be written out.
  [in]	rings	Row-ordered vector of vectors for rings of a single type.
  [in]	nList	Row-ordered neighbour list by index.
  [in]	yCloud	The input PointCloud.
  [in]	maxDepth	The maximum possible size of the primitive rings.
  [in]	maxDepth	The first frame.

Definition at line 49 of file topo_one_dim.cpp.

Code

                                                                         {
  //
  std::vector<std::vector<int>>
      ringsOneType;           // Vector of vectors of rings of a single size
  std::vector<int> listPrism; // Vector for ring indices of n-sided prism
  std::vector<ring::strucType>
      ringType; // This vector will have a value for each ring inside
  std::vector<ring::strucType>
      atomState; // This vector will have a value for each atom, depending on
                 // the ring type
  int nPerfectPrisms;          // Number of perfect prisms of each type
  int nImperfectPrisms;        // Number of deformed prisms of each type
  std::vector<int> nPrismList; // Vector of the values of the number of perfect
                               // prisms for a particular frame
  std::vector<int> nDefPrismList; // Vector of the values of the number of
                                  // deformed prisms for a particular frame
  std::vector<double>
      heightPercent; // Height percent for a particular n and frame
  std::vector<int>
      atomTypes; // contains int values for each prism type considered
  double avgPrismHeight = 2.845; // A value of 2.7-2.85 Angstrom is reasonable
  // Qualifier for the RMSD per atom:
  std::vector<double> rmsdPerAtom;
  // -------------------------------------------------------------------------------
  // Init
  nPrismList.resize(
      maxDepth -
      2); // Has a value for every value of ringSize from 3, upto maxDepth
  nDefPrismList.resize(maxDepth - 2);
  heightPercent.resize(maxDepth - 2);
  // The atomTypes vector is the same size as the pointCloud atoms
  atomTypes.resize(yCloud->nop, 1); // The dummy or unclassified value is 1
  // The rmsdPerAtom vector is the same size as the pointCloud atoms and has an
  // RMSD value for every atom
  rmsdPerAtom.resize(yCloud->nop, -1);
  // Resize the atom state vector
  atomState.resize(yCloud->nop); // Dummy or unclassified
  // -------------------------------------------------------------------------------
  // Run this loop for rings of sizes upto maxDepth
  // The smallest possible ring is of size 3
  for (int ringSize = 3; ringSize <= maxDepth; ringSize++) {
    // Clear ringsOneType
    ring::clearRingList(ringsOneType);
    // Get rings of the current ring size
    ringsOneType = ring::getSingleRingSize(rings, ringSize);
    //
    // Continue if there are zero rings of ringSize
    if (ringsOneType.size() == 0) {
      nPrismList[ringSize - 3] = 0;      // Update the number of prisms
      nDefPrismList[ringSize - 3] = 0;   // Update the number of deformed prisms
      heightPercent[ringSize - 3] = 0.0; // Update the height%
      continue;
    } // skip if there are no rings
    //
    // -------------
    // Init of variables specific to ringSize prisms
    listPrism.resize(0);
    ringType.resize(0);
    nPerfectPrisms = 0;
    nImperfectPrisms = 0;
    ringType.resize(
        ringsOneType.size()); // Has a value for each ring. init to zero.
    // -------------
    // Now that you have rings of a certain size:
    // Find prisms, saving the ring indices to listPrism
    listPrism = ring::findPrisms(ringsOneType, &ringType, &nPerfectPrisms,
                                 &nImperfectPrisms, nList, yCloud, &rmsdPerAtom,
                                 doShapeMatching);
    // -------------
    nPrismList[ringSize - 3] =
        nPerfectPrisms; // Update the number of perfect prisms
    nDefPrismList[ringSize - 3] =
        nImperfectPrisms; // Update the number of defromed prisms
    int totalPrisms = nPerfectPrisms + nImperfectPrisms;
    // Update the height% for the phase
    heightPercent[ringSize - 3] =
        topoparam::normHeightPercent(yCloud, totalPrisms, avgPrismHeight);
    // Continue if there are no prism units
    if (nPerfectPrisms + nImperfectPrisms == 0) {
      continue;
    } // skip for no prisms
    // Do a bunch of write-outs and calculations
    // TODO: Write out each individual prism as data files (maybe with an
    // option)
    // Get the atom types for a particular prism type
    ring::assignPrismType(ringsOneType, listPrism, ringSize, ringType,
                          &atomTypes, &atomState);
    // -------------
  } // end of loop through every possible ringSize
 
  // Calculate the height%
 
  // Write out the prism information
  sout::writePrismNum(path, nPrismList, nDefPrismList, heightPercent, maxDepth,
                      yCloud->currentFrame, firstFrame);
  // Reassign the prism block types for the deformed prisms
  if (doShapeMatching) {
    ring::deformedPrismTypes(atomState, &atomTypes, maxDepth);
  } // reassign prism block types for deformed prisms
 
  // Get rid of translations along the axial dimension
  ring::rmAxialTranslations(yCloud, atomID, firstFrame, currentFrame);
 
  // Write out dump files with the RMSD per atom for the shape-matching
  // criterion
  if (doShapeMatching) {
    sout::writeLAMMPSdumpINT(yCloud, rmsdPerAtom, atomTypes, maxDepth, path);
  } // reassign prism block types for deformed prisms
 
  // Write out the lammps data file for the particular frame
  sout::writeLAMMPSdataAllPrisms(yCloud, nList, atomTypes, maxDepth, path,
                                 doShapeMatching);
 
  return 0;
}

relaxedPrismConditions()

bool ring::relaxedPrismConditions	(	std::vector< std::vector< int > >	nList,
		std::vector< int > *	basal1,
		std::vector< int > *	basal2
	)

Reduced criterion: Two candidate basal rings of a prism block should have at least one bond between them

Relaxed criteria for deformed prism blocks: at least one bond should exist between the basal rings.

Definition at line 436 of file topo_one_dim.cpp.

Code

                                                          {
  int ringSize =
      (*basal1).size();     // Size of the ring; each ring contains n elements
  int m_k;                  // Atom ID of element in basal2
  bool isNeighbour = false; // This is true if there is at least one bond
                            // between the basal rings
  int l_k;                  // Atom ID of element in basal1
 
  // ---------------------------------------------
  // COMPARISON OF basal2 ELEMENTS (m_k) WITH basal1 ELEMENTS (l_k)
  // Loop through all the elements of basal1
  for (int l = 0; l < ringSize; l++) {
    l_k = (*basal1)[l];
    // Search for the nearest neighbour of l_k in basal2
    // Loop through basal2 elements
    for (int m = 0; m < ringSize; m++) {
      m_k = (*basal2)[m];
      // Find m_k inside l_k neighbour list
      auto it = std::find(nList[l_k].begin() + 1, nList[l_k].end(), m_k);
 
      // If the element has been found, for l1
      if (it != nList[l_k].end()) {
        isNeighbour = true;
        break;
      } // found element
    }   // end of loop through all the elements of basal2
 
    // If a neighbour has been found then
    if (isNeighbour) {
      return true;
    }
  } // end of loop through all the elements of basal1
 
  // If a neighbour has not been found, return false
  return false;
}

rmAxialTranslations()

int ring::rmAxialTranslations	(	molSys::PointCloud< molSys::Point< double >, double > *	yCloud,
		int *	atomID,
		int	firstFrame,
		int	currentFrame
	)

Shift the entire ice nanotube and remove axial translations.

Assign an atomType value for atoms belonging to deformed prisms.

Parameters

[in]	rings	The vector of vectors containing the primitive rings, of a particular ring size.
[in]	listPrism	The list of prism blocks found.
[in]	ringSize	The current ring size or number of nodes in each ring.
[in,out]	atomTypes	A vector which contains a type for each atom, depending on it's type as classified by the prism identification scheme.

Definition at line 791 of file topo_one_dim.cpp.

Code

                                      {
  //
  int atomIndex;          // Index of the atom to be shifted first
  double shiftDistance;   // Value by which all atoms will be shifted
  double distFrmBoundary; // Distance from either the upper or lower boundary
                          // along the axial dimension
  // Get the axial dimension
  int axialDim = std::max_element(yCloud->box.begin(), yCloud->box.end()) -
                 yCloud->box.begin();
  // Check: (default z)
  if (axialDim < 0 || axialDim > 2) {
    axialDim = 2;
  } // end of check to set the axial dimension
  // Lower and higher limits of the box in the axial dimension
  double boxLowAxial = yCloud->boxLow[axialDim];
  double boxHighAxial = boxLowAxial + yCloud->box[axialDim];
  //
  // -----------------------------------
  // Update the atomID of the 'first' atom in yCloud if the current frame is the
  // first frame.
  // Get the atom index of the first atom to be shifted down or up
  if (currentFrame == firstFrame) {
    *atomID = yCloud->pts[0].atomID;
    atomIndex = 0;
  } // end of update of the atom ID to be shifted for the first frame
  else {
    //
    int iatomID = *atomID; // Atom ID of the first particle to be moved down
    // Find the index given the atom ID
    auto it = yCloud->idIndexMap.find(iatomID);
 
    if (it != yCloud->idIndexMap.end()) {
      atomIndex = it->second;
    } // found jatom
    else {
      std::cerr << "Lost atoms.\n";
      return 1;
    } // error handling
    //
  } // Update of atomIndex for all other frames
  // -----------------------------------
  // Calculate the distance by which the atoms must be shifted (negative value)
  if (axialDim == 0) {
    shiftDistance = boxLowAxial - yCloud->pts[atomIndex].x;
  } // x coordinate
  else if (axialDim == 1) {
    shiftDistance = boxLowAxial - yCloud->pts[atomIndex].y;
  } // y coordinate
  else {
    shiftDistance = boxLowAxial - yCloud->pts[atomIndex].z;
  } // z coordinate
  // -----------------------------------
  // Shift all the particles
  for (int iatom = 0; iatom < yCloud->nop; iatom++) {
    // Shift the particles along the axial dimension only
    if (axialDim == 0) {
      yCloud->pts[iatom].x += shiftDistance;
      // Wrap-around
      if (yCloud->pts[iatom].x < boxLowAxial) {
        distFrmBoundary = boxLowAxial - yCloud->pts[iatom].x; // positive value
        yCloud->pts[iatom].x = boxHighAxial - distFrmBoundary;
      } // end of wrap-around
    }   // x coordinate
    else if (axialDim == 1) {
      yCloud->pts[iatom].y += shiftDistance;
      // Wrap-around
      if (yCloud->pts[iatom].y < boxLowAxial) {
        distFrmBoundary = boxLowAxial - yCloud->pts[iatom].y; // positive value
        yCloud->pts[iatom].y = boxHighAxial - distFrmBoundary;
      } // end of wrap-around
    }   // y coordinate
    else {
      yCloud->pts[iatom].z += shiftDistance;
      // Wrap-around
      if (yCloud->pts[iatom].z < boxLowAxial) {
        distFrmBoundary = boxLowAxial - yCloud->pts[iatom].z; // positive value
        yCloud->pts[iatom].z = boxHighAxial - distFrmBoundary;
      } // end of wrap-around
    }   // z coordinate
  }     // end of shifting all paritcles downward
  // -----------------------------------
  return 0;
} // end of function

setAtomsWithSameMolID()

void gen::setAtomsWithSameMolID	(	molSys::PointCloud< molSys::Point< double >, double > *	yCloud,
		std::unordered_multimap< int, int >	molIDAtomIDmap,
		int	molID,
		bool	inSliceValue = true
	)

Given a particular molecule ID and a pointCloud set the inSlice bool for all atoms, with that molecule ID

Function that loops through a given input PointCloud and sets the inSlice bool for every Point according to whether the molecule
is in the specified (single) slice or not. If even one atom of a molecule is inside the region, then all atoms belonging to that molecule should be inside the slice as well (therefore, inSlice would be set to true)

Parameters

[in]	yCloud	The given input PointCloud
[in]	clearPreviousSliceSelection	sets all inSlice bool values to false before adding Points to the slice
[in]	coordLow	Contains the lower limits of the slice, if a slice is to be created
[in]	coordHigh	Contains the upper limits of the slice, if a slice is to be created

Definition at line 225 of file selection.cpp.

Code

                                  {
  //
  int jatomID;    // atom ID of the current jatom
  int jatomIndex; // index of jatom
 
  // Find all atoms with molID and set the inSlice bool to inSliceValue
  auto range = molIDAtomIDmap.equal_range(molID); 
 
  // Loop through all atoms with molID
  for (auto it = range.first; it != range.second; it++){
    // it->second gives the value (in this case, the atom ID)
    jatomID = it->second; // Atom ID with molecule ID equal to iatomMolID
    auto gotJ = yCloud->idIndexMap.find(jatomID);
    jatomIndex = gotJ->second;
    // Set the jatom inSlice bool to true
    yCloud->pts[jatomIndex].inSlice = inSliceValue; // jatomIndex is assigned inSliceValue
  } // end of loop through all atoms with molID
 
  // IC(yCloud->pts[jatomIndex]);
 
  return;
}

shapeMatchDDC()

int tum3::shapeMatchDDC	(	molSys::PointCloud< molSys::Point< double >, double > *	yCloud,
		const Eigen::MatrixXd &	refPoints,
		std::vector< cage::Cage >	cageList,
		int	cageIndex,
		std::vector< std::vector< int > >	rings,
		std::vector< double > *	quat,
		double *	rmsd
	)

Shape-matching for a target DDC.

Match the input cage with a perfect DDC. The quaternion for the rotation and the RMSD is outputted.

Parameters

[in]	refPoints	The point set of the reference system (or right system). This is a $(n\times 3)$ Eigen matrix. Here, $n$ is the number of particles.
[in]	targetPoints	$(n\times 3)$ Eigen matrix of the candidate/test system (or left system).
[in,out]	quat	The quaternion for the optimum rotation of the left system into the right system.
[in,out]	scale	The scale factor obtained from Horn's algorithm.

Definition at line 323 of file bulkTUM.cpp.

Code

                                           {
  //
  std::vector<int> ddcOrder;             // Connectivity of the DDC
  int ringSize = 6;                      // Each ring has 6 nodes
  Eigen::MatrixXd targetPointSet(14, 3); // Target point set (Eigen matrix)
  //
  std::vector<double> rmsdList; // List of RMSD per atom
  // Variables for looping through possible permutations
  //
  std::vector<double> currentQuat; // quaternion rotation
  double currentRmsd;              // least RMSD
  double currentScale;             // scale
  double scale;                    // Final scale
  int index; // Int for describing which permutation matches the best
 
  // ----------------
  // Save the order of the DDC in a vector
  ddcOrder = pntToPnt::relOrderDDC(cageIndex, rings, cageList);
  // ----------------
  // Loop through all possible permutations
  //
  for (int i = 0; i < ringSize; i++) {
    // Change the order of the target points somehow!
    //
    targetPointSet = pntToPnt::changeDiaCageOrder(yCloud, ddcOrder, i);
    // Shape-matching
    absor::hornAbsOrientation(refPoints, targetPointSet, &currentQuat,
                              &currentRmsd, &rmsdList, &currentScale);
    if (i == 0) {
      *quat = currentQuat;
      *rmsd = currentRmsd;
      scale = currentScale;
      index = 0;
    } else {
      if (currentRmsd < *rmsd) {
        *quat = currentQuat;
        *rmsd = currentRmsd;
        scale = currentScale;
        index = i;
      } // update
    }   // Update if this is a better match
  }     // Loop through possible permutations
  // ----------------
  return 0;
}

shapeMatchHC()

int tum3::shapeMatchHC	(	molSys::PointCloud< molSys::Point< double >, double > *	yCloud,
		const Eigen::MatrixXd &	refPoints,
		cage::Cage	cageUnit,
		std::vector< std::vector< int > >	rings,
		std::vector< std::vector< int > >	nList,
		std::vector< double > *	quat,
		double *	rmsd
	)

Shape-matching for a target HC.

Match the input cage with a perfect HC. The quaternion for the rotation and the RMSD is outputted.

Parameters

[in]	refPoints	The point set of the reference system (or right system). This is a $(n\times 3)$ Eigen matrix. Here, $n$ is the number of particles.
[in]	targetPoints	$(n\times 3)$ Eigen matrix of the candidate/test system (or left system).
[in,out]	quat	The quaternion for the optimum rotation of the left system into the right system.
[in,out]	scale	The scale factor obtained from Horn's algorithm.

Definition at line 251 of file bulkTUM.cpp.

Code

                                           {
  //
  int iring,
      jring; // Indices in the ring vector of vectors for the basal rings
  std::vector<int> basal1,
      basal2;       // Re-ordered basal rings matched to each other
  int ringSize = 6; // Each ring has 6 nodes
  Eigen::MatrixXd targetPointSet(12, 3); // Target point set (Eigen matrix)
  //
  std::vector<double> rmsdList; // List of RMSD per atom
  // Variables for looping through possible permutations
  //
  std::vector<double> currentQuat; // quaternion rotation
  double currentRmsd;              // least RMSD
  double currentScale;             // scale
  double scale;                    // Final scale
  int index; // Int for describing which permutation matches the best
 
  // Init
  iring = cageUnit.rings[0];    // Index of basal1
  jring = cageUnit.rings[1];    // Index of basal2
  rmsdList.resize(yCloud->nop); // Not actually updated here
  //
  // ----------------
  // Re-order the basal rings so that they are matched
  pntToPnt::relOrderHC(yCloud, rings[iring], rings[jring], nList, &basal1,
                       &basal2);
  // ----------------
  // Loop through all possible permutations
  //
  for (int i = 0; i < ringSize; i++) {
    // Change the order of the target points somehow!
    //
    targetPointSet = pntToPnt::changeHexCageOrder(yCloud, basal1, basal2, i);
    // Shape-matching
    absor::hornAbsOrientation(refPoints, targetPointSet, &currentQuat,
                              &currentRmsd, &rmsdList, &currentScale);
    if (i == 0) {
      *quat = currentQuat;
      *rmsd = currentRmsd;
      scale = currentScale;
      index = 0;
    } else {
      if (currentRmsd < *rmsd) {
        *quat = currentQuat;
        *rmsd = currentRmsd;
        scale = currentScale;
        index = i;
      } // update
    }   // Update if this is a better match
  }     // Loop through possible permutations
  // ----------------
  return 0;
}

topoBulkCriteria()

std::vector< cage::Cage > tum3::topoBulkCriteria	(	std::string	path,
		std::vector< std::vector< int > >	rings,
		std::vector< std::vector< int > >	nList,
		molSys::PointCloud< molSys::Point< double >, double > *	yCloud,
		int	firstFrame,
		int *	numHC,
		int *	numDDC,
		std::vector< ring::strucType > *	ringType
	)

Topological network methods Finds the HCs and DDCs for the system

Finds out if rings constitute double-diamond cages or hexagonal cages. Requires a neighbour list (by index) and a vector of vectors of primitive rings which should also be by index. This is only for rings of size 6! Internally, this function calls the following functions:

• ring::getSingleRingSize (Saves rings of a single ring size into a new vector of vectors, which is subsequently used for finding DDCs, HCs etc).
• ring::findDDC (Finds the DDCs).
• ring::findHC (Finds the HCs).
• ring::findMixedRings (Finds the mixed rings, which are shared by DDCs and HCs both).
• ring::getStrucNumbers (Gets the number of structures (DDCs, HCs, mixed rings, basal rings, prismatic rings, to be used for write-outs).
• sout::writeTopoBulkData (Writes out the numbers and data obtained for the current frame).
• ring::getAtomTypesTopoBulk (Gets the atom type for every atom, to be used for printing out the ice types found).
• sout::writeLAMMPSdataTopoBulk (Writes out the atoms, with the classified types, into a LAMMPS data file, which can be visualized in OVITO).

  Parameters

  [in]	path	The file path of the output directory to which output files will be written.
  [in]	rings	Vector of vectors containing the primitive rings. This should contain rings of only size 6.
  [in]	nList	Row-ordered neighbour list, by index.
  [in]	yCloud	The input PointCloud, with respect to which the indices in the rings and nList vector of vectors have been saved.
  [in]	firstFrame	First frame to be analyzed
  [in]	onlyTetrahedral	Flag for only finding DDCs and HCs (true) or also finding PNCs (false)

Definition at line 581 of file bulkTUM.cpp.

Code

                                                                 {
  //
  // Ring IDs of each type will be saved in these vectors
  std::vector<int> listDDC; // Vector for ring indices of DDC
  std::vector<int> listHC;  // Vector for ring indices of HC
  std::vector<int>
      listMixed; // Vector for ring indices of rings that are both DDC and HC
                 // ringList, of type enum strucType in gen.hpp
  // Make a list of all the DDCs and HCs
  std::vector<cage::Cage> cageList;
  // Number of types
  int mixedRings, prismaticRings, basalRings;
 
  // ----------------------------------------------
  // Init
  //
  *numHC = 0;  // Number of hexagonal cages
  *numDDC = 0; // Init the number of DDCs
  // Quit the function for zero rings
  if (rings.size() == 0) {
    return cageList;
  } // skip for no rings of ringSize
  //
  // Init the ringType vector
  (*ringType).resize(rings.size()); // Has a value for each ring. init to zero.
  // ----------------------------------------------
  // Get the cages
 
  // Find HC rings, saving the ring IDs (starting from 0) to listHC
  listHC = ring::findHC(rings, ringType, nList, &cageList);
 
  // Find DDC rings, saving the IDs to listDDC
  listDDC = ring::findDDC(rings, ringType, listHC, &cageList);
 
  // Find rings which are both DDCs and HCs (mixed)
  // A dummy value of -10 in the listDDC and listHC vectors for mixed rings
  listMixed = ring::findMixedRings(rings, ringType, &listDDC, &listHC);
 
  // Get the number of structures (DDCs, HCs, mixed rings, basal rings,
  // prismatic rings)
  ring::getStrucNumbers(*ringType, cageList, numHC, numDDC, &mixedRings,
                        &prismaticRings, &basalRings);
 
  // Write out to a file
  sout::writeTopoBulkData(path, yCloud->currentFrame, *numHC, *numDDC,
                          mixedRings, basalRings, prismaticRings, firstFrame);
 
  return cageList;
}

topoUnitMatchingBulk()

int tum3::topoUnitMatchingBulk	(	std::string	path,
		std::vector< std::vector< int > >	rings,
		std::vector< std::vector< int > >	nList,
		molSys::PointCloud< molSys::Point< double >, double > *	yCloud,
		int	firstFrame,
		bool	printClusters,
		bool	onlyTetrahedral
	)

Topological unit matching for bulk water. If printClusters is true, individual clusters of connected cages are printed.

Finds out if rings constitute double-diamond cages or hexagonal cages. Requires a neighbour list (by index) and a vector of vectors of primitive rings which should also be by index. This is registered as a Lua function and is accessible to the user. Internally, this function calls the following functions:

• ring::getSingleRingSize (Saves rings of a single ring size into a new vector of vectors, which is subsequently used for finding DDCs, HCs etc).
• ring::findDDC (Finds the DDCs).
• ring::findHC (Finds the HCs).
• ring::findMixedRings (Finds the mixed rings, which are shared by DDCs and HCs both).
• ring::getStrucNumbers (Gets the number of structures (DDCs, HCs, mixed rings, basal rings, prismatic rings, to be used for write-outs).
• sout::writeTopoBulkData (Writes out the numbers and data obtained for the current frame).
• ring::getAtomTypesTopoBulk (Gets the atom type for every atom, to be used for printing out the ice types found).
• sout::writeLAMMPSdataTopoBulk (Writes out the atoms, with the classified types, into a LAMMPS data file, which can be visualized in OVITO).

  Parameters

  [in]	path	The file path of the output directory to which output files will be written.
  [in]	rings	Vector of vectors containing the primitive rings. This contains rings of all sizes.
  [in]	nList	Row-ordered neighbour list, by index.
  [in]	yCloud	The input PointCloud, with respect to which the indices in the rings and nList vector of vectors have been saved.
  [in]	firstFrame	First frame to be analyzed
  [in]	onlyTetrahedral	Flag for only finding DDCs and HCs (true) or also finding PNCs (false)

Definition at line 52 of file bulkTUM.cpp.

Code

                                              {
  // The input rings vector has rings of all sizes
 
  // ringType has a value for rings of a particular size
  std::vector<ring::strucType>
      ringType; // This vector will have a value for each ring inside
                // ringList, of type enum strucType in gen.hpp
  // Make a list of all the DDCs and HCs
  std::vector<cage::Cage> cageList;
  std::vector<std::vector<int>>
      ringsOneType;    // Vector of vectors of rings of a single size
  int initRingSize;    // Todo or not: calculate the PNCs or not
  int maxRingSize = 6; // DDCs and HCs are for 6-membered rings
  std::vector<cage::iceType>
      atomTypes; // This vector will have a value for every atom
  // Number of types
  int numHC, numDDC, mixedRings, prismaticRings, basalRings;
  // Shape-matching variables ----
  double rmsd;              // RMSD value for a particular cage type
  std::vector<double> quat; // Quaternion obtained from shape-matching
  std::vector<std::vector<double>> quatList; // List of quaternions
  // Reference points
  Eigen::MatrixXd refPntsDDC(14,
                             3); // Reference point set (Eigen matrix) for a DDC
  Eigen::MatrixXd refPntsHC(12,
                            3); // Reference point set (Eigen matrix) for a DDC
 
  // Vector for the RMSD per atom and the RMSD per ring:
  std::vector<double> rmsdPerAtom;
  std::vector<int>
      noOfCommonElements; // An atom can belong to more than one ring or shape
  //
 
  if (onlyTetrahedral) {
    initRingSize = 6;
  } else {
    initRingSize = 5;
  }
 
  // Init the atom type vector
  atomTypes.resize(yCloud->nop); // Has a value for each atom
  // Init the rmsd per atom
  rmsdPerAtom.resize(yCloud->nop, -1);    // Has a value for each atom
  noOfCommonElements.resize(yCloud->nop); // Has a value for each atom
 
  // ----------------------------------------------
  // Init
  // Get rings of size 5 or 6.
  for (int ringSize = initRingSize; ringSize <= maxRingSize; ringSize++) {
    // Clear ringsOneType
    ring::clearRingList(ringsOneType);
    // Get rings of the current ring size
    ringsOneType = ring::getSingleRingSize(rings, ringSize);
    // Skip for zero rings
    if (ringsOneType.size() == 0) {
      continue;
    } // skip for no rings of ringSize
    //
    // Init the ringType vector
    ringType.resize(
        ringsOneType.size()); // Has a value for each ring. init to zero.
    // ----------------------------------------------
    if (ringSize == 6) {
      // Get the cages
 
      // Find the DDCs and HCs
      cageList = tum3::topoBulkCriteria(path, ringsOneType, nList, yCloud,
                                        firstFrame, &numHC, &numDDC, &ringType);
 
      // Gets the atom type for every atom, to be used for printing out the ice
      // types found
      ring::getAtomTypesTopoBulk(ringsOneType, ringType, &atomTypes);
    }
    // Finding the 5-membered prismatic blocks
    else {
      // Get the prism block classifications
      prism3::findBulkPrisms(ringsOneType, &ringType, nList, yCloud,
                             &rmsdPerAtom);
      // Gets the atom type for every atom, to be used for printing out the ice
      // types found
      ring::getAtomTypesTopoBulk(ringsOneType, ringType, &atomTypes);
    }
  }
 
  // --------------------------------------------------
  // SHAPE-MATCHING
  //
  // Get the reference point sets
  //
  std::string filePathHC = "templates/hc.xyz";
  std::string filePathDDC = "templates/ddc.xyz";
  refPntsHC = tum3::buildRefHC(filePathHC);    // HC
  refPntsDDC = tum3::buildRefDDC(filePathDDC); // DDC
  //
  // Init
  //
  // Loop through all the atoms and set commonElements to 1 for PNC atoms
  for (int i = 0; i < yCloud->nop; i++) {
    if (rmsdPerAtom[i] != -1) {
      noOfCommonElements[i] = 1;
    } // for atoms which have been updated
  }   // end of updating commonElements
  //
  // Loop through the entire cageList vector of cages to match the HCs and DDCs
  //
  // HCs are first, followed by DDCs.
  //
  // Go through all the HCs
  for (int icage = 0; icage < numHC; icage++) {
    // Match against a perfect HC
    tum3::shapeMatchHC(yCloud, refPntsHC, cageList[icage], ringsOneType, nList,
                       &quat, &rmsd);
    // Update the vector of quaternions
    quatList.push_back(quat);
    // Update the RMSD per ring
    tum3::updateRMSDatom(ringsOneType, cageList[icage], rmsd, &rmsdPerAtom,
                         &noOfCommonElements, atomTypes);
  } // end of looping through all HCs
  // --------------------------------------------------
  // Go through all the DDCs
  for (int icage = numHC; icage < cageList.size(); icage++) {
    // Match against a perfect DDC
    tum3::shapeMatchDDC(yCloud, refPntsDDC, cageList, icage, ringsOneType,
                        &quat, &rmsd);
    // Update the vector of quaternions
    quatList.push_back(quat);
    // Update the RMSD per ring
    tum3::updateRMSDatom(ringsOneType, cageList[icage], rmsd, &rmsdPerAtom,
                         &noOfCommonElements, atomTypes);
  } // end of looping through all HCs
 
  // --------------------------------------------------
  // Getting the RMSD per atom
  // Average the RMSD per atom
  tum3::averageRMSDatom(&rmsdPerAtom, &noOfCommonElements);
  // --------------------------------------------------
 
  // Print out the lammps data file with the bonds and types
  sout::writeLAMMPSdataTopoBulk(yCloud, nList, atomTypes, path);
  // To output the bonds between dummy atoms, uncomment the following line
  // sout::writeLAMMPSdataTopoBulk(yCloud, nList, atomTypes, path, true);
 
  // --------------------------------------------------
  // Writes out the dumpfile
  //
  std::vector<int> dumpAtomTypes; // Must be typecast to int
  dumpAtomTypes.resize(atomTypes.size());
  // Change from enum to int
  for (int i = 0; i < atomTypes.size(); i++) {
    if (atomTypes[i] == cage::iceType::hc) {
      dumpAtomTypes[i] = 1;
    } // HC
    else if (atomTypes[i] == cage::iceType::ddc) {
      dumpAtomTypes[i] = 2;
    } // DDC
    else if (atomTypes[i] == cage::iceType::mixed) {
      dumpAtomTypes[i] = 3;
    } // mixed rings
    else if (atomTypes[i] == cage::iceType::pnc) {
      dumpAtomTypes[i] = 4;
    } // pnc
    else if (atomTypes[i] == cage::iceType::mixed2) {
      dumpAtomTypes[i] = 5;
    } // shared by pnc and ddc/hc
    else {
      dumpAtomTypes[i] = 0;
    } // dummy atoms
  }   // end of changing the type to int
  sout::writeLAMMPSdumpCages(yCloud, rmsdPerAtom, dumpAtomTypes, path,
                             firstFrame);
  // --------------------------------------------------
  //
  // PRINT THE CLUSTERS
  //
  if (printClusters) {
    // Print the clusters
    tum3::clusterCages(yCloud, path, ringsOneType, cageList, numHC, numDDC);
  } // end of printing individual clusters
  // --------------------------------------------------
  return 0;
}

updateRMSDatom()

int tum3::updateRMSDatom	(	std::vector< std::vector< int > >	rings,
		cage::Cage	cageUnit,
		double	rmsd,
		std::vector< double > *	rmsdPerAtom,
		std::vector< int > *	noOfCommonAtoms,
		std::vector< cage::iceType >	atomTypes
	)

Calulate the RMSD for each ring, using RMSD values (rmsd) obtained from the shape-matching of each cage

Update the calculated RMSD per ring using the RMSD values of each cage, and also update the values in the noOfCommonRings vector, which will be used for averaging the RMSD per atom depending on the number of cages that share that particular ring.

Definition at line 486 of file bulkTUM.cpp.

Code

                                                           {
  //
  int nRings = cageUnit.rings.size(); // Number of rings in the current cage
  int iring; // Index according to the rings vector of vector, for the current
             // ring inside the cage
  int iatom; // Current atom index
  int ringSize = rings[0].size(); // Number of nodes in each ring (6)
 
  // Loop through the rings in each cage
  for (int i = 0; i < nRings; i++) {
    iring = cageUnit.rings[i]; // Current ring
 
    // Loop through all the atoms in the ring
    for (int j = 0; j < ringSize; j++) {
      iatom = rings[iring][j]; // Current atom index
 
      // Skip for PNC atoms
      if (atomTypes[iatom] == cage::iceType::pnc || atomTypes[iatom] == cage::iceType::mixed2) {
        continue;
      } // Do not update if the atom is a PNC
      //
      // UPDATE
      if ((*rmsdPerAtom)[iatom] == -1) {
        (*rmsdPerAtom)[iatom] = rmsd;
        (*noOfCommonAtoms)[iatom] = 1;
      } else {
        (*rmsdPerAtom)[iatom] += rmsd;
        (*noOfCommonAtoms)[iatom] += 1;
      }
    } // end of looping through all the atoms in the ring
 
  } // end of looping through the rings in the cage
 
  return 0;
}