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