1 //////////////////////////////////////////////////////////////////////////////// 2 // 3 // Copyright (c) 2009 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 // AtomicExtForce.cpp 16 // 17 //////////////////////////////////////////////////////////////////////////////// 18 19 #include "AtomicExtForce.h" 20 #include "AtomSet.h" 21 #include "Atom.h" 22 #include "Species.h" 23 #include <cassert> 24 #include <cmath> 25 #include <iostream> 26 #include <iomanip> 27 using namespace std; 28 29 //////////////////////////////////////////////////////////////////////////////// setup(const AtomSet & atoms)30void AtomicExtForce::setup(const AtomSet& atoms) 31 { 32 // find position in tau array corresponding to atom name1 33 is1_ = atoms.is(name1_); 34 ia1_ = atoms.ia(name1_); 35 assert(is1_>=0); 36 assert(ia1_>=0); 37 } 38 39 //////////////////////////////////////////////////////////////////////////////// energy(const vector<vector<double>> & r,vector<vector<double>> & f)40double AtomicExtForce::energy(const vector<vector<double> > &r, 41 vector<vector<double> > &f) 42 { 43 f[is1_][3*ia1_+0] += force_.x; 44 f[is1_][3*ia1_+1] += force_.y; 45 f[is1_][3*ia1_+2] += force_.z; 46 return -force_.x*r[is1_][3*ia1_+0] + 47 force_.y*r[is1_][3*ia1_+1] + 48 force_.z*r[is1_][3*ia1_+2]; 49 } 50 51 //////////////////////////////////////////////////////////////////////////////// print(ostream & os)52ostream& AtomicExtForce::print( ostream &os ) 53 { 54 os.setf(ios::left,ios::adjustfield); 55 os << " <extforce name=\"" << name(); 56 os << "\" type=\"" << type(); 57 os << "\" atoms=\"" << name1_ << "\"\n"; 58 os.setf(ios::fixed,ios::floatfield); 59 os.setf(ios::right,ios::adjustfield); 60 os << " value=\"" << setprecision(6) << force_.x << " " 61 << setprecision(6) << force_.y << " " 62 << setprecision(6) << force_.z << "\"/>"; 63 return os; 64 } 65