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: mariano $       $Locker:  $             $State: Exp $
00006  *      $Revision: 1.7 $       $Date: 2020/04/13 21:18:28 $
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(1, 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   char versionstr[BUFSIZ];
00674   int mainVersion, secondDigit, thirdDigit;
00675   buffer[0] = '\0';
00676   programLine = goto_keyline(data->file, "O   R   C   A", NULL);
00677   if (programLine != 1) {
00678     return FALSE;
00679   }
00680   strcpy(data->version_string, "ORCA ");
00681   versionLine = goto_keyline(data->file, "Program Version", NULL);
00682   // thisline(data->file);
00683   GET_LINE(buffer, data->file);
00684   if (strstr(buffer,"Version") != NULL) {
00685     sscanf(buffer, "%*s %*s %d.%d.%d", &mainVersion, &secondDigit, &thirdDigit);
00686     #ifdef DEBUGGING_ORCA
00687       printf("DEBUG: build: %d.%d.%d\n", mainVersion, secondDigit, thirdDigit);
00688     #endif
00689     int build[3] = { mainVersion, secondDigit, thirdDigit };
00690     sprintf(versionstr,"%d.%d.%d",mainVersion, secondDigit, thirdDigit);
00691     strcat(data->version_string,versionstr);
00692     for (size_t i = 0; i < 3; i++) {
00693         orca->digits[i] = build[i];
00694     }
00695     switch (mainVersion) {
00696       case ORCA4:
00697         orca->version = 2;
00698         break;
00699       case ORCA3:
00700         orca->version = 1;
00701         break;
00702       default:
00703         orca->version = 0;
00704         break;
00705     }
00706   } else {
00707     PRINTERR;
00708     return FALSE;
00709   }
00710 
00711   return TRUE;
00712 }
00713 
00714 
00715 
00716 static int parse_static_data(qmdata_t *data, int* natoms) {
00717   orcadata *orca = (orcadata *)data->format_specific_data;
00718   if (!get_job_info(data)) return FALSE;
00719 
00720   if (!get_input_structure(data, orca)) return FALSE;
00721 
00722   if (!get_basis(data)) return FALSE;
00723 
00724   if (!analyze_traj(data, orca)) {
00725     printf("orcaplugin) WARNING: Truncated or abnormally terminated file!\n\n");
00726   }
00727 
00728   *natoms = data->numatoms;
00729 
00730   read_first_frame(data);
00731 
00732   // print_input_data(data);
00733 
00734   return TRUE;
00735 }
00736 
00737 
00738 /**********************************************************
00739  *
00740  * Read the input atom definitions and geometry
00741  *
00742  **********************************************************/
00743 static int get_input_structure(qmdata_t *data, orcadata *orca) {
00744   char buffer[BUFSIZ];
00745   int numatoms = -1;
00746   int bohr;
00747   long filepos;
00748   filepos = ftell(data->file);
00749 
00750   if (goto_keyline(data->file, "CARTESIAN COORDINATES (ANGSTROEM)", NULL)) {
00751     GET_LINE(buffer, data->file);
00752     // thisline(data->file);
00753     // UNITS ARE ANGSTROEM
00754     bohr = 0;
00755     // sscanf()
00756   } else {
00757     printf("orcaplugin) No cartesian coordinates in ANGSTROEM found.\n");
00758     return FALSE;
00759   }
00760 
00761   // skip the ---- line
00762   eatline(data->file, 1);
00763   /* Read the coordinate block */
00764   if (get_coordinates(data->file, &data->atoms, bohr, &numatoms))
00765     data->num_frames_read = 0;
00766   else {
00767     printf("orcaplugin) Bad atom coordinate block!\n");
00768     return FALSE;
00769   }
00770 
00771   data->numatoms = numatoms;
00772   return TRUE;
00773 }
00774 
00775 
00776 static int get_basis(qmdata_t *data) {
00777 //  orcadata *orca = (orcadata *)data->format_specific_data; // unused?!?
00778   char buffer[BUFSIZ];
00779   char word[4][BUFSIZ];
00780   int i = 0;
00781   int numread, numshells;
00782   shell_t *shell;
00783   long filepos;
00784   int semiempirical = 0;
00785 
00786   if (!strcmp(data->gbasis, "MNDO") ||
00787       !strcmp(data->gbasis, "AM1")  ||
00788       !strcmp(data->gbasis, "PM3")) {
00789     semiempirical = 1;
00790   }
00791 
00792   /* Search for "ATOMIC BASIS SET" line */
00793   if ((pass_keyline(data->file, "BASIS SET IN INPUT FORMAT", NULL) != FOUND) &&
00794       (pass_keyline(data->file, "GAUSSIAN BASIS SET", NULL) != FOUND)) {
00795     printf("orcaplugin) No basis set found!\n");
00796     return FALSE;
00797   }
00798 
00799   /* initialize buffers */
00800   buffer[0] = '\0';
00801   for (i=0; i<3; i++)
00802     word[i][0] = '\0';
00803 
00804 
00805   /* skip the next 3 lines */
00806   // eatline(data->file, 2);
00807 
00808   /* Allocate space for the basis for all atoms */
00809   data->basis_set = (basis_atom_t*)calloc(data->numatoms, sizeof(basis_atom_t));
00810 
00811   filepos = ftell(data->file);
00812   i = 0; /* basis atom counter */
00813   char elementName[11];
00814   int finished = FALSE;
00815 
00816   basis_atom_t* tempBasis = (basis_atom_t*)calloc(1, sizeof(basis_atom_t));
00817 
00818   prim_t *prim;
00819 
00820   std::vector<std::string> readElements;
00821 #if 0
00822   for (size_t atom = 0; atom < data->numatoms; atom++) {
00823 // WTF?!?!?!
00824   }
00825 #endif
00826 
00827   while (!finished && !semiempirical) {
00828     // printf("Trying to read bf. \n");
00829     if (pass_keyline(data->file, "Basis set for element", NULL) == FOUND ) {
00830       GET_LINE(buffer, data->file);
00831       numread = sscanf(buffer,"%s %s",&word[0][0], &word[1][0]);
00832       strcpy(elementName, &word[1][0]);
00833       printf("New element found: %s\n", &word[1][0]);
00834       std::string atype(&word[1][0]);
00835       if (std::find(readElements.begin(), readElements.end(), atype) != readElements.end()) {
00836         break;
00837       } else {
00838         readElements.push_back(atype);
00839       }
00840       int elementCompleted = 0;
00841 
00842 
00843       shell = (shell_t*)calloc(1, sizeof(shell_t));
00844       numshells = 0;
00845 //      int readingShell = 0;
00846       int primcounter = 0;
00847 
00848       // this is very sloppy at the moment...
00849       // for PM3 etc. Orca prints the bf per atom...
00850 
00851       while(!elementCompleted) {
00852         GET_LINE(buffer, data->file);
00853         numread = sscanf(buffer,"%s %s %s",&word[0][0], &word[1][0],&word[2][0]);
00854         // printf("numread: %d -- %s %s %s \n",numread, &word[0][0], &word[1][0],&word[2][0]);
00855         switch (numread) {
00856           case 1:
00857             if (strcmp(trimleft(trimright(&word[0][0])), "end")) {
00858               // printf("Section ended. \n");
00859               elementCompleted = 1;
00860               break;
00861             }
00862           case 2:
00863             shell[numshells].numprims = atoi(trimleft(trimright(&word[1][0])));
00864             shell[numshells].type = shelltype_int(word[0][0]);
00865             // printf("orcaplugin) Type: %d NPrims: %d\n", shell[numshells].type, shell[numshells].numprims);
00866             primcounter = 0;
00867             prim = (prim_t*)calloc(shell[numshells].numprims, sizeof(prim_t));
00868             // printf("!! Address of prim is %p\n", (void *)prim);
00869             shell[numshells].prim = prim;
00870             numshells++;
00871             if (numshells) {
00872               shell = (shell_t*)realloc(shell, (numshells+1)*sizeof(shell_t));
00873             }
00874             break;
00875           case 3:
00876             prim[primcounter].exponent = atof(&word[1][0]);
00877             prim[primcounter].contraction_coeff = atof(&word[2][0]);
00878             // printf("%f - %f\n", prim[primcounter].exponent, prim[primcounter].contraction_coeff);
00879             primcounter++;
00880             break;
00881           default:
00882             printf("orcaplugin) Unknown line in basis functions.\n");
00883             elementCompleted = 1;
00884             break;
00885         }
00886       }
00887       // printf("Number of shells: %d \n", numshells);
00888       strcpy(tempBasis[i].name, elementName);
00889       tempBasis[i].numshells = numshells;
00890       tempBasis[i].shell = shell;
00891       i++;
00892       if (i) {
00893         tempBasis = (basis_atom_t*)realloc(tempBasis, (i+1)*sizeof(basis_atom_t));
00894       }
00895       // set prim to nullpointer!
00896       prim = NULL;
00897     } else {
00898       finished = TRUE;
00899       printf("orcaplugin) Reading basis set finished! \n");
00900     }
00901   }
00902 
00903   finished = 0;
00904   // semiempirical sto-3g
00905   int ngauss = 3;
00906   int atomCounter = 0;
00907   int shellNumber;
00908   while(!finished && semiempirical) {
00909     if(goto_keyline(data->file,"shells",NULL) == FOUND && atomCounter <(data->numatoms-1)) {
00910       // thisline(data->file);
00911       GET_LINE(buffer, data->file);
00912       std::string lineString(buffer);
00913       std::vector<std::string> elements = split(reduce(lineString), ' ');
00914 
00915       //stoi is not supported
00916       //shellNumber = stoi(elements[2]);
00917       shellNumber = atoi(elements[2].c_str());
00918 
00919       // std::cout << "shell number: " << shellNumber << std::endl;
00920       shell = (shell_t*) calloc(shellNumber, sizeof(shell_t));
00921       data->basis_set[atomCounter].shell = shell;
00922       data->basis_set[atomCounter].numshells = shellNumber;
00923       for (size_t shellIdx = 0; shellIdx < shellNumber; ++shellIdx) {
00924         GET_LINE(buffer, data->file);
00925         prim = (prim_t*) calloc(3, sizeof(prim_t));
00926         shell[shellIdx].prim = prim;
00927         shell[shellIdx].numprims = 3;
00928         shell[shellIdx].type = shellIdx;
00929         for (size_t nbas = 0; nbas < ngauss; ++nbas) {
00930           GET_LINE(buffer, data->file);
00931           std::string l(buffer);
00932           std::vector<std::string> coeff = split(reduce(l), ' ');
00933 #if 1
00934           prim[nbas].exponent = atof(coeff[0].c_str());
00935           prim[nbas].contraction_coeff = atof(coeff[1].c_str());
00936 #else
00937           prim[nbas].exponent = stod(coeff[0]);
00938           prim[nbas].contraction_coeff = stod(coeff[1]);
00939 #endif
00940         }
00941       }
00942 
00943       data->num_shells += shellNumber;
00944       data->num_basis_atoms++;
00945       data->num_basis_funcs += 3;
00946       strncpy(data->basis_set[atomCounter].name, data->atoms[atomCounter].type, 11);
00947       atomCounter++;
00948     } else {
00949       finished = TRUE;
00950       prim = NULL;
00951       std::cout << "orcaplugin) Reading STO-3G basis for semiempirical method finished." << std::endl;
00952       // free unused stuff
00953       free(tempBasis);
00954       // We return here without further ado.
00955       return fill_basis_arrays(data);
00956     }
00957   }
00958 
00959   // As we read GTOs from the Orca output file, we need to
00960   // loop over all atoms and assign the basis functions
00961 
00962   // printf("orcaplugin) Parsed basis set of %d elements. \n", i);
00963   // for (size_t j = 0; j < i; j++) {
00964   //   printf("Element: %s\n", tempBasis[j].name);
00965   //   printf("- NShells: %d\n", tempBasis[j].numshells);
00966   //   for (size_t k = 0; k < tempBasis[j].numshells; k++) {
00967   //     printf("--- NPrims: %d \n", tempBasis[j].shell[k].numprims);
00968   //     for (size_t o = 0; o < tempBasis[j].shell[k].numprims; o++) {
00969   //       float expo = tempBasis[j].shell[k].prim[o].exponent;
00970   //       float cont = tempBasis[j].shell[k].prim[o].contraction_coeff;
00971   //       printf("----- E= %f , C= %f \n", expo, cont);
00972   //     }
00973   //   }
00974   // }
00975 
00976   // Allocate an array of zeros to store whether the tempBasis
00977   // of the same index was actually used.
00978   int *tempBasisUsed = (int*) calloc(i, sizeof(int));
00979 
00980   const char* currentElement;
00981   basis_atom_t* currentBasis;
00982   for (size_t n = 0; n < data->numatoms; n++) {
00983     currentElement = data->atoms[n].type;
00984     for (size_t j = 0; j < i; j++) {
00985       if (!strcmp(currentElement, tempBasis[j].name)) {
00986         // printf("orcaplugin) found basis for element %s\n", currentElement);
00987         currentBasis = &tempBasis[j];
00988         tempBasisUsed[j] = 1;
00989       }
00990     }
00991 
00992     // printf("orcaplugin) Basis for element %s has %d shells.\n", currentElement, currentBasis->numshells);
00993     data->basis_set[n].shell = (shell_t *) calloc(currentBasis->numshells, sizeof(shell_t));
00994     memcpy(data->basis_set[n].shell, currentBasis->shell, currentBasis->numshells * sizeof(shell_t));
00995     data->basis_set[n].numshells = currentBasis->numshells;
00996     data->num_shells += currentBasis->numshells;
00997     for (size_t p = 0; p < currentBasis->numshells; ++p) {
00998       data->basis_set[n].shell[p].prim = (prim_t *) calloc(currentBasis->shell[p].numprims, sizeof(prim_t));
00999       memcpy(data->basis_set[n].shell[p].prim, currentBasis->shell[p].prim, currentBasis->shell[p].numprims * sizeof(prim_t));
01000       data->num_basis_funcs += currentBasis->shell[p].numprims;
01001     }
01002     data->num_basis_atoms++;
01003     strncpy(data->basis_set[n].name, currentBasis->name, 11);
01004     currentBasis = NULL;
01005     currentElement  = NULL;
01006   }
01007 
01008   // cleaning up memory for tempBasis and its shells
01009   for (size_t idx = 0; idx < i; ++idx) {
01010     for (size_t shellIdx = 0; shellIdx < tempBasis[idx].numshells; shellIdx++) {
01011       shell_t* cshell = &tempBasis[idx].shell[shellIdx];
01012       free(cshell->prim);
01013       cshell->prim = NULL;
01014     }
01015     free(tempBasis[idx].shell);
01016     tempBasis[idx].shell = NULL;
01017   }
01018   free(tempBasis);
01019   tempBasis = NULL;
01020   shell = NULL;
01021   printf("orcaplugin) Parsed %d uncontracted basis functions.\n", data->num_basis_funcs);
01022   free(tempBasisUsed);
01023 
01024   // /* allocate and populate flat arrays needed for molfileplugin */
01025   return fill_basis_arrays(data);
01026 }
01027 
01028 
01029 static int get_coordinates(FILE *file, qm_atom_t **atoms, int unit,
01030                            int *numatoms) {
01031   int i = 0;
01032   int growarray = 0;
01033 
01034   if (*numatoms<0) {
01035     *atoms = (qm_atom_t*)calloc(1, sizeof(qm_atom_t));
01036     growarray = 1;
01037   }
01038 
01039   /* Read in the coordinates until an empty line is reached.
01040    * We expect 5 entries per line */
01041   while (1) {
01042     char buffer[BUFSIZ];
01043     char atname[BUFSIZ];
01044     float atomicnum;
01045     float x,y,z;
01046     int n;
01047     qm_atom_t *atm;
01048 
01049     GET_LINE(buffer, file);
01050     // thisline(file);
01051 
01052     /* For FMO there is an additional atom index in the
01053      * second column. Try both variants: */
01054     n = sscanf(buffer,"%s %f %f %f",atname,&x,&y,&z);
01055     // printf("%s\n", atname);
01056     if (n!=4) {
01057       // n = sscanf(buffer,"%s %f %f %f %f",atname,&atomicnum,&x,&y,&z);
01058       break;
01059     }
01060     // if (n!=5 && n!=6) break;
01061 
01062     if (growarray && i>0) {
01063       *atoms = (qm_atom_t*)realloc(*atoms, (i+1)*sizeof(qm_atom_t));
01064     }
01065     atm = (*atoms)+i;
01066 
01067     // just get the atomic number from periodic_table.h
01068     atomicnum = get_pte_idx(atname);
01069 
01070     strncpy(atm->type, atname, sizeof(atm->type));
01071     atm->atomicnum = floor(atomicnum+0.5); /* nuclear charge */
01072     // printf("coor: %s %d %f %f %f\n", atm->type, atm->atomicnum, x, y, z);
01073 
01074     /* if coordinates are in Bohr convert them to Angstrom */
01075     if (unit==BOHR) {
01076       x *= BOHR_TO_ANGS;
01077       y *= BOHR_TO_ANGS;
01078       z *= BOHR_TO_ANGS;
01079     }
01080 
01081     atm->x = x;
01082     atm->y = y;
01083     atm->z = z;
01084     i++;
01085   }
01086 
01087   /* If file is broken off in the middle of the coordinate block
01088    * we cannot use this frame. */
01089   if (*numatoms>=0 && *numatoms!=i) {
01090     (*numatoms) = i;
01091     return FALSE;
01092   }
01093 
01094   (*numatoms) = i;
01095   return TRUE;
01096 }
01097 
01098 
01099 
01100 /* Read the first trajectory frame. */
01101 static int read_first_frame(qmdata_t *data) {
01102   /* Try reading the first frame.
01103    * If there is only one frame then also read the
01104    * final wavefunction. */
01105   if (!get_traj_frame(data, data->atoms, data->numatoms)) {
01106     return FALSE;
01107   }
01108 
01109   return TRUE;
01110 }
01111 
01112 
01113 
01114 /******************************************************
01115  *
01116  * this function extracts the trajectory information
01117  * from the output file
01118  *
01119  * *****************************************************/
01120 static int get_traj_frame(qmdata_t *data, qm_atom_t *atoms,
01121                           int natoms) {
01122 //  orcadata *orca = (orcadata *)data->format_specific_data;
01123   qm_timestep_t *cur_ts;
01124   char buffer[BUFSIZ];
01125   char word[BUFSIZ];
01126   int units=ANGSTROM;
01127   buffer[0] = '\0';
01128   word[0]   = '\0';
01129 
01130   printf("orcaplugin) Timestep %d:\n", data->num_frames_read);
01131   printf("orcaplugin) ============\n");
01132 
01133   cur_ts = data->qm_timestep + data->num_frames_read;
01134 
01135   // debugging the trajectory reading file positions
01136   // printf("nfread: %d \n", data->num_frames_read);
01137   if (!data->filepos_array) {
01138     printf("filepos array empty!!!\n");
01139     return FALSE;
01140   }
01141 
01142   fseek(data->file, data->filepos_array[data->num_frames_read], SEEK_SET);
01143 
01144   /*
01145   * distinguish between job types
01146   * at the moment, only Single Points will work
01147   * lines 2840 - 3122 in gamessplugin.c
01148    */
01149 
01150   // reading geometries...
01151 
01152   if ((data->runtype == MOLFILE_RUNTYPE_OPTIMIZE || data->runtype == MOLFILE_RUNTYPE_GRADIENT) && data->num_frames > 1) {
01153     if (goto_keyline(data->file, "CARTESIAN COORDINATES (ANGSTROEM)", NULL)) {
01154       GET_LINE(buffer, data->file);
01155       // thisline(data->file);
01156       // UNITS ARE ANGSTROEM
01157       // bohr = 0;
01158       // sscanf()
01159     } else {
01160       printf("orcaplugin) No cartesian coordinates in ANGSTROEM found.\n");
01161     }
01162     // skip the ---- line
01163     eatline(data->file, 1);
01164     if (!get_coordinates(data->file, &data->atoms, units, &natoms)) {
01165       printf("orcaplugin) Couldn't find coordinates for timestep %d\n", data->num_frames_read);
01166     }
01167   }
01168 
01169   if (get_scfdata(data, cur_ts) == FALSE) {
01170     printf("orcaplugin) Couldn't find SCF iterations for timestep %d\n",
01171            data->num_frames_read);
01172   }
01173 
01174   /* Try reading canonical alpha/beta wavefunction */
01175   check_add_wavefunctions(data, cur_ts);
01176 
01177   /* Read point charged */
01178   if (!cur_ts->have_mulliken && get_population(data, cur_ts)) {
01179     printf("orcaplugin) Mulliken charges found\n");
01180   }
01181 
01182   /* Read the energy gradients (= -forces on atoms) */
01183   // if (get_gradient(data, cur_ts)) {
01184   //   printf("orcaplugin) Energy gradient found.\n");
01185   // }
01186 
01187   /* If this is the last frame of the trajectory and the file
01188    * wasn't truncated and the program didn't terminate
01189    * abnormally then read the final wavefunction. */
01190   if ((data->runtype == MOLFILE_RUNTYPE_OPTIMIZE ||
01191        data->runtype == MOLFILE_RUNTYPE_SADPOINT) &&
01192       (data->num_frames_read+1 == data->num_frames &&
01193        (data->status == MOLFILE_QMSTATUS_UNKNOWN ||
01194         data->status == MOLFILE_QMSTATUS_OPT_CONV ||
01195         data->status == MOLFILE_QMSTATUS_OPT_NOT_CONV))) {
01196 
01197     /* We need to jump over the end of the trajectory because
01198      * this is also the keystring for get_wavefunction() to
01199      * bail out. */
01200     if (data->status == MOLFILE_QMSTATUS_OPT_CONV ||
01201         data->status == MOLFILE_QMSTATUS_OPT_NOT_CONV) {
01202       fseek(data->file, data->end_of_traj, SEEK_SET);
01203       std::cout << "orcaplugin) Finished trajectory." << std::endl;
01204     }
01205 
01206     /* Try to read final wavefunction and orbital energies
01207      * A preexisting canonical wavefunction for this timestep
01208      * with the same characteristics (spin, exci, info) will
01209      * be overwritten by the final wavefuntion if it has more
01210      * orbitals. */
01211     check_add_wavefunctions(data, cur_ts);
01212   }
01213 
01214   data->num_frames_read++;
01215   #ifdef DEBUGGING_ORCA
01216     std::cout << "orcaplugin) Frames read: " << data->num_frames_read << std::endl;
01217   #endif
01218 
01219   return TRUE;
01220 }
01221 
01222 
01223 static bool only_numbers(const std::vector<std::string>& l) {
01224   //for (auto s : l) {
01225   for (int i = 0; i < l.size(); i++){
01226     if (l[i].find_first_not_of("0123456789-.") != std::string::npos) {
01227       return false;
01228     }
01229   }
01230   return true;
01231 }
01232 
01233 
01234 static int get_scfdata(qmdata_t *data, qm_timestep_t *ts) {
01235   char buffer[BUFSIZ];
01236   long filepos;
01237   filepos = ftell(data->file);
01238   if (!goto_keyline(data->file, "SCF ITERATIONS", NULL)) {
01239     fseek(data->file, filepos, SEEK_SET);
01240     ts->num_scfiter = 0;
01241     return FALSE;
01242   }
01243   eatline(data->file, 2);
01244   GET_LINE(buffer, data->file);
01245   std::string scfIterLine(buffer);
01246   int ncols = (split(reduce(scfIterLine), ' ')).size();
01247   int currentCols, numiter = 0;
01248   filepos = ftell(data->file); // save position for later
01249   while (scfIterLine.find("SUCCESS") == std::string::npos && scfIterLine.find("ERROR") == std::string::npos) {
01250     GET_LINE(buffer, data->file);
01251     scfIterLine = buffer;
01252     std::vector<std::string> currentCol = split(reduce(scfIterLine), ' ');
01253     currentCols = currentCol.size();
01254     if (currentCols == ncols && only_numbers(currentCol)) {
01255       numiter++;
01256     }
01257   }
01258   ts->num_scfiter = numiter;
01259   std::cout << "orcaplugin) Number of SCF iterations: " << numiter << std::endl;
01260   // NEW: Added SCF energies
01261   ts->scfenergies = (double *)calloc(numiter,sizeof(double));
01262   numiter = 0; // reset to check consistency
01263   int numread, dum;
01264   fseek(data->file, filepos, SEEK_SET); // go back to read scfenergies
01265   GET_LINE(buffer, data->file);
01266   scfIterLine = buffer;
01267   while (scfIterLine.find("SUCCESS") == std::string::npos && scfIterLine.find("ERROR") == std::string::npos) {
01268     GET_LINE(buffer, data->file);
01269     scfIterLine = buffer;
01270     std::vector<std::string> currentCol = split(reduce(scfIterLine), ' ');
01271     currentCols = currentCol.size();
01272     if (currentCols == ncols && only_numbers(currentCol)) {
01273       numread = sscanf(buffer,"%i %lf", &dum, &ts->scfenergies[numiter]);
01274       numiter++;
01275     }
01276   }  
01277   
01278   return TRUE;
01279 }
01280 
01281 static int get_population(qmdata_t *data, qm_timestep_t *ts) {
01282   int i;
01283   char buffer[BUFSIZ];
01284   long filepos;
01285   ts->have_mulliken = FALSE;
01286   ts->have_lowdin   = FALSE;
01287   ts->have_esp      = FALSE;
01288   filepos = ftell(data->file);
01289 
01290   if (pass_keyline(data->file, "MULLIKEN ATOMIC CHARGES", NULL) != FOUND) {
01291     fseek(data->file, filepos, SEEK_SET);
01292     return FALSE;
01293   }
01294 
01295   /* Read Mulliken charges if present */
01296   ts->mulliken_charges = (double *)calloc(data->numatoms, sizeof(double));
01297 
01298   if (!ts->mulliken_charges) {
01299     PRINTERR;
01300     return FALSE;
01301   }
01302 
01303   // ts->lowdin_charges = (double *)calloc(data->numatoms, sizeof(double));
01304   //
01305   // if (!ts->lowdin_charges) {
01306   //   free(ts->mulliken_charges);
01307   //   ts->mulliken_charges = NULL;
01308   //   PRINTERR;
01309   //   return FALSE;
01310   // }
01311   eatline(data->file, 1);
01312 
01313   for (i=0; i<data->numatoms; i++) {
01314     int n;
01315     float mullcharge;
01316     GET_LINE(buffer, data->file);
01317     std::string currentLine(buffer);
01318     std::vector<std::string> elements = split(reduce(currentLine), ' ');
01319     n = elements.size();
01320     if (n!=4) {
01321       free(ts->mulliken_charges);
01322       // free(ts->lowdin_charges);
01323       ts->mulliken_charges = NULL;
01324       // ts->lowdin_charges   = NULL;
01325       return FALSE;
01326     }
01327     mullcharge = (float)atof((*(elements.end()-1)).c_str());
01328     ts->mulliken_charges[i] = mullcharge;
01329     // ts->lowdin_charges[i]   = lowcharge;
01330   }
01331 
01332   if (i!=data->numatoms) {
01333     free(ts->mulliken_charges);
01334     free(ts->lowdin_charges);
01335     ts->mulliken_charges = NULL;
01336     ts->lowdin_charges   = NULL;
01337     return FALSE;
01338   }
01339 
01340   ts->have_mulliken = TRUE;
01341   // ts->have_lowdin   = TRUE;
01342   return TRUE;
01343 }
01344 
01345 /*********************************************************
01346  *
01347  * Reads a set of wavefunctions for the current timestep.
01348  * These are typically the alpha and beta spin wavefunctions
01349  * or the MCSCF natural and optimized orbitals or the GVB
01350  * canonical orbitals and geminal pairs.
01351  *
01352  **********************************************************/
01353  // 3461
01354 static int check_add_wavefunctions(qmdata_t *data, qm_timestep_t *ts) {
01355   qm_wavefunction_t *wavef;
01356   int i, n=1;
01357 
01358   if (data->scftype==MOLFILE_SCFTYPE_UHF ||
01359       data->scftype==MOLFILE_SCFTYPE_GVB ||
01360       data->scftype==MOLFILE_SCFTYPE_MCSCF) {
01361     /* Try to read second wavefunction, e.g. spin beta */
01362     n = 2;
01363   }
01364 
01365   for (i=0; i<n; i++) {
01366     /* Allocate memory for new wavefunction */
01367     wavef = add_wavefunction(ts);
01368 
01369     /* Try to read wavefunction and orbital energies */
01370     if (get_wavefunction(data, ts, wavef) == FALSE) {
01371       /* Free the last wavefunction again. */
01372       del_wavefunction(ts);
01373 #ifdef DEBUGGING_ORCA
01374       printf("orcaplugin) No canonical wavefunction present for timestep %d\n", data->num_frames_read);
01375 #endif
01376       break;
01377 
01378     } else {
01379       char action[32];
01380       char spinstr[32];
01381       strcpy(spinstr, "");
01382       if (data->scftype==MOLFILE_SCFTYPE_UHF) {
01383         if (wavef->spin==SPIN_BETA) {
01384           strcat(spinstr, "spin  beta, ");
01385         } else {
01386           strcat(spinstr, "spin alpha, ");
01387         }
01388       }
01389 
01390       /* The last SCF energy is the energy of this electronic state */
01391       if (ts->scfenergies) {
01392         wavef->energy = ts->scfenergies[ts->num_scfiter-1];
01393       } else {
01394         wavef->energy = 0.f;
01395       }
01396 
01397       /* Multiplicity */
01398       wavef->mult = data->multiplicity;
01399 
01400 
01401       /* String telling wether wavefunction was added, updated
01402        * or ignored. */
01403       strcpy(action, "added");
01404 
01405       /* If there exists a canonical wavefunction of the same spin
01406        * we'll replace it */
01407       if (ts->numwave>1 && wavef->type==MOLFILE_WAVE_CANON) {
01408         int i, found =-1;
01409         for (i=0; i<ts->numwave-1; i++) {
01410           if (ts->wave[i].type==wavef->type &&
01411               ts->wave[i].spin==wavef->spin &&
01412               ts->wave[i].exci==wavef->exci &&
01413               !strncmp(ts->wave[i].info, wavef->info, MOLFILE_BUFSIZ)) {
01414             found = i;
01415             break;
01416           }
01417         }
01418         if (found>=0) {
01419           /* If the new wavefunction has more orbitals we
01420            * replace the old one for this step. */
01421           if (wavef->num_orbitals >
01422               ts->wave[found].num_orbitals) {
01423             /* Replace existing wavefunction for this step */
01424             replace_wavefunction(ts, found);
01425             sprintf(action, "%d updated", found);
01426           } else {
01427             /* Delete last wavefunction again */
01428             del_wavefunction(ts);
01429             sprintf(action, "matching %d ignored", found);
01430           }
01431           wavef = &ts->wave[ts->numwave-1];
01432         }
01433       }
01434 
01435       printf("orcaplugin) Wavefunction %s (%s):\n", action, wavef->info);
01436       printf("orcaplugin)   %d orbitals, %sexcitation %d, multiplicity %d\n",
01437              wavef->num_orbitals, spinstr, wavef->exci, wavef->mult);
01438     }
01439   }
01440 
01441   return i;
01442 }
01443 
01444 //std::vector<std::vector<int> > dAngMom{{1,0,1},{0,1,1},{1,1,0},{2,0,0},{0,2,0},{0,0,2}};
01445 //std::vector<std::vector<int> > fAngMom{{1,1,1},{2,1,0},{2,0,1},{1,2,0},{0,2,1},{1,0,2},
01446 //  {0,1,2},{3,0,0},{0,3,0},{0,0,3}};
01447 int dAngMom[18] = {1,0,1, 0,1,1, 1,1,0, 2,0,0, 0,2,0, 0,0,2};
01448 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};
01449 
01450 typedef enum {
01451   dshell, fshell
01452 } Shelltype;
01453 
01454 static CoeffRowBlock convertPure(CoeffRowBlock pureBlock) {
01455   Shelltype stype;
01456   CoeffRowBlock resultBlock;
01457   switch (pureBlock.size()) {
01458     case 5:
01459       stype = dshell;
01460       break;
01461     case 7:
01462       stype = fshell;
01463       break;
01464     default:
01465       return resultBlock;
01466   }
01467 
01468   Matrix *m = new Matrix(pureBlock);
01469   Matrix *multiplication;
01470   if (stype == dshell) {
01471     multiplication = Matrix::multiply(convD, m);
01472   } else if (stype == fshell) {
01473     multiplication = Matrix::multiply(convF, m);
01474   } else {
01475     return resultBlock;
01476   }
01477   // m->printMatrix();
01478   // multiplication->printMatrix();
01479   resultBlock = multiplication->toVector();
01480   delete multiplication;
01481   delete m;
01482 
01483   // OLD CODE!
01484   // first get the unchanged d-orbitals xz, yz, xy
01485   // std::vector<std::string> unchangedOrbitals{"xz", "yz", "xy"};
01486   // std::vector<int> unchangedIndices{1, 2, 4};
01487   // for (auto idx : unchangedIndices) {
01488   //   resultBlock.push_back(pureBlock[idx]);
01489   // }
01490   //
01491   // int alpha = 3;
01492   // int beta = 0;
01493   // // now convert the others
01494   //
01495   // int rows = pureBlock[0].size(); // should be 5 in this case
01496   //
01497   // std::vector<float> x2;
01498   // for (size_t i = 0; i < pureBlock[0].size(); i++) {
01499   //   x2.push_back(pureBlock[alpha][i]-pureBlock[beta][i]);
01500   // }
01501   //
01502   // std::vector<float> y2;
01503   // for (size_t i = 0; i < pureBlock[0].size(); i++) {
01504   //   y2.push_back(-pureBlock[beta][i]-pureBlock[alpha][i]);
01505   // }
01506   //
01507   // std::vector<float> z2;
01508   // for (size_t i = 0; i < pureBlock[0].size(); i++) {
01509   //   z2.push_back(2*pureBlock[beta][i]);
01510   // }
01511   //
01512   // resultBlock.push_back(x2);
01513   // resultBlock.push_back(y2);
01514   // resultBlock.push_back(z2);
01515 
01516   return resultBlock;
01517 }
01518 
01519 static int get_wavefunction(qmdata_t *data, qm_timestep_t *ts, qm_wavefunction_t *wf) {
01520   std::vector<float> orbitalEnergies;
01521   std::vector<int> orbitalOccupancies;
01522   std::vector<float> wavefunctionCoeffs;
01523   int num_orbitals = 0;
01524 
01525   char buffer[BUFSIZ];
01526   char word[6][BUFSIZ];
01527   long filepos;
01528   char *line;
01529 
01530   buffer[0] = '\0';
01531   int i = 0;
01532   for (i=0; i<6; i++) word[i][0] = '\0';
01533 
01534   if (wf == NULL) {
01535     PRINTERR;
01536     return FALSE;
01537   }
01538 
01539   wf->has_occup = FALSE;
01540   wf->has_orben = FALSE;
01541   wf->type = MOLFILE_WAVE_UNKNOWN;
01542   wf->spin = SPIN_ALPHA;
01543   wf->exci = 0;
01544   strncpy(wf->info, "unknown", MOLFILE_BUFSIZ);
01545 
01546   filepos = ftell(data->file);
01547 
01548   do {
01549     GET_LINE(buffer, data->file);
01550     line = trimleft(trimright(buffer));
01551     if(!strcmp(line, "MOLECULAR ORBITALS")) {
01552       wf->type = MOLFILE_WAVE_CANON;
01553       strncpy(wf->info, "canonical", MOLFILE_BUFSIZ);
01554     }
01555   } while(wf->type == MOLFILE_WAVE_UNKNOWN && strcmp(line, "FINAL SINGLE POINT ENERGY"));
01556 
01557   if(wf->type == MOLFILE_WAVE_UNKNOWN) {
01558     #ifdef DEBUGGING_ORCA
01559         printf("orcaplugin) get_wavefunction(): No wavefunction found!\n");
01560     #endif
01561         fseek(data->file, filepos, SEEK_SET);
01562         return FALSE;
01563   } else {
01564     #ifdef DEBUGGING_ORCA
01565       printf("orcaplugin) Found wavefunction of type %d.\n", wf->type);
01566     #endif
01567   }
01568 
01569   eatline(data->file, 1);
01570 
01571   // number of read values from line;
01572   int numReadOrbitalIndices = 0;
01573   int numReadEnergies = 0;
01574   int numReadOccupancies = 0;
01575   int numReadCoefficients = 0;
01576   float numberOfElectrons = 0;
01577   float occupiedOrbitals = 0;
01578   std::vector<int> numberContractedBf;
01579   std::vector<int> wfAngMoment;
01580   std::vector<std::string> orbitalNames;
01581   MoCoeff allCoefficients;
01582 
01583   // orbital indices
01584   int n[6];
01585   int wavefunctionRead = 0;
01586   int firstRead = 1;
01587   int haveAngMom = 0;
01588   while(!wavefunctionRead) {
01589     float coeff[6], energies[6];
01590     float occ[6];
01591     char dumpName[BUFSIZ];
01592     char dumpBasisFunc[BUFSIZ];
01593     filepos = ftell(data->file);
01594 
01595     // reads the orbital indices
01596     if (firstRead == 1) {
01597       GET_LINE(buffer, data->file);
01598       firstRead++;
01599     }
01600 
01601     numReadOrbitalIndices = sscanf(buffer, "%d %d %d %d %d %d", &n[0], &n[1], &n[2], &n[3], &n[4], &n[5]);
01602     if (!numReadOrbitalIndices || numReadOrbitalIndices == -1) {
01603       /* If there are no orbital indexes then this must be the
01604        * end of the wavefunction coefficient table. */
01605       fseek(data->file, filepos, SEEK_SET);
01606       wavefunctionRead = 1;
01607       break;
01608     }
01609 
01610     // reads the orbital energies
01611     GET_LINE(buffer, data->file);
01612     numReadEnergies = sscanf(buffer, "%f %f %f %f %f %f", &energies[0], &energies[1], &energies[2],
01613      &energies[3], &energies[4], &energies[5]);
01614     if (numReadEnergies != numReadOrbitalIndices) {
01615       printf("orcaplugin) Molecular Orbital section corrupted! energies.\n");
01616       break;
01617     }
01618 
01619     // store the energies in vector
01620     if (numReadEnergies != -1) {
01621       for (size_t c = 0; c < numReadEnergies; c++) {
01622         orbitalEnergies.push_back(energies[c]);
01623         // std::cout << "Energy: " <<energies[c]<< std::endl;
01624       }
01625     }
01626     // reads the orbital occupancies
01627     GET_LINE(buffer, data->file);
01628     numReadOccupancies = sscanf(buffer, "%f %f %f %f %f %f", &occ[0], &occ[1], &occ[2],
01629       &occ[3], &occ[4], &occ[5]);
01630     if (numReadOccupancies != numReadOrbitalIndices) {
01631       printf("orcaplugin) Molecular Orbital section corrupted! %d\n", numReadOccupancies);
01632       break;
01633     }
01634 
01635     // stores the occupancies in vector
01636     if(numReadOccupancies) {
01637       for (size_t c = 0; c < numReadOccupancies; c++) {
01638         orbitalOccupancies.push_back((int)occ[c]);
01639         // std::cout << "Occupancy: " << occ[c] << std::endl;
01640         numberOfElectrons += occ[c];
01641         if (occ[c]) {
01642           occupiedOrbitals++;
01643         }
01644       }
01645       num_orbitals += numReadOccupancies;
01646     }
01647 
01648     // skip --- line
01649     eatline(data->file, 1);
01650 
01651     std::vector<std::vector<float> > moCoefficients;
01652     // we expect as many coefficients as numReadOccupancies, numReadEnergies, numReadOrbitalIndices!
01653     // read them as long as we find new coefficients
01654     int readingBlock = 1;
01655     int coefficientNumber = 0;
01656     while(readingBlock) {
01657       GET_LINE(buffer, data->file);
01658       numReadCoefficients = sscanf(buffer, "%s %s %f %f %f %f %f %f", &dumpName, &dumpBasisFunc,
01659       &coeff[0], &coeff[1], &coeff[2],&coeff[3], &coeff[4], &coeff[5]);
01660       // the coefficient number is the number of read elements minus 2 bc. of the atom and bf name
01661       coefficientNumber = (numReadCoefficients - 2);
01662       if (coefficientNumber == numReadOrbitalIndices) {
01663         // std::cout << "found coeffs: " << dumpName << "," << dumpBasisFunc << std::endl;
01664         if (firstRead == 2 && !haveAngMom) {
01665           std::string bfn = dumpBasisFunc;
01666           size_t found = bfn.find_first_not_of("-0123456789 ");
01667           if (found!=std::string::npos) {
01668             std::string orbital =  bfn.substr(found);
01669             orbitalNames.push_back(orbital);
01670             // std::cout << orbital << std::endl;
01671           } else {
01672             printf("orcaplugin) Could not determine orbital description.\n");
01673             return FALSE;
01674           }
01675         }
01676         // reading coefficients
01677         std::vector<float> currentMoCoeffs;
01678         for (size_t cidx = 0; cidx < coefficientNumber; cidx++) {
01679           currentMoCoeffs.push_back(coeff[cidx]);
01680           // std::cout << coeff[cidx] << std::endl;
01681         }
01682         moCoefficients.push_back(currentMoCoeffs);
01683       } else {
01684         // block seems to be finished
01685         readingBlock = 0;
01686         haveAngMom = 1;
01687       }
01688     }
01689     allCoefficients.push_back(moCoefficients);
01690   }
01691 
01692   if ( std::adjacent_find( numberContractedBf.begin(), numberContractedBf.end(), std::not_equal_to<int>() ) != numberContractedBf.end() ) {
01693     printf("orcaplugin) Molecular orbital section corrupted. Did not read consistent number of contracted basis functions!\n");
01694     return FALSE;
01695   }
01696 
01697   // now loop over MO blocks and convert coefficients if needed
01698   int readingPureFunction = 0;
01699   std::vector<std::vector<float> > pureFunction;
01700   std::vector<std::string> pureFunctionName;
01701   int expectedNumberOfPureFunctions = 0;
01702   //std::vector<std::string> orbList{"s","px","py","pz","dz2","dxz","dyz","dx2y2","dxy",
01703   //"f0","f+1","f-1","f+2","f-2","f+3","f-3"};
01704   const char* orbitalList[] = {"s","px","py","pz","dz2","dxz","dyz","dx2y2","dxy",
01705   "f0","f+1","f-1","f+2","f-2","f+3","f-3"};
01706   std::vector<std::string> orbList(orbitalList, orbitalList+15);
01707   int orbRowIndex = 0;
01708   int blockIdx = 0;
01709 
01710   // for (auto name : orbitalNames) {
01711   //   std::cout << name << std::endl;
01712   // }
01713 
01714   MoCoeff newAllCoefficients;
01715   //MoCoeff is a vector<vector<vecto<float>r>>
01716   //for (auto moBlock : allCoefficients) {
01717   for (size_t i = 0; i < allCoefficients.size(); i++) {
01718     CoeffRowBlock newRows;
01719     int blockNumberOfContracted = 0;
01720     //for (auto moRow : moBlock) {
01721     for (size_t j = 0; j < allCoefficients[i].size(); j++) {
01722       int angX = 0, angY = 0, angZ = 0;
01723       std::string orbital = orbitalNames[orbRowIndex];
01724       //std::begin and std::end are not supported on c++98
01725       //std::vector<std::string>::iterator orbIndex = find(std::begin(orbList), std::end(orbList), orbital);
01726       std::vector<std::string>::iterator orbIndex = find(orbList.begin(), orbList.end(), orbital);
01727       //if (std::end(orbList) != orbIndex) {
01728       if (orbList.end() != orbIndex) {
01729         int listIndex = orbIndex-orbList.begin();
01730         switch (listIndex) {
01731           case 0:
01732             break;
01733           case 1:
01734             angX = 1;
01735             break;
01736           case 2:
01737             angY = 1;
01738             break;
01739           case 3:
01740             angZ = 1;
01741             break;
01742           default:
01743             break;
01744         }
01745         // d-shell found
01746         // std::cout << "orbital index: " << orbIndex-orbList.begin() << " " << orbital << std::endl;
01747         if (listIndex > 3) {
01748           if (!readingPureFunction) {
01749             pureFunction.clear();
01750             pureFunctionName.clear();
01751             if (orbital.compare(0, 1, "d") == 0) {
01752               expectedNumberOfPureFunctions = 5;
01753               // std::cout << "Trying to read d-functions." << std::endl;
01754             } else if (orbital.compare(0, 1, "f") == 0) {
01755               expectedNumberOfPureFunctions = 7;
01756               // std::cout << "Trying to read f-functions." << std::endl;
01757             }
01758             readingPureFunction = 1;
01759            // pureFunction.push_back(moRow);
01760             pureFunction.push_back(allCoefficients[i][j]);
01761             pureFunctionName.push_back(orbital);
01762           } else {
01763             //pureFunction.push_back(moRow);
01764             pureFunction.push_back(allCoefficients[i][j]);
01765             pureFunctionName.push_back(orbital);
01766             if (pureFunction.size() == expectedNumberOfPureFunctions) {
01767               // std::cout << "found complete pure function set." << std::endl;
01768               //CoeffRowBlock is a vector<vector<float>>
01769               CoeffRowBlock newBlock = convertPure(pureFunction);
01770               blockNumberOfContracted+= newBlock.size();
01771               //for (auto r : newBlock) {
01772               for (size_t k = 0; k < newBlock.size(); k++) {
01773                 newRows.push_back(newBlock[k]);
01774               }
01775               if (!blockIdx) {
01776                 for (size_t i = 0; i < newBlock.size(); i++) {
01777                   if (expectedNumberOfPureFunctions == 5) {
01778                     //for (auto angMom : dAngMom[i]) {
01779                     //wfAngMoment.push_back(dAngMom[i][0]);
01780                     //wfAngMoment.push_back(dAngMom[i][1]);
01781                     //wfAngMoment.push_back(dAngMom[i][2]);
01782                     wfAngMoment.push_back(dAngMom[i*3 + 0]);
01783                     wfAngMoment.push_back(dAngMom[i*3 + 1]);
01784                     wfAngMoment.push_back(dAngMom[i*3 + 2]);
01785                     //}
01786                   } else if (expectedNumberOfPureFunctions == 7) {
01787                     //for (auto angMom : fAngMom[i]) {
01788                     //wfAngMoment.push_back(fAngMom[i][0]);
01789                     //wfAngMoment.push_back(fAngMom[i][1]);
01790                     //wfAngMoment.push_back(fAngMom[i][2]);
01791                     wfAngMoment.push_back(fAngMom[i*3 + 0]);
01792                     wfAngMoment.push_back(fAngMom[i*3 + 1]);
01793                     wfAngMoment.push_back(fAngMom[i*3 + 2]);
01794                     //}
01795                   }
01796                 }
01797               }
01798               readingPureFunction = 0;
01799             }
01800           }
01801         } else {
01802           newRows.push_back(allCoefficients[i][j]);
01803           blockNumberOfContracted++;
01804           if (!blockIdx) {
01805             wfAngMoment.push_back(angX);
01806             wfAngMoment.push_back(angY);
01807             wfAngMoment.push_back(angZ);
01808           }
01809         }
01810       } else {
01811         std::cout << "orcaplugin) ERROR. Only s/p/d/f-shells are supported." << std::endl;
01812         return FALSE;
01813       }
01814 
01815       orbRowIndex++;
01816     }
01817     numberContractedBf.push_back(blockNumberOfContracted);
01818     newAllCoefficients.push_back(newRows);
01819     orbRowIndex = 0;
01820     blockIdx++;
01821   }
01822 
01823 
01824   // assign the number of contracted functions to wavefunction size
01825   data->wavef_size = numberContractedBf[0];
01826   wf->num_orbitals  = num_orbitals;
01827   wf->orb_energies = (float *) calloc(num_orbitals, sizeof(float));
01828   wf->orb_occupancies = (float *) calloc(num_orbitals, sizeof(float));
01829   wf->wave_coeffs = (float *) calloc(num_orbitals * data->wavef_size, sizeof(float));
01830   wf->has_occup = TRUE;
01831   wf->has_orben = TRUE;
01832 
01833   int cnt = 0;
01834   
01835   //orbitalEnergies is a std::vector<float>
01836   //Ranges do not work...
01837   //for (auto en : orbitalEnergies) {
01838   for (int i = 0; i < orbitalEnergies.size(); i++) {
01839     //wf->orb_energies[cnt] = en;
01840     //cnt++;
01841     wf->orb_energies[i] = orbitalEnergies[i];
01842   }
01843   cnt = 0;
01844   //for (auto occ : orbitalOccupancies) {
01845   for (int i = 0; i < orbitalOccupancies.size();i++) {
01846     //wf->orb_occupancies[cnt] = occ;
01847     //cnt++;
01848     wf->orb_occupancies[i] = orbitalOccupancies[i];
01849   }
01850 
01851   int rowIndex = 0, columnIndex = 0;
01852   blockIdx = 0;
01853   int moBlockSize = 0;
01854   int moBlockIdx = 0;
01855   //for (auto moBlock : newAllCoefficients) {
01856   for (int i = 0; i < newAllCoefficients.size();i++) {
01857     //for (auto moRow : moBlock) {
01858     for (int j = 0; j < newAllCoefficients[i].size(); j++) {
01859       // std::cout << rowIndex << std::endl;
01860       //for (auto moCo : moRow) {
01861       for (int k = 0; k < newAllCoefficients[i][j].size(); k++) {
01862         if ((columnIndex * data->wavef_size + rowIndex) > num_orbitals * data->wavef_size) {
01863           std::cout << "something went wrong:" << columnIndex << std::endl;
01864           std::cout << "something went wrong:" << (columnIndex * data->wavef_size + rowIndex) << " vs. " << num_orbitals * data->wavef_size << std::endl;
01865           return FALSE;
01866         }
01867         // std::cout << orbitalNames[rowIndex] << std::endl;
01868         //wf->wave_coeffs[columnIndex * data->wavef_size + rowIndex] = moCo;
01869         wf->wave_coeffs[columnIndex * data->wavef_size + rowIndex] = newAllCoefficients[i][j][k];
01870         columnIndex++;
01871       }
01872       columnIndex = moBlockSize;
01873       rowIndex++;
01874     }
01875     rowIndex = 0;
01876     // 0-based!!!
01877     //moBlockSize += moBlock[moBlockIdx].size();
01878     moBlockSize += newAllCoefficients[i][moBlockIdx].size();
01879     columnIndex = moBlockSize;
01880     moBlockIdx++;
01881     // std::cout << "bs: " << moBlockSize << std::endl;
01882   }
01883 
01884   // LOGGING for MO coefficients
01885   /*
01886   float coeff2 = 0;
01887   for (size_t t = 0; t < (num_orbitals * data->wavef_size); t++) {
01888     if (t % data->wavef_size == 0) {
01889       std::cout << "---------- " << t/num_orbitals << " c2: " << coeff2 << std::endl;
01890       coeff2 = 0;
01891     }
01892     coeff2 += wf->wave_coeffs[t]*wf->wave_coeffs[t];
01893     std::cout << wf->wave_coeffs[t] << std::endl;
01894   }
01895   */
01896 
01897   if (data->num_frames_read < 1) {
01898     data->angular_momentum = (int*)calloc(wfAngMoment.size(), sizeof(int));
01899   }
01900 
01901   // std::cout << "wfang: " << wfAngMoment.size() <<  " " << 3*data->wavef_size <<std::endl;
01902   for (size_t ang = 0; ang < wfAngMoment.size(); ang++) {
01903     data->angular_momentum[ang] = wfAngMoment[ang];
01904   }
01905 
01906   // TODO: This is just a workaround and might give wrong
01907   // results when reading unrestricted jobs
01908   data->num_occupied_A = occupiedOrbitals;
01909   data->num_occupied_B = occupiedOrbitals;
01910   // data->num_electrons = numberOfElectrons;
01911 
01912   // if (data->num_electrons != numberOfElectrons) {
01913   //
01914   // }
01915 
01916   std::cout << "----------------------------------------" << std::endl;
01917   std::cout << "Total number of orbitals: " << num_orbitals << std::endl;
01918   std::cout << "Number of electrons in read wf: " << numberOfElectrons<< std::endl;
01919   std::cout << "Number of occupied orbitals: " << occupiedOrbitals << std::endl;
01920   std::cout << "Number of contracted bf: " << numberContractedBf[0] << std::endl;
01921   std::cout << "----------------------------------------" << std::endl;
01922 
01923   return TRUE;
01924 }
01925 
01926 /******************************************************
01927  *
01928  * Populate the flat arrays containing the basis
01929  * set data.
01930  *
01931  ******************************************************/
01932 static int fill_basis_arrays(qmdata_t *data) {
01933   //orcadata *orca = (orcadata *)data->format_specific_data;
01934   int i, j, k;
01935   int shellcount = 0;
01936   int primcount = 0;
01937 
01938   float *basis;
01939   int *num_shells_per_atom;
01940   int *num_prim_per_shell;
01941   int *shell_types;
01942   int *atomicnum_per_basisatom;
01943 
01944   /* Count the total number of primitives which
01945    * determines the size of the basis array. */
01946   for (i=0; i<data->num_basis_atoms; i++) {
01947     for (j=0; j<data->basis_set[i].numshells; j++) {
01948       primcount += data->basis_set[i].shell[j].numprims;
01949     }
01950   }
01951 
01952   /* reserve space for pointer to array containing basis
01953    * info, i.e. contraction coeficients and expansion
01954    * coefficients; need 2 entries per basis function, i.e.
01955    * exponent and contraction coefficient; also,
01956    * allocate space for the array holding the orbital symmetry
01957    * information per primitive Gaussian.
01958    * Finally, initialize the arrays holding the number of
01959    * shells per atom and the number of primitives per shell*/
01960   basis = (float *)calloc(2*primcount,sizeof(float));
01961 
01962   /* make sure memory was allocated properly */
01963   if (basis == NULL) {
01964     PRINTERR;
01965     return FALSE;
01966   }
01967 
01968   shell_types = (int *)calloc(data->num_shells, sizeof(int));
01969 
01970   /* make sure memory was allocated properly */
01971   if (shell_types == NULL) {
01972     PRINTERR;
01973     return FALSE;
01974   }
01975 
01976   num_shells_per_atom = (int *)calloc(data->num_basis_atoms, sizeof(int));
01977 
01978   /* make sure memory was allocated properly */
01979   if (num_shells_per_atom == NULL) {
01980     PRINTERR;
01981     return FALSE;
01982   }
01983 
01984   num_prim_per_shell = (int *)calloc(data->num_shells, sizeof(int));
01985 
01986   /* make sure memory was allocated properly */
01987   if (num_prim_per_shell == NULL) {
01988     PRINTERR;
01989     return FALSE;
01990   }
01991 
01992   atomicnum_per_basisatom = (int *)calloc(data->num_basis_atoms, sizeof(int));
01993 
01994   /* make sure memory was allocated properly */
01995   if (atomicnum_per_basisatom == NULL) {
01996     PRINTERR;
01997     return FALSE;
01998   }
01999 
02000 
02001   /* store pointers in struct qmdata_t */
02002   data->basis = basis;
02003   data->shell_types = shell_types;
02004   data->num_shells_per_atom = num_shells_per_atom;
02005   data->num_prim_per_shell = num_prim_per_shell;
02006   data->atomicnum_per_basisatom = atomicnum_per_basisatom;
02007 
02008   /* Go through all basis set atoms and try to assign the
02009    * atomic numbers. The basis set atoms are specified by
02010    * name strings (the same as in the coordinate section,
02011    * except for FMO calcs.) and we try to match the names
02012    * from the two lists. The basis set atom list is symmetry
02013    * unique while the coordinate atom list is complete.*/
02014   primcount = 0;
02015   for (i=0; i<data->num_basis_atoms; i++) {
02016     int found = 0;
02017 
02018     /* For this basis atom find a matching atom from the
02019      * coordinate atom list. */
02020     for (j=0; j<data->numatoms; j++) {
02021       char basisname[BUFSIZ];
02022       strcpy(basisname, data->basis_set[i].name);
02023 
02024       /* for FMO calculations we have to strip the "-n" tail
02025        * of the basis atom name. */
02026       // if (gms->have_fmo) {
02027       //   *strchr(basisname, '-') = '\0';
02028       // }
02029 
02030       if (!strcmp(data->atoms[j].type, basisname)) {
02031         found = 1;
02032         break;
02033       }
02034     }
02035     if (!found) {
02036       printf("orcaplugin) WARNING: Couldn't find atomic number for basis set atom %s\n",
02037              data->basis_set[i].name);
02038       data->basis_set[i].atomicnum = 0;
02039       atomicnum_per_basisatom[i] = 0;
02040     } else {
02041       /* assign atomic number */
02042       data->basis_set[i].atomicnum = data->atoms[j].atomicnum;
02043       atomicnum_per_basisatom[i]   = data->atoms[j].atomicnum;
02044     }
02045     num_shells_per_atom[i] = data->basis_set[i].numshells;
02046 
02047     for (j=0; j<data->basis_set[i].numshells; j++) {
02048       shell_types[shellcount] = data->basis_set[i].shell[j].type;
02049       num_prim_per_shell[shellcount] = data->basis_set[i].shell[j].numprims;
02050 
02051       for (k=0; k<data->basis_set[i].shell[j].numprims; k++) {
02052         basis[2*primcount  ] = data->basis_set[i].shell[j].prim[k].exponent;
02053         basis[2*primcount+1] = data->basis_set[i].shell[j].prim[k].contraction_coeff;
02054         primcount++;
02055       }
02056       shellcount++;
02057     }
02058   }
02059   printf("orcaplugin) Filled basis arrays.\n");
02060 
02061   return TRUE;
02062 }
02063 
02064 
02065 /**************************************************
02066  *
02067  * Convert shell type from char to int.
02068  *
02069  ************************************************ */
02070 static int shelltype_int(char type) {
02071   int shelltype;
02072 
02073   switch (type) {
02074     case 'L':
02075       /* SP_P shells are assigned in get_basis() */
02076       shelltype = SP_S_SHELL;
02077       break;
02078     case 'S':
02079       shelltype = S_SHELL;
02080       break;
02081     case 'P':
02082       shelltype = P_SHELL;
02083       break;
02084     case 'D':
02085       shelltype = D_SHELL;
02086       break;
02087     case 'F':
02088       shelltype = F_SHELL;
02089       break;
02090     case 'G':
02091       shelltype = G_SHELL;
02092       break;
02093     default:
02094       shelltype = UNK_SHELL;
02095       break;
02096   }
02097 
02098   return shelltype;
02099 }
02100 
02101 /* Analyze the trajectory.
02102  * Read the parameters controlling geometry search and
02103  * find the end of the trajectory, couinting the frames
02104  * on the way. Store the filepointer for the beginning of
02105  * each frame in *filepos_array. */
02106 static int analyze_traj(qmdata_t *data, orcadata *orca) {
02107   char buffer[BUFSIZ];
02108   char *line;
02109   long filepos;
02110   filepos = ftell(data->file);
02111 
02112   data->filepos_array = (long* )calloc(1, sizeof(long ));
02113 
02114   /* currently, only one frame is supported!
02115    * lines 3130-3348 in gamessplugin.c
02116    */
02117   if (data->runtype == MOLFILE_RUNTYPE_ENERGY) {
02118     /* We have just one frame */
02119     data->num_frames = 1;
02120     pass_keyline(data->file, "Single Point Calculation", NULL);
02121     data->filepos_array[0] = ftell(data->file);
02122 
02123     /* Check wether SCF has converged */
02124     // if (pass_keyline(data->file,
02125     //                  "SCF IS UNCONVERGED, TOO MANY ITERATIONS",
02126     //                  "ENERGY COMPONENTS")==FOUND) {
02127     //   printf("orcaplugin) SCF IS UNCONVERGED, TOO MANY ITERATIONS\n");
02128     //   data->status = MOLFILE_QMSTATUS_SCF_NOT_CONV;
02129     // } else {
02130     //   data->status = MOLFILE_QMSTATUS_OPT_CONV;
02131     //   fseek(data->file, data->filepos_array[0], SEEK_SET);
02132     // }
02133 
02134     pass_keyline(data->file, "FINAL SINGLE POINT ENERGY", NULL);
02135     data->end_of_traj = ftell(data->file);
02136 
02137     /* Allocate memory for the frame */
02138     data->qm_timestep = (qm_timestep_t *)calloc(1, sizeof(qm_timestep_t));
02139     memset(data->qm_timestep, 0, sizeof(qm_timestep_t));
02140 
02141     return TRUE;
02142   } else if (data->runtype == MOLFILE_RUNTYPE_GRADIENT) {
02143     int appendedCalculations = 0;
02144     rewind(data->file);
02145     pass_keyline(data->file, "Energy+Gradient Calculation", NULL);
02146     data->filepos_array[0] = ftell(data->file);
02147     data->num_frames = 1;
02148 
02149     while(TRUE) {
02150       if (!fgets(buffer, sizeof(buffer), data->file)) break;
02151       line = trimleft(buffer);
02152 
02153       std::string l(line);
02154       if (l.find("Energy+Gradient Calculation") != std::string::npos && data->runtype==MOLFILE_RUNTYPE_GRADIENT) {
02155         appendedCalculations++;
02156         // std::cout << l << std::endl;
02157         if (data->num_frames > 0) {
02158           data->filepos_array = (long*)realloc(data->filepos_array, (data->num_frames+1)*sizeof(long));
02159         }
02160         data->filepos_array[data->num_frames] = ftell(data->file);
02161         data->num_frames++;
02162       }
02163     }
02164 
02165     if (appendedCalculations) {
02166       std::cout << "orcaplugin) Found multiple appended gradient calculations: " << data->num_frames << std::endl;
02167       pass_keyline(data->file, "FINAL SINGLE POINT ENERGY", NULL);
02168       data->end_of_traj = ftell(data->file);
02169       fseek(data->file, filepos, SEEK_SET);
02170 
02171       data->qm_timestep = (qm_timestep_t *)calloc(data->num_frames,
02172                                                   sizeof(qm_timestep_t));
02173       memset(data->qm_timestep, 0, data->num_frames*sizeof(qm_timestep_t));
02174     } else {
02175       data->num_frames = 1;
02176       pass_keyline(data->file, "FINAL SINGLE POINT ENERGY", NULL);
02177       data->end_of_traj = ftell(data->file);
02178 
02179       /* Allocate memory for the frame */
02180       data->qm_timestep = (qm_timestep_t *)calloc(data->num_frames, sizeof(qm_timestep_t));
02181       memset(data->qm_timestep, 0, sizeof(qm_timestep_t));
02182     }
02183     return TRUE;
02184   }
02185   else if (data->runtype == MOLFILE_RUNTYPE_OPTIMIZE) {
02186     std::cout << "orcaplugin) Reading trajectory of optimization." << std::endl;
02187     rewind(data->file);
02188     goto_keyline(data->file, "Geometry Optimization Run", NULL);
02189   }
02190   else {
02191     std::cout << "orcaplugin) Jobtype not supported for trajectory reading." << std::endl;
02192     return FALSE;
02193   }
02194 
02195   // printf("orcaplugin) Analyzing trajectory...\n");
02196   data->status = MOLFILE_QMSTATUS_UNKNOWN;
02197 
02198   while (TRUE) {
02199     if (!fgets(buffer, sizeof(buffer), data->file)) break;
02200     line = trimleft(buffer);
02201 
02202     std::string l(line);
02203     if (l.find("GEOMETRY OPTIMIZATION CYCLE") != std::string::npos && data->runtype==MOLFILE_RUNTYPE_OPTIMIZE) {
02204       // std::cout << l << std::endl;
02205       if (data->num_frames > 0) {
02206         data->filepos_array = (long*)realloc(data->filepos_array, (data->num_frames+1)*sizeof(long));
02207       }
02208       data->filepos_array[data->num_frames] = ftell(data->file);
02209       data->num_frames++;
02210     }
02211     else if (l.find("THE OPTIMIZATION HAS CONVERGED") != std::string::npos) {
02212       printf("orcaplugin) ==== End of trajectory (%d frames) ====\n", data->num_frames);
02213       data->status = MOLFILE_QMSTATUS_OPT_CONV;
02214     }
02215     else if (data->status == MOLFILE_QMSTATUS_OPT_CONV) {
02216       if(l.find("FINAL ENERGY EVALUATION AT THE STATIONARY POINT") != std::string::npos) {
02217         if (data->num_frames > 0) {
02218           data->filepos_array = (long*)realloc(data->filepos_array, (data->num_frames+1)*sizeof(long));
02219         }
02220         std::cout << "orcaplugin) found equilibrium geometry." << std::endl;
02221         data->filepos_array[data->num_frames] = ftell(data->file);
02222         data->num_frames++;
02223         goto_keyline(data->file, "TOTAL RUN TIME", NULL);
02224         break;
02225       }
02226     }
02227   }
02228 
02229   data->end_of_traj = ftell(data->file);
02230   fseek(data->file, filepos, SEEK_SET);
02231 
02232   if (data->status == MOLFILE_QMSTATUS_UNKNOWN) {
02233     /* We didn't find any of the regular key strings,
02234      * the run was most likely broken off and we have an
02235      * incomplete file. */
02236     data->status = MOLFILE_QMSTATUS_FILE_TRUNCATED;
02237   }
02238 
02239 
02240   /* Allocate memory for all frames */
02241   data->qm_timestep = (qm_timestep_t *)calloc(data->num_frames,
02242                                               sizeof(qm_timestep_t));
02243   memset(data->qm_timestep, 0, data->num_frames*sizeof(qm_timestep_t));
02244 
02245 
02246   if (data->status == MOLFILE_QMSTATUS_SCF_NOT_CONV ||
02247       data->status == MOLFILE_QMSTATUS_FILE_TRUNCATED) {
02248     return FALSE;
02249   }
02250 
02251   return TRUE;
02252 }
02253 
02254 
02255 
02256 /***********************************************************
02257  * Provide QM metadata for next timestep.
02258  * This actually triggers reading the entire next timestep
02259  * since we have to parse the whole timestep anyway in order
02260  * to get the metadata. So we store the read data locally
02261  * and hand them to VMD when requested by read_timestep().
02262  *
02263  ***********************************************************/
02264 static int read_qm_timestep_metadata(void *mydata,
02265                                     molfile_qm_timestep_metadata_t *meta) {
02266   int have = 0;
02267 
02268   qmdata_t *data = (qmdata_t *)mydata;
02269 
02270   meta->count = -1; /* Don't know the number of frames yet */
02271 
02272   if (data->num_frames_read > data->num_frames_sent) {
02273     have = 1;
02274   }
02275   else if (data->num_frames_read < data->num_frames) {
02276     printf("orcaplugin) Probing timestep %d\n", data->num_frames_read);
02277 
02278     have = get_traj_frame(data, data->atoms, data->numatoms);
02279   }
02280 
02281   if (have) {
02282     // std::cout << "have frame" << std::endl;
02283     int i;
02284     qm_timestep_t *cur_ts;
02285 
02286     /* get a pointer to the current qm timestep */
02287     cur_ts = data->qm_timestep+data->num_frames_sent;
02288 
02289     // std::cout << "numwave: " << cur_ts->numwave << std::endl;
02290 
02291     for (i=0; (i<MOLFILE_MAXWAVEPERTS && i<cur_ts->numwave); i++) {
02292       meta->num_orbitals_per_wavef[i] = cur_ts->wave[i].num_orbitals;
02293       meta->has_occup_per_wavef[i]    = cur_ts->wave[i].has_occup;
02294       meta->has_orben_per_wavef[i]    = cur_ts->wave[i].has_orben;
02295       // std::cout << "occ: " << cur_ts->wave[i].has_occup << std::endl;
02296       // std::cout << "energy: " << cur_ts->wave[i].has_orben << std::endl;
02297     }
02298     meta->wavef_size      = data->wavef_size;
02299     meta->num_wavef       = cur_ts->numwave;
02300     meta->num_scfiter     = cur_ts->num_scfiter;
02301     meta->num_charge_sets = cur_ts->have_mulliken +
02302       cur_ts->have_lowdin + cur_ts->have_esp;
02303     if (cur_ts->gradient) meta->has_gradient = TRUE;
02304 
02305   } else {
02306     // std::cout << "not have frame" << std::endl;
02307     meta->has_gradient = FALSE;
02308     meta->num_scfiter  = 0;
02309     meta->num_orbitals_per_wavef[0] = 0;
02310     meta->has_occup_per_wavef[0] = FALSE;
02311     meta->num_wavef = 0;
02312     meta->wavef_size = 0;
02313     meta->num_charge_sets = 0;
02314     data->trajectory_done = TRUE;
02315   }
02316 
02317   return MOLFILE_SUCCESS;
02318 }
02319 
02320 
02321 
02322 
02323 /*************************************************************
02324  *
02325  * plugin registration
02326  *
02327  **************************************************************/
02328 static molfile_plugin_t plugin;
02329 
02330 VMDPLUGIN_API int VMDPLUGIN_init(void) {
02331   memset(&plugin, 0, sizeof(molfile_plugin_t));
02332   plugin.abiversion = vmdplugin_ABIVERSION;
02333   plugin.type = MOLFILE_PLUGIN_TYPE;
02334   plugin.name = "orca";
02335   plugin.prettyname = "Orca";
02336   plugin.author = "Maximilian Scheurer, Michael F. Herbst, Marcelo Melo, Julio Maia, John Stone, M Spivak";
02337   plugin.majorv = 0;
02338   plugin.minorv = 2;
02339   plugin.is_reentrant = VMDPLUGIN_THREADUNSAFE;
02340   plugin.filename_extension = "orca";
02341   plugin.open_file_read = open_orca_read;
02342   plugin.read_structure = read_orca_structure;
02343   plugin.close_file_read = close_orca_read;
02344   //
02345   plugin.read_qm_metadata = read_orca_metadata;
02346   plugin.read_qm_rundata  = read_orca_rundata;
02347 
02348 #if vmdplugin_ABIVERSION > 11
02349   plugin.read_timestep_metadata    = read_timestep_metadata;
02350   plugin.read_qm_timestep_metadata = read_qm_timestep_metadata;
02351   plugin.read_timestep = read_timestep;
02352 #endif
02353 
02354   return VMDPLUGIN_SUCCESS;
02355 }
02356 
02357 VMDPLUGIN_API int VMDPLUGIN_register(void *v, vmdplugin_register_cb cb) {
02358   (*cb)(v, (vmdplugin_t *)&plugin);
02359   return VMDPLUGIN_SUCCESS;
02360 }
02361 
02362 VMDPLUGIN_API int VMDPLUGIN_fini(void) {
02363   return VMDPLUGIN_SUCCESS;
02364 }
02365 
02366 
02367 #if vmdplugin_ABIVERSION > 11
02368 
02369 /***********************************************************
02370  * Provide non-QM metadata for next timestep.
02371  * Required by the plugin interface.
02372  ***********************************************************/
02373 static int read_timestep_metadata(void *mydata,
02374                                   molfile_timestep_metadata_t *meta) {
02375   meta->count = -1;
02376   meta->has_velocities = 0;
02377 
02378   return MOLFILE_SUCCESS;
02379 }
02380 /***********************************************************
02381  *
02382  * This function provides the data of the next timestep.
02383  * Here we actually don't read the data from file, that had
02384  * to be done already upon calling read_timestep_metadata().
02385  * Instead we copy the stuff from the local data structure
02386  * into the one's provided by VMD.
02387  *
02388  ***********************************************************/
02389 static int read_timestep(void *mydata, int natoms,
02390        molfile_timestep_t *ts, molfile_qm_metadata_t *qm_metadata,
02391        molfile_qm_timestep_t *qm_ts)
02392 {
02393   qmdata_t *data = (qmdata_t *)mydata;
02394   qm_timestep_t *cur_ts;
02395   int offset;
02396   int i = 0;
02397   int num_charge_sets = 0;
02398 
02399   if (data->trajectory_done == TRUE) {
02400     printf("orcaplugin) Trajectory done.\n");
02401     return MOLFILE_ERROR;
02402   }
02403 
02404   /* copy the coordinates */
02405   for (i=0; i<natoms; i++) {
02406     ts->coords[3*i  ] = data->atoms[i].x;
02407     ts->coords[3*i+1] = data->atoms[i].y;
02408     ts->coords[3*i+2] = data->atoms[i].z;
02409     // printf("x: %f y: %f z: %f\n", data->atoms[i].x, data->atoms[i].y, data->atoms[i].z);
02410   }
02411 
02412   /* get a convenient pointer to the current qm timestep */
02413   cur_ts = data->qm_timestep+data->num_frames_sent;
02414   //
02415   // /* store the SCF energies */
02416   for (i=0; i<cur_ts->num_scfiter; i++) {
02417      qm_ts->scfenergies[i] = cur_ts->scfenergies[i];
02418   }
02419   //
02420   // /* store gradients */
02421   // if (cur_ts->gradient) {
02422   //   for (i=0; i<3*natoms; i++) {
02423   //     qm_ts->gradient[i] = cur_ts->gradient[i];
02424   //   }
02425   // }
02426   //
02427   // /* store charge sets*/
02428   if (cur_ts->have_mulliken) {
02429     offset = num_charge_sets*data->numatoms;
02430     for (i=0; i<data->numatoms; i++) {
02431       qm_ts->charges[offset+i] = cur_ts->mulliken_charges[i];
02432     }
02433     qm_ts->charge_types[num_charge_sets] = MOLFILE_QMCHARGE_MULLIKEN;
02434     num_charge_sets++;
02435   }
02436   //
02437   // if (cur_ts->have_lowdin) {
02438   //   offset = num_charge_sets*data->numatoms;
02439   //   for (i=0; i<data->numatoms; i++) {
02440   //     qm_ts->charges[offset+i] = cur_ts->lowdin_charges[i];
02441   //   }
02442   //   qm_ts->charge_types[num_charge_sets] = MOLFILE_QMCHARGE_LOWDIN;
02443   //   num_charge_sets++;
02444   // }
02445   // if (cur_ts->have_esp) {
02446   //   offset = num_charge_sets*data->numatoms;
02447   //   for (i=0; i<data->numatoms; i++) {
02448   //     qm_ts->charges[offset+i] = cur_ts->esp_charges[i];
02449   //   }
02450   //   qm_ts->charge_types[num_charge_sets] = MOLFILE_QMCHARGE_ESP;
02451   //   num_charge_sets++;
02452   // }
02453   //
02454   //
02455   /* store the wave function and orbital energies */
02456   if (cur_ts->wave) {
02457     std::cout << "orcaplugin) Have wavefunctions: " << cur_ts->numwave << " in frame: " << data->num_frames_sent << std::endl;
02458     for (i=0; i<cur_ts->numwave; i++) {
02459       qm_wavefunction_t *wave = &cur_ts->wave[i];
02460       qm_ts->wave[i].type         = wave->type;
02461       qm_ts->wave[i].spin         = wave->spin;
02462       qm_ts->wave[i].excitation   = wave->exci;
02463       qm_ts->wave[i].multiplicity = wave->mult;
02464       qm_ts->wave[i].energy       = wave->energy;
02465       strncpy(qm_ts->wave[i].info, wave->info, MOLFILE_BUFSIZ);
02466 
02467       if (wave->wave_coeffs) {
02468         memcpy(qm_ts->wave[i].wave_coeffs, wave->wave_coeffs, wave->num_orbitals*data->wavef_size*sizeof(float));
02469       }
02470       if (wave->orb_energies) {
02471         memcpy(qm_ts->wave[i].orbital_energies, wave->orb_energies,
02472                wave->num_orbitals*sizeof(float));
02473       }
02474       if (wave->has_occup) {
02475         memcpy(qm_ts->wave[i].occupancies, wave->orb_occupancies,
02476                wave->num_orbitals*sizeof(float));
02477       }
02478     }
02479   }
02480   //
02481   if (data->runtype == MOLFILE_RUNTYPE_ENERGY ||
02482       data->runtype == MOLFILE_RUNTYPE_HESSIAN) {
02483     /* We have only a single point */
02484     data->trajectory_done = TRUE;
02485   }
02486 
02487   data->num_frames_sent++;
02488 
02489   return MOLFILE_SUCCESS;
02490 }
02491 #endif
02492 
02493 /*****************************************************
02494  *
02495  * provide VMD with the sizes of the QM related
02496  * data structure arrays that need to be made
02497  * available
02498  *
02499  *****************************************************/
02500 static int read_orca_metadata(void *mydata,
02501     molfile_qm_metadata_t *metadata) {
02502 
02503   qmdata_t *data = (qmdata_t *)mydata;
02504 
02505   if (data->runtype == MOLFILE_RUNTYPE_HESSIAN) {
02506     metadata->ncart = (3*data->numatoms);
02507     metadata->nimag = data->nimag;
02508 
02509     if (data->have_internals) {
02510       metadata->nintcoords = data->nintcoords;
02511     } else {
02512       metadata->nintcoords = 0;
02513     }
02514   }
02515   else {
02516     metadata->ncart = 0;
02517     metadata->nimag = 0;
02518     metadata->nintcoords = 0;
02519   }
02520 
02521   /* orbital data */
02522   metadata->num_basis_funcs = data->num_basis_funcs;
02523   metadata->num_basis_atoms = data->num_basis_atoms;
02524   metadata->num_shells      = data->num_shells;
02525   metadata->wavef_size      = data->wavef_size;
02526 
02527 #if vmdplugin_ABIVERSION > 11
02528   /* system and run info */
02529   metadata->have_sysinfo = 1;
02530 
02531   /* hessian info */
02532   metadata->have_carthessian = data->have_cart_hessian;
02533   metadata->have_inthessian  = data->have_int_hessian;
02534 
02535   /* normal mode info */
02536   metadata->have_normalmodes = data->have_normal_modes;
02537 #endif
02538 
02539   return MOLFILE_SUCCESS;
02540 }
02541 
02542 
02543 /******************************************************
02544  *
02545  * Provide VMD with the static (i.e. non-trajectory)
02546  * data. That means we are filling the molfile_plugin
02547  * data structures.
02548  *
02549  ******************************************************/
02550 static int read_orca_rundata(void *mydata,
02551                                molfile_qm_t *qm_data) {
02552 
02553   qmdata_t *data = (qmdata_t *)mydata;
02554   int i, j;
02555   int ncart;
02556   molfile_qm_hessian_t *hessian_data = &qm_data->hess;
02557   molfile_qm_basis_t   *basis_data   = &qm_data->basis;
02558   molfile_qm_sysinfo_t *sys_data     = &qm_data->run;
02559 
02560   /* fill in molfile_qm_hessian_t */
02561   if (data->runtype == MOLFILE_RUNTYPE_HESSIAN) {
02562     ncart = 3*data->numatoms;
02563 
02564     /* Hessian matrix in cartesian coordinates */
02565     if (data->have_cart_hessian) {
02566       for (i=0; i<ncart; i++) {
02567         for (j=0; j<=i; j++) {
02568           hessian_data->carthessian[ncart*i+j] = data->carthessian[ncart*i+j];
02569           hessian_data->carthessian[ncart*j+i] = data->carthessian[ncart*i+j];
02570         }
02571       }
02572     }
02573 
02574     /* Hessian matrix in internal coordinates */
02575     if (data->have_int_hessian) {
02576       for (i=0; i<(data->nintcoords)*(data->nintcoords); i++) {
02577         hessian_data->inthessian[i] = data->inthessian[i];
02578       }
02579     }
02580 
02581     /* wavenumbers, intensities, normal modes */
02582     if (data->have_normal_modes) {
02583       for (i=0; i<ncart*ncart; i++) {
02584         hessian_data->normalmodes[i] = data->normal_modes[i];
02585       }
02586       for (i=0; i<ncart; i++) {
02587         hessian_data->wavenumbers[i] = data->wavenumbers[i];
02588         hessian_data->intensities[i] = data->intensities[i];
02589       }
02590     }
02591 
02592     /* imaginary modes */
02593     for (i=0; i<data->nimag; i++) {
02594       /*printf("imag_modes[%d]=%d\n", i, data->imag_modes[i]);*/
02595       hessian_data->imag_modes[i] = data->imag_modes[i];
02596     }
02597   }
02598 
02599   /* fill in molfile_qm_sysinfo_t */
02600   sys_data->runtype = data->runtype;
02601   sys_data->scftype = data->scftype;
02602   sys_data->nproc   = data->nproc;
02603   sys_data->num_electrons  = data->num_electrons;
02604   sys_data->totalcharge    = data->totalcharge;
02605   sys_data->num_occupied_A = data->num_occupied_A;
02606   sys_data->num_occupied_B = data->num_occupied_B;
02607   sys_data->status         = data->status;
02608 
02609 
02610   strncpy(sys_data->basis_string, data->basis_string,
02611           sizeof(sys_data->basis_string));
02612 
02613   sys_data->memory = 0; /* XXX fixme */
02614 
02615   strncpy(sys_data->runtitle, data->runtitle, sizeof(sys_data->runtitle));
02616   strncpy(sys_data->geometry, data->geometry, sizeof(sys_data->geometry));
02617   strncpy(sys_data->version_string, data->version_string,
02618           sizeof(sys_data->version_string));
02619 
02620 #if vmdplugin_ABIVERSION > 11
02621   /* fill in molfile_qm_basis_t */
02622   if (data->num_basis_funcs) {
02623     for (i=0; i<data->num_basis_atoms; i++) {
02624       basis_data->num_shells_per_atom[i] = data->num_shells_per_atom[i];
02625       basis_data->atomic_number[i] = data->atomicnum_per_basisatom[i];
02626     }
02627 
02628     for (i=0; i<data->num_shells; i++) {
02629       basis_data->num_prim_per_shell[i] = data->num_prim_per_shell[i];
02630       basis_data->shell_types[i] = data->shell_types[i];
02631     }
02632 
02633     for (i=0; i<2*data->num_basis_funcs; i++) {
02634       basis_data->basis[i] = data->basis[i];
02635     }
02636 
02637     for (i=0; i<3*data->wavef_size; i++) {
02638       basis_data->angular_momentum[i] = data->angular_momentum[i];
02639     }
02640     // std::cout << "free data->angular_momentum" << std::endl;
02641     // free(data->angular_momentum);
02642   }
02643 #endif
02644 
02645   return MOLFILE_SUCCESS;
02646 }
02647 
02648 
02649 
02650 
02651 /**********************************************************
02652  *
02653  * close file and free memory
02654  *
02655  **********************************************************/
02656  // move to top later...
02657 static void close_orca_read(void *mydata) {
02658   printf("Freeing memory.\n");
02659 
02660   qmdata_t *data = (qmdata_t *)mydata;
02661   int i, j;
02662   fclose(data->file);
02663 
02664   free(data->atoms);
02665   free(data->basis);
02666   free(data->shell_types);
02667   free(data->atomicnum_per_basisatom);
02668   free(data->num_shells_per_atom);
02669   free(data->num_prim_per_shell);
02670   free(data->bonds);
02671   free(data->angles);
02672   free(data->dihedrals);
02673   free(data->impropers);
02674   free(data->internal_coordinates);
02675   free(data->bond_force_const);
02676   free(data->angle_force_const);
02677   free(data->dihedral_force_const);
02678   free(data->improper_force_const);
02679   free(data->inthessian);
02680   free(data->carthessian);
02681   free(data->wavenumbers);
02682   free(data->intensities);
02683   free(data->normal_modes);
02684   free(data->imag_modes);
02685   free(data->angular_momentum);
02686   data->angular_momentum = NULL;
02687   free(data->filepos_array);
02688 
02689   if (data->basis_set) {
02690     for (i=0; i<data->num_basis_atoms; i++) {
02691       // printf("Freeing basis set of atom %d\n", i);
02692       for (j=0; j<data->basis_set[i].numshells; j++) {
02693         // printf("Freeing shell primitives %d\n", j);
02694         // printf("--- Address of prim is %p\n", (void *)data->basis_set[i].shell[j].prim);
02695         free(data->basis_set[i].shell[j].prim);
02696         data->basis_set[i].shell[j].prim = NULL;
02697       }
02698       // printf("- Address of shell is %p\n", (void *)data->basis_set[i].shell);
02699       free(data->basis_set[i].shell);
02700       data->basis_set[i].shell = NULL;
02701       // printf("- Address of shell is %p\n", (void *)data->basis_set[i].shell);
02702     }
02703     free(data->basis_set);
02704     data->basis_set = NULL;
02705   }
02706 
02707   for (i=0; i<data->num_frames; i++) {
02708     free(data->qm_timestep[i].scfenergies);
02709     free(data->qm_timestep[i].gradient);
02710     free(data->qm_timestep[i].mulliken_charges);
02711     free(data->qm_timestep[i].lowdin_charges);
02712     free(data->qm_timestep[i].esp_charges);
02713     for (j=0; j<data->qm_timestep[i].numwave; j++) {
02714       free(data->qm_timestep[i].wave[j].wave_coeffs);
02715       free(data->qm_timestep[i].wave[j].orb_energies);
02716       free(data->qm_timestep[i].wave[j].orb_occupancies);
02717     }
02718     free(data->qm_timestep[i].wave);
02719   }
02720   free(data->qm_timestep);
02721   free(data->format_specific_data);
02722   free(data);
02723 }
02724 
02725 
02726 
02727 #ifdef TEST_ORCAPLUGIN
02728 
02729 static void print_input_data(qmdata_t *data) {
02730   int i, j, k;
02731   int primcount=0;
02732   int shellcount=0;
02733 
02734   printf("\nDATA READ FROM FILE:\n\n");
02735   printf(" %10s WORDS OF MEMORY AVAILABLE\n", data->memory);
02736   printf("\n");
02737   printf("     BASIS OPTIONS\n");
02738   printf("     -------------\n");
02739   printf("%s\n", data->basis_string);
02740   printf("\n\n\n");
02741   printf("     RUN TITLE\n");
02742   printf("     ---------\n");
02743   printf(" %s\n", data->runtitle);
02744   printf("\n");
02745   printf(" THE POINT GROUP OF THE MOLECULE IS %s\n", "XXX");
02746   printf(" THE ORDER OF THE PRINCIPAL AXIS IS %5i\n", 0);
02747   printf("\n");
02748   printf(" YOUR FULLY SUBSTITUTED Z-MATRIX IS\n");
02749   printf("\n");
02750   printf(" THE MOMENTS OF INERTIA ARE (AMU-ANGSTROM**2)\n");
02751   printf(" IXX=%10.3f   IYY=%10.3f   IZZ=%10.3f\n", 0.0, 0.0, 0.0);
02752   printf("\n");
02753   printf(" ATOM      ATOMIC                      COORDINATES (BOHR)\n");
02754   printf("           CHARGE         X                   Y                   Z\n");
02755   for (i=0; i<data->numatoms; i++) {
02756     printf(" %-8s %6d", data->atoms[i].type, data->atoms[i].atomicnum);
02757 
02758     printf("%17.10f",   ANGS_TO_BOHR*data->atoms[i].x);
02759     printf("%20.10f",   ANGS_TO_BOHR*data->atoms[i].y);
02760     printf("%20.10f\n", ANGS_TO_BOHR*data->atoms[i].z);
02761   }
02762   printf("\n");
02763   printf("     ATOMIC BASIS SET\n");
02764   printf("     ----------------\n");
02765   printf(" THE CONTRACTED PRIMITIVE FUNCTIONS HAVE BEEN UNNORMALIZED\n");
02766   printf(" THE CONTRACTED BASIS FUNCTIONS ARE NOW NORMALIZED TO UNITY\n");
02767   printf("\n");
02768   printf("  SHELL TYPE  PRIMITIVE        EXPONENT          CONTRACTION COEFFICIENT(S)\n");
02769   printf("\n");
02770 
02771 #if 0
02772   for (i=0; i<data->numatoms; i++) {
02773     printf("%-8s\n\n", data->atoms[i].type);
02774     printf("\n");
02775     printf("nshells=%d\n", data->num_shells_per_atom[i]);
02776 
02777     for (j=0; j<data->num_shells_per_atom[i]; j++) {
02778       printf("nprim=%d\n", data->num_prim_per_shell[shellcount]);
02779 
02780       for (k=0; k<data->num_prim_per_shell[shellcount]; k++) {
02781         printf("%6d   %d %7d %22f%22f\n", j, data->shell_types[shellcount],
02782                primcount+1, data->basis[2*primcount], data->basis[2*primcount+1]);
02783         primcount++;
02784       }
02785 
02786       printf("\n");
02787       shellcount++;
02788     }
02789   }
02790 #endif
02791   printf("orcaplugin) =================================================================\n");
02792   for (i=0; i<data->num_basis_atoms; i++) {
02793     printf("%-8s (%10s)\n\n", data->atoms[i].type, data->basis_set[i].name);
02794     printf("\n");
02795 
02796     for (j=0; j<data->basis_set[i].numshells; j++) {
02797 
02798       for (k=0; k<data->basis_set[i].shell[j].numprims; k++) {
02799         printf("%6d   %d %7d %22f%22f\n", j,
02800                data->basis_set[i].shell[j].type,
02801                primcount+1,
02802                data->basis_set[i].shell[j].prim[k].exponent,
02803                data->basis_set[i].shell[j].prim[k].contraction_coeff);
02804         primcount++;
02805       }
02806 
02807       printf("\n");
02808       shellcount++;
02809     }
02810   }
02811   printf("\n");
02812   printf(" TOTAL NUMBER OF BASIS SET SHELLS             =%5d\n", data->num_shells);
02813   printf(" NUMBER OF CARTESIAN GAUSSIAN BASIS FUNCTIONS =%5d\n", data->wavef_size);
02814   printf(" NUMBER OF ELECTRONS                          =%5d\n", data->num_electrons);
02815   printf(" CHARGE OF MOLECULE                           =%5d\n", data->totalcharge);
02816   printf(" SPIN MULTIPLICITY                            =%5d\n", data->multiplicity);
02817   printf(" NUMBER OF OCCUPIED ORBITALS (ALPHA)          =%5d\n", data->num_occupied_A);
02818   printf(" NUMBER OF OCCUPIED ORBITALS (BETA )          =%5d\n", data->num_occupied_B);
02819   printf(" TOTAL NUMBER OF ATOMS                        =%5i\n", data->numatoms);
02820   printf("\n");
02821 }
02822 
02823 #endif
02824 
02825 

Generated on Fri Dec 6 03:09:50 2024 for VMD Plugins (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002