/src/openbabel/src/formats/pwscfformat.cpp

Source
/**********************************************************************
Copyright (C) 2004 by Chris Morley for template
Copyright (C) 2009 by David C. Lonie for PWscf

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>
#include <openbabel/mol.h>
#include <openbabel/atom.h>
#include <openbabel/bond.h>
#include <openbabel/obiter.h>
#include <openbabel/elements.h>
#include <openbabel/generic.h>
#include <cstdlib>


#define RYDBERG_TO_KCAL_PER_MOL 313.755026
#define RYDBERG_TO_ELECTRON_VOLT 13.60569193
#define BOHR_TO_ANGSTROM .529177
#define EV_TO_KCAL_PER_MOL 23.060538

using namespace std;
namespace OpenBabel {
  class PWscfFormat : public OBMoleculeFormat
  {
  public:

    PWscfFormat()
    {
      OBConversion::RegisterFormat("pwscf",this);
    }

    const char* Description() override
    {
      return "PWscf format\n"
             "The format used by PWscf, part of Quantum Espresso.\n\n";
    }

    const char* SpecificationURL() override {
      return "https://www.quantum-espresso.org/";
    }

    /* Flags() can return be any of the following combined by |
       or be omitted if none apply
       NOTREADABLE  READONEONLY  NOTWRITABLE  WRITEONEONLY  DEFAULTFORMAT
       READBINARY  WRITEBINARY  READXML  ZEROATOMSOK */
    unsigned int Flags() override
    {
      return READONEONLY | NOTWRITABLE;
    }

    int SkipObjects(int /*n*/, OBConversion* /*pConv*/) override
    {
      return 0;
    }

    ////////////////////////////////////////////////////
    /// Declarations for the "API" interface functions. Definitions are below
    bool ReadMolecule(OBBase* pOb, OBConversion* pConv) override;
    //    virtual bool WriteMolecule(OBBase* pOb, OBConversion* pConv);

  private:
    /* Add declarations for any local function or member variables used.
       Generally only a single instance of a format class is used. Keep this in
       mind if you employ member variables. */
  };
  ////////////////////////////////////////////////////

  //Make an instance of the format class
  PWscfFormat thePWscfFormat;

  /////////////////////////////////////////////////////////////////

  bool PWscfFormat::ReadMolecule(OBBase* pOb, OBConversion* pConv)
  {
    OBMol* pmol = pOb->CastAndClear<OBMol>();
    if (pmol == nullptr)
      return false;

    //Define some references so we can use the old parameter names
    istream &ifs = *pConv->GetInStream();

    char buffer[BUFF_SIZE], tag[BUFF_SIZE];
    double x,y,z;
    double alat = 1.0;
    vector<string> vs;
    matrix3x3 ortho;
    int atomicNum;
    OBUnitCell *cell = new OBUnitCell();
    bool hasEnthalpy=false;
    double enthalpy, pv;

    pmol->BeginModify();

    while (ifs.getline(buffer,BUFF_SIZE)) {

      // Older version of pwscf may use this for alat
      if (strstr(buffer, "lattice parameter (a_0)")) {
        tokenize(vs, buffer);
        alat = atof(vs.at(4).c_str());
      }

      // Newer versions will use this for alat instead
      if (strstr(buffer, "lattice parameter (alat)")) {
        tokenize(vs, buffer);
        alat = atof(vs.at(4).c_str());
      }

      // Unit cell info
      // Newer versions will also say "CELL_PARAMETERS" to complain that no
      // units were specified
      if (strstr(buffer, "CELL_PARAMETERS") &&
          !strstr(buffer, "no units specified in CELL_PARAMETERS card")) {
        // Discover units
        double conv = 1.0;
        tokenize(vs, buffer);

        if (vs.size() > 1) {
          if (strstr(vs[1].c_str(), "alat")) {
            conv = alat * BOHR_TO_ANGSTROM;
          }
          else if (strstr(vs[1].c_str(), "bohr")) {
            conv = BOHR_TO_ANGSTROM;
          }
        }
        // Add others if needed

        double v11, v12, v13,
          v21, v22, v23,
          v31, v32, v33;

        ifs.getline(buffer,BUFF_SIZE); // v1
        tokenize(vs, buffer);
        v11 = atof(vs.at(0).c_str()) * conv;
        v12 = atof(vs.at(1).c_str()) * conv;
        v13 = atof(vs.at(2).c_str()) * conv;

        ifs.getline(buffer,BUFF_SIZE); // v2
        tokenize(vs, buffer);
        v21 = atof(vs.at(0).c_str()) * conv;
        v22 = atof(vs.at(1).c_str()) * conv;
        v23 = atof(vs.at(2).c_str()) * conv;

        ifs.getline(buffer,BUFF_SIZE); // v3
        tokenize(vs, buffer);
        v31 = atof(vs.at(0).c_str()) * conv;
        v32 = atof(vs.at(1).c_str()) * conv;
        v33 = atof(vs.at(2).c_str()) * conv;

        // Build unit cell
        cell->SetData(vector3(v11,v12,v13),
                      vector3(v21,v22,v23),
                      vector3(v31,v32,v33));
      }

      // Unit cell info (for non-variable cell calcs)
      if (strstr(buffer, "crystal axes: (cart. coord. in units of a_0)") ||
          strstr(buffer, "crystal axes: (cart. coord. in units of alat)")) {
        double conv = alat * BOHR_TO_ANGSTROM;
        double v11, v12, v13,
          v21, v22, v23,
          v31, v32, v33;

        ifs.getline(buffer,BUFF_SIZE); // v1
        tokenize(vs, buffer);
        v11 = atof(vs.at(3).c_str()) * conv;
        v12 = atof(vs.at(4).c_str()) * conv;
        v13 = atof(vs.at(5).c_str()) * conv;

        ifs.getline(buffer,BUFF_SIZE); // v2
        tokenize(vs, buffer);
        v21 = atof(vs.at(3).c_str()) * conv;
        v22 = atof(vs.at(4).c_str()) * conv;
        v23 = atof(vs.at(5).c_str()) * conv;

        ifs.getline(buffer,BUFF_SIZE); // v3
        tokenize(vs, buffer);
        v31 = atof(vs.at(3).c_str()) * conv;
        v32 = atof(vs.at(4).c_str()) * conv;
        v33 = atof(vs.at(5).c_str()) * conv;

        // Build unit cell
        cell->SetData(vector3(v11,v12,v13),
                      vector3(v21,v22,v23),
                      vector3(v31,v32,v33));
      }

      // Atoms info
      if (strstr(buffer, "ATOMIC_POSITIONS")) {
        // Clear old atoms from pmol
        vector<OBAtom*> toDelete;
        FOR_ATOMS_OF_MOL(a, *pmol)
          toDelete.push_back(&*a);
        for (size_t i = 0; i < toDelete.size(); i++)
          pmol->DeleteAtom(toDelete.at(i));


        // Discover units
        matrix3x3 conv (1);
        tokenize(vs, buffer);

        if (vs.size() > 1) {
          if (strstr(vs[1].c_str(), "alat")) {
            conv *= (alat * BOHR_TO_ANGSTROM);
          }
          else if (strstr(vs[1].c_str(), "crystal")) {
            // Set to the zero matrix and test below.
            conv = matrix3x3 (0.0);
          }
        }
        // Add others if needed

        // Load new atoms from molecule
        ifs.getline(buffer,BUFF_SIZE); // First entry
        tokenize(vs, buffer);
        int size = vs.size();
        while (size == 4) {
          atomicNum = OBElements::GetAtomicNum(vs[0].c_str());
          x = atof((char*)vs[1].c_str());
          y = atof((char*)vs[2].c_str());
          z = atof((char*)vs[3].c_str());
          // Add atom
          OBAtom *atom = pmol->NewAtom();
          atom->SetAtomicNum(atomicNum);
          vector3 coords (x,y,z);
          if (conv.determinant() == 0.0) { // Fractional coords
            atom->SetVector(cell->FractionalToCartesian(coords));
          }
          else {
            atom->SetVector(conv * coords);
          }

          // Reset vars
          ifs.getline(buffer,BUFF_SIZE); // First entry
          tokenize(vs, buffer);
          size = vs.size();
        }
      }

      // Free energy
      if (strstr(buffer, "Final energy =")) {
        tokenize(vs, buffer);
        if (vs.size() > 3) {
          pmol->SetEnergy(atof(vs[3].c_str()) * RYDBERG_TO_KCAL_PER_MOL);
        }
      }

      // H - PV = U energy
      if (strstr(buffer, "!    total energy              =")) {
        tokenize(vs, buffer);
        if (vs.size() > 4) {
          pmol->SetEnergy(atof(vs[4].c_str()) * RYDBERG_TO_KCAL_PER_MOL);
        }
      }

      // Enthalphy
      if (strstr(buffer, "Final enthalpy =")) {
        tokenize(vs, buffer);

        hasEnthalpy = true;
        enthalpy = atof(vs.at(3).c_str()) * RYDBERG_TO_KCAL_PER_MOL;
        pv = enthalpy - pmol->GetEnergy();
      }
    }

    // set final unit cell
    pmol->SetData(cell);

    // Set enthalpy
    if (hasEnthalpy) {
      OBPairData *enthalpyPD = new OBPairData();
      OBPairData *enthalpyPD_pv = new OBPairData();
      OBPairData *enthalpyPD_eV = new OBPairData();
      OBPairData *enthalpyPD_pv_eV = new OBPairData();
      enthalpyPD->SetAttribute("Enthalpy (kcal/mol)");
      enthalpyPD_pv->SetAttribute("Enthalpy PV term (kcal/mol)");
      enthalpyPD_eV->SetAttribute("Enthalpy (eV)");
      enthalpyPD_pv_eV->SetAttribute("Enthalpy PV term (eV)");
      double en_kcal_per_mole = enthalpy;
      double pv_kcal_per_mole = pv;
      double en_eV = enthalpy / EV_TO_KCAL_PER_MOL;
      double pv_eV = pv / EV_TO_KCAL_PER_MOL;
      snprintf(tag, BUFF_SIZE, "%f", en_kcal_per_mole);
      enthalpyPD->SetValue(tag);
      snprintf(tag, BUFF_SIZE, "%f", pv_kcal_per_mole);
      enthalpyPD_pv->SetValue(tag);
      snprintf(tag, BUFF_SIZE, "%f", en_eV);
      enthalpyPD_eV->SetValue(tag);
      snprintf(tag, BUFF_SIZE, "%f", pv_eV);
      enthalpyPD_pv_eV->SetValue(tag);
      pmol->SetData(enthalpyPD);
      pmol->SetData(enthalpyPD_pv);
      pmol->SetData(enthalpyPD_eV);
      pmol->SetData(enthalpyPD_pv_eV);
    }

    pmol->EndModify();

    return true;
  }

} //namespace OpenBabel

Coverage Report

Created: 2026-08-31 06:40

Line	Count	Source
1		/**********************************************************************
2		Copyright (C) 2004 by Chris Morley for template
3		Copyright (C) 2009 by David C. Lonie for PWscf
4
5		This program is free software; you can redistribute it and/or modify
6		it under the terms of the GNU General Public License as published by
7		the Free Software Foundation version 2 of the License.
8
9		This program is distributed in the hope that it will be useful,
10		but WITHOUT ANY WARRANTY; without even the implied warranty of
11		MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12		GNU General Public License for more details.
13		***********************************************************************/
14
15		#include <openbabel/babelconfig.h>
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	0	#define RYDBERG_TO_KCAL_PER_MOL 313.755026
27		#define RYDBERG_TO_ELECTRON_VOLT 13.60569193
28	0	#define BOHR_TO_ANGSTROM .529177
29	0	#define EV_TO_KCAL_PER_MOL 23.060538
30
31		using namespace std;
32		namespace OpenBabel {
33		class PWscfFormat : public OBMoleculeFormat
34		{
35		public:
36
37		PWscfFormat()
38	12	{
39	12	OBConversion::RegisterFormat("pwscf",this);
40	12	}
41
42		const char* Description() override
43	0	{
44	0	return "PWscf format\n"
45	0	"The format used by PWscf, part of Quantum Espresso.\n\n";
46	0	}
47
48	0	const char* SpecificationURL() override {
49	0	return "https://www.quantum-espresso.org/";
50	0	}
51
52		/* Flags() can return be any of the following combined by \|
53		or be omitted if none apply
54		NOTREADABLE READONEONLY NOTWRITABLE WRITEONEONLY DEFAULTFORMAT
55		READBINARY WRITEBINARY READXML ZEROATOMSOK */
56		unsigned int Flags() override
57	13	{
58	13	return READONEONLY \| NOTWRITABLE;
59	13	}
60
61		int SkipObjects(int /n/, OBConversion* /pConv/) override
62	0	{
63	0	return 0;
64	0	}
65
66		////////////////////////////////////////////////////
67		/// Declarations for the "API" interface functions. Definitions are below
68		bool ReadMolecule(OBBase* pOb, OBConversion* pConv) override;
69		// virtual bool WriteMolecule(OBBase* pOb, OBConversion* pConv);
70
71		private:
72		/* Add declarations for any local function or member variables used.
73		Generally only a single instance of a format class is used. Keep this in
74		mind if you employ member variables. */
75		};
76		////////////////////////////////////////////////////
77
78		//Make an instance of the format class
79		PWscfFormat thePWscfFormat;
80
81		/////////////////////////////////////////////////////////////////
82
83		bool PWscfFormat::ReadMolecule(OBBase* pOb, OBConversion* pConv)
84	0	{
85	0	OBMol* pmol = pOb->CastAndClear<OBMol>();
86	0	if (pmol == nullptr)
87	0	return false;
88
89		//Define some references so we can use the old parameter names
90	0	istream &ifs = *pConv->GetInStream();
91
92	0	char buffer[BUFF_SIZE], tag[BUFF_SIZE];
93	0	double x,y,z;
94	0	double alat = 1.0;
95	0	vector<string> vs;
96	0	matrix3x3 ortho;
97	0	int atomicNum;
98	0	OBUnitCell *cell = new OBUnitCell();
99	0	bool hasEnthalpy=false;
100	0	double enthalpy, pv;
101
102	0	pmol->BeginModify();
103
104	0	while (ifs.getline(buffer,BUFF_SIZE)) {
105
106		// Older version of pwscf may use this for alat
107	0	if (strstr(buffer, "lattice parameter (a_0)")) {
108	0	tokenize(vs, buffer);
109	0	alat = atof(vs.at(4).c_str());
110	0	}
111
112		// Newer versions will use this for alat instead
113	0	if (strstr(buffer, "lattice parameter (alat)")) {
114	0	tokenize(vs, buffer);
115	0	alat = atof(vs.at(4).c_str());
116	0	}
117
118		// Unit cell info
119		// Newer versions will also say "CELL_PARAMETERS" to complain that no
120		// units were specified
121	0	if (strstr(buffer, "CELL_PARAMETERS") &&
122	0	!strstr(buffer, "no units specified in CELL_PARAMETERS card")) {
123		// Discover units
124	0	double conv = 1.0;
125	0	tokenize(vs, buffer);
126
127	0	if (vs.size() > 1) {
128	0	if (strstr(vs[1].c_str(), "alat")) {
129	0	conv = alat * BOHR_TO_ANGSTROM;
130	0	}
131	0	else if (strstr(vs[1].c_str(), "bohr")) {
132	0	conv = BOHR_TO_ANGSTROM;
133	0	}
134	0	}
135		// Add others if needed
136
137	0	double v11, v12, v13,
138	0	v21, v22, v23,
139	0	v31, v32, v33;
140
141	0	ifs.getline(buffer,BUFF_SIZE); // v1
142	0	tokenize(vs, buffer);
143	0	v11 = atof(vs.at(0).c_str()) * conv;
144	0	v12 = atof(vs.at(1).c_str()) * conv;
145	0	v13 = atof(vs.at(2).c_str()) * conv;
146
147	0	ifs.getline(buffer,BUFF_SIZE); // v2
148	0	tokenize(vs, buffer);
149	0	v21 = atof(vs.at(0).c_str()) * conv;
150	0	v22 = atof(vs.at(1).c_str()) * conv;
151	0	v23 = atof(vs.at(2).c_str()) * conv;
152
153	0	ifs.getline(buffer,BUFF_SIZE); // v3
154	0	tokenize(vs, buffer);
155	0	v31 = atof(vs.at(0).c_str()) * conv;
156	0	v32 = atof(vs.at(1).c_str()) * conv;
157	0	v33 = atof(vs.at(2).c_str()) * conv;
158
159		// Build unit cell
160	0	cell->SetData(vector3(v11,v12,v13),
161	0	vector3(v21,v22,v23),
162	0	vector3(v31,v32,v33));
163	0	}
164
165		// Unit cell info (for non-variable cell calcs)
166	0	if (strstr(buffer, "crystal axes: (cart. coord. in units of a_0)") \|\|
167	0	strstr(buffer, "crystal axes: (cart. coord. in units of alat)")) {
168	0	double conv = alat * BOHR_TO_ANGSTROM;
169	0	double v11, v12, v13,
170	0	v21, v22, v23,
171	0	v31, v32, v33;
172
173	0	ifs.getline(buffer,BUFF_SIZE); // v1
174	0	tokenize(vs, buffer);
175	0	v11 = atof(vs.at(3).c_str()) * conv;
176	0	v12 = atof(vs.at(4).c_str()) * conv;
177	0	v13 = atof(vs.at(5).c_str()) * conv;
178
179	0	ifs.getline(buffer,BUFF_SIZE); // v2
180	0	tokenize(vs, buffer);
181	0	v21 = atof(vs.at(3).c_str()) * conv;
182	0	v22 = atof(vs.at(4).c_str()) * conv;
183	0	v23 = atof(vs.at(5).c_str()) * conv;
184
185	0	ifs.getline(buffer,BUFF_SIZE); // v3
186	0	tokenize(vs, buffer);
187	0	v31 = atof(vs.at(3).c_str()) * conv;
188	0	v32 = atof(vs.at(4).c_str()) * conv;
189	0	v33 = atof(vs.at(5).c_str()) * conv;
190
191		// Build unit cell
192	0	cell->SetData(vector3(v11,v12,v13),
193	0	vector3(v21,v22,v23),
194	0	vector3(v31,v32,v33));
195	0	}
196
197		// Atoms info
198	0	if (strstr(buffer, "ATOMIC_POSITIONS")) {
199		// Clear old atoms from pmol
200	0	vector<OBAtom*> toDelete;
201	0	FOR_ATOMS_OF_MOL(a, *pmol)
202	0	toDelete.push_back(&*a);
203	0	for (size_t i = 0; i < toDelete.size(); i++)
204	0	pmol->DeleteAtom(toDelete.at(i));
205
206
207		// Discover units
208	0	matrix3x3 conv (1);
209	0	tokenize(vs, buffer);
210
211	0	if (vs.size() > 1) {
212	0	if (strstr(vs[1].c_str(), "alat")) {
213	0	conv = (alat BOHR_TO_ANGSTROM);
214	0	}
215	0	else if (strstr(vs[1].c_str(), "crystal")) {
216		// Set to the zero matrix and test below.
217	0	conv = matrix3x3 (0.0);
218	0	}
219	0	}
220		// Add others if needed
221
222		// Load new atoms from molecule
223	0	ifs.getline(buffer,BUFF_SIZE); // First entry
224	0	tokenize(vs, buffer);
225	0	int size = vs.size();
226	0	while (size == 4) {
227	0	atomicNum = OBElements::GetAtomicNum(vs[0].c_str());
228	0	x = atof((char*)vs[1].c_str());
229	0	y = atof((char*)vs[2].c_str());
230	0	z = atof((char*)vs[3].c_str());
231		// Add atom
232	0	OBAtom *atom = pmol->NewAtom();
233	0	atom->SetAtomicNum(atomicNum);
234	0	vector3 coords (x,y,z);
235	0	if (conv.determinant() == 0.0) { // Fractional coords
236	0	atom->SetVector(cell->FractionalToCartesian(coords));
237	0	}
238	0	else {
239	0	atom->SetVector(conv * coords);
240	0	}
241
242		// Reset vars
243	0	ifs.getline(buffer,BUFF_SIZE); // First entry
244	0	tokenize(vs, buffer);
245	0	size = vs.size();
246	0	}
247	0	}
248
249		// Free energy
250	0	if (strstr(buffer, "Final energy =")) {
251	0	tokenize(vs, buffer);
252	0	if (vs.size() > 3) {
253	0	pmol->SetEnergy(atof(vs[3].c_str()) * RYDBERG_TO_KCAL_PER_MOL);
254	0	}
255	0	}
256
257		// H - PV = U energy
258	0	if (strstr(buffer, "! total energy =")) {
259	0	tokenize(vs, buffer);
260	0	if (vs.size() > 4) {
261	0	pmol->SetEnergy(atof(vs[4].c_str()) * RYDBERG_TO_KCAL_PER_MOL);
262	0	}
263	0	}
264
265		// Enthalphy
266	0	if (strstr(buffer, "Final enthalpy =")) {
267	0	tokenize(vs, buffer);
268
269	0	hasEnthalpy = true;
270	0	enthalpy = atof(vs.at(3).c_str()) * RYDBERG_TO_KCAL_PER_MOL;
271	0	pv = enthalpy - pmol->GetEnergy();
272	0	}
273	0	}
274
275		// set final unit cell
276	0	pmol->SetData(cell);
277
278		// Set enthalpy
279	0	if (hasEnthalpy) {
280	0	OBPairData *enthalpyPD = new OBPairData();
281	0	OBPairData *enthalpyPD_pv = new OBPairData();
282	0	OBPairData *enthalpyPD_eV = new OBPairData();
283	0	OBPairData *enthalpyPD_pv_eV = new OBPairData();
284	0	enthalpyPD->SetAttribute("Enthalpy (kcal/mol)");
285	0	enthalpyPD_pv->SetAttribute("Enthalpy PV term (kcal/mol)");
286	0	enthalpyPD_eV->SetAttribute("Enthalpy (eV)");
287	0	enthalpyPD_pv_eV->SetAttribute("Enthalpy PV term (eV)");
288	0	double en_kcal_per_mole = enthalpy;
289	0	double pv_kcal_per_mole = pv;
290	0	double en_eV = enthalpy / EV_TO_KCAL_PER_MOL;
291	0	double pv_eV = pv / EV_TO_KCAL_PER_MOL;
292	0	snprintf(tag, BUFF_SIZE, "%f", en_kcal_per_mole);
293	0	enthalpyPD->SetValue(tag);
294	0	snprintf(tag, BUFF_SIZE, "%f", pv_kcal_per_mole);
295	0	enthalpyPD_pv->SetValue(tag);
296	0	snprintf(tag, BUFF_SIZE, "%f", en_eV);
297	0	enthalpyPD_eV->SetValue(tag);
298	0	snprintf(tag, BUFF_SIZE, "%f", pv_eV);
299	0	enthalpyPD_pv_eV->SetValue(tag);
300	0	pmol->SetData(enthalpyPD);
301	0	pmol->SetData(enthalpyPD_pv);
302	0	pmol->SetData(enthalpyPD_eV);
303	0	pmol->SetData(enthalpyPD_pv_eV);
304	0	}
305
306	0	pmol->EndModify();
307
308	0	return true;
309	0	}
310
311		} //namespace OpenBabel