QMTimestep.C

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 }