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 // GlobalExtForce.cpp 16 // 17 //////////////////////////////////////////////////////////////////////////////// 18 19 #include "GlobalExtForce.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 GlobalExtForce::setup(const AtomSet& atoms) 31 { 32 // All atoms are affected, no setup needed 33 na_ = atoms.size(); 34 } 35 36 //////////////////////////////////////////////////////////////////////////////// energy(const vector<vector<double>> & r,vector<vector<double>> & f)37double GlobalExtForce::energy(const vector<vector<double> > &r, 38 vector<vector<double> > &f) 39 { 40 double sum = 0.0; 41 for ( int is = 0; is < f.size(); is++ ) 42 { 43 const int nais = f[is].size() / 3; 44 for ( int ia = 0; ia < nais; ia++ ) 45 { 46 f[is][3*ia+0] += force_.x / na_; 47 f[is][3*ia+1] += force_.y / na_; 48 f[is][3*ia+2] += force_.z / na_; 49 sum -= force_.x * r[is][3*ia+0] / na_ + 50 force_.y * r[is][3*ia+1] / na_ + 51 force_.z * r[is][3*ia+2] / na_; 52 } 53 } 54 return sum; 55 } 56 57 //////////////////////////////////////////////////////////////////////////////// print(ostream & os)58ostream& GlobalExtForce::print( ostream &os ) 59 { 60 os.setf(ios::left,ios::adjustfield); 61 os << " <extforce name=\"" << name(); 62 os << "\" type=\"" << type(); 63 os << "\" atoms=\"_all_\"\n"; 64 os.setf(ios::fixed,ios::floatfield); 65 os.setf(ios::right,ios::adjustfield); 66 os << " value=\"" << setprecision(6) << force_.x << " " 67 << setprecision(6) << force_.y << " " 68 << setprecision(6) << force_.z << "\"/>"; 69 return os; 70 } 71