libqt: The Quantum-Trio Miscellaneous Library

Files
file	3d_array.cc
	Routines for 3d arrays.
file	blas_intfc.cc
	The PSI3 BLAS interface routines.
file	cc_excited.cc
	Determine if a wavefunction is a CC-excited wavefunction type.
file	cc_wfn.cc
	Check if wavefunction is coupled-cluster type.
file	ci_wfn.cc
	Check if wavefunction is CI-type.
file	david.cc
	In-core Davidson-Liu diagonalization of symm matrices.
file	dirprd_block.cc
	Take a direct product of two matrices.
file	dot_block.cc
	Take dot product of two block matrices.
file	eri.cc
	Inefficient electron repulsion integral evaluation By Edward Valeev.
file	fill_sym_matrix.cc
	Fill a symmetric matrix from a lower triangle.
file	filter.cc
	Filter out unneeded frozen core/virt integrals.
file	get_frzpi.cc
	Get frozen core/virtuals per irrep.
file	invert.cc
	Invert a small matrix.
file	lapack_intfc.cc
	Interface to LAPACK routines.
file	mat_in.cc
	read in a matrix from an input stream (deprecated)
file	mat_print.cc
	Print a matrix to a file in a formatted style.
file	newmm_rking.cc
	Matrix multiplication routine (deprecated).
file	normalize.cc
	Normalize a set of vectors.
file	pople.cc
	Pople's method for solving linear equations.
file	probabil.cc
	Contains some probability functions.
file	qt.h
	Header file for the Quantum Trio Library.
file	ras_set.cc
	Obtain orbital space and reordering for CI/MCSCF wavefunctions.
file	reorder_qt.cc
	Obtain the QT orbital reordering array between Pitzer and correlated order.
file	rootfind.cc
	Simple root finding methods.
file	schmidt.cc
	Gram-Schmidt orthogonalize a set of vectors.
file	lib/libqt/schmidt_add.cc
	Gram-Schmidt orthogonalize vector and add to set.
file	slaterdset.cc
	Utility functions for importing/exporting sets of Slater determinants from the CI codes.
file	slaterdset.h
	Header file for SlaterDetSets.
file	solve_pep.cc
	Solve a 2x2 pseudo-eigenvalue problem.
file	lib/libqt/sort.cc
	Sort eigenvalues and eigenvectors into ascending order.
file	strncpy.cc
	Override strncpy to ensure strings are terminated.
file	timer.cc
	Obtain user and system timings for blocks of code.
Functions
double ***	init_3d_array (int p, int q, int r)
void	free_3d_array (double ***A, int p, int q)
void	C_DAXPY (int length, double a, double *x, int inc_x, double *y, int inc_y)
void	C_DCOPY (int length, double *x, int inc_x, double *y, int inc_y)
void	C_DSCAL (int length, double alpha, double *vec, int inc)
void	C_DROT (int length, double *x, int inc_x, double *y, int inc_y, double costheta, double sintheta)
void	C_DGEMM (char transa, char transb, int m, int n, int k, double alpha, double *A, int nca, double *B, int ncb, double beta, double *C, int ncc)
void	C_DGEMV (char transa, int m, int n, double alpha, double *A, int nca, double *X, int inc_x, double beta, double *Y, int inc_y)
void	C_DSPMV (char uplo, int n, double alpha, double *A, double *X, int inc_x, double beta, double *Y, int inc_y)
double	C_DDOT (int n, double *x, int inc_x, double *y, int inc_y)
int	cc_excited (char *wfn)
int	cc_wfn (char *wfn)
int	ci_wfn (char *wfn)
int	david (double **A, int N, int M, double *eps, double **v, double cutoff, int print)
void	dirprd_block (double **A, double **B, int rows, int cols)
double	dot_block (double **A, double **B, int rows, int cols, double alpha)
double	eri (unsigned int l1, unsigned int m1, unsigned int n1, double alpha1, double A[3], unsigned int l2, unsigned int m2, unsigned int n2, double alpha2, double B[3], unsigned int l3, unsigned int m3, unsigned int n3, double alpha3, double C[3], unsigned int l4, unsigned int m4, unsigned int n4, double alpha4, double D[3], int norm_flag)
double	norm_const (unsigned int l1, unsigned int m1, unsigned int n1, double alpha1, double A[3])
void	fill_sym_matrix (double **A, int size)
void	filter (double *input, double *output, int *ioff, int norbs, int nfzc, int nfzv)
int *	get_frzcpi ()
int *	get_frzvpi ()
double	invert_matrix (double **a, double **y, int N, FILE *outfile)
int	C_DGEEV (int n, double **a, int lda, double *wr, double *wi, double **vl, int ldvl, double **vr, int ldvr, double *work, int lwork, int info)
int	C_DGESV (int n, int nrhs, double *a, int lda, int *ipiv, double *b, int ldb)
int	C_DGETRF (int nrow, int ncol, double *a, int lda, int *ipiv)
int	C_DGETRI (int n, double *a, int lda, int *ipiv, double *work, int lwork)
int	C_DGESVD (char jobu, char jobvt, int m, int n, double *A, int lda, double *s, double *u, int ldu, double *vt, int ldvt, double *work, int lwork)
int	C_DSYEV (char jobz, char uplo, int n, double *A, int lda, double *w, double *work, int lwork)
int	mat_in (FILE *fp, double **array, int width, int max_length, int *stat)
int	mat_print (double **matrix, int rows, int cols, FILE *outfile)
void	normalize (double **A, int rows, int cols)
int	pople (double **A, double *x, int dimen, int num_vecs, double tolerance, FILE *outfile, int print_lvl)
double	combinations (int n, int k)
double	factorial (int n)
int	ras_set (int nirreps, int nbfso, int freeze_core, int *orbspi, int *docc, int *socc, int *frdocc, int *fruocc, int **ras_opi, int *order, int ras_type)
int	ras_set2 (int nirreps, int nmo, int delete_fzdocc, int delete_restrdocc, int *orbspi, int *docc, int *socc, int *frdocc, int *fruocc, int *restrdocc, int *restruocc, int **ras_opi, int *order, int ras_type, int hoffmann)
void	reorder_qt (int *docc_in, int *socc_in, int *frozen_docc_in, int *frozen_uocc_in, int *order, int *orbs_per_irrep, int nirreps)
void	reorder_qt_uhf (int *docc, int *socc, int *frozen_docc, int *frozen_uocc, int *order_alpha, int *order_beta, int *orbspi, int nirreps)
double	bisect (double(*function)(double), double low, double high, double tolerance, int maxiter, int printflag)
double	newton (double(*F)(double), double(*dF)(double), double x, double tolerance, int maxiter, int printflag)
double	secant (double(*F)(double), double x0, double x1, double tolerance, int maxiter, int printflag)
void	schmidt (double **A, int rows, int cols, FILE *outfile)
int	schmidt_add (double **A, int rows, int cols, double *v)
void	stringset_init (StringSet *sset, int size, int nelec, int nfzc, short int *frozen_occ)
void	stringset_delete (StringSet *sset)
void	stringset_add (StringSet *sset, int index, unsigned char *Occ)
void	stringset_reindex (StringSet *sset, short int *mo_map)
void	stringset_write (ULI unit, char *prefix, StringSet *sset)
void	stringset_read (ULI unit, char *prefix, StringSet **stringset)
void	slaterdetset_init (SlaterDetSet *sdset, int size, StringSet *alphastrings, StringSet *betastrings)
void	slaterdetset_delete (SlaterDetSet *sdset)
void	slaterdetset_delete_full (SlaterDetSet *sdset)
void	slaterdetset_add (SlaterDetSet *sdset, int index, int alphastring, int betastring)
void	slaterdetset_write (ULI unit, char *prefix, SlaterDetSet *sdset)
void	slaterdetset_read (ULI unit, char *prefix, SlaterDetSet **slaterdetset)
void	slaterdetvector_init (SlaterDetVector *sdvector, SlaterDetSet *sdset)
void	slaterdetvector_delete (SlaterDetVector *sdvector)
void	slaterdetvector_delete_full (SlaterDetVector *sdvector)
void	slaterdetvector_add (SlaterDetVector *sdvector, int index, double coeff)
void	slaterdetvector_set (SlaterDetVector *sdvector, double *coeffs)
void	slaterdetvector_write (ULI unit, char *prefix, SlaterDetVector *vector)
void	slaterdetset_write_vect (ULI unit, char *prefix, double *coeffs, int size, int vectnum)
void	slaterdetvector_read (ULI unit, char *prefix, SlaterDetVector **sdvector)
void	slaterdetset_read_vect (ULI unit, char *prefix, double *coeffs, int size, int vectnum)
void	solve_2x2_pep (double **H, double S, double *evals, double **evecs)
void	sort_vector (double *A, int n)
void	timer_init (void)
void	timer_done (void)
struct timer *	timer_scan (char *key)
void	timer_on (char *key)
void	timer_off (char *key)

---

Detailed Description

---

Function Documentation

double bisect	(	double(*)(double)	function,
		double	low,
		double	high,
		double	tolerance,
		int	maxiter,
		int	printflag
	)

bisect(): Finds the root of a function between two points to within a given tolerance. Iterations are limited to a given maximum, and a print flag specifies whether this function should print the results of each iteration to stdout. Note that the values of the function at the two endpoints of the interval must have _different_ signs for the bisection method to work! This routine checks for this initially.

Parameters:

	function	= pointer to function we want to examine (must return double)
	low	= lower bound of interval to search for root
	high	= upper bound of interval
	tolerance	= how small is the maximum allowable error
	maxiter	= maximum number of iterations
	printflag	= whether or not to print results for each iteration (1 or 0)

Returns: the value of the root

Definition at line 56 of file rootfind.cc.

{
   int iter = 1 ;
   double tmpval ;
   double center ;      /* half way between low and high, i.e. current guess */
   double fl, fc, fh ;  /* values of function at the three points */
   double errorval ; 

   /* check the input parameters */
   if (low > high) {  /* make sure low is the smaller of the 2 limits */
      tmpval = low ;
      low = high ; 
      high = tmpval ;
      }
 
   fl = function(low) ;
   fh = function(high) ;
   if (fl * fh > 0.0) {
      printf("(bisect): Product of function values is greater than 0\n") ;
      return(0.0) ;
      }

   if (printflag) {
      printf("Bisection method root finder\n") ;
      printf("%4s %12s  %12s\n", "iter", "x", "max error") ;
      }

   /* now iterate until we hit maxiter or until we get within the tolerance */
   while (iter < maxiter) {  
      center = (low + high) / 2.0 ;
      fc = function(center) ;
      errorval = (high - low) / 2.0 ;
      if (printflag) 
         printf("%4d %12.8f  %12.8f\n", iter, center, errorval) ;
      if (errorval < tolerance) return(center) ;
      if (fl * fc > 0.0) {
         low = center ;
         fl = fc ;
         }
      else {
         high = center ;
         fh = fc ;
         }
      iter++ ;
      }

   printf("(bisect): Failed to converge within %d iterations\n", iter) ;
   return(center) ;
} 

void C_DAXPY	(	int	length,
		double	a,
		double *	x,
		int	inc_x,
		double *	y,
		int	inc_y
	)

C_DAXPY(): This function performs y = a * x + y.

Steps every inc_x in x and every inc_y in y (normally both 1).

Parameters:

	length	= length of arrays
	a	= scalar a to multiply vector x
	x	= vector x
	inc_x	= how many places to skip to get to next element in x
	y	= vector y
	inc_y	= how many places to skip to get to next element in y

Returns:none

Definition at line 93 of file blas_intfc.cc.

{
  F_DAXPY(&length, &a, x, &inc_x, y, &inc_y);
}

void C_DCOPY	(	int	length,
		double *	x,
		int	inc_x,
		double *	y,
		int	inc_y
	)

C_DCOPY(): This function copies x into y.

Steps every inc_x in x and every inc_y in y (normally both 1).

Parameters:

	length	= length of array
	x	= vector x
	inc_x	= how many places to skip to get to next element in x
	y	= vector y
	inc_y	= how many places to skip to get to next element in y

Returns: none

Definition at line 115 of file blas_intfc.cc.

{
  F_DCOPY(&length, x, &inc_x, y, &inc_y);
}

double C_DDOT	(	int	n,
		double *	x,
		int	inc_x,
		double *	y,
		int	inc_y
	)

C_DDOT(): This function returns the dot product of two vectors, x and y.

Parameters:

	length	= Number of elements in x and y.
	x	= A pointer to the beginning of the data in x. Must be of at least length (1+(N-1)*abs(inc_x).
	inc_x	= how many places to skip to get to next element in x
	y	= A pointer to the beginning of the data in y.
	inc_y	= how many places to skip to get to next element in y

Returns: the dot product

Interface written by ST Brown. July 2000

Definition at line 365 of file blas_intfc.cc.

{
   if(n == 0) return 0.0;

   return F_DDOT(&n,x,&inc_x,y,&inc_y);
}

int C_DGEEV	(	int	n,
		double **	a,
		int	lda,
		double *	wr,
		double *	wi,
		double **	vl,
		int	ldvl,
		double **	vr,
		int	ldvr,
		double *	work,
		int	lwork,
		int	info
	)

C_DGEEV()

This function computes the eigenvalues and the left and right right eigenvectors of a real, nonsymmetric matrix A. For symmetric matrices, refer to C_DSYEV().

Parameters:

	n	= The order of the matrix A. n >= 0.
	a	= The n-by-n matrix A. As with all other lapack routines, must be allocated by contiguous memory (i.e., block_matrix()).
	lda	= The leading dimension of matrix A. lda >= max(1,n).
	wr	= array of length n containing the real parts of computed eigenvalues.
	wi	= array of length n containing the imaginary parts of computed eigenvalues.
	vl	= matrix of dimensions ldvl*n. The columns store the left eigenvectors u(j). If the j-th eigenvalues is real, then u(j) = vl(:,j), the j-th column of vl. If the j-th and (j+1)-st eigenvalues form a complex conjugate pair, then u(j) = vl(:,j) + i*vl(:,j+1) and u(j+1) = vl(:,j) - i*vl(:,j+1). Note: this is the Fortran documentation, may need to change cols <-> rows.
	ldvl	= The leading dimension of matrix vl.
	vr	= matrix of dimensions ldvr*n. The columns store the right eigenvectors v(j). If the j-th eigenvalues is real, then v(j) = vr(:,j), the j-th column of vr. If the j-th and (j+1)-st eigenvalues form a complex conjugate pair, then v(j) = vr(:,j) + i*vr(:,j+1) and v(j+1) = vr(:,j) - i*vr(:,j+1). Note: this is the Fortran documentation, may need to change cols <-> rows.
	ldvr	= The leading dimension of matrix vr.
	work	= Array for scratch computations, of dimension lwork. On successful exit, work[0] returns the optimal value of lwork.
	lwork	= The dimension of the array work. lwork >= max(1,3*n), and if eigenvectors are required (default for this wrapper at present) then actually lwork >= 4*n. For good performance, lwork must generally be larger. If lwork = -1, then a workspace query is assumed. The routine only calculates the optimal size of the work array, returns this value ans the first entry of the work array, and no error message related to lword is issued by xerbla.
	info	= On output (returned by C_DGEEV), a status flag. info = 0 for successful exit, if info = -i, the ith argument had an illegal value. If info = i, the QR algorithm failed to compute all the eigenvalues, and no eigenvectors have been computed. Elements i+1:n of wr and wi contain eigenvalues which have converged.

Definition at line 119 of file lapack_intfc.cc.

{
  char jobvl, jobvr;
  jobvl = 'V';
  jobvr = 'V';
  F_DGEEV(&jobvl, &jobvr, &n, &(a[0][0]), &lda, &(wr[0]), &(wi[0]),
       &(vl[0][0]), &ldvl, &(vr[0][0]), &ldvr, &(work[0]), &lwork, &info);

  return info;
}

void C_DGEMM	(	char	transa,
		char	transb,
		int	m,
		int	n,
		int	k,
		double	alpha,
		double *	A,
		int	nca,
		double *	B,
		int	ncb,
		double	beta,
		double *	C,
		int	ncc
	)

C_DGEMM(): This function calculates C(m,n)=alpha*(opT)A(m,k)*(opT)B(k,n)+ beta*C(m,n)

These arguments mimic their Fortran conterparts; parameters have been reversed (nca, ncb, ncc, A, B, C), to make it correct for C.

Parameters:

	transa	= On entry, specifies the form of (op)A used in the matrix multiplication: If transa = 'N' or 'n', (op)A = A. If transa = 'T' or 't', (op)A = transp(A). If transa = 'R' or 'r', (op)A = conjugate(A). If transa = 'C' or 'c', (op)A = conjug_transp(A). On exit, transa is unchanged.
	transb	= On entry, specifies the form of (op)B used in the matrix multiplication: If transb = 'N' or 'n', (op)B = B. If transb = 'T' or 't', (op)B = transp(B). If transb = 'R' or 'r', (op)B = conjugate(B)
	m	= On entry, the number of rows of the matrix (op)A and of the matrix C; m >= 0. On exit, m is unchanged.
	n	= On entry, the number of columns of the matrix (op)B and of the matrix C; n >= 0. On exit, n is unchanged.
	k	= On entry, the number of columns of the matrix (op)A and the number of rows of the matrix (op)B; k >= 0. On exit, k is unchanged.
	alpha	= On entry, specifies the scalar alpha. On exit, alpha is unchanged.
	A	= On entry, a two-dimensional array A with dimensions ka by nca. For (op)A = A or conjugate(A), nca >= k and the leading m by k portion of the array A contains the matrix A. For (op)A = transp(A) or conjug_transp(A), nca >= m and the leading k by m part of the array A contains the matrix A. On exit, a is unchanged.
	nca	= On entry, the second dimension of array A. For (op)A = A or conjugate(A), nca >= MAX(1,k). For (op)A=transp(A) or conjug_transp(A), nca >= MAX(1,m). On exit, nca is unchanged.
	B	= On entry, a two-dimensional array B with dimensions kb by ncb. For (op)B = B or conjugate(B), kb >= k and the leading k by n portion of the array contains the matrix B. For (op)B = transp(B) or conjug_transp(B), ncb >= k and the leading n by k part of the array contains the matrix B. On exit, B is unchanged.
	ncb	= On entry, the second dimension of array B. For (op)B = B or <conjugate(B), ncb >= MAX(1,n). For (op)B = transp(B) or conjug_transp(B), ncb >= MAX(1,k). On exit, ncb is unchanged.
	beta	= On entry, specifies the scalar beta. On exit, beta is unchanged.
	C	= On entry, a two-dimensional array with the dimension at least m by ncc. On exit, the leading m by n part of array C is overwritten by the matrix alpha*(op)A*(op)B + beta*C.
	ncc	= On entry, the second dimension of array C; ncc >=MAX(1,n). On exit, ncc is unchanged.

Returns: none

Definition at line 235 of file blas_intfc.cc.

Referenced by ael(), david(), dpd_contract244(), and dpd_contract424().

{

  /* the only strange thing we need to do is reverse everything
     since the stride runs differently in C vs. Fortran
   */
  
  /* also, do nothing if a dimension is 0 */
  if (m == 0 || n == 0 || k == 0) return;

  F_DGEMM(&transb,&transa,&n,&m,&k,&alpha,B,&ncb,A,&nca,&beta,C,&ncc);

}

void C_DGEMV	(	char	transa,
		int	m,
		int	n,
		double	alpha,
		double *	A,
		int	nca,
		double *	X,
		int	inc_x,
		double	beta,
		double *	Y,
		int	inc_y
	)

C_DGEMV(): This function calculates the matrix-vector product.

Y = alpha * A * X + beta * Y

where X and Y are vectors, A is a matrix, and alpha and beta are constants.

Parameters:

	transa	= Indicates whether the matrix A should be transposed ('t') or left alone ('n').
	m	= The row dimension of A (regardless of transa).
	n	= The column dimension of A (regardless of transa).
	alpha	= The scalar alpha.
	A	= A pointer to the beginning of the data in A.
	nca	= The number of columns *actually* in A. This is useful if one only wishes to multiply the first n columns of A times X even though A contains nca columns.
	X	= A pointer to the beginning of the data in X.
	inc_x	= The desired stride for X. Useful for skipping sections of data to treat only one column of a complete matrix. Usually 1, though.
	beta	= The scalar beta.
	Y	= A pointer to the beginning of the data in Y.
	inc_y	= The desired stride for Y.

Returns: none

Interface written by TD Crawford and EF Valeev. June 1999

Definition at line 286 of file blas_intfc.cc.

{
  if (m == 0 || n == 0) return;

  if(transa == 'n') transa = 't';
  else transa = 'n';

  F_DGEMV(&transa,&n,&m,&alpha,A,&nca,X,&inc_x,&beta,Y,&inc_y);

}

int C_DGESV	(	int	n,
		int	nrhs,
		double *	a,
		int	lda,
		int *	ipiv,
		double *	b,
		int	ldb
	)

C_DGESV()

This function solves a system of linear equations A * X = B, where A is an n x n matrix and X and B are n x nrhs matrices.

Parameters:

	n	= The number of linear equations, i.e., the order of the matrix A. n >= 0.
	nrhs	= The number of right hand sides, i.e., the number of columns of the matrix B. nrhs >= 0.
	A	= On entry, the n-by-n coefficient matrix A. On exit, the factors L and U from the factorization A = P*L*U; the unit diagonal elements of L are not stored.

lda = The leading dimension of the array A. lda >= max(1,n).

Parameters:

	int	ipiv = An integer array of length n. The pivot indices that define the permutation matrix P; row i of the matrix was interchanged with row ipiv(i).
	B	= On entry, the n-by-nrhs matrix of right hand side matrix B. On exit, if info = 0, the n-by-nrhs solution matrix X.
	ldb	= The leading dimension of the array B. ldb >= max(1,n).

Returns: int info. info = 0: successful exit. info < 0: if info = -i, the i-th argument had an illegal value. info > 0: if info = i, U(i,i) is exactly zero. The factorization has been completed, but the factor U is exactly singular, so the solution could not be computed.

Definition at line 169 of file lapack_intfc.cc.

{
  int info;

  F_DGESV(&n, &nrhs, &(a[0]), &lda, &(ipiv[0]), &(b[0]), &ldb, &info);

  return info;
}

int C_DGESVD	(	char	jobu,
		char	jobvt,
		int	m,
		int	n,
		double *	A,
		int	lda,
		double *	s,
		double *	u,
		int	ldu,
		double *	vt,
		int	ldvt,
		double *	work,
		int	lwork
	)

C_DGESVD() This function computes the singular value decomposition (SVD) of a real mxn matrix A, ** optionally computing the left and/or right singular vectors. The SVD is written

A = U * S * transpose(V)

where S is an mxn matrix which is zero except for its min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and V is an N-by-N orthogonal matrix. The diagonal elements of SIGMA are the singular values of A; they are real and non-negative, and are returned in descending order. The first min(m,n) columns of U and V are the left and right singular vectors of A.

Note that the routine returns V^t, not V;

These arguments mimic their Fortran counterparts. See the LAPACK manual for additional information.

Parameters:

	jobu	= 'A' = all m columns of U are returned 'S' = the first min(m,n) columns of U are returned 'O' = the first min(m,n) columns of U are returned in the input matrix A 'N' = no columns of U are returned
	jobvt	= 'A' = all n rows of VT are returned 'S' = the first min(m,n) rows of VT are returned 'O' = the first min(m,n) rows of VT are returned in the input matrix A 'N' = no rows of VT are returned

Obviously jobu and jobvt cannot both be 'O'.

Parameters:

	m	= The row dimension of the matrix A.
	n	= The column dimension of the matrix A.
	A	= On entry, the mxn matrix A with dimensions m by lda. On exit, if jobu='O' the first min(m,n) columns of A are overwritten with the left singular vectors; if jobvt='O' the first min(m,n) rows of A are overwritten with the right singular vectors; otherwise, the contents of A are destroyed.
	lda	= The number of columns allocated for A.
	s	= The singular values of A.
	u	= The right singular vectors of A, returned as column of the matrix u if jobu='A'. If jobu='N' or 'O', u is not referenced.
	ldu	= The number of columns allocated for u.
	vt	= The left singular vectors of A, returned as rows of the matrix VT is jobvt='A'. If jobvt='N' or 'O', vt is not referenced.
	ldvt	= The number of columns allocated for vt.
	work	= Workspace array of length lwork.
	lwork	= The length of the workspace array work, which should be at least as large as max(3*min(m,n)+max(m,n),5*min(m,n)). For good performance, lwork should generally be larger. If lwork=-1, a workspace query is assumed, and the value of work[0] upon return is the optimal size of the work array.

Return value: int info. If info=0, successful exit. If info = -i, the i-th argument had an illegal value. If info > 0, Related to failure in the convergence of the upper bidiagonal matrix B. See the LAPACK manual for additiona information.

Interface written by TDC, July 2001, updated April 2004

Definition at line 339 of file lapack_intfc.cc.

{
  int info;

  F_DGESVD(&jobvt, &jobu, &n, &m, A, &lda, s, vt, &ldvt, u, &ldu, work, 
    &lwork, &info);

  return info;
}

int C_DGETRF	(	int	nrow,
		int	ncol,
		double *	a,
		int	lda,
		int *	ipiv
	)

C_DGETRF(): Compute an LU factorization of a general M-by-N matrix a using partial pivoting with row interchanges.

The factorization has the form A = P * L * U where P is a permutation matrix, L is lower triangular with unit diagonal elements (lower trapezoidal if nrow > ncol), and U is upper triangular (upper trapezoidal if nrow < ncol).

Parameters:

	nrow	= number of rows
	ncol	= number of colums
	a	= matrix to factorize
	lda	= leading dimension of a, lda >= max(1,ncol)
	ipiv	= output integer array of the pivot indices; for 1 <= i <= min(nrow, ncol), col i of the matrix was interchanged with col ipiv(i).

Returns: int info. info = 0: successful exit. info < 0: if info=-i, the ith-argument had an illegal vale. info>0: if info=i, U(i,i) is exactly zero. The factorization has been completed, but the factor U is exactly singular, and division by zero will occur if it is used to solve a system of equations.

Definition at line 205 of file lapack_intfc.cc.

{
  int info;

  F_DGETRF(&ncol, &nrow, &(a[0]), &lda, &(ipiv[0]), &info);

  return info;
}

int C_DGETRI	(	int	n,
		double *	a,
		int	lda,
		int *	ipiv,
		double *	work,
		int	lwork
	)

C_DGETRI(): computes the inverse of a matrix using the LU factorization computed by DGETRF

This method inverts U and then computes inv(A) by solving the system inv(A)*L = inv(U) for inv(A).

Parameters:

	n	= order of the matrix a. n >= 0.
	a	= array of dimension (n,lda). On entry, the factors L and U from the factorization A = P*L*U as computed by DGETRF. On exit, if info=0, the inverse of hte original matrix A.
	lda	= the leading dimension of the array a. lda >= max(1,n).
	ipiv	= The pivot indices from DGETRF. For 1<=i<=n, row i of the matrix was interchanged with row ipiv(i)
	work	= workspace of dimension max(1,lwork). On exit, if info=0, then work(0) returns the optimal lwork.
	lwork	= the dimension of the array work. lwork >= max(1,n). For optimal performance, lwork > n*nb, where nb is the optimal blocksize returned by ilaenv. If lwork=-1, then a workspace query is assumed; the routine only calculates the optimal size of the work array, returns this value as the first entry of the work array, and no error message related to lwork is issued by xerbla.

Returns: int info. info=0: successful exit. info<0: if info=-i, the i-th argument had an illegal value. info>0: if info=i, U(i,i) is exactly zero; the matrix is singular and its inverse could not be computed.

Definition at line 246 of file lapack_intfc.cc.

{
  int info;

  F_DGETRI(&n, &(a[0]), &lda, &(ipiv[0]), &(work[0]), &lwork, &info);

  return info;
}

void C_DROT	(	int	length,
		double *	x,
		int	inc_x,
		double *	y,
		int	inc_y,
		double	costheta,
		double	sintheta
	)

C_DROT(): Calculates a plane Givens rotation for vectors x, y and angle theta. x = x*cos + y*sin, y = -x*sin + y*cos.

Parameters:

	x	= vector x
	y	= vector Y
	length	= length of x,y
	inc_x	= how many places to skip to get to the next element of x
	inc_y	= how many places to skip to get to the next element of y

Returns: none

Definition at line 154 of file blas_intfc.cc.

{

  F_DROT(&length,x,&inc_x,y,&inc_y,&costheta,&sintheta);
}

void C_DSCAL	(	int	length,
		double	alpha,
		double *	vec,
		int	inc
	)

C_DSCAL(): This function scales a vector by a real scalar.

Parameters:

	length	= length of array
	alpha	= scale factor
	vec	= vector to scale
	inc	= how many places to skip to get to next element in vec

Returns: none

Definition at line 134 of file blas_intfc.cc.

Referenced by psi::detcas::calc_gradient(), and dpd_buf4_scm().

{
  F_DSCAL(&length, &alpha, vec, &inc);
}

void C_DSPMV	(	char	uplo,
		int	n,
		double	alpha,
		double *	A,
		double *	X,
		int	inc_x,
		double	beta,
		double *	Y,
		int	inc_y
	)

C_DSPMV(): This function calculates the matrix-vector product

Y = alpha * A * X + beta * Y

where X and Y are vectors, A is a matrix, and alpha and beta are constants.

Parameters:

	uplo	= Indicates whether the matrix A is packed in upper ('U' or 'u') or lower ('L' or 'l') triangular form. We reverse what is passed before sending it on to Fortran because of the different Fortran/C conventions
	n	= The order of the matrix A (number of rows/columns)
	alpha	= The scalar alpha.
	A	= A pointer to the beginning of the data in A.
	X	= A pointer to the beginning of the data in X.
	inc_x	= The desired stride for X. Useful for skipping sections of data to treat only one column of a complete matrix. Usually 1, though.
	beta	= The scalar beta.
	Y	= A pointer to the beginning of the data in Y.
	inc_y	= The desired stride for Y.

Returns: none

Interface written by CD Sherrill July 2003

Definition at line 331 of file blas_intfc.cc.

{
  if (n == 0) return;

  if (uplo != 'U' && uplo != 'u' && uplo != 'L' && uplo != 'l')
    fprintf(stderr, "C_DSPMV: called with unrecognized option for uplo!\n");

  if (uplo == 'U' || uplo == 'u') uplo = 'L';
  else uplo = 'U';

  F_DSPMV(&uplo,&n,&alpha,A,X,&inc_x,&beta,Y,&inc_y);

}

int C_DSYEV	(	char	jobz,
		char	uplo,
		int	n,
		double *	A,
		int	lda,
		double *	w,
		double *	work,
		int	lwork
	)

C_DSYEV(): Computes all eigenvalues and, optionally, eigenvectors of a real symmetric matrix A.

These arguments mimic their Fortran counterparts.

Parameters:

	jobz	= 'N' or 'n' = compute eigenvalues only; 'V' or 'v' = compute both eigenvalues and eigenvectors.
	uplo	= 'U' or 'u' = A contains the upper triangular part of the matrix; 'L' or 'l' = A contains the lower triangular part of the matrix.
	n	= The order of the matrix A.
	A	= On entry, the two-dimensional array with dimensions n by lda. On exit, if jobz = 'V', the columns of the matrix contain the eigenvectors of A, but if jobz = 'N', the contents of the matrix are destroyed.
	lda	= The second dimension of A (i.e., the number of columns allocated for A).
	w	= The computed eigenvalues in ascending order.
	work	= An array of length lwork. On exit, if the return value is 0, work[0] contains the optimal value of lwork.
	lwork	= The length of the array work. A useful value of lwork seems to be 3*N.

Returns: 0 = successful exit <0 = the value of the i-th argument to the function was illegal >0 = the algorithm failed to converge.

Interface written by TDC, 10/2002

Definition at line 392 of file lapack_intfc.cc.

{
  int info;

  F_DSYEV(&jobz, &uplo, &n, A, &lda, w, work, &lwork, &info);

  return info;
}

int cc_excited

(

char *

wfn

)

cc_excited(): This function takes a WFN string and returns 1 if the WFN is an excited-state method and 0 if the WFN is a ground-state method.

Parameters:

*wfn

= wavefunction string

Returns: 1 if an excited state method, else 0

Definition at line 23 of file cc_excited.cc.

Referenced by main().

{

  if ( !strcmp(wfn, "CCSD")   || !strcmp(wfn, "CCSD_T") || !strcmp(wfn, "BCCD") 
    || !strcmp(wfn, "BCCD_T") || !strcmp(wfn, "CC2")    || !strcmp(wfn, "CC3") 
    || !strcmp(wfn, "CCSD_MVD") ) {
    return 0;
  }
  else if ( !strcmp(wfn, "EOM_CCSD") || !strcmp(wfn, "LEOM_CCSD") ||
            !strcmp(wfn, "EOM_CC2")  || !strcmp(wfn, "EOM_CC3") ) {
    return 1;
  }
  else {
    printf("Invalid value of input keyword WFN: %s\n", wfn);
    exit(PSI_RETURN_FAILURE);
  }
}

int cc_wfn

(

char *

wfn

)

cc_wfn(): Checks if the given wavefunction string is a coupled-cluster type and returns 1 if yes and 0 if no.

Note: "coupled-cluster type" means it is handled by PSI like the coupled-cluster codes, not necessarily that it is literally a coupled-cluster wavefunction

Parameters:

*wfn

= wavefunction string

Returns: 1 if the WFN is a CC method, 0 otherwise

Definition at line 28 of file cc_wfn.cc.

{

  if ( !strcmp(wfn, "CCSD")     || !strcmp(wfn, "CCSD_T") || 
       !strcmp(wfn, "BCCD")     || !strcmp(wfn, "BCCD_T") || 
       !strcmp(wfn, "CC2")      || !strcmp(wfn, "CC3")    ||
       !strcmp(wfn, "EOM_CCSD") || !strcmp(wfn, "LEOM_CCSD") ||
       !strcmp(wfn, "EOM_CC2")  || !strcmp(wfn, "EOM_CC3") ||
       !strcmp(wfn, "CIS") ) {
    return 1;
  }
  else {
    return 0;
  }
}

int ci_wfn

(

char *

wfn

)

ci_wfn(): Examine the wavefunction type and return 1 if a CI/MCSCF-type, otherwise 0

Parameters:

*wfn

= wavefunction string

Returns: 1 if a CI/MCSCF-type wavefunction, otherwise 0

Definition at line 24 of file ci_wfn.cc.

{

  if (strcmp(wfn, "CI")==0     || strcmp(wfn, "DETCAS")==0 || 
      strcmp(wfn, "CASSCF")==0 || strcmp(wfn, "RASSCF")==0 ||
      strcmp(wfn, "DETCI")==0 || strcmp(wfn, "MCSCF")==0 ||
      strcmp(wfn, "OOCCD")==0 || strcmp(wfn,"ZAPTN")==0) 
  {
    return(1);
  }
  else  {
    return 0;
  }
}

double combinations	(	int	n,
		int	k
	)

combinations() : Calculates the number of ways to choose k objects from n objects, or "n choose k"

Parameters:

Parameters:

	n	= number of objects in total
	k	= number of objects taken at a time

Returns: number of combinations of n objects taken k at a time ("n choose k") (returned as a double).

Definition at line 23 of file probabil.cc.

References factorial().

{
   double factorial(int) ;
   double comb ;

   if (n == k) return (1.0) ;
   else if (k > n) return(0.0) ;
   else if (k == 0) return(1.0) ; 
   comb = factorial(n) / (factorial(k) * factorial(n-k)) ;
 
   return(comb) ;
}

int david	(	double **	A,
		int	N,
		int	M,
		double *	eps,
		double **	v,
		double	cutoff,
		int	print
	)

david(): Computes the lowest few eigenvalues and eigenvectors of a symmetric matrix, A, using the Davidson-Liu algorithm.

The matrix must be small enough to fit entirely in core. This algorithm is useful if one is interested in only a few roots of the matrix rather than the whole spectrum.

NB: This implementation will keep up to eight guess vectors for each root desired before collapsing to one vector per root. In addition, if smart_guess=1 (the default), guess vectors are constructed by diagonalization of a sub-matrix of A; otherwise, unit vectors are used.

TDC, July-August 2002

Parameters:

	A	= matrix to diagonalize
	N	= dimension of A
	M	= number of roots desired
	eps	= eigenvalues
	v	= eigenvectors
	cutoff	= tolerance for convergence of eigenvalues
	print	= Boolean for printing additional information

Returns: number of converged roots

Definition at line 47 of file david.cc.

References block_matrix(), C_DGEMM(), init_array(), init_int_array(), schmidt_add(), sq_rsp(), and zero_int_array().

{
  int i, j, k, L, I;
  double minimum;
  int min_pos, numf, iter, *conv, converged, maxdim, skip_check;
  int *small2big, init_dim;
  int smart_guess = 1;
  double *Adiag, **b, **bnew, **sigma, **G;
  double *lambda, **alpha, **f, *lambda_old;
  double norm, denom, diff;

  maxdim = 8 * M;

  b = block_matrix(maxdim, N);  /* current set of guess vectors,
                                   stored by row */
  bnew = block_matrix(M, N); /* guess vectors formed from old vectors,
                                stored by row*/
  sigma = block_matrix(N, maxdim); /* sigma vectors, stored by column */
  G = block_matrix(maxdim, maxdim); /* Davidson mini-Hamitonian */
  f = block_matrix(maxdim, N); /* residual eigenvectors, stored by row */
  alpha = block_matrix(maxdim, maxdim); /* eigenvectors of G */
  lambda = init_array(maxdim); /* eigenvalues of G */
  lambda_old = init_array(maxdim); /* approximate roots from previous
                                      iteration */

  if(smart_guess) { /* Use eigenvectors of a sub-matrix as initial guesses */

    if(N > 7*M) init_dim = 7*M;
    else init_dim = M;
    Adiag = init_array(N);
    small2big = init_int_array(7*M);
    for(i=0; i < N; i++) { Adiag[i] = A[i][i]; }
    for(i=0; i < init_dim; i++) {
      minimum = Adiag[0];
      min_pos = 0;
      for(j=1; j < N; j++)
        if(Adiag[j] < minimum) {        
          minimum = Adiag[j]; 
          min_pos = j; 
          small2big[i] = j; 
        }

      Adiag[min_pos] = BIGNUM;
      lambda_old[i] = minimum;
    }
    for(i=0; i < init_dim; i++) {
      for(j=0; j < init_dim; j++)
        G[i][j] = A[small2big[i]][small2big[j]];
    }

    sq_rsp(init_dim, init_dim, G, lambda, 1, alpha, 1e-12);

    for(i=0; i < init_dim; i++) {
      for(j=0; j < init_dim; j++)
        b[i][small2big[j]] = alpha[j][i];
    }

    free(Adiag);
    free(small2big);
  }
  else { /* Use unit vectors as initial guesses */
    Adiag = init_array(N);
    for(i=0; i < N; i++) { Adiag[i] = A[i][i]; }
    for(i=0; i < M; i++) {
      minimum = Adiag[0];
      min_pos = 0;
      for(j=1; j < N; j++)
        if(Adiag[j] < minimum) { minimum = Adiag[j]; min_pos = j; }

      b[i][min_pos] = 1.0; 
      Adiag[min_pos] = BIGNUM; 
      lambda_old[i] = minimum;
    }
    free(Adiag);
  }

  L = init_dim;
  iter =0;
  converged = 0;
  conv = init_int_array(M); /* boolean array for convergence of each
                               root */
  while(converged < M && iter < MAXIT) {

    skip_check = 0;
    if(print) printf("\niter = %d\n", iter); 

    /* form mini-matrix */
    C_DGEMM('n','t', N, L, N, 1.0, &(A[0][0]), N, &(b[0][0]), N,
            0.0, &(sigma[0][0]), maxdim);
    C_DGEMM('n','n', L, L, N, 1.0, &(b[0][0]), N,
            &(sigma[0][0]), maxdim, 0.0, &(G[0][0]), maxdim);

    /* diagonalize mini-matrix */
    sq_rsp(L, L, G, lambda, 1, alpha, 1e-12);

    /* form preconditioned residue vectors */
    for(k=0; k < M; k++)
      for(I=0; I < N; I++) {
        f[k][I] = 0.0;
        for(i=0; i < L; i++) {
          f[k][I] += alpha[i][k] * (sigma[I][i] - lambda[k] * b[i][I]);
        }
        denom = lambda[k] - A[I][I];
        if(fabs(denom) > 1e-6) f[k][I] /= denom;
        else f[k][I] = 0.0;
      }

    /* normalize each residual */
    for(k=0; k < M; k++) {
      norm = 0.0;
      for(I=0; I < N; I++) {
        norm += f[k][I] * f[k][I];
      }
      norm = sqrt(norm);
      for(I=0; I < N; I++) {
        if(norm > 1e-6) f[k][I] /= norm;
        else f[k][I] = 0.0;
      }
    }

    /* schmidt orthogonalize the f[k] against the set of b[i] and add
       new vectors */
    for(k=0,numf=0; k < M; k++)
      if(schmidt_add(b, L, N, f[k])) { L++; numf++; }

    /* If L is close to maxdim, collapse to one guess per root */
    if(maxdim - L < M) {
      if(print) {
        printf("Subspace too large: maxdim = %d, L = %d\n", maxdim, L);
        printf("Collapsing eigenvectors.\n");
      }
      for(i=0; i < M; i++) {
        memset((void *) bnew[i], 0, N*sizeof(double));
        for(j=0; j < L; j++) {
          for(k=0; k < N; k++) {
            bnew[i][k] += alpha[j][i] * b[j][k];
          }
        }
      }

      /* copy new vectors into place */
      for(i=0; i < M; i++) 
        for(k=0; k < N; k++)
          b[i][k] = bnew[i][k];

      skip_check = 1;

      L = M;
    }

    /* check convergence on all roots */
    if(!skip_check) {
      converged = 0;
      zero_int_array(conv, M);
      if(print) {
        printf("Root      Eigenvalue       Delta  Converged?\n");
        printf("---- -------------------- ------- ----------\n");
      }
      for(k=0; k < M; k++) {
        diff = fabs(lambda[k] - lambda_old[k]);
        if(diff < cutoff) {
          conv[k] = 1;
          converged++;
        }
        lambda_old[k] = lambda[k];
        if(print) {
          printf("%3d  %20.14f %4.3e    %1s\n", k, lambda[k], diff,
                 conv[k] == 1 ? "Y" : "N");
        }
      }
    }

    iter++;
  }

  /* generate final eigenvalues and eigenvectors */
  if(converged == M) {
    for(i=0; i < M; i++) {
      eps[i] = lambda[i];
      for(j=0; j < L; j++) {
        for(I=0; I < N; I++) {
          v[I][i] += alpha[j][i] * b[j][I];
        }
      }
    }
    if(print) printf("Davidson algorithm converged in %d iterations.\n", iter);
  }

  free(conv);
  free_block(b);
  free_block(bnew);
  free_block(sigma);
  free_block(G);
  free_block(f);
  free_block(alpha);
  free(lambda);
  free(lambda_old);

  return converged;
}

void dirprd_block	(	double **	A,
		double **	B,
		int	rows,
		int	cols
	)

dirprd_block()

This function takes two block matrices A and B and multiplies each element of B by the corresponding element of A

Parameters:

	A	= block matrix A
	B	= block matrix B
	nrows	= number of rows of A and B
	ncols	= number of columns of A and B

Returns: none

Definition at line 25 of file dirprd_block.cc.

{
  register long int i;
  double *a, *b;
  long size;

  size = ((long) rows) * ((long) cols);

  if(!size) return;

  a = A[0]; b= B[0];

  for(i=0; i < size; i++, a++, b++) (*b) = (*a) * (*b);
}

double dot_block	(	double **	A,
		double **	B,
		int	rows,
		int	cols,
		double	alpha
	)

dot_block(): Find dot product of two block matrices

Parameters:

	A	= block matrix A
	B	= block matrix B
	nrows	= number of rows of A and B
	ncols	= number of columns of A and B
	alpha	= scale factor by which the dot product is multiplied

Returns: dot product

Definition at line 21 of file dot_block.cc.

{
  register long int i;
  double *a, *b;
  double value;
  long int size;

  size = ((long) rows) * ((long) cols);

  if(!size) return 0.0;

  a = A[0]; b = B[0];

  value = 0.0;
  for(i=0; i < size; i++,a++,b++) value += (*a) * (*b);

  return alpha*value;
}

double eri	(	unsigned int	l1,
		unsigned int	m1,
		unsigned int	n1,
		double	alpha1,
		double	A[3],
		unsigned int	l2,
		unsigned int	m2,
		unsigned int	n2,
		double	alpha2,
		double	B[3],
		unsigned int	l3,
		unsigned int	m3,
		unsigned int	n3,
		double	alpha3,
		double	C[3],
		unsigned int	l4,
		unsigned int	m4,
		unsigned int	n4,
		double	alpha4,
		double	D[3],
		int	norm_flag
	)

eri()

This is a very inefficient function for computing ERIs over primitive Gaussian functions. The argument list is self-explanatory, except for norm_flag:

Parameters:

norm_flag

= tells what kind of normalization to use, 0 - no normalization, >0 - normalized ERI

Definition at line 39 of file eri.cc.

References block_matrix(), init_array(), and norm_const().

{

  const double gammap = alpha1 + alpha2;
  const double Px = (alpha1*A[0] + alpha2*B[0])/gammap;
  const double Py = (alpha1*A[1] + alpha2*B[1])/gammap;
  const double Pz = (alpha1*A[2] + alpha2*B[2])/gammap;
  const double PAx = Px - A[0];
  const double PAy = Py - A[1];
  const double PAz = Pz - A[2];
  const double PBx = Px - B[0];
  const double PBy = Py - B[1];
  const double PBz = Pz - B[2];
  const double AB2 = (A[0]-B[0])*(A[0]-B[0]) + (A[1]-B[1])*(A[1]-B[1]) + (A[2]-B[2])*(A[2]-B[2]); 
  
  const double gammaq = alpha3 + alpha4;
  const double gammapq = gammap*gammaq/(gammap+gammaq);
  const double Qx = (alpha3*C[0] + alpha4*D[0])/gammaq;
  const double Qy = (alpha3*C[1] + alpha4*D[1])/gammaq;
  const double Qz = (alpha3*C[2] + alpha4*D[2])/gammaq;
  const double QCx = Qx - C[0];
  const double QCy = Qy - C[1];
  const double QCz = Qz - C[2];
  const double QDx = Qx - D[0];
  const double QDy = Qy - D[1];
  const double QDz = Qz - D[2];
  const double CD2 = (C[0]-D[0])*(C[0]-D[0]) + (C[1]-D[1])*(C[1]-D[1]) + (C[2]-D[2])*(C[2]-D[2]);

  const double PQx = Px - Qx;
  const double PQy = Py - Qy;
  const double PQz = Pz - Qz;
  const double PQ2 = PQx*PQx + PQy*PQy + PQz*PQz;

  int u1,u2,v1,v2,w1,w2,tx,ty,tz,txmax,tymax,tzmax;
  int i,j,k;
  int lp,lq,mp,mq,np,nq;
  int zeta;
  double *flp, *flq, *fmp, *fmq, *fnp, *fnq;
  double *F;
  double K1, K2;
  double Gx,Gy,Gz;
  double pfac;
  double result = 0.0;
  double tmp;
  int u1max,u2max,v1max,v2max,w1max,w2max;

  K1 = exp(-alpha1*alpha2*AB2/gammap);
  K2 = exp(-alpha3*alpha4*CD2/gammaq);
  pfac = 2*pow(M_PI,2.5)*K1*K2/(gammap*gammaq*sqrt(gammap+gammaq));


  if (fac == NULL) {
    fac = init_array(MAXFAC);
    fac[0] = 1.0;
    for(i=1;i<MAXFAC;i++)
      fac[i] = fac[i-1]*i;
    bc = block_matrix(MAXFAC,MAXFAC);
    for(i=0;i<MAXFAC;i++)
      for(j=0;j<=i;j++)
        bc[i][j] = fac[i]/(fac[i-j]*fac[j]);
  }

  if (norm_flag > 0) {
    pfac *= norm_const(l1,m1,n1,alpha1,A);
    pfac *= norm_const(l2,m2,n2,alpha2,B);
    pfac *= norm_const(l3,m3,n3,alpha3,C);
    pfac *= norm_const(l4,m4,n4,alpha4,D);
  }
  

  F = init_array(l1+l2+l3+l4+m1+m2+m3+m4+n1+n2+n3+n4+1);
  calc_f(F,l1+l2+l3+l4+m1+m2+m3+m4+n1+n2+n3+n4,PQ2*gammapq);

  
  flp = init_array(l1+l2+1);
  for(k=0;k<=l1+l2;k++)
    for(i=0;i<=MIN(k,l1);i++) {
      j = k-i;
      if (j > l2) continue;
      tmp = bc[l1][i]*bc[l2][j];
      if (l1-i > 0)
        tmp *= pow(PAx,l1-i);
      if (l2-j > 0)
        tmp *= pow(PBx,l2-j);
      flp[k] += tmp;
    }
  fmp = init_array(m1+m2+1);
  for(k=0;k<=m1+m2;k++)
    for(i=0;i<=MIN(k,m1);i++) {
      j = k-i;
      if (j > m2) continue;
      tmp = bc[m1][i]*bc[m2][j];
      if (m1-i > 0)
        tmp *= pow(PAy,m1-i);
      if (m2-j > 0)
        tmp *= pow(PBy,m2-j);
      fmp[k] += tmp;
    }
  fnp = init_array(n1+n2+1);
  for(k=0;k<=n1+n2;k++)
    for(i=0;i<=MIN(k,n1);i++) {
      j = k-i;
      if (j > n2) continue;
      tmp = bc[n1][i]*bc[n2][j];
      if (n1-i > 0)
        tmp *= pow(PAz,n1-i);
      if (n2-j > 0)
        tmp *= pow(PBz,n2-j);
      fnp[k] += tmp;
    }
  flq = init_array(l3+l4+1);
  for(k=0;k<=l3+l4;k++)
    for(i=0;i<=MIN(k,l3);i++) {
      j = k-i;
      if (j > l4) continue;
      tmp = bc[l3][i]*bc[l4][j];
      if (l3-i > 0)
        tmp *= pow(QCx,l3-i);
      if (l4-j > 0)
        tmp *= pow(QDx,l4-j);
      flq[k] += tmp;
    }
  fmq = init_array(m3+m4+1);
  for(k=0;k<=m3+m4;k++)
    for(i=0;i<=MIN(k,m3);i++) {
      j = k-i;
      if (j > m4) continue;
      tmp = bc[m3][i]*bc[m4][j];
      if (m3-i > 0)
        tmp *= pow(QCy,m3-i);
      if (m4-j > 0)
        tmp *= pow(QDy,m4-j);
      fmq[k] += tmp;
    }
  fnq = init_array(n3+n4+1);
  for(k=0;k<=n3+n4;k++)
    for(i=0;i<=MIN(k,n3);i++) {
      j = k-i;
      if (j > n4) continue;
      tmp = bc[n3][i]*bc[n4][j];
      if (n3-i > 0)
        tmp *= pow(QCz,n3-i);
      if (n4-j > 0)
        tmp *= pow(QDz,n4-j);
      fnq[k] += tmp;
    }


  for(lp=0;lp<=l1+l2;lp++)
    for(lq=0;lq<=l3+l4;lq++) {
      u1max = lp/2;
      u2max = lq/2;
      for(u1=0;u1<=u1max;u1++)
        for(u2=0;u2<=u2max;u2++) {
          Gx = pow(-1,lp)*flp[lp]*flq[lq]*fac[lp]*fac[lq]*pow(gammap,u1-lp)*pow(gammaq,u2-lq)*fac[lp+lq-2*u1-2*u2]*
               pow(gammapq,lp+lq-2*u1-2*u2)/(fac[u1]*fac[u2]*fac[lp-2*u1]*fac[lq-2*u2]);
          for(mp=0;mp<=m1+m2;mp++)
            for(mq=0;mq<=m3+m4;mq++) {
              v1max = mp/2;
              v2max = mq/2;
              for(v1=0;v1<=v1max;v1++)
                for(v2=0;v2<=v2max;v2++) {
                  Gy = pow(-1,mp)*fmp[mp]*fmq[mq]*fac[mp]*fac[mq]*pow(gammap,v1-mp)*pow(gammaq,v2-mq)*fac[mp+mq-2*v1-2*v2]*
                       pow(gammapq,mp+mq-2*v1-2*v2)/(fac[v1]*fac[v2]*fac[mp-2*v1]*fac[mq-2*v2]);
                  for(np=0;np<=n1+n2;np++)
                    for(nq=0;nq<=n3+n4;nq++) {
                      w1max = np/2;
                      w2max = nq/2;
                      for(w1=0;w1<=w1max;w1++)
                        for(w2=0;w2<=w2max;w2++) {
                          Gz = pow(-1,np)*fnp[np]*fnq[nq]*fac[np]*fac[nq]*pow(gammap,w1-np)*pow(gammaq,w2-nq)*fac[np+nq-2*w1-2*w2]*
                               pow(gammapq,np+nq-2*w1-2*w2)/(fac[w1]*fac[w2]*fac[np-2*w1]*fac[nq-2*w2]);
                          txmax = (lp+lq-2*u1-2*u2)/2;
                          tymax = (mp+mq-2*v1-2*v2)/2;
                          tzmax = (np+nq-2*w1-2*w2)/2;
                          for(tx=0;tx<=txmax;tx++)
                          for(ty=0;ty<=tymax;ty++)
                          for(tz=0;tz<=tzmax;tz++) {
                            zeta = lp+lq+mp+mq+np+nq-2*u1-2*u2-2*v1-2*v2-2*w1-2*w2-tx-ty-tz;
                            result += Gx*Gy*Gz*F[zeta]*
                                      pow(-1,tx+ty+tz)*
                                      pow(PQx,lp+lq-2*u1-2*u2-2*tx)*
                                      pow(PQy,mp+mq-2*v1-2*v2-2*ty)*
                                      pow(PQz,np+nq-2*w1-2*w2-2*tz)/
                                      (pow(4,u1+u2+tx+v1+v2+ty+w1+w2+tz)*
                                      pow(gammapq,tx)*
                                      pow(gammapq,ty)*
                                      pow(gammapq,tz)*
                                      fac[lp+lq-2*u1-2*u2-2*tx]*fac[tx]*
                                      fac[mp+mq-2*v1-2*v2-2*ty]*fac[ty]*
                                      fac[np+nq-2*w1-2*w2-2*tz]*fac[tz]);
                          }
                        }
                    }
                }
            }
        }
    }

  free(F);
  free(flp);
  free(fmp);
  free(fnp);
  free(flq);
  free(fmq);
  free(fnq);

  return result*pfac;
}

double factorial

(

int

n

)

factorial(): Returns n!

Parameters:

Parameters:

n

= number to take factorial of

Returns: n factorial, as a double word (since n! can get very large).

Definition at line 47 of file probabil.cc.

References factorial().

Referenced by combinations(), and factorial().

{

   if (n == 0 || n == 1) return(1.0);
   if (n < 0) return(0.0) ;
   else {
      return ((double) n * factorial(n-1)) ;
      }
}

void fill_sym_matrix	(	double **	A,
		int	size
	)

fill_sym_matrix(): Fills a symmetric matrix by placing the elements of the lower triangle into the upper triangle.

Parameters:

	A	= matrix to symmetrize
	size	= number of rows or columns (assume square)

Returns: none

Definition at line 19 of file fill_sym_matrix.cc.

{
   double **row, *col; 
   int rc, cc;

   row = A;
   for (rc = 0; rc < (size-1); rc++) {
     col = *row;
     for (cc = 0; cc < size; cc++) {
       if (cc > rc) {
         *col = A[cc][rc];
       }
       col++;
     }
   row++;
  }
}

void filter	(	double *	input,
		double *	output,
		int *	ioff,
		int	norbs,
		int	nfzc,
		int	nfzv
	)

filter(): Filter out undesired (frozen core/virt) integrals

Given a lower-triangle array of integrals in the full space of orbitals as well as numbers of frozen core and virtual orbitals, this function returns a list of integrals involving only active orbitals.

TDC, June 2001

Note: Based on the code written by CDS in the original iwl_rd_one_all_act() function in LIBIWL.

Definition at line 25 of file filter.cc.

{
  int i, j, ij, IJ;
  int nact;

  nact = norbs - nfzc - nfzv;

  for(i=0,ij=0; i < nact; i++) {
    for(j=0; j <= i; j++,ij++) {
      IJ = ioff[i+nfzc] + (j + nfzc);
      output[ij] = input[IJ];
    }
  }
}

void free_3d_array	(	double ***	A,
		int	p,
		int	q
	)

free_3d_array(): Free a (non-contiguous) 3D array

Parameters:

	A	= triple-star pointer to the 3D array
	p	= size of first dimension
	q	= size of second dimension

Returns: none

Definition at line 55 of file 3d_array.cc.

{
  int i,j;

  for(i=0; i < p; i++)
    for(j=0; j < q; j++)
      free(A[i][j]);

  for(i=0; i < p; i++) free(A[i]);

  free(A);
}

int* get_frzcpi

(

)

get_frzcpi(): Get frozen core per irrep array.

Parameters: none

Returns: pointer to int array

Definition at line 42 of file get_frzpi.cc.

References chkpt_rd_frzcpi(), chkpt_rd_nirreps(), init_int_array(), ip_exist(), and ip_int_array().

Referenced by ras_set(), and ras_set2().

{
  int errcod;
  int nirreps, nirr;
  int need_to_init_psio = 0;
  int need_to_init_chkpt = 0;
  int if_exists;
  int* frzcpi;

PSIO_INIT
CHKPT_INIT(PSIO_OPEN_OLD);
  nirreps = chkpt_rd_nirreps();
  
  frzcpi = init_int_array(nirreps);
  
  if_exists = ip_exist("FROZEN_DOCC",0);
  if (if_exists) {
    errcod = ip_int_array("FROZEN_DOCC",frzcpi,nirreps);
  }
  else {
    free(frzcpi);
    frzcpi = chkpt_rd_frzcpi();
  }

CHKPT_DONE
PSIO_DONE

  return frzcpi;
}

int* get_frzvpi

(

)

get_frzvpi(): Get frozen virtuals per irrep array.

Parameters: none

Returns: pointer to int array

Definition at line 82 of file get_frzpi.cc.

References chkpt_rd_frzvpi(), chkpt_rd_nirreps(), init_int_array(), ip_exist(), and ip_int_array().

Referenced by ras_set(), and ras_set2().

{
  int errcod;
  int nirreps, nirr;
  int need_to_init_psio = 0;
  int need_to_init_chkpt = 0;
  int if_exists;
  int* frzvpi;

PSIO_INIT
CHKPT_INIT(PSIO_OPEN_OLD);
  nirreps = chkpt_rd_nirreps();
  
  frzvpi = init_int_array(nirreps);
  
  if_exists = ip_exist("FROZEN_UOCC",0);
  if (if_exists) {
    errcod = ip_int_array("FROZEN_UOCC",frzvpi,nirreps);
  }
  else {
    free(frzvpi);
    frzvpi = chkpt_rd_frzvpi();
  }

CHKPT_DONE
PSIO_DONE

  return frzvpi;
}

double*** init_3d_array	(	int	p,
		int	q,
		int	r
	)

3d_array(): Initialize a (non-contiguous) 3D array

Parameters:

	p	= size of first dimension
	q	= size of second dimension
	r	= size of third dimension

Returns: triple-star pointer to 3D array

Definition at line 24 of file 3d_array.cc.

{
  double ***A;
  int i,j,k;

  A = (double ***) malloc(p * sizeof(double **));
  for(i=0; i < p; i++) {
    A[i] = (double **) malloc(q * sizeof(double *));
    for(j=0; j < q; j++) {
      A[i][j] = (double *) malloc(r * sizeof(double));
      for(k=0; k < r; k++) {
        A[i][j][k] = 0.0;
      }
    }
  }

  return A;
}

double invert_matrix	(	double **	a,
		double **	y,
		int	N,
		FILE *	outfile
	)

INVERT_MATRIX(): The function takes the inverse of a matrix using the C routines in Numerical Recipes in C.

Matt Leininger, Summer 1994

Parameters:

Parameters:

	a	= matrix to take the inverse of (is modified by invert_matrix())
	y	= the inverse matrix
	N	= the size of the matrices
	outfile	= file for error messages

Other variables: col and indx are temporary arrays d is 1 or -1

Returns: double (determinant) Note: The original matrix is modified by invert_matrix()

Definition at line 38 of file invert.cc.

References init_array(), and init_int_array().

Referenced by psi::detcas::bfgs_hessian().

{
   double  d, *col, *colptr;
   register int i, j;
   int *indx ;

   col = init_array(N) ;
   indx = init_int_array(N) ;

   ludcmp(a,N,indx,&d) ;
   for (j=0; j<N; j++) d *= a[j][j];
  
   /* fprintf(outfile,"detH0 in invert = %lf\n", fabs(d));
   fflush(outfile); */

    if (fabs(d) < SMALL_DET) {
      fprintf(outfile,"Warning (invert_matrix): Determinant is %g\n", d);
      printf("Warning (invert_matrix): Determinant is %g\n", d);
      fflush(outfile);
      }

   for (j=0; j<N; j++) {
       bzero(col,sizeof(double)*N);
       col[j] = 1.0 ;
       lubksb(a,N,indx,col) ;
       colptr = col;
       for (i=0; i<N; i++) y[i][j] = *colptr++;
       }

   d = fabs(d);
   return(d);
}

int mat_in	(	FILE *	fp,
		double **	array,
		int	width,
		int	max_length,
		int *	stat
	)

MAT_IN(): Function to read in a matrix. Simple version for now.

Parameters:

Parameters:

	fp	= file pointer to input stream
	array	= matrix to hold data
	width	= number of columns to read
	max_length	= maximum number of rows to read
	stat	= pointer to int to hold status flag (0=read ok, 1=error)

Returns: number of rows read Also modifies stat to = error code (0 = ok, 1 = error)

Definition at line 28 of file mat_in.cc.

{
   int i=0, j, errcod=0 ;
   int nr ;
   double data ;

   while ( (i < max_length) && (!errcod) ) {
      for (j=0; j<width; j++) {
         nr = fscanf(fp, "%lf", &data) ;
         if (feof(fp)) break ;
         if (nr != 1) {
            errcod = 1 ;
            break ;
            }
         else {
            array[i][j] = data ;
            }
         }
      if (feof(fp)) break ;
      if (!errcod) i++ ;
      }

   *stat = errcod ;
   return(i) ;
}

int mat_print	(	double **	matrix,
		int	rows,
		int	cols,
		FILE *	outfile
	)

mat_print(): Prints a matrix to a file in a formatted style

Parameters:

	matrix	= matrix to print
	rows	= number of rows
	cols	= number of columns
	outfile	= output file pointer for printing

Returns: Always returns zero...

Definition at line 24 of file mat_print.cc.

{
  div_t fraction;
  int i,j;
  int cols_per_page, num_pages, last_page, page, first_col;

  cols_per_page = 5;

  /* Determine the number of cols_per_page groups */
  fraction = div(cols,cols_per_page);
  num_pages = fraction.quot;  /* Number of complete column groups */
  last_page = fraction.rem;  /* Number of columns in last group */

  /* Loop over the complete column groups */
  for(page=0; page < num_pages; page++) {
      first_col = page*cols_per_page;

      fprintf(outfile,"\n      ");
      for(i=first_col; i < first_col+cols_per_page; i++) 
          fprintf(outfile,"      %5d        ",i);

      fprintf (outfile,"\n");
      for(i=0; i < rows; i++) {
          fprintf(outfile,"\n%5d ",i);

          for(j=first_col; j < first_col+cols_per_page; j++)        
              fprintf (outfile,"%19.15f",matrix[i][j]);
        }

      fprintf (outfile,"\n");
    }

  /* Now print the remaining columns */
  if(last_page) {
      first_col = page*cols_per_page;

      fprintf(outfile,"\n      ");
      for(i=first_col; i < first_col+last_page; i++) 
          fprintf(outfile,"      %5d        ",i);
      
      fprintf (outfile,"\n");
      for(i=0; i < rows; i++) {
          fprintf(outfile,"\n%5d ",i);

          for(j=first_col; j < first_col+last_page; j++)
              fprintf (outfile,"%19.15f",matrix[i][j]);
        }

      fprintf (outfile,"\n");
    }

  return 0;

}

double newton	(	double(*)(double)	F,
		double(*)(double)	dF,
		double	x,
		double	tolerance,
		int	maxiter,
		int	printflag
	)

newton(): Find the root of a function by Newton's method. Iterations are limited to a maximum value. The algorithm stops when the difference between successive estimates of the root is less than the specified tolerance. An initial guess for the root must be given, as well as the function AND it's derivative.

Parameters:

	F	= pointer to function we want to examine (must return double)
	dF	= pointer to _derivative_ of function F
	x	= initial guess for root
	tolerance	= how close successive guesses must get before convergence
	maxiter	= maximum number of iterations
	printflag	= whether or not to print results for each iteration (1 or 0)

Returns: the value of the root

Definition at line 127 of file rootfind.cc.

{
   int iter = 1 ;
   double newx ;

   if (printflag) {
      printf("Newton's method root finder\n") ;
      printf("%4s %12s\n", "iter", "x") ;
      printf("%4d %12.8f\n", iter, x) ;
      }

   while (iter < maxiter) {
      newx = x - F(x) / dF(x) ;
      iter++ ;
      if (printflag) printf("%4d %12.8f\n", iter, newx) ;
      if (fabs(newx - x) < tolerance) return (newx) ;
      x = newx ;
      }

   printf("(newton): Failed to converge within %d iterations\n", iter) ;
   return(newx) ;
}

double norm_const	(	unsigned int	l1,
		unsigned int	m1,
		unsigned int	n1,
		double	alpha1,
		double	A[3]
	)

norm_const()

Returns: normalization constant

Definition at line 319 of file eri.cc.

References init_array().

Referenced by eri().

{
  int i;
  
  if (df == NULL) {
    df = init_array(2*MAXFAC);
    df[0] = 1.0;
    df[1] = 1.0;
    df[2] = 1.0;
    for(i=3; i<MAXFAC*2; i++) {
      df[i] = (i-1)*df[i-2];
    }
  }
  
  return pow(2*alpha1/M_PI,0.75)*pow(4*alpha1,0.5*(l1+m1+n1))/sqrt(df[2*l1]*df[2*m1]*df[2*n1]);
}

void normalize	(	double **	A,
		int	rows,
		int	cols
	)

normalize(): Normalize a set of vectors

Assume we're normalizing the ROWS

Parameters:

	A	= matrix holding vectors to normalize
	rows	= number of rows in A
	cols	= number of columns in A

Returns: none

David Sherrill, Feb 1994

Definition at line 30 of file normalize.cc.

References dot_arr().

{
  double normval;
  register int i, j;

  /* divide each row by the square root of its norm */
  for (i=0; i<rows; i++) {
    dot_arr(A[i], A[i], cols, &normval);
    normval = sqrt(normval);
    for (j=0; j<cols; j++) A[i][j] /= normval;
  }

}

int pople	(	double **	A,
		double *	x,
		int	dimen,
		int	num_vecs,
		double	tolerance,
		FILE *	outfile,
		int	print_lvl
	)

POPLE(): Uses Pople's method to iteratively solve linear equations Ax = b

Matt Leininger, April 1998

Parameters:

	A	= matrix
	x	= initially has vector b, but returns vector x.
	dimen	= dimension of vector x.
	num_vecs	= number of vectors x to obtain.
	tolerance	= cutoff threshold for norm of expansion vector.

Returns: 0 for success, 1 for failure

Definition at line 32 of file pople.cc.

References block_matrix(), dot_arr(), flin(), init_array(), print_mat(), zero_arr(), and zero_mat().

{
   double det, tval;
   double **Bmat; /* Matrix of expansion vectors */ 
   double **Ab;   /* Matrix of A x expansion vectors */ 
   double **M;    /* Pople subspace matrix */ 
   double *n;     /* overlap of b or transformed b in Ax = b 
                     with expansion vector zero */
   double *r;     /* residual vector */ 
   double *b;     /* b vector in Ax = b, or transformed b vector */  
   double *sign;  /* sign array  to insure diagonal element of A are positive */
   double **Mtmp; /* tmp M matrix passed to flin */
   register int i, j, L=0, I;
   double norm, rnorm=1.0, norm_Ab=1.0, *dotprod, *alpha;
   double *dvec;
   int llast=0, last=0, maxdimen;
   int quit=0;
   
   maxdimen = 200;
   /* initialize working array */
   Bmat = block_matrix(maxdimen, dimen);
   alpha = init_array(maxdimen);
   M = block_matrix(maxdimen, maxdimen);
   Mtmp = block_matrix(maxdimen, maxdimen);
   dvec = init_array(dimen);
   Ab = block_matrix(maxdimen, dimen);
   r = init_array(dimen);
   b = init_array(dimen);
   n = init_array(dimen);
   dotprod = init_array(maxdimen);
   sign = init_array(dimen);

   if (print_lvl > 6) {
   fprintf(outfile,"\n\n Using Pople's Method for solving linear equations.\n");
   fprintf(outfile,"     --------------------------------------------------\n");
   fprintf(outfile,"         Iter             Norm of Residual Vector      \n");
   fprintf(outfile,"        ------           -------------------------     \n");
   }

   norm = 0.0;
   for (i=0; i<dimen; i++) norm += x[i]*x[i];
   if (norm < ZERO) quit = 1;

   if (!quit) {
       for (i=0; i<dimen; i++) {
           if (A[i][i] > ZERO) sign[i] = 1.0;
           else sign[i] = -1.0;
           x[i] *= sign[i];
         }

       for (i=0; i<dimen; i++) { 
           Bmat[0][i] = x[i]; 
           b[i] = x[i];
           dvec[i] = sqrt(fabs(A[i][i]));
           b[i] /= dvec[i];
           /*   fprintf(outfile,"A[%d][%d] = %lf\n",i,i, A[i][i]);
                fprintf(outfile,"dvec[%d] = %lf\n",i, dvec[i]);
                fprintf(outfile,"x[%d] = %lf\n",i, x[i]);
           */
         }

       if (print_lvl > 8) {
           fprintf(outfile," A matrix in POPLE(LIBQT):\n");
           print_mat(A, dimen, dimen, outfile);
         }

       /* Constructing P matrix */
       for (i=0; i<dimen; i++) { 
           for (j=0; j<dimen; j++) {
               if (i==j) A[i][i] = 0.0;
               else A[i][j] = -A[i][j]*sign[i];
             }
         }
       if (print_lvl > 8) {
           fprintf(outfile," P matrix in POPLE(LIBQT):\n");
           print_mat(A, dimen, dimen, outfile);
         }

       /* Precondition P matrix with diagonal elements of A */
       for (i=0; i<dimen; i++)
           for (j=0; j<dimen; j++) {
               tval = 1.0/(dvec[j]*dvec[i]);
               A[i][j] *= tval;
             }

       if (print_lvl > 8) {
           fprintf(outfile," Preconditioned P matrix in POPLE(LIBQT):\n");
           print_mat(A, dimen, dimen, outfile);
         }

       /* 
          for (i=0; i<dimen; i++) { 
          for (j=0; j<dimen; j++) {
          if (i==j) A[i][i] = 1.0 - A[i][i]; 
          else A[i][j] = -A[i][j];
          }
          }
       */

       for (i=0; i<dimen; i++) Bmat[0][i] /= dvec[i];


       dot_arr(Bmat[0],Bmat[0],dimen,&norm);
       norm = sqrt(norm);
       for (i=0; i<dimen; i++) {
           x[i] = Bmat[0][i];
           Bmat[0][i] /= norm;
         }

       dot_arr(Bmat[0],x,dimen,&(n[0]));
 
       while (!last) {
           /* form A*b_i */
           for (i=0; i<dimen; i++) 
               dot_arr(A[i], Bmat[L], dimen, &(Ab[L][i]));


           /* Construct M matrix */
           /* M = delta_ij - <b_new|Ab_j> */
           zero_mat(M, maxdimen, maxdimen);
           for (i=0; i<=L; i++) {
               for (j=0; j<=L; j++) {
                   dot_arr(Bmat[i], Ab[j], dimen, &(dotprod[i]));
                   if (i==j) M[i][j] = 1.0 - dotprod[i]; 
                   else M[i][j] = -dotprod[i];
                 }
             }

           if (llast) last = llast;

           for (i=0; i<=L; i++) { 
               alpha[i] = n[i];
               for (j=0; j<=L; j++)
                   Mtmp[i][j] = M[i][j];
             }

           flin(Mtmp, alpha, L+1, 1, &det);
 
 
           /* Need to form and backtransform x to orig basis x = D^(-1/2) x */
           zero_arr(x, dimen);
           for (i=0; i<=L; i++)
               for (I=0; I<dimen; I++)
                   x[I] += alpha[i]*Bmat[i][I]/dvec[I];

           /* Form residual vector Ax - b = r */
           zero_arr(r, dimen);
           for (I=0; I<dimen; I++)
               for (i=0; i<=L; i++) 
                   r[I] += alpha[i]*(Bmat[i][I] - Ab[i][I]);

           for (I=0; I<dimen; I++) r[I] -= b[I]; 
 
           dot_arr(r, r, dimen, &rnorm);
           rnorm = sqrt(rnorm);
           if (print_lvl > 6) {
               fprintf(outfile,
                 "        %3d                     %10.3E\n",L+1,rnorm);
               fflush(outfile);
             }

           if (L+1>dimen) {
               fprintf(outfile,"POPLE: Too many vectors in expansion space.\n");
               return 1;
             }

           /* place residual in b vector space */
           if (L+1>= maxdimen) {
               fprintf(outfile,
                 "POPLE (LIBQT): Number of expansion vectors exceeds"
                       " maxdimen (%d)\n", L+1);
               return 1;
             }
           for (i=0; i<dimen; i++) Bmat[L+1][i] = Ab[L][i];  
           /* 
              for (i=0; i<dimen; i++) Bmat[L+1][i] = r[i]; 
              for (i=0; i<dimen; i++) Bmat[L+1][i] = r[i]/(dvec[i]*dvec[i]); 
           */


           /* Schmidt orthonormalize new b vec to expansion space */
           for (j=0; j<=L; j++) {
               dot_arr(Bmat[j], Bmat[L+1], dimen, &(dotprod[j]));
               for (I=0; I<dimen; I++)
                   Bmat[L+1][I] -= dotprod[j] * Bmat[j][I];
             }

           /* Normalize new expansion vector */
           dot_arr(Bmat[L+1], Bmat[L+1], dimen, &norm);  
           norm = sqrt(norm);
           for (I=0; I<dimen; I++) Bmat[L+1][I] /= norm;

           /* check orthogonality of subspace */
           if (0) {
               for (i=0; i<=L+1; i++) { 
                   for (j=0; j<=i; j++) { 
                       dot_arr(Bmat[i], Bmat[j], dimen, &tval);
                       fprintf(outfile, "Bvec[%d] * Bvec[%d] = %lf\n",i,j,tval); 
                     }
                 }
             } 


           if (rnorm<tolerance) llast = last = 1;
           L++;
         }

       zero_arr(x, dimen);
       for (i=0; i<=L; i++) 
           for (I=0; I<dimen; I++)
               x[I] += alpha[i]*Bmat[i][I];
   
       /* Need to backtransform x to orginal basis x = D^(-1/2) x */
       for (I=0; I<dimen; I++)
           x[I] /= dvec[I]; 

     }
   free_block(Bmat);
   free(alpha);
   free_block(M);
   free_block(Mtmp);
   free(n);
   free(dvec);
   free_block(Ab);
   free(r);
   free(b);
   free(dotprod);

   return 0;
}

int ras_set	(	int	nirreps,
		int	nbfso,
		int	freeze_core,
		int *	orbspi,
		int *	docc,
		int *	socc,
		int *	frdocc,
		int *	fruocc,
		int **	ras_opi,
		int *	order,
		int	ras_type
	)

ras_set(): Deprecated

This function sets up the number of orbitals per irrep for each of the RAS subspaces [frozen core, RAS I, RAS II, RAS III, RAS IV, frozen virts]. It also obtains the appropriate orbital reordering array. The reordering array takes a basis function in Pitzer ordering (orbitals grouped according to irrep) and gives the corresponding index in the RAS numbering scheme. Orbitals are numbered according to irrep within each of the subspaces.

Formerly, docc, socc, frdocc, and fruocc were read in this function. Now docc and socc will be left as-is if they are not present in input.

C. David Sherrill Center for Computational Quantum Chemistry University of Georgia, 25 June 1995

Parameters:

	nirreps	= num of irreps in computational point group
	nbfso	= num of basis functions in symmetry orbitals (num MOs)
	freeze_core	= 1 to remove frozen core orbitals from ras_opi
	orbspi	= array giving num symmetry orbitals (or MOs) per irrep
	docc	= array of doubly occupied orbitals per irrep
	socc	= array of singly occupied orbitals per irrep
	frdocc	= array of frozen core per irrep
	fruocc	= array of frozen virtuals per irrep
	ras_opi	= matrix giving num of orbitals per irrep per ras space, addressed as ras_opi[ras_space][irrep]
	order	= array nbfso big which maps Pitzer to Correlated order
	ras_type	= if 1, put docc and socc together in same RAS space (RAS I), as appropriate for DETCI. If 0, put socc in its own RAS space (RAS II), as appropriate for CC.

Returns: 1 for success, 0 otherwise

Definition at line 54 of file ras_set.cc.

References free_int_matrix(), get_frzcpi(), get_frzvpi(), init_int_array(), init_int_matrix(), ip_exist(), ip_int_array(), and zero_int_array().

{
  int i, irrep, point, tmpi, cnt=0;
  int errcod, errbad, parsed_ras1=0, parsed_ras2=0, do_ras4;
  int *used, *offset, **tras;
  int *tmp_frdocc, *tmp_fruocc;

  used = init_int_array(nirreps);
  offset = init_int_array(nirreps);
  
  /* if we have trouble reading DOCC and SOCC, we'll take them as
   * provided to this routine.  Zero out everything else
   */
  for (i=0; i<4; i++) {
    zero_int_array(ras_opi[i], nirreps);
  }
  zero_int_array(order, nbfso);
  
  
  /* now use the parser to get the arrays we require */
  tmp_frdocc = get_frzcpi();
  tmp_fruocc = get_frzvpi();
  for (i=0; i<nirreps; i++) {
    frdocc[i] = tmp_frdocc[i];
    fruocc[i] = tmp_fruocc[i];
  }
  free(tmp_frdocc);
  free(tmp_fruocc);

  /* replace DOCC and SOCC only if they are in input */
  if (ip_exist("DOCC",0)) 
    errcod = ip_int_array("DOCC",docc,nirreps); 
  if (ip_exist("DOCC",0))
    errcod = ip_int_array("SOCC",socc,nirreps);

  errbad=0; do_ras4=1;
  errcod = ip_int_array("RAS1", ras_opi[0], nirreps);
  if (errcod == IPE_OK) parsed_ras1 = 1;
  else if (errcod == IPE_KEY_NOT_FOUND) {
    errcod = ip_int_array("ACTIVE_CORE", ras_opi[0], nirreps);
    if (errcod == IPE_OK) parsed_ras1 = 1;
    else if (errcod != IPE_KEY_NOT_FOUND) errbad = 1;
  }
  else errbad = 1;
  errcod = ip_int_array("RAS2", ras_opi[1], nirreps);
  if (errcod == IPE_OK) parsed_ras2 = 1;
  else if (errcod == IPE_KEY_NOT_FOUND) {
    errcod = ip_int_array("MODEL_SPACE", ras_opi[1], nirreps);
    if (errcod == IPE_OK) parsed_ras2 = 1;
    else if (errcod != IPE_KEY_NOT_FOUND) errbad = 1;
  }
  else errbad = 1;
  errcod = ip_int_array("RAS3", ras_opi[2], nirreps);
  if (errcod != IPE_OK && errcod != IPE_KEY_NOT_FOUND) errbad=1;
  if (errcod == IPE_KEY_NOT_FOUND) do_ras4=0;
  
  if (errbad == 1) {
    fprintf(stderr, "(ras_set): trouble parsing RAS keyword\n");
    return(0);
  }
  
  /* if the user has not specified RAS I, we must deduce it.
   * RAS I does not include any FZC orbs but does include COR orbs
   */
  
  if (!parsed_ras1) {
    for (irrep=0; irrep<nirreps; irrep++) {
      if (ras_type==1) ras_opi[0][irrep] = docc[irrep] + socc[irrep];
      if (ras_type==0) ras_opi[0][irrep] = docc[irrep];
      ras_opi[0][irrep] -= frdocc[irrep]; /* add back later for COR */ 
    }
  }
  
  
  /* if the user hasn't specified RAS II, look for val_orb             */
  /* val_orb should typically be RAS I + RAS II, so subtract out RAS I */ 
  if (!parsed_ras2) {
    errcod = ip_int_array("VAL_ORB",ras_opi[1],nirreps);
    if (errcod != IPE_OK) {
      for (irrep=0; irrep<nirreps; irrep++) {
        if (ras_type==1) ras_opi[1][irrep] = 0;
        if (ras_type==0) ras_opi[1][irrep] = socc[irrep];
      }
    }
    else {
      for (irrep = 0; irrep<nirreps; irrep++) { 
         ras_opi[1][irrep] -= ras_opi[0][irrep];
         if (ras_opi[1][irrep] < 0) {
           fprintf(stderr, "(ras_set): val_orb must be larger than RAS I\n");
           return(0);
         }
      }
    }
  }
  
  /* set up the RAS III or IV array: if RAS IV is used, RAS III must
   * be specified and then RAS IV is deduced.
   */
  
  for (irrep=0; irrep<nirreps; irrep++) { 
    tmpi = frdocc[irrep] + fruocc[irrep] + ras_opi[0][irrep] + 
      ras_opi[1][irrep];
    if (do_ras4) tmpi += ras_opi[2][irrep];
    if (tmpi > orbspi[irrep]) {
      fprintf(stderr, "(ras_set): orbitals don't add up for irrep %d\n",
              irrep);
      return(0);
    }
    if (do_ras4) ras_opi[3][irrep] = orbspi[irrep] - tmpi;
    else ras_opi[2][irrep] = orbspi[irrep] - tmpi;
  } 
  
  /* copy everything to the temporary RAS arrays: */
  /* add subspaces for frozen orbitals            */
  tras = init_int_matrix(6, nirreps);
  for (irrep=0; irrep<nirreps; irrep++) {
    tras[0][irrep] = frdocc[irrep];
    tras[5][irrep] = fruocc[irrep];
  }
  for (i=0; i<4; i++) {
    for (irrep=0; irrep<nirreps; irrep++) {
      tras[i+1][irrep] = ras_opi[i][irrep];
    }
  }
 
  /* construct the offset array */
  offset[0] = 0;
  for (irrep=1; irrep<nirreps; irrep++) {
    offset[irrep] = offset[irrep-1] + orbspi[irrep-1];
  }
  
  for (i=0; i<6; i++) {
    for (irrep=0; irrep<nirreps; irrep++) { 
      while (tras[i][irrep]) {
        point = used[irrep] + offset[irrep];
        if (point < 0 || point >= nbfso) {
          fprintf(stderr, "(ras_set): Invalid point value\n");
          exit(PSI_RETURN_FAILURE);
        }
        order[point] = cnt++;
        used[irrep]++;
        tras[i][irrep]--;
      }
    }
  }
  

  /* do a final check */

  for (irrep=0; irrep<nirreps; irrep++) {
    if (used[irrep] > orbspi[irrep]) {
      fprintf(stderr, "(ras_set): on final check, used more orbitals");
      fprintf(stderr, "   than were available (%d vs %d) for irrep %d\n",
              used[irrep], orbspi[irrep], irrep);
      errbad = 1;
    } 
    if (used[irrep] < orbspi[irrep]) {
      fprintf(stderr, "(ras_set): on final check, used fewer orbitals");
      fprintf(stderr, "   than were available (%d vs %d) for irrep %d\n",
              used[irrep], orbspi[irrep], irrep);
      errbad = 1;
    }
  }
  
  /* for restricted COR orbitals */
  if (!freeze_core) {
    for (irrep=0; irrep<nirreps; irrep++) {
      ras_opi[0][irrep] += frdocc[irrep];
    }
  }
  
  free(used);  free(offset);
  free_int_matrix(tras);

  return(!errbad);

}

int ras_set2	(	int	nirreps,
		int	nmo,
		int	delete_fzdocc,
		int	delete_restrdocc,
		int *	orbspi,
		int *	docc,
		int *	socc,
		int *	frdocc,
		int *	fruocc,
		int *	restrdocc,
		int *	restruocc,
		int **	ras_opi,
		int *	order,
		int	ras_type,
		int	hoffmann
	)

ras_set2()

This function sets up the number of orbitals per irrep for each of the RAS subspaces [frozen core (FZC), restricted core (COR), RAS I, RAS II, RAS III, RAS IV, restricted virts (VIR), frozen virts (FZV)]. It also obtains the appropriate orbital reordering array. The reordering array takes a basis function in Pitzer ordering (orbitals grouped according to irrep) and gives the corresponding index in the RAS numbering scheme. Orbitals are numbered according to irrep within each of the subspaces.

Formerly, docc, socc, frdocc, and fruocc were read in this function. Now docc and socc will be left as-is if they are not present in input.

Assume we always want integrals (at least some of them...) involving restricted orbitals, but we may not need them for frozen orbitals unless perhaps it's a gradient. The frozen orbitals will never enter explicitly in DETCI, but restricted orbitals may or may not, depending on the type of computation (CI vs CAS vs CASPT2, etc). For conventional CI, FCI, MRCI, RASCI, CASSCF, restricted orbitals will not participate explicitly in DETCI. For CASPT2, perhaps they will. CLAG and DETCAS will still need to have some indices in the restricted (and possibly frozen) orbital subspaces?

C. David Sherrill

Updated June 2002 to distinguish between frozen and restricted spaces

Parameters:

	nirreps	= num of irreps in computational point group
	nmo	= number of MO's
	delete_fzdocc	= 1 to remove frozen core orbitals from ras_opi[0]
	delete_restrdocc	= 1 to remove restricted core orbs from ras_opi[0]
	orbspi	= array giving num symmetry orbitals (or MOs) per irrep
	docc	= array of doubly occupied orbitals per irrep (guess should be provided)
	socc	= array of singly occupied orbitals per irrep (guess should be provided)
	frdocc	= array of frozen core per irrep (returned by function, but allocate before call)
	fruocc	= array of frozen virtuals per irrep (returned by function, but allocate before call)
	rstrdocc	= array of restricted core per irrep (returned by function, but allocate before call)
	rstruocc	= array of restricted core per irrep (returned by function, but allocate before call)
	ras_opi	= matrix giving num of orbitals per irrep per ras space, addressed as ras_opi[ras_space][irrep] (returned by function, but allocate before call)
	order	= array nmo big which maps Pitzer to Correlated order (returned by function, but allocate before call)
	ras_type	= if 1, put docc and socc together in same RAS space (RAS I), as appropriate for DETCI. If 0, put socc in its own RAS space (RAS II), as appropriate for CC.
	hoffmann	= if > 0, order orbitals according to Mark Hoffmann. hoffmann==1: ras1, ras2, ..., rasn, COR, FZC, VIR, FZV. hoffmann==2: VIR, ras1, ras2, ..., rasn, COR, FZC, FZV. Note odd placement of FZC in middle!

Returns: 1 for success, 0 otherwise

Definition at line 299 of file ras_set.cc.

References free_int_matrix(), get_frzcpi(), get_frzvpi(), init_int_array(), init_int_matrix(), ip_exist(), ip_int_array(), and zero_int_array().

{
  int i, irrep, point, tmpi, cnt=0;
  int errcod, errbad, parsed_ras1=0, parsed_ras2=0, do_ras4;
  int parsed_restr_uocc=0;
  int *used, *offset, **tras;
  int *tmp_frdocc, *tmp_fruocc;

  /* Hoffmann's code never wants FZC or COR lumped in with RAS I,
     even if those orbitals need to be transformed also */
  if (hoffmann > 0) {
    delete_fzdocc = 1;
    delete_restrdocc = 1;
  }

  used = init_int_array(nirreps);
  offset = init_int_array(nirreps);
  
  /* if we have trouble reading DOCC and SOCC, we'll take them as
   * provided to this routine.  Zero out everything else
   */
  zero_int_array(restrdocc, nirreps);
  zero_int_array(restruocc, nirreps);

  for (i=0; i<MAX_RAS_SPACES; i++) {
    zero_int_array(ras_opi[i], nirreps);
  }
  zero_int_array(order, nmo);
  
  /* now use the parser to get the arrays we require */
  tmp_frdocc = get_frzcpi();
  tmp_fruocc = get_frzvpi();
  for (i=0; i<nirreps; i++) {
    frdocc[i] = tmp_frdocc[i];
    fruocc[i] = tmp_fruocc[i];
  }
  free(tmp_frdocc);
  free(tmp_fruocc);


  /* replace DOCC and SOCC only if they are in input */
  if (ip_exist("DOCC",0))
    errcod = ip_int_array("DOCC",docc,nirreps); 
  if (ip_exist("SOCC",0))
    errcod = ip_int_array("SOCC",socc,nirreps);
  
  /* now use the parser to get the arrays we require */
  errcod = ip_int_array("RESTRICTED_DOCC",restrdocc,nirreps);
  
  errbad=0; do_ras4=1;
  errcod = ip_int_array("RAS1", ras_opi[0], nirreps);
  if (errcod == IPE_OK) parsed_ras1 = 1;
  else if (errcod != IPE_KEY_NOT_FOUND) errbad = 1;
  errcod = ip_int_array("RAS2", ras_opi[1], nirreps);
  if (errcod == IPE_OK) parsed_ras2 = 1;
  else if (errcod != IPE_KEY_NOT_FOUND) errbad = 1;
  errcod = ip_int_array("RAS3", ras_opi[2], nirreps);
  if (errcod != IPE_OK && errcod != IPE_KEY_NOT_FOUND) errbad=1;
  if (errcod == IPE_KEY_NOT_FOUND) do_ras4=0;
  
  if (errbad == 1) {
    fprintf(stderr, "(ras_set): trouble parsing RAS keyword\n");
    return(0);
  }
  
  /*
   * If the user has not specified RAS I, we must deduce it.
   * As of 2004, RAS I will not include any FZC ("frozen core")
   * orbitals *OR* any COR ("restricted core") orbitals.  However,
   * we'll add these back into RAS I later if the user requests it
   * by setting delete_fzdocc or delete_rstrdocc = 0.
   */
  
  if (!parsed_ras1) {
    for (irrep=0; irrep<nirreps; irrep++) {
      if (ras_type==1) ras_opi[0][irrep] = docc[irrep] + socc[irrep];
      if (ras_type==0) ras_opi[0][irrep] = docc[irrep];
      ras_opi[0][irrep] -= (frdocc[irrep] + restrdocc[irrep]);
    }
  }
  
  
  /* 
   * if the user hasn't specified RAS II, look for ACTIVE (replaces
   * VAL_ORB).  Assume this always goes after FZC + COR. 
   */ 
  if (!parsed_ras2) {
    errcod = ip_int_array("ACTIVE",ras_opi[1],nirreps);
    if (errcod != IPE_OK) { /* we didn't find ACTIVE keyword */
      for (irrep=0; irrep<nirreps; irrep++) {
        if (ras_type==1) ras_opi[1][irrep] = 0;
        if (ras_type==0) ras_opi[1][irrep] = socc[irrep];
      }
    }
    else { /* we found a ACTIVE keyword */
      if (parsed_ras1) 
        fprintf(outfile, "ras_set(): Warning: ACTIVE overrides RAS 1\n");

      for (irrep=0; irrep<nirreps; irrep++) 
        ras_opi[0][irrep] = 0; /* ACTIVE overrides RAS 1 */

      if (!ip_exist("RESTRICTED_UOCC",0)) { /* default restrict other virs */
        for (irrep=0; irrep<nirreps; irrep++) {
          ras_opi[0][irrep] = 0; /* ACTIVE overrides RAS 1 */
          restruocc[irrep] = orbspi[irrep] - frdocc[irrep] - 
                             restrdocc[irrep] - ras_opi[0][irrep] - 
                             ras_opi[1][irrep] - fruocc[irrep];
        }
        parsed_restr_uocc = 1;
      } 

      /* do a quick check */
      for (irrep=0; irrep<nirreps; irrep++) {
        if (frdocc[irrep]+restrdocc[irrep]+ras_opi[0][irrep]+ras_opi[1][irrep] 
            < docc[irrep] + socc[irrep])
          fprintf(outfile, 
                  "ras_set():Warning:Occupied electrons beyond ACTIVE orbs!\n");
      }
    }
  } /* end looking for "ACTIVE" */

  if (!parsed_restr_uocc) 
    errcod = ip_int_array("RESTRICTED_UOCC",restruocc,nirreps);
  
  /* set up the RAS III or IV array: if RAS IV is used, RAS III must
   * be specified and then RAS IV is deduced.
   */
  
  for (irrep=0; irrep<nirreps; irrep++) { 
    tmpi = frdocc[irrep] + restrdocc[irrep] + fruocc[irrep] + restruocc[irrep]
           + ras_opi[0][irrep] + ras_opi[1][irrep];
    if (do_ras4) tmpi += ras_opi[2][irrep];
    if (tmpi > orbspi[irrep]) {
      fprintf(stderr, "(ras_set): orbitals don't add up for irrep %d\n",
              irrep);
      return(0);
    }
    if (do_ras4) ras_opi[3][irrep] = orbspi[irrep] - tmpi;
    else ras_opi[2][irrep] = orbspi[irrep] - tmpi;
  } 
  
  /* copy everything to the temporary RAS arrays: */
  /* add subspaces for frozen orbitals            */
  tras = init_int_matrix(MAX_RAS_SPACES+4, nirreps);

  /* for usual DETCI, DETCAS, etc */
  if (hoffmann==0) {
    for (irrep=0; irrep<nirreps; irrep++) {
      tras[0][irrep] = frdocc[irrep];
      tras[1][irrep] = restrdocc[irrep];
      tras[MAX_RAS_SPACES+2][irrep] = restruocc[irrep];
      tras[MAX_RAS_SPACES+3][irrep] = fruocc[irrep];
    }
    for (i=0; i<MAX_RAS_SPACES; i++) {
      for (irrep=0; irrep<nirreps; irrep++) {
        tras[i+2][irrep] = ras_opi[i][irrep];
      }
    }
  }
  /* for Mark Hoffmann's UND order for GVVPT2, etc. */
  else if (hoffmann == 1) {
    for (i=0; i<MAX_RAS_SPACES; i++) {
      for (irrep=0; irrep<nirreps; irrep++) {
        tras[i][irrep] = ras_opi[i][irrep];
      }
    }
    for (irrep=0; irrep<nirreps; irrep++) {
      tras[MAX_RAS_SPACES+0][irrep] = restrdocc[irrep];
      tras[MAX_RAS_SPACES+1][irrep] = frdocc[irrep];
      tras[MAX_RAS_SPACES+2][irrep] = restruocc[irrep];
      tras[MAX_RAS_SPACES+3][irrep] = fruocc[irrep];
    }
  }
  /* for Mark Hoffmann's UND GUGA code */
  else if (hoffmann == 2) {
    for (i=0; i<MAX_RAS_SPACES; i++) {
      for (irrep=0; irrep<nirreps; irrep++) {
        tras[i+1][irrep] = ras_opi[i][irrep];
      }
    }
    for (irrep=0; irrep<nirreps; irrep++) {
      tras[0][irrep] = restruocc[irrep];
      tras[1+MAX_RAS_SPACES][irrep] = restrdocc[irrep];
      tras[2+MAX_RAS_SPACES][irrep] = frdocc[irrep];
      tras[MAX_RAS_SPACES+3][irrep] = fruocc[irrep];
    }
  }
   
  /* construct the offset array */
  offset[0] = 0;
  for (irrep=1; irrep<nirreps; irrep++) {
    offset[irrep] = offset[irrep-1] + orbspi[irrep-1];
  }
  
  for (i=0; i<MAX_RAS_SPACES+4; i++) {
    for (irrep=0; irrep<nirreps; irrep++) { 
      while (tras[i][irrep]) {
        point = used[irrep] + offset[irrep];
        if (point < 0 || point >= nmo) {
          fprintf(stderr, "(ras_set): Invalid point value\n");
          abort();
        }
        order[point] = cnt++;
        used[irrep]++;
        tras[i][irrep]--;
      }
    }
  }
  

  /* do a final check */
  for (irrep=0; irrep<nirreps; irrep++) {
    if (used[irrep] > orbspi[irrep]) {
      fprintf(stderr, "(ras_set): on final check, used more orbitals");
      fprintf(stderr, "   than were available (%d vs %d) for irrep %d\n",
              used[irrep], orbspi[irrep], irrep);
      errbad = 1;
    } 
    if (used[irrep] < orbspi[irrep]) {
      fprintf(stderr, "(ras_set): on final check, used fewer orbitals");
      fprintf(stderr, "   than were available (%d vs %d) for irrep %d\n",
              used[irrep], orbspi[irrep], irrep);
      errbad = 1;
    }
  }
  
  /* add back FZC and COR orbitals to RAS I if they are to be included */
  for (irrep=0; irrep<nirreps; irrep++) {
    if (!delete_restrdocc) ras_opi[0][irrep] += restrdocc[irrep];
    if (!delete_fzdocc)    ras_opi[0][irrep] += frdocc[irrep];
  }
  
  free(used);  free(offset);
  free_int_matrix(tras);

  return(!errbad);

}

void reorder_qt	(	int *	docc_in,
		int *	socc_in,
		int *	frozen_docc_in,
		int *	frozen_uocc_in,
		int *	order,
		int *	orbs_per_irrep,
		int	nirreps
	)

reorder_qt()

This function constructs a reordering array according to the "Quantum Trio" standard ordering, in which the orbitals are divided into the following sets: frozen core, then doubly occupied, then singly occupied, then virtuals, then deleted (frozen) virtuals. The reordering array takes a basis function in Pitzer ordering (orbitals grouped according to irrep) and gives the corresponding index in the Quantum Trio numbering scheme.

Should give the same reordering array as in the old libread30 routines.

C. David Sherrill Center for Computational Quantum Chemistry University of Georgia, 1995

Parameters:

	docc_in	= doubly occupied orbitals per irrep
	socc_in	= singly occupied orbitals per irrep
	frozen_docc_in	= frozen occupied orbitals per irrep
	frozen_uocc_in	= frozen unoccupied orbitals per irrep
	order	= reordering array (Pitzer->QT order)
	nirreps	= number of irreducible representations

Definition at line 40 of file reorder_qt.cc.

References init_int_array().

{

   int cnt=0, irrep, point, tmpi;
   int *used, *offset; 
   int *docc, *socc, *frozen_docc, *frozen_uocc; 
   int *uocc;

   used = init_int_array(nirreps);
   offset = init_int_array(nirreps);

   docc = init_int_array(nirreps);
   socc = init_int_array(nirreps);
   frozen_docc = init_int_array(nirreps);
   frozen_uocc = init_int_array(nirreps);
   uocc = init_int_array(nirreps);

   for (irrep=0; irrep<nirreps; irrep++) {
      docc[irrep] = docc_in[irrep];
      socc[irrep] = socc_in[irrep];
      frozen_docc[irrep] = frozen_docc_in[irrep];
      frozen_uocc[irrep] = frozen_uocc_in[irrep];
      }
 
   /* construct the offset array */
   offset[0] = 0;
   for (irrep=1; irrep<nirreps; irrep++) {
      offset[irrep] = offset[irrep-1] + orbs_per_irrep[irrep-1];
      }
   
   /* construct the uocc array */
   for (irrep=0; irrep<nirreps; irrep++) {
      tmpi = frozen_uocc[irrep] + docc[irrep] + socc[irrep];
      if (tmpi > orbs_per_irrep[irrep]) {
         fprintf(stderr, "(reorder_qt): orbitals don't add up for irrep %d\n",
            irrep);
         return;
         }
      else
         uocc[irrep] = orbs_per_irrep[irrep] - tmpi;
      }

   /* do the frozen core */
   for (irrep=0; irrep<nirreps; irrep++) { 
      while (frozen_docc[irrep]) {
         point = used[irrep] + offset[irrep];
         order[point] = cnt++;
         used[irrep]++;
         frozen_docc[irrep]--;
         docc[irrep]--;
         }
      }

   /* do doubly occupied orbitals */
   for (irrep=0; irrep<nirreps; irrep++) {
      while (docc[irrep]) {
         point = used[irrep] + offset[irrep];
         order[point] = cnt++;
         used[irrep]++;
         docc[irrep]--;
         }
      }

   /* do singly-occupied orbitals */
   for (irrep=0; irrep<nirreps; irrep++) {
      while (socc[irrep]) {
         point = used[irrep] + offset[irrep];
         order[point] = cnt++;
         used[irrep]++;
         socc[irrep]--;
         }
      }

   /* do virtual orbitals */
   for (irrep=0; irrep<nirreps; irrep++) {
      while (uocc[irrep]) {
         point = used[irrep] + offset[irrep];
         order[point] = cnt++;
         used[irrep]++;
         uocc[irrep]--;
         }
      }

   /* do frozen uocc */
   for (irrep=0; irrep<nirreps; irrep++) {
      while (frozen_uocc[irrep]) {
         point = used[irrep] + offset[irrep];
         order[point] = cnt++;
         used[irrep]++;
         frozen_uocc[irrep]--;
         }
      }


   /* do a final check */
   for (irrep=0; irrep<nirreps; irrep++) {
      if (used[irrep] > orbs_per_irrep[irrep]) {
         fprintf(stderr, "(reorder_qt): on final check, used more orbitals");
         fprintf(stderr, "   than were available (%d vs %d) for irrep %d\n",
            used[irrep], orbs_per_irrep[irrep], irrep);
         } 
      }

   free(used);  free(offset);
   free(docc);  free(socc);  free(frozen_docc);  free(frozen_uocc);
   free(uocc);
}

void reorder_qt_uhf	(	int *	docc,
		int *	socc,
		int *	frozen_docc,
		int *	frozen_uocc,
		int *	order_alpha,
		int *	order_beta,
		int *	orbspi,
		int	nirreps
	)

reorder_qt_uhf()

Generalization of reorder_qt() for UHF case

Parameters:

	docc	= doubly occupied orbitals per irrep
	socc	= singly occupied orbitals per irrep
	frozen_docc	= frozen occupied orbitals per irrep
	frozen_uocc	= frozen unoccupied orbitals per irrep
	order_alpha	= reordering array for alpha (Pitzer->QT order)
	order_beta	= reordering array for beta (Pitzer->QT order)
	nirreps	= number of irreducible representations

Definition at line 164 of file reorder_qt.cc.

References init_int_array().

{
  int p, nmo;
  int cnt_alpha, cnt_beta, irrep, point, point_alpha, point_beta, tmpi;
  int *offset, this_offset; 
  int *uocc;

  offset = init_int_array(nirreps);

  uocc = init_int_array(nirreps);

  /* construct the offset array */
  offset[0] = 0;
  for (irrep=1; irrep<nirreps; irrep++) {
    offset[irrep] = offset[irrep-1] + orbspi[irrep-1];
  }
   
  /* construct the uocc array */
  nmo = 0;
  for (irrep=0; irrep<nirreps; irrep++) {
    nmo += orbspi[irrep];
    tmpi = frozen_uocc[irrep] + docc[irrep] + socc[irrep];
    if (tmpi > orbspi[irrep]) {
      fprintf(stderr, "(reorder_qt_uhf): orbitals don't add up for irrep %d\n",
              irrep);
      return;
    }
    else
      uocc[irrep] = orbspi[irrep] - tmpi;
  }

  cnt_alpha = cnt_beta = 0;

  /* do the frozen core */
  for(irrep=0; irrep<nirreps; irrep++) { 
    this_offset = offset[irrep];
    for(p=0; p < frozen_docc[irrep]; p++) {
      order_alpha[this_offset+p] = cnt_alpha++;
      order_beta[this_offset+p] = cnt_beta++;
    }
  }

  /* alpha occupied orbitals */
  for(irrep=0; irrep<nirreps; irrep++) {
    this_offset = offset[irrep] + frozen_docc[irrep];
    for(p=0; p < docc[irrep] + socc[irrep] - frozen_docc[irrep]; p++) {
      order_alpha[this_offset + p] = cnt_alpha++;
    }
  }

  /* beta occupied orbitals */
  for(irrep=0; irrep<nirreps; irrep++) {
    this_offset = offset[irrep] + frozen_docc[irrep];
    for(p=0; p < docc[irrep] - frozen_docc[irrep]; p++) {
      order_beta[this_offset + p] = cnt_beta++;
    }
  }

  /* alpha unoccupied orbitals */
  for(irrep=0; irrep<nirreps; irrep++) {
    this_offset = offset[irrep] + docc[irrep] + socc[irrep];
    for(p=0; p < uocc[irrep]; p++) {
      order_alpha[this_offset + p] = cnt_alpha++;
    }
  }

  /* beta unoccupied orbitals */
  for(irrep=0; irrep<nirreps; irrep++) {
    this_offset = offset[irrep] + docc[irrep];
    for(p=0; p < uocc[irrep] + socc[irrep]; p++) {
      order_beta[this_offset + p] = cnt_beta++;
    }
  }

  /* do the frozen uocc */
  for (irrep=0; irrep<nirreps; irrep++) {
    this_offset = offset[irrep] + docc[irrep] + socc[irrep] + uocc[irrep];
    for(p=0; p < frozen_uocc[irrep]; p++) {
      order_alpha[this_offset + p] = cnt_alpha++;
      order_beta[this_offset + p] = cnt_beta++;
    }
  }

  /* do a final check */
  for (irrep=0; irrep<nirreps; irrep++) {
    if (cnt_alpha > nmo) {
      fprintf(stderr, "(reorder_qt_uhf): on final check, used more orbitals");
      fprintf(stderr, "   than were available (%d vs %d) for irrep %d\n",
              cnt_alpha, nmo, irrep);
    } 
    if (cnt_beta > nmo) {
      fprintf(stderr, "(reorder_qt_uhf): on final check, used more orbitals");
      fprintf(stderr, "   than were available (%d vs %d) for irrep %d\n",
              cnt_beta, nmo, irrep);
    } 
  }

  free(offset);
  free(uocc);
}

void schmidt	(	double **	A,
		int	rows,
		int	cols,
		FILE *	outfile
	)

SCHMIDT(): Gram-Schmidt orthogonalize a set of vectors

Assume we're orthogonalizing the ROWS, since in C a vector is usually a row more often than a column.

David Sherrill, Feb 1994

Parameters:

	A	= matrix to orthogonalize (matrix of doubles)
	rows	= rows of A
	cols	= columns of A

Returns: none

Definition at line 31 of file schmidt.cc.

References dot_arr(), and init_array().

{
   double *tmp;
   double normval, dotval;
   int i, j;
   register int I;

   /* initialize working array */
   tmp = init_array(cols);
  
   /* always take the first vector (normalized) as given */
   dot_arr(A[0], A[0], cols, &normval) ; /* normval = dot (A0 * A0) */
   normval = sqrt(normval) ;

   for (i=0; i<cols; i++) A[0][i] /= normval; 

   /* now, one at a time, get the new rows */
   for (i=1; i<rows; i++) {
      for (I=0; I<cols; I++) tmp[I] = A[i][I] ;
      for (j=0; j<i; j++) {
         dot_arr(A[i], A[j], cols, &dotval) ; 
         for (I=0; I<cols; I++) tmp[I] -= dotval * A[j][I];
         }
      dot_arr(tmp, tmp, cols, &normval);
      normval = sqrt(normval);
      /* fprintf(outfile,"\n norm[%d] = %20.15f\n",i, (1.0/normval));
      fflush(outfile); */
      for (I=0; I<cols; I++) A[i][I] = tmp[I] / normval; 
      }

   free(tmp);

}

int schmidt_add	(	double **	A,
		int	rows,
		int	cols,
		double *	v
	)

SCHMIDT_ADD(): Assume A is a orthogonal matrix. This function Gram-Schmidt orthogonalizes a new vector v and adds it to matrix A. A must contain a free row pointer for a new row. Don't add orthogonalized v' if norm(v') < NORM_TOL.

David Sherrill, Feb 1994

Parameters:

	A	= matrix to add new vector to
	rows	= current number of rows in A (A must have ptr for 'rows+1' row.)
	cols	= columns in A v = vector to add to A after it has been made orthogonal to rest of A

Returns: 1 if a vector is added to A, 0 otherwise

Definition at line 33 of file lib/libqt/schmidt_add.cc.

References dot_arr(), and init_array().

Referenced by david().

{
   double dotval, normval ;
   int i, I ;

   for (i=0; i<rows; i++) {
      dot_arr(A[i], v, cols, &dotval) ;
      for (I=0; I<cols; I++) v[I] -= dotval * A[i][I] ;
      }                 

   dot_arr(v, v, cols, &normval) ;
   normval = sqrt(normval) ;

   if (normval < NORM_TOL) 
      return(0) ;
   else {
      if (A[rows] == NULL) A[rows] = init_array(cols) ;
      for (I=0; I<cols; I++) A[rows][I] = v[I] / normval ;
      return(1) ;
      }
}

double secant	(	double(*)(double)	F,
		double	x0,
		double	x1,
		double	tolerance,
		int	maxiter,
		int	printflag
	)

secant(): Find the root of a function by the Secant Method. Iterations are limited to a maximum value. The algorithm stops when the relative difference between successive guesses is less than the specified tolerance. An initial TWO guesses for the root must be given, as well as the function itself.

Parameters:

	F	= pointer to function we want to examine (must return double)
	x0	= 1st guess for root
	x1	= 2nd guess for root
	tolerance	= how close successive guesses must get before convergence
	maxiter	= maximum number of iterations
	printflag	= whether or not to print results for each iteration (1 or 0)

Returns: the value of the root

Definition at line 171 of file rootfind.cc.

{
   int iter = 2 ;
   double x2=0.0 ;
   double fx1, fx0 ;

   if (printflag) {
      printf("Secant Method root finder\n") ;
      printf("%4s %12s\n", "iter", "x") ;
      printf("%4d %12.8f\n", 1, x0) ;
      printf("%4d %12.8f\n", 2, x1) ;
      }


   while (iter < maxiter) {
      fx1 = F(x1) ; 
      fx0 = F(x0) ;
      x2 = x1 - fx1 * (x1 - x0) / (fx1 - fx0) ;
      iter++ ;
      if (printflag) printf("%4d %12.8f\n", iter, x2) ;
      if (fabs((x2 - x1)/x1) < tolerance) return(x2) ;
      x0 = x1 ;
      x1 = x2 ;
      }

   printf("(secant): Failed to converge within %d iterations\n", iter) ;
   return(x2) ;
}

void slaterdetset_add	(	SlaterDetSet *	sdset,
		int	index,
		int	alphastring,
		int	betastring
	)

slaterdetset_add(): Add info for a particular Slater determinant to a SlaterDetSet.

Parameters:

	sdset	= pointer to SlaterDetSet to add to
	index	= location in the set to add to
	alphastring	= alpha string ID for the new determinant
	betastring	= beta string ID for the new determinant

Returns: none

Definition at line 358 of file slaterdset.cc.

References SlaterDet::alphastring, SlaterDetSet::alphastrings, SlaterDet::betastring, SlaterDetSet::betastrings, SlaterDetSet::dets, and SlaterDet::index.

{
  SlaterDet *det;
  StringSet *alphaset = sdset->alphastrings;
  StringSet *betaset = sdset->betastrings;

  if (index < sdset->size && index >= 0) {
    det = sdset->dets + index;
  }
  det->index = index;
  if (alphastring < alphaset->size && alphastring >= 0)
    det->alphastring = alphastring;
  if (betastring < betaset->size && betastring >= 0)
    det->betastring = betastring;
}

void slaterdetset_delete

(

SlaterDetSet *

sdset

)

slaterdetset_delete(): Delete a Slater Determinant Set.

Does not free the members alphastrings and betastrings. See also: slaterdetset_delete_full() which does this.

Parameters:

sdset

= pointer to SlaterDetSet to be de-allocated

Returns: none

Definition at line 306 of file slaterdset.cc.

References SlaterDetSet::alphastrings, SlaterDetSet::betastrings, SlaterDetSet::dets, and SlaterDetSet::size.

{
  sdset->size = 0;
  if (sdset->dets) {
    free(sdset->dets);
    sdset->dets = NULL;
  }
  sdset->alphastrings = NULL;
  sdset->betastrings = NULL;
}

void slaterdetset_delete_full

(

SlaterDetSet *

sdset

)

slaterdetset_delete_full(): De-allocate a Slater Determinant Set.

Frees memory including alpha and beta strings. See slaterdetset_delete() for a similar version which does not free the alpha and beta strings.

Parameters:

sdset

= pointer to SlaterDetSet to be de-allocated

Returns: none

Definition at line 329 of file slaterdset.cc.

References SlaterDetSet::alphastrings, SlaterDetSet::betastrings, SlaterDetSet::dets, SlaterDetSet::size, and stringset_delete().

Referenced by slaterdetvector_delete_full().

{
  sdset->size = 0;
  if (sdset->dets) {
    free(sdset->dets);
    sdset->dets = NULL;
  }
  if (sdset->alphastrings) {
    stringset_delete(sdset->alphastrings);
    sdset->alphastrings = NULL;
  }
  if (sdset->betastrings) {
    stringset_delete(sdset->betastrings);
    sdset->betastrings = NULL;
  }
}

void slaterdetset_init	(	SlaterDetSet *	sdset,
		int	size,
		StringSet *	alphastrings,
		StringSet *	betastrings
	)

slaterdetset_init(): Initialize a Slater Determinant Set

Parameters:

	sdset	= pointer to SlaterDetSet being initialized
	size	= number of SlaterDets to be held
	alphastrings	= pointer to StringSet of alpha strings
	betastrings	= pointer to StringSet of beta strings

Returns: none

Definition at line 285 of file slaterdset.cc.

References SlaterDetSet::alphastrings, SlaterDetSet::betastrings, SlaterDetSet::dets, and SlaterDetSet::size.

Referenced by slaterdetset_read().

{
  sdset->size = size;
  sdset->dets = (SlaterDet *) malloc(size*sizeof(SlaterDet));
  memset(sdset->dets,0,size*sizeof(SlaterDet));
  sdset->alphastrings = alphastrings;
  sdset->betastrings = betastrings;
}

void slaterdetset_read	(	ULI	unit,
		char *	prefix,
		SlaterDetSet **	slaterdetset
	)

slaterdetset_read(): Read a Slater Determinant Set

Parameters:

	unit	= file number of the PSIO file
	prefix	= prefix string to come before libpsio entry keys
	sdset	= pointer to SlaterDetSet to read into

Returns: none

Definition at line 437 of file slaterdset.cc.

References SlaterDetSet::dets, psio_read_entry(), SlaterDetSet::size, slaterdetset_init(), and stringset_read().

Referenced by slaterdetvector_read().

{
  int i, size;
  int need_to_init_psio = 0;
  int unit_opened = 1;
  char *size_key, *set_key;
  char *alphaprefix, *betaprefix;
  psio_address ptr;
  StringSet *alphastrings, *betastrings;
  SlaterDetSet *sdset = (SlaterDetSet *) malloc(sizeof(SlaterDetSet));

PSIO_INIT
PSIO_OPEN(unit,PSIO_OPEN_OLD)

  alphaprefix = (char *) malloc( strlen(prefix) + 
    strlen(SDSET_KEY_ALPHASTRINGS) + 2);
  sprintf(alphaprefix,"%s:%s",prefix,SDSET_KEY_ALPHASTRINGS);
  betaprefix = (char *) malloc( strlen(prefix) + 
    strlen(SDSET_KEY_BETASTRINGS) + 2);
  sprintf(betaprefix,"%s:%s",prefix,SDSET_KEY_BETASTRINGS);

  stringset_read( unit, alphaprefix, &alphastrings);
  stringset_read( unit, betaprefix, &betastrings);
  
  free(alphaprefix);
  free(betaprefix);

  size_key = (char *) malloc( strlen(prefix) + strlen(SDSET_KEY_SIZE) + 3);
  sprintf(size_key,":%s:%s",prefix,SDSET_KEY_SIZE);
  set_key = (char *) malloc( strlen(prefix) + strlen(SDSET_KEY_DETERMINANTS) 
    + 3);
  sprintf(set_key,":%s:%s",prefix,SDSET_KEY_DETERMINANTS);

  psio_read_entry( unit, size_key, (char *)&size, sizeof(int));
  slaterdetset_init(sdset,size,alphastrings,betastrings);
  psio_read_entry( unit, set_key, (char *)sdset->dets, 
    sdset->size*sizeof(SlaterDet));

PSIO_CLOSE(unit)
PSIO_DONE

  free(size_key);
  free(set_key);

  *slaterdetset = sdset;
}

void slaterdetset_read_vect	(	ULI	unit,
		char *	prefix,
		double *	coeffs,
		int	size,
		int	vectnum
	)

slaterdetset_read_vect(): Read in the coefficients for a single vector associated with a SlaterDetSet.

This function already assumes we've already called slaterdetset_read() to read in the string and determinant information. This is only going to read in the coefficients. This has been split out because we might want to read several roots for a given determinant setup. This does not actually depend on the presence of a SlaterDetVector object and is called a SlaterDetSet function.

Parameters:

	unit	= file number of the UNINITIALIZED PSIO file
	prefix	= prefix string to come before libpsio entry keys
	coeffs	= array to hold coefficients read
	size	= number of elements in coeffs array
	vectnum	= the vector number (for the PSIO key). Start from 0.

Returns: none

CDS 8/03

Definition at line 733 of file slaterdset.cc.

References psio_read_entry().

Referenced by slaterdetvector_read().

{
  int need_to_init_psio = 0;
  int unit_opened = 1;
  char *vector_key;

PSIO_INIT
PSIO_OPEN(unit,PSIO_OPEN_OLD)

  vector_key = (char *) malloc(strlen(prefix)+strlen(SDVECTOR_KEY_VECTOR)+5);
  sprintf(vector_key,":%s:%s%2d",prefix,SDVECTOR_KEY_VECTOR,vectnum);

  psio_read_entry(unit, vector_key, (char *)coeffs, size*sizeof(double));

PSIO_CLOSE(unit)
PSIO_DONE
  
  free(vector_key);

}

void slaterdetset_write	(	ULI	unit,
		char *	prefix,
		SlaterDetSet *	sdset
	)

slaterdetset_write(): Write a Slater Determinant Set to disk.

Parameters:

	unit	= file number to write to
	prefix	= prefix string to come before libpsio entry keys
	sdset	= pointer to SlaterDetSet to write

Returns: none

Definition at line 385 of file slaterdset.cc.

References SlaterDetSet::alphastrings, SlaterDetSet::betastrings, SlaterDetSet::dets, psio_write_entry(), SlaterDetSet::size, and stringset_write().

Referenced by slaterdetvector_write().

{
  int i;
  int need_to_init_psio = 0;
  int unit_opened = 1;
  char *size_key, *set_key;
  char *alphaprefix, *betaprefix;
  psio_address ptr;

PSIO_INIT
PSIO_OPEN(unit,PSIO_OPEN_OLD)

  alphaprefix = (char *) malloc( strlen(prefix) + 
    strlen(SDSET_KEY_ALPHASTRINGS) + 2);
  sprintf(alphaprefix,"%s:%s",prefix,SDSET_KEY_ALPHASTRINGS);
  betaprefix = (char *) malloc( strlen(prefix) + 
    strlen(SDSET_KEY_BETASTRINGS) + 2);
  sprintf(betaprefix,"%s:%s",prefix,SDSET_KEY_BETASTRINGS);

  stringset_write( unit, alphaprefix, sdset->alphastrings);
  stringset_write( unit, betaprefix, sdset->betastrings);
  
  free(alphaprefix);
  free(betaprefix);

  size_key = (char *) malloc( strlen(prefix) + strlen(SDSET_KEY_SIZE) + 3);
  sprintf(size_key,":%s:%s",prefix,SDSET_KEY_SIZE);
  set_key = (char *) malloc( strlen(prefix) + strlen(SDSET_KEY_DETERMINANTS) 
    + 3);
  sprintf(set_key,":%s:%s",prefix,SDSET_KEY_DETERMINANTS);

  psio_write_entry( unit, size_key, (char *)&sdset->size, sizeof(int));
  psio_write_entry( unit, set_key, (char *)sdset->dets, 
    sdset->size*sizeof(SlaterDet));

PSIO_CLOSE(unit)
PSIO_DONE

  free(size_key);
  free(set_key);
}

void slaterdetset_write_vect	(	ULI	unit,
		char *	prefix,
		double *	coeffs,
		int	size,
		int	vectnum
	)

slaterdetset_write_vect(): Write to disk the coefficients for a single vector associated with a set of Slater determinants.

This function already assumes we've already called slaterdetset_write() to write out the string and determinant information. This is only going to write out the coefficients. This has been split out because we might want to write several roots for a given determinant setup. This does not actually dpend on the presence of a SlaterDetVector object so it is called a SlaterDetSet function.

Parameters:

	unit	= file number of the UNINITIALIZED PSIO file
	prefix	= prefix string to come before libpsio entry keys
	coeffs	= array of coefficients to write
	size	= number of elements in coeffs array
	vectnum	= the vector number (to make a PSIO key). Start numbering from zero.

Returns: none

CDS 8/03

Definition at line 647 of file slaterdset.cc.

References psio_write_entry().

Referenced by slaterdetvector_write().

{
  int need_to_init_psio = 0;
  int unit_opened = 1;
  char *vector_key;

PSIO_INIT
PSIO_OPEN(unit,PSIO_OPEN_OLD)

  if (vectnum < 0 || vectnum > 99) {
    fprintf(stderr, "(slaterdetset_write_vect): vectnum out of bounds\n");
    abort();
  }

  vector_key = (char *) malloc(strlen(prefix)+strlen(SDVECTOR_KEY_VECTOR)+5);
  sprintf(vector_key,":%s:%s%2d",prefix,SDVECTOR_KEY_VECTOR,vectnum);

  psio_write_entry(unit, vector_key, (char *)coeffs, size*sizeof(double));

PSIO_CLOSE(unit)
PSIO_DONE

  free(vector_key);
}

void slaterdetvector_add	(	SlaterDetVector *	sdvector,
		int	index,
		double	coeff
	)

slaterdetvector_add(): Add a coefficient to a SlaterDetVector

Parameters:

	sdvector	= Pointer to SlaterDetVector to add to
	index	= location in vector for writing the coefficient
	coeff	= the coefficient to write to location index

Returns: none

Definition at line 561 of file slaterdset.cc.

References SlaterDetVector::coeffs.

{
  if (index < sdvector->size && index >= 0) {
    sdvector->coeffs[index] = coeff;
  }
}

void slaterdetvector_delete

(

SlaterDetVector *

sdvector

)

slaterdetvector_delete(): De-allocate a SlaterDetVector

Parameters:

sdvector

= pointer to SlaterDetVector to de-allocate

Note: does NOT fully free the associated SlaterDetSet. For that, see function slaterdetvector_delete_full()

Returns: none

Definition at line 514 of file slaterdset.cc.

References SlaterDetVector::coeffs, SlaterDetVector::sdset, and SlaterDetVector::size.

{
  sdvector->size = 0;
  sdvector->sdset = NULL;
  if (sdvector->coeffs) {
    free(sdvector->coeffs);
    sdvector->coeffs = NULL;
  }
}

void slaterdetvector_delete_full

(

SlaterDetVector *

sdvector

)

slaterdetvector_delete_full(): De-allocate a SlaterDetVector and its associated SlaterDetSet.

To keep the SlaterDetSet itself, use similar function slaterdetvector_delete().

Parameters:

sdvector

= pointer to SlaterDetVector to delete

Returns: none

Definition at line 537 of file slaterdset.cc.

References SlaterDetVector::coeffs, SlaterDetVector::sdset, SlaterDetVector::size, and slaterdetset_delete_full().

{
  sdvector->size = 0;
  if (sdvector->sdset) {
    slaterdetset_delete_full(sdvector->sdset);
    sdvector->sdset = NULL;
  }
  if (sdvector->coeffs) {
    free(sdvector->coeffs);
    sdvector->coeffs = NULL;
  }
}

void slaterdetvector_init	(	SlaterDetVector *	sdvector,
		SlaterDetSet *	sdset
	)

slaterdetvector_init(): Initialize a vector of coefficients corresponding to a Slater Determinant set

Parameters:

	sdvector	= pointer to SlaterDetVector to initialize (coeffs member will be allocated)
	sdset	= pointer to SlaterDetSet the vector is associated with

Returns: none

Definition at line 496 of file slaterdset.cc.

References SlaterDetVector::coeffs, init_array(), SlaterDetVector::sdset, SlaterDetSet::size, and SlaterDetVector::size.

Referenced by slaterdetvector_read().

{
  sdvector->size = sdset->size;
  sdvector->sdset = sdset;
  sdvector->coeffs = init_array(sdvector->size);
}

void slaterdetvector_read	(	ULI	unit,
		char *	prefix,
		SlaterDetVector **	sdvector
	)

slaterdetvector_read(): Read a SlaterDetVector from disk

Use this if we only need to read a single vector. Otherwise, call slaterdetset_read(); slaterdetset_read_vect(); to allow for multiple vectors per slaterdetset to be read from disk.

Parameters:

	unit	= file number to read from
	prefix	= prefix string to come before libpsio entry keys
	sdvector	= pointer to hold pointer to SlaterDetVector allocated by this function

Returns: none

Definition at line 690 of file slaterdset.cc.

References SlaterDetVector::coeffs, SlaterDetVector::size, slaterdetset_read(), slaterdetset_read_vect(), and slaterdetvector_init().

{
  int need_to_init_psio = 0;
  int unit_opened = 1;
  SlaterDetSet *sdset;
  SlaterDetVector *vector = (SlaterDetVector *) malloc(sizeof(SlaterDetVector));

PSIO_INIT
PSIO_OPEN(unit,PSIO_OPEN_OLD)

  slaterdetset_read(unit, prefix, &sdset);
  slaterdetvector_init(vector,sdset);
  slaterdetset_read_vect(unit, prefix, vector->coeffs, vector->size, 0);

PSIO_CLOSE(unit)
PSIO_DONE
  
  *sdvector = vector;
}

void slaterdetvector_set	(	SlaterDetVector *	sdvector,
		double *	coeffs
	)

slaterdetvector_set(): Set a SlaterDetVector's vector to a set of coefficients supplied by array coeffs

Parameters:

	sdvector	= pointer to SlaterDetVector for writing coefficients
	coeffs	= array of coefficients to write to sdvector

Definition at line 578 of file slaterdset.cc.

References SlaterDetVector::coeffs, and SlaterDetVector::size.

{
  int i;
  const int size = sdvector->size;
  double *v = sdvector->coeffs;
  if (v) {
    for(i=0; i<size; i++)
      v[i] = coeffs[i];
  }
}

void slaterdetvector_write	(	ULI	unit,
		char *	prefix,
		SlaterDetVector *	vector
	)

slaterdetvector_write(): Write a SlaterDetVector to disk.

This writes a vector in the space of Slater determinants, along with the set of determinants itself, to a PSIO file.

Use this if we only need to write a single vector. Otherwise, call slaterdetset_write(); slaterdetset_write_vect(); to allow for multiple vectors per slaterdetset to be written to disk.

Parameters:

	unit	= file number of the UNINITIALIZED PSIO file
	prefix	= prefix string to come before libpsio entry keys
	vector	= SlaterDetVector to write to disk

Returns: none

Definition at line 607 of file slaterdset.cc.

References SlaterDetVector::coeffs, SlaterDetVector::sdset, SlaterDetVector::size, slaterdetset_write(), and slaterdetset_write_vect().

{
  int need_to_init_psio = 0;
  int unit_opened = 1;

PSIO_INIT
PSIO_OPEN(unit,PSIO_OPEN_OLD)

  slaterdetset_write(unit, prefix, vector->sdset);
  slaterdetset_write_vect(unit, prefix, vector->coeffs, vector->size, 0);

PSIO_CLOSE(unit)
PSIO_DONE

}

void solve_2x2_pep	(	double **	H,
		double	S,
		double *	evals,
		double **	evecs
	)

solve_2x2_pep(): Solve a 2x2 pseudo-eigenvalue problem of the form [ H11 - E H12 - E*S ] [c1] [ H12 - E*S H22 - E ] [c2] = 0

Parameters:

	H	= matrix to get eigenvalues of
	S	= overlap between states 1 & 2
	evals	= pointer to array to hold 2 eigenvalues
	evecs	= matrix to hold 2 eigenvectors

Returns: none

Definition at line 28 of file solve_pep.cc.

{
   int i ;
   double a, b, c, p, q, x ;
   double norm, tval, tval2; 

   /* put in quadratic form */
   a = 1.0 - S * S ;
   b = 2.0 * S * H[0][1] - H[0][0] - H[1][1] ;
   c = H[0][0] * H[1][1] - H[0][1] * H[0][1] ; 
  
   /* solve the quadratic equation for E0 and E1 */

   tval = b*b ;
   tval -= 4.0 * a * c ;
   tval2 = sqrt(tval) ;
   if (tval < 0.0) {
      fprintf(stderr, "(solve_2x2_pep): radical less than 0\n") ;
      return ;
      }
   else if (fabs(a) < A_MIN) {
      printf("min a reached\n") ;
      evals[0] = evals[1] = H[1][1] ;
      }
   else { 
      evals[0] = -b / (2.0 * a) ;
      evals[0] -= tval2 / (2.0 * a) ; 
      evals[1] = -b / (2.0 * a) ;
      evals[1] += tval2 / (2.0 * a) ;
      }
 
   /* Make sure evals[0] < evals[1] */
   if (evals[1] < evals[0]) { 
      tval = evals[0] ;  
      evals[0] = evals[1] ;
      evals[1] = tval ;
      } 

   /* Make sure evals[0] < H[1][1] */
   if (evals[0] > H[1][1]) {
      printf("Warning: using H11 as eigenvalues\n") ;
      evals[0] = evals[1] = H[1][1] ;
      printf("Got evals[0] = H[1][1] = %12.7f\n", evals[0]) ;
      }

   /* get the eigenvectors */
   for (i=0; i<2; i++) {
      p = H[0][0] - evals[i] ;
      q = H[0][1] - S * evals[i] ;
      x = -p/q ; 
      norm = 1.0 + x * x ;
      norm = sqrt(norm) ;
      evecs[i][0] = 1.0 / norm ;
      evecs[i][1] = x / norm ;
      }

   /* test 
   for (i=0; i<2; i++) {
      p = H[i][0] * evecs[0][0] + H[i][1] * evecs[0][1] ;
      if (i==0) q = evecs[0][0] + S * evecs[0][1] ;
      else q = S * evecs[0][0] + evecs[0][1] ;
      q *= evals[0] ;
      printf("2x2 check %d: LHS = %12.6f  RHS = %12.6f\n", i, p, q) ;
      }
    */

}

void sort_vector	(	double *	A,
		int	n
	)

sort_vector(): Sort the elements of a vector into increasing order

Parameters:

	A	= array to sort
	n	= length of array

Returns: none

Definition at line 55 of file lib/libqt/sort.cc.

{
  int i, j, k;
  double val;

  for(i=0; i < n-1; i++) {
    val = A[k=i];

    for(j=i+1; j < n; j++)
      if(A[j] <= val) val = A[k=j];

    if(k != i) {
      A[k] = A[i];
      A[i] = val;
    }
  }
}

void stringset_add	(	StringSet *	sset,
		int	index,
		unsigned char *	Occ
	)

stringset_add(): Add a string (in Pitzer order, given by Occ) to the StringSet, writing to position index.

Parameters:

	sset	= StringSet to add to
	index	= location in StringSet to add to
	Occ	= orbital occupations (Pitzer order) of the string to add

Returns: none

Definition at line 98 of file slaterdset.cc.

References String::index, StringSet::nelec, StringSet::nfzc, String::occ, and StringSet::strings.

{
  int i;
  int nact = sset->nelec - sset->nfzc;
  String *s;

  if (index < sset->size && index >= 0) {
    s = sset->strings + index;
  }
  s->index = index;
  s->occ = (short int*) malloc(nact*sizeof(short int));
  for(i=0;i<nact;i++)
    s->occ[i] = Occ[i];
}

void stringset_delete

(

StringSet *

sset

)

stringset_delete(): Delete a StringSet

Parameters:

sset

= pointer to StringSet to delete

Returns: none

Definition at line 77 of file slaterdset.cc.

References StringSet::fzc_occ, StringSet::nelec, StringSet::nfzc, StringSet::size, and StringSet::strings.

Referenced by slaterdetset_delete_full().

{
  if (sset->nfzc > 0) free(sset->fzc_occ);
  sset->size = 0;
  sset->nelec = 0;
  sset->nfzc = 0;
  if (sset->strings) free(sset->strings);
  sset->strings = NULL;
}

void stringset_init	(	StringSet *	sset,
		int	size,
		int	nelec,
		int	nfzc,
		short int *	frozen_occ
	)

stringset_init(): Initialize a set of alpha/beta strings

Parameters:

	sset	= pointer to StringSet (contains occs in Pitzer order)
	size	= number of strings
	nelec	= number of electrons
	nfzc	= number of frozen core orbitals
	frozen_occ	= array of frozen occupied orbitals (Pitzer numbering!)

Returns: none

Definition at line 49 of file slaterdset.cc.

References StringSet::fzc_occ, StringSet::nelec, StringSet::nfzc, StringSet::size, and StringSet::strings.

Referenced by stringset_read().

{
  int i;

  sset->size = size;
  sset->nelec = nelec;
  sset->nfzc = nfzc;
  sset->strings = (String *) malloc(size*sizeof(String));
  memset(sset->strings,0,size*sizeof(String));
  if (nfzc > 0) {
    sset->fzc_occ = (short int *) malloc(nfzc * sizeof(short int));
    for (i=0; i<nfzc; i++) {
      sset->fzc_occ[i] = frozen_occ[i];  
    }
  }
}

void stringset_read	(	ULI	unit,
		char *	prefix,
		StringSet **	stringset
	)

stringset_read(): Read a StringSet from disk

Parameters:

	unit	= file number to read from
	prefix	= prefix string to come before libpsio entry keys
	stringset	= double pointer to StringSet allocated by this function

Returns: none

Definition at line 213 of file slaterdset.cc.

References String::index, String::occ, psio_read(), psio_read_entry(), PSIO_ZERO, StringSet::strings, and stringset_init().

Referenced by slaterdetset_read().

{
  int i, size, nelec, nfzc, nact;
  int need_to_init_psio = 0;
  int unit_opened = 1;
  char *size_key, *nelec_key, *nfzc_key, *fzc_occ_key, *strings_key;
  short int *fzc_occ;
  psio_address ptr;
  StringSet *sset = (StringSet *) malloc(sizeof(StringSet));

PSIO_INIT
PSIO_OPEN(unit,PSIO_OPEN_OLD)

  size_key = (char *) malloc( strlen(prefix) + strlen(STRINGSET_KEY_SIZE) + 3);
  sprintf(size_key,":%s:%s",prefix,STRINGSET_KEY_SIZE);
  nelec_key = (char *) malloc( strlen(prefix) + strlen(STRINGSET_KEY_NELEC)+ 3);
  sprintf(nelec_key,":%s:%s",prefix,STRINGSET_KEY_NELEC);
  nfzc_key = (char *) malloc( strlen(prefix) + strlen(STRINGSET_KEY_NFZC) + 3);
  sprintf(nfzc_key,":%s:%s",prefix,STRINGSET_KEY_NFZC);
  fzc_occ_key = (char *) malloc(strlen(prefix) +
    strlen(STRINGSET_KEY_FZC_OCC) + 3);
  sprintf(fzc_occ_key,":%s:%s",prefix,STRINGSET_KEY_FZC_OCC);
  strings_key = (char *) malloc( strlen(prefix) + strlen(STRINGSET_KEY_STRINGS)
    + 3);
  sprintf(strings_key,":%s:%s",prefix,STRINGSET_KEY_STRINGS);

  psio_read_entry( unit, size_key, (char *)&size, sizeof(int));
  psio_read_entry( unit, nelec_key, (char *)&nelec, sizeof(int));
  psio_read_entry( unit, nfzc_key, (char *)&nfzc, sizeof(int));
  if (nfzc > 0) {
    fzc_occ = (short int *) malloc(nfzc*sizeof(short int));
    psio_read_entry( unit, fzc_occ_key, (char *)fzc_occ, 
      nfzc*sizeof(short int));
  }
  else fzc_occ = NULL;

  stringset_init(sset, size, nelec, nfzc, fzc_occ);

  nact = nelec - nfzc;
  ptr = PSIO_ZERO;
  for(i=0; i<size; i++) {
    psio_read( unit, strings_key, (char *) &(sset->strings[i].index), 
      sizeof(int), ptr, &ptr);
    sset->strings[i].occ = (short int*) malloc(nact*sizeof(short int));
    psio_read( unit, strings_key, (char *) sset->strings[i].occ, 
      nact*sizeof(short int), ptr, &ptr);
  }

PSIO_CLOSE(unit)
PSIO_DONE

  free(size_key);
  free(nelec_key);
  free(nfzc_key);
  free(fzc_occ_key);
  free(strings_key);
  if (nfzc > 0) free(fzc_occ);
  *stringset = sset;
}

void stringset_reindex	(	StringSet *	sset,
		short int *	mo_map
	)

stringset_reindex(): Remap orbital occupations from one ordering to another.

Parameters:

	sset	= pointer to StringSet
	mo_map	= mapping array from original orbital order to new order

Returns: none

Definition at line 123 of file slaterdset.cc.

References StringSet::fzc_occ, StringSet::nelec, StringSet::nfzc, StringSet::size, and StringSet::strings.

{
  int s, mo, core;
  short int* occ;
  int nstrings = sset->size;
  int nact = sset->nelec - sset->nfzc;

  for (core=0; core<sset->nfzc; core++) {
    sset->fzc_occ[core] = mo_map[sset->fzc_occ[core]];
  }

  for(s=0; s<nstrings; s++) {
    occ = (sset->strings + s)->occ;
    for(mo=0; mo<nact; mo++)
      occ[mo] = mo_map[occ[mo]];
  }
}

void stringset_write	(	ULI	unit,
		char *	prefix,
		StringSet *	sset
	)

stringset_write(): Write a stringset to a PSI file

Parameters:

	unit	= file number to write to
	prefix	= prefix string to come before libpsio entry keys
	sset	= pointer to StringSet to write

Returns: none

Definition at line 151 of file slaterdset.cc.

References StringSet::fzc_occ, String::index, StringSet::nelec, StringSet::nfzc, String::occ, psio_write(), psio_write_entry(), PSIO_ZERO, StringSet::size, and StringSet::strings.

Referenced by slaterdetset_write().

{
  int i, size, nact;
  int need_to_init_psio = 0;
  int unit_opened = 1;
  char *size_key, *nelec_key, *nfzc_key, *strings_key, *fzc_occ_key;
  psio_address ptr;

PSIO_INIT
PSIO_OPEN(unit,PSIO_OPEN_OLD)

  size_key = (char *) malloc(strlen(prefix) + strlen(STRINGSET_KEY_SIZE) + 3);
  sprintf(size_key,":%s:%s",prefix,STRINGSET_KEY_SIZE);
  nelec_key = (char *) malloc(strlen(prefix) + strlen(STRINGSET_KEY_NELEC) + 3);
  sprintf(nelec_key,":%s:%s",prefix,STRINGSET_KEY_NELEC);
  nfzc_key = (char *) malloc(strlen(prefix) + strlen(STRINGSET_KEY_NFZC) + 3);
  sprintf(nfzc_key,":%s:%s",prefix,STRINGSET_KEY_NFZC);
  fzc_occ_key = (char *) malloc(strlen(prefix) +
    strlen(STRINGSET_KEY_FZC_OCC) + 3);
  sprintf(fzc_occ_key,":%s:%s",prefix,STRINGSET_KEY_FZC_OCC);
  strings_key = (char *) malloc(strlen(prefix) + strlen(STRINGSET_KEY_STRINGS) + 3);
  sprintf(strings_key,":%s:%s",prefix,STRINGSET_KEY_STRINGS);

  psio_write_entry( unit, size_key, (char *)&sset->size, sizeof(int));
  psio_write_entry( unit, nelec_key, (char *)&sset->nelec, sizeof(int));
  psio_write_entry( unit, nfzc_key, (char *)&sset->nfzc, sizeof(int));
  if (sset->nfzc) {
    psio_write_entry( unit, fzc_occ_key, (char *)sset->fzc_occ, 
      sset->nfzc*sizeof(short int));
  }

  ptr = PSIO_ZERO;
  size = sset->size;
  nact = sset->nelec - sset->nfzc;
  for(i=0; i<size; i++) {
    psio_write( unit, strings_key, (char *) &(sset->strings[i].index), 
      sizeof(int), ptr, &ptr);
    psio_write( unit, strings_key, (char *) sset->strings[i].occ, 
      nact*sizeof(short int), ptr, &ptr);
  }

PSIO_CLOSE(unit)
PSIO_DONE

  free(size_key);
  free(nelec_key);
  free(nfzc_key);
  free(strings_key);
  free(fzc_occ_key);
}

void timer_done

(

void

)

timer_done(): Close down all timers and write results to timer.dat

Definition at line 91 of file timer.cc.

References ffile().

Referenced by main().

{
  FILE *timer_out;
  char *host;
  extern time_t timer_start, timer_end;
  struct timer *this_timer, *next_timer;
  extern struct timer *global_timer;

  timer_end = time(NULL);

  host = (char *) malloc(40 * sizeof(char));
  gethostname(host, 40);

  /* Dump the timing data to timer.dat and free the timers */
  ffile(&timer_out, "timer.dat", 1);
  fprintf(timer_out, "\n");
  fprintf(timer_out, "Host: %s\n", host);
  fprintf(timer_out, "\n");
  fprintf(timer_out, "Timers On : %s", ctime(&timer_start));
  fprintf(timer_out, "Timers Off: %s", ctime(&timer_end));
  fprintf(timer_out, "\nWall Time:  %10.2f seconds\n\n",
          (double) timer_end - timer_start);

  this_timer = global_timer;
  while(this_timer != NULL) {
      if(this_timer->calls > 1) 
          fprintf(timer_out, "%-12s: %10.2fu %10.2fs %10.2fw %6d calls\n",
                  this_timer->key, this_timer->utime, this_timer->stime,
                  this_timer->wtime, this_timer->calls);
      else if(this_timer->calls == 1)
          fprintf(timer_out, "%-12s: %10.2fu %10.2fs %10.2fw %6d call\n",
                  this_timer->key, this_timer->utime, this_timer->stime,
                  this_timer->wtime, this_timer->calls);
      next_timer = this_timer->next;
      free(this_timer);
      this_timer = next_timer;
    }

  fprintf(timer_out, 
          "\n***********************************************************\n");
  fclose(timer_out);

  free(host);
}

void timer_init

(

void

)

timer_init(): Initialize the linked list of timers

Definition at line 76 of file timer.cc.

Referenced by main().

{
  extern struct timer *global_timer;
  extern time_t timer_start;

  timer_start = time(NULL);

  global_timer = NULL;
}

void timer_off

(

char *

key

)

timer_off(): Turn off the timer with the name given as an argument. Can be turned on and off, time will accumulate while on.

Parameters:

key

= Name of timer

Definition at line 229 of file timer.cc.

References timer_scan().

Referenced by dpd_buf4_scm(), and main().

{
  struct tms ontime, offtime;
  struct timer *this_timer;
  time_t wall_stop;

  this_timer = timer_scan(key);

  if(this_timer == NULL) {
      fprintf(stderr, "Bad timer key: %s\n", key);
      exit(PSI_RETURN_FAILURE);
    }

  if(this_timer->status == TIMER_OFF) {
     fprintf(stderr, "Timer %s is already off.\n", this_timer->key);
     exit(PSI_RETURN_FAILURE);
    }

  ontime = this_timer->ontime;

  times(&offtime);

  this_timer->utime += ((double) (offtime.tms_utime-ontime.tms_utime))/HZ;
  this_timer->stime += ((double) (offtime.tms_stime-ontime.tms_stime))/HZ;

  wall_stop = time(NULL);
  this_timer->wtime += ((double) (wall_stop - this_timer->wall_start));

  this_timer->status = TIMER_OFF;
}

void timer_on

(

char *

key

)

timer_on(): Turn on the timer with the name given as an argument. Can be turned on and off, time will accumulate while on.

Parameters:

key

= Name of timer

Definition at line 189 of file timer.cc.

References timer_scan().

Referenced by dpd_buf4_scm(), and main().

{
  struct timer *this_timer;

  this_timer = timer_scan(key);

  if(this_timer == NULL) { /* New timer */
      this_timer = (struct timer *) malloc(sizeof(struct timer));
      strcpy(this_timer->key,key);
      this_timer->calls = 0;
      this_timer->utime = 0;
      this_timer->stime = 0;
      this_timer->wtime = 0;
      this_timer->next = NULL;
      this_timer->last = timer_last();
      if(this_timer->last != NULL) this_timer->last->next = this_timer;
      else global_timer = this_timer;
    }
  else {
    if((this_timer->status == TIMER_ON) && (this_timer->calls)) {
        fprintf(stderr, "Timer %s is already on.\n", key);
        exit(PSI_RETURN_FAILURE);
    }
  }

  this_timer->status = TIMER_ON;
  this_timer->calls++;
  
  times(&(this_timer->ontime));
  this_timer->wall_start = time(NULL);
}

struct timer* timer_scan

(

char *

key

)

[read]

timer_scan(): Return a timer structure whose name matches that given by supplied string

Parameters:

key

= name of timer to search for

Returns: the timer structure with the given name, else NULL

Definition at line 145 of file timer.cc.

Referenced by timer_off(), and timer_on().

{
  extern struct timer *global_timer;
  struct timer *this_timer;

  this_timer = global_timer;
  
  while(this_timer != NULL) {
      if(!strcmp(this_timer->key,key)) return(this_timer);
      this_timer = this_timer->next;
    }

  return(this_timer);
}

---