openbabel/tools/obxyz.cpp

/**********************************************************************
obxyz - Open Babel XYZ transformations

Copyright (C) 2001-2006 Geoffrey R. Hutchison
 
This file is part of the Open Babel project.
For more information, see <http://openbabel.sourceforge.net/>

This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation version 2 of the License.
 
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
***********************************************************************/

#include "mol.h"
#include "math/matrix3x3.h"
#include "obconversion.h"

using namespace std;
using namespace OpenBabel;

int main(int argc,char **argv)
{
    char c;
    int transX, transY, transZ, centerIdx;
    transX = transY = transZ = centerIdx = -1;

    char *FileIn = NULL;
    char *program_name = argv[0];
    char *iext;

    OBConversion conv(&cin,&cout);
    OBFormat *pFormat = conv.FindFormat("smi"); // default format is SMILES

    // Still need to add:
    // rotate X, Y, or Z by set angle
    // translate X, Y, or Z by set amount
    
    // Parse options
    while ((c = getopt(argc, argv, "x:y:z:c:")) != -1)
    {
        switch (c)
        {
        case 'c': /// atom to be centered

            c = sscanf(optarg, "%d", &centerIdx);
            if (c != 1 )
            {
                cerr << program_name << ": unable to parse -c option" << endl;
                exit (-1);
            }
            break;

        case 'x': /// atom to be centered

            c = sscanf(optarg, "%d", &transX);
            if (c != 1 )
            {
                cerr << program_name << ": unable to parse -x option" << endl;
                exit (-1);
            }
            break;

        case 'y': /// atom to y-axis

            c = sscanf(optarg, "%d", &transY);
            if (c != 1 )
            {
                cerr << program_name << ": unable to parse -y option" << endl;
                exit (-1);
            }
            break;

        case 'z': /// atom to z-axis

            c = sscanf(optarg, "%d", &transZ);
            if (c != 1 )
            {
                cerr << program_name << ": unable to parse -z option" << endl;
                exit (-1);
            }
            break;

        }
    }
    
    ifstream ifs;
    FileIn  = argv[optind];
    if (FileIn != NULL)
      {
	// Read the file
	ifs.open(FileIn);
	if (!ifs)
	  {
	    cerr << program_name << ": cannot read input file!" << endl;
	    exit (-1);
	  }
	conv.SetInStream(&ifs);
	
	// Find Input filetype
	pFormat = conv.FormatFromExt(FileIn);
	if (pFormat == NULL)
	  {
	    cerr << program_name << ": cannot read input format!" << endl;
	    return (-1);
	  }
      }

    if (! conv.SetInAndOutFormats(pFormat, pFormat))
      {
	cerr << program_name << ": cannot read or write to this file format" << endl;
	return (-1);
      }

    OBMol mol;
    mol.Clear();
    conv.Read(&mol);
    if (mol.Empty())
      return(1);

    vector3 v;
    if (centerIdx != -1)
      {
	v = (mol.GetAtom(centerIdx))->GetVector();
	mol.Translate(-1.0f*v);
      }

    double alpha, beta, gamma;
    alpha = beta = gamma = 0.0f;
    if (transX != -1)
      {
	v = (mol.GetAtom(transX))->GetVector();

	// angle to rotate vector into XZ plane
	// (i.e., angle between x and y components)
	gamma = PI/2.0 - atan(v.x() / v.y());

	// angle to rotate vector from XZ plane to X axis 
	// (by rotation along Y-axis)
	beta = atan(v.z() / sqrt(v.y()*v.y() + v.x()*v.x()) );
      }
    else if (transY != -1)
      {
	v = (mol.GetAtom(transY))->GetVector();

	// angle to rotate vector into YZ plane
	// (i.e., angle between x and y components)
	gamma = -1.0 * atan(v.x() / v.y());

	// angle to rotate vector from YZ plane to Y axis 
	// (by rotation along X-axis)
	alpha = -1.0 * atan(v.z() / sqrt(v.y()*v.y() + v.x()*v.x()) );
      }
    else if (transZ != -1)
      {
	v = (mol.GetAtom(transZ))->GetVector();

	// angle to rotate vector into YZ plane
	// (i.e., angle between x and y components)
	gamma = -1.0 * atan(v.x() / v.y());

	// angle to rotate vector from YZ plane to Z axis
	// (by rotation along X axis)
	alpha = PI/2.0f - atan(v.z() / sqrt(v.y()*v.y() + v.x()*v.x()) );
      }

    if ( !IsNearZero(alpha) ||
	 !IsNearZero(beta)  ||
	 !IsNearZero(gamma) )
      {
	matrix3x3 mat;
	mat.SetupRotMat(alpha*RAD_TO_DEG, beta*RAD_TO_DEG, gamma*RAD_TO_DEG);
	double array[9];
	mat.GetArray(array);
	mol.Rotate(array);
      }

    conv.Write(&mol, &cout);

    return(0);
}