sqs2poscar.cpp

//
// sqs2poscar.cpp
//
// Usage: sqs2poscar FileName
//
// The input file must be a standard "bestsqs.out" file
// generated by mcsqs (ATAT 3.0). The output file will
// be named as "FileName-POSCAR".
//
// Written by Changning Niu (Oct. 10, 2013)
//

#include <iostream>
#include <fstream>
#include <sstream>
#include <string>
#include <cmath>
#include <iomanip>

using namespace std;

/// read_bestsqs:
///     return true if successful; false otherwise.
bool read_bestsqs (char *file, double dArrVec1[3][3], double dArrVec2[3][3],
                   double dArrAtom[500][3], string strArrAtom[500], int& nAtom,
                   int& nElem, int nArrElem[10], string strArrElem[10]) {
    istringstream iss;
    ifstream is;
    string line, elem;
    double x, y, z;
    int num = 0;
    bool signal;

    is.open(file);
    if (!is) {
        return false;
    }
    for (int i = 0; i < 3; i++) {   // line 1-3: basis coordinate system
        getline(is, line);
        iss.str(line);
        if (iss >> x >> y >> z) {
            iss.clear();
            dArrVec1[i][0] = x;
            dArrVec1[i][1] = y;
            dArrVec1[i][2] = z;
        }
        else
            return false;
    }
    for (int i = 0; i < 3; i++) {   // line 4-6: lattice vectors scaled by line 1-3
        getline(is, line);
        iss.str(line);
        if (iss >> x >> y >> z) {
            iss.clear();
            dArrVec2[i][0] = x;
            dArrVec2[i][1] = y;
            dArrVec2[i][2] = z;
        }
        else
            return false;
    }
    while (getline(is, line)) {     // line 7-: atomic positions scaled by line 1-3
        iss.str(line);
        if (iss >> x >> y >> z >> elem) {
            iss.clear();
            dArrAtom[num][0] = x;
            dArrAtom[num][1] = y;
            dArrAtom[num][2] = z;
            strArrAtom[num] = elem;
            num++;
        }
        else
            return false;
    }
    nAtom = num;
    /// count number of elements
    for (int i = 0; i < 10; i++) {
        nArrElem[i] = 0;
    }
    strArrElem[0] = strArrAtom[0];      // this is the first atom
    nElem = 1; nArrElem[0] = 1;
    for (int i = 1; i < nAtom; i++) {
        signal = false;
        for (int j = 0; j < nElem; j++) {
            if (strArrAtom[i] == strArrElem[j]) {
                nArrElem[j]++;
                signal = true;
                break;
            }
        }
        if (signal == false) {
            nArrElem[nElem] = 1;
            strArrElem[nElem] = strArrAtom[i];
            nElem++;
        }
    }
    return true;
}

/// calculate inverse matrix
bool inverse_matrix (double m1[3][3], double m2[3][3]) {
    double det;
    det = m1[0][0] * m1[1][1] * m1[2][2]
        - m1[0][0] * m1[1][2] * m1[2][1]
        - m1[0][1] * m1[1][0] * m1[2][2]
        + m1[0][1] * m1[1][2] * m1[2][0]
        + m1[0][2] * m1[1][0] * m1[2][1]
        - m1[0][2] * m1[1][1] * m1[2][0];
    if (det == 0)
        return false;
    m2[0][0] = m1[1][1] * m1[2][2] - m1[1][2] * m1[2][1];
    m2[0][1] = m1[0][2] * m1[2][1] - m1[0][1] * m1[2][2];
    m2[0][2] = m1[0][1] * m1[1][2] - m1[0][2] * m1[1][1];
    m2[1][0] = m1[1][2] * m1[2][0] - m1[1][0] * m1[2][2];
    m2[1][1] = m1[0][0] * m1[2][2] - m1[0][2] * m1[2][0];
    m2[1][2] = m1[0][2] * m1[1][0] - m1[0][0] * m1[1][2];
    m2[2][0] = m1[1][0] * m1[2][1] - m1[1][1] * m1[2][0];
    m2[2][1] = m1[0][1] * m1[2][0] - m1[0][0] * m1[2][1];
    m2[2][2] = m1[0][0] * m1[1][1] - m1[0][1] * m1[1][0];
    for (int i = 0; i < 3; i++) {
        for (int j = 0; j < 3; j++) {
            m2[i][j] /= det;
        }
    }
    return true;
}

int main (int argc, char* argv[]) {
    /// Variables
    char *file;                 // file name given by user
    double dArrVec1[3][3];      // basis vectors in bestsqs
    double dArrVec2[3][3];      // lattice vectors in bestsqs
    double dArrAtom[500][3];    // atomic coordinates in bestsqs
    string strArrAtom[500];     // atomic species in bestsqs
    int nAtom;                  // # of atoms
    int nElem;                  // # of elements
    int nArrElem[10];           // # of atoms of each element
    string strArrElem[10];      // name of elements
    double dArrLatVec[3][3];    // lattice vectors in POSCAR
    double dArrLatVecInv[3][3]; // inverse matrix of vectors
    double dArrAtom2[500][3];   // atomic coordinates in POSCAR
    string strTemp;
    int count;
    double x, y, z;

    if (argc != 2) {
        cout << "Usage: koala_sqs2poscar FileName\n";
        exit(1);
    }
    file = argv[1];
    if (!read_bestsqs(file, dArrVec1, dArrVec2, dArrAtom, strArrAtom, nAtom,
                      nElem, nArrElem, strArrElem)) {
        cout << "Errors occurred when reading " << file << "!\n";
        exit(1);
    }
    /// calculate lattice vectors in cartesian
    for (int i = 0; i < 3; i++) {
        dArrLatVec[i][0] = dArrVec2[i][0] * dArrVec1[0][0]
                         + dArrVec2[i][1] * dArrVec1[1][0]
                         + dArrVec2[i][2] * dArrVec1[2][0];
        dArrLatVec[i][1] = dArrVec2[i][0] * dArrVec1[0][1]
                         + dArrVec2[i][1] * dArrVec1[1][1]
                         + dArrVec2[i][2] * dArrVec1[2][1];
        dArrLatVec[i][2] = dArrVec2[i][0] * dArrVec1[0][2]
                         + dArrVec2[i][1] * dArrVec1[1][2]
                         + dArrVec2[i][2] * dArrVec1[2][2];
    }
    /// Calculate atomic positions in cartesian
    for (int i = 0; i < nAtom; i++) {
        x = dArrAtom[i][0];
        y = dArrAtom[i][1];
        z = dArrAtom[i][2];
        dArrAtom[i][0] = x * dArrVec1[0][0]
                       + y * dArrVec1[1][0]
                       + z * dArrVec1[2][0];
        dArrAtom[i][1] = x * dArrVec1[0][1]
                       + y * dArrVec1[1][1]
                       + z * dArrVec1[2][1];
        dArrAtom[i][2] = x * dArrVec1[0][2]
                       + y * dArrVec1[1][2]
                       + z * dArrVec1[2][2];
    }
    /// Sort atomic coordinates based on its symbol
    count = 0;
    for (int i = 0; i < nElem; i++) {
        for (int j = 0; j < nAtom; j++) {
            if (strArrAtom[j] == strArrElem[i]) {
                for (int k = 0; k < 3; k++) {
                    dArrAtom2[count][k] = dArrAtom[j][k];
                }
                count++;
            }
        }
    }
    if (count != nAtom) {
        cout << "Error occurred when sorting!\n";
        exit(1);
    }
    /// Calculate fractional coordinates
    if (inverse_matrix(dArrLatVec, dArrLatVecInv)) {
        for (int i = 0; i < nAtom; i++) {
            x = dArrAtom2[i][0];
            y = dArrAtom2[i][1];
            z = dArrAtom2[i][2];
            dArrAtom2[i][0] = x * dArrLatVecInv[0][0]
                            + y * dArrLatVecInv[1][0]
                            + z * dArrLatVecInv[2][0];
            dArrAtom2[i][1] = x * dArrLatVecInv[0][1]
                            + y * dArrLatVecInv[1][1]
                            + z * dArrLatVecInv[2][1];
            dArrAtom2[i][2] = x * dArrLatVecInv[0][2]
                            + y * dArrLatVecInv[1][2]
                            + z * dArrLatVecInv[2][2];
        }
    }
    else {
        cout << "Errors occurred when converting to direct.\n";
        exit (1);
    }
    /// Output
    ofstream os;
    os.open(strcat(file, "-POSCAR"));
    if (!os) {
        cout << "Cannot write to file!\n";
        exit(1);
    }
    os << fixed;
    os << "POSCAR\n";   // line 1: head
    os << "xxx\n";      // line 2: universal scaler
    for (int i = 0; i < 3; i++) {   // line 3-5: lattice vectors
        for (int j = 0; j < 3; j++) {
            os << setw(12) << setprecision(8) << dArrLatVec[i][j];
        }
        os << endl;
    }
    for (int i = 0; i < nElem; i++) {   // line 6: element names
        os << setw(4) << strArrElem[i];
    }
    os << endl;
    for (int i = 0; i < nElem; i++) {   // line 7: # of atoms of each elemt
        os << setw(4) << nArrElem[i];
    }
    os << "\nDirect\n";
    for (int i = 0; i < nAtom; i++) {
        for (int j = 0; j < 3; j++) {
            os << setw(12) << setprecision(8) << dArrAtom2[i][j];
        }
        os << endl;
    }
    os.close();
    exit(0);
}