Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Copy bonds and bond orders when generating super cells #1898

Merged
merged 2 commits into from
Dec 28, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
29 changes: 22 additions & 7 deletions avogadro/core/crystaltools.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@ struct WrapAtomsToCellFunctor

void operator()(Vector3& pos) { unitCell.wrapCartesian(pos, pos); }
};
}
} // namespace

bool CrystalTools::wrapAtomsToUnitCell(Molecule& molecule)
{
Expand Down Expand Up @@ -182,7 +182,7 @@ T niggliRound(T v, T dec)
const T shifted = v * shift;
return std::floor(shifted + 0.5) / shift;
}
}
} // namespace

bool CrystalTools::niggliReduce(Molecule& molecule, Options opts)
{
Expand Down Expand Up @@ -433,7 +433,7 @@ bool CrystalTools::niggliReduce(Molecule& molecule, Options opts)

// fix coordinates with COB matrix:
const Matrix3 invCob(cob.inverse());
for (auto & fcoord : fcoords) {
for (auto& fcoord : fcoords) {
fcoord = invCob * fcoord;
}

Expand Down Expand Up @@ -558,6 +558,10 @@ bool CrystalTools::buildSupercell(Molecule& molecule, unsigned int a,
Vector3 newB = oldB * b;
Vector3 newC = oldC * c;

// archive the old bond pairs and orders
Array<std::pair<Index, Index>> bondPairs = molecule.bondPairs();
Array<unsigned char> bondOrders = molecule.bondOrders();

// Add in the atoms to the new subcells of the supercell
Index numAtoms = molecule.atomCount();
Array<Vector3> atoms = molecule.atomPositions3d();
Expand All @@ -578,6 +582,17 @@ bool CrystalTools::buildSupercell(Molecule& molecule, unsigned int a,
}
}

// now we need to add the bonds
unsigned copies = molecule.atomCount() / numAtoms;
// we loop through the original bonds to add copies
for (Index i = 0; i < bondPairs.size(); ++i) {
std::pair<Index, Index> bond = bondPairs.at(i);
for (unsigned j = 0; j < copies; ++j) {
molecule.addBond(bond.first + j * numAtoms, bond.second + j * numAtoms,
bondOrders.at(i));
}
}

// Now set the unit cell
molecule.unitCell()->setAVector(newA);
molecule.unitCell()->setBVector(newB);
Expand All @@ -595,7 +610,7 @@ struct TransformAtomsFunctor

void operator()(Vector3& pos) { pos = transform * pos; }
};
}
} // namespace

bool CrystalTools::setCellMatrix(Molecule& molecule,
const Matrix3& newCellColMatrix, Options opt)
Expand Down Expand Up @@ -625,7 +640,7 @@ struct FractionalCoordinatesFunctor

void operator()(Vector3& pos) { unitCell.toFractional(pos, pos); }
};
}
} // namespace

bool CrystalTools::fractionalCoordinates(const UnitCell& unitCell,
const Array<Vector3>& cart,
Expand Down Expand Up @@ -664,7 +679,7 @@ struct SetFractionalCoordinatesFunctor

Vector3 operator()(const Vector3& pos) { return unitCell.toCartesian(pos); }
};
}
} // namespace

bool CrystalTools::setFractionalCoordinates(Molecule& molecule,
const Array<Vector3>& coords)
Expand All @@ -684,4 +699,4 @@ bool CrystalTools::setFractionalCoordinates(Molecule& molecule,
return true;
}

} // namespace Avogadro
} // namespace Avogadro::Core
9 changes: 6 additions & 3 deletions avogadro/core/spacegroups.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@ namespace Avogadro::Core {

unsigned short SpaceGroups::hallNumber(const std::string& spaceGroup)
{
unsigned short hall = 0; // can't find anything
unsigned short hall = 0; // can't find anything
const unsigned short hall_count = 531; // 530 but first one is empty
// some files use " instead of = for the space group symbol
std::string sg = spaceGroup;
Expand Down Expand Up @@ -275,6 +275,9 @@ void SpaceGroups::fillUnitCell(Molecule& mol, unsigned short hallNumber,
for (Index j = 1; j < newAtoms.size(); ++j) {
// The new atoms are in fractional coordinates. Convert to cartesian.
Vector3 newCandidate = uc->toCartesian(newAtoms[j]);
// if we are wrapping to the cell, we need to wrap the new atom
if (wrapToCell)
newCandidate = uc->wrapCartesian(newCandidate);

// If there is already an atom in this location within a
// certain tolerance, do not add the atom.
Expand All @@ -298,8 +301,8 @@ void SpaceGroups::fillUnitCell(Molecule& mol, unsigned short hallNumber,
}
}

if (wrapToCell)
CrystalTools::wrapAtomsToUnitCell(mol);
// if (wrapToCell)
// CrystalTools::wrapAtomsToUnitCell(mol);

// Now we need to generate any copies on the unit boundary
// We need to loop through all the atoms again
Expand Down
Loading