libdiffpy/src/diffpy/srreal/ObjCrystStructureAdapter.cpp at develop · cfarrow/libdiffpy · GitHub

Repository navigation

History

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

239

240

241

242

243

244

245

246

247

248

249

250

251

252

253

254

255

256

257

258

259

260

261

262

263

264

/*****************************************************************************

* diffpy.srreal by DANSE Diffraction group

* Simon J. L. Billinge

* File coded by: Chris Farrow

* See AUTHORS.txt for a list of people who contributed.

* See LICENSE.txt for license information.

******************************************************************************

* class ObjCrystStructureAdapter

* -- adapter to the Crystal class from ObjCryst++.

* class ObjCrystBondGenerator

* -- Generate bonds from periodic ObjCrystStructureAdapter.

* Factory function createStructureAdapter(ObjCryst::Molecule)

* -- builds AtomicStructureAdapter from the ObjCryst Molecule object

*****************************************************************************/

#include <diffpy/PythonInterface.hpp>

#include <algorithm>

#include <cassert>

#include <cmath>

#include <set>

#include <diffpy/srreal/PythonStructureAdapter.hpp>

#include <diffpy/srreal/ObjCrystStructureAdapter.hpp>

using namespace std;

namespace diffpy {

namespace srreal {

namespace {

// Utility functions ---------------------------------------------------------

R3::Matrix

getUij(const ObjCryst::ScatteringPower* sp)

{

const double BtoU = 1.0 / (8 * M_PI * M_PI);

R3::Matrix Uij = R3::zeromatrix();

if (sp->IsIsotropic())

{

Uij(0,0) = Uij(1,1) = Uij(2,2) = sp->GetBiso() * BtoU;

}

else

{

Uij(0,0) = sp->GetBij(1,1) * BtoU;

Uij(1,1) = sp->GetBij(2,2) * BtoU;

Uij(2,2) = sp->GetBij(3,3) * BtoU;

Uij(0,1) = Uij(1,0) = sp->GetBij(1,2) * BtoU;

Uij(0,2) = Uij(2,0) = sp->GetBij(1,3) * BtoU;

Uij(1,2) = Uij(2,1) = sp->GetBij(2,3) * BtoU;

}

return Uij;

}

CrystalStructureAdapter::SymOpVector

fetchSymmetryOperations(const ObjCryst::SpaceGroup& spacegroup)

{

CrystalStructureAdapter::SymOpVector rv;

const int nbsym = spacegroup.GetNbSymmetrics();

rv.resize(nbsym);

const vector<ObjCryst::SpaceGroup::SMx>& sgsymops =

spacegroup.GetSymmetryOperations();

// copy convert the ObjCryst symmetry operations

int last = sgsymops.size();

assert(last <= nbsym);

for (int i = 0; i != last; ++i)

{

CrystalStructureAdapter::SymOpVector::value_type& op = rv[i];

copy(sgsymops[i].mx, sgsymops[i].mx + 9, op.R.data().begin());

copy(sgsymops[i].tr, sgsymops[i].tr + 3, op.t.data().begin());

}

// apply all origin translations

const vector<ObjCryst::SpaceGroup::TRx>& sgtrans =

spacegroup.GetTranslationVectors();

const int nbtran = spacegroup.GetNbTranslationVectors();

assert(int(sgtrans.size()) == nbtran);

assert(nbtran * last <= nbsym);

for (int nt = 0; nt < nbtran; ++nt)

{

const double* pt = sgtrans[nt].tr;

R3::Vector sgt(pt[0], pt[1], pt[2]);

for (int i = 0; i < last; ++i)

{

int k = i + nt * last;

assert(k < nbsym);

rv[k] = rv[i];

rv[k].t += sgt;

}

last = sgsymops.size() * nbtran;

// apply the inversion center if necessary

if (last < nbsym)

{

assert(spacegroup.IsCentrosymmetric());

assert(nbsym == 2 * last);

bool noCenter = false, noTransl = true, noIdentical = true;

CrystMatrix_REAL xyzinv = spacegroup.GetAllSymmetrics(

0.0, 0.0, 0.0, noCenter, noTransl, noIdentical);

// check if symmetry center is away from the origin

diffpy::mathutils::EpsilonEqual allclose;

R3::Vector dxyz = R3::zerovector;

for (int i = 0; i < xyzinv.rows(); ++i)

{

using mathutils::eps_eq;

R3::Vector xyzi(

round(xyzinv(i, 0) / 0.5),

round(xyzinv(i, 1) / 0.5),

round(xyzinv(i, 2) / 0.5));

if (eps_eq(0.0, R3::norm(xyzi))) continue;

// here we found inversion center at non-origin location

dxyz = xyzi - R3::floor(xyzi);

break;

}

for (int i = 0, j = last; i < last; ++i, ++j)

{

rv[j] = rv[i];

rv[j].R *= -1;

rv[j].t = dxyz - rv[j].t;

}

last *= 2;

}

assert(nbsym == last);

return rv;

}

} // namespace

//////////////////////////////////////////////////////////////////////////////

// StructureAdapter factory for ObjCryst::Crystal

//////////////////////////////////////////////////////////////////////////////

StructureAdapterPtr

createStructureAdapter(const ObjCryst::Crystal& cryst)

{

const double radtodeg = 180 / M_PI;

CrystalStructureAdapterPtr adpt(new CrystalStructureAdapter);

adpt->setLatPar(

cryst.GetLatticePar(0),

cryst.GetLatticePar(1),

cryst.GetLatticePar(2),

radtodeg * cryst.GetLatticePar(3),

radtodeg * cryst.GetLatticePar(4),

radtodeg * cryst.GetLatticePar(5));

// find out number of scatterers in the asymmetric unit

const ObjCryst::ScatteringComponentList& scl =

cryst.GetScatteringComponentList();

size_t nbComponent = scl.GetNbComponent();

adpt->reserve(nbComponent);

Atom ai;

const Lattice& L = adpt->getLattice();

for (size_t i = 0; i < nbComponent; ++i)

{

const ObjCryst::ScatteringComponent& sc = scl(i);

const ObjCryst::ScatteringPower* sp = sc.mpScattPow;

// Skip over this if it is a dummy atom. A dummy atom has no

// mpScattPow, and therefore no type. It's just in a structure as a

// reference position.

if (sp == NULL) continue;

ai.occupancy = sc.mOccupancy;

ai.anisotropy = !(sp->IsIsotropic());

ai.atomtype = sp->GetSymbol();

R3::Vector xyz(sc.mX, sc.mY, sc.mZ);

ai.cartesianposition = L.cartesian(xyz);

// Store Uij

R3::Matrix uijl = getUij(sp);

ai.cartesianuij = ai.anisotropy ? L.cartesianMatrix(uijl) : uijl;

adpt->append(ai);

}

const ObjCryst::SpaceGroup& spacegroup = cryst.GetSpaceGroup();

CrystalStructureAdapter::SymOpVector symops =

fetchSymmetryOperations(spacegroup);

CrystalStructureAdapter::SymOpVector::const_iterator op;

for (op = symops.begin(); op != symops.end(); ++op) adpt->addSymOp(*op);

adpt->updateSymmetryPositions();

return adpt;

}

//////////////////////////////////////////////////////////////////////////////

// StructureAdapter factory for ObjCryst::Molecule

//////////////////////////////////////////////////////////////////////////////

StructureAdapterPtr

createStructureAdapter(const ObjCryst::Molecule& molecule)

{

AtomicStructureAdapterPtr adpt(new AtomicStructureAdapter);

size_t nbComponent = molecule.GetNbComponent();

adpt->reserve(nbComponent);

Atom ai;

for (size_t i = 0; i < nbComponent; ++i)

{

const ObjCryst::MolAtom& atom = molecule.GetAtom(i);

if (atom.IsDummy()) continue;

ai.occupancy = atom.GetOccupancy();

const ObjCryst::ScatteringPower* sp = &(atom.GetScatteringPower());

ai.anisotropy = !(sp->IsIsotropic());

ai.atomtype = sp->GetSymbol();

ai.cartesianposition[0] = atom.X();

ai.cartesianposition[1] = atom.Y();

ai.cartesianposition[2] = atom.Z();

// Store Uij

ai.cartesianuij = getUij(sp);

adpt->append(ai);

}

return adpt;

}

// Factory Function and its Registration -------------------------------------

StructureAdapterPtr

createPyObjCrystStructureAdapter(boost::python::object stru)

{

using diffpy::importFromPyModule;

boost::python::object cls_Crystal, None;

cls_Crystal = importFromPyModule("pyobjcryst.crystal", "Crystal", None);

StructureAdapterPtr rv;

if (cls_Crystal.ptr() != Py_None &&

PyObject_IsInstance(stru.ptr(), cls_Crystal.ptr()) == 1)

{

ObjCryst::Crystal* pcryst =

boost::python::extract<ObjCryst::Crystal*>(stru);

rv = createStructureAdapter(*pcryst);

}

return rv;

}

StructureAdapterPtr

createPyObjCrystMoleculeAdapter(boost::python::object stru)

{

using diffpy::importFromPyModule;

boost::python::object cls_Molecule, None;

cls_Molecule = importFromPyModule("pyobjcryst.molecule", "Molecule", None);

StructureAdapterPtr rv;

if (cls_Molecule.ptr() != Py_None &&

PyObject_IsInstance(stru.ptr(), cls_Molecule.ptr()) == 1)

{

ObjCryst::Molecule* pmolecule =

boost::python::extract<ObjCryst::Molecule*>(stru);

rv = createStructureAdapter(*pmolecule);

}

return rv;

}

bool reg_PyObjCrystStructureAdapter =

registerPythonStructureAdapterFactory(createPyObjCrystStructureAdapter);

bool reg_PyObjCrystMoleculeAdapter =

registerPythonStructureAdapterFactory(createPyObjCrystMoleculeAdapter);

} // namespace srreal

} // namespace diffpy

// End of file

libdiffpy/src/diffpy/srreal/ObjCrystStructureAdapter.cpp at develop · cfarrow/libdiffpy · GitHub

Repository navigation

Expand file tree

Breadcrumbs

ObjCrystStructureAdapter.cpp

Latest commit

History

Breadcrumbs

ObjCrystStructureAdapter.cpp

File metadata and controls