I think this works OK. Ron

On Mon, Apr 29, 2013 at 9:21 PM, Cohen, Ron <rco...@carnegiescience.edu> wrote:
> I will be sending a working version to you by tomorrow, I hope. Ron
>
>
>
> On Mon, Apr 29, 2013 at 4:44 PM, Cohen, Ron <rco...@carnegiescience.edu> 
> wrote:
>> The atom positions are wrong too. I have been working on it, but not
>> sure how to fix the multiple line thing, and not sure how
>> mol>Addconformer works. I can get correct atomPositions but something
>> is still messed up. Ron
>>
>>
>> On Mon, Apr 29, 2013 at 2:52 PM, Cohen, Ron <rco...@carnegiescience.edu> 
>> wrote:
>>> It still has some problems. Converts to (in cif format):
>>>
>>> _chemical_name_common 'BTO_Mn_B_V_z_PAW.abinit'
>>> _cell_length_a 0
>>> _cell_length_b nan
>>> _cell_length_c nan
>>> _cell_angle_alpha nan
>>> _cell_angle_beta nan
>>> _cell_angle_gamma nan
>>> _space_group_name_H-M_alt 'P 4 m m'
>>> _space_group_name_Hall 'P 4 -2'
>>> ...
>>> whole cif file attached.
>>>
>>> Thanks!
>>>
>>> Ron
>>>
>>>
>>>
>>> On Mon, Apr 29, 2013 at 2:44 PM, Cohen, Ron <rco...@carnegiescience.edu> 
>>> wrote:
>>>> Sure. It is attached. The problem is in tokenize, I think, so it seems
>>>> it might be more general. If I take out the newline to
>>>>            typat      1  1  1  1  1  1  1  1  4  2  2  2  2  2  2  2
>>>> 3  3  3  3 3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3
>>>> It seems to run OK.
>>>>
>>>> Ron
>>>>
>>>>
>>>> On Mon, Apr 29, 2013 at 2:41 PM, David Lonie <lonieda...@gmail.com> wrote:
>>>>> (Moving this to the openbabel developers list)
>>>>>
>>>>> Can you provide a full sample file that demonstrates the issue?
>>>>>
>>>>> Thanks,
>>>>> Dave
>>>>>
>>>>> On Mon, Apr 29, 2013 at 2:39 PM, Cohen, Ron <rco...@carnegiescience.edu> 
>>>>> wrote:
>>>>>> It fails when the typeat goes over two lines, as in:
>>>>>>
>>>>>>           typat      1  1  1  1  1  1  1  1  4  2  2  2  2  2  2  2  3  
>>>>>> 3  3  3
>>>>>>                        3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3 
>>>>>>  3  3
/**********************************************************************
Copyright (C) 2011 by Geoffrey Hutchison

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 <openbabel/babelconfig.h>

#include <openbabel/obmolecformat.h>

using namespace std;
namespace OpenBabel
{

#define BOHR_TO_ANGSTROM 0.5291772108
#define ANGSTROM_TO_BOHR 1.0 / BOHR_TO_ANGSTROM


  class ABINITFormat : public OBMoleculeFormat
  {
  public:
    //Register this format type ID
    ABINITFormat()
    {
      OBConversion::RegisterFormat("abinit",this);
    }

    virtual const char* Description() //required
    {
      return
        "ABINIT Output Format\n"
        "Read Options e.g. -as\n"
        "  s  Output single bonds only\n"
        "  b  Disable bonding entirely\n\n";
    };

    virtual const char* SpecificationURL()
    { return "http://abinit.org/"; ;}; //optional

    //Flags() can return be any the following combined by | or be omitted if none apply
    // NOTREADABLE  READONEONLY  NOTWRITABLE  WRITEONEONLY
    virtual unsigned int Flags()
    {
      return READONEONLY | NOTWRITABLE;
    };

    //*** This section identical for most OBMol conversions ***
    ////////////////////////////////////////////////////
    /// The "API" interface functions
    virtual bool ReadMolecule(OBBase* pOb, OBConversion* pConv);
  };
  //***

  //Make an instance of the format class
  ABINITFormat theABINITFormat;

  /////////////////////////////////////////////////////////////////
  bool ABINITFormat::ReadMolecule(OBBase* pOb, OBConversion* pConv)
  {

    OBMol* pmol = pOb->CastAndClear<OBMol>();
    if(pmol==NULL)
      return false;

    //Define some references so we can use the old parameter names
    istream &ifs = *pConv->GetInStream();
    OBMol &mol = *pmol;
    const char* title = pConv->GetTitle();

    char buffer[BUFF_SIZE];
    string str;
    vector<string> vs;

    OBAtom *atom;
    int natom;
    vector<int> atomicNumbers, atomTypes;
    double x, y, z;
    vector<vector3> atomPositions;
    vector<double> energies;

    // Translation vectors (if present)
    // aka rprim
    vector3 translationVectors[3];
    double acell[3]; // scale of lattice vectors
    int numTranslationVectors = 0;
    int symmetryCode = 0;
    bool last = false;
    mol.BeginModify();

    while (ifs.getline(buffer,BUFF_SIZE))
      {
	if(strstr(buffer,"outvars")){
	     last = true;
	     continue;
	   }
	if(! last){
	  continue;
	}
        // tokens are listed in alphabetical order for code clarity
        if (strstr(buffer, "acell")) {
          tokenize(vs, buffer);
          if (vs.size() < 4)
            continue; // invalid line

          // acell=  7.6967369631E+00  7.6967369631E+00  7.6967369631E+00
          for (int i = 0; i < 3; ++i) {
            acell[i] = atof(vs[i+1].c_str());
          }
        }
        else if (strstr(buffer, " xcart ")) {
          double unit = BOHR_TO_ANGSTROM;
          if (strstr(buffer, "ngstrom"))
            unit = 1.0; // no conversion needed

	  //          ifs.getline(buffer,BUFF_SIZE);
          tokenize(vs, buffer);

          // first line, rprim takes up a token
          x = atof((char*)vs[1].c_str()) * unit;
          y = atof((char*)vs[2].c_str()) * unit;
	  z = atof((char*)vs[3].c_str()) * unit;
	  atomPositions.push_back(vector3(x, y, z));
	  // get next line
	  ifs.getline(buffer,BUFF_SIZE);
	  tokenize(vs, buffer);
	  //rest of lines
	  while(vs.size() == 3) {
	    x = atof(vs[0].c_str()) * unit;
	    y = atof(vs[1].c_str()) * unit;
	    z = atof(vs[2].c_str()) * unit;
            atomPositions.push_back(vector3(x, y, z));
            // get next line
            ifs.getline(buffer,BUFF_SIZE);
            tokenize(vs, buffer);
          }
      }
        else if (strstr(buffer, "natom")) {
          tokenize(vs, buffer);
          if (vs.size() != 2)
            continue;
          natom = atoi(vs[1].c_str());
        }
        else if (strstr(buffer, "rprim")) {
          numTranslationVectors = 0;
          int column;
	  ifs.getline(buffer,BUFF_SIZE);
          for (int i = 1; i < 4; ++i) {
            tokenize(vs, buffer);
            if (vs.size() < 3)
              break;

            // first line, rprim takes up a token
	    //            if (i == 0)
	    // column = 1;
            //else
              column = 0;

            x = atof((char*)vs[column].c_str()) * BOHR_TO_ANGSTROM;
            y = atof((char*)vs[column+1].c_str()) * BOHR_TO_ANGSTROM;
            z = atof((char*)vs[column+2].c_str()) * BOHR_TO_ANGSTROM;
            translationVectors[numTranslationVectors++].Set(x, y,z);
            ifs.getline(buffer,BUFF_SIZE);
          }
        }
        else if (strstr(buffer, "Symmetries")) {
          tokenize(vs, buffer, "()");
          // Should be something like (#160)
          symmetryCode = atoi(vs[1].substr(1).c_str());
        }
        else if (strstr(buffer, "typat")) {
          atomTypes.clear();
	  int n=0;
	    while( n<=natom){
	      tokenize(vs, buffer);
	      for (unsigned int i = 1; i < vs.size(); ++i) {
		atomTypes.push_back(atoi(vs[i].c_str()));
	      }
	      n+=vs.size();
	      ifs.getline(buffer,BUFF_SIZE);
	    }
        }
        else if (strstr(buffer, "znucl")) {
          tokenize(vs, buffer);
          // make sure znucl is first token
          if (vs[0] != "znucl")
            continue;
          // push back the remaining tokens into atomicNumbers
          atomicNumbers.clear();
          atomicNumbers.push_back(0); // abinit starts typat with 1
          for (unsigned int i = 1; i < vs.size(); ++i)
            atomicNumbers.push_back(int(atof(vs[i].c_str())));
        }
        // xangst
        // forces
      }

    for (int i = 0; i < natom; ++i) {
      atom = mol.NewAtom();
      //set atomic number
      int idx = atom->GetIdx();
      int type = atomTypes[idx - 1];
      atom->SetAtomicNum(atomicNumbers[type]);
      // we set the coordinates by conformers in another loop
    }

    mol.EndModify();

    int numConformers = atomPositions.size() / natom;
    for (int i = 0; i < numConformers; ++i) {
      double *coordinates = new double[natom * 3];
      for (int j = 0; j < natom; ++j) {
        vector3 currentPosition = atomPositions[i*natom + j];
        coordinates[j*3] = currentPosition.x();
        coordinates[j*3 + 1] = currentPosition.y();
        coordinates[j*3 + 2] = currentPosition.z();
      }
      mol.AddConformer(coordinates);
    }
    // Delete first conformer, created by EndModify, bunch of 0s
    mol.DeleteConformer(0);
    // Set geometry to last one
    mol.SetConformer(mol.NumConformers() - 1);

    if (!pConv->IsOption("b",OBConversion::INOPTIONS))
      mol.ConnectTheDots();
    if (!pConv->IsOption("s",OBConversion::INOPTIONS) && !pConv->IsOption("b",OBConversion::INOPTIONS))
      mol.PerceiveBondOrders();

    // Attach unit cell translation vectors if found
    if (numTranslationVectors > 0) {
      OBUnitCell* uc = new OBUnitCell;
	    //      uc->SetData(acell[0] * translationVectors[0], acell[1] * translationVectors[1], acell[2] * translationVectors[2]);
      uc->SetData(translationVectors[0], translationVectors[1], translationVectors[2]);
      uc->SetOrigin(fileformatInput);
      if (symmetryCode)
        uc->SetSpaceGroup(symmetryCode);
      mol.SetData(uc);
    }

    mol.SetTitle(title);
    return(true);
  }

} //namespace OpenBabel
------------------------------------------------------------------------------
Introducing AppDynamics Lite, a free troubleshooting tool for Java/.NET
Get 100% visibility into your production application - at no cost.
Code-level diagnostics for performance bottlenecks with <2% overhead
Download for free and get started troubleshooting in minutes.
http://p.sf.net/sfu/appdyn_d2d_ap1
_______________________________________________
OpenBabel-Devel mailing list
OpenBabel-Devel@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/openbabel-devel

Reply via email to