PSCF v1.2
UnitCell3.cpp
1/*
2* PSCF - Polymer Self-Consistent Field Theory
3*
4* Copyright 2016 - 2022, The Regents of the University of Minnesota
5* Distributed under the terms of the GNU General Public License.
6*/
7
8#include "UnitCell.h"
9#include <util/math/Constants.h>
10
11namespace Pscf {
12namespace Prdc {
13
14 using namespace Util;
15
16 /*
17 * Constructor.
18 */
19 UnitCell<3>::UnitCell()
20 : lattice_(Null)
21 {}
22
23 /*
24 * Read the lattice system and set nParameter.
25 */
26 void UnitCell<3>::setNParameter()
27 {
28 UTIL_CHECK(lattice_ != UnitCell<3>::Null);
29 if (lattice_ == UnitCell<3>::Cubic) {
30 nParameter_ = 1;
31 } else
32 if (lattice_ == UnitCell<3>::Tetragonal) {
33 nParameter_ = 2;
34 } else
35 if (lattice_ == UnitCell<3>::Orthorhombic) {
36 nParameter_ = 3;
37 } else
38 if (lattice_ == UnitCell<3>::Monoclinic) {
39 nParameter_ = 4;
40 } else
41 if (lattice_ == UnitCell<3>::Triclinic) {
42 nParameter_ = 6;
43 } else
44 if (lattice_ == UnitCell<3>::Rhombohedral) {
45 nParameter_ = 2;
46 } else
47 if (lattice_ == UnitCell<3>::Hexagonal) {
48 nParameter_ = 2;
49 } else {
50 UTIL_THROW("Invalid value");
51 }
52 }
53
54 /*
55 * Set Bravais and reciprocal lattice parameters, and drBasis.
56 *
57 * Function UnitCellBase::initializeToZero() must be called before
58 * setBasis() to initialize all elements of vectors and tensors to
59 * zero. Only nonzero elements need to be set here.
60 */
61 void UnitCell<3>::setBasis()
62 {
63 UTIL_CHECK(lattice_ != UnitCell<3>::Null);
64 UTIL_CHECK(nParameter_ > 0);
65
66 // Set elements for specific lattice types
67 double twoPi = 2.0*Constants::Pi;
68 int i;
69 if (lattice_ == UnitCell<3>::Cubic) {
70 UTIL_CHECK(nParameter_ == 1);
71 for (i=0; i < 3; ++i) {
72 rBasis_[i][i] = parameters_[0];
73 kBasis_[i][i] = twoPi/parameters_[0];
74 drBasis_[0](i,i) = 1.0;
75 }
76 } else
77 if (lattice_ == UnitCell<3>::Tetragonal) {
78 UTIL_CHECK(nParameter_ == 2);
79
80 rBasis_[0][0] = parameters_[0];
81 rBasis_[1][1] = parameters_[0];
82 rBasis_[2][2] = parameters_[1];
83
84 drBasis_[0](0,0) = 1.0;
85 drBasis_[0](1,1) = 1.0;
86 drBasis_[1](2,2) = 1.0;
87
88 kBasis_[0][0] = twoPi/parameters_[0];
89 kBasis_[1][1] = kBasis_[0][0];
90 kBasis_[2][2] = twoPi/parameters_[1];
91
92 }else
93 if (lattice_ == UnitCell<3>::Orthorhombic) {
94 UTIL_CHECK(nParameter_ == 3);
95 for (i=0; i < 3; ++i) {
96 rBasis_[i][i] = parameters_[i];
97 drBasis_[i](i,i) = 1.0;
98 kBasis_[i][i] = twoPi/parameters_[i];
99 }
100 } else
101 if (lattice_ == UnitCell<3>::Hexagonal) {
102 UTIL_CHECK(nParameter_ == 2);
103 double a = parameters_[0];
104 double c = parameters_[1];
105 double rt3 = sqrt(3.0);
106
107 rBasis_[0][0] = a;
108 rBasis_[1][0] = -0.5*a;
109 rBasis_[1][1] = 0.5*rt3*a;
110 rBasis_[2][2] = c;
111
112 drBasis_[0](0, 0) = 1.0;
113 drBasis_[0](1, 0) = -0.5;
114 drBasis_[0](1, 1) = 0.5*rt3;
115 drBasis_[1](2, 2) = 1.0;
116
117 kBasis_[0][0] = twoPi/a;
118 kBasis_[0][1] = twoPi/(rt3*a);
119 kBasis_[1][0] = 0.0;
120 kBasis_[1][1] = twoPi/(0.5*rt3*a);
121 kBasis_[2][2] = twoPi/(c);
122
123 } else
124 if (lattice_ == UnitCell<3>::Rhombohedral) {
125 UTIL_CHECK(nParameter_ == 2);
126
127 // Set parameters
128 double a, beta, theta;
129 a = parameters_[0]; // magnitude of Bravais basis vectors
130 beta = parameters_[1]; // angle between basis vectors
131 theta = acos( sqrt( (2.0*cos(beta) + 1.0)/3.0 ) );
132 // theta is the angle of all basis vectors from the z axis
133
134 /*
135 * Bravais lattice vectora a_i with i=0, 1, or 2 have form:
136 *
137 * a_i = a * (z * cos(theta) + u_i * sin(theta) )
138 *
139 * where z denotes a unit vector parallel to the z axis, and
140 * the u_i are three unit vectors in the x-y separated by
141 * angles of 2 pi /3 (or 120 degrees), for which u_0 is
142 * parallel to the x axis. Note that u_i . u_j = -1/2 for any
143 * i not equal to j.
144 *
145 * The angle between any two such Bravais basis vectors is easily
146 * computed from the dot product, a_i . a_j for unequal i and j.
147 * This gives:
148 *
149 * cos beta = cos^2(theta) - 0.5 sin^2 (theta)
150 * = (3cos^{2}(theta) - 1)/2
151 *
152 * The angle theta is computed above by solving this equation
153 * for theta as a function of the input parameter beta.
154 *
155 * The corresponding reciprocal lattice vectors have the form
156 *
157 * b_i = (2*pi/3*a)( z/cos(theta) + 2 * u_i / sin(theta) )
158 *
159 * It is straightforward to confirm that these reciprocal
160 * vectors satisfy the condition a_i . b_j = 2 pi delta_{ij}
161 */
162
163 // Trig function values
164 double cosT = cos(theta);
165 double sinT = sin(theta);
166 double sqrt3 = sqrt(3.0);
167 double sinPi3 = 0.5*sqrt3;
168 double cosPi3 = 0.5;
169
170 // Bravais basis vectors
171 rBasis_[0][0] = sinT * a;
172 rBasis_[0][2] = cosT * a;
173 rBasis_[1][0] = -cosPi3 * sinT * a;
174 rBasis_[1][1] = sinPi3 * sinT * a;
175 rBasis_[1][2] = cosT * a;
176 rBasis_[2][0] = -cosPi3 * sinT * a;
177 rBasis_[2][1] = -sinPi3 * sinT * a;
178 rBasis_[2][2] = cosT * a;
179
180 // Reciprocal lattice basis vectors
181 kBasis_[0][0] = 2.0 * twoPi /( 3.0 * sinT * a);
182 kBasis_[0][2] = twoPi /( 3.0 * cosT * a);
183 kBasis_[1][0] = -2.0 * cosPi3 * twoPi /( 3.0 * sinT * a);
184 kBasis_[1][1] = 2.0 * sinPi3 * twoPi /( 3.0 * sinT * a);
185 kBasis_[1][2] = twoPi /( 3.0 * cosT * a);
186 kBasis_[2][0] = -2.0 * cosPi3 * twoPi /( 3.0 * sinT * a);
187 kBasis_[2][1] = -2.0 * sinPi3 * twoPi /( 3.0 * sinT * a);
188 kBasis_[2][2] = twoPi /( 3.0 * cosT * a);
189
190 // Derivatives with respect to length a
191 drBasis_[0](0,0) = sinT;
192 drBasis_[0](0,2) = cosT;
193 drBasis_[0](1,0) = -cosPi3 * sinT;
194 drBasis_[0](1,1) = sinPi3 * sinT;
195 drBasis_[0](1,2) = cosT;
196 drBasis_[0](2,0) = -cosPi3 * sinT;
197 drBasis_[0](2,1) = -sinPi3* sinT;
198 drBasis_[0](2,2) = cosT;
199
200 // Define alpha = d(theta)/d(beta)
201 double alpha = -sin(beta)/(3.0*cosT*sinT);
202
203 // Derivatives with respect to beta
204 drBasis_[1](0,0) = cosT * alpha * a;
205 drBasis_[1](0,2) = -sinT * alpha * a;
206 drBasis_[1](1,0) = -cosPi3 * cosT * alpha * a;
207 drBasis_[1](1,1) = sinPi3 * cosT * alpha * a;
208 drBasis_[1](1,2) = -sinT * alpha * a;
209 drBasis_[1](2,0) = -cosPi3 * cosT * alpha * a;
210 drBasis_[1](2,1) = -sinPi3* cosT * alpha * a;
211 drBasis_[1](2,2) = -sinT * alpha * a;
212
213 } else
214 if (lattice_ == UnitCell<3>::Monoclinic) {
215 UTIL_CHECK(nParameter_ == 4);
216
217 /*
218 * Description: Bravais basis vectors A, B, and C
219 * A (vector 0) is along x axis
220 * B (vector 1) is along y axis
221 * C (vector 2) is in the x-z plane, an angle beta from A/x
222 * B is the unique axis (perpendicular to A and C)
223 * Angle beta is stored in radians.
224 */
225
226 // Parameters
227 double a = parameters_[0]; // length of A
228 double b = parameters_[1]; // length of B
229 double c = parameters_[2]; // length of C
230 double beta = parameters_[3]; // angle between C and A/x
231
232 double cb = cos(beta);
233 double sb = sin(beta);
234
235 // Nonzero components of basis vectors
236 // For rBasis_[i][j], i:basis vector, j:Cartesian component
237 rBasis_[0][0] = a;
238 rBasis_[1][1] = b;
239 rBasis_[2][0] = c*cb;
240 rBasis_[2][2] = c*sb;
241
242 // Nonzero derivatives of basis vectors with respect to parameters
243 // For drBasis_[k](i,j), k:parameter, i:basis vector, j:Cartesian
244 drBasis_[0](0,0) = 1.0;
245 drBasis_[1](1,1) = 1.0;
246 drBasis_[2](2,0) = cb;
247 drBasis_[2](2,2) = sb;
248 drBasis_[3](2,0) = c*sb;
249 drBasis_[3](2,2) = -c*cb;
250
251 // Reciprocal lattice vectors
252 kBasis_[0][0] = twoPi/a;
253 kBasis_[0][2] = -twoPi*cb/(a*sb);
254 kBasis_[1][1] = twoPi / b;
255 kBasis_[2][2] = twoPi/(c*sb);
256
257 } else
258 if (lattice_ == UnitCell<3>::Triclinic) {
259 UTIL_CHECK(nParameter_ == 6);
260
261 /*
262 * Description: Bravais basis vectors A, B, and C
263 * A (vector 0) is along x axis
264 * B (vector 1) is in the x-y plane, an angle gamma from x axis
265 * C (vector 2) is tilted by an angle theta from z axis
266 * phi is the angle beween c-z an x-z (or a-z) planes
267 * All three angles are stored in radians
268 */
269 double a = parameters_[0]; // length of A
270 double b = parameters_[1]; // length of B
271 double c = parameters_[2]; // length of C
272 double phi = parameters_[3]; // angle between c-z and a-z planes
273 double theta = parameters_[4]; // angle between c and z axis
274 double gamma = parameters_[5]; // angle between a and b
275
276 // sine and cosine of all angles
277 double cosPhi = cos(phi);
278 double sinPhi = sin(phi);
279 double cosTheta = cos(theta);
280 double sinTheta = sin(theta);
281 double cosGamma = cos(gamma);
282 double sinGamma = sin(gamma);
283
284 // Nonzero components of Bravais basis vectors
285 rBasis_[0][0] = a;
286 rBasis_[1][0] = b * cosGamma;
287 rBasis_[1][1] = b * sinGamma;
288 rBasis_[2][0] = c * sinTheta * cosPhi;
289 rBasis_[2][1] = c * sinTheta * sinPhi;
290 rBasis_[2][2] = c * cosTheta;
291
292 // Nonzero derivatives of basis vectors with respect to parameters
293 drBasis_[0](0,0) = 1.0;
294 drBasis_[1](1,0) = cosGamma;
295 drBasis_[1](1,1) = sinGamma;
296 drBasis_[2](2,0) = sinTheta * cosPhi;
297 drBasis_[2](2,1) = sinTheta * sinPhi;
298 drBasis_[2](2,2) = cosTheta;
299 drBasis_[3](2,0) = -c * sinTheta * sinPhi;
300 drBasis_[3](2,1) = c * sinTheta * cosPhi;
301 drBasis_[4](2,0) = c * cosTheta * cosPhi;
302 drBasis_[4](2,1) = c * cosTheta * sinPhi;
303 drBasis_[4](2,2) = -c * sinTheta;
304 drBasis_[5](1,0) = -b * sinGamma;
305 drBasis_[5](1,1) = b * cosGamma;
306
307 // Reciprocal lattice vectors
308 kBasis_[0][0] = twoPi / a;
309 kBasis_[0][1] = -twoPi * cosGamma / (a * sinGamma);
310 kBasis_[0][2] = sinTheta * (cosGamma * sinPhi - sinGamma * cosPhi);
311 kBasis_[0][2] = twoPi * kBasis_[0][2] / (a * sinGamma * cosTheta);
312
313 kBasis_[1][1] = twoPi/(b*sinGamma);
314 kBasis_[1][2] = -twoPi*sinPhi*sinTheta/(b*cosTheta*sinGamma);
315
316 kBasis_[2][2] = twoPi/(c*cosTheta);
317
318 } else {
319 UTIL_THROW("Unimplemented 3D lattice type");
320 }
321 }
322
323 /*
324 * Extract a UnitCell<3>::LatticeSystem from an istream as a string.
325 */
326 std::istream& operator >> (std::istream& in,
327 UnitCell<3>::LatticeSystem& lattice)
328 {
329
330 std::string buffer;
331 in >> buffer;
332 if (buffer == "Cubic" || buffer == "cubic") {
333 lattice = UnitCell<3>::Cubic;
334 } else
335 if (buffer == "Tetragonal" || buffer == "tetragonal") {
336 lattice = UnitCell<3>::Tetragonal;
337 } else
338 if (buffer == "Orthorhombic" || buffer == "orthorhombic") {
339 lattice = UnitCell<3>::Orthorhombic;
340 } else
341 if (buffer == "Monoclinic" || buffer == "monoclinic") {
342 lattice = UnitCell<3>::Monoclinic;
343 } else
344 if (buffer == "Triclinic" || buffer == "triclinic") {
345 lattice = UnitCell<3>::Triclinic;
346 } else
347 if (buffer == "Rhombohedral" || buffer == "rhombohedral") {
348 lattice = UnitCell<3>::Rhombohedral;
349 } else
350 if (buffer == "Hexagonal" || buffer == "hexagonal") {
351 lattice = UnitCell<3>::Hexagonal;
352 } else {
353 UTIL_THROW("Invalid UnitCell<3>::LatticeSystem value input");
354 }
355 return in;
356 }
357
358 /*
359 * Insert a UnitCell<3>::LatticeSystem to an ostream as a string.
360 */
361 std::ostream& operator << (std::ostream& out,
362 UnitCell<3>::LatticeSystem lattice)
363 {
364 if (lattice == UnitCell<3>::Cubic) {
365 out << "cubic";
366 } else
367 if (lattice == UnitCell<3>::Tetragonal) {
368 out << "tetragonal";
369 } else
370 if (lattice == UnitCell<3>::Orthorhombic) {
371 out << "orthorhombic";
372 } else
373 if (lattice == UnitCell<3>::Monoclinic) {
374 out << "monoclinic";
375 } else
376 if (lattice == UnitCell<3>::Triclinic) {
377 out << "triclinic";
378 } else
379 if (lattice == UnitCell<3>::Rhombohedral) {
380 out << "rhombohedral";
381 } else
382 if (lattice == UnitCell<3>::Hexagonal) {
383 out << "hexagonal";
384 } else
385 if (lattice == UnitCell<3>::Null) {
386 out << "Null";
387 } else {
388 UTIL_THROW("This should never happen");
389 }
390 return out;
391 }
392
393 /*
394 * Assignment operator.
395 */
396 UnitCell<3>& UnitCell<3>::operator = (UnitCell<3> const & other)
397 {
398 if (lattice_ != UnitCell<3>::Null) {
399 UTIL_CHECK(other.lattice_ == lattice_);
400 }
401 isInitialized_ = false;
402 lattice_ = other.lattice_;
403 setNParameter();
404 UTIL_CHECK(nParameter_ == other.nParameter_);
405 for (int i = 0; i < nParameter_; ++i) {
406 parameters_[i] = other.parameters_[i];
407 }
408 setLattice();
409 return *this;
410 }
411
412 /*
413 * Set lattice system of the unit cell (but not the parameters).
414 */
415 void UnitCell<3>::set(UnitCell<3>::LatticeSystem lattice)
416 {
417 UTIL_CHECK(lattice != UnitCell<3>::Null);
418 if (lattice_ != UnitCell<3>::Null) {
419 UTIL_CHECK(lattice == lattice_);
420 }
421 isInitialized_ = false;
422 lattice_ = lattice;
423 setNParameter();
424 }
425
426 /*
427 * Set state of the unit cell.
428 */
429 void UnitCell<3>::set(UnitCell<3>::LatticeSystem lattice,
430 FSArray<double, 6> const & parameters)
431 {
432 set(lattice);
433 UTIL_CHECK(parameters.size() == nParameter_);
434 for (int i = 0; i < nParameter_; ++i) {
435 parameters_[i] = parameters[i];
436 }
437 setLattice();
438 }
439
440 /*
441 * Get the length of the unit cell.
442 */
443 double UnitCell<3>::volume() const
444 {
445 double v = 0.0;
446 v += rBasis_[0][0]*rBasis_[1][1]*rBasis_[2][2];
447 v -= rBasis_[0][0]*rBasis_[1][2]*rBasis_[2][1];
448 v += rBasis_[0][1]*rBasis_[1][2]*rBasis_[2][0];
449 v -= rBasis_[0][1]*rBasis_[1][0]*rBasis_[2][2];
450 v += rBasis_[0][2]*rBasis_[1][0]*rBasis_[2][1];
451 v -= rBasis_[0][2]*rBasis_[1][1]*rBasis_[2][0];
452 return v;
453 }
454
455}
456}
Pscf::Prdc::UnitCellBase::parameters_
FArray< double, 6 > parameters_
Parameters used to describe the unit cell.
Definition UnitCellBase.h:268
Pscf::Prdc::UnitCellBase::nParameter_
int nParameter_
Number of parameters required to specify unit cell.
Definition UnitCellBase.h:273
Pscf::Prdc::UnitCell
Base template for UnitCell<D> classes, D=1, 2 or 3.
Definition rpg/System.h:34
Util::Constants::Pi
static const double Pi
Trigonometric constant Pi.
Definition Constants.h:35
Util::FSArray
A fixed capacity (static) contiguous array with a variable logical size.
Definition rpg/System.h:28
Util::FSArray::size
int size() const
Return logical size of this array (i.e., number of elements).
Definition FSArray.h:207
UTIL_CHECK
#define UTIL_CHECK(condition)
Assertion macro suitable for serial or parallel production code.
Definition global.h:68
UTIL_THROW
#define UTIL_THROW(msg)
Macro for throwing an Exception, reporting function, file and line number.
Definition global.h:51
Pscf
PSCF package top-level namespace.
Definition param_pc.dox:1
Util
Utility classes for scientific computation.
Definition accumulators.mod:1
Util::operator>>
std::istream & operator>>(std::istream &in, Pair< Data > &pair)
Input a Pair from an istream.
Definition Pair.h:44
Util::operator<<
std::ostream & operator<<(std::ostream &out, const Pair< Data > &pair)
Output a Pair to an ostream, without line breaks.
Definition Pair.h:57
Generated on Mon Apr 14 2025 16:29:24 for PSCF by doxygen 1.12.0