1 /*
2 * CifCrystal.cpp
3 *
4 * Copyright Notice: see copyright.txt
5 *
6 * Date: 3/8/2012 (D/M/Y)
7 * Author: Dmitrij Lioubartsev
8 * Email: dmitrijl42@gmail.com
9 *
10 * For class info, see CifCrystal.h
11 */
12
13 #include "CifCrystal.h"
14 #include "Atom.h"
15 #include "SymmetryCalculator.h"
16
17
18 #include <iostream>
19 #include <cmath>
20 using namespace std;
21
22 #define _PI__ 3.14159265
23
CifCrystal()24 CifCrystal::CifCrystal() {
25 //constants
26 a = b = c = 0;
27 alpha = beta = gamma = 0;
28 aNew = bNew = cNew = 0;
29 aDup = bDup = cDup = 0;
30 for(int i = 0; i < 6; i++) {
31 isSet[i] = false;
32 }
33 //symmetry
34 vector<string> symmetries(0);
35 //atoms
36 vector<Atom> atoms(0);
37 numEl = numx = numy = numz = 0;
38 vector<Bond> bonds(0);
39
40 }
41
~CifCrystal()42 CifCrystal::~CifCrystal() {
43 // do nothing
44 }
45
seta(double value)46 bool CifCrystal::seta(double value) {
47 if(!isSet[0]) {
48 a = value;
49 isSet[0] = true;
50 return true;
51 } else {
52 return false;
53 }
54 }
55
setb(double value)56 bool CifCrystal::setb(double value) {
57 if(!isSet[1]) {
58 b = value;
59 isSet[1] = true;
60 return true;
61 } else {
62 return false;
63 }
64 }
65
setc(double value)66 bool CifCrystal::setc(double value) {
67 if(!isSet[2]) {
68 c = value;
69 isSet[2] = true;
70 return true;
71 } else {
72 return false;
73 }
74 }
75
setalpha(double value)76 bool CifCrystal::setalpha(double value) {
77 if(!isSet[3]) {
78 alpha = value;
79 isSet[3] = true;
80 return true;
81 } else {
82 return false;
83 }
84 }
85
setbeta(double value)86 bool CifCrystal::setbeta(double value) {
87 if(!isSet[4]) {
88 beta = value;
89 isSet[4] = true;
90 return true;
91 } else {
92 return false;
93 }
94 }
95
setgamma(double value)96 bool CifCrystal::setgamma(double value) {
97 if(!isSet[5]) {
98 gamma = value;
99 isSet[5] = true;
100 return true;
101 } else {
102 return false;
103 }
104 }
105
106 //returns an array of size 6 to indicate which values has been set.
107 // order is {a, b, c, alpha, beta, gamma}
areSet()108 bool* CifCrystal::areSet() { return isSet;}
109
get_a()110 double CifCrystal::get_a() { return a;}
get_b()111 double CifCrystal::get_b() { return b;}
get_c()112 double CifCrystal::get_c() { return c;}
get_alpha()113 double CifCrystal::get_alpha() { return alpha;}
get_beta()114 double CifCrystal::get_beta() { return beta;}
get_gamma()115 double CifCrystal::get_gamma() { return gamma;}
getNewa()116 double CifCrystal::getNewa() { return aNew;}
getNewb()117 double CifCrystal::getNewb() { return bNew;}
getNewc()118 double CifCrystal::getNewc() { return cNew;}
getDupa()119 int CifCrystal::getDupa() { return aDup;}
getDupb()120 int CifCrystal::getDupb() { return bDup;}
getDupc()121 int CifCrystal::getDupc() { return cDup;}
122
123 /*
124 * Adds a symmetry
125 */
addSymmetry(string symmetry)126 void CifCrystal::addSymmetry(string symmetry) {
127 symmetries.push_back(symmetry);
128 }
129
130 //atoms
131
132 //setters, see CifCrystal.h for more info
setAx(double x)133 void CifCrystal::setAx(double x) {
134 if(numx < atoms.size()) {
135 atoms[numx].setx(x);
136 } else {
137 Atom a = Atom();
138 a.setx(x);
139 atoms.push_back(a);
140 }
141 numx++;
142 }
143
setAy(double y)144 void CifCrystal::setAy(double y) {
145 if(numy < atoms.size()) {
146 atoms[numy].sety(y);
147 } else {
148 Atom a = Atom();
149 a.sety(y);
150 atoms.push_back(a);
151 }
152 numy++;
153 }
154
setAz(double z)155 void CifCrystal::setAz(double z) {
156 if(numz < atoms.size()) {
157 atoms[numz].setz(z);
158 } else {
159 Atom a = Atom();
160 a.setz(z);
161 atoms.push_back(a);
162 }
163 numz++;
164 }
165
setAElement(string e)166 void CifCrystal::setAElement(string e) {
167 if(numEl < atoms.size()) {
168 atoms[numEl].setElement(e);
169 } else {
170 Atom a = Atom();
171 a.setElement(e);
172 atoms.push_back(a);
173 }
174 numEl++;
175 }
176
177 /*
178 * Number of atoms currently in the crystal structure
179 */
numAtoms()180 unsigned int CifCrystal::numAtoms() {
181 return atoms.size();
182 }
183
184 /*
185 * Returns a vector with all atoms
186 */
getAtomList()187 vector<Atom>* CifCrystal::getAtomList() {
188 return &atoms;
189 }
190
191 /*
192 * Returns a vector with atomic bonds.
193 */
getBondList()194 vector<CifCrystal::Bond>* CifCrystal::getBondList() {
195 return &bonds;
196 }
197
198 /*
199 * A check so that all atoms have all properties assigned.
200 */
okData()201 bool CifCrystal::okData() {
202 return (numEl == numx) && (numEl == numy) && (numEl == numz)
203 && (numEl == atoms.size());
204 }
205
206 /*
207 * Apply each symmetry operation on each atom. If the newly created atom
208 * does not already exist (compare coordinates and account for symmetry,
209 * this is what Atom::equals does), add it to current Atom list.
210 * Afterwards, clear symmetry list.
211 * For information about what a symmetry is and how its calculated, see
212 * class SymmetryCalculator.
213 */
processSymmetries()214 void CifCrystal::processSymmetries() {
215
216 SymmetryCalculator sc = SymmetryCalculator();
217 Atom newAtom;
218 vector<Atom> newAtoms(0);
219
220 for(unsigned int i = 0; i < symmetries.size(); i++) {//for each symmetry
221 //string& sym = symmetries[i];
222 if(!sc.setSymmetry(symmetries[i])) { //load the symmetry
223 cout << "WARNING: Symmetry number " << i + 1 << " could not be parsed. "
224 "Ignoring symmetry." << endl;
225 continue;
226 }
227
228
229 for(unsigned int j = 0; j < atoms.size(); j++) { //for each atom
230 sc.loadAtom(atoms[j].getx(),atoms[j].gety(),atoms[j].getz()); //load atom values
231 if(!sc.calculate()) { //and calculate
232 //something went wrong with the calculations
233 cout << "WARNING: Symmetry number " << i + 1 << " could not be calculated."
234 " Ignoring symmetry." << endl;
235 break; //skip the rest of atoms, because an atom can't be "wrong", only the symmetry
236 }
237 //now create atom
238 newAtom = Atom(atoms[j].getElement(),sc.getx(),sc.gety(),sc.getz());
239 newAtom.toCubic();
240 bool isNew = true;
241 //now cheack if atom already exists
242 for(unsigned int k = 0; k < atoms.size();k++) { //check current atoms
243 if(newAtom.equals(atoms[k]) ) {
244 isNew = false;
245 break;
246 }
247 }
248 if(isNew)
249 for(unsigned int k = 0; k < newAtoms.size();k++) { //check newly added atoms
250 if(newAtom.equals(newAtoms[k]) ) {
251 isNew = false;
252 break;
253 }
254 }
255
256 if(isNew) {
257 newAtoms.push_back(newAtom);
258 }
259 }
260 }
261
262 //for all new atoms, add to the atom list.
263 for(unsigned int k = 0; k < newAtoms.size(); k++) {
264 atoms.push_back(newAtoms[k]);
265 }
266
267 //clear all the processed symmetries.
268 symmetries.clear();
269
270 }
271
272 /*
273 * Creates a supercell. This means, duplicate atom coordinates
274 * in x,y and z directions respectivly, so that the new cell
275 * (called supercell), has sides a,b,c that are at minimum
276 * r_min.
277 * E.g. r_min = 20, a = 9. Duplicate into x-direction 3 times,
278 * so that new a becomes 27 > 20.
279 */
makeSuperCell(double r_min)280 void CifCrystal::makeSuperCell(double r_min) {
281
282 int i = 1; //number of copies.
283 vector<Atom> newAtoms; //a holder vector for new atoms.
284 while((i*a) < r_min) { //duplicate in x-direction
285 for(unsigned int k = 0;k < atoms.size();k++) {
286 Atom& at = atoms[k]; //for easier-to-read code
287 //now add a new atom, with the same properties as the current one ("at"),
288 //but with its x-coordinate incremented with i.
289 newAtoms.push_back(Atom(at.getElement(),at.getx()+i,at.gety(),at.getz()));
290 }
291 i++;
292 }
293 /*for(unsigned int j = 0;j < newAtoms.size();j++) {
294 atoms.push_back(newAtoms[j]);
295 }*/ //this is the alternative way to do this
296 for(vector<Atom>::iterator it = newAtoms.begin(); it != newAtoms.end();it++) {
297 atoms.push_back(*it); //add all atoms from newAtoms to atoms.
298 }
299 aNew = a*i; //the new a.
300 aDup = i;
301 newAtoms.clear(); //clear
302 i=1; //reset i for next iteration
303 while((i*b) < r_min) { //duplicate in y-direction
304 for(unsigned int k = 0;k < atoms.size();k++) {
305 Atom& at = atoms[k]; //for easier-to-read code
306 newAtoms.push_back(Atom(at.getElement(),at.getx(),at.gety()+i,at.getz()));
307 }
308 i++;
309 }
310 for(vector<Atom>::iterator it = newAtoms.begin(); it != newAtoms.end();it++) {
311 atoms.push_back(*it);
312 }
313 bNew = b*i; //the new b.
314 bDup = i;
315 newAtoms.clear();
316 i=1;
317 while((i*c) < r_min) { //duplicate in z-direction
318 for(unsigned int k = 0;k < atoms.size();k++) {
319 Atom& at = atoms[k]; //for easier-to-read code
320 newAtoms.push_back(Atom(at.getElement(),at.getx(),at.gety(),at.getz()+i));
321 }
322 i++;
323 }
324 for(vector<Atom>::iterator it = newAtoms.begin(); it != newAtoms.end();it++) {
325 atoms.push_back(*it);
326 }
327 cNew = c*i; //the new c.
328 cDup = i;
329 }
330
makeTransformations()331 void CifCrystal::makeTransformations() {
332 /*
333 * First, a quick lesson in linear algebra and how this is done, and how
334 * stuff was calculated, because having it only written down on a paper
335 * that I later may or may not find can be bad, in case something will
336 * need to be altered and I won't understand anything here.
337 */
338
339 /*
340 * The atoms in the file are given in relative coordinates, in the
341 * coordinate system of a,b,c, which are vectors in the euclidian system
342 * with lengths a,b,c, and with angles alpha (between b,c), beta (between a,c)
343 * and gamma (between a,b) between them. To not think of the relative thingy
344 * on the coordinates, I see the current coordinates as being written in
345 * the basis of {x',y',z'}, where |x'| == a, |y'| == b, |z'| == c, and with
346 * the given angles. (basically the basis are the vectors a,b,c but technically
347 * a,b,c represents scalars (lengths) and not vectors, so I can't really write
348 * that).
349 *
350 * If I set the x-axis to be oriented along the vector of a, then x' will be
351 * (a,0,0). That is one of the basis vectors for the current coordinate system
352 * (or vector room). Now to y'. See triangle below, which is drawn in x-y plane:
353 * /|
354 * / |
355 * /__|
356 *
357 * Bot-left corner is angle gamma, and the hypotenuse length is b. Lets call the
358 * bot side f and right side h. Then we will have b*cos(gamma) = f, b*sin(gamma) = h.
359 * This will give y' = (b*cos(gamma),b*sin(gamma),0).
360 *
361 * With similar geometry, z'.x will be c*cos(beta), and z'.y == c*cos(alpha). It isn't
362 * that easy to calculate z'.z, so I will use the formula (z'.x)²+(z'.y)²+(z'.z)² == c².
363 * Solving for z'.z, I get, in the end, that
364 * z' = (c*cos(beta),c*cos(alpha),c*sqrt(1-cos²(beta)-cos²(alpha))).
365 *
366 * So, in the .cif file the coordinates are given in the basis of
367 * (a,0,0) , (b*cos(gamma), b*sin(gamma),0), (c*cos(beta),c*cos(alpha),c*sqrt(1-cos²(beta)-cos²(alpha))).
368 *
369 * I want to transform it to euclidian basis (1,0,0),(0,1,0),(0,0,1). To do that, I use
370 * the formula B = T*B'. B' is the vector given in the "old" basis, and B are the coordinates
371 * in the new euclidian basis. T is the transformation matrix where its columns are given by
372 * x',y' and z'. (not gonna write it out, too large).
373 *
374 * Esentially, using simple matrix multiplication, this formulas are given, where Ax,Ay,Az are the atom
375 * coordinates as written in the file:
376 *
377 * Nx = Ax*a + Ay*b*cos(gamma) + Ax*c*cos(beta)
378 * Ny = Ay*b*sin(gamma) + Az*c*cos(alpha)
379 * Nz = Az*c*sqrt(1-cos²(beta)-cos²(alpha))
380 *
381 *IMPORTANT IMPORTANT IMPORTANT
382 * NOTE: Transformation matrix is slightly different, the z'.y and z'.z are not the same.
383 * See http://en.wikipedia.org/wiki/Fractional_coordinates for more info. I got a bit
384 * overzealous with these. Using the wiki Transformation matrix.
385 *
386 * NOTE: AFTER more testing, it appears my own transformation matrix gives exact same
387 * result as wikipedia's one... Using my own one just because. And they are simpler.
388 *
389 * Apply these on every atom.
390 */
391
392 double alphar, betar, gammar; //angles in radian.
393
394 //cmath's trigonometric functions need radian angles.
395 alphar = alpha*_PI__/180;
396 betar = beta*_PI__/180;
397 gammar = gamma*_PI__/180;
398
399 //these are some stuff that needs to be calculated many times, might as well
400 //do it once here and keep the result.
401 const double cosg = cos(gammar);
402 const double cosb = cos(betar);
403 const double sing = sin(gammar);
404 const double cosa = cos(alphar);
405
406 const double zeta = sqrt(1-pow(cosb,2)-pow(cosa,2));
407 //const double u = ((cos(alphar)-(cosb*cosg))/(sing));
408 //const double v = sqrt(1-pow(cosb,2)-pow(cos(alphar),2)-pow(cosg,2) + (2*cos(alphar)*cosb*cosg));
409
410 //now make the cell cubic, instead of a parallelepiped, so new values
411 //of aNew, bNew and cNew are needed, for the makeSmall() function.
412 //aNew = aNew; stays the same
413 //bNew = bNew * sing; //the y part of the b vector.
414 //cNew = cNew * zeta; //the z-part of the c-vector
415 //cNew = cNew * v/sing; //the wiki way of doing it.
416
417
418 for(unsigned int i = 0; i < atoms.size(); i++) {
419 const double& Ax = atoms[i].getx();
420 const double& Ay = atoms[i].gety();
421 const double& Az = atoms[i].getz();
422 double newx, newy, newz;
423 newx = (Ax*a + Ay*b*cosg + Az*c*cosb);
424 newy = (Ay*b*sing + Az*c*cosa);
425 newz = (Az*c*zeta);
426 //newy = (Ay*b*sing + Az*c*u);
427 //newz = (Az*c*v/sing);
428 atoms[i].setx(newx);
429 atoms[i].sety(newy);
430 atoms[i].setz(newz);
431
432 }
433 }
434
435 /*
436 * Create atomic bonds. An atomic bond should be created between two atoms
437 * that have distance less than 2 Å between each other. However, since the
438 * cell is repeated many times, atoms may be bound to other atoms outside
439 * of the cell. To account for that, bind to the corresponding atom in
440 * the current cell.
441 */
makeBonds()442 void CifCrystal::makeBonds() {
443 //first initialize some variables
444
445 double alphar, betar, gammar;
446
447 //cmath's trigonometric functions need radian angles.
448 alphar = alpha*_PI__/180;
449 betar = beta*_PI__/180;
450 gammar = gamma*_PI__/180;
451
452 double bx = bNew*cos(gammar);
453 double by = bNew*sin(gammar);
454 double cx = cNew*cos(betar);
455 double cy = cNew*cos(alphar);
456 double cz = cNew*sqrt(1-pow(cos(betar),2)-pow(cos(alphar),2));
457
458
459
460 for(unsigned int mainInt = 0; mainInt < atoms.size(); mainInt++) { //for each atom
461 /*
462 * For each atom, compute the distance to every other atom. If it is
463 * less than 2, create a bond.
464 */
465 Atom& a1 = atoms[mainInt];
466 for(unsigned int otherInt = 0; otherInt < atoms.size(); otherInt++) {
467 if(mainInt!=otherInt) { //no point comparing atom to itself.
468 Atom& a2 = atoms[otherInt];
469 //now get distance between a1 and a2
470 double minDistance = 100;
471 for(int i = -1; i < 2; i++) {//from -1 to 1
472 for(int j = -1; j < 2; j++) {
473 for(int k = -1; k < 2; k++) {
474 double tmp = sqrt(
475 pow(a1.getx()-a2.getx()+i*aNew+j*bx+k*cx,2) +
476 pow(a1.gety()-a2.gety()+j*by+k*cy,2) +
477 pow(a1.getz()-a2.getz()+k*cz,2));
478 if(tmp < minDistance)
479 minDistance = tmp;
480 }
481 }
482 }
483 if(minDistance <= 2) {
484 //add bond to atom
485 a1.addBond(&a2);
486 if(mainInt < otherInt) {//dont add duplicate bonds.
487 Bond b = {mainInt,otherInt};
488 bonds.push_back(b);
489 }
490 }
491 }
492 }
493 }
494 }
495
496 //for testing purposes
printThis()497 void CifCrystal::printThis() {
498 cout.setf(ios_base::fixed);
499 cout.precision(4);
500 cout << "----- CELL VALUES -----" << endl
501 << "a: ";
502 if(isSet[0])
503 cout << a << endl;
504 else
505 cout << "NOT SET" << endl;
506 cout << "b: ";
507 if(isSet[1])
508 cout << b << endl;
509 else
510 cout << "NOT SET" << endl;
511 cout << "c: ";
512 if(isSet[2])
513 cout << c << endl;
514 else
515 cout << "NOT SET" << endl;
516 cout << "alpha: ";
517 if(isSet[3])
518 cout << alpha << endl;
519 else
520 cout << "NOT SET" << endl;
521 cout << "beta: ";
522 if(isSet[4])
523 cout << beta << endl;
524 else
525 cout << "NOT SET" << endl;
526 cout << "gamma: ";
527 if(isSet[5])
528 cout << gamma << endl;
529 else
530 cout << "NOT SET" << endl;
531 cout << endl;
532 cout << "----- SYMMETRIES -----" << endl;
533 for(unsigned int i=0;i<symmetries.size();i++) {
534 cout << symmetries[i] << endl;
535 }
536 cout << endl;
537 int numO = 0;
538 int numSi = 0;
539 cout << "----- ATOMS -----" << endl;
540 if(!okData()) {
541 cout << "BLÄH: " << numEl << ' '
542 << numx << ' ' << numy << ' '
543 << numz << endl;
544 }
545
546 for(unsigned int i=0;i<atoms.size();i++) {
547 cout << atoms[i].getElement() << i
548 << "\t\t" << atoms[i].getx()
549 << "\t\t" << atoms[i].gety()
550 << "\t\t" << atoms[i].getz() << endl;
551
552 if(atoms[i].getElement() == "O")
553 numO++;
554 else if(atoms[i].getElement() == "Si")
555 numSi++;
556
557 }
558 cout << "NumO: " << numO << "\tnumSi: " << numSi;
559
560 cout << endl << endl;
561 unsigned int i;
562 for(i = 0; i < bonds.size();i++) {
563 cout << bonds[i].a1 << ' ' << bonds[i].a2 << endl;
564 }
565 cout << "NumBounds: " << i << endl;
566
567 }
568
569
570