Skip to content

Commit

Permalink
Merge pull request #1853 from bbucior/periodic-boundaries
Browse files Browse the repository at this point in the history 
Closes #1536: implement periodic boundary conditions
  • Loading branch information
@ghutchis
ghutchis committed Jan 25, 2020
2 parents e030d2e + 884b164 commit 3d6384b
Show file tree
Hide file tree
Showing 17 changed files with 656 additions and 79 deletions.
2 changes: 2 additions & 0 deletions include/openbabel/atom.h
View file
Expand Up @@ -410,6 +410,8 @@ namespace OpenBabel
bool IsAromaticNOxide();
//! \return Is this atom chiral?
bool IsChiral();
//! \return Is the atom part of a periodic unit cell?
bool IsPeriodic() const;
//! \return Is this atom an axial atom in a ring
bool IsAxial();
//! \return Is this atom a hydrogen-bond acceptor (considering also atom surrounding)
Expand Down
2 changes: 2 additions & 0 deletions include/openbabel/bond.h
View file
Expand Up @@ -217,6 +217,8 @@ namespace OpenBabel
bool IsRotor(bool includeRingBonds=false);
/** \return Is the bond an amide link (i.e., between a carbonyl C and a N)?
No distinction is made between primary, secondary, and tertiary amides. **/
bool IsPeriodic() const;
//! \return Is the bond within a periodic unit cell?
bool IsAmide();
/** \return Is the bond a primary amide (i.e., between carbonyl C and a NH2)?
In versions prior to 2.3, this function incorrectly identified secondary amides. **/
Expand Down
23 changes: 23 additions & 0 deletions include/openbabel/generic.h
View file
Expand Up @@ -487,6 +487,29 @@ namespace OpenBabel
//! \todo Make OBUnitCell::WrapFractionalCoordinate static in the next ABI break
vector3 WrapFractionalCoordinate(vector3 frac);
vector3 WrapFractionalCoordinate(vector3 frac) const;
//! Unwraps Cartesian coordinates near a reference location
//! \param new_loc Cartesian coordinates of target
//! \param ref_loc Cartesian coordinates of the reference location
//! \return Unwrapped coordinates of new_loc near ref_loc
vector3 UnwrapCartesianNear(vector3 new_loc, vector3 ref_loc);
vector3 UnwrapCartesianNear(vector3 new_loc, vector3 ref_loc) const;
//! Unwraps fractional coordinates near a reference location.
//! \param new_loc Fractional coordinates of target
//! \param ref_loc Fractional coordinates of the reference location
//! \return Unwrapped coordinates of new_loc near ref_loc
//! \todo Add a simple test case/example, like unwrapNear(<0.9, 0.2, 0.2>, <0.3, 0.9, 0.2>) -> <-0.1, 1.2, 0.2>
vector3 UnwrapFractionalNear(vector3 new_loc, vector3 ref_loc);
vector3 UnwrapFractionalNear(vector3 new_loc, vector3 ref_loc) const;
//! Applies the minimum image convention to a Cartesian displacement vector
//! \param cart Displacement vector between two atoms in Cartesian coordinates
//! \return Cartesian difference, wrapped within half the unit cell
vector3 MinimumImageCartesian(vector3 cart);
vector3 MinimumImageCartesian(vector3 cart) const;
//! Applies the minimum image convention to a fractional displacement vector
//! \param cart Displacement vector between two atoms in fractional coordinates
//! \return Fractional difference, wrapped within half the unit cell (-0.5 to 0.5)
vector3 MinimumImageFractional(vector3 frac);
vector3 MinimumImageFractional(vector3 frac) const;

//! \return The numeric value of the given spacegroup
int GetSpaceGroupNumber( std::string name = "" );
Expand Down
14 changes: 12 additions & 2 deletions include/openbabel/mol.h
View file
Expand Up @@ -104,7 +104,9 @@ namespace OpenBabel
#define OB_ATOMSPIN_MOL (1<<21)
//! Treat as reaction
#define OB_REACTION_MOL (1<<22)
// flags 22-32 unspecified
//! Molecule is repeating in a periodic unit cell
#define OB_PERIODIC_MOL (1<<23)
// flags 24-32 unspecified

#define SET_OR_UNSET_FLAG(X) \
if (value) SetFlag(X); \
Expand Down Expand Up @@ -387,16 +389,21 @@ enum HydrogenType { AllHydrogen, PolarHydrogen, NonPolarHydrogen };
//! Mark that ring closure bonds have been assigned by graph traversal
void SetClosureBondsPerceived(bool value = true) { SET_OR_UNSET_FLAG(OB_CLOSURE_MOL); }
//! Mark that explicit hydrogen atoms have been added

void SetHydrogensAdded(bool value = true) { SET_OR_UNSET_FLAG(OB_H_ADDED_MOL); }
void SetCorrectedForPH(bool value = true) { SET_OR_UNSET_FLAG(OB_PH_CORRECTED_MOL); }
void SetSpinMultiplicityAssigned(bool value = true) { SET_OR_UNSET_FLAG(OB_ATOMSPIN_MOL); }
//! The OBMol is a pattern, not a complete molecule. Left unchanged by Clear().
void SetIsPatternStructure(bool value = true) { SET_OR_UNSET_FLAG(OB_PATTERN_STRUCTURE); }
void SetIsReaction(bool value = true) { SET_OR_UNSET_FLAG(OB_REACTION_MOL) };
void SetIsReaction(bool value = true) { SET_OR_UNSET_FLAG(OB_REACTION_MOL); }
//! Mark that distance calculations, etc., should apply periodic boundary conditions through the minimimum image convention.
//! Does not automatically recalculate bonding.
void SetPeriodicMol(bool value = true){ SET_OR_UNSET_FLAG(OB_PERIODIC_MOL); }
bool HasFlag(int flag) { return (_flags & flag) ? true : false; }
void SetFlag(int flag) { _flags |= flag; }
void UnsetFlag(int flag) { _flags &= (~(flag)); }
void SetFlags(int flags) { _flags = flags; }

//@}

//! \name Molecule modification methods
Expand Down Expand Up @@ -587,6 +594,9 @@ enum HydrogenType { AllHydrogen, PolarHydrogen, NonPolarHydrogen };
bool HasSpinMultiplicityAssigned() { return(HasFlag(OB_ATOMSPIN_MOL)); }
//! Does this OBMol represent a reaction?
bool IsReaction() { return HasFlag(OB_REACTION_MOL); }
//! Is this molecule periodic? Should periodic boundary conditions be applied?
bool IsPeriodic() { return(HasFlag(OB_PERIODIC_MOL)); }

//! Are there any atoms in this molecule?
bool Empty() { return(_natoms == 0); }
//@}
Expand Down
46 changes: 27 additions & 19 deletions src/atom.cpp
View file
Expand Up @@ -23,6 +23,7 @@ GNU General Public License for more details.
#include <openbabel/bond.h>
#include <openbabel/mol.h>
#include <openbabel/obiter.h>
#include <openbabel/generic.h>
#include <openbabel/molchrg.h>
#include <openbabel/ring.h>
#include <openbabel/phmodel.h>
Expand Down Expand Up @@ -802,6 +803,12 @@ namespace OpenBabel
return stereoFacade.HasTetrahedralStereo(_id);
}

bool OBAtom::IsPeriodic() const
{
OBMol *mol = (OBMol*)((OBAtom*)this)->GetParent();
return mol->IsPeriodic();
}

bool OBAtom::IsInRingSize(int size) const
{
vector<OBRing*> rlist;
Expand Down Expand Up @@ -898,9 +905,7 @@ namespace OpenBabel
k = j;
for (c = NextNbrAtom(k); c; c = NextNbrAtom(k))
{
v1 = b->GetVector() - GetVector();
v2 = c->GetVector() - GetVector();
degrees = vectorAngle(v1, v2);
degrees = b->GetAngle((OBAtom*)this, c);
if (degrees < minDegrees)
minDegrees = degrees;
}
Expand All @@ -924,9 +929,7 @@ namespace OpenBabel
k = j;
for (c = NextNbrAtom(k); c; c = NextNbrAtom(k))
{
v1 = b->GetVector() - GetVector();
v2 = c->GetVector() - GetVector();
degrees = vectorAngle(v1, v2);
degrees = b->GetAngle((OBAtom*)this, c);
avgDegrees += degrees;
n++;
}
Expand Down Expand Up @@ -1102,13 +1105,21 @@ namespace OpenBabel

double OBAtom::GetDistance(OBAtom *b)
{
return(( this->GetVector() - b->GetVector() ).length());
if (!IsPeriodic())
{
return(( this->GetVector() - b->GetVector() ).length());
}
else
{
OBUnitCell *box = (OBUnitCell*)GetParent()->GetData(OBGenericDataType::UnitCell);
return (box->MinimumImageCartesian(this->GetVector() - b->GetVector())).length();
}
}

double OBAtom::GetDistance(int b)
{
OBMol *mol = (OBMol*)GetParent();
return(( this->GetVector() - mol->GetAtom(b)->GetVector() ).length());
return( this->GetDistance(mol->GetAtom(b)) );
}

double OBAtom::GetDistance(vector3 *v)
Expand All @@ -1122,6 +1133,13 @@ namespace OpenBabel

v1 = this->GetVector() - b->GetVector();
v2 = c->GetVector() - b->GetVector();
if (IsPeriodic())
{
OBUnitCell *box = (OBUnitCell*)GetParent()->GetData(OBGenericDataType::UnitCell);
v1 = box->MinimumImageCartesian(v1);
v2 = box->MinimumImageCartesian(v2);
}

if (IsNearZero(v1.length(), 1.0e-3)
|| IsNearZero(v2.length(), 1.0e-3)) {
return(0.0);
Expand All @@ -1133,17 +1151,7 @@ namespace OpenBabel
double OBAtom::GetAngle(int b, int c)
{
OBMol *mol = (OBMol*)GetParent();
vector3 v1,v2;

v1 = this->GetVector() - mol->GetAtom(b)->GetVector();
v2 = mol->GetAtom(c)->GetVector() - mol->GetAtom(b)->GetVector();

if (IsNearZero(v1.length(), 1.0e-3)
|| IsNearZero(v2.length(), 1.0e-3)) {
return(0.0);
}

return(vectorAngle(v1, v2));
return(this->GetAngle(mol->GetAtom(b), mol->GetAtom(c)));
}

bool OBAtom::GetNewBondVector(vector3 &v,double length)
Expand Down
31 changes: 22 additions & 9 deletions src/bond.cpp
View file
Expand Up @@ -24,6 +24,7 @@ GNU General Public License for more details.
#include <openbabel/ring.h>
#include <openbabel/bond.h>
#include <openbabel/mol.h>
#include <openbabel/generic.h>
#include <climits>

using namespace std;
Expand Down Expand Up @@ -96,6 +97,7 @@ namespace OpenBabel
_order = (char)order;
}

// TODO: Figure out how to consider periodicity, etc.
void OBBond::SetLength(OBAtom *fixed, double length)
{
unsigned int i;
Expand Down Expand Up @@ -181,6 +183,12 @@ namespace OpenBabel
return (_bgn->GetHvyDegree() > 1 && _end->GetHvyDegree() > 1);
}

bool OBBond::IsPeriodic() const
{
OBMol *mol = (OBMol*)((OBBond*)this)->GetParent();
return mol->IsPeriodic();
}

bool OBBond::IsAmide()
{
OBAtom *c,*n;
Expand Down Expand Up @@ -491,10 +499,7 @@ namespace OpenBabel
{
if (nbrEnd != _bgn)
{
torsion=fabs(CalcTorsionAngle(nbrStart->GetVector(),
static_cast<OBAtom*>(_bgn)->GetVector(),
static_cast<OBAtom*>(_end)->GetVector(),
nbrEnd->GetVector()));
torsion=fabs(_parent->GetTorsion(nbrStart, _bgn, _end, nbrEnd));

// >12&&<168 not enough
if (torsion > 15.0 && torsion < 160.0)
Expand Down Expand Up @@ -595,11 +600,19 @@ namespace OpenBabel
begin = GetBeginAtom();
end = GetEndAtom();

d2 = SQUARE(begin->GetX() - end->GetX());
d2 += SQUARE(begin->GetY() - end->GetY());
d2 += SQUARE(begin->GetZ() - end->GetZ());

return(sqrt(d2));
if (!IsPeriodic())
{
d2 = SQUARE(begin->GetX() - end->GetX());
d2 += SQUARE(begin->GetY() - end->GetY());
d2 += SQUARE(begin->GetZ() - end->GetZ());
return(sqrt(d2));
}
else
{
OBMol *mol = (OBMol*)((OBBond*)this)->GetParent();
OBUnitCell *box = (OBUnitCell*)mol->GetData(OBGenericDataType::UnitCell);
return (box->MinimumImageCartesian(begin->GetVector() - end->GetVector())).length();
}
}

/*Now in OBBase
Expand Down
1 change: 1 addition & 0 deletions src/formats/cacaoformat.cpp
View file
Expand Up @@ -269,6 +269,7 @@ namespace OpenBabel
a = vit[i]->_a;
b = vit[i]->_b;
c = vit[i]->_c;
// TODO: fix explicit calculations for PBC?
v1 = atom->GetVector() - a->GetVector();
v2 = b->GetVector() - a->GetVector();
vit[i]->_ang = vectorAngle(v1,v2);
Expand Down
97 changes: 96 additions & 1 deletion src/formats/cifformat.cpp
View file
Expand Up @@ -80,7 +80,10 @@ namespace OpenBabel
"Read Options e.g. -ab:\n"
" s Output single bonds only\n"
" b Disable bonding entirely\n"
" B Use bonds listed in CIF file from _geom_bond_etc records (overrides option b) \n\n";
" B Use bonds listed in CIF file from _geom_bond_etc records (overrides option b)\n\n"

"Write Options e.g. -xg:\n"
" g Write bonds using _geom_bond_etc fields \n\n";
};

virtual const char* SpecificationURL()
Expand Down Expand Up @@ -1610,6 +1613,7 @@ namespace OpenBabel
<< " _atom_site_fract_y" << endl
<< " _atom_site_fract_z" << endl
<< " _atom_site_occupancy" << endl;
std::map<OBAtom*,std::string> label_table;
unsigned int i = 0;
FOR_ATOMS_OF_MOL(atom, *pmol)
{
Expand Down Expand Up @@ -1641,13 +1645,104 @@ namespace OpenBabel
label_str = OBElements::GetSymbol(atom->GetAtomicNum()) + to_string(i);
i++;
}
// Save the existing or generated label for optional bonding
label_table[&*atom] = label_str;

snprintf(buffer, BUFF_SIZE, " %-8s%-5s%.5f%10.5f%10.5f%8.3f\n",
label_str.c_str(), OBElements::GetSymbol(atom->GetAtomicNum()),
X, Y, Z, occup);

ofs << buffer;
}

if (pConv->IsOption("g", OBConversion::OUTOPTIONS))
{
if (pmol->NumBonds() > 0) {
obErrorLog.ThrowError(__FUNCTION__, "Writing bonds to output CIF", obDebug);
ofs << "loop_" << endl
<< " _geom_bond_atom_site_label_1" << endl
<< " _geom_bond_atom_site_label_2" << endl
<< " _geom_bond_distance" << endl
<< " _geom_bond_site_symmetry_2" << endl
<< " _ccdc_geom_bond_type" << endl;
} else {
obErrorLog.ThrowError(__FUNCTION__, "No bonds defined in molecule for CIF export", obDebug);
}

FOR_BONDS_OF_MOL(bond, *pmol)
{
std::string label_1 = label_table[bond->GetBeginAtom()];
std::string label_2 = label_table[bond->GetEndAtom()];

std::string sym_key;
int symmetry_num = 555;
if (bond->IsPeriodic()) {
OBUnitCell *box = (OBUnitCell*)pmol->GetData(OBGenericDataType::UnitCell);
vector3 begin, end_orig, end_expected, uc_direction;
// Use consistent coordinates with the X Y Z written in the _atom_site_* loop earlier
begin = box->CartesianToFractional(bond->GetBeginAtom()->GetVector());
begin = box->WrapFractionalCoordinate(begin);
end_orig = box->CartesianToFractional(bond->GetEndAtom()->GetVector());
end_orig = box->WrapFractionalCoordinate(end_orig);
end_expected = box->UnwrapFractionalNear(end_orig, begin);

// To get the signs right, consider the example {0, 0.7}. We want -1 as the periodic direction.
// TODO: Think about edge cases, particularly atoms on the border of the unit cell.
uc_direction = end_expected - end_orig;

std:vector<int> uc;
for (int i = 0; i < 3; ++i) {
double raw_cell = uc_direction[i];
uc.push_back(static_cast<int>(lrint(raw_cell)));
}
symmetry_num += 100*uc[0] + 10*uc[1] + 1*uc[2]; // Unit cell directionality vs. 555, per CIF spec
}
if (symmetry_num == 555)
{
sym_key = ".";
}
else
{
stringstream ss;
ss << "1_" << symmetry_num;
sym_key = ss.str();
}

std::string bond_type;
int bond_order = bond->GetBondOrder();
switch (bond_order)
{
case 1:
bond_type = "S";
break;
case 2:
bond_type = "D";
break;
case 3:
bond_type = "T";
break;
case 5: // FIXME: this will be different in upstream code
bond_type = "A"; // aromatic, per OBBond::_order
break;
default:
stringstream ss;
ss << "Unexpected bond order " << bond_order
<< " for bond" << label_1 << "-" << label_2 << std::endl
<< "Defaulting to single bond.";
obErrorLog.ThrowError(__FUNCTION__, ss.str(), obWarning);
bond_type = "S";
}


//printf("%p: %s\n", &*atom, label_table[&*atom].c_str());
snprintf(buffer, BUFF_SIZE, " %-7s%-7s%10.5f%7s%4s\n",
label_1.c_str(), label_2.c_str(),
bond->GetLength(), sym_key.c_str(),
bond_type.c_str());

ofs << buffer;
}
}
return true;
}//WriteMolecule
} //namespace OpenBabel

0 comments on commit 3d6384b

Please sign in to comment.