1 /********************************************************************** 2 Copyright (C) 2011 by Geoffrey Hutchison 3 4 This program is free software; you can redistribute it and/or modify 5 it under the terms of the GNU General Public License as published by 6 the Free Software Foundation version 2 of the License. 7 8 This program is distributed in the hope that it will be useful, 9 but WITHOUT ANY WARRANTY; without even the implied warranty of 10 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 11 GNU General Public License for more details. 12 ***********************************************************************/ 13 14 #include <openbabel/babelconfig.h> 15 16 #include <openbabel/obmolecformat.h> 17 #include <openbabel/mol.h> 18 #include <openbabel/atom.h> 19 #include <openbabel/bond.h> 20 #include <openbabel/obiter.h> 21 #include <openbabel/elements.h> 22 #include <openbabel/generic.h> 23 #include <cstdlib> 24 25 26 using namespace std; 27 namespace OpenBabel 28 { 29 30 class XSFFormat : public OBMoleculeFormat 31 { 32 public: 33 //Register this format type ID XSFFormat()34 XSFFormat() 35 { 36 OBConversion::RegisterFormat("xsf",this); 37 // animation variant 38 OBConversion::RegisterFormat("axsf",this); 39 } 40 Description()41 virtual const char* Description() //required 42 { 43 return 44 "XCrySDen Structure Format\n" 45 "Read Options e.g. -as\n" 46 " s Output single bonds only\n" 47 " b Disable bonding entirely\n\n"; 48 }; 49 SpecificationURL()50 virtual const char* SpecificationURL() 51 { return "http://www.xcrysden.org/doc/XSF.html/" ;}; //optional 52 53 //Flags() can return be any the following combined by | or be omitted if none apply 54 // NOTREADABLE READONEONLY NOTWRITABLE WRITEONEONLY Flags()55 virtual unsigned int Flags() 56 { 57 return READONEONLY | NOTWRITABLE; 58 }; 59 60 //*** This section identical for most OBMol conversions *** 61 //////////////////////////////////////////////////// 62 /// The "API" interface functions 63 virtual bool ReadMolecule(OBBase* pOb, OBConversion* pConv); 64 }; 65 //*** 66 67 //Make an instance of the format class 68 XSFFormat theXSFFormat; 69 70 ///////////////////////////////////////////////////////////////// ReadMolecule(OBBase * pOb,OBConversion * pConv)71 bool XSFFormat::ReadMolecule(OBBase* pOb, OBConversion* pConv) 72 { 73 74 OBMol* pmol = pOb->CastAndClear<OBMol>(); 75 if (pmol == nullptr) 76 return false; 77 78 //Define some references so we can use the old parameter names 79 istream &ifs = *pConv->GetInStream(); 80 OBMol &mol = *pmol; 81 const char* title = pConv->GetTitle(); 82 83 char buffer[BUFF_SIZE]; 84 string str; 85 double x,y,z; 86 OBAtom *atom; 87 vector3 translationVectors[3]; 88 int numTranslationVectors = 0; 89 vector<string> vs; 90 vector<vector3> atomPositions; 91 bool createdAtoms = false; 92 int atomicNum; 93 94 mol.BeginModify(); 95 96 while (ifs.getline(buffer, BUFF_SIZE)) 97 { 98 if (buffer[0] == '#') 99 continue; // comment 100 if (strstr(buffer, "ATOMS") != nullptr) { 101 // Minimum of 4 columns -- AtNum, x, y, z (forces) 102 // where AtNum stands for atomic number (or symbol), while X Y Z are 103 ifs.getline(buffer, BUFF_SIZE); 104 tokenize(vs, buffer); 105 while (vs.size() >= 4) { 106 if (!createdAtoms) { 107 atom = mol.NewAtom(); 108 //set atomic number 109 atomicNum = OBElements::GetAtomicNum(vs[0].c_str()); 110 if (atomicNum == 0) { 111 atomicNum = atoi(vs[0].c_str()); 112 } 113 atom->SetAtomicNum(atomicNum); 114 } 115 x = atof((char*)vs[1].c_str()); 116 y = atof((char*)vs[2].c_str()); 117 z = atof((char*)vs[3].c_str()); 118 atomPositions.push_back(vector3(x, y, z)); // we may have a movie or animation 119 120 ifs.getline(buffer, BUFF_SIZE); 121 tokenize(vs, buffer); 122 } 123 createdAtoms = true; // don't run NewAtom() anymore 124 } 125 else if ( strstr(buffer, "PRIMVEC") 126 || strstr(buffer, "CONVVEC") ) { 127 // translation vectors 128 numTranslationVectors = 0; // if we have an animation 129 while (numTranslationVectors < 3 && ifs.getline(buffer,BUFF_SIZE)) { 130 tokenize(vs,buffer); // we really need to check that it's 3 entries only 131 if (vs.size() < 3) return false; // timvdm 18/06/2008 132 x = atof((char*)vs[0].c_str()); 133 y = atof((char*)vs[1].c_str()); 134 z = atof((char*)vs[2].c_str()); 135 translationVectors[numTranslationVectors++].Set(x, y, z); 136 } 137 } 138 else if (strstr(buffer, "PRIMCOORD") != nullptr) { 139 // read the coordinates 140 ifs.getline(buffer, BUFF_SIZE); 141 tokenize(vs, buffer); 142 if (vs.size() < 2) return false; 143 int numAtoms = atoi(vs[0].c_str()); 144 for (int a = 0; a < numAtoms; ++a) { 145 if (!ifs.getline(buffer,BUFF_SIZE)) 146 break; 147 tokenize(vs,buffer); 148 if (vs.size() < 4) 149 break; 150 151 if (!createdAtoms) { 152 atom = mol.NewAtom(); 153 //set atomic number 154 atomicNum = OBElements::GetAtomicNum(vs[0].c_str()); 155 if (atomicNum == 0) { 156 atomicNum = atoi(vs[0].c_str()); 157 } 158 atom->SetAtomicNum(atomicNum); 159 } 160 x = atof((char*)vs[1].c_str()); 161 y = atof((char*)vs[2].c_str()); 162 z = atof((char*)vs[3].c_str()); 163 atomPositions.push_back(vector3(x, y, z)); 164 } 165 } 166 } 167 168 mol.EndModify(); 169 170 int natom = mol.NumAtoms(); 171 int numConformers = atomPositions.size() / natom; 172 for (int i = 0; i < numConformers; ++i) { 173 double *coordinates = new double[natom * 3]; 174 for (int j = 0; j < natom; ++j) { 175 vector3 currentPosition = atomPositions[i*natom + j]; 176 coordinates[j*3] = currentPosition.x(); 177 coordinates[j*3 + 1] = currentPosition.y(); 178 coordinates[j*3 + 2] = currentPosition.z(); 179 } 180 mol.AddConformer(coordinates); 181 } 182 // Delete first conformer, created by EndModify, bunch of 0s 183 mol.DeleteConformer(0); 184 // Set geometry to last one 185 mol.SetConformer(mol.NumConformers() - 1); 186 187 if (!pConv->IsOption("b",OBConversion::INOPTIONS)) 188 mol.ConnectTheDots(); 189 if (!pConv->IsOption("s",OBConversion::INOPTIONS) && !pConv->IsOption("b",OBConversion::INOPTIONS)) 190 mol.PerceiveBondOrders(); 191 192 // Add final properties 193 mol.SetTitle(title); 194 if (numTranslationVectors == 3) { 195 OBUnitCell *uc = new OBUnitCell; 196 uc->SetOrigin(fileformatInput); 197 uc->SetData(translationVectors[0], 198 translationVectors[1], 199 translationVectors[2]); 200 mol.SetData(uc); 201 } 202 203 return(true); 204 } 205 206 } //namespace OpenBabel 207