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 // AngleCmd.h:
16 //
17 ////////////////////////////////////////////////////////////////////////////////
18 
19 #ifndef ANGLECMD_H
20 #define ANGLECMD_H
21 
22 #include <iostream>
23 #include "UserInterface.h"
24 #include "Sample.h"
25 #include <cstdlib>
26 
27 class AngleCmd : public Cmd
28 {
29   public:
30 
31   Sample *s;
32 
AngleCmd(Sample * sample)33   AngleCmd(Sample *sample) : s(sample) {};
34 
name(void)35   const char *name(void) const { return "angle"; }
help_msg(void)36   const char *help_msg(void) const
37   {
38     return
39     "\n angle\n\n"
40     " syntax: angle [-pbc] name1 name2 name3\n\n"
41     "   The angle command prints the angle defined by three atoms.\n"
42     "   If the -pbc option is used, the angle is computed using the\n"
43     "   nearest atoms taking into account periodic boundary conditions.\n\n";
44   }
45 
action(int argc,char ** argv)46   int action(int argc, char **argv)
47   {
48     if ( ! ( argc == 4 || argc == 5 ) )
49     {
50       if ( ui->onpe0() )
51       {
52         cout << " use: angle [-pbc] name1 name2 name3" << endl;
53       }
54       return 1;
55     }
56 
57     string name1, name2, name3;
58     bool use_pbc = false;
59 
60     if ( argc == 4 )
61     {
62       name1 = argv[1];
63       name2 = argv[2];
64       name3 = argv[3];
65     }
66     if ( argc == 5 )
67     {
68       if ( strcmp(argv[1],"-pbc") )
69       {
70         if ( ui->onpe0() )
71         {
72           cout << " use: angle [-pbc] name1 name2 name3" << endl;
73         }
74         return 1;
75       }
76       use_pbc = true;
77       name1 = argv[2];
78       name2 = argv[3];
79       name3 = argv[4];
80     }
81 
82     Atom* a1 = s->atoms.findAtom(name1);
83     Atom* a2 = s->atoms.findAtom(name2);
84     Atom* a3 = s->atoms.findAtom(name3);
85     if ( a1 == 0 || a2 == 0 || a3 == 0 )
86     {
87       if ( ui->onpe0() )
88       {
89         if ( a1 == 0 )
90           cout << " AngleCmd: atom " << name1 << " not defined" << endl;
91         if ( a2 == 0 )
92           cout << " AngleCmd: atom " << name2 << " not defined" << endl;
93         if ( a3 == 0 )
94           cout << " AngleCmd: atom " << name3 << " not defined" << endl;
95       }
96       return 1;
97     }
98 
99     if ( a1 == a2 || a2 == a3 || a3 == a1 )
100     {
101       if ( ui->onpe0() )
102       {
103         cout << " AngleCmd: replicated atoms in " << name1
104              << " " << name2 << " " << name3 << endl;
105       }
106       return 1;
107     }
108 
109     if ( ui->onpe0() )
110     {
111       D3vector r12(a1->position()-a2->position());
112       D3vector r32(a3->position()-a2->position());
113       if ( norm2(r12) == 0.0 || norm2(r32) == 0.0 )
114       {
115         cout << " AngleCmd: atoms are too close" << endl;
116         return 1;
117       }
118 
119       if ( use_pbc )
120       {
121         const UnitCell& cell = s->wf.cell();
122         cell.fold_in_ws(r12);
123         cell.fold_in_ws(r32);
124       }
125       const double sp = normalized(r12) * normalized(r32);
126       const double c = max(-1.0,min(1.0,sp));
127       const double a = (180.0/M_PI)*acos(c);
128       cout.setf(ios::fixed,ios::floatfield);
129       cout << " angle " << name1 << "-" << name2  << "-" << name3
130            << ": "
131            << setprecision(3)
132            << a << " (deg)" << endl;
133     }
134     return 0;
135   }
136 };
137 #endif
138