New functions for AtomSet manipulation and measurement.

git-svn-id: http://qboxcode.org/svn/qb/trunk@415 cba15fb0-1239-40c8-b417-11db7ca47a34

src/AngleCmd.h

+////////////////////////////////////////////////////////////////////////////////
+//
+// AngleCmd.h:
+//
+////////////////////////////////////////////////////////////////////////////////
+// $Id: AngleCmd.h,v 1.1 2005-06-27 22:35:26 fgygi Exp $
+
+#ifndef ANGLECMD_H
+#define ANGLECMD_H
+
+#include <iostream>
+#include "UserInterface.h"
+#include "Sample.h"
+#include <cstdlib>
+
+class AngleCmd : public Cmd
+{
+  public:
+
+  Sample *s;
+
+  AngleCmd(Sample *sample) : s(sample) {};
+
+  char *name(void) const { return "angle"; }
+  char *help_msg(void) const
+  {
+    return 
+    "\n angle\n\n"
+    " syntax: angle name1 name2 name3\n\n"
+    "   The angle command prints the angle defined by three atoms.\n\n";
+  }
+
+  int action(int argc, char **argv)
+  {
+    if ( argc != 4 )
+    {
+      if ( ui->onpe0() )
+      {
+        cout << " use: angle name1 name2 name3" << endl;
+      }
+      return 1;
+    }
+    
+    string name1(argv[1]);                
+    string name2(argv[2]);                
+    string name3(argv[3]);                
+    Atom* a1 = s->atoms.findAtom(name1);  
+    Atom* a2 = s->atoms.findAtom(name2);  
+    Atom* a3 = s->atoms.findAtom(name3);  
+    if ( a1 == 0 || a2 == 0 || a3 == 0 )             
+    {                                     
+      if ( ui->onpe0() )
+      {
+        if ( a1 == 0 )
+          cout << " <!-- AngleCmd: atom " << name1 << " not defined -->"
+               << endl;
+        if ( a2 == 0 )
+          cout << " <!-- AngleCmd: atom " << name2 << " not defined -->"
+               << endl;
+        if ( a3 == 0 )
+          cout << " <!-- AngleCmd: atom " << name3 << " not defined -->"
+               << endl;
+      }
+      return 1;
+    }
+    
+    if ( a1 == a2 || a2 == a3 || a3 == a1 )             
+    {                                     
+      if ( ui->onpe0() )
+      {
+        cout << " <!-- AngleCmd: replicated atoms in " << name1   
+             << " " << name2 << " " << name3 << " -->" << endl;   
+      }
+      return 1;
+    }
+    
+    D3vector r12(a1->position()-a2->position());
+    D3vector r32(a3->position()-a2->position());
+    if ( norm(r12) == 0.0 || norm(r32) == 0.0 )
+    {                                     
+      if ( ui->onpe0() )
+      {
+        cout << " <!-- AngleCmd: atoms are too close -->" << endl;
+      }
+      return 1;
+    }
+    
+    const double sp = normalized(r12) * normalized(r32);
+    const double c = max(-1.0,min(1.0,sp));
+    const double a = (180.0/M_PI)*acos(c);
+    if ( ui->onpe0() )
+    {
+      cout.setf(ios::fixed,ios::floatfield);
+      cout << " <!-- angle " << name1 << "-" << name2  << "-" << name3
+           << ": "  
+           << setprecision(3)                                    
+           << a << " (deg) -->" << endl;     
+    }
+    return 0;
+  }
+};
+#endif

src/DistanceCmd.h

+////////////////////////////////////////////////////////////////////////////////
+//
+// DistanceCmd.h:
+//
+////////////////////////////////////////////////////////////////////////////////
+// $Id: DistanceCmd.h,v 1.1 2005-06-27 22:35:26 fgygi Exp $
+
+#ifndef DISTANCECMD_H
+#define DISTANCECMD_H
+
+#include <iostream>
+#include "UserInterface.h"
+#include "Sample.h"
+#include <cstdlib>
+
+class DistanceCmd : public Cmd
+{
+  public:
+
+  Sample *s;
+
+  DistanceCmd(Sample *sample) : s(sample) {};
+
+  char *name(void) const { return "distance"; }
+  char *help_msg(void) const
+  {
+    return 
+    "\n distance\n\n"
+    " syntax: distance name1 name2\n\n"
+    "   The distance command prints the distance between two atoms.\n\n";
+  }
+
+  int action(int argc, char **argv)
+  {
+    if ( argc != 3 )
+    {
+      if ( ui->onpe0() )
+      {
+        cout << " use: distance name1 name2" << endl;
+      }
+      return 1;
+    }
+    
+    string name1(argv[1]);                
+    string name2(argv[2]);                
+    Atom* a1 = s->atoms.findAtom(name1);  
+    Atom* a2 = s->atoms.findAtom(name2);  
+    if ( a1 == 0 || a2 == 0 )             
+    {                                     
+      // either a1 or a2 was not found
+      if ( ui->onpe0() )
+      {
+        if ( a1 == 0 )
+          cout << " <!-- DistanceCmd: atom " << name1 << " not defined -->"
+               << endl;
+        if ( a2 == 0 )
+          cout << " <!-- DistanceCmd: atom " << name2 << " not defined -->"
+               << endl;
+      }
+      return 1;
+    }
+    
+    if ( ui->onpe0() )
+    {
+      const double d = length(a1->position()-a2->position());
+      cout.setf(ios::fixed,ios::floatfield);
+      cout << " <!-- distance " << name1 << "-" << name2 << ": "  
+           << setprecision(3)                                    
+           << d << " (a.u.) / " << 0.529177*d << " (Ang)" << endl;     
+    }
+    return 0;
+  }
+};
+#endif

src/MoveCmd.h

+////////////////////////////////////////////////////////////////////////////////
+//
+// MoveCmd.h:
+//
+////////////////////////////////////////////////////////////////////////////////
+// $Id: MoveCmd.h,v 1.1 2005-06-27 22:35:26 fgygi Exp $
+
+#ifndef MOVECMD_H
+#define MOVECMD_H
+
+#include <string>
+#include <cstdlib>
+#include <iostream>
+using namespace std;
+
+#include "UserInterface.h"
+#include "Sample.h"
+
+class MoveCmd : public Cmd
+{
+  public:
+
+  Sample *s;
+
+  MoveCmd(Sample *sample) : s(sample) {};
+
+  char *name(void) const { return "move"; }
+  char *help_msg(void) const
+  {
+    return 
+    "\n move\n\n"
+    " syntax: move atom_name {to|by} x y z \n\n"
+    "   The move command displaces an atom to a new position.\n"
+    "   The new position is defined by absolute coordinates (to) or\n"
+    "   by a relative displacement (by).\n\n";
+  }
+
+  int action(int argc, char **argv)
+  {
+    // move must have 6 arguments including the command name
+    if ( argc != 6 )
+    {
+      if ( ui->onpe0() )
+        cout << " use: move atom_name {to|by} x y z " << endl;
+      return 1;
+    }
+  
+    const string atom_name = argv[1];
+    const string mode = argv[2];
+    const string xs = argv[3];
+    const string ys = argv[4];
+    const string zs = argv[5];
+
+    Atom* pa = s->atoms.findAtom(atom_name);
+    if ( !pa )
+    {
+      if ( ui->onpe0() )
+        cout << " <!-- MoveCmd: could not find atom " << atom_name << " -->"
+             << endl;
+      return 1;
+    }
+
+    D3vector pos = pa->position();
+
+    double x = pos.x;
+    double y = pos.y;
+    double z = pos.z;
+ 
+    if ( mode == "to" )
+    {
+      if ( xs != "*" )
+        x = atof(argv[3]);
+      if ( ys != "*" )
+        y = atof(argv[4]);
+      if ( zs != "*" )
+        z = atof(argv[5]);
+      pos = D3vector(x,y,z);      
+    }
+    else if ( mode == "by" )
+    {
+      x += atof(argv[3]);
+      y += atof(argv[4]);
+      z += atof(argv[5]);
+      pos = D3vector(x,y,z);      
+    }
+    else
+    {
+      if ( ui->onpe0() )
+        cout << " <!-- MoveCmd: unknown mode -->" << endl;
+      return 1;
+    }
+
+    pa->set_position(pos);
+    if ( ui->onpe0() )
+      cout << " <!-- MoveCmd: atom " << atom_name << " moved to "
+           << pos << " -->" << endl;
+    
+    return 0;
+  }
+};
+#endif

src/TorsionCmd.h

+////////////////////////////////////////////////////////////////////////////////
+//
+// TorsionCmd.h:
+//
+////////////////////////////////////////////////////////////////////////////////
+// $Id: TorsionCmd.h,v 1.1 2005-06-27 22:35:26 fgygi Exp $
+
+#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) {};
+
+  char *name(void) const { return "torsion"; }
+  char *help_msg(void) const
+  {
+    return 
+    "\n torsion\n\n"
+    " syntax: torsion name1 name2 name3 name4\n\n"
+    "   The torsion command prints the dihedral defined by four atoms.\n\n";
+  }
+
+  int action(int argc, char **argv)
+  {
+    if ( argc != 5 )
+    {
+      if ( ui->onpe0() )
+      {
+        cout << " use: torsion name1 name2 name3 name4" << endl;
+      }
+      return 1;
+    }
+    
+    string name1(argv[1]);                
+    string name2(argv[2]);                
+    string name3(argv[3]);                
+    string name4(argv[4]);                
+    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 )             
+    {                                     
+      if ( ui->onpe0() )
+      {
+        if ( a1 == 0 )
+          cout << " <!-- TorsionCmd: atom " << name1 << " not defined -->"
+               << endl;
+        if ( a2 == 0 )
+          cout << " <!-- TorsionCmd: atom " << name2 << " not defined -->"
+               << endl;
+        if ( a3 == 0 )
+          cout << " <!-- TorsionCmd: atom " << name3 << " not defined -->"
+               << endl;
+        if ( a4 == 0 )
+          cout << " <!-- TorsionCmd: atom " << name4 << " not defined -->"
+               << endl;
+      }
+      return 1;
+    }
+    
+    if ( a1 == a2 || a1 == a3 || a1 == a4 || 
+         a2 == a3 || a2 == a4 || a3 == a4 )             
+    {                                     
+      if ( ui->onpe0() )
+      {
+        cout << " <!-- TorsionCmd: replicated atoms in " 
+             << name1 << " " << name2 << " "   
+             << name3 << " " << name4 << " -->" << endl;   
+      }
+      return 1;
+    }
+    
+    D3vector r12(a1->position()-a2->position());
+    D3vector r32(a3->position()-a2->position());
+    D3vector r43(a4->position()-a3->position());
+    if ( norm(r12) == 0.0 || norm(r32) == 0.0 || norm(r43) == 0.0 )
+    {                                     
+      if ( ui->onpe0() )
+      {
+        cout << " <!-- TorsionCmd: atoms are too close -->" << endl;
+      }
+      return 1;
+    }
+    
+    D3vector e12(normalized(r12));
+    D3vector e32(normalized(r32));
+    D3vector e23(-e32);
+    D3vector e43(normalized(r43));
+    
+    const double sin123 = length(e12^e32);
+    const double sin234 = length(e23^e43);
+    
+    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);
+      double ss = max(min((e123^e234)*e32,1.0),-1.0);
+      a = (180.0/M_PI) * atan2(ss,cc);
+    }
+    
+    if ( ui->onpe0() )
+    {
+      cout.setf(ios::fixed,ios::floatfield);
+      cout << " <!-- torsion " 
+           << name1 << "-" << name2 << "-"  
+           << name3 << "-" << name4 << ": "  
+           << setprecision(3) << a << " (deg) -->" << endl;                                    
+    }
+    return 0;
+  }
+};
+#endif