Main Page   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Namespace Members   Compound Members   File Members   Related Pages  

orcaplugin.C

Go to the documentation of this file.
00001 /***************************************************************************
00002  * RCS INFORMATION:
00003  *
00004  *      $RCSfile: orcaplugin.C,v $
00005  *      $Author: johns $       $Locker:  $             $State: Exp $
00006  *      $Revision: 1.5 $       $Date: 2019/01/29 22:39:44 $
00007  *
00008  ***************************************************************************
00009  * DESCRIPTION:
00010  *   Orca VMD plugin
00011  *   Authors: Maximilian Scheurer, Marcelo Melo, April 2017
00012  *   Contributions from: Michael F. Herbst
00013  *   Inspired from gamessplugin.c
00014  *
00015  ***************************************************************************/
00016 
00017 #include <stdlib.h>
00018 #include <stdarg.h>
00019 #include <string.h>
00020 #include <ctype.h>
00021 #include <errno.h>
00022 #include <time.h>
00023 #include <iostream>
00024 #include <sstream>
00025 #include <algorithm>
00026 #include <string>
00027 #include <vector>
00028 #include <cctype>
00029 #include "qmplugin.h"
00030 #include "unit_conversion.h"
00031 #include "periodic_table.h"
00032 // #include "Matrix.h" // trivial matrix routines
00033 
00034 using namespace std; 
00035 
00036 class Matrix {
00037 public:
00038   double **value;
00039   int columns, rows;
00040 
00041   Matrix() {
00042     value = NULL;
00043     rows = 0;
00044     columns = 0;
00045   };
00046 
00047   Matrix(std::string input) {
00048     value = parseInput(input);
00049   }
00050 
00051   // template <typename T>
00052   Matrix(std::vector<std::vector<float> > input) {
00053     rows = input.size();
00054     columns = input[0].size();
00055     value = new double *[rows];
00056     for (int i = 0; i < rows; i++) {
00057       value[i] = new double [columns];
00058       if (columns != input[i].size()) {
00059 //        cout << "Matrix rows have different number of columns." << endl;
00060         return;
00061       }
00062     }
00063 
00064     int row = 0, col = 0;
00065     // for (auto r : input) {
00066     //   for (auto v : r) {
00067     for (int lr=0; lr<input.size(); lr++) {
00068       for (int lv=0; lv<input[0].size(); lv++) {
00069         float v=input[lr][lv];
00070         value[row][col] = v;
00071         col++;
00072       }
00073       row++;
00074       col = 0;
00075     }
00076   }
00077 
00078 
00079   ~Matrix() {
00080     for (int i=0; i < rows; i++) {
00081       delete[] value[i];
00082     }
00083     delete[] value;
00084   }
00085 
00086 
00087   double** parseInput(std::string input) {
00088     rows = 0;
00089     columns = 0;
00090     vector<double> tempValues;
00091     string tempString;
00092     bool firstBracket = false;
00093     bool closedBracket = false;
00094     char previousChar = '\0';
00095 
00096 //    for (char &c : input) {
00097     for (int lv=0; lv<input.size(); lv++) {
00098       char c = input[lv]; 
00099 
00100       if (firstBracket && c == '{') {
00101         rows++;
00102       } else if(c == '{') {
00103         firstBracket = true;
00104       } else if(c == ',') {
00105         tempValues.push_back(stringToDouble(tempString));
00106 //        tempString.clear();
00107         tempString = ""; // for very old revs of MSVC
00108         if (closedBracket == false) {
00109           columns++;
00110         }
00111       } else if(c == '}' && closedBracket == false) {
00112         closedBracket = true;
00113         columns++;
00114       } else if (c != ',' && c != '{' && c != '}') {
00115         tempString.append(&c);
00116       } else if (c == '}' && previousChar == '}') {
00117         tempValues.push_back(stringToDouble(tempString));
00118       }
00119       previousChar = c;
00120     }
00121 
00122     double **matrix = new double *[rows];
00123     for (int i = 0; i < rows; i++) {
00124       matrix[i] = new double[columns];
00125     }
00126 
00127 #if 1
00128     // support ancient versions of MSVC
00129     int n;
00130     for (n=0; n<tempValues.size(); n++) {
00131       div_t mod = div(n, columns);
00132       matrix[mod.quot][mod.rem] = tempValues[n];
00133     }
00134 #else
00135     int n = 0;
00136     for (vector<double>::iterator i = tempValues.begin(); i < tempValues.end(); i++) {
00137       div_t mod = div(n, columns);
00138       matrix[mod.quot][mod.rem] = *i;
00139       n++;
00140     }
00141 #endif
00142 
00143     return matrix;
00144   }
00145 
00146 #if 1
00147   static double stringToDouble(std::string& s) {
00148     return atof(s.c_str());
00149   }
00150 #else
00151   static double stringToDouble(std::string& s) {
00152     istringstream i(s);
00153     double x;
00154     if (!(i >> x))
00155         return 0;
00156     return x;
00157   }
00158 #endif
00159 
00160   static void printMatrix(Matrix *matrix) {
00161     if (matrix != NULL) {
00162       for (int y = 0; y < matrix->rows; y++) {
00163         for (int x = 0; x < matrix->columns; x++) {
00164           if (fabs(matrix->value[y][x]) == 0) {
00165             matrix->value[y][x] = fabs(matrix->value[y][x]);
00166           }
00167           cout << matrix->value[y][x] << " ";
00168         }
00169         cout << endl;
00170       }
00171       cout << endl;
00172     } else {
00173       cout << "Matrix empty" << endl << endl;
00174     }
00175   }
00176 
00177 
00178   void printMatrix() {
00179     if (value != NULL) {
00180       for (int y = 0; y < rows; y++) {
00181         for (int x = 0; x < columns; x++) {
00182           if (fabs(value[y][x]) == 0) {
00183             value[y][x] = fabs(value[y][x]);
00184           }
00185           cout << value[y][x] << " ";
00186         }
00187         cout << endl;
00188       }
00189     cout << endl;
00190     } else {
00191       cout << "Matrix empty" << endl << endl;
00192     }
00193   }
00194 
00195 
00196   std::vector< std::vector<float> > toVector() {
00197     std::vector<std::vector<float> > result;
00198     for (size_t i = 0; i < rows; i++) {
00199         std::vector<float> tmpRow;
00200         for (size_t j = 0; j < columns; j++) {
00201             tmpRow.push_back(value[i][j]);
00202         }
00203         result.push_back(tmpRow);
00204     }
00205     return result;
00206   }
00207 
00208 
00209   static std::vector<double>* rowAtIndex(Matrix *input, unsigned int index) {
00210     vector<double> *result = new vector<double>;
00211     for (int i = 0; i < input->columns; i++) {
00212         result->push_back(input->value[index][i]);
00213     }
00214     return result;
00215   }
00216 
00217 
00218   static std::vector<double>* columnAtIndex(Matrix *input, unsigned int index) {
00219     vector<double> *result = new vector<double>;
00220     for (int i = 0; i < input->rows; i++) {
00221         result->push_back(input->value[i][index]);
00222     }
00223     return result;
00224   }
00225 
00226 
00227   static double dotProduct(std::vector<double> *firstVector, std::vector<double> *secondVector) {
00228     double result = 0.0;
00229     if (firstVector == NULL || secondVector == NULL){
00230       cout << "Nullpointer Exception. \n";
00231     } else if (firstVector->size() == secondVector->size()) {
00232       for (int i = 0; i < firstVector->size(); i++) {
00233         result += firstVector->at(i) * secondVector->at(i);
00234       }
00235     }
00236     return result;
00237   }
00238 
00239 
00240   static Matrix* multiply(Matrix *firstMatrix, Matrix *secondMatrix) {
00241     if (firstMatrix->columns == secondMatrix->rows) {
00242       double **matrix = 0;
00243       matrix = new double *[firstMatrix->rows];
00244       for (int i = 0; i < firstMatrix->rows; i++) {
00245         matrix[i] = new double [secondMatrix->columns];
00246       }
00247 
00248       for (int y = 0; y < firstMatrix->rows; y++) {
00249         for (int x = 0; x < secondMatrix->columns; x++) {
00250           std::vector<double, std::allocator<double> > *r = rowAtIndex(firstMatrix, y);
00251           std::vector<double, std::allocator<double> > *c = columnAtIndex(secondMatrix, x);
00252           matrix[y][x] = dotProduct(r, c);
00253           delete r;
00254           delete c;
00255         }
00256       }
00257       Matrix *result = new Matrix();
00258       result->rows = firstMatrix->rows;
00259       result->columns = secondMatrix->columns;
00260       result->value = matrix;
00261       return result;
00262     } else {
00263       return NULL;
00264     }
00265   }
00266 };
00267 
00268 // end of Matrix class
00269 
00270 
00271 
00272 typedef std::vector<std::vector<std::vector<float> > > MoCoeff;
00273 
00274 #define DEBUGGING_ORCA 1
00275 #ifdef DEBUGGING_ORCA
00276 #define PRINTERR fprintf(stderr, "\n In file %s, line %d: \n %s \n \n", \
00277                            __FILE__, __LINE__, strerror(errno))
00278 #else
00279 #define PRINTERR (void)(0)
00280 #endif
00281 
00282 #define ANGSTROM 0
00283 #define BOHR     1
00284 
00285 /*
00286  * Error reporting macro for the multiple fgets calls in
00287  * the code
00288  */
00289 #define GET_LINE(x,y) if (!fgets(x, sizeof(x), y)) return FALSE
00290 
00291 #define NOTFOUND 0
00292 #define FOUND    1
00293 #define STOPPED  2
00294 
00295 #define ORCA4 4
00296 #define ORCA3 3
00297 
00298 /* ######################################################## */
00299 /*                 orca specific struct                     */
00300 /*
00301   might be extended in the future...
00302  */
00303 /* ######################################################## */
00304 
00305 typedef struct {
00306   int version; /* version = 0 : not supported, do not read the file!
00307                 * version = 1 : version 3.0.0, 3.0.1, 3.0.3
00308                 * version = 2 : version 4.0.0
00309                 */
00310 
00311   int digits[3]; /* The three digits of the orca version number
00312                   * e.g. 4.0.0 (i.e. digits[0].digits[1].digits[2])
00313                   * */
00314 } orcadata;
00315 
00316 typedef std::vector<std::vector<float> > CoeffRowBlock;
00317 
00318 // Matrix for conversion of pure coefficients to cartesian coefficients, d- and f-shell
00319 Matrix *convD = new Matrix("{{ 0, 1, 0, 0, 0},"
00320                             "{ 0, 0, 1, 0, 0},"
00321                             "{ 0, 0, 0, 0, 1},"
00322                             "{-1, 0, 0, 1, 0},"
00323                             "{-1, 0, 0,-1, 0},"
00324                             "{ 2, 0, 0, 0, 0}}");
00325 
00326 Matrix *convF = new Matrix("{{ 0, 0, 0, 0, 1, 0, 0},"
00327                             "{ 0, 0,-1, 0, 0, 0, 3},"
00328                             "{-3, 0, 0, 1, 0, 0, 0},"
00329                             "{ 0,-1, 0, 0, 0,-3, 0},"
00330                             "{-3, 0, 0,-1, 0, 0, 0},"
00331                             "{ 0, 4, 0, 0, 0, 0, 0},"
00332                             "{ 0, 0, 4, 0, 0, 0, 0},"
00333                             "{ 0,-1, 0, 0, 0, 1, 0},"
00334                             "{ 0, 0,-1, 0, 0, 0,-1},"
00335                             "{ 2, 0, 0, 0, 0, 0, 0}}");
00336 
00337 
00338 /* ######################################################## */
00339 /*                    static functions                      */
00340 /* ######################################################## */
00341 
00342 // Checks if loaded file is really an Orca output file
00343 static int have_orca(qmdata_t *data, orcadata* orca);
00344 
00345 static int read_orca_structure(void *mydata, int *optflags, molfile_atom_t *atoms);
00346 
00347 // Freeing memory
00348 static void close_orca_read(void *mydata);
00349 
00350 // Function for reading timestep independent information: Main Parser
00351 static int parse_static_data(qmdata_t *, int *);
00352 
00353 
00354 // collect information about the jobtype
00355 static int get_job_info(qmdata_t *data);
00356 
00357 // analyze trajectory, getting number of frames, file positions etc.
00358 static int analyze_traj(qmdata_t *data, orcadata *orca);
00359 
00360 // atom definitions & geometry
00361 static int get_input_structure(qmdata_t *data, orcadata *orca);
00362 
00363 // reading coord block
00364 static int get_coordinates(FILE *file, qm_atom_t **atoms, int unit, int *numatoms);
00365 
00366 // reading first trajectory frame
00367 static int read_first_frame(qmdata_t *data);
00368 
00369 static int get_basis(qmdata_t *data);
00370 
00371 static int get_wavefunction(qmdata_t *data, qm_timestep_t *ts, qm_wavefunction_t *wf);
00372 
00373 // main routine for extracting "step"-dependent info from the output file
00374 static int get_traj_frame(qmdata_t *data, qm_atom_t *atoms, int natoms);
00375 
00376 static int get_scfdata(qmdata_t *data, qm_timestep_t *ts);
00377 
00378 static int get_population(qmdata_t *data, qm_timestep_t *ts);
00379 
00380 static int check_add_wavefunctions(qmdata_t *data, qm_timestep_t *ts);
00381 
00382 static int fill_basis_arrays(qmdata_t *data);
00383 
00384 static int shelltype_int(char type);
00385 
00386 static CoeffRowBlock convertPure(CoeffRowBlock pureBlock);
00387 
00388 // for VMD
00389 static int read_orca_metadata(void *mydata, molfile_qm_metadata_t *metadata);
00390 static int read_orca_rundata(void *mydata, molfile_qm_t *qm_data);
00391 static int read_qm_timestep_metadata(void *mydata, molfile_qm_timestep_metadata_t *meta);
00392 
00393 static int read_timestep(void *mydata, int natoms,
00394        molfile_timestep_t *ts, molfile_qm_metadata_t *qm_metadata,
00395        molfile_qm_timestep_t *qm_ts);
00396 
00397 static int read_timestep_metadata(void *mydata, molfile_timestep_metadata_t *meta);
00398 
00399 static void print_input_data(qmdata_t *data);
00400 
00401 /*
00402 String manipulation
00403 */
00404 template<typename Out>
00405 void split(const std::string &s, char delim, Out result);
00406 
00407 /*************************************************************
00408  *
00409  * MAIN ORCA CODE PART
00410  *
00411  **************************************************************/
00412 
00413 static void* open_orca_read(const char* filename, const char* filetype, int *natoms) {
00414   FILE* fd;
00415   qmdata_t *data = NULL;
00416 
00417   #ifdef DEBUGGING_ORCA
00418     printf("DEBUG: Open Orca Read called: %s\n", filename);
00419   #endif
00420 
00421   orcadata* orca;
00422 
00423   // open the orca output files
00424   fd = fopen(filename, "rb");
00425   if (!fd) {
00426     PRINTERR;
00427     return NULL;
00428   }
00429 
00430   // initialize/allocate qm data (qmplugin.h)
00431   data = init_qmdata();
00432   if (data == NULL) {
00433     PRINTERR;
00434     return NULL;
00435   }
00436 
00437   // orca specific information alloc
00438   orca = (orcadata *) calloc(1, sizeof(orcadata));
00439   orca->version = 0;
00440   data->format_specific_data = orca;
00441 
00442   // filename in qm data struct
00443   data->file = fd;
00444 
00445   if(have_orca(data, orca)) {
00446     if (orca->version != 0) {
00447       printf("orcaplugin) Orca version: %d.%d.%d \n", orca->digits[0],
00448                                             orca->digits[1],orca->digits[2]);
00449     } else {
00450       printf("orcaplugin) Orca version not supported: %d.%d.%d \n", orca->digits[0],
00451                                             orca->digits[1],orca->digits[2]);
00452       return NULL;
00453     }
00454     if (parse_static_data(data, natoms) == FALSE) {
00455       return NULL;
00456     }
00457   } else {
00458     printf("orcaplugin) This is not an Orca output file!\n");
00459     return NULL;
00460   }
00461 
00462   return data;
00463 }
00464 
00465 
00466 
00467 static int read_orca_structure(void *mydata, int *optflags, molfile_atom_t *atoms) {
00468   qmdata_t *data = (qmdata_t *)mydata;
00469   qm_atom_t *cur_atom;
00470   molfile_atom_t *atom;
00471   int i = 0;
00472   *optflags = MOLFILE_ATOMICNUMBER;
00473 
00474   cur_atom = data->atoms;
00475 
00476   for (i=0; i<data->numatoms; i++) {
00477     atom = atoms+i;
00478     strncpy(atom->name, cur_atom->type, sizeof(atom->name));
00479     strncpy(atom->type, cur_atom->type, sizeof(atom->type));
00480     strncpy(atom->resname,"", sizeof(atom->resname));
00481     atom->resid = 1;
00482     atom->chain[0] = '\0';
00483     atom->segid[0] = '\0';
00484     atom->atomicnumber = cur_atom->atomicnum;
00485     #ifdef DEBUGGING_ORCA
00486       printf("orcaplugin) atomicnum[%d] = %d\n", i, atom->atomicnumber);
00487     #endif
00488 
00489     /* if (data->have_mulliken)
00490     atom->charge = data->qm_timestep->mulliken_charges[i];
00491     */
00492     cur_atom++;
00493   }
00494 
00495   return MOLFILE_SUCCESS;
00496 }
00497 
00498 
00499 
00500 template<typename Out>
00501 void split(const std::string &s, char delim, Out result) {
00502     std::stringstream ss;
00503     ss.str(s);
00504     std::string item;
00505     while (std::getline(ss, item, delim)) {
00506         *(result++) = item;
00507     }
00508 }
00509 
00510 
00511 
00512 std::vector<std::string> split(const std::string &s, char delim) {
00513     std::vector<std::string> elems;
00514     split(s, delim, std::back_inserter(elems));
00515     return elems;
00516 }
00517 
00518 
00519 
00520 std::string trim(const std::string& str, 
00521                  const std::string& whitespace = " \t") {
00522   const int strBegin = str.find_first_not_of(whitespace);
00523   if (strBegin == std::string::npos)
00524     return ""; // no content
00525 
00526   const int strEnd = str.find_last_not_of(whitespace);
00527   const int strRange = strEnd - strBegin + 1;
00528 
00529   return str.substr(strBegin, strRange);
00530 }
00531 
00532 
00533 
00534 std::string reduce(const std::string& str,
00535                    const std::string& fill = " ",
00536                    const std::string& whitespace = " \t\n") {
00537   // trim first
00538   std::string result = trim(str, whitespace);
00539 
00540   // replace sub ranges
00541   size_t beginSpace = result.find_first_of(whitespace);
00542   while (beginSpace != std::string::npos) {
00543     size_t endSpace = result.find_first_not_of(whitespace, beginSpace);
00544     const size_t range = endSpace - beginSpace;
00545 
00546     result.replace(beginSpace, range, fill);
00547 
00548     const size_t newStart = beginSpace + fill.length();
00549     beginSpace = result.find_first_of(whitespace, newStart);
00550   }
00551 
00552   return result;
00553 }
00554 
00555 
00556 
00557 static int get_job_info(qmdata_t *data) {
00558   long filepos;
00559   char buffer[BUFSIZ];
00560   filepos = ftell(data->file);
00561   if (goto_keyline(data->file, "INPUT FILE", NULL)) {
00562     eatline(data->file, 3);
00563   }
00564   std::vector<std::string> inputFile;
00565   int endOfInput = 0;
00566   while(!endOfInput) {
00567     GET_LINE(buffer, data->file);
00568     std::string test(buffer);
00569     size_t ln = test.find_first_of("123456789");
00570     size_t beforeLine = test.find_first_of(">");
00571     //stoi is not supported.
00572     //int lineNumber = stoi(test.substr(ln,beforeLine-ln));
00573     int lineNumber = atoi(test.substr(ln,beforeLine-ln).c_str());
00574     std::string lContent = trim(test.substr(beforeLine+1));
00575     inputFile.push_back(lContent);
00576     if (lineNumber == 1) {
00577       if (lContent[0] != '!') {
00578         std::cout << "orcaplugin) File corrupted in"
00579                      " input section. Exiting" << std::endl;
00580         return FALSE;
00581       } else {
00582         std::cout << "Found commands." << std::endl;
00583       }
00584 
00585       const char *semiempiricals[3] = {"MNDO","PM3","AM1"};
00586 
00587       //for (auto method : semiempiricals) {
00588       for (unsigned int i = 0; i < 3 /* semiempiricals.size()*/ ; i++){
00589         std::string method(semiempiricals[i]);
00590         if (lContent.find(method) != std::string::npos) {
00591           const char *m = method.c_str();
00592           strncpy(data->gbasis, m, sizeof(char)*strlen(m));
00593           strncpy(data->basis_string, "STO-3G", sizeof(char)*strlen("STO-3G"));
00594           std::cout << "orcaplugin) semiempirical Method used." << std::endl;
00595           break;
00596         }
00597       }
00598     }
00599 
00600     if (test.find("END OF INPUT") !=std::string::npos) {
00601       endOfInput = 1;
00602     }
00603   }
00604 
00605   data->runtype = MOLFILE_RUNTYPE_ENERGY;
00606   std::string lower = inputFile[0];
00607   std::transform(lower.begin(), lower.end(), lower.begin(), ::tolower);
00608   if (lower.find("opt") != std::string::npos) {
00609     data->runtype = MOLFILE_RUNTYPE_OPTIMIZE;
00610     std::cout << "orcaplugin) Optimization loaded." << std::endl;
00611   } else if (lower.find("engrad") != std::string::npos) {
00612     data->runtype = MOLFILE_RUNTYPE_GRADIENT;
00613   }
00614 
00615   int totalCharge;
00616   if (goto_keyline(data->file, "Total Charge", NULL)) {
00617     GET_LINE(buffer,data->file);
00618     std::string chargeLine(buffer);
00619     std::vector<std::string> chargeVec = split(reduce(chargeLine), ' ');
00620     //stoi is not supported
00621     //totalCharge = stoi(*(chargeVec.end()-1));
00622     totalCharge = (atoi((*(chargeVec.end()-1)).c_str()));
00623     std::cout << "orcaplugin) Found molecule charge: " << totalCharge << std::endl;
00624     data->totalcharge = totalCharge;
00625   } else {
00626     std::cout << "orcaplugin) No molecule charge found. Exiting" << std::endl;
00627     data->totalcharge = 0;
00628     // return FALSE;
00629   }
00630 
00631   int multiplicity;
00632   if (goto_keyline(data->file, "Multiplicity", NULL)) {
00633     GET_LINE(buffer,data->file);
00634     std::string multLine(buffer);
00635     std::vector<std::string> multVec = split(reduce(multLine), ' ');
00636     //stoi is not supported
00637     //multiplicity = stoi(*(multVec.end()-1));
00638     multiplicity = atoi((*(multVec.end()-1)).c_str());
00639     std::cout << "orcaplugin) Found molecule multiplicity: " << multiplicity << std::endl;
00640     data->multiplicity = multiplicity;
00641   } else {
00642     std::cout << "orcaplugin) No molecule multiplicity found. Exiting" << std::endl;
00643     data->multiplicity = -1;
00644     // return FALSE;
00645   }
00646 
00647   int nEl;
00648   if (goto_keyline(data->file, "Number of Electrons", NULL)) {
00649     GET_LINE(buffer,data->file);
00650     std::string nElLine(buffer);
00651     std::vector<std::string> nElVec = split(reduce(nElLine), ' ');
00652     //stoi is not supported
00653     //nEl = stoi(*(nElVec.end()-1));
00654     nEl = atoi((*(nElVec.end()-1)).c_str());
00655     std::cout << "orcaplugin) Found number of electrons: " << nEl << std::endl;
00656     data->num_electrons = nEl;
00657   } else {
00658     std::cout << "orcaplugin) Number of electrons not found. Exiting" << std::endl;
00659     data->num_electrons = -1;
00660     // return FALSE;
00661   }
00662 
00663   rewind(data->file);
00664   return TRUE;
00665 }
00666 
00667 
00668 
00669 static int have_orca(qmdata_t *data, orcadata* orca) {
00670   int programLine;
00671   int versionLine;
00672   char buffer[BUFSIZ];
00673   int mainVersion, secondDigit, thirdDigit;
00674   buffer[0] = '\0';
00675   programLine = goto_keyline(data->file, "O   R   C   A", NULL);
00676   if (programLine != 1) {
00677     return FALSE;
00678   }
00679 
00680   versionLine = goto_keyline(data->file, "Program Version", NULL);
00681   // thisline(data->file);
00682   GET_LINE(buffer, data->file);
00683   if (strstr(buffer,"Version") != NULL) {
00684     sscanf(buffer, "%*s %*s %d.%d.%d", &mainVersion, &secondDigit, &thirdDigit);
00685     #ifdef DEBUGGING_ORCA
00686       printf("DEBUG: build: %d.%d.%d\n", mainVersion, secondDigit, thirdDigit);
00687     #endif
00688     int build[3] = { mainVersion, secondDigit, thirdDigit };
00689     for (size_t i = 0; i < 3; i++) {
00690         orca->digits[i] = build[i];
00691     }
00692     switch (mainVersion) {
00693       case ORCA4:
00694         orca->version = 2;
00695         break;
00696       case ORCA3:
00697         orca->version = 1;
00698         break;
00699       default:
00700         orca->version = 0;
00701         break;
00702     }
00703   } else {
00704     PRINTERR;
00705     return FALSE;
00706   }
00707 
00708   return TRUE;
00709 }
00710 
00711 
00712 
00713 static int parse_static_data(qmdata_t *data, int* natoms) {
00714   orcadata *orca = (orcadata *)data->format_specific_data;
00715   if (!get_job_info(data)) return FALSE;
00716 
00717   if (!get_input_structure(data, orca)) return FALSE;
00718 
00719   if (!get_basis(data)) return FALSE;
00720 
00721   if (!analyze_traj(data, orca)) {
00722     printf("orcaplugin) WARNING: Truncated or abnormally terminated file!\n\n");
00723   }
00724 
00725   *natoms = data->numatoms;
00726 
00727   read_first_frame(data);
00728 
00729   // print_input_data(data);
00730 
00731   return TRUE;
00732 }
00733 
00734 
00735 /**********************************************************
00736  *
00737  * Read the input atom definitions and geometry
00738  *
00739  **********************************************************/
00740 static int get_input_structure(qmdata_t *data, orcadata *orca) {
00741   char buffer[BUFSIZ];
00742   int numatoms = -1;
00743   int bohr;
00744   long filepos;
00745   filepos = ftell(data->file);
00746 
00747   if (goto_keyline(data->file, "CARTESIAN COORDINATES (ANGSTROEM)", NULL)) {
00748     GET_LINE(buffer, data->file);
00749     // thisline(data->file);
00750     // UNITS ARE ANGSTROEM
00751     bohr = 0;
00752     // sscanf()
00753   } else {
00754     printf("orcaplugin) No cartesian coordinates in ANGSTROEM found.\n");
00755     return FALSE;
00756   }
00757 
00758   // skip the ---- line
00759   eatline(data->file, 1);
00760   /* Read the coordinate block */
00761   if (get_coordinates(data->file, &data->atoms, bohr, &numatoms))
00762     data->num_frames_read = 0;
00763   else {
00764     printf("orcaplugin) Bad atom coordinate block!\n");
00765     return FALSE;
00766   }
00767 
00768   data->numatoms = numatoms;
00769   return TRUE;
00770 }
00771 
00772 
00773 static int get_basis(qmdata_t *data) {
00774 //  orcadata *orca = (orcadata *)data->format_specific_data; // unused?!?
00775   char buffer[BUFSIZ];
00776   char word[4][BUFSIZ];
00777   int i = 0;
00778   int numread, numshells;
00779   shell_t *shell;
00780   long filepos;
00781   int semiempirical = 0;
00782 
00783   if (!strcmp(data->gbasis, "MNDO") ||
00784       !strcmp(data->gbasis, "AM1")  ||
00785       !strcmp(data->gbasis, "PM3")) {
00786     semiempirical = 1;
00787   }
00788 
00789   /* Search for "ATOMIC BASIS SET" line */
00790   if ((pass_keyline(data->file, "BASIS SET IN INPUT FORMAT", NULL) != FOUND) &&
00791       (pass_keyline(data->file, "GAUSSIAN BASIS SET", NULL) != FOUND)) {
00792     printf("orcaplugin) No basis set found!\n");
00793     return FALSE;
00794   }
00795 
00796   /* initialize buffers */
00797   buffer[0] = '\0';
00798   for (i=0; i<3; i++)
00799     word[i][0] = '\0';
00800 
00801 
00802   /* skip the next 3 lines */
00803   // eatline(data->file, 2);
00804 
00805   /* Allocate space for the basis for all atoms */
00806   data->basis_set = (basis_atom_t*)calloc(data->numatoms, sizeof(basis_atom_t));
00807 
00808   filepos = ftell(data->file);
00809   i = 0; /* basis atom counter */
00810   char elementName[11];
00811   int finished = FALSE;
00812 
00813   basis_atom_t* tempBasis = (basis_atom_t*)calloc(1, sizeof(basis_atom_t));
00814 
00815   prim_t *prim;
00816 
00817   std::vector<std::string> readElements;
00818 #if 0
00819   for (size_t atom = 0; atom < data->numatoms; atom++) {
00820 // WTF?!?!?!
00821   }
00822 #endif
00823 
00824   while (!finished && !semiempirical) {
00825     // printf("Trying to read bf. \n");
00826     if (pass_keyline(data->file, "Basis set for element", NULL) == FOUND ) {
00827       GET_LINE(buffer, data->file);
00828       numread = sscanf(buffer,"%s %s",&word[0][0], &word[1][0]);
00829       strcpy(elementName, &word[1][0]);
00830       printf("New element found: %s\n", &word[1][0]);
00831       std::string atype(&word[1][0]);
00832       if (std::find(readElements.begin(), readElements.end(), atype) != readElements.end()) {
00833         break;
00834       } else {
00835         readElements.push_back(atype);
00836       }
00837       int elementCompleted = 0;
00838 
00839 
00840       shell = (shell_t*)calloc(1, sizeof(shell_t));
00841       numshells = 0;
00842 //      int readingShell = 0;
00843       int primcounter = 0;
00844 
00845       // this is very sloppy at the moment...
00846       // for PM3 etc. Orca prints the bf per atom...
00847 
00848       while(!elementCompleted) {
00849         GET_LINE(buffer, data->file);
00850         numread = sscanf(buffer,"%s %s %s",&word[0][0], &word[1][0],&word[2][0]);
00851         // printf("numread: %d -- %s %s %s \n",numread, &word[0][0], &word[1][0],&word[2][0]);
00852         switch (numread) {
00853           case 1:
00854             if (strcmp(trimleft(trimright(&word[0][0])), "end")) {
00855               // printf("Section ended. \n");
00856               elementCompleted = 1;
00857               break;
00858             }
00859           case 2:
00860             shell[numshells].numprims = atoi(trimleft(trimright(&word[1][0])));
00861             shell[numshells].type = shelltype_int(word[0][0]);
00862             // printf("orcaplugin) Type: %d NPrims: %d\n", shell[numshells].type, shell[numshells].numprims);
00863             primcounter = 0;
00864             prim = (prim_t*)calloc(shell[numshells].numprims, sizeof(prim_t));
00865             // printf("!! Address of prim is %p\n", (void *)prim);
00866             shell[numshells].prim = prim;
00867             numshells++;
00868             if (numshells) {
00869               shell = (shell_t*)realloc(shell, (numshells+1)*sizeof(shell_t));
00870             }
00871             break;
00872           case 3:
00873             prim[primcounter].exponent = atof(&word[1][0]);
00874             prim[primcounter].contraction_coeff = atof(&word[2][0]);
00875             // printf("%f - %f\n", prim[primcounter].exponent, prim[primcounter].contraction_coeff);
00876             primcounter++;
00877             break;
00878           default:
00879             printf("orcaplugin) Unknown line in basis functions.\n");
00880             elementCompleted = 1;
00881             break;
00882         }
00883       }
00884       // printf("Number of shells: %d \n", numshells);
00885       strcpy(tempBasis[i].name, elementName);
00886       tempBasis[i].numshells = numshells;
00887       tempBasis[i].shell = shell;
00888       i++;
00889       if (i) {
00890         tempBasis = (basis_atom_t*)realloc(tempBasis, (i+1)*sizeof(basis_atom_t));
00891       }
00892       // set prim to nullpointer!
00893       prim = NULL;
00894     } else {
00895       finished = TRUE;
00896       printf("orcaplugin) Reading basis set finished! \n");
00897     }
00898   }
00899 
00900   finished = 0;
00901   // semiempirical sto-3g
00902   int ngauss = 3;
00903   int atomCounter = 0;
00904   int shellNumber;
00905   while(!finished && semiempirical) {
00906     if(goto_keyline(data->file,"shells",NULL) == FOUND && atomCounter <(data->numatoms-1)) {
00907       // thisline(data->file);
00908       GET_LINE(buffer, data->file);
00909       std::string lineString(buffer);
00910       std::vector<std::string> elements = split(reduce(lineString), ' ');
00911 
00912       //stoi is not supported
00913       //shellNumber = stoi(elements[2]);
00914       shellNumber = atoi(elements[2].c_str());
00915 
00916       // std::cout << "shell number: " << shellNumber << std::endl;
00917       shell = (shell_t*) calloc(shellNumber, sizeof(shell_t));
00918       data->basis_set[atomCounter].shell = shell;
00919       data->basis_set[atomCounter].numshells = shellNumber;
00920       for (size_t shellIdx = 0; shellIdx < shellNumber; ++shellIdx) {
00921         GET_LINE(buffer, data->file);
00922         prim = (prim_t*) calloc(3, sizeof(prim_t));
00923         shell[shellIdx].prim = prim;
00924         shell[shellIdx].numprims = 3;
00925         shell[shellIdx].type = shellIdx;
00926         for (size_t nbas = 0; nbas < ngauss; ++nbas) {
00927           GET_LINE(buffer, data->file);
00928           std::string l(buffer);
00929           std::vector<std::string> coeff = split(reduce(l), ' ');
00930 #if 1
00931           prim[nbas].exponent = atof(coeff[0].c_str());
00932           prim[nbas].contraction_coeff = atof(coeff[1].c_str());
00933 #else
00934           prim[nbas].exponent = stod(coeff[0]);
00935           prim[nbas].contraction_coeff = stod(coeff[1]);
00936 #endif
00937         }
00938       }
00939 
00940       data->num_shells += shellNumber;
00941       data->num_basis_atoms++;
00942       data->num_basis_funcs += 3;
00943       strncpy(data->basis_set[atomCounter].name, data->atoms[atomCounter].type, 11);
00944       atomCounter++;
00945     } else {
00946       finished = TRUE;
00947       prim = NULL;
00948       std::cout << "orcaplugin) Reading STO-3G basis for semiempirical method finished." << std::endl;
00949       // free unused stuff
00950       free(tempBasis);
00951       // We return here without further ado.
00952       return fill_basis_arrays(data);
00953     }
00954   }
00955 
00956   // As we read GTOs from the Orca output file, we need to
00957   // loop over all atoms and assign the basis functions
00958 
00959   // printf("orcaplugin) Parsed basis set of %d elements. \n", i);
00960   // for (size_t j = 0; j < i; j++) {
00961   //   printf("Element: %s\n", tempBasis[j].name);
00962   //   printf("- NShells: %d\n", tempBasis[j].numshells);
00963   //   for (size_t k = 0; k < tempBasis[j].numshells; k++) {
00964   //     printf("--- NPrims: %d \n", tempBasis[j].shell[k].numprims);
00965   //     for (size_t o = 0; o < tempBasis[j].shell[k].numprims; o++) {
00966   //       float expo = tempBasis[j].shell[k].prim[o].exponent;
00967   //       float cont = tempBasis[j].shell[k].prim[o].contraction_coeff;
00968   //       printf("----- E= %f , C= %f \n", expo, cont);
00969   //     }
00970   //   }
00971   // }
00972 
00973   // Allocate an array of zeros to store whether the tempBasis
00974   // of the same index was actually used.
00975   int *tempBasisUsed = (int*) calloc(i, sizeof(int));
00976 
00977   const char* currentElement;
00978   basis_atom_t* currentBasis;
00979   for (size_t n = 0; n < data->numatoms; n++) {
00980     currentElement = data->atoms[n].type;
00981     for (size_t j = 0; j < i; j++) {
00982       if (!strcmp(currentElement, tempBasis[j].name)) {
00983         // printf("orcaplugin) found basis for element %s\n", currentElement);
00984         currentBasis = &tempBasis[j];
00985         tempBasisUsed[j] = 1;
00986       }
00987     }
00988 
00989     // printf("orcaplugin) Basis for element %s has %d shells.\n", currentElement, currentBasis->numshells);
00990     data->basis_set[n].shell = (shell_t *) calloc(currentBasis->numshells, sizeof(shell_t));
00991     memcpy(data->basis_set[n].shell, currentBasis->shell, currentBasis->numshells * sizeof(shell_t));
00992     data->basis_set[n].numshells = currentBasis->numshells;
00993     data->num_shells += currentBasis->numshells;
00994     for (size_t p = 0; p < currentBasis->numshells; ++p) {
00995       data->basis_set[n].shell[p].prim = (prim_t *) calloc(currentBasis->shell[p].numprims, sizeof(prim_t));
00996       memcpy(data->basis_set[n].shell[p].prim, currentBasis->shell[p].prim, currentBasis->shell[p].numprims * sizeof(prim_t));
00997       data->num_basis_funcs += currentBasis->shell[p].numprims;
00998     }
00999     data->num_basis_atoms++;
01000     strncpy(data->basis_set[n].name, currentBasis->name, 11);
01001     currentBasis = NULL;
01002     currentElement  = NULL;
01003   }
01004 
01005   // cleaning up memory for tempBasis and its shells
01006   for (size_t idx = 0; idx < i; ++idx) {
01007     for (size_t shellIdx = 0; shellIdx < tempBasis[idx].numshells; shellIdx++) {
01008       shell_t* cshell = &tempBasis[idx].shell[shellIdx];
01009       free(cshell->prim);
01010       cshell->prim = NULL;
01011     }
01012     free(tempBasis[idx].shell);
01013     tempBasis[idx].shell = NULL;
01014   }
01015   free(tempBasis);
01016   tempBasis = NULL;
01017   shell = NULL;
01018   printf("orcaplugin) Parsed %d uncontracted basis functions.\n", data->num_basis_funcs);
01019   free(tempBasisUsed);
01020 
01021   // /* allocate and populate flat arrays needed for molfileplugin */
01022   return fill_basis_arrays(data);
01023 }
01024 
01025 
01026 static int get_coordinates(FILE *file, qm_atom_t **atoms, int unit,
01027                            int *numatoms) {
01028   int i = 0;
01029   int growarray = 0;
01030 
01031   if (*numatoms<0) {
01032     *atoms = (qm_atom_t*)calloc(1, sizeof(qm_atom_t));
01033     growarray = 1;
01034   }
01035 
01036   /* Read in the coordinates until an empty line is reached.
01037    * We expect 5 entries per line */
01038   while (1) {
01039     char buffer[BUFSIZ];
01040     char atname[BUFSIZ];
01041     float atomicnum;
01042     float x,y,z;
01043     int n;
01044     qm_atom_t *atm;
01045 
01046     GET_LINE(buffer, file);
01047     // thisline(file);
01048 
01049     /* For FMO there is an additional atom index in the
01050      * second column. Try both variants: */
01051     n = sscanf(buffer,"%s %f %f %f",atname,&x,&y,&z);
01052     // printf("%s\n", atname);
01053     if (n!=4) {
01054       // n = sscanf(buffer,"%s %f %f %f %f",atname,&atomicnum,&x,&y,&z);
01055       break;
01056     }
01057     // if (n!=5 && n!=6) break;
01058 
01059     if (growarray && i>0) {
01060       *atoms = (qm_atom_t*)realloc(*atoms, (i+1)*sizeof(qm_atom_t));
01061     }
01062     atm = (*atoms)+i;
01063 
01064     // just get the atomic number from periodic_table.h
01065     atomicnum = get_pte_idx(atname);
01066 
01067     strncpy(atm->type, atname, sizeof(atm->type));
01068     atm->atomicnum = floor(atomicnum+0.5); /* nuclear charge */
01069     // printf("coor: %s %d %f %f %f\n", atm->type, atm->atomicnum, x, y, z);
01070 
01071     /* if coordinates are in Bohr convert them to Angstrom */
01072     if (unit==BOHR) {
01073       x *= BOHR_TO_ANGS;
01074       y *= BOHR_TO_ANGS;
01075       z *= BOHR_TO_ANGS;
01076     }
01077 
01078     atm->x = x;
01079     atm->y = y;
01080     atm->z = z;
01081     i++;
01082   }
01083 
01084   /* If file is broken off in the middle of the coordinate block
01085    * we cannot use this frame. */
01086   if (*numatoms>=0 && *numatoms!=i) {
01087     (*numatoms) = i;
01088     return FALSE;
01089   }
01090 
01091   (*numatoms) = i;
01092   return TRUE;
01093 }
01094 
01095 
01096 
01097 /* Read the first trajectory frame. */
01098 static int read_first_frame(qmdata_t *data) {
01099   /* Try reading the first frame.
01100    * If there is only one frame then also read the
01101    * final wavefunction. */
01102   if (!get_traj_frame(data, data->atoms, data->numatoms)) {
01103     return FALSE;
01104   }
01105 
01106   return TRUE;
01107 }
01108 
01109 
01110 
01111 /******************************************************
01112  *
01113  * this function extracts the trajectory information
01114  * from the output file
01115  *
01116  * *****************************************************/
01117 static int get_traj_frame(qmdata_t *data, qm_atom_t *atoms,
01118                           int natoms) {
01119 //  orcadata *orca = (orcadata *)data->format_specific_data;
01120   qm_timestep_t *cur_ts;
01121   char buffer[BUFSIZ];
01122   char word[BUFSIZ];
01123   int units=ANGSTROM;
01124   buffer[0] = '\0';
01125   word[0]   = '\0';
01126 
01127   printf("orcaplugin) Timestep %d:\n", data->num_frames_read);
01128   printf("orcaplugin) ============\n");
01129 
01130   cur_ts = data->qm_timestep + data->num_frames_read;
01131 
01132   // debugging the trajectory reading file positions
01133   // printf("nfread: %d \n", data->num_frames_read);
01134   if (!data->filepos_array) {
01135     printf("filepos array empty!!!\n");
01136     return FALSE;
01137   }
01138 
01139   fseek(data->file, data->filepos_array[data->num_frames_read], SEEK_SET);
01140 
01141   /*
01142   * distinguish between job types
01143   * at the moment, only Single Points will work
01144   * lines 2840 - 3122 in gamessplugin.c
01145    */
01146 
01147   // reading geometries...
01148 
01149   if ((data->runtype == MOLFILE_RUNTYPE_OPTIMIZE || data->runtype == MOLFILE_RUNTYPE_GRADIENT) && data->num_frames > 1) {
01150     if (goto_keyline(data->file, "CARTESIAN COORDINATES (ANGSTROEM)", NULL)) {
01151       GET_LINE(buffer, data->file);
01152       // thisline(data->file);
01153       // UNITS ARE ANGSTROEM
01154       // bohr = 0;
01155       // sscanf()
01156     } else {
01157       printf("orcaplugin) No cartesian coordinates in ANGSTROEM found.\n");
01158     }
01159     // skip the ---- line
01160     eatline(data->file, 1);
01161     if (!get_coordinates(data->file, &data->atoms, units, &natoms)) {
01162       printf("orcaplugin) Couldn't find coordinates for timestep %d\n", data->num_frames_read);
01163     }
01164   }
01165 
01166   if (get_scfdata(data, cur_ts) == FALSE) {
01167     printf("orcaplugin) Couldn't find SCF iterations for timestep %d\n",
01168            data->num_frames_read);
01169   }
01170 
01171   /* Try reading canonical alpha/beta wavefunction */
01172   check_add_wavefunctions(data, cur_ts);
01173 
01174   /* Read point charged */
01175   if (!cur_ts->have_mulliken && get_population(data, cur_ts)) {
01176     printf("orcaplugin) Mulliken charges found\n");
01177   }
01178 
01179   /* Read the energy gradients (= -forces on atoms) */
01180   // if (get_gradient(data, cur_ts)) {
01181   //   printf("orcaplugin) Energy gradient found.\n");
01182   // }
01183 
01184   /* If this is the last frame of the trajectory and the file
01185    * wasn't truncated and the program didn't terminate
01186    * abnormally then read the final wavefunction. */
01187   if ((data->runtype == MOLFILE_RUNTYPE_OPTIMIZE ||
01188        data->runtype == MOLFILE_RUNTYPE_SADPOINT) &&
01189       (data->num_frames_read+1 == data->num_frames &&
01190        (data->status == MOLFILE_QMSTATUS_UNKNOWN ||
01191         data->status == MOLFILE_QMSTATUS_OPT_CONV ||
01192         data->status == MOLFILE_QMSTATUS_OPT_NOT_CONV))) {
01193 
01194     /* We need to jump over the end of the trajectory because
01195      * this is also the keystring for get_wavefunction() to
01196      * bail out. */
01197     if (data->status == MOLFILE_QMSTATUS_OPT_CONV ||
01198         data->status == MOLFILE_QMSTATUS_OPT_NOT_CONV) {
01199       fseek(data->file, data->end_of_traj, SEEK_SET);
01200       std::cout << "orcaplugin) Finished trajectory." << std::endl;
01201     }
01202 
01203     /* Try to read final wavefunction and orbital energies
01204      * A preexisting canonical wavefunction for this timestep
01205      * with the same characteristics (spin, exci, info) will
01206      * be overwritten by the final wavefuntion if it has more
01207      * orbitals. */
01208     check_add_wavefunctions(data, cur_ts);
01209   }
01210 
01211   data->num_frames_read++;
01212   #ifdef DEBUGGING_ORCA
01213     std::cout << "orcaplugin) Frames read: " << data->num_frames_read << std::endl;
01214   #endif
01215 
01216   return TRUE;
01217 }
01218 
01219 
01220 static bool only_numbers(const std::vector<std::string>& l) {
01221   //for (auto s : l) {
01222   for (int i = 0; i < l.size(); i++){
01223     if (l[i].find_first_not_of("0123456789-.") != std::string::npos) {
01224       return false;
01225     }
01226   }
01227   return true;
01228 }
01229 
01230 
01231 static int get_scfdata(qmdata_t *data, qm_timestep_t *ts) {
01232   char buffer[BUFSIZ];
01233   long filepos;
01234   filepos = ftell(data->file);
01235   if (!goto_keyline(data->file, "SCF ITERATIONS", NULL)) {
01236     fseek(data->file, filepos, SEEK_SET);
01237     ts->num_scfiter = 0;
01238     return FALSE;
01239   }
01240   eatline(data->file, 2);
01241   GET_LINE(buffer, data->file);
01242   std::string scfIterLine(buffer);
01243   int ncols = (split(reduce(scfIterLine), ' ')).size();
01244   int currentCols, numiter = 0;
01245   while (scfIterLine.find("SUCCESS") == std::string::npos && scfIterLine.find("ERROR") == std::string::npos) {
01246     GET_LINE(buffer, data->file);
01247     scfIterLine = buffer;
01248     std::vector<std::string> currentCol = split(reduce(scfIterLine), ' ');
01249     currentCols = currentCol.size();
01250     if (currentCols == ncols && only_numbers(currentCol)) {
01251       numiter++;
01252     }
01253   }
01254   ts->num_scfiter = numiter;
01255   std::cout << "orcaplugin) Number of SCF iterations: " << numiter << std::endl;
01256   return TRUE;
01257 }
01258 
01259 static int get_population(qmdata_t *data, qm_timestep_t *ts) {
01260   int i;
01261   char buffer[BUFSIZ];
01262   long filepos;
01263   ts->have_mulliken = FALSE;
01264   ts->have_lowdin   = FALSE;
01265   ts->have_esp      = FALSE;
01266   filepos = ftell(data->file);
01267 
01268   if (pass_keyline(data->file, "MULLIKEN ATOMIC CHARGES", NULL) != FOUND) {
01269     fseek(data->file, filepos, SEEK_SET);
01270     return FALSE;
01271   }
01272 
01273   /* Read Mulliken charges if present */
01274   ts->mulliken_charges = (double *)calloc(data->numatoms, sizeof(double));
01275 
01276   if (!ts->mulliken_charges) {
01277     PRINTERR;
01278     return FALSE;
01279   }
01280 
01281   // ts->lowdin_charges = (double *)calloc(data->numatoms, sizeof(double));
01282   //
01283   // if (!ts->lowdin_charges) {
01284   //   free(ts->mulliken_charges);
01285   //   ts->mulliken_charges = NULL;
01286   //   PRINTERR;
01287   //   return FALSE;
01288   // }
01289   eatline(data->file, 1);
01290 
01291   for (i=0; i<data->numatoms; i++) {
01292     int n;
01293     float mullcharge;
01294     GET_LINE(buffer, data->file);
01295     std::string currentLine(buffer);
01296     std::vector<std::string> elements = split(reduce(currentLine), ' ');
01297     n = elements.size();
01298     if (n!=4) {
01299       free(ts->mulliken_charges);
01300       // free(ts->lowdin_charges);
01301       ts->mulliken_charges = NULL;
01302       // ts->lowdin_charges   = NULL;
01303       return FALSE;
01304     }
01305     mullcharge = (float)atof((*(elements.end()-1)).c_str());
01306     ts->mulliken_charges[i] = mullcharge;
01307     // ts->lowdin_charges[i]   = lowcharge;
01308   }
01309 
01310   if (i!=data->numatoms) {
01311     free(ts->mulliken_charges);
01312     free(ts->lowdin_charges);
01313     ts->mulliken_charges = NULL;
01314     ts->lowdin_charges   = NULL;
01315     return FALSE;
01316   }
01317 
01318   ts->have_mulliken = TRUE;
01319   // ts->have_lowdin   = TRUE;
01320   return TRUE;
01321 }
01322 
01323 /*********************************************************
01324  *
01325  * Reads a set of wavefunctions for the current timestep.
01326  * These are typically the alpha and beta spin wavefunctions
01327  * or the MCSCF natural and optimized orbitals or the GVB
01328  * canonical orbitals and geminal pairs.
01329  *
01330  **********************************************************/
01331  // 3461
01332 static int check_add_wavefunctions(qmdata_t *data, qm_timestep_t *ts) {
01333   qm_wavefunction_t *wavef;
01334   int i, n=1;
01335 
01336   if (data->scftype==MOLFILE_SCFTYPE_UHF ||
01337       data->scftype==MOLFILE_SCFTYPE_GVB ||
01338       data->scftype==MOLFILE_SCFTYPE_MCSCF) {
01339     /* Try to read second wavefunction, e.g. spin beta */
01340     n = 2;
01341   }
01342 
01343   for (i=0; i<n; i++) {
01344     /* Allocate memory for new wavefunction */
01345     wavef = add_wavefunction(ts);
01346 
01347     /* Try to read wavefunction and orbital energies */
01348     if (get_wavefunction(data, ts, wavef) == FALSE) {
01349       /* Free the last wavefunction again. */
01350       del_wavefunction(ts);
01351 #ifdef DEBUGGING_ORCA
01352       printf("orcaplugin) No canonical wavefunction present for timestep %d\n", data->num_frames_read);
01353 #endif
01354       break;
01355 
01356     } else {
01357       char action[32];
01358       char spinstr[32];
01359       strcpy(spinstr, "");
01360       if (data->scftype==MOLFILE_SCFTYPE_UHF) {
01361         if (wavef->spin==SPIN_BETA) {
01362           strcat(spinstr, "spin  beta, ");
01363         } else {
01364           strcat(spinstr, "spin alpha, ");
01365         }
01366       }
01367 
01368       /* The last SCF energy is the energy of this electronic state */
01369       if (ts->scfenergies) {
01370         wavef->energy = ts->scfenergies[ts->num_scfiter-1];
01371       } else {
01372         wavef->energy = 0.f;
01373       }
01374 
01375       /* Multiplicity */
01376       wavef->mult = data->multiplicity;
01377 
01378 
01379       /* String telling wether wavefunction was added, updated
01380        * or ignored. */
01381       strcpy(action, "added");
01382 
01383       /* If there exists a canonical wavefunction of the same spin
01384        * we'll replace it */
01385       if (ts->numwave>1 && wavef->type==MOLFILE_WAVE_CANON) {
01386         int i, found =-1;
01387         for (i=0; i<ts->numwave-1; i++) {
01388           if (ts->wave[i].type==wavef->type &&
01389               ts->wave[i].spin==wavef->spin &&
01390               ts->wave[i].exci==wavef->exci &&
01391               !strncmp(ts->wave[i].info, wavef->info, MOLFILE_BUFSIZ)) {
01392             found = i;
01393             break;
01394           }
01395         }
01396         if (found>=0) {
01397           /* If the new wavefunction has more orbitals we
01398            * replace the old one for this step. */
01399           if (wavef->num_orbitals >
01400               ts->wave[found].num_orbitals) {
01401             /* Replace existing wavefunction for this step */
01402             replace_wavefunction(ts, found);
01403             sprintf(action, "%d updated", found);
01404           } else {
01405             /* Delete last wavefunction again */
01406             del_wavefunction(ts);
01407             sprintf(action, "matching %d ignored", found);
01408           }
01409           wavef = &ts->wave[ts->numwave-1];
01410         }
01411       }
01412 
01413       printf("orcaplugin) Wavefunction %s (%s):\n", action, wavef->info);
01414       printf("orcaplugin)   %d orbitals, %sexcitation %d, multiplicity %d\n",
01415              wavef->num_orbitals, spinstr, wavef->exci, wavef->mult);
01416     }
01417   }
01418 
01419   return i;
01420 }
01421 
01422 //std::vector<std::vector<int> > dAngMom{{1,0,1},{0,1,1},{1,1,0},{2,0,0},{0,2,0},{0,0,2}};
01423 //std::vector<std::vector<int> > fAngMom{{1,1,1},{2,1,0},{2,0,1},{1,2,0},{0,2,1},{1,0,2},
01424 //  {0,1,2},{3,0,0},{0,3,0},{0,0,3}};
01425 int dAngMom[18] = {1,0,1, 0,1,1, 1,1,0, 2,0,0, 0,2,0, 0,0,2};
01426 int fAngMom[30] = {1,1,1, 2,1,0, 2,0,1, 1,2,0, 0,2,1, 1,0,2, 0,1,2, 3,0,0, 0,3,0, 0,0,3};
01427 
01428 typedef enum {
01429   dshell, fshell
01430 } Shelltype;
01431 
01432 static CoeffRowBlock convertPure(CoeffRowBlock pureBlock) {
01433   Shelltype stype;
01434   CoeffRowBlock resultBlock;
01435   switch (pureBlock.size()) {
01436     case 5:
01437       stype = dshell;
01438       break;
01439     case 7:
01440       stype = fshell;
01441       break;
01442     default:
01443       return resultBlock;
01444   }
01445 
01446   Matrix *m = new Matrix(pureBlock);
01447   Matrix *multiplication;
01448   if (stype == dshell) {
01449     multiplication = Matrix::multiply(convD, m);
01450   } else if (stype == fshell) {
01451     multiplication = Matrix::multiply(convF, m);
01452   } else {
01453     return resultBlock;
01454   }
01455   // m->printMatrix();
01456   // multiplication->printMatrix();
01457   resultBlock = multiplication->toVector();
01458   delete multiplication;
01459   delete m;
01460 
01461   // OLD CODE!
01462   // first get the unchanged d-orbitals xz, yz, xy
01463   // std::vector<std::string> unchangedOrbitals{"xz", "yz", "xy"};
01464   // std::vector<int> unchangedIndices{1, 2, 4};
01465   // for (auto idx : unchangedIndices) {
01466   //   resultBlock.push_back(pureBlock[idx]);
01467   // }
01468   //
01469   // int alpha = 3;
01470   // int beta = 0;
01471   // // now convert the others
01472   //
01473   // int rows = pureBlock[0].size(); // should be 5 in this case
01474   //
01475   // std::vector<float> x2;
01476   // for (size_t i = 0; i < pureBlock[0].size(); i++) {
01477   //   x2.push_back(pureBlock[alpha][i]-pureBlock[beta][i]);
01478   // }
01479   //
01480   // std::vector<float> y2;
01481   // for (size_t i = 0; i < pureBlock[0].size(); i++) {
01482   //   y2.push_back(-pureBlock[beta][i]-pureBlock[alpha][i]);
01483   // }
01484   //
01485   // std::vector<float> z2;
01486   // for (size_t i = 0; i < pureBlock[0].size(); i++) {
01487   //   z2.push_back(2*pureBlock[beta][i]);
01488   // }
01489   //
01490   // resultBlock.push_back(x2);
01491   // resultBlock.push_back(y2);
01492   // resultBlock.push_back(z2);
01493 
01494   return resultBlock;
01495 }
01496 
01497 static int get_wavefunction(qmdata_t *data, qm_timestep_t *ts, qm_wavefunction_t *wf) {
01498   std::vector<float> orbitalEnergies;
01499   std::vector<int> orbitalOccupancies;
01500   std::vector<float> wavefunctionCoeffs;
01501   int num_orbitals = 0;
01502 
01503   char buffer[BUFSIZ];
01504   char word[6][BUFSIZ];
01505   long filepos;
01506   char *line;
01507 
01508   buffer[0] = '\0';
01509   int i = 0;
01510   for (i=0; i<6; i++) word[i][0] = '\0';
01511 
01512   if (wf == NULL) {
01513     PRINTERR;
01514     return FALSE;
01515   }
01516 
01517   wf->has_occup = FALSE;
01518   wf->has_orben = FALSE;
01519   wf->type = MOLFILE_WAVE_UNKNOWN;
01520   wf->spin = SPIN_ALPHA;
01521   wf->exci = 0;
01522   strncpy(wf->info, "unknown", MOLFILE_BUFSIZ);
01523 
01524   filepos = ftell(data->file);
01525 
01526   do {
01527     GET_LINE(buffer, data->file);
01528     line = trimleft(trimright(buffer));
01529     if(!strcmp(line, "MOLECULAR ORBITALS")) {
01530       wf->type = MOLFILE_WAVE_CANON;
01531       strncpy(wf->info, "canonical", MOLFILE_BUFSIZ);
01532     }
01533   } while(wf->type == MOLFILE_WAVE_UNKNOWN && strcmp(line, "FINAL SINGLE POINT ENERGY"));
01534 
01535   if(wf->type == MOLFILE_WAVE_UNKNOWN) {
01536     #ifdef DEBUGGING_ORCA
01537         printf("orcaplugin) get_wavefunction(): No wavefunction found!\n");
01538     #endif
01539         fseek(data->file, filepos, SEEK_SET);
01540         return FALSE;
01541   } else {
01542     #ifdef DEBUGGING_ORCA
01543       printf("orcaplugin) Found wavefunction of type %d.\n", wf->type);
01544     #endif
01545   }
01546 
01547   eatline(data->file, 1);
01548 
01549   // number of read values from line;
01550   int numReadOrbitalIndices = 0;
01551   int numReadEnergies = 0;
01552   int numReadOccupancies = 0;
01553   int numReadCoefficients = 0;
01554   float numberOfElectrons = 0;
01555   float occupiedOrbitals = 0;
01556   std::vector<int> numberContractedBf;
01557   std::vector<int> wfAngMoment;
01558   std::vector<std::string> orbitalNames;
01559   MoCoeff allCoefficients;
01560 
01561   // orbital indices
01562   int n[6];
01563   int wavefunctionRead = 0;
01564   int firstRead = 1;
01565   int haveAngMom = 0;
01566   while(!wavefunctionRead) {
01567     float coeff[6], energies[6];
01568     float occ[6];
01569     char dumpName[BUFSIZ];
01570     char dumpBasisFunc[BUFSIZ];
01571     filepos = ftell(data->file);
01572 
01573     // reads the orbital indices
01574     if (firstRead == 1) {
01575       GET_LINE(buffer, data->file);
01576       firstRead++;
01577     }
01578 
01579     numReadOrbitalIndices = sscanf(buffer, "%d %d %d %d %d %d", &n[0], &n[1], &n[2], &n[3], &n[4], &n[5]);
01580     if (!numReadOrbitalIndices || numReadOrbitalIndices == -1) {
01581       /* If there are no orbital indexes then this must be the
01582        * end of the wavefunction coefficient table. */
01583       fseek(data->file, filepos, SEEK_SET);
01584       wavefunctionRead = 1;
01585       break;
01586     }
01587 
01588     // reads the orbital energies
01589     GET_LINE(buffer, data->file);
01590     numReadEnergies = sscanf(buffer, "%f %f %f %f %f %f", &energies[0], &energies[1], &energies[2],
01591      &energies[3], &energies[4], &energies[5]);
01592     if (numReadEnergies != numReadOrbitalIndices) {
01593       printf("orcaplugin) Molecular Orbital section corrupted! energies.\n");
01594       break;
01595     }
01596 
01597     // store the energies in vector
01598     if (numReadEnergies != -1) {
01599       for (size_t c = 0; c < numReadEnergies; c++) {
01600         orbitalEnergies.push_back(energies[c]);
01601         // std::cout << "Energy: " <<energies[c]<< std::endl;
01602       }
01603     }
01604     // reads the orbital occupancies
01605     GET_LINE(buffer, data->file);
01606     numReadOccupancies = sscanf(buffer, "%f %f %f %f %f %f", &occ[0], &occ[1], &occ[2],
01607       &occ[3], &occ[4], &occ[5]);
01608     if (numReadOccupancies != numReadOrbitalIndices) {
01609       printf("orcaplugin) Molecular Orbital section corrupted! %d\n", numReadOccupancies);
01610       break;
01611     }
01612 
01613     // stores the occupancies in vector
01614     if(numReadOccupancies) {
01615       for (size_t c = 0; c < numReadOccupancies; c++) {
01616         orbitalOccupancies.push_back((int)occ[c]);
01617         // std::cout << "Occupancy: " << occ[c] << std::endl;
01618         numberOfElectrons += occ[c];
01619         if (occ[c]) {
01620           occupiedOrbitals++;
01621         }
01622       }
01623       num_orbitals += numReadOccupancies;
01624     }
01625 
01626     // skip --- line
01627     eatline(data->file, 1);
01628 
01629     std::vector<std::vector<float> > moCoefficients;
01630     // we expect as many coefficients as numReadOccupancies, numReadEnergies, numReadOrbitalIndices!
01631     // read them as long as we find new coefficients
01632     int readingBlock = 1;
01633     int coefficientNumber = 0;
01634     while(readingBlock) {
01635       GET_LINE(buffer, data->file);
01636       numReadCoefficients = sscanf(buffer, "%s %s %f %f %f %f %f %f", &dumpName, &dumpBasisFunc,
01637       &coeff[0], &coeff[1], &coeff[2],&coeff[3], &coeff[4], &coeff[5]);
01638       // the coefficient number is the number of read elements minus 2 bc. of the atom and bf name
01639       coefficientNumber = (numReadCoefficients - 2);
01640       if (coefficientNumber == numReadOrbitalIndices) {
01641         // std::cout << "found coeffs: " << dumpName << "," << dumpBasisFunc << std::endl;
01642         if (firstRead == 2 && !haveAngMom) {
01643           std::string bfn = dumpBasisFunc;
01644           size_t found = bfn.find_first_not_of("-0123456789 ");
01645           if (found!=std::string::npos) {
01646             std::string orbital =  bfn.substr(found);
01647             orbitalNames.push_back(orbital);
01648             // std::cout << orbital << std::endl;
01649           } else {
01650             printf("orcaplugin) Could not determine orbital description.\n");
01651             return FALSE;
01652           }
01653         }
01654         // reading coefficients
01655         std::vector<float> currentMoCoeffs;
01656         for (size_t cidx = 0; cidx < coefficientNumber; cidx++) {
01657           currentMoCoeffs.push_back(coeff[cidx]);
01658           // std::cout << coeff[cidx] << std::endl;
01659         }
01660         moCoefficients.push_back(currentMoCoeffs);
01661       } else {
01662         // block seems to be finished
01663         readingBlock = 0;
01664         haveAngMom = 1;
01665       }
01666     }
01667     allCoefficients.push_back(moCoefficients);
01668   }
01669 
01670   if ( std::adjacent_find( numberContractedBf.begin(), numberContractedBf.end(), std::not_equal_to<int>() ) != numberContractedBf.end() ) {
01671     printf("orcaplugin) Molecular orbital section corrupted. Did not read consistent number of contracted basis functions!\n");
01672     return FALSE;
01673   }
01674 
01675   // now loop over MO blocks and convert coefficients if needed
01676   int readingPureFunction = 0;
01677   std::vector<std::vector<float> > pureFunction;
01678   std::vector<std::string> pureFunctionName;
01679   int expectedNumberOfPureFunctions = 0;
01680   //std::vector<std::string> orbList{"s","px","py","pz","dz2","dxz","dyz","dx2y2","dxy",
01681   //"f0","f+1","f-1","f+2","f-2","f+3","f-3"};
01682   const char* orbitalList[] = {"s","px","py","pz","dz2","dxz","dyz","dx2y2","dxy",
01683   "f0","f+1","f-1","f+2","f-2","f+3","f-3"};
01684   std::vector<std::string> orbList(orbitalList, orbitalList+15);
01685   int orbRowIndex = 0;
01686   int blockIdx = 0;
01687 
01688   // for (auto name : orbitalNames) {
01689   //   std::cout << name << std::endl;
01690   // }
01691 
01692   MoCoeff newAllCoefficients;
01693   //MoCoeff is a vector<vector<vecto<float>r>>
01694   //for (auto moBlock : allCoefficients) {
01695   for (size_t i = 0; i < allCoefficients.size(); i++) {
01696     CoeffRowBlock newRows;
01697     int blockNumberOfContracted = 0;
01698     //for (auto moRow : moBlock) {
01699     for (size_t j = 0; j < allCoefficients[i].size(); j++) {
01700       int angX = 0, angY = 0, angZ = 0;
01701       std::string orbital = orbitalNames[orbRowIndex];
01702       //std::begin and std::end are not supported on c++98
01703       //std::vector<std::string>::iterator orbIndex = find(std::begin(orbList), std::end(orbList), orbital);
01704       std::vector<std::string>::iterator orbIndex = find(orbList.begin(), orbList.end(), orbital);
01705       //if (std::end(orbList) != orbIndex) {
01706       if (orbList.end() != orbIndex) {
01707         int listIndex = orbIndex-orbList.begin();
01708         switch (listIndex) {
01709           case 0:
01710             break;
01711           case 1:
01712             angX = 1;
01713             break;
01714           case 2:
01715             angY = 1;
01716             break;
01717           case 3:
01718             angZ = 1;
01719             break;
01720           default:
01721             break;
01722         }
01723         // d-shell found
01724         // std::cout << "orbital index: " << orbIndex-orbList.begin() << " " << orbital << std::endl;
01725         if (listIndex > 3) {
01726           if (!readingPureFunction) {
01727             pureFunction.clear();
01728             pureFunctionName.clear();
01729             if (orbital.compare(0, 1, "d") == 0) {
01730               expectedNumberOfPureFunctions = 5;
01731               // std::cout << "Trying to read d-functions." << std::endl;
01732             } else if (orbital.compare(0, 1, "f") == 0) {
01733               expectedNumberOfPureFunctions = 7;
01734               // std::cout << "Trying to read f-functions." << std::endl;
01735             }
01736             readingPureFunction = 1;
01737            // pureFunction.push_back(moRow);
01738             pureFunction.push_back(allCoefficients[i][j]);
01739             pureFunctionName.push_back(orbital);
01740           } else {
01741             //pureFunction.push_back(moRow);
01742             pureFunction.push_back(allCoefficients[i][j]);
01743             pureFunctionName.push_back(orbital);
01744             if (pureFunction.size() == expectedNumberOfPureFunctions) {
01745               // std::cout << "found complete pure function set." << std::endl;
01746               //CoeffRowBlock is a vector<vector<float>>
01747               CoeffRowBlock newBlock = convertPure(pureFunction);
01748               blockNumberOfContracted+= newBlock.size();
01749               //for (auto r : newBlock) {
01750               for (size_t k = 0; k < newBlock.size(); k++) {
01751                 newRows.push_back(newBlock[k]);
01752               }
01753               if (!blockIdx) {
01754                 for (size_t i = 0; i < newBlock.size(); i++) {
01755                   if (expectedNumberOfPureFunctions == 5) {
01756                     //for (auto angMom : dAngMom[i]) {
01757                     //wfAngMoment.push_back(dAngMom[i][0]);
01758                     //wfAngMoment.push_back(dAngMom[i][1]);
01759                     //wfAngMoment.push_back(dAngMom[i][2]);
01760                     wfAngMoment.push_back(dAngMom[i*3 + 0]);
01761                     wfAngMoment.push_back(dAngMom[i*3 + 1]);
01762                     wfAngMoment.push_back(dAngMom[i*3 + 2]);
01763                     //}
01764                   } else if (expectedNumberOfPureFunctions == 7) {
01765                     //for (auto angMom : fAngMom[i]) {
01766                     //wfAngMoment.push_back(fAngMom[i][0]);
01767                     //wfAngMoment.push_back(fAngMom[i][1]);
01768                     //wfAngMoment.push_back(fAngMom[i][2]);
01769                     wfAngMoment.push_back(fAngMom[i*3 + 0]);
01770                     wfAngMoment.push_back(fAngMom[i*3 + 1]);
01771                     wfAngMoment.push_back(fAngMom[i*3 + 2]);
01772                     //}
01773                   }
01774                 }
01775               }
01776               readingPureFunction = 0;
01777             }
01778           }
01779         } else {
01780           newRows.push_back(allCoefficients[i][j]);
01781           blockNumberOfContracted++;
01782           if (!blockIdx) {
01783             wfAngMoment.push_back(angX);
01784             wfAngMoment.push_back(angY);
01785             wfAngMoment.push_back(angZ);
01786           }
01787         }
01788       } else {
01789         std::cout << "orcaplugin) ERROR. Only s/p/d/f-shells are supported." << std::endl;
01790         return FALSE;
01791       }
01792 
01793       orbRowIndex++;
01794     }
01795     numberContractedBf.push_back(blockNumberOfContracted);
01796     newAllCoefficients.push_back(newRows);
01797     orbRowIndex = 0;
01798     blockIdx++;
01799   }
01800 
01801 
01802   // assign the number of contracted functions to wavefunction size
01803   data->wavef_size = numberContractedBf[0];
01804   wf->num_orbitals  = num_orbitals;
01805   wf->orb_energies = (float *) calloc(num_orbitals, sizeof(float));
01806   wf->orb_occupancies = (float *) calloc(num_orbitals, sizeof(float));
01807   wf->wave_coeffs = (float *) calloc(num_orbitals * data->wavef_size, sizeof(float));
01808   wf->has_occup = TRUE;
01809   wf->has_orben = TRUE;
01810 
01811   int cnt = 0;
01812   
01813   //orbitalEnergies is a std::vector<float>
01814   //Ranges do not work...
01815   //for (auto en : orbitalEnergies) {
01816   for (int i = 0; i < orbitalEnergies.size(); i++) {
01817     //wf->orb_energies[cnt] = en;
01818     //cnt++;
01819     wf->orb_energies[i] = orbitalEnergies[i];
01820   }
01821   cnt = 0;
01822   //for (auto occ : orbitalOccupancies) {
01823   for (int i = 0; i < orbitalOccupancies.size();i++) {
01824     //wf->orb_occupancies[cnt] = occ;
01825     //cnt++;
01826     wf->orb_occupancies[i] = orbitalOccupancies[i];
01827   }
01828 
01829   int rowIndex = 0, columnIndex = 0;
01830   blockIdx = 0;
01831   int moBlockSize = 0;
01832   int moBlockIdx = 0;
01833   //for (auto moBlock : newAllCoefficients) {
01834   for (int i = 0; i < newAllCoefficients.size();i++) {
01835     //for (auto moRow : moBlock) {
01836     for (int j = 0; j < newAllCoefficients[i].size(); j++) {
01837       // std::cout << rowIndex << std::endl;
01838       //for (auto moCo : moRow) {
01839       for (int k = 0; k < newAllCoefficients[i][j].size(); k++) {
01840         if ((columnIndex * data->wavef_size + rowIndex) > num_orbitals * data->wavef_size) {
01841           std::cout << "something went wrong:" << columnIndex << std::endl;
01842           std::cout << "something went wrong:" << (columnIndex * data->wavef_size + rowIndex) << " vs. " << num_orbitals * data->wavef_size << std::endl;
01843           return FALSE;
01844         }
01845         // std::cout << orbitalNames[rowIndex] << std::endl;
01846         //wf->wave_coeffs[columnIndex * data->wavef_size + rowIndex] = moCo;
01847         wf->wave_coeffs[columnIndex * data->wavef_size + rowIndex] = newAllCoefficients[i][j][k];
01848         columnIndex++;
01849       }
01850       columnIndex = moBlockSize;
01851       rowIndex++;
01852     }
01853     rowIndex = 0;
01854     // 0-based!!!
01855     //moBlockSize += moBlock[moBlockIdx].size();
01856     moBlockSize += newAllCoefficients[i][moBlockIdx].size();
01857     columnIndex = moBlockSize;
01858     moBlockIdx++;
01859     // std::cout << "bs: " << moBlockSize << std::endl;
01860   }
01861 
01862   // LOGGING for MO coefficients
01863   /*
01864   float coeff2 = 0;
01865   for (size_t t = 0; t < (num_orbitals * data->wavef_size); t++) {
01866     if (t % data->wavef_size == 0) {
01867       std::cout << "---------- " << t/num_orbitals << " c2: " << coeff2 << std::endl;
01868       coeff2 = 0;
01869     }
01870     coeff2 += wf->wave_coeffs[t]*wf->wave_coeffs[t];
01871     std::cout << wf->wave_coeffs[t] << std::endl;
01872   }
01873   */
01874 
01875   if (data->num_frames_read < 1) {
01876     data->angular_momentum = (int*)calloc(wfAngMoment.size(), sizeof(int));
01877   }
01878 
01879   // std::cout << "wfang: " << wfAngMoment.size() <<  " " << 3*data->wavef_size <<std::endl;
01880   for (size_t ang = 0; ang < wfAngMoment.size(); ang++) {
01881     data->angular_momentum[ang] = wfAngMoment[ang];
01882   }
01883 
01884   // TODO: This is just a workaround and might give wrong
01885   // results when reading unrestricted jobs
01886   data->num_occupied_A = occupiedOrbitals;
01887   data->num_occupied_B = occupiedOrbitals;
01888   // data->num_electrons = numberOfElectrons;
01889 
01890   // if (data->num_electrons != numberOfElectrons) {
01891   //
01892   // }
01893 
01894   std::cout << "----------------------------------------" << std::endl;
01895   std::cout << "Total number of orbitals: " << num_orbitals << std::endl;
01896   std::cout << "Number of electrons in read wf: " << numberOfElectrons<< std::endl;
01897   std::cout << "Number of occupied orbitals: " << occupiedOrbitals << std::endl;
01898   std::cout << "Number of contracted bf: " << numberContractedBf[0] << std::endl;
01899   std::cout << "----------------------------------------" << std::endl;
01900 
01901   return TRUE;
01902 }
01903 
01904 /******************************************************
01905  *
01906  * Populate the flat arrays containing the basis
01907  * set data.
01908  *
01909  ******************************************************/
01910 static int fill_basis_arrays(qmdata_t *data) {
01911   //orcadata *orca = (orcadata *)data->format_specific_data;
01912   int i, j, k;
01913   int shellcount = 0;
01914   int primcount = 0;
01915 
01916   float *basis;
01917   int *num_shells_per_atom;
01918   int *num_prim_per_shell;
01919   int *shell_types;
01920   int *atomicnum_per_basisatom;
01921 
01922   /* Count the total number of primitives which
01923    * determines the size of the basis array. */
01924   for (i=0; i<data->num_basis_atoms; i++) {
01925     for (j=0; j<data->basis_set[i].numshells; j++) {
01926       primcount += data->basis_set[i].shell[j].numprims;
01927     }
01928   }
01929 
01930   /* reserve space for pointer to array containing basis
01931    * info, i.e. contraction coeficients and expansion
01932    * coefficients; need 2 entries per basis function, i.e.
01933    * exponent and contraction coefficient; also,
01934    * allocate space for the array holding the orbital symmetry
01935    * information per primitive Gaussian.
01936    * Finally, initialize the arrays holding the number of
01937    * shells per atom and the number of primitives per shell*/
01938   basis = (float *)calloc(2*primcount,sizeof(float));
01939 
01940   /* make sure memory was allocated properly */
01941   if (basis == NULL) {
01942     PRINTERR;
01943     return FALSE;
01944   }
01945 
01946   shell_types = (int *)calloc(data->num_shells, sizeof(int));
01947 
01948   /* make sure memory was allocated properly */
01949   if (shell_types == NULL) {
01950     PRINTERR;
01951     return FALSE;
01952   }
01953 
01954   num_shells_per_atom = (int *)calloc(data->num_basis_atoms, sizeof(int));
01955 
01956   /* make sure memory was allocated properly */
01957   if (num_shells_per_atom == NULL) {
01958     PRINTERR;
01959     return FALSE;
01960   }
01961 
01962   num_prim_per_shell = (int *)calloc(data->num_shells, sizeof(int));
01963 
01964   /* make sure memory was allocated properly */
01965   if (num_prim_per_shell == NULL) {
01966     PRINTERR;
01967     return FALSE;
01968   }
01969 
01970   atomicnum_per_basisatom = (int *)calloc(data->num_basis_atoms, sizeof(int));
01971 
01972   /* make sure memory was allocated properly */
01973   if (atomicnum_per_basisatom == NULL) {
01974     PRINTERR;
01975     return FALSE;
01976   }
01977 
01978 
01979   /* store pointers in struct qmdata_t */
01980   data->basis = basis;
01981   data->shell_types = shell_types;
01982   data->num_shells_per_atom = num_shells_per_atom;
01983   data->num_prim_per_shell = num_prim_per_shell;
01984   data->atomicnum_per_basisatom = atomicnum_per_basisatom;
01985 
01986   /* Go through all basis set atoms and try to assign the
01987    * atomic numbers. The basis set atoms are specified by
01988    * name strings (the same as in the coordinate section,
01989    * except for FMO calcs.) and we try to match the names
01990    * from the two lists. The basis set atom list is symmetry
01991    * unique while the coordinate atom list is complete.*/
01992   primcount = 0;
01993   for (i=0; i<data->num_basis_atoms; i++) {
01994     int found = 0;
01995 
01996     /* For this basis atom find a matching atom from the
01997      * coordinate atom list. */
01998     for (j=0; j<data->numatoms; j++) {
01999       char basisname[BUFSIZ];
02000       strcpy(basisname, data->basis_set[i].name);
02001 
02002       /* for FMO calculations we have to strip the "-n" tail
02003        * of the basis atom name. */
02004       // if (gms->have_fmo) {
02005       //   *strchr(basisname, '-') = '\0';
02006       // }
02007 
02008       if (!strcmp(data->atoms[j].type, basisname)) {
02009         found = 1;
02010         break;
02011       }
02012     }
02013     if (!found) {
02014       printf("orcaplugin) WARNING: Couldn't find atomic number for basis set atom %s\n",
02015              data->basis_set[i].name);
02016       data->basis_set[i].atomicnum = 0;
02017       atomicnum_per_basisatom[i] = 0;
02018     } else {
02019       /* assign atomic number */
02020       data->basis_set[i].atomicnum = data->atoms[j].atomicnum;
02021       atomicnum_per_basisatom[i]   = data->atoms[j].atomicnum;
02022     }
02023     num_shells_per_atom[i] = data->basis_set[i].numshells;
02024 
02025     for (j=0; j<data->basis_set[i].numshells; j++) {
02026       shell_types[shellcount] = data->basis_set[i].shell[j].type;
02027       num_prim_per_shell[shellcount] = data->basis_set[i].shell[j].numprims;
02028 
02029       for (k=0; k<data->basis_set[i].shell[j].numprims; k++) {
02030         basis[2*primcount  ] = data->basis_set[i].shell[j].prim[k].exponent;
02031         basis[2*primcount+1] = data->basis_set[i].shell[j].prim[k].contraction_coeff;
02032         primcount++;
02033       }
02034       shellcount++;
02035     }
02036   }
02037   printf("orcaplugin) Filled basis arrays.\n");
02038 
02039   return TRUE;
02040 }
02041 
02042 
02043 /**************************************************
02044  *
02045  * Convert shell type from char to int.
02046  *
02047  ************************************************ */
02048 static int shelltype_int(char type) {
02049   int shelltype;
02050 
02051   switch (type) {
02052     case 'L':
02053       /* SP_P shells are assigned in get_basis() */
02054       shelltype = SP_S_SHELL;
02055       break;
02056     case 'S':
02057       shelltype = S_SHELL;
02058       break;
02059     case 'P':
02060       shelltype = P_SHELL;
02061       break;
02062     case 'D':
02063       shelltype = D_SHELL;
02064       break;
02065     case 'F':
02066       shelltype = F_SHELL;
02067       break;
02068     case 'G':
02069       shelltype = G_SHELL;
02070       break;
02071     default:
02072       shelltype = UNK_SHELL;
02073       break;
02074   }
02075 
02076   return shelltype;
02077 }
02078 
02079 /* Analyze the trajectory.
02080  * Read the parameters controlling geometry search and
02081  * find the end of the trajectory, couinting the frames
02082  * on the way. Store the filepointer for the beginning of
02083  * each frame in *filepos_array. */
02084 static int analyze_traj(qmdata_t *data, orcadata *orca) {
02085   char buffer[BUFSIZ];
02086   char *line;
02087   long filepos;
02088   filepos = ftell(data->file);
02089 
02090   data->filepos_array = (long* )calloc(1, sizeof(long ));
02091 
02092   /* currently, only one frame is supported!
02093    * lines 3130-3348 in gamessplugin.c
02094    */
02095   if (data->runtype == MOLFILE_RUNTYPE_ENERGY) {
02096     /* We have just one frame */
02097     data->num_frames = 1;
02098     pass_keyline(data->file, "Single Point Calculation", NULL);
02099     data->filepos_array[0] = ftell(data->file);
02100 
02101     /* Check wether SCF has converged */
02102     // if (pass_keyline(data->file,
02103     //                  "SCF IS UNCONVERGED, TOO MANY ITERATIONS",
02104     //                  "ENERGY COMPONENTS")==FOUND) {
02105     //   printf("orcaplugin) SCF IS UNCONVERGED, TOO MANY ITERATIONS\n");
02106     //   data->status = MOLFILE_QMSTATUS_SCF_NOT_CONV;
02107     // } else {
02108     //   data->status = MOLFILE_QMSTATUS_OPT_CONV;
02109     //   fseek(data->file, data->filepos_array[0], SEEK_SET);
02110     // }
02111 
02112     pass_keyline(data->file, "FINAL SINGLE POINT ENERGY", NULL);
02113     data->end_of_traj = ftell(data->file);
02114 
02115     /* Allocate memory for the frame */
02116     data->qm_timestep = (qm_timestep_t *)calloc(1, sizeof(qm_timestep_t));
02117     memset(data->qm_timestep, 0, sizeof(qm_timestep_t));
02118 
02119     return TRUE;
02120   } else if (data->runtype == MOLFILE_RUNTYPE_GRADIENT) {
02121     int appendedCalculations = 0;
02122     rewind(data->file);
02123     pass_keyline(data->file, "Energy+Gradient Calculation", NULL);
02124     data->filepos_array[0] = ftell(data->file);
02125     data->num_frames = 1;
02126 
02127     while(TRUE) {
02128       if (!fgets(buffer, sizeof(buffer), data->file)) break;
02129       line = trimleft(buffer);
02130 
02131       std::string l(line);
02132       if (l.find("Energy+Gradient Calculation") != std::string::npos && data->runtype==MOLFILE_RUNTYPE_GRADIENT) {
02133         appendedCalculations++;
02134         // std::cout << l << std::endl;
02135         if (data->num_frames > 0) {
02136           data->filepos_array = (long*)realloc(data->filepos_array, (data->num_frames+1)*sizeof(long));
02137         }
02138         data->filepos_array[data->num_frames] = ftell(data->file);
02139         data->num_frames++;
02140       }
02141     }
02142 
02143     if (appendedCalculations) {
02144       std::cout << "orcaplugin) Found multiple appended gradient calculations: " << data->num_frames << std::endl;
02145       pass_keyline(data->file, "FINAL SINGLE POINT ENERGY", NULL);
02146       data->end_of_traj = ftell(data->file);
02147       fseek(data->file, filepos, SEEK_SET);
02148 
02149       data->qm_timestep = (qm_timestep_t *)calloc(data->num_frames,
02150                                                   sizeof(qm_timestep_t));
02151       memset(data->qm_timestep, 0, data->num_frames*sizeof(qm_timestep_t));
02152     } else {
02153       data->num_frames = 1;
02154       pass_keyline(data->file, "FINAL SINGLE POINT ENERGY", NULL);
02155       data->end_of_traj = ftell(data->file);
02156 
02157       /* Allocate memory for the frame */
02158       data->qm_timestep = (qm_timestep_t *)calloc(data->num_frames, sizeof(qm_timestep_t));
02159       memset(data->qm_timestep, 0, sizeof(qm_timestep_t));
02160     }
02161     return TRUE;
02162   }
02163   else if (data->runtype == MOLFILE_RUNTYPE_OPTIMIZE) {
02164     std::cout << "orcaplugin) Reading trajectory of optimization." << std::endl;
02165     rewind(data->file);
02166     goto_keyline(data->file, "Geometry Optimization Run", NULL);
02167   }
02168   else {
02169     std::cout << "orcaplugin) Jobtype not supported for trajectory reading." << std::endl;
02170     return FALSE;
02171   }
02172 
02173   // printf("orcaplugin) Analyzing trajectory...\n");
02174   data->status = MOLFILE_QMSTATUS_UNKNOWN;
02175 
02176   while (TRUE) {
02177     if (!fgets(buffer, sizeof(buffer), data->file)) break;
02178     line = trimleft(buffer);
02179 
02180     std::string l(line);
02181     if (l.find("GEOMETRY OPTIMIZATION CYCLE") != std::string::npos && data->runtype==MOLFILE_RUNTYPE_OPTIMIZE) {
02182       // std::cout << l << std::endl;
02183       if (data->num_frames > 0) {
02184         data->filepos_array = (long*)realloc(data->filepos_array, (data->num_frames+1)*sizeof(long));
02185       }
02186       data->filepos_array[data->num_frames] = ftell(data->file);
02187       data->num_frames++;
02188     }
02189     else if (l.find("THE OPTIMIZATION HAS CONVERGED") != std::string::npos) {
02190       printf("orcaplugin) ==== End of trajectory (%d frames) ====\n", data->num_frames);
02191       data->status = MOLFILE_QMSTATUS_OPT_CONV;
02192     }
02193     else if (data->status == MOLFILE_QMSTATUS_OPT_CONV) {
02194       if(l.find("FINAL ENERGY EVALUATION AT THE STATIONARY POINT") != std::string::npos) {
02195         if (data->num_frames > 0) {
02196           data->filepos_array = (long*)realloc(data->filepos_array, (data->num_frames+1)*sizeof(long));
02197         }
02198         std::cout << "orcaplugin) found equilibrium geometry." << std::endl;
02199         data->filepos_array[data->num_frames] = ftell(data->file);
02200         data->num_frames++;
02201         goto_keyline(data->file, "TOTAL RUN TIME", NULL);
02202         break;
02203       }
02204     }
02205   }
02206 
02207   data->end_of_traj = ftell(data->file);
02208   fseek(data->file, filepos, SEEK_SET);
02209 
02210   if (data->status == MOLFILE_QMSTATUS_UNKNOWN) {
02211     /* We didn't find any of the regular key strings,
02212      * the run was most likely broken off and we have an
02213      * incomplete file. */
02214     data->status = MOLFILE_QMSTATUS_FILE_TRUNCATED;
02215   }
02216 
02217 
02218   /* Allocate memory for all frames */
02219   data->qm_timestep = (qm_timestep_t *)calloc(data->num_frames,
02220                                               sizeof(qm_timestep_t));
02221   memset(data->qm_timestep, 0, data->num_frames*sizeof(qm_timestep_t));
02222 
02223 
02224   if (data->status == MOLFILE_QMSTATUS_SCF_NOT_CONV ||
02225       data->status == MOLFILE_QMSTATUS_FILE_TRUNCATED) {
02226     return FALSE;
02227   }
02228 
02229   return TRUE;
02230 }
02231 
02232 
02233 
02234 /***********************************************************
02235  * Provide QM metadata for next timestep.
02236  * This actually triggers reading the entire next timestep
02237  * since we have to parse the whole timestep anyway in order
02238  * to get the metadata. So we store the read data locally
02239  * and hand them to VMD when requested by read_timestep().
02240  *
02241  ***********************************************************/
02242 static int read_qm_timestep_metadata(void *mydata,
02243                                     molfile_qm_timestep_metadata_t *meta) {
02244   int have = 0;
02245 
02246   qmdata_t *data = (qmdata_t *)mydata;
02247 
02248   meta->count = -1; /* Don't know the number of frames yet */
02249 
02250   if (data->num_frames_read > data->num_frames_sent) {
02251     have = 1;
02252   }
02253   else if (data->num_frames_read < data->num_frames) {
02254     printf("orcaplugin) Probing timestep %d\n", data->num_frames_read);
02255 
02256     have = get_traj_frame(data, data->atoms, data->numatoms);
02257   }
02258 
02259   if (have) {
02260     // std::cout << "have frame" << std::endl;
02261     int i;
02262     qm_timestep_t *cur_ts;
02263 
02264     /* get a pointer to the current qm timestep */
02265     cur_ts = data->qm_timestep+data->num_frames_sent;
02266 
02267     // std::cout << "numwave: " << cur_ts->numwave << std::endl;
02268 
02269     for (i=0; (i<MOLFILE_MAXWAVEPERTS && i<cur_ts->numwave); i++) {
02270       meta->num_orbitals_per_wavef[i] = cur_ts->wave[i].num_orbitals;
02271       meta->has_occup_per_wavef[i]    = cur_ts->wave[i].has_occup;
02272       meta->has_orben_per_wavef[i]    = cur_ts->wave[i].has_orben;
02273       // std::cout << "occ: " << cur_ts->wave[i].has_occup << std::endl;
02274       // std::cout << "energy: " << cur_ts->wave[i].has_orben << std::endl;
02275     }
02276     meta->wavef_size      = data->wavef_size;
02277     meta->num_wavef       = cur_ts->numwave;
02278     meta->num_scfiter     = cur_ts->num_scfiter;
02279     meta->num_charge_sets = cur_ts->have_mulliken +
02280       cur_ts->have_lowdin + cur_ts->have_esp;
02281     if (cur_ts->gradient) meta->has_gradient = TRUE;
02282 
02283   } else {
02284     // std::cout << "not have frame" << std::endl;
02285     meta->has_gradient = FALSE;
02286     meta->num_scfiter  = 0;
02287     meta->num_orbitals_per_wavef[0] = 0;
02288     meta->has_occup_per_wavef[0] = FALSE;
02289     meta->num_wavef = 0;
02290     meta->wavef_size = 0;
02291     meta->num_charge_sets = 0;
02292     data->trajectory_done = TRUE;
02293   }
02294 
02295   return MOLFILE_SUCCESS;
02296 }
02297 
02298 
02299 
02300 
02301 /*************************************************************
02302  *
02303  * plugin registration
02304  *
02305  **************************************************************/
02306 static molfile_plugin_t plugin;
02307 
02308 VMDPLUGIN_API int VMDPLUGIN_init(void) {
02309   memset(&plugin, 0, sizeof(molfile_plugin_t));
02310   plugin.abiversion = vmdplugin_ABIVERSION;
02311   plugin.type = MOLFILE_PLUGIN_TYPE;
02312   plugin.name = "orca";
02313   plugin.prettyname = "Orca";
02314   plugin.author = "Maximilian Scheurer, Michael F. Herbst, Marcelo Melo, Julio Maia, John STone";
02315   plugin.majorv = 0;
02316   plugin.minorv = 1;
02317   plugin.is_reentrant = VMDPLUGIN_THREADUNSAFE;
02318   plugin.filename_extension = "orca";
02319   plugin.open_file_read = open_orca_read;
02320   plugin.read_structure = read_orca_structure;
02321   plugin.close_file_read = close_orca_read;
02322   //
02323   plugin.read_qm_metadata = read_orca_metadata;
02324   plugin.read_qm_rundata  = read_orca_rundata;
02325 
02326 #if vmdplugin_ABIVERSION > 11
02327   plugin.read_timestep_metadata    = read_timestep_metadata;
02328   plugin.read_qm_timestep_metadata = read_qm_timestep_metadata;
02329   plugin.read_timestep = read_timestep;
02330 #endif
02331 
02332   return VMDPLUGIN_SUCCESS;
02333 }
02334 
02335 VMDPLUGIN_API int VMDPLUGIN_register(void *v, vmdplugin_register_cb cb) {
02336   (*cb)(v, (vmdplugin_t *)&plugin);
02337   return VMDPLUGIN_SUCCESS;
02338 }
02339 
02340 VMDPLUGIN_API int VMDPLUGIN_fini(void) {
02341   return VMDPLUGIN_SUCCESS;
02342 }
02343 
02344 
02345 #if vmdplugin_ABIVERSION > 11
02346 
02347 /***********************************************************
02348  * Provide non-QM metadata for next timestep.
02349  * Required by the plugin interface.
02350  ***********************************************************/
02351 static int read_timestep_metadata(void *mydata,
02352                                   molfile_timestep_metadata_t *meta) {
02353   meta->count = -1;
02354   meta->has_velocities = 0;
02355 
02356   return MOLFILE_SUCCESS;
02357 }
02358 /***********************************************************
02359  *
02360  * This function provides the data of the next timestep.
02361  * Here we actually don't read the data from file, that had
02362  * to be done already upon calling read_timestep_metadata().
02363  * Instead we copy the stuff from the local data structure
02364  * into the one's provided by VMD.
02365  *
02366  ***********************************************************/
02367 static int read_timestep(void *mydata, int natoms,
02368        molfile_timestep_t *ts, molfile_qm_metadata_t *qm_metadata,
02369        molfile_qm_timestep_t *qm_ts)
02370 {
02371   qmdata_t *data = (qmdata_t *)mydata;
02372   qm_timestep_t *cur_ts;
02373   int offset;
02374   int i = 0;
02375   int num_charge_sets = 0;
02376 
02377   if (data->trajectory_done == TRUE) {
02378     printf("orcaplugin) Trajectory done.\n");
02379     return MOLFILE_ERROR;
02380   }
02381 
02382   /* copy the coordinates */
02383   for (i=0; i<natoms; i++) {
02384     ts->coords[3*i  ] = data->atoms[i].x;
02385     ts->coords[3*i+1] = data->atoms[i].y;
02386     ts->coords[3*i+2] = data->atoms[i].z;
02387     // printf("x: %f y: %f z: %f\n", data->atoms[i].x, data->atoms[i].y, data->atoms[i].z);
02388   }
02389 
02390   /* get a convenient pointer to the current qm timestep */
02391   cur_ts = data->qm_timestep+data->num_frames_sent;
02392   //
02393   // /* store the SCF energies */
02394   // for (i=0; i<cur_ts->num_scfiter; i++) {
02395   //   qm_ts->scfenergies[i] = cur_ts->scfenergies[i];
02396   // }
02397   //
02398   // /* store gradients */
02399   // if (cur_ts->gradient) {
02400   //   for (i=0; i<3*natoms; i++) {
02401   //     qm_ts->gradient[i] = cur_ts->gradient[i];
02402   //   }
02403   // }
02404   //
02405   // /* store charge sets*/
02406   if (cur_ts->have_mulliken) {
02407     offset = num_charge_sets*data->numatoms;
02408     for (i=0; i<data->numatoms; i++) {
02409       qm_ts->charges[offset+i] = cur_ts->mulliken_charges[i];
02410     }
02411     qm_ts->charge_types[num_charge_sets] = MOLFILE_QMCHARGE_MULLIKEN;
02412     num_charge_sets++;
02413   }
02414   //
02415   // if (cur_ts->have_lowdin) {
02416   //   offset = num_charge_sets*data->numatoms;
02417   //   for (i=0; i<data->numatoms; i++) {
02418   //     qm_ts->charges[offset+i] = cur_ts->lowdin_charges[i];
02419   //   }
02420   //   qm_ts->charge_types[num_charge_sets] = MOLFILE_QMCHARGE_LOWDIN;
02421   //   num_charge_sets++;
02422   // }
02423   // if (cur_ts->have_esp) {
02424   //   offset = num_charge_sets*data->numatoms;
02425   //   for (i=0; i<data->numatoms; i++) {
02426   //     qm_ts->charges[offset+i] = cur_ts->esp_charges[i];
02427   //   }
02428   //   qm_ts->charge_types[num_charge_sets] = MOLFILE_QMCHARGE_ESP;
02429   //   num_charge_sets++;
02430   // }
02431   //
02432   //
02433   /* store the wave function and orbital energies */
02434   if (cur_ts->wave) {
02435     std::cout << "orcaplugin) Have wavefunctions: " << cur_ts->numwave << " in frame: " << data->num_frames_sent << std::endl;
02436     for (i=0; i<cur_ts->numwave; i++) {
02437       qm_wavefunction_t *wave = &cur_ts->wave[i];
02438       qm_ts->wave[i].type         = wave->type;
02439       qm_ts->wave[i].spin         = wave->spin;
02440       qm_ts->wave[i].excitation   = wave->exci;
02441       qm_ts->wave[i].multiplicity = wave->mult;
02442       qm_ts->wave[i].energy       = wave->energy;
02443       strncpy(qm_ts->wave[i].info, wave->info, MOLFILE_BUFSIZ);
02444 
02445       if (wave->wave_coeffs) {
02446         memcpy(qm_ts->wave[i].wave_coeffs, wave->wave_coeffs, wave->num_orbitals*data->wavef_size*sizeof(float));
02447       }
02448       if (wave->orb_energies) {
02449         memcpy(qm_ts->wave[i].orbital_energies, wave->orb_energies,
02450                wave->num_orbitals*sizeof(float));
02451       }
02452       if (wave->has_occup) {
02453         memcpy(qm_ts->wave[i].occupancies, wave->orb_occupancies,
02454                wave->num_orbitals*sizeof(float));
02455       }
02456     }
02457   }
02458   //
02459   if (data->runtype == MOLFILE_RUNTYPE_ENERGY ||
02460       data->runtype == MOLFILE_RUNTYPE_HESSIAN) {
02461     /* We have only a single point */
02462     data->trajectory_done = TRUE;
02463   }
02464 
02465   data->num_frames_sent++;
02466 
02467   return MOLFILE_SUCCESS;
02468 }
02469 #endif
02470 
02471 /*****************************************************
02472  *
02473  * provide VMD with the sizes of the QM related
02474  * data structure arrays that need to be made
02475  * available
02476  *
02477  *****************************************************/
02478 static int read_orca_metadata(void *mydata,
02479     molfile_qm_metadata_t *metadata) {
02480 
02481   qmdata_t *data = (qmdata_t *)mydata;
02482 
02483   if (data->runtype == MOLFILE_RUNTYPE_HESSIAN) {
02484     metadata->ncart = (3*data->numatoms);
02485     metadata->nimag = data->nimag;
02486 
02487     if (data->have_internals) {
02488       metadata->nintcoords = data->nintcoords;
02489     } else {
02490       metadata->nintcoords = 0;
02491     }
02492   }
02493   else {
02494     metadata->ncart = 0;
02495     metadata->nimag = 0;
02496     metadata->nintcoords = 0;
02497   }
02498 
02499   /* orbital data */
02500   metadata->num_basis_funcs = data->num_basis_funcs;
02501   metadata->num_basis_atoms = data->num_basis_atoms;
02502   metadata->num_shells      = data->num_shells;
02503   metadata->wavef_size      = data->wavef_size;
02504 
02505 #if vmdplugin_ABIVERSION > 11
02506   /* system and run info */
02507   metadata->have_sysinfo = 1;
02508 
02509   /* hessian info */
02510   metadata->have_carthessian = data->have_cart_hessian;
02511   metadata->have_inthessian  = data->have_int_hessian;
02512 
02513   /* normal mode info */
02514   metadata->have_normalmodes = data->have_normal_modes;
02515 #endif
02516 
02517   return MOLFILE_SUCCESS;
02518 }
02519 
02520 
02521 /******************************************************
02522  *
02523  * Provide VMD with the static (i.e. non-trajectory)
02524  * data. That means we are filling the molfile_plugin
02525  * data structures.
02526  *
02527  ******************************************************/
02528 static int read_orca_rundata(void *mydata,
02529                                molfile_qm_t *qm_data) {
02530 
02531   qmdata_t *data = (qmdata_t *)mydata;
02532   int i, j;
02533   int ncart;
02534   molfile_qm_hessian_t *hessian_data = &qm_data->hess;
02535   molfile_qm_basis_t   *basis_data   = &qm_data->basis;
02536   molfile_qm_sysinfo_t *sys_data     = &qm_data->run;
02537 
02538   /* fill in molfile_qm_hessian_t */
02539   if (data->runtype == MOLFILE_RUNTYPE_HESSIAN) {
02540     ncart = 3*data->numatoms;
02541 
02542     /* Hessian matrix in cartesian coordinates */
02543     if (data->have_cart_hessian) {
02544       for (i=0; i<ncart; i++) {
02545         for (j=0; j<=i; j++) {
02546           hessian_data->carthessian[ncart*i+j] = data->carthessian[ncart*i+j];
02547           hessian_data->carthessian[ncart*j+i] = data->carthessian[ncart*i+j];
02548         }
02549       }
02550     }
02551 
02552     /* Hessian matrix in internal coordinates */
02553     if (data->have_int_hessian) {
02554       for (i=0; i<(data->nintcoords)*(data->nintcoords); i++) {
02555         hessian_data->inthessian[i] = data->inthessian[i];
02556       }
02557     }
02558 
02559     /* wavenumbers, intensities, normal modes */
02560     if (data->have_normal_modes) {
02561       for (i=0; i<ncart*ncart; i++) {
02562         hessian_data->normalmodes[i] = data->normal_modes[i];
02563       }
02564       for (i=0; i<ncart; i++) {
02565         hessian_data->wavenumbers[i] = data->wavenumbers[i];
02566         hessian_data->intensities[i] = data->intensities[i];
02567       }
02568     }
02569 
02570     /* imaginary modes */
02571     for (i=0; i<data->nimag; i++) {
02572       /*printf("imag_modes[%d]=%d\n", i, data->imag_modes[i]);*/
02573       hessian_data->imag_modes[i] = data->imag_modes[i];
02574     }
02575   }
02576 
02577   /* fill in molfile_qm_sysinfo_t */
02578   sys_data->runtype = data->runtype;
02579   sys_data->scftype = data->scftype;
02580   sys_data->nproc   = data->nproc;
02581   sys_data->num_electrons  = data->num_electrons;
02582   sys_data->totalcharge    = data->totalcharge;
02583   sys_data->num_occupied_A = data->num_occupied_A;
02584   sys_data->num_occupied_B = data->num_occupied_B;
02585   sys_data->status         = data->status;
02586 
02587 
02588   strncpy(sys_data->basis_string, data->basis_string,
02589           sizeof(sys_data->basis_string));
02590 
02591   sys_data->memory = 0; /* XXX fixme */
02592 
02593   strncpy(sys_data->runtitle, data->runtitle, sizeof(sys_data->runtitle));
02594   strncpy(sys_data->geometry, data->geometry, sizeof(sys_data->geometry));
02595   strncpy(sys_data->version_string, data->version_string,
02596           sizeof(sys_data->version_string));
02597 
02598 #if vmdplugin_ABIVERSION > 11
02599   /* fill in molfile_qm_basis_t */
02600   if (data->num_basis_funcs) {
02601     for (i=0; i<data->num_basis_atoms; i++) {
02602       basis_data->num_shells_per_atom[i] = data->num_shells_per_atom[i];
02603       basis_data->atomic_number[i] = data->atomicnum_per_basisatom[i];
02604     }
02605 
02606     for (i=0; i<data->num_shells; i++) {
02607       basis_data->num_prim_per_shell[i] = data->num_prim_per_shell[i];
02608       basis_data->shell_types[i] = data->shell_types[i];
02609     }
02610 
02611     for (i=0; i<2*data->num_basis_funcs; i++) {
02612       basis_data->basis[i] = data->basis[i];
02613     }
02614 
02615     for (i=0; i<3*data->wavef_size; i++) {
02616       basis_data->angular_momentum[i] = data->angular_momentum[i];
02617     }
02618     // std::cout << "free data->angular_momentum" << std::endl;
02619     // free(data->angular_momentum);
02620   }
02621 #endif
02622 
02623   return MOLFILE_SUCCESS;
02624 }
02625 
02626 
02627 
02628 
02629 /**********************************************************
02630  *
02631  * close file and free memory
02632  *
02633  **********************************************************/
02634  // move to top later...
02635 static void close_orca_read(void *mydata) {
02636   printf("Freeing memory.\n");
02637 
02638   qmdata_t *data = (qmdata_t *)mydata;
02639   int i, j;
02640   fclose(data->file);
02641 
02642   free(data->atoms);
02643   free(data->basis);
02644   free(data->shell_types);
02645   free(data->atomicnum_per_basisatom);
02646   free(data->num_shells_per_atom);
02647   free(data->num_prim_per_shell);
02648   free(data->bonds);
02649   free(data->angles);
02650   free(data->dihedrals);
02651   free(data->impropers);
02652   free(data->internal_coordinates);
02653   free(data->bond_force_const);
02654   free(data->angle_force_const);
02655   free(data->dihedral_force_const);
02656   free(data->improper_force_const);
02657   free(data->inthessian);
02658   free(data->carthessian);
02659   free(data->wavenumbers);
02660   free(data->intensities);
02661   free(data->normal_modes);
02662   free(data->imag_modes);
02663   free(data->angular_momentum);
02664   data->angular_momentum = NULL;
02665   free(data->filepos_array);
02666 
02667   if (data->basis_set) {
02668     for (i=0; i<data->num_basis_atoms; i++) {
02669       // printf("Freeing basis set of atom %d\n", i);
02670       for (j=0; j<data->basis_set[i].numshells; j++) {
02671         // printf("Freeing shell primitives %d\n", j);
02672         // printf("--- Address of prim is %p\n", (void *)data->basis_set[i].shell[j].prim);
02673         free(data->basis_set[i].shell[j].prim);
02674         data->basis_set[i].shell[j].prim = NULL;
02675       }
02676       // printf("- Address of shell is %p\n", (void *)data->basis_set[i].shell);
02677       free(data->basis_set[i].shell);
02678       data->basis_set[i].shell = NULL;
02679       // printf("- Address of shell is %p\n", (void *)data->basis_set[i].shell);
02680     }
02681     free(data->basis_set);
02682     data->basis_set = NULL;
02683   }
02684 
02685   for (i=0; i<data->num_frames; i++) {
02686     free(data->qm_timestep[i].scfenergies);
02687     free(data->qm_timestep[i].gradient);
02688     free(data->qm_timestep[i].mulliken_charges);
02689     free(data->qm_timestep[i].lowdin_charges);
02690     free(data->qm_timestep[i].esp_charges);
02691     for (j=0; j<data->qm_timestep[i].numwave; j++) {
02692       free(data->qm_timestep[i].wave[j].wave_coeffs);
02693       free(data->qm_timestep[i].wave[j].orb_energies);
02694       free(data->qm_timestep[i].wave[j].orb_occupancies);
02695     }
02696     free(data->qm_timestep[i].wave);
02697   }
02698   free(data->qm_timestep);
02699   free(data->format_specific_data);
02700   free(data);
02701 }
02702 
02703 
02704 
02705 #ifdef TEST_ORCAPLUGIN
02706 
02707 static void print_input_data(qmdata_t *data) {
02708   int i, j, k;
02709   int primcount=0;
02710   int shellcount=0;
02711 
02712   printf("\nDATA READ FROM FILE:\n\n");
02713   printf(" %10s WORDS OF MEMORY AVAILABLE\n", data->memory);
02714   printf("\n");
02715   printf("     BASIS OPTIONS\n");
02716   printf("     -------------\n");
02717   printf("%s\n", data->basis_string);
02718   printf("\n\n\n");
02719   printf("     RUN TITLE\n");
02720   printf("     ---------\n");
02721   printf(" %s\n", data->runtitle);
02722   printf("\n");
02723   printf(" THE POINT GROUP OF THE MOLECULE IS %s\n", "XXX");
02724   printf(" THE ORDER OF THE PRINCIPAL AXIS IS %5i\n", 0);
02725   printf("\n");
02726   printf(" YOUR FULLY SUBSTITUTED Z-MATRIX IS\n");
02727   printf("\n");
02728   printf(" THE MOMENTS OF INERTIA ARE (AMU-ANGSTROM**2)\n");
02729   printf(" IXX=%10.3f   IYY=%10.3f   IZZ=%10.3f\n", 0.0, 0.0, 0.0);
02730   printf("\n");
02731   printf(" ATOM      ATOMIC                      COORDINATES (BOHR)\n");
02732   printf("           CHARGE         X                   Y                   Z\n");
02733   for (i=0; i<data->numatoms; i++) {
02734     printf(" %-8s %6d", data->atoms[i].type, data->atoms[i].atomicnum);
02735 
02736     printf("%17.10f",   ANGS_TO_BOHR*data->atoms[i].x);
02737     printf("%20.10f",   ANGS_TO_BOHR*data->atoms[i].y);
02738     printf("%20.10f\n", ANGS_TO_BOHR*data->atoms[i].z);
02739   }
02740   printf("\n");
02741   printf("     ATOMIC BASIS SET\n");
02742   printf("     ----------------\n");
02743   printf(" THE CONTRACTED PRIMITIVE FUNCTIONS HAVE BEEN UNNORMALIZED\n");
02744   printf(" THE CONTRACTED BASIS FUNCTIONS ARE NOW NORMALIZED TO UNITY\n");
02745   printf("\n");
02746   printf("  SHELL TYPE  PRIMITIVE        EXPONENT          CONTRACTION COEFFICIENT(S)\n");
02747   printf("\n");
02748 
02749 #if 0
02750   for (i=0; i<data->numatoms; i++) {
02751     printf("%-8s\n\n", data->atoms[i].type);
02752     printf("\n");
02753     printf("nshells=%d\n", data->num_shells_per_atom[i]);
02754 
02755     for (j=0; j<data->num_shells_per_atom[i]; j++) {
02756       printf("nprim=%d\n", data->num_prim_per_shell[shellcount]);
02757 
02758       for (k=0; k<data->num_prim_per_shell[shellcount]; k++) {
02759         printf("%6d   %d %7d %22f%22f\n", j, data->shell_types[shellcount],
02760                primcount+1, data->basis[2*primcount], data->basis[2*primcount+1]);
02761         primcount++;
02762       }
02763 
02764       printf("\n");
02765       shellcount++;
02766     }
02767   }
02768 #endif
02769   printf("orcaplugin) =================================================================\n");
02770   for (i=0; i<data->num_basis_atoms; i++) {
02771     printf("%-8s (%10s)\n\n", data->atoms[i].type, data->basis_set[i].name);
02772     printf("\n");
02773 
02774     for (j=0; j<data->basis_set[i].numshells; j++) {
02775 
02776       for (k=0; k<data->basis_set[i].shell[j].numprims; k++) {
02777         printf("%6d   %d %7d %22f%22f\n", j,
02778                data->basis_set[i].shell[j].type,
02779                primcount+1,
02780                data->basis_set[i].shell[j].prim[k].exponent,
02781                data->basis_set[i].shell[j].prim[k].contraction_coeff);
02782         primcount++;
02783       }
02784 
02785       printf("\n");
02786       shellcount++;
02787     }
02788   }
02789   printf("\n");
02790   printf(" TOTAL NUMBER OF BASIS SET SHELLS             =%5d\n", data->num_shells);
02791   printf(" NUMBER OF CARTESIAN GAUSSIAN BASIS FUNCTIONS =%5d\n", data->wavef_size);
02792   printf(" NUMBER OF ELECTRONS                          =%5d\n", data->num_electrons);
02793   printf(" CHARGE OF MOLECULE                           =%5d\n", data->totalcharge);
02794   printf(" SPIN MULTIPLICITY                            =%5d\n", data->multiplicity);
02795   printf(" NUMBER OF OCCUPIED ORBITALS (ALPHA)          =%5d\n", data->num_occupied_A);
02796   printf(" NUMBER OF OCCUPIED ORBITALS (BETA )          =%5d\n", data->num_occupied_B);
02797   printf(" TOTAL NUMBER OF ATOMS                        =%5i\n", data->numatoms);
02798   printf("\n");
02799 }
02800 
02801 #endif
02802 
02803 

Generated on Tue Jan 21 02:56:18 2020 for VMD Plugins (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002