ael.c

Go to the documentation of this file.

#include <stdio.h>
#include <libqt/qt.h>
#include <psifiles.h>
#define EXTERN
#include "globals.h"
#include <math.h>

void ael(struct RHO_Params *rho_params)
{
  int dim,i,j,k;
  double **rho_g, *evals, **evects, **tmat, **rho_x, ael, **rho_diff, trace;

  dim = moinfo.nmo - moinfo.nfzv;
  rho_g = block_matrix(dim,dim);
  rho_x = block_matrix(dim,dim);
  rho_diff = block_matrix(dim,dim);
  evals = init_array(dim);
  evects = block_matrix(dim,dim);
  tmat = block_matrix(dim,dim);

  /* read and diagonalize the ground-state rho */
  psio_open(PSIF_MO_OPDM, PSIO_OPEN_OLD);
  psio_read_entry(PSIF_MO_OPDM, rho_params[0].opdm_lbl, (char *) &(rho_g[0][0]), sizeof(double)*dim*dim);
  psio_close(PSIF_MO_OPDM, 1);

  sq_rsp(dim, dim, rho_g, evals, 3, evects, 1.0E-14);
  C_DGEMM('t', 'n', dim, dim, dim, 1.0, &(evects[0][0]), dim, &(rho_g[0][0]), dim, 0.0, &(tmat[0][0]), dim);
  C_DGEMM('n', 'n', dim, dim, dim, 1.0,   &(tmat[0][0]), dim, &(evects[0][0]), dim, 0.0, &(rho_g[0][0]), dim);

  for (i=1; i<params.nstates; ++i) {
    /* read in the excited state density */
    psio_open(PSIF_MO_OPDM, PSIO_OPEN_OLD);
    psio_read_entry(PSIF_MO_OPDM, rho_params[i].opdm_lbl, (char *) &(rho_x[0][0]),
                     sizeof(double)*dim*dim);
    psio_close(PSIF_MO_OPDM, 1);

    /* transform the excited-state density */
    C_DGEMM('t', 'n', dim, dim, dim, 1.0, &(evects[0][0]), dim, &(rho_x[0][0]), dim, 0.0, &(tmat[0][0]), dim);
    C_DGEMM('n', 'n', dim, dim, dim, 1.0,   &(tmat[0][0]), dim, &(evects[0][0]), dim, 0.0, &(rho_x[0][0]), dim);

    /* compute the ith AEL */
    ael = 0.0;
    for (j=0; j<dim; ++j) {
      ael += 0.5 * fabs( rho_x[j][k] - rho_g[j][k] );
    }
    fprintf(outfile, "\tAEL %d: %10.7lf\n", i, ael);
  }

  free(evals);
  free_block(evects);
  free_block(tmat);
  free_block(rho_g);
  free_block(rho_x);
  free_block(rho_diff);
  return;
}

---