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)30 void 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)40 double 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)52 ostream& 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