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)28SampleStepper::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)44SampleStepper::~SampleStepper(void) {} 45 46 //////////////////////////////////////////////////////////////////////////////// print_stress(void)47void 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)114void 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