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)30 void 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)37 double 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)58 ostream& 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