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

QMTimestep.C

Go to the documentation of this file.
00001 /***************************************************************************
00002  *cr                                                                       
00003  *cr            (C) Copyright 1995-2019 The Board of Trustees of the           
00004  *cr                        University of Illinois                       
00005  *cr                         All Rights Reserved                        
00006  *cr                                                                   
00007  ***************************************************************************/
00008 
00009 /***************************************************************************
00010  * RCS INFORMATION:
00011  *
00012  *      $RCSfile: QMTimestep.C,v $
00013  *      $Author: johns $        $Locker:  $             $State: Exp $
00014  *      $Revision: 1.74 $       $Date: 2019/11/07 03:52:53 $
00015  *
00016  ***************************************************************************
00017  * DESCRIPTION:
00018  *
00019  * The QMTimestep and Wavefunction classes.
00020  * The QMTimestep class stores and manages all quantm chemistry related
00021  * data that are timestep dependent. These include gradients, charges,
00022  * SCF energies. Most importantly we also consider wave function
00023  * trajectories which enables the user to visualize orbital dynamics.
00024  * (None of the other visualization tools I'm aware of can do that :-)
00025  * 
00026  * Moreover, each timestep can have multiple different wave functions.
00027  * This is useful for UHF calculations where we have two sets of orbitals
00028  * with opposite spins or for calculations with multiple excited states. 
00029  * The number of existing wave function may vary between frames.
00030  * So we could have, for instance, a coordinate trajectory with orbitals
00031  * defined only for the last frame. Another typical case would be to have
00032  * a canonical wavefunction for each frame but for the last frame there
00033  * exists an additional set of localized orbitals.
00034  * The number of orbitals present for a given wavefunction may also
00035  * differ from frame to frame. 
00036  *
00037  * In order to identify a wavefunction throughout the trajectory we
00038  * assign a unique wavefunction ID that is based on its 'signature'
00039  * (type, spin, excitation, multiplicity and an optional info string).
00040  * A list with the signatures of all different wavefunction occuring
00041  * over the course of the trajectory is kept in the QMData class.
00042  *
00043  ***************************************************************************/
00044 
00045 #include <math.h>
00046 #include <stdio.h>
00047 #include "Inform.h"
00048 #include "QMTimestep.h"
00049 #include "QMData.h"
00050 #include "Molecule.h"
00051 
00052 #define ANGMOM_X  0
00053 #define ANGMOM_Y  1
00054 #define ANGMOM_Z  2
00055 
00056 //#define DEBUG 1
00057 
00058 Wavefunction::Wavefunction()
00059 {
00060   idtag      = 0;
00061   type       = 0;
00062   spin       = 0;
00063   excitation = 0;
00064   multiplicity = 0;
00065   num_orbitals = 0;
00066   num_coeffs = 0;
00067   energy = 0.0;
00068   wave_coeffs  = NULL;
00069   orb_energies = NULL;
00070   occupancies  = NULL;
00071   orb_ids      = NULL;
00072   orb_id2index = NULL;
00073 }
00074 
00075 Wavefunction::Wavefunction(const Wavefunction& wf)
00076 {
00077   wave_coeffs  = NULL;
00078   orb_energies = NULL;
00079   occupancies  = NULL;
00080   orb_ids      = NULL;
00081   orb_id2index = NULL;
00082   *this = wf;
00083 }
00084 
00085 Wavefunction::Wavefunction(int ncoeffs,
00086                            int norbitals, 
00087                            const float *coeffs,
00088                            const float *energies,
00089                            const float *occ,
00090                            const int   *orbids,
00091                            double _energy,
00092                            int _idtag,
00093                            int _type,
00094                            int _spin,
00095                            int _excitation,
00096                            int _multiplicity,
00097                            char *infostr) :
00098   idtag     (_idtag),
00099   type      (_type),
00100   excitation(_excitation),
00101   multiplicity(_multiplicity),
00102   num_orbitals(norbitals),
00103   num_coeffs  (ncoeffs),
00104   energy      (_energy),
00105   wave_coeffs (NULL),
00106   orb_energies(NULL),
00107   occupancies (NULL),
00108   orb_ids     (NULL),
00109   orb_id2index(NULL)
00110 {
00111   strncpy(info, infostr, QMDATA_BUFSIZ - 1);
00112 
00113   set_coeffs(coeffs, norbitals, ncoeffs);
00114   set_orbenergies(energies, norbitals);
00115   set_occupancies(occ, norbitals);
00116   set_orbids(orbids, norbitals);
00117 }
00118 
00119 Wavefunction& Wavefunction::operator=(const Wavefunction& wf) {
00120   idtag        = wf.idtag;
00121   type         = wf.type;
00122   spin         = wf.spin;
00123   excitation   = wf.excitation;
00124   multiplicity = wf.multiplicity;
00125   strncpy(info, wf.info, QMDATA_BUFSIZ);
00126 
00127   num_orbitals = wf.num_orbitals;
00128   num_coeffs   = wf.num_coeffs;
00129   energy       = wf.energy;
00130 
00131   if (orb_energies) delete [] orb_energies;
00132   if (orb_ids)      delete [] orb_ids;
00133   if (orb_id2index) delete [] orb_id2index;
00134   if (wave_coeffs)  delete [] wave_coeffs;
00135   if (occupancies)  delete [] occupancies;
00136   wave_coeffs  = NULL;
00137   orb_energies = NULL;
00138   occupancies  = NULL;
00139   orb_ids      = NULL;
00140   orb_id2index = NULL;
00141 
00142   if (wf.orb_energies) {
00143     orb_energies = new float[num_orbitals];
00144     memcpy(orb_energies, wf.orb_energies, num_orbitals * sizeof(float));
00145   }
00146 
00147   if (wf.wave_coeffs) {
00148     wave_coeffs = new float[num_orbitals * num_coeffs];
00149     memcpy(wave_coeffs, wf.wave_coeffs, num_orbitals * num_coeffs * sizeof(float));
00150   }
00151 
00152   if (wf.occupancies) {
00153     occupancies = new float[num_orbitals];
00154     memcpy(occupancies, wf.occupancies, num_orbitals * sizeof(int));
00155   }
00156 
00157   if (wf.orb_ids) {
00158     orb_ids = new int[num_orbitals];
00159     memcpy(orb_ids, wf.orb_ids, num_orbitals * sizeof(int));
00160   }
00161 
00162   if (wf.orb_id2index) {
00163     orb_id2index = new int[num_coeffs];
00164     memcpy(orb_id2index, wf.orb_id2index, num_coeffs * sizeof(int));
00165   }
00166 
00167   return *this;
00168 }
00169 
00170 // Move the data over from the given wavefunction wf
00171 // and set the pointers in wf to NULL.
00172 // This avoids copying the arrays.
00173 void Wavefunction::movefrom(Wavefunction& wf) {
00174   idtag        = wf.idtag;
00175   type         = wf.type;
00176   spin         = wf.spin;
00177   excitation   = wf.excitation;
00178   multiplicity = wf.multiplicity;
00179   strncpy(info, wf.info, QMDATA_BUFSIZ);
00180 
00181   num_orbitals = wf.num_orbitals;
00182   num_coeffs   = wf.num_coeffs;
00183   energy       = wf.energy;
00184 
00185   wave_coeffs  = wf.wave_coeffs;
00186   orb_energies = wf.orb_energies;
00187   occupancies  = wf.occupancies;
00188   orb_ids      = wf.orb_ids;
00189   orb_id2index = wf.orb_id2index;
00190   wf.wave_coeffs  = NULL;
00191   wf.orb_energies = NULL;
00192   wf.occupancies  = NULL;
00193   wf.orb_ids      = NULL;
00194   wf.orb_id2index = NULL;
00195 }
00196 
00197 Wavefunction::~Wavefunction()
00198 {
00199   if (wave_coeffs)  delete [] wave_coeffs;
00200   if (orb_energies) delete [] orb_energies;
00201   if (occupancies)  delete [] occupancies;
00202   if (orb_ids)      delete [] orb_ids;
00203   if (orb_id2index) delete [] orb_id2index;
00204 }
00205 
00206 #if 0
00207 const float* Wavefunction::get_coeffs(int orb)
00208 {
00209   if (!wave_coeffs || orb<0 || orb>=num_orbitals) return NULL;
00210   return wave_coeffs + orb*num_coeffs;
00211 }
00212 #endif
00213 
00214 
00215 // Return energy of the given orbital or zero if that orbital
00216 // doesn't exist.
00217 float Wavefunction::get_orbitalenergy(int orb) const
00218 {
00219   if (orb_energies && orb>=0 && orb<num_orbitals)
00220     return orb_energies[orb];
00221   else
00222     return 0.f;
00223 }
00224 
00225 
00226 // Set the wavefunction coefficients
00227 void Wavefunction::set_coeffs(const float *wfn, int norbitals, int wavef_size) 
00228 {
00229   if (!wfn || !norbitals || !wavef_size) return;
00230   num_orbitals = norbitals;
00231   num_coeffs = wavef_size;
00232   
00233   wave_coeffs = new float[num_orbitals*num_coeffs];
00234   memcpy(wave_coeffs, wfn, num_orbitals*num_coeffs*sizeof(float));
00235 }
00236 
00237 
00238 // Set the orbital energies
00239 void Wavefunction::set_orbenergies(const float *energies, int norbitals) 
00240 {
00241   if (!energies || !norbitals) return;
00242 
00243   if (num_orbitals < 1) 
00244     num_orbitals = norbitals;
00245 
00246   if (num_orbitals != norbitals) 
00247     msgWarn << "Mismatched number of orbitals in " << "QMTimestep::set_orbenergies()" << ": " 
00248             << norbitals << " != " << num_orbitals << sendmsg;
00249   
00250   orb_energies = new float[norbitals];
00251   memcpy(orb_energies, energies, norbitals*sizeof(float));
00252 }
00253 
00254 
00255 // Set the orbital occupancies
00256 void Wavefunction::set_occupancies(const float *occ, int norbitals) 
00257 {
00258   if (!occ || !norbitals) return;
00259 
00260   if (num_orbitals < 1) 
00261     num_orbitals = norbitals;
00262 
00263   if (num_orbitals != norbitals) 
00264     msgWarn << "Mismatched number of orbitals in " << "QMTimestep::set_occupancies()" << ": " 
00265             << norbitals << " != " << num_orbitals << sendmsg;
00266   
00267   occupancies = new float[norbitals];
00268   memcpy(occupancies, occ, norbitals*sizeof(float));
00269 }
00270 
00271 // Set orbital ID number array.
00272 // Assumed 1,2,3,...num_orbitals if orbids==NULL.      
00273 void Wavefunction::set_orbids(const int *orbids, int norbitals) 
00274 {
00275   if (!norbitals) return;
00276 
00277   if (num_orbitals < 1) 
00278     num_orbitals = norbitals;
00279 
00280   if (num_orbitals != norbitals) 
00281     msgWarn << "Mismatched number of orbitals in " << "QMTimestep::set_orbids()" << ": " 
00282             << norbitals << " != " << num_orbitals << sendmsg;
00283   
00284   int i;
00285   orb_ids = new int[norbitals];
00286   if (orbids) {
00287     memcpy(orb_ids, orbids, norbitals*sizeof(int));
00288   } else {
00289     for (i=0; i<num_orbitals; i++) {
00290       orb_ids[i] = i+1;
00291     }
00292   }
00293 
00294   orb_id2index = new int[num_coeffs+1];
00295   for (i=0; i<num_coeffs+1; i++) {
00296     orb_id2index[i] = -1;
00297   }
00298   for (i=0; i<num_orbitals; i++) {
00299     orb_id2index[orb_ids[i]] = i;
00300 //    printf("orb_id2index[%d]=%d\n", orb_ids[i], orb_id2index[orb_ids[i]]);
00301   }
00302 }
00303 
00304 
00305 // Get a string describing the wavefunction type
00306 void Wavefunction::get_typestr(char *&typestr) const {
00307   typestr = new char[64];
00308   switch (type) {
00309   case WAVE_CANON:
00310     strcpy(typestr, "Canonical");
00311     break;
00312   case WAVE_CINATUR:
00313     strcpy(typestr, "CI natural");
00314     break;
00315   case WAVE_GEMINAL:
00316     strcpy(typestr, "GVB geminal pairs");
00317     break;
00318   case WAVE_BOYS:
00319     strcpy(typestr, "Boys localized");
00320     break;
00321   case WAVE_RUEDEN:
00322     strcpy(typestr, "Ruedenberg localized");
00323     break;
00324   case WAVE_PIPEK:
00325     strcpy(typestr, "Pipek-Mezey localized");
00326     break;
00327   case WAVE_MCSCFOPT:
00328     strcpy(typestr, "MCSCF optimized");
00329     break;
00330   case WAVE_MCSCFNAT:
00331     strcpy(typestr, "MCSCF natural");
00332     break;
00333   default:
00334     strcpy(typestr, "Unknown");
00335   }
00336 }
00337 
00338 
00339 // Get orbital index for HOMO
00340 int Wavefunction::get_homo() const {
00341   if (!occupancies) return -1;
00342   if (!orb_energies && type!=WAVE_CANON) return -1;
00343 
00344   int i;
00345   int homo = -1;
00346   float maxorben = 0.f;
00347   if (orb_energies) maxorben = orb_energies[0];
00348   for (i=0; i<num_orbitals; i++) {
00349     int intocc = (int)floor(0.5f+occupancies[i]);
00350     if (intocc>0) {
00351       if (type==WAVE_CANON ||
00352           (orb_energies && orb_energies[i]>maxorben)) {
00353         homo = i;
00354       }
00355     }
00356   }
00357 
00358   return homo;
00359 }
00360 
00361 
00362 // Get number of double occupied orbitals
00363 int Wavefunction::get_num_occupied_double() const {
00364   if (!occupancies) return -1;
00365 
00366   int i;
00367   int numocc = 0;
00368   for (i=0; i<num_orbitals; i++) {
00369     int intocc = (int)floor(0.5f+occupancies[i]);
00370     if (intocc==2) numocc ++;
00371   }
00372   return numocc;
00373 }
00374 
00375 
00376 // Get number of single occupied orbitals
00377 int Wavefunction::get_num_occupied_single() const {
00378   if (!occupancies) return -1;
00379 
00380   int i;
00381   int numocc = 0;
00382   for (i=0; i<num_orbitals; i++) {
00383     int intocc = (int)floor(0.5f+occupancies[i]);
00384     if (intocc==1) numocc++;
00385   }
00386   return numocc;
00387 }
00388 
00389 
00390 // The array angular_momentum consists of 3*num_wave_f flags
00391 // determining which cartesian component of the angular
00392 // momentum each wavefunction coefficient corresponds to. 
00393 // Each triplet represents the exponents of the x, y, and z
00394 // components. I.e. (1, 0, 2) means xzz for an F shell.
00395 // Our inner loop in the orbital computation assumes an order
00396 // with the z-component changing fastest and x slowest, i.e
00397 // for a D- shell the order would be xx xy xz yy yz zz
00398 // (I hope I did that right :-)
00399 // This is a routine to sort the wavefunction coefficients
00400 // of each shell according to that scheme.
00401 void Wavefunction::sort_wave_coefficients(QMData *qmdata) {
00402   if (!wave_coeffs || !num_orbitals ||
00403       !qmdata->num_wave_f || !qmdata->num_basis) {
00404     return;
00405   }
00406 
00407   int atom, j;
00408   for (atom=0; atom<qmdata->get_num_atoms(); atom++) {
00409     for (j=0; j<qmdata->get_num_shells_per_atom()[atom]; j++) {
00410 
00411       const shell_t *shell = qmdata->get_basis(atom, j);
00412       if (!shell) {
00413         printf("sort_wave_coefficients(): NO SHELL %d %d\n", atom, j);
00414       }
00415       int shelltype = shell->type;
00416 
00417       // Sort the wavefunction coefficients of this shell
00418       // according to the angular momentum.
00419       // First we must sort by increasing exponent of the
00420       // z-component:
00421       sort_incr(qmdata, atom, j, ANGMOM_Z, 0, shell->num_cart_func);
00422 
00423       // Then we sort the coefficients for each z-component
00424       // by increasing exponent of the y-component.
00425       // The x-component needs no sorting since it is
00426       // dependent on y and z.
00427       int k, first=0;
00428       for (k=0; k<=shelltype; k++) {
00429         sort_incr(qmdata, atom, j, ANGMOM_Y, first, shelltype-k+1);
00430         first += shelltype-k+1;
00431       }
00432     }
00433   }
00434 }
00435 
00436 
00437 
00438 // Get density matrix D
00439 void Wavefunction::density_matrix(float *(&D)) const {
00440   D = new float[num_coeffs*num_coeffs];
00441 
00442   // XXX This won't work for unsorted orbitals
00443   int num_occ = get_num_occupied_double();
00444   printf("numocc=%d\n", num_occ);
00445   int i, j, k;
00446   for (i=0; i<num_coeffs; i++) {
00447     for (j=0; j<num_coeffs; j++) {
00448       float Dij = 0.f;
00449       for (k=0; k<num_occ; k++) {
00450         Dij += wave_coeffs[k*num_coeffs+i]*wave_coeffs[k*num_coeffs+j];
00451       }
00452       D[i*num_coeffs+j] = 2.f*Dij;
00453     }
00454   }
00455 }
00456 
00457 // Get density matrix D for atom <atomid>
00458 void Wavefunction::density_matrix(const QMData *qmdata, int atomid,
00459                                   float *(&D)) const {
00460   int first_coeff = qmdata->get_wave_offset(atomid, 0);
00461   int num_atom_coeffs = qmdata->get_num_wavecoeff_per_atom(atomid);
00462 
00463   D = new float[num_coeffs*num_atom_coeffs];
00464 
00465   // XXX This won't work for unsorted orbitals
00466   int num_occ = get_num_occupied_double();
00467   printf("numocc=%d\n", num_occ);
00468   int i, j, k;
00469   for (i=0; i<num_atom_coeffs; i++) {
00470     for (j=0; j<num_coeffs; j++) {
00471       float Dij = 0.f;
00472       for (k=0; k<num_occ; k++) {
00473         Dij += wave_coeffs[k*num_coeffs+first_coeff+i]*wave_coeffs[k*num_coeffs+first_coeff+j];
00474       }
00475       D[i*num_atom_coeffs+j] = 2.f*Dij;
00476     }
00477   }
00478 }
00479 
00480 
00481 // Get Population matrix P
00482 void Wavefunction::population_matrix(const float *S, float *(&P)) const {
00483   float *D;
00484   density_matrix(D);
00485 
00486   P = new float[num_coeffs*num_coeffs];
00487   int i;
00488   for (i=0; i<num_coeffs*num_coeffs; i++) {
00489     P[i] = S[i]*D[i];
00490   }
00491 
00492   delete [] D;
00493 }
00494 
00495 
00496 // Get Population matrix P for atom <atomid>
00497 void Wavefunction::population_matrix(const QMData *qmdata, int atomid,
00498                                      const float *S, float *(&P)) const {
00499   float *D;
00500   density_matrix(D);
00501 
00502   int first_coeff = qmdata->get_wave_offset(atomid, 0);
00503   int num_atom_coeffs = qmdata->get_num_wavecoeff_per_atom(atomid);
00504   P = new float[num_coeffs*num_atom_coeffs];
00505 
00506   int i, j;
00507   for (i=0; i<num_atom_coeffs; i++) {
00508     for (j=0; j<num_coeffs; j++) {
00509       P[i*num_coeffs+j] = S[first_coeff+i*num_coeffs+j]*D[first_coeff+i*num_coeffs+j];
00510     }
00511   }
00512 
00513   delete [] D;
00514 }
00515 
00516 
00517 // Sort array of indexes *idx according to the order obtained 
00518 // by sorting the integers in array *tag. Array *idx can then
00519 // be used to reorder any array according to the string tags.
00520 static void quicksort(const int *tag, int *A, int p, int r);
00521 static int  quickpart(const int *tag, int *A, int p, int r);
00522 
00523 
00524 // Sort the wavefunction coefficients of the specified shell
00525 // according to the increasing exponent of requested angular
00526 // momentum component.
00527 // Parameters:
00528 // comp:  0,1,2; sort according to the x,y,z component 
00529 // first: the array element where the sorting shall begin
00530 //        counted from the beginning of the shell
00531 // num:   the number of consecutive array elements to be
00532 //        sorted
00533 void Wavefunction::sort_incr(QMData *qmdata, int atom, int ishell,
00534                            int comp, int first, int num) {
00535   int i, orb;
00536   int wave_offset = qmdata->get_wave_offset(atom, ishell);
00537 
00538   // Initialize the index array;
00539   int *index = new int[num];
00540   int *powz  = new int[num];
00541   for (i=0; i<num; i++) {
00542     index[i] = i;
00543     powz[i] = qmdata->get_angular_momentum(atom, ishell, first+i, comp);
00544   }
00545 
00546   float *wave_f = wave_coeffs + wave_offset + first;
00547 #if DEBUG
00548   for (i=0; i<num; i++) {
00549     printf("unsorted %i: %i %f\n", i, powz[i], wave_f[i]);
00550   }
00551 #endif
00552 
00553   // Sort index array according to power of z
00554   quicksort(powz, index, 0, num-1);
00555 
00556   // Copy angular moments over into sorted array
00557   int   *sorted_ang    = new int[3*num];
00558   for (i=0; i<num; i++) {
00559     sorted_ang[3*i+(comp+1)%3] = qmdata->get_angular_momentum(atom, ishell, first+index[i], (comp+1)%3);
00560     sorted_ang[3*i+(comp+2)%3] = qmdata->get_angular_momentum(atom, ishell, first+index[i], (comp+2)%3);
00561     sorted_ang[3*i+comp] = powz[index[i]];
00562   }
00563 
00564   // Copy sorted angular moments back into original arrays
00565   for (i=0; i<num; i++) {
00566     qmdata->set_angular_momentum(atom, ishell, first+i, &sorted_ang[3*i]);
00567   }
00568 
00569   float *sorted_wave_f = new float[num];
00570 
00571   // Sort the wavefunction coefficients for each orbital
00572   for (orb=0; orb<num_orbitals; orb++) {
00573     wave_f = wave_coeffs + (num_coeffs*orb) + wave_offset + first;
00574 
00575     // Create sorted array
00576     for (i=0; i<num; i++) {
00577       sorted_wave_f[i] = wave_f[index[i]];
00578     }
00579 
00580     // Copy sorted coeffs back to original array
00581     for (i=0; i<num; i++) {
00582       wave_f[i] = sorted_wave_f[i];
00583     }
00584   }
00585 
00586 #if DEBUG
00587   orb=0;
00588   for (i=0; i<qmdata->get_basis(atom, ishell)->num_cart_func; i++) {
00589     printf("sorted %i: %i %i %i %f\n", i, 
00590            qmdata->get_angular_momentum(atom, ishell, i, 0),
00591            qmdata->get_angular_momentum(atom, ishell, i, 1),
00592            qmdata->get_angular_momentum(atom, ishell, i, 2),
00593            wave_coeffs[(num_coeffs*orb) + wave_offset + i]);
00594   }
00595 #endif
00596 
00597   delete [] powz;
00598   delete [] sorted_wave_f;
00599   delete [] sorted_ang;
00600   delete [] index;
00601 }
00602 
00603 
00604 // **************************************************
00605 // ************** QMTimestep class ******************
00606 // **************************************************
00607 
00609 QMTimestep::QMTimestep(int numatoms) :
00610   num_scfiter(0),
00611   num_atoms  (numatoms),
00612   num_wavef  (0),
00613   num_idtags (0),
00614   num_charge_sets(0)
00615 {
00616   wavef_id_map     = NULL;
00617   wavef            = NULL;
00618   scfenergies      = NULL;
00619   gradients        = NULL;
00620   charges          = NULL;
00621   chargetypes      = NULL;
00622 }
00623 
00624 
00626 QMTimestep::QMTimestep(const QMTimestep& ts) 
00627 {
00628   num_scfiter = ts.num_scfiter;
00629   num_atoms   = ts.num_atoms;
00630   num_wavef   = ts.num_wavef;
00631   //wavef_size = ts.wavef_size;
00632   num_idtags = ts.num_idtags;
00633   num_charge_sets = ts.num_charge_sets;
00634 
00635   wavef            = NULL;
00636   wavef_id_map     = NULL;
00637   scfenergies      = NULL;
00638   gradients        = NULL;
00639   charges          = NULL;
00640   chargetypes      = NULL;
00641 
00642   if (ts.wavef) {
00643     int i;
00644     wavef = new Wavefunction[num_wavef];
00645     for (i=0; i<num_wavef; i++) {
00646       wavef[i] = ts.wavef[i];
00647     }
00648   }
00649 
00650   if (ts.wavef_id_map) {
00651     wavef_id_map = (int *)calloc(num_idtags, sizeof(int));
00652     memcpy(wavef_id_map, ts.wavef_id_map, num_idtags*sizeof(int));
00653   }
00654 
00655   if (ts.scfenergies) {
00656     scfenergies = new double[num_scfiter];
00657     memcpy(scfenergies, ts.scfenergies, num_scfiter*sizeof(double));
00658   }
00659 
00660   if (ts.gradients) {
00661     gradients = new float[3*num_atoms];
00662     memcpy(gradients, ts.gradients, 3*num_atoms*sizeof(float));
00663   }
00664 
00665   if (ts.charges) {
00666     charges = new double[num_atoms*num_charge_sets];
00667     memcpy(charges, ts.charges, num_atoms*num_charge_sets*sizeof(double));
00668   }
00669 
00670   if (ts.chargetypes) {
00671     chargetypes = new int[num_charge_sets];
00672     memcpy(chargetypes, ts.chargetypes, num_charge_sets*sizeof(int));
00673   }
00674 }
00675 
00676 
00678 QMTimestep::~QMTimestep() 
00679 {
00680   free(wavef_id_map);
00681   delete [] wavef;
00682   delete [] gradients;
00683   delete [] scfenergies;
00684   delete [] charges;
00685   delete [] chargetypes;
00686 }
00687 
00689 Wavefunction* QMTimestep::get_wavefunction(int iwave)
00690 {
00691   if (iwave<0 || iwave>=num_wavef) return NULL;
00692   return &wavef[iwave];
00693 }
00694 
00696 const float* QMTimestep::get_wavecoeffs(int iwave)
00697 {
00698   if (iwave<0 || iwave>=num_wavef) return NULL;
00699   return wavef[iwave].get_coeffs();
00700 }
00701 
00703 const float* QMTimestep::get_orbitalenergy(int iwave)
00704 {
00705   if (iwave<0 || iwave>=num_wavef) return NULL;
00706   return wavef[iwave].get_orbenergies();
00707 }
00708  
00710 const float* QMTimestep::get_occupancies(int iwave)
00711 {
00712   if (iwave<0 || iwave>=num_wavef) return NULL;
00713   return wavef[iwave].occupancies;
00714 }
00715 
00717 const int* QMTimestep::get_orbitalids(int iwave)
00718 {
00719   if (iwave<0 || iwave>=num_wavef) return NULL;
00720   return wavef[iwave].orb_ids;
00721 }
00722 
00724 const double* QMTimestep::get_charge_set(int iset)
00725 {
00726   if (iset<0 || iset>=num_charge_sets) return NULL;
00727   return &charges[iset*num_atoms];
00728 }
00729 
00731 int QMTimestep::get_charge_type(int iset)
00732 {
00733   if (iset<0 || iset>=num_charge_sets) return -1;
00734   return chargetypes[iset];
00735 }
00736 
00738 const char* QMTimestep::get_charge_type_str(int iset)
00739 {
00740   if (iset<0 || iset>=num_charge_sets) return "";
00741 
00742   switch (chargetypes[iset]) {
00743   case QMCHARGE_MULLIKEN:
00744     return "Mulliken";
00745     break;
00746   case QMCHARGE_LOWDIN:
00747     return "Loewdin";
00748     break;
00749   case QMCHARGE_ESP:
00750     return "ESP";
00751     break;
00752   case QMCHARGE_NPA:
00753     return "NPA";
00754     break;
00755   default:
00756     return "Unknown";
00757   }
00758 }
00759 
00761 int QMTimestep::get_num_coeffs(int iwave) {
00762   if (iwave<0 || iwave>=num_wavef) return -1;
00763   return wavef[iwave].num_coeffs;
00764 }
00765 
00767 int QMTimestep::get_num_orbitals(int iwave) {
00768   if (iwave<0 || iwave>=num_wavef) return -1;
00769   return wavef[iwave].num_orbitals;
00770 }
00771 
00773 int QMTimestep::get_waveid(int iwave) {
00774   if (iwave<0 || iwave>=num_wavef) return -1;
00775   return wavef[iwave].idtag;
00776 }
00777 
00779 int QMTimestep::get_spin(int iwave) {
00780   if (iwave<0 || iwave>=num_wavef) return -1;
00781   return wavef[iwave].spin;
00782 }
00783 
00785 int QMTimestep::get_excitation(int iwave) {
00786   if (iwave<0 || iwave>=num_wavef) return -1;
00787   return wavef[iwave].excitation;
00788 }
00789 
00791 int QMTimestep::get_multiplicity(int iwave) {
00792   if (iwave<0 || iwave>=num_wavef) return -1;
00793   return wavef[iwave].multiplicity;
00794 }
00795 
00796 
00799 double QMTimestep::get_wave_energy(int iwave) {
00800   if (iwave<0 || iwave>=num_wavef) return -1;
00801   return wavef[iwave].energy;
00802 }
00803 
00804 // Translate the wavefunction ID tag into the index the
00805 // wavefunction has in this timestep
00806 int QMTimestep::get_wavef_index(int idtag) {
00807   if (idtag<0 || idtag>=num_idtags) return -1;
00808   return wavef_id_map[idtag];
00809 }
00810 
00812 int QMTimestep::add_wavefunction(QMData *qmdata,
00813                                  int numcoeffs,
00814                                  int numorbitals, 
00815                                  const float *coeffs,
00816                                  const float *orbenergies,
00817                                  float *occupancies,
00818                                  const int *orbids,
00819                                  double energy,
00820                                  int type,
00821                                  int spin,
00822                                  int excitation,
00823                                  int multiplicity,
00824                                  const char *info,
00825                                  wavef_signa_t *(&signa_ts),
00826                                  int &num_signa_ts) 
00827 {
00828   if (!numcoeffs || !numorbitals) return 0;
00829   int iwave = num_wavef;
00830 
00831   Wavefunction *newwavef = new Wavefunction[num_wavef+1];
00832   memset(static_cast<void *>(newwavef), 0, (num_wavef+1)*sizeof(Wavefunction));
00833   int i;
00834   for (i=0; i<num_wavef; i++) {
00835     newwavef[i].movefrom(wavef[i]);
00836   }
00837   delete [] wavef;
00838   wavef = newwavef;
00839   num_wavef++;
00840 
00841   wavef[iwave].energy       = energy;
00842   wavef[iwave].type         = type;
00843   wavef[iwave].spin         = spin;
00844   wavef[iwave].excitation   = excitation;
00845   wavef[iwave].multiplicity = multiplicity;
00846   if (info) 
00847     strncpy(wavef[iwave].info, info, QMDATA_BUFSIZ-1);
00848   else 
00849     wavef[iwave].info[0] = '\0';
00850   
00851   wavef[iwave].set_coeffs(coeffs, numorbitals, numcoeffs);
00852   wavef[iwave].set_orbenergies(orbenergies, numorbitals);
00853   wavef[iwave].set_orbids(orbids, numorbitals);
00854 
00855   if (occupancies) {
00856     wavef[iwave].set_occupancies(occupancies, numorbitals);
00857   } else {
00858     // Assign the MO occupancies depending on RHF, ROHF, UHF
00859     vmd_set_default_occ(wavef[iwave].occupancies,
00860                         qmdata->scftype,
00861                         qmdata->num_electrons,
00862                         numorbitals,
00863                         wavef[iwave].multiplicity);
00864   }
00865 
00866   // Sort the wavefunction coefficients of each shell in
00867   // each orbital according to the angular momenta
00868   wavef[iwave].sort_wave_coefficients(qmdata);
00869 
00870   // Determine a unique ID for each wavefuntion
00871   // based on type, spin, excitation and info string.
00872   int idtag = qmdata->assign_wavef_id(
00873                                        wavef[i].type,
00874                            wavef[i].spin,
00875                            wavef[i].excitation,
00876                            wavef[i].info,
00877                            signa_ts, num_signa_ts);
00878 
00879   // Set the wavefunction ID
00880   set_wavef_idtag(iwave, idtag);
00881 
00882   // Update the list of available orbitals
00883   qmdata->update_avail_orbs(iwave, numorbitals,
00884                             get_orbitalids(i),
00885                             get_occupancies(i));     
00886 
00887   return 1;
00888 }
00889 
00890 // Set timestep independent IDtag for a wavefunction
00891 void QMTimestep::set_wavef_idtag(int iwave, int idtag) {
00892   if (iwave<0 || iwave>=num_wavef) return;
00893   wavef[iwave].idtag = idtag;
00894   if (idtag>=num_idtags) {
00895     if (!wavef_id_map) {
00896       wavef_id_map = (int *)calloc(1, sizeof(int));
00897     } else {
00898       wavef_id_map = (int *)realloc(wavef_id_map,
00899                                     (num_idtags+1)*sizeof(int));
00900     }
00901     num_idtags++;
00902   }
00903   wavef_id_map[idtag] = iwave;
00904 }
00905 
00906 // Initialize SCF energies
00907 void QMTimestep::set_scfenergies(const double *energies, int numscfiter) 
00908 {
00909   if (!energies || !numscfiter) return;
00910 
00911   num_scfiter = numscfiter;
00912   scfenergies = new double[numscfiter];
00913   memcpy(scfenergies, energies, numscfiter*sizeof(double));
00914 }
00915 
00916 // Initialize energy gradient for each atom
00917 void QMTimestep::set_gradients(const float *grad, int natoms) {
00918   if (!grad || !natoms || natoms!=num_atoms) return;
00919 
00920   gradients = new float[3*num_atoms];
00921   memcpy(gradients, grad, 3*num_atoms*sizeof(float));
00922 }
00923 
00924 // Initialize the sets of atom charges.
00925 void QMTimestep::set_charges(const double *q, const int *qtype,
00926                              int natoms, int numqsets) {
00927   if (!q || !natoms || natoms!=num_atoms || !numqsets || !qtype)
00928     return;
00929   num_charge_sets = numqsets;
00930   charges = new double[num_atoms*num_charge_sets];
00931   memcpy(charges, q, num_atoms*num_charge_sets*sizeof(double));
00932 
00933   chargetypes = new int[num_charge_sets];
00934   memcpy(chargetypes, qtype, num_charge_sets*sizeof(int));
00935 }
00936 
00937 
00940 int QMTimestep::get_homo(int iwave) {
00941   if (iwave<0 || iwave>=num_wavef) return -1;
00942 
00943   return wavef[iwave].get_homo();
00944 }
00945 
00946 
00951 int QMTimestep::get_lumo(int iwave) {
00952   if (iwave<0 || iwave>=num_wavef) return -1;
00953   if (!wavef[iwave].get_occupancies()) return -1;
00954   return get_homo(iwave)+1;
00955 }
00956 
00958 void QMTimestep::get_orbital_occ_energy(int iwave, int orb, float &occ, float &energy) {
00959   if (wavef[iwave].get_occupancies() && orb>=0 && orb<wavef[iwave].get_num_orbitals())
00960     occ = wavef[iwave].get_occupancies()[orb];
00961   else
00962     occ = -1.f;      // undefined, return a sentinel value for now
00963 
00964   if (wavef[iwave].get_orbenergies())
00965     energy = wavef[iwave].get_orbenergies()[orb];
00966   else
00967     energy = -666.f; // undefined, return a sentinel value for now
00968 }
00969 
00970 
00975 int QMTimestep::get_orbital_id_from_index(int iwave, int index) {
00976   if (iwave<0 || iwave>=num_wavef) return -1;
00977   if (index<0 || index>=wavef[iwave].num_orbitals) return -1;
00978   return wavef[iwave].orb_ids[index];
00979 }
00980 
00981 
00986 int QMTimestep::get_orbital_index_from_id(int iwave, int id) {
00987   if (iwave<0 || iwave>=num_wavef) return -1;
00988   if (id<1 || id>wavef[iwave].num_orbitals) return -1;
00989   return wavef[iwave].orb_id2index[id];
00990 }
00991 
00992 
00993 // Generate mapping that sorts the orbitals by similarity
00994 // throughout the trajectory (rather than by energy).
00995 void Wavefunction::sort_orbitals(Wavefunction *previous_wavef) {
00996   if (previous_wavef->num_orbitals!=num_orbitals) {
00997     msgInfo << "sort_orbitals(): Number of orbitals differs. Ignoring."
00998             << sendmsg;
00999     return;
01000   }
01001   int *orb_used = new int[num_orbitals];
01002   memset(orb_used, 0, num_orbitals*sizeof(int));
01003   float *scores = new float[num_orbitals];
01004   float *orb_sort_score = new float[num_orbitals];
01005   orb_sort_map = new int[num_orbitals];
01006   int i, j, k;
01007   for (i=0; i<num_orbitals; i++) {
01008     printf("Orbital %d\n", i);
01009     for (j=0; j<previous_wavef->num_orbitals; j++) {
01010       float dot = 0.f;
01011       // Compute the mean square difference between the two
01012       // orbitals.
01013       for (k=0; k<num_coeffs; k++) {
01014         float d = wave_coeffs[i*num_orbitals+k] - 
01015           previous_wavef->wave_coeffs[j*previous_wavef->num_orbitals+k] ;
01016 
01017         dot += d*d;
01018       }
01019       scores[j] = dot/num_coeffs;
01020       printf("score[%d] = % .3f\n", j, scores[j]);
01021     }
01022 
01023     int bestorbital = 0;
01024     float minscore = 1.f;
01025     for (j=0; j<num_orbitals; j++) {
01026       if (scores[j]<minscore) {
01027         minscore = scores[j];
01028         bestorbital = j;
01029       }
01030     }
01031     if (sqrtf(minscore)>0.2) {
01032       msgWarn << "Less than 80% similarity between orbitals."
01033               << sendmsg;
01034       for (k=0; k<num_coeffs; k++) {
01035         printf("coeff[%d]:  % f % f\n", k,
01036                previous_wavef->wave_coeffs[bestorbital*previous_wavef->num_orbitals+k],
01037                wave_coeffs[i*num_orbitals+k]);
01038       }
01039     }
01040     orb_sort_map[i] = bestorbital;
01041     orb_sort_score[i] = sqrtf(minscore);
01042   }
01043 #if 0
01044   float *sorted_coeffs = new float[num_orbitals*num_coeffs];
01045   for (i=0; i<num_orbitals; i++) {
01046     printf("orb_sort_map[%d] = %d, score = %.2f\n", i, orb_sort_map[i], orb_sort_score[i]);
01047     memcpy(&sorted_coeffs[i*num_coeffs],
01048            &wave_coeffs[orb_sort_map[i]*num_coeffs],
01049            num_coeffs*sizeof(float));
01050   }
01051   delete [] wave_coeffs;
01052   wave_coeffs = sorted_coeffs;
01053 #endif  
01054   delete [] scores;
01055   delete [] orb_sort_score;
01056 }
01057 
01058 
01059 // Generate mapping that sorts the orbitals by similarity
01060 // throughout the trajectory (rather than by energy).
01061 void QMTimestep::sort_orbitals(QMTimestep *prev_qmts) {
01062   if (!prev_qmts) return;
01063 
01064   int waveid, i, j;
01065   for (waveid=0; waveid<num_wavef; waveid++) {
01066     i = get_wavef_index(waveid);
01067     j = prev_qmts->get_wavef_index(waveid);
01068 
01069     if (j<0) {
01070       msgInfo << "No wavefunction in timestep matches waveid " 
01071               << waveid << " from previous timestep." << sendmsg;
01072     } else {
01073       wavef[i].sort_orbitals(&prev_qmts->wavef[j]);
01074     }
01075   }
01076 
01077 }
01078 
01079 
01080 
01081 // The standard quicksort algorithm except for it doesn't
01082 // sort the data itself but rather sorts array of ints *idx
01083 // in the same order as it would sort the integers in array
01084 // *tag. Array *idx can then be used to reorder any array
01085 // according to the string tags.
01086 // Example:
01087 // tag:   BC DD BB AA  -->  AA BB BC DD
01088 // index:  0  1  2  3  -->   3  2  0  1
01089 //
01090 static void quicksort(const int* tag, int *idx, int p, int r) {
01091   int q;
01092   if (p < r) {
01093     q=quickpart(tag, idx, p, r);
01094     quicksort(tag, idx, p, q);
01095     quicksort(tag, idx, q+1, r);
01096   }
01097 }
01098 
01099 
01100 // Partitioning for quicksort.
01101 static int quickpart(const int *tag, int *idx, int p, int r) {
01102   int i, j;
01103   int tmp;
01104   int x = tag[idx[p]];
01105   i = p-1;
01106   j = r+1;
01107 
01108   while (1) {
01109     // Find highest element smaller than idx[p]
01110     do j--; while (tag[idx[j]] > x);
01111 
01112     // Find lowest element larger than idx[p]
01113     do i++; while (tag[idx[i]] < x);
01114 
01115     if (i < j) {
01116       tmp    = idx[i];
01117       idx[i] = idx[j];
01118       idx[j] = tmp;
01119     } else {
01120       return j;
01121     }
01122   }
01123 }
01124 
01125 
01126 #if 0
01127 static void quicksort(const char **tag, int *A, int p, int r);
01128 static int  quickpart(const char **tag, int *A, int p, int r);
01129 
01130 // Sorts the indexlist such that the indexes refer to elements in
01131 // tag in an increasing order.
01132 void QMTimestep::sort_shell(QMData *qmdata, int atom, int ishell) {
01133   int i;
01134   const shell_t *shell = qmdata->get_basis(atom, ishell);
01135 
01136   // Initialize the index array;
01137   int *index = new int[shell->num_cart_func];
01138   for (i=0; i<shell->num_cart_func; i++) index[i] = i;
01139 
01140   // Initialize array of sortable tag strings
01141   const char **tag = new const char*[shell->num_cart_func];
01142   for (i=0; i<shell->num_cart_func; i++) {
01143     tag[i] = qmdata->get_angular_momentum(atom, ishell, i);
01144   }
01145 
01146   float *wave_f = wave_function + shell->wave_offset;
01147 #if DEBUG
01148   for (i=0; i<shell->num_cart_func; i++) {
01149     printf("unsorted %i: %s %f\n", i, tag[i], wave_f[i]);
01150   }
01151 #endif
01152 
01153   // Sort index array according to tags.
01154   quicksort(tag, index, 0, shell->num_cart_func-1);
01155 
01156 #if DEBUG
01157   for (i=0; i<shell->num_cart_func; i++) {
01158     printf("sorted %i: %s %f\n", i, tag[index[i]], wave_f[i]);
01159   }
01160 #endif
01161 
01162   // Copy data over into sorted arrays
01163   float *sorted_wave_f = new float[shell->num_cart_func];
01164   char **sorted_tag = new char*[shell->num_cart_func];
01165   for (i=0; i<shell->num_cart_func; i++) {
01166     sorted_wave_f[i] = wave_f[index[i]];
01167     sorted_tag[i] = new char[strlen(tag[index[i]])+1];
01168     strcpy(sorted_tag[i], tag[index[i]]);
01169   }
01170 
01171   // Copy sorted data back into original arrays
01172   for (i=0; i<shell->num_cart_func; i++) {
01173     wave_f[i] = sorted_wave_f[i];
01174     qmdata->set_angular_momentum_str(atom, ishell, i, sorted_tag[i]);
01175   }
01176 
01177   for (i=0; i<shell->num_cart_func; i++) {
01178     delete [] tag[i];
01179     delete [] sorted_tag[i];
01180   }
01181   delete [] sorted_wave_f;
01182   delete [] sorted_tag;
01183   delete [] index;
01184 }
01185 
01186 
01187 // The standard quicksort algorithm except for it doesn't
01188 // sort the data itself but rather sorts array of ints *A
01189 // in the same order as it would sort the strings in array
01190 // **tag. Array *A can then be used to reorder any array
01191 // according to the string tags.
01192 // Example:
01193 // tag:   BC DD BB AA  -->  AA BB BC DD
01194 // index:  0  1  2  3  -->   3  2  0  1
01195 //
01196 static void quicksort(const char **tag, int *idx, int p, int r) {
01197   int q;
01198   if (p < r) {
01199     q=quickpart(tag, idx, p, r);
01200     quicksort(tag, idx, p, q);
01201     quicksort(tag, idx, q+1, r);
01202   }
01203 }
01204 
01205 
01206 // Partitioning for quicksort.
01207 static int quickpart(const char **tag, int *idx,
01208                              int p, int r) {
01209   int i, j;
01210   int tmp;
01211   const char *x = tag[idx[p]];
01212   i = p-1;
01213   j = r+1;
01214 
01215   while (1) {
01216     // Find highest element smaller than idx[p]
01217     do j--; while (strcmp(tag[idx[j]], x) > 0);
01218 
01219     // Find lowest element larger than idx[p]
01220     do i++; while (strcmp(tag[idx[i]], x) < 0);
01221 
01222     if (i < j) {
01223       tmp    = idx[i];
01224       idx[i] = idx[j];
01225       idx[j] = tmp;
01226 
01227     } else {
01228       return j;
01229     }
01230   }
01231 }
01232 
01233 #endif
01234 
01235 
01236 // Assign default occupancies depending on calculation method,
01237 // number of electrons and multiplicity.
01238 // Memory for array *occupancies will be allocated.
01239 void vmd_set_default_occ(float *(&occupancies), int scftype, int numelec, 
01240                          int numorbitals, int multiplicity) 
01241 {
01242   if (!numorbitals) return;
01243   if (occupancies) {
01244     delete [] occupancies;
01245   }
01246   occupancies = new float[numorbitals];
01247 
01248   int i;
01249   for (i=0; i<numorbitals; i++) {
01250     switch(scftype) {
01251       case SCFTYPE_RHF:
01252         if (i<numelec/2) occupancies[i] = 2;
01253         else             occupancies[i] = 0;
01254         break;
01255         
01256       case SCFTYPE_UHF:
01257         if (i<(numelec/2+multiplicity/2)) {
01258           occupancies[i] = 1.f;
01259         } else if ((i >= numorbitals/2) && 
01260                    (i < (numorbitals/2 + numelec/2+multiplicity/2))) {
01261           occupancies[i] = 1.f;
01262         } else {
01263           occupancies[i] = 0.f;
01264         }
01265         break;
01266 
01267       case SCFTYPE_ROHF:
01268         if (i<(numelec/2-multiplicity/2))      occupancies[i] = 2.f;
01269         else if (i<(numelec/2+multiplicity/2)) occupancies[i] = 1.f;
01270         else                                   occupancies[i] = 0.f;
01271         break;
01272 
01273       default:
01274         occupancies[i] = -1.f;
01275         break;
01276     }
01277   }
01278 }

Generated on Thu Mar 28 02:44:06 2024 for VMD (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002