TorsionCmd.h 4.27 KB
Newer Older
1 2
////////////////////////////////////////////////////////////////////////////////
//
Francois Gygi committed
3 4 5 6
// Copyright (c) 2008 The Regents of the University of California
//
// This file is part of Qbox
//
Francois Gygi committed
7 8
// Qbox is distributed under the terms of the GNU General Public License
// as published by the Free Software Foundation, either version 2 of
Francois Gygi committed
9 10 11 12 13 14
// the License, or (at your option) any later version.
// See the file COPYING in the root directory of this distribution
// or <http://www.gnu.org/licenses/>.
//
////////////////////////////////////////////////////////////////////////////////
//
15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34
// TorsionCmd.h:
//
////////////////////////////////////////////////////////////////////////////////

#ifndef TORSIONCMD_H
#define TORSIONCMD_H

#include <iostream>
#include "UserInterface.h"
#include "Sample.h"
#include <cstdlib>

class TorsionCmd : public Cmd
{
  public:

  Sample *s;

  TorsionCmd(Sample *sample) : s(sample) {};

35 36
  const char *name(void) const { return "torsion"; }
  const char *help_msg(void) const
37
  {
38
    return
39
    "\n torsion\n\n"
40 41 42 43
    " syntax: torsion [-pbc] name1 name2 name3 name4\n\n"
    "   The torsion command prints the dihedral defined by four atoms.\n"
    "   If the -pbc option is used, the dihedral is computed using the\n"
    "   nearest atoms taking into account periodic boundary conditions.\n\n";
44 45 46 47
  }

  int action(int argc, char **argv)
  {
48
    if ( ! ( argc == 5 || argc == 6 ) )
49 50 51
    {
      if ( ui->onpe0() )
      {
52
        cout << " use: torsion [-pbc] name1 name2 name3 name4" << endl;
53 54 55
      }
      return 1;
    }
56

57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83
    string name1, name2, name3, name4;
    bool use_pbc = false;

    if ( argc == 5 )
    {
      name1 = argv[1];
      name2 = argv[2];
      name3 = argv[3];
      name4 = argv[4];
    }
    if ( argc == 6 )
    {
      if ( strcmp(argv[1],"-pbc") )
      {
        if ( ui->onpe0() )
        {
          cout << " use: torsion [-pbc] name1 name2 name3 name4" << endl;
        }
        return 1;
      }
      use_pbc = true;
      name1 = argv[2];
      name2 = argv[3];
      name3 = argv[4];
      name4 = argv[5];
    }

84 85 86 87 88 89
    Atom* a1 = s->atoms.findAtom(name1);
    Atom* a2 = s->atoms.findAtom(name2);
    Atom* a3 = s->atoms.findAtom(name3);
    Atom* a4 = s->atoms.findAtom(name4);
    if ( a1 == 0 || a2 == 0 || a3 == 0 || a4 == 0 )
    {
90 91 92
      if ( ui->onpe0() )
      {
        if ( a1 == 0 )
93
          cout << " TorsionCmd: atom " << name1 << " not defined" << endl;
94
        if ( a2 == 0 )
95
          cout << " TorsionCmd: atom " << name2 << " not defined" << endl;
96
        if ( a3 == 0 )
97
          cout << " TorsionCmd: atom " << name3 << " not defined" << endl;
98
        if ( a4 == 0 )
99
          cout << " TorsionCmd: atom " << name4 << " not defined" << endl;
100 101 102
      }
      return 1;
    }
103 104 105 106

    if ( a1 == a2 || a1 == a3 || a1 == a4 ||
         a2 == a3 || a2 == a4 || a3 == a4 )
    {
107 108
      if ( ui->onpe0() )
      {
109
        cout << " TorsionCmd: replicated atoms in "
110
             << name1 << " " << name2 << " "
111
             << name3 << " " << name4 << endl;
112 113 114
      }
      return 1;
    }
115

116
    if ( ui->onpe0() )
117
    {
118 119 120 121
      D3vector r12(a1->position()-a2->position());
      D3vector r32(a3->position()-a2->position());
      D3vector r43(a4->position()-a3->position());
      if ( norm2(r12) == 0.0 || norm2(r32) == 0.0 || norm2(r43) == 0.0 )
122
      {
123
        cout << " TorsionCmd: atoms are too close" << endl;
124
        return 1;
125
      }
126

127 128 129 130 131 132 133 134 135 136 137
      if ( use_pbc )
      {
        const UnitCell& cell = s->wf.cell();
        cell.fold_in_ws(r12);
        cell.fold_in_ws(r32);
        cell.fold_in_ws(r43);
      }
      D3vector e12(normalized(r12));
      D3vector e32(normalized(r32));
      D3vector e23(-e32);
      D3vector e43(normalized(r43));
138

139 140
      const double sin123 = length(e12^e32);
      const double sin234 = length(e23^e43);
141

142 143 144 145 146 147
      double a = 0;
      if ( sin123 != 0.0 && sin234 != 0.0 )
      {
        D3vector e123 = normalized(e12^e32);
        D3vector e234 = normalized(e23^e43);
        double cc = max(min(e123*e234,1.0),-1.0);
148
        double ss = max(min((e123^e234)*e32,1.0),-1.0);
149 150
        a = (180.0/M_PI) * atan2(ss,cc);
      }
151

152
      cout.setf(ios::fixed,ios::floatfield);
153
      cout << " torsion "
154 155
           << name1 << "-" << name2 << "-"
           << name3 << "-" << name4 << ": "
156
           << setprecision(3) << a << " (deg)" << endl;
157 158 159 160 161
    }
    return 0;
  }
};
#endif