1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (c) 2008 The Regents of the University of California
4 //
5 // This file is part of Qbox
6 //
7 // Qbox is distributed under the terms of the GNU General Public License
8 // as published by the Free Software Foundation, either version 2 of
9 // the License, or (at your option) any later version.
10 // See the file COPYING in the root directory of this distribution
11 // or <http://www.gnu.org/licenses/>.
12 //
13 ////////////////////////////////////////////////////////////////////////////////
14 //
15 // SampleStepper.cpp
16 //
17 ////////////////////////////////////////////////////////////////////////////////
18 
19 #include "SampleStepper.h"
20 #include "Sample.h"
21 #include "Species.h"
22 
23 #include <iostream>
24 #include <iomanip>
25 using namespace std;
26 
27 ////////////////////////////////////////////////////////////////////////////////
SampleStepper(Sample & s)28 SampleStepper::SampleStepper(Sample& s) : s_(s)
29 {
30   fion.resize(s_.atoms.nsp());
31   for ( int is = 0; is < fion.size(); is++ )
32     fion[is].resize(3*s_.atoms.na(is));
33 
34   sigma_eks.resize(6);
35   sigma_kin.resize(6);
36   sigma_ext.resize(6);
37   sigma.resize(6);
38   sigma_ext = valarray<double>(s.ctrl.ext_stress,6);
39   const double gpa = 29421.5;
40   sigma_ext /= gpa;
41 }
42 
43 ////////////////////////////////////////////////////////////////////////////////
~SampleStepper(void)44 SampleStepper::~SampleStepper(void) {}
45 
46 ////////////////////////////////////////////////////////////////////////////////
print_stress(void)47 void SampleStepper::print_stress(void)
48 {
49   const double gpa = 29421.5;
50   if ( MPIdata::onpe0() )
51   {
52     cout.setf(ios::fixed,ios::floatfield);
53     cout.setf(ios::right,ios::adjustfield);
54     cout << setprecision(8);
55     cout << " <stress_tensor unit=\"GPa\">" << endl;
56 
57     cout << "   <sigma_eks_xx> " << setw(12) << sigma_eks[0]*gpa
58          << " </sigma_eks_xx>\n"
59          << "   <sigma_eks_yy> " << setw(12) << sigma_eks[1]*gpa
60          << " </sigma_eks_yy>\n"
61          << "   <sigma_eks_zz> " << setw(12) << sigma_eks[2]*gpa
62          << " </sigma_eks_zz>\n"
63          << "   <sigma_eks_xy> " << setw(12) << sigma_eks[3]*gpa
64          << " </sigma_eks_xy>\n"
65          << "   <sigma_eks_yz> " << setw(12) << sigma_eks[4]*gpa
66          << " </sigma_eks_yz>\n"
67          << "   <sigma_eks_xz> " << setw(12) << sigma_eks[5]*gpa
68          << " </sigma_eks_xz>\n\n";
69 
70     cout << "   <sigma_kin_xx> " << setw(12) << sigma_kin[0]*gpa
71          << " </sigma_kin_xx>\n"
72          << "   <sigma_kin_yy> " << setw(12) << sigma_kin[1]*gpa
73          << " </sigma_kin_yy>\n"
74          << "   <sigma_kin_zz> " << setw(12) << sigma_kin[2]*gpa
75          << " </sigma_kin_zz>\n"
76          << "   <sigma_kin_xy> " << setw(12) << sigma_kin[3]*gpa
77          << " </sigma_kin_xy>\n"
78          << "   <sigma_kin_yz> " << setw(12) << sigma_kin[4]*gpa
79          << " </sigma_kin_yz>\n"
80          << "   <sigma_kin_xz> " << setw(12) << sigma_kin[5]*gpa
81          << " </sigma_kin_xz>\n\n";
82 
83     cout << "   <sigma_ext_xx> " << setw(12) << sigma_ext[0]*gpa
84          << " </sigma_ext_xx>\n"
85          << "   <sigma_ext_yy> " << setw(12) << sigma_ext[1]*gpa
86          << " </sigma_ext_yy>\n"
87          << "   <sigma_ext_zz> " << setw(12) << sigma_ext[2]*gpa
88          << " </sigma_ext_zz>\n"
89          << "   <sigma_ext_xy> " << setw(12) << sigma_ext[3]*gpa
90          << " </sigma_ext_xy>\n"
91          << "   <sigma_ext_yz> " << setw(12) << sigma_ext[4]*gpa
92          << " </sigma_ext_yz>\n"
93          << "   <sigma_ext_xz> " << setw(12) << sigma_ext[5]*gpa
94          << " </sigma_ext_xz>\n\n";
95 
96     cout << "   <sigma_xx> " << setw(12) << sigma[0]*gpa
97          << " </sigma_xx>\n"
98          << "   <sigma_yy> " << setw(12) << sigma[1]*gpa
99          << " </sigma_yy>\n"
100          << "   <sigma_zz> " << setw(12) << sigma[2]*gpa
101          << " </sigma_zz>\n"
102          << "   <sigma_xy> " << setw(12) << sigma[3]*gpa
103          << " </sigma_xy>\n"
104          << "   <sigma_yz> " << setw(12) << sigma[4]*gpa
105          << " </sigma_yz>\n"
106          << "   <sigma_xz> " << setw(12) << sigma[5]*gpa
107          << " </sigma_xz>\n";
108 
109     cout << " </stress_tensor>" << endl;
110   }
111 }
112 
113 ////////////////////////////////////////////////////////////////////////////////
compute_sigma(void)114 void SampleStepper::compute_sigma(void)
115 {
116   sigma_kin = 0.0;
117   // compute kinetic contribution to stress using velocities at time t0
118   for ( int is = 0; is < s_.atoms.atom_list.size(); is++ )
119   {
120     int i = 0;
121     double mass = s_.atoms.species_list[is]->mass() * 1822.89;
122     for ( int ia = 0; ia < s_.atoms.atom_list[is].size(); ia++ )
123     {
124       Atom* pa = s_.atoms.atom_list[is][ia];
125       D3vector v = pa->velocity();
126       const double vx = v.x;
127       const double vy = v.y;
128       const double vz = v.z;
129 
130       sigma_kin[0] += mass * vx * vx;
131       sigma_kin[1] += mass * vy * vy;
132       sigma_kin[2] += mass * vz * vz;
133       sigma_kin[3] += mass * vx * vy;
134       sigma_kin[4] += mass * vy * vz;
135       sigma_kin[5] += mass * vx * vz;
136 
137       i += 3;
138     }
139   }
140   sigma_kin /= s_.wf.cell().volume();
141 
142   sigma = sigma_eks + sigma_kin - sigma_ext;
143 }
144