/* Copyright (C) 1989, California Institute of Technology */ /* All rights reserved. U. S. Government Sponsorship under */ /* NASA Contract NAS7-918 is acknowledged. */ /* Herbert M. Pickett, 20 March 1989 */ /* Revised version in c, 22 March 1999 */ /* 24 March 1999: Fixed sin scaling in dircos */ /* 25 March 1999: Fixed diag = 2 error in call to ordham*/ /* 27 March 1999: Changed diag = 3 to take care of l-doubled states */ /* 30 March 1999: Fixed sin scaling bug */ /* 21 April 1999: Fixed bug at end of hamx for KROLL case */ /* 22 April 1999: Fixed bug in Itot option */ /* 25 May 1999: fix bug when largest number of spins is not default */ /* 23 July 1999: fix bug in SETOPT that returned incorrect IQNFMT */ /* 7 Sept. 1999: fix bug in GETQQ for that filled ivs with short */ /* 29 Dec. 1999: improve ORDHAM */ /* 12 April 2000: negate Hamiltonian for kappa > 0 on first pass */ /* 12 April 2000: separate ORDHAM and FIXHAM */ /* 12 April 2000: subtract constant from energy for more precision*/ /* 1 August 2000: special code for unit matrix */ /* 12 Sept. 2000: fix up phase for cross-multiplet interaction */ /* 9 Oct. 2000: change dipole type 10 to double commutator with N^2 */ /* 29 Nov. 2000: redefine idpars and spar to put more calc. in initial */ /* 5 Feb. 2001: define augmented Euler and Fourier series */ /* 7 July 2001: fix bug for l doublet pairs in setopt */ /* 26 Sept. 2001: fix bug in idpari for oblate 6x and 2x interactions */ /* 12 March 2002: improve selections for off-diagonal l in dircos */ /* 27 June 2002: increase vibration limit to 359 */ /* 29 July 2002: fix bug in fixham (DIAG = 2,3)*/ /* 18 Aug. 2003: code cleanup, @comment@ is for splint */ /* 24 Nov. 2003: implement code for higher symmetry IAX options */ /* 3 Feb 2004: implement code for higher symmetry spins */ /* 19 May 2004: implement code for bcd storage of idpar and idip */ /* 20 July 2004: fix bug in code for bcd storage of idip */ /* 19 Jan. 2005: extended Euler operator */ /* 28 May 2005: fix bug in getmask for operator between l-doubled state */ /* 27 Sept. 2005: fix bug for 6-fold symmetry B-weights in getqn & setopt */ /* 17 Oct. 2005: fix bug in getmask for operator within l-doubled state */ /* 12 Dec. 2005: fix bug Itot basis and include mixed basis */ #include #include #include #include #include #include "calpgm.h" #include "spinit.h" #define ODD(i) (((int)(i) & 1) != 0) #define EVEN(i) (((int)(i) & 1) == 0) #define ODD2(i) (((int)(i) & 2) != 0) #define MOD(i,n) ((int)(i) % (int)(n)) #define C0 '\0' #define MAXINT 12 /* maximum number of dipole types */ #define MAXVIB 999 #define MAXSPIN 9 #define MAXN_DIRCOS 359 /* MAX K in DIRCOS */ #define MAXVDEC 3 #define NDECPAR 5 /* number of digit pairs for idpar not including vib */ #define NSBCD (2 * (NDECPAR + MAXVDEC + 1)) #define MAXII 20 /* I <= 10 */ #define MAXNS 20 /* 2 * NSPINS + 2 */ #define MNSQ 2 /* bit for operator which is different from previous by only N*(N+1) */ #define MCOS_OK 4 /* bit for operator which uses cosine from previous */ #define MNOUNIT 8 /* bit for operator which is not a unit matrix */ #define MMASK (0x0fe0) /* clear bits 0-4 */ #define MIDEN 16 /* bit for Identity operator under alpha */ #define MODD 0x040 /* bit for imaginary operator */ #define MSYM2 0x080 /* bit for I_alpha - I_-alpha */ #define MLZ 0x100 /* bit for Lz operator */ #define MINOQ 4 /* bit for NO Delta N = 0 (commutator) */ #define MELEC 2 /* bit for electric dipole (must be 2)*/ #define NDXCOM 8 /* XCOM size, XCOM enum: */ #define XDIM 0 /* 0: DIMENSION OF SUB-BLOCK */ #define XSYM 1 /* 1: SYMMETRY CODE FOR BLOCK (0=A,1=Bx,2=By,3=Bz) */ #define XNVAL 2 /* 2: 'N' QUANTUM NUMBER */ #define XKBGN 3 /* 3: BEGINNING K */ #define XLVAL 4 /* 4: L VALUE */ #define XVIB 5 /* 5: VIBRATIONAL QUANTA */ #define XISYM 6 /* 6: SPIN SYMMETRY */ #define XIQN 7 /* 6: SPIN SYMMETRY QUANTUM NUMBER */ /* Common Declarations */ typedef struct { /* vibrational information */ /*@dependent@*/ short *spt; int nspstat; short knmin[4], knmax, wt[5], ewt[2], lvqn, lvupper, gsym, nqn; } SVIB; static SVIB vinfo1; /*@owned@*/ static SVIB *vinfo = &vinfo1; typedef struct str_spar { /* local parameter information */ /*@null@*/ /*@owned@*/ struct str_spar *next; double zfac; int ip; unsigned int ipsym; short flags; signed char ksq, fc, msi1, msi2, mldel; unsigned char njq, mln, mld, mkdel, mins, msznz, euler, alpha; } SPAR; typedef /*@null@*/ /*@owned@*/ SPAR *PSPAR; static PSPAR spar_head[MAXVIB]; static short sptzero[] = { 1, 0 }; /* default for no spins */ typedef struct struct_ssp { /*@null@*/ /*@owned@*/ /*@reldef@*/ struct struct_ssp *next; /*@notnull@*//*@owned@*/ short *sspt; int ssize, nitot; } SSP; static SSP ssp_head = { NULL, sptzero, 1 , 0}; static struct { /* save data for getqn */ int cblk, cnblk, csblk, cff, cwt[5]; } *cgetq, cgetv[2]; typedef struct { /* local dipole information */ double fac; short flg; signed char kd, ld, ldel, fc; } SDIP; /*@owned@*/ static SDIP *dipinfo; static double zero = 0.; static double zwk = 0.; static double spfac[MAXII]; static double spfac2[MAXII]; static int zmoldv = 0; static int zblkptr = 0; static int zivs = 0; static int zipder = 0; static int revsym[] = { 0, 3, 2, 1}; static int isoddk[] = { 0, 1, 1, 0}; static int is_esym[MAXITOT]; static int lscom[MAXNS], iscom[MAXNS], jscom[MAXNS], ismap[MAXNS]; static int ixcom[NDXCOM], jxcom[NDXCOM]; static int itptr, itsym, nspin, nsqmax, ndmx, ndmax, nddip; static short szero = 0; static short zidx = 0; static short zjdx = 0; static short zibkptr = 0; static short zikmin = 0; static struct { int mxspin, idiag, nvib, nfit, nbkpj, ixz, nqnn, nqn, maxqn, vibfac; int parinit, maxblk, nitot, vibdec, esymdec, msmask, nqn0, iqfmt0; unsigned int msshft; BOOL lsym, esym, oblate, vibfmt; } glob; static char sbcd[NSBCD]; /* pointers to dynamically allocated arrays */ /*@owned@*/ static int *moldv; /*@owned@*/ static int *blkptr; /*@owned@*/ static int *ipder, *ivs; /*@owned@*/ static double *wk; /*@owned@*/ static short *idx, *jdx, *ibkptr, *ikmin; int dclr(const int n1, const int n2, double *vec, const int ix); int specop(const int neuler, BOOL * newblk, int *nsqj, int *ikq, const int ksi, const int ksj, const int ni, const int nj, const int ncos, double *wk, const short *ix, const short *jx, const double par); int specfc(const int ifc, const int iv, const int jv, const int kdel, const int ksi, const int ksj, const int ncos, double *wk, const short *ix, const short *jx); int sznzfix(const int sznz, const int ni, const int nj, int *ixcom, int *jxcom, int *iscom, int *jscom); int sznzop(const int ni, const int nj, const int ksi, const int ksj, const int *iscom, const int *jscom, const int ncos, double *wk, const short *ix, const short *jx); unsigned int blksym(const int *ixcom, const int *jxcom); int ordham(const int nn, double *egy, const short *isblk, short *iswap); int fixham(const int ndm, const int nn, double *t, double *egy, double *p, const short *iswap); BOOL kroll(const int nsizd, double *t, const int nsblk, const short *sbkptr, const short *kmin); int getqs(const int mvs, const int iff, const int nsiz, const int kbgn, /*@out@*/ int *ixcom, /*@out@*/ int *iscom, /*@out@*/ int *iv); int idpars(SPAR * pspar, /*@out@*/ int *ksq, /*@out@*/ int *itp, /*@out@*/ int *l, /*@out@*/ int *ld, /*@out@*/ int *kdel, /*@out@*/ int *ins, /*@out@*/ int *si1, /*@out@*/ int *si2, /*@out@*/ int *sznz, /*@out@*/ int *ifc, /*@out@*/ int *alpha, /*@out@*/ int *ldel); int getll(const int llf, const int ld, const int kd, const int lt, const int si1, const int si2, int *lscom, const int *iscom, const int *jscom); int getmask(const int *xbra, const int *xket, const int kd, const int ldel, const int loff, const int alpha); double rmatrx(const int ld, const int lv, const int *ixcom, const int *jxcom); int symnsq(const int inq, const int ins, const int *iscom, const int *jscom, double *z); int symksq(const int ikq, const int ksi, const int ksj, const int n, double *wk, short *ix, short *jx); int dpmake(const int nsize, double *dp, const double *t, const int n, const double *wk, const short *ix, const short *jx, const BOOL isunit); int pasort(FILE * lu, const int npar, bcd_t *idpar, const double *par); int idpari(bcd_t *idval, int itp, /*@out@*/ SPAR * pspar); int checksp(const BOOL first, int si1, int si2, const short *iiv1, const short *iiv2, double *zfac); int setwt(SVIB * pvinfo, const int ivib, const int iax, const int iwtpl, const int iwtmn, double vsym); int getwt(SVIB * pvinfo, const int isym, const int iispin, /*@out@*/ int *ivwt); int setgsym(const int gsym); int getsp(const bcd_t *ispnx, SVIB *pvinfo); void setsp(void); int getqq(const int iblk, /*@out@*/ int *f, /*@out@*/ int *iwtb, /*@out@*/ short *sbkptr, /*@out@*/ short *kmin, /*@out@*/ int *vs); int dircos(const int *xbra, const int *xket, const int ld, const int kd, const int ncmax, /*@out@*/ double *direl, /*@out@*/ short *ibra, /*@out@*/ short *iket, const int ifup, const int loff, int mask, /*@out@*/ BOOL *isunit); int ffcal(const int nff, const int kff, /*@out@*/ double *ff); /*****************************************************************************/ int hamx(iblk, nsize, npar, idpar, par, egy, t, dedp, pmix, ifdump) const int iblk, nsize, npar; const BOOL ifdump; const bcd_t *idpar; const double *par; double *egy, *t, *dedp, *pmix; { /* .. PACKAGE FOR 98 INTERACTING VIBRATIONAL STATES WITH MULTI-SPIN */ /* calculate energies and derivatives for block IBLK */ /* on entry: */ /* IBLK= block number */ /* NSIZE= block size and dimension */ /* IDPAR= list of parameter identifiers ( element 0 is length ) */ /* PAR = list of parameter values */ /* NPARD= number of parameters with positive IDPAR */ /* on return: */ /* EGY = energy list */ /* T = eigenvector matrix */ /* DEDP = derivative of energy with respect to parameter */ static double sqj[10]; static int nold = -1; SPAR *spar_now; short *itau; double *pt, *dbar, pbar, sqnn, zpar, dn, ele, egy0; long ndm; int ipar, ierr, ncos, ispn, jspn, nsqj, iiwt[5], ndmd, alpha, mkd; int i, ii, lt, n, ibase, jbase, kbgni, kbgnj, nsblk, idflags; int kd, ld, nd, ni, nj, ivbase, jvbase, ivmin, ifc, iz, jz, npair, ldel; int si1, si2, iff, ijd, ikq, ins, neuler, sznz, ixx, jxx, kl, isgn; unsigned int ivsym, lastvv, ivcmp; BOOL oldpar, isegy0, roll, parskp, isunit, isneg, first; BOOL newblk, firstpar; if (ndmx <= 0) { puts("working vectors not allocated"); exit(EXIT_FAILURE); } dbar = &wk[ndmx]; ndm = nsize; ndmd = nsize + 1; isneg = (glob.oblate && !ifdump); cgetv[0].cblk = 0; cgetv[1].cblk = 0; /* get F and sub-block structure */ itau = &ikmin[glob.maxblk]; nsblk = getqq(iblk, &iff, iiwt, ibkptr, ikmin, ivs); /* zero hamiltonian matrix */ dclr(nsize, nsize, t, 1); /* zero derivative matrix */ n = glob.nfit; dclr(nsize, n, dedp, 1); /* set up hamiltonian */ oldpar = parskp = roll = isunit = firstpar = FALSE; ncos = sznz = 0; isgn = 1; pbar = zpar = 0.; first = TRUE; do { /* loop on wang blocks */ lastvv = 0; newblk = TRUE; egy0 = 0.; for (ixx = 0; ixx < nsblk; ++ixx) { ibase = ibkptr[ixx]; n = ibkptr[ixx + 1] - ibase; kbgni = ikmin[ixx]; ispn = ivs[ixx]; getqs(ispn, iff, n, kbgni, ixcom, iscom, &ivbase); ni = ixcom[XNVAL]; /* set up itau */ ld = ixcom[XLVAL]; if (ld == 0) ld = (ni + ixcom[XSYM]) & 1; else ld = (ld > 0)? 1: 0; kd = kbgni + kbgni - ld; kd = kd + kd; if (ivbase != ixcom[XVIB]) ++kd; /* upper l-doubled state */ itau[ixx] = (short) kd; if (ni != nold) { /* set up powers of N(N+1) */ nold = ni; dn = (double) ni; sqnn = dn * (ni + 1); sqj[0] = 1; for (i = 1; i <= nsqmax; ++i) { sqj[i] = sqj[i - 1] * sqnn; } } for (jxx = 0; jxx <= ixx; ++jxx) { /* set up sub-block quanta */ ijd = ixx - jxx; jbase = ibkptr[jxx]; n = ibkptr[jxx + 1] - jbase; kbgnj = ikmin[jxx]; jspn = ivs[jxx]; getqs(jspn, iff, n, kbgnj, jxcom, jscom, &jvbase); nj = jxcom[XNVAL]; nd = ni - nj; if (nd > ndmax) continue; ivmin = (ivbase < jvbase) ? ivbase : jvbase; ivsym = (unsigned int) (ivbase + jvbase + ivmin * glob.vibfac) << 2; ivcmp = ivsym + 3; ivsym += blksym(ixcom, jxcom); if (lastvv != ivcmp) { newblk = TRUE; lastvv = ivcmp; } isegy0 = FALSE; firstpar = TRUE; for (spar_now = spar_head[ivmin]; TRUE; spar_now = spar_now->next) { /* loop over parameters */ idflags = 0; if (spar_now != NULL) { ivcmp = spar_now->ipsym; if (ODD(spar_now->euler)) { if (ivcmp > lastvv) continue; if (ivcmp == lastvv) idflags = (int) spar_now->flags; } else { if (ivcmp > ivsym) continue; if (ivcmp == ivsym) idflags = (int) spar_now->flags; } } if ((idflags & MCOS_OK) == 0) { /* new cosine or end */ if (oldpar) { /* add diagonal or sub-diagonal to Hamiltonian */ /* using composite parameter value PBAR */ if (isegy0) { egy0 += pbar; } else { if (isneg) pbar = -pbar; if (isunit) { iz = idx[0]; jz = jdx[0]; pt = &t[jz * ndm]; for (i = 0; i < ncos; ++i) { pt[iz] += pbar; iz += ndmd; } } else { for (i = 0; i < ncos; ++i) { iz = idx[i]; jz = jdx[i]; t[iz + jz * ndm] += pbar * wk[i]; } } } oldpar = FALSE; } if (idflags == 0) break; /* end loop over parameters */ ncos = -1; } else { /* check if old cosine operator was zero */ if (ncos == 0) continue; } if (spar_now == NULL) break; ipar = spar_now->ip; nsqj = (int) spar_now->njq; if ((idflags & MNSQ) != 0 && nd == 0) { /* sum parameters which are different only by N*(N+1) */ if (parskp) continue; /* jump if mother operator was zero */ } else { /* new type of operator */ parskp = TRUE; if (sznz < 0) sznzfix(sznz, ni, nj, ixcom, jxcom, iscom, jscom); kl = idpars(spar_now, &ikq, &neuler, <, &ld, &kd, &ins, &si1, &si2, &sznz, &ifc, &alpha, &ldel); if (sznz > 0) { sznz = sznzfix(sznz, ni, nj, ixcom, jxcom, iscom, jscom); if (sznz == 0) continue; /* get next parameter */ } mkd = getmask(ixcom, jxcom, kd, ldel, kl, alpha); if (mkd == 0) continue; /* test tensor compatability */ npair = getll(0, ld, lt, kd, si1, si2, lscom, iscom, jscom); if (npair < 0) continue; /* get direction cosines and K dependent part of operator */ if (ncos < 0) { ncos = dircos(ixcom, jxcom, ld, kd, ndmx, wk, idx, jdx, ijd, kl, mkd, &isunit); if (ncos <= 0) { if (ncos == 0) continue; /* get next parameter */ printf("DIRCOS WORKING VECTOR TOO SHORT IN HAMX BY %d\n", ncos); exit(EXIT_FAILURE); } if ((idflags & MNOUNIT) != 0) isunit = FALSE; isgn = 1; if (isunit && wk[0] < 0.) isgn = -1; if (neuler != 0) { /*Euler operator */ ncos = specop(neuler, &newblk, &nsqj, &ikq, kbgni, kbgnj, ni, nj, ncos, wk, idx, jdx, par[ipar]); if (ncos == 0) continue; } if (ifc != 0) specfc(ifc, ivbase, jvbase, kd, kbgni, kbgnj, ncos, wk, idx, jdx); if (sznz != 0) sznzop(ni, nj, kbgni, kbgnj, iscom, jscom, ncos, wk, idx, jdx); if (ikq != 0) { /* find K*K dependence (anti-commutated) */ ncos = symksq(ikq, kbgni, kbgnj, ncos, wk, idx, jdx); if (ncos == 0) continue; } if (ibase != 0) { for (i = 0; i < ncos; ++i) idx[i] = (short) (idx[i] + ibase); } if (jbase != 0) { for (i = 0; i < ncos; ++i) jdx[i] = (short) (jdx[i] + jbase); } if (first) { isegy0 = FALSE; pbar = 0.; if (isunit && ixx == jxx) { /* diagonal unit operator */ oldpar = TRUE; if (ixx == 0) { isegy0 = TRUE; } else if (firstpar){ pbar = -egy0; firstpar = FALSE; } } if (ipar < 0) continue; /* get next parameter */ } else { /* begin by setting up reduced expectation values */ if (ipar < 0) continue; /* get next parameter */ dpmake(nsize, dbar, t, ncos, wk, idx, jdx, isunit); } } /* end for: get direction cosines and K dependent parts */ if (si2 < 0 && nd == 0) continue; /* commutator with N * N */ /* find K independent parts */ zpar = rmatrx(ld, lt, ixcom, jxcom); /* spherical tensor transformations */ tensor(&zpar, iscom, jscom, lscom, ismap, npair, alpha); if (sznz < 0) sznz = sznzfix(sznz, ni, nj, ixcom, jxcom, iscom, jscom); /* find N*(N+1) and N.S dependence */ symnsq(nsqj, ins, iscom, jscom, &zpar); parskp = (fabs(zpar) < 1e-30); if (parskp) continue; /* get next parameter */ if (si2 < 0) { /* commutator with N * N / 2 */ zpar *= ni; } } if (ni == 0 && nsqj != 0) continue; /* get next parameter */ ele = zpar * sqj[nsqj] * spar_now->zfac; if (isgn < 0) ele = -ele; if (first) { /* sum into composite parameter */ pbar += ele * par[ipar] ; oldpar = TRUE; } else { /* sum into derivative matrix */ i = ipder[ipar]; if (i < 0) { i = -1 - i; ele *= par[ipar]; i = ipder[i]; } daxpy(nsize, ele, dbar, 1, &dedp[i * ndm], 1); } } /* end parameter loop */ } /* end JXX loop */ } /* end IXX loop */ if (first) { if (ifdump) { for (i = 0; i < nsize; ++i) { pt = &t[i + i * ndm]; *pt += egy0; egy[i] = (*pt); pmix[i] = 0.; } return 0; } if (glob.idiag < 0 || nsize <= 1) { nd = nsize + 1; dcopy(nsize, t, nd, egy, 1); for (i = 0; i < nsize; ++i) { pt = &t[i * ndm]; dcopy(nsize, &zero, 0, pt, 1); pmix[i] = pt[i] = 1.; } } else { if (glob.idiag == 0) { roll = kroll(nsize, t, nsblk, ibkptr, ikmin); } else if (glob.idiag == 2) { nd = nsize + 1; dcopy(nsize, t, nd, egy, 1); ordham(nsize, egy, ibkptr, jdx); } /* diagonalize */ ierr = hdiag(nsize, nsize, t, egy, pmix); if (ierr != 0) { printf(" diagonalization failure, err = %d in block %d, %d\n", ierr, iblk, nsize); exit(EXIT_FAILURE); } if (glob.idiag == 0) { /* energy sort of subblock */ ordblk(nsize, nsize, t, egy, ibkptr, pmix, idx); } else if (glob.idiag == 1) { /* projection sort of full block */ for (i = 0; i <= nsize; ++i) { jdx[i] = (short) i; } ordblk(nsize, nsize, t, egy, jdx, pmix, idx); } else if (glob.idiag == 2) { ordblk(nsize, nsize, t, egy, ibkptr, pmix, idx); fixham(nsize, nsize, t, egy, pmix, jdx); } else if (glob.idiag == 3) { jxx = 0; ivbase = -1; for (i = 0; i < nsblk; ++i) { ixx = ibkptr[i]; jvbase = (int) ((unsigned) ivs[i] >> 2) - (itau[i] & 1); if (ivbase != jvbase) { jdx[jxx++] = (short) ixx; ivbase = jvbase; } wk[ixx] = (double) (itau[i] / 2); nd = ibkptr[i + 1] - 1; for (ii = ixx; ii < nd; ++ii) { wk[ii + 1] = wk[ii] + 4; } } jdx[jxx] = ibkptr[nsblk]; ordblk(nsize, nsize, t, egy, jdx, pmix, idx); ordham(nsize, wk, jdx, idx); fixham(nsize, nsize, t, egy, pmix, idx); } } if (isneg) { /* correct for oblate ordering */ for (i = 0; i < nsize; ++i) egy[i] = -egy[i]; } for (i = 0; i < nsize; ++i) /* add energy offset */ egy[i] += egy0; first = FALSE; } else { /* not first */ if (roll) { /* go back and get expectation values */ dcopy(nsize, &zero, 0, egy, 1); for (ipar = npar - 1; ipar >= 0; --ipar) { i = ipder[ipar]; if (i >= 0) { ele = par[ipar]; daxpy(nsize, ele, &dedp[i * ndm], 1, egy, 1); } } } break; } } while (!first); /* first loop */ i = 1; if (idpar[0] == (bcd_t)0) i = -1; /* this makes idpar referenced */ return i; } /* hamx */ int dclr(n1, n2, vec, ix) const int n1, n2, ix; double *vec; { /* clear a N1*N2 block */ static long nbig = 0x7ff0; long nsq; double *pvec; pvec = vec; nsq = n1 * (long) n2; while (nsq > nbig) { dcopy((int) nbig, &zero, 0, pvec, ix); pvec += nbig; nsq -= nbig; } dcopy((int) nsq, &zero, 0, pvec, ix); return 0; } /* dclr */ int specop(neuler, newblk, nsqj, ikq, ksi, ksj, ni, nj, ncos, wk, ix, jx, par) const int neuler, ksi, ksj, ni, nj, ncos; int *nsqj, *ikq; BOOL *newblk; double *wk; const double par; const short *ix, *jx; { /* subroutine to make Euler transform operators */ /* on entry: */ /* neuler= type of operator */ /* oldblk= flag for new block (=0 initially) */ /* NSQJ=power of N*N+N +1 */ /* IKQ= power of K*K */ /* KSI,KSJ= starting values of K */ /* NCOS= length of vector WK */ /* WK= vector of a sub-diagonal */ /* IX,JX= relative indices of elements in WK */ /* on return: */ /* WK= modified vector */ #define NSPOP 5 static double aden[NSPOP], bden[NSPOP]; double akisq, akjsq, anisq, anjsq, da, dd, di, dj, xi, xj; double sqi, sqj, di0, dj0; int i, nkq, nkq0, nnq, nnq0, ki, kj, ist; unsigned int ipwr; if (*newblk) { for (i = 0; i < NSPOP; ++i) { aden[i] = 0.; bden[i] = 0.; } *newblk = FALSE; } ist = (neuler - 2) >> 1; i = ist; if (i >= NSPOP) i = 0; nkq0 = (*ikq); nnq0 = (*nsqj); if (ODD(neuler)) { if (nnq0 == 0 && nkq0 == 1) { aden[i] = par; } else if (nnq0 == 1 && nkq0 == 0) { bden[i] = par; } return 0; } dd = bden[i]; da = aden[i] - dd; sqi = (double)(ni * (ni + 1)); di0 = 1. + dd * sqi; if (ni == nj) { sqj = sqi; dj0 = di0; }else { sqj = (double)(nj * (nj + 1)); dj0 = 1. + dd * sqj; } for (i = 0; i < ncos; ++i) { nkq = nkq0; nnq = nnq0; xi = wk[i]; ki = ksi + ((int) ix[i] << 1); akisq = (double) (ki * ki); di = di0 + da * akisq; akisq /= di; kj = ksj + ((int) jx[i] << 1); if (ki == kj && ni == nj) { if (nnq > 0) { /* N*N+N-K*K part */ ipwr = (unsigned int) nnq; anisq = sqi / di - akisq; for (;;) { if (ODD(ipwr)) xi *= anisq; if ((ipwr >>= 1) == 0) break; anisq *= anisq; } } if (nkq > 0) { /* K*K part */ ipwr = (unsigned int) nkq; for (;;) { if (ODD(ipwr)) xi *= akisq; if ((ipwr >>= 1) == 0) break; akisq *= akisq; } } if (ist != 0) xi /= di; /* extra part */ wk[i] = xi; } else { akjsq = (double) (kj * kj); dj = dj0 + da * akjsq; akjsq /= dj; xj = xi; if (nnq > 0) { /* N*N+N-K*K part */ ipwr = (unsigned int) nnq; anisq = sqi / di - akisq; anjsq = sqj / dj - akjsq; for (;;) { if (ODD(ipwr)) { xi *= anisq; xj *= anjsq; } if ((ipwr >>= 1) == 0) break; anisq *= anisq; anjsq *= anjsq; } } if (nkq > 0) { /* K*K part */ ipwr = (unsigned int) nkq; for (;;) { if (ODD(ipwr)) { xi *= akisq; xj *= akjsq; } if ((ipwr >>= 1) == 0) break; akisq *= akisq; akjsq *= akjsq; } } if (ist != 0) { xi /= di; xj /= dj; } /* extra part */ wk[i] = 0.5 * (xi + xj); } } *nsqj = 0; *ikq = 0; return ncos; } /* specop */ int specfc(ifc, iv, jv, kdel, ksi, ksj, ncos, wk, ix, jx) const int ifc, iv, jv, kdel, ksi, ksj, ncos; double *wk; const short *ix, *jx; { /* subroutine to make Fourier coefficient operators */ /* on entry: */ /* ITY= type of operator */ /* IV= v', JV=v" */ /* KDEL= change in K */ /* KSI,KSJ= starting values of K */ /* NCOS= length of vector WK */ /* WK= vector of a sub-diagonal */ /* IX,JX= relative indices of elements in WK */ /* on return: */ /* WK= modified vector */ static double pfac[MAXVIB]; static short ipfac[MAXVIB]; double di, dj, ang; int i, n, ki, kj, kd; if (ifc == 0) { if (ncos == 0) { for (i = 0; i < MAXVIB; ++i) { pfac[i] = 0.; ipfac[i] = 0; } } else { /* multiply by pi/3 */ pfac[iv] = wk[0] * acos(0.5); if (pfac[iv] < 0.) { pfac[iv] = -pfac[iv]; ipfac[iv] = 1; } } return 0; } n = ifc; if (n < 0) n = -1 - n; ang = (double) n; if (ODD(ipfac[iv] + ipfac[jv])) ang -= 0.5; di = ang * pfac[iv]; dj = ang * pfac[jv]; for (i = 0; i < ncos; ++i) { ki = ksi + (ix[i] << 1); kj = ksj + (jx[i] << 1); kd = ki - kj; if (kd < 0) kd = -kd; if (kdel != kd) kj = -kj; if (n != 0) { ang = di * ki + dj * kj; if (ifc < 0) { wk[i] *= sin(ang); } else { wk[i] *= cos(ang); } } else { wk[i] *= 0.5 * (ki + kj); } } return 0; } /* specfc */ int sznzfix(sznz, ni, nj, ixcom, jxcom, iscom, jscom) const int sznz, ni, nj; int *ixcom, *jxcom, *iscom, *jscom; { if (sznz <= 0) { if (sznz <= -4) { jscom[nspin] = nj << 1; jxcom[XNVAL] = nj; } else { iscom[nspin] = ni << 1; ixcom[XNVAL] = ni; } return 0; } else { switch (sznz) { case 1: if (ni + nj == 0) return 0; break; case 2: if (ni <= 1) return 0; iscom[nspin] -= 2; --ixcom[XNVAL]; break; case 3: if (ni <= 0) return 0; iscom[nspin] += 2; ++ixcom[XNVAL]; break; case 4: if (nj <= 1) return 0; jscom[nspin] -= 2; --jxcom[XNVAL]; break; case 5: if (nj <= 0) return 0; jscom[nspin] += 2; ++jxcom[XNVAL]; break; default: return 0; } return -sznz; } } /* sznzfix */ int sznzop(ni, nj, ksi, ksj, iscom, jscom, ncos, wk, ix, jx) const int ni, nj, ksi, ksj, ncos; const int *iscom, *jscom; double *wk; const short *ix, *jx; { /* subroutine to add SzNz */ /* on entry: */ /* ISCOM,JSCOM= spins */ /* NCOS= length of vector WK */ /* WK= vector of a sub-diagonal */ /* IX,JX= relative indices of elements in WK */ /* on return: */ /* WK= modified vector */ int iflg, i, k, kk, nni, nnj, nni0, nnj0, jj, iis; double vali, valj, twk, sum; if (ncos <= 0) return 0; nni = ni << 1; nni0 = iscom[nspin]; vali = (double) (nni0 + 1); nnj = nj << 1; nnj0 = jscom[nspin]; valj = (double) (nnj0 + 1); if (nni != nni0) { iflg = 1; vali = sqrt(vali * (nni + 1)); } else if (nnj != nnj0) { iflg = -1; valj = sqrt(valj * (nnj + 1)); } else { iflg = 0; } if (iflg <= 0) { jj = jscom[0]; iis = jscom[nspin + 1]; k = ksj + ksj + nnj0 + nnj + jj + iis; if ((k & 3) != 0) valj = -valj; valj *= c6jj(nnj0, 2, nnj, iis, jj, iis); } if (iflg >= 0) { jj = iscom[0]; iis = iscom[nspin + 1]; k = ksi + ksi + nni0 + nni + jj + iis; if ((k & 3) != 0) vali = -vali; vali *= c6jj(nni, 2, nni0, iis, jj, iis); } for (i = 0; i < ncos; ++i) { twk = wk[i]; sum = 0.; if (iflg <= 0) { k = ksj + ((int) jx[i] << 1); kk = k << 1; sum = valj * twk * k * c3jj(nnj0, 2, nnj, -kk, 0, kk); } if (iflg >= 0) { k = ksi + ((int) ix[i] << 1); kk = k << 1; sum += vali * twk * k * c3jj(nni, 2, nni0, -kk, 0, kk); } wk[i] = sum; } return 0; } /* sznzop */ unsigned int blksym(ixcom, jxcom) const int *ixcom, *jxcom; { /* get block symmetry from state symmetry */ return (unsigned int) (3 & (ixcom[XSYM] ^ jxcom[XSYM])); } /* blksym */ int ordham(nn, egy, isblk, iswap) const int nn; double *egy; const short *isblk; short *iswap; { /* subroutine to order eigenvalues within sub-block like diagonal of */ /* Hamiltonian */ double tmp; int iblk, i, ii, iq, is, inext; iblk = 1; i = 0; inext = isblk[1]; for (is = 1; is < nn; ++is) { if (is == inext) { inext = isblk[++iblk]; } else { /* find min energy */ tmp = egy[i]; iq = i; for (ii = is; ii < inext; ++ii) { if (egy[ii] < tmp) { tmp = egy[ii]; iq = ii; } } iswap[i] = (short) iq; if (iq > i) { egy[iq] = egy[i]; egy[i] = tmp; } } i = is; } return 0; } /* ordham */ int fixham(ndm, nn, t, egy, p, iswap) const int ndm, nn; double *t, *egy, *p; const short *iswap; { /* subroutine to order eigenvalues within sub-block like diagonal of */ /* Hamiltonian using permutation found with ordham */ double tmp, *tp; long ndml; int i, iq, is; is = nn - 2; ndml = ndm; if (is >= 0) { tp = &t[is * ndml]; for (i = is; i >= 0; --i) { iq = iswap[i]; if (iq > i) { tmp = egy[iq]; egy[iq] = egy[i]; egy[i] = tmp; tmp = p[iq]; p[iq] = p[i]; p[i] = tmp; dswap(nn, tp, 1, &t[iq * ndml], 1); } tp -= ndm; } } return 0; } /* fixham */ BOOL kroll(nsizd, t, nsblk, sbkptr, kmin) const int nsizd, nsblk; double *t; const short *sbkptr, *kmin; { /* subroutine to make sure that diagonal elements of the Hamiltonian */ /* are monotonically increasing (decreasing for oblate) */ /* after K=KTHRSH */ double vall, tlast, tmp, *ptmp; long ndmt; int ibgn, iend, i, k, n, ixx; BOOL roll; ndmt = nsizd + 1; roll = FALSE; for (ixx = 0; ixx < nsblk; ++ixx) { i = 5 - kmin[ixx]; ibgn = sbkptr[ixx]; if (i >= 2) ibgn += i >> 1; iend = sbkptr[ixx + 1] - 1; if (ibgn >= iend) continue; ptmp = &t[ibgn * ndmt]; tmp = tlast = *ptmp; ++ibgn; for (i = ibgn; i <= iend; ++i) { ptmp += ndmt; tlast = tmp; tmp = (*ptmp); if (tmp < tlast) break; } if (i > iend) continue; k = kmin[ixx] + ((i - sbkptr[ixx]) << 1); printf(" roll-over at K = %d\n", k); roll = TRUE; vall = 1e15; for (; i <= iend; ++i) { n = i - 1; dcopy(n, &zero, 0, &t[i], nsizd); n = nsizd - i; dcopy(n, &zero, 0, &t[i * ndmt], 1); tlast += vall; *ptmp = tlast; ptmp += ndmt; } } return roll; } /* kroll */ int getqs(im, iff, nsiz, kbgn, ixcom, iscom, iv) const int im, iff, nsiz, kbgn; int *ixcom, *iscom, *iv; { /* subroutine to get quanta corresponding to sub-block */ /* on entry: */ /* IM= position in sub-block mold */ /* IFF = F quantum number *2 */ /* NSIZ= size of sub-block */ /* KBGN= first K value */ /* on return: */ /* IXCOM = list of rotational containing */ /* 0: DIMENSION OF SUB-BLOCK */ /* 1: SYMMETRY CODE FOR BLOCK (0=A,1=Bx,2=By,3=Bz) */ /* 2: 'N' QUANTUM NUMBER */ /* 3: BEGINNING K */ /* 4: L VALUE */ /* 5: VIBRATIONAL QUANTA */ /* 5: ITOT SYMMETRY */ /* 5: EXTRA ITOT QUANTA */ /* ISCOM = list of angular momenta */ /* IV = vibrational quantum number if L=0,1 */ /* IV = vibrational quantum number-1 if L=-1 */ /* find list of angular momenta */ SVIB *pvinfo; int *isscom; short *iis, *jjs; int i, iiv, nset, nspinv, ispx, k, lupper; unsigned int mvs; mvs = (unsigned int) im; ispx = (int) (mvs >> glob.msshft); iiv = (int) (mvs >> 2) & glob.msmask; pvinfo = &vinfo[iiv]; iis = pvinfo->spt; nset = iis[0]; jjs = &iis[nset * ispx]; i = (int)(mvs & 3); iscom[0] = iff + (*jjs); ixcom[XNVAL] = iscom[0] >> 1; /* find N */ ixcom[XDIM] = nsiz; ixcom[XSYM] = i; /* find symmetry */ if (nsiz < 0) setgsym((int) pvinfo->gsym); k = kbgn; if (k < 0) { k = pvinfo->knmin[i]; i = pvinfo->knmin[3 - i]; if (k > i) k = i; if (k < 0) k = 0; } ixcom[XKBGN] = k; ixcom[XVIB] = iiv; ixcom[XISYM] = 0; lupper = pvinfo->lvupper; if (ODD(lupper)) --iiv; *iv = iiv; i = pvinfo->lvqn; ixcom[XLVAL] = i; if (i != 0 || k != 0) i = 1; ixcom[XIQN] = i; if (nspin == 0) return 1; nspinv = nset - 1; if (nspinv > nspin) nspinv = nspin; isscom = &iscom[nspin]; isscom[0] = iscom[0]; for (i = 1; i < nspinv; ++i) { iscom[i - 1] = iff + jjs[i]; isscom[i] = iis[i]; } iscom[i - 1] = iff; isscom[i] = iis[i]; while (i < nspin) { ++i; iscom[i - 1] = iff; isscom[i] = 0; } if (glob.nitot > 0) { if (glob.nitot >= 3) { if (itsym < nspinv) { i = jjs[itsym + 1]; iscom[itsym] = i; k = MOD(i >> 2, glob.nitot); i -= k << 2; if (is_esym[k] != 0) ++i; ixcom[XISYM] = k; ixcom[XIQN] |= i; } else { iscom[itsym] = 0; } } iscom[itptr] = jjs[itptr + 1]; } return ispx; } /* getqs */ int idpars(pspar, ksq, itp, ln, ld, kdel, ins, si1, si2, sznz, ifc, alpha, ldel) SPAR *pspar; int *ksq, *itp, *ln, *ld, *kdel, *ins, *si1, *si2, *sznz, *ifc, *alpha, *ldel; { /* subroutine to parse parameter identifiers to subfields */ /* on entry: */ /* IDPAR= 10*( parameter ID/1000) + ITP */ /* on return: */ /* KSQ= power of K*K */ /* ITP= parameter subtype */ /* LN= tensor order */ /* LD= order of direction cosine ( L <= LD ) */ /* KDEL= change in K */ /* INS=power of N.S */ /* SI1,SI2=spin identifiers */ int iret; *ksq = (int) pspar->ksq; *itp = (int) pspar->euler; *ln = (int) pspar->mln; *ld = (int) pspar->mld; *kdel = (int) pspar->mkdel; *ins = (int) pspar->mins; *si1 = (int) pspar->msi1; *si2 = (int) pspar->msi2; *sznz = (int) pspar->msznz; *ifc = (int) pspar->fc; *alpha = (int)pspar->alpha; *ldel = (int)pspar->mldel; iret = (int)pspar->flags; return iret; } /* idpars */ int getll(llf, ld, ln, kd, si1, si2, lscom, iscom, jscom) const int llf, ld, ln, kd, si1, si2; int lscom[]; const int iscom[], jscom[]; { /* function to find sequence of L values for coupling of spins */ /* on entry: */ /* LLF = 2 * tensor order for F */ /* LD = order of direction cosine */ /* LN = order of N dependent tensor */ /* SI1, SI2 = position of spin values */ /* ISCOM,JSCOM = list of 2 * angular momenta */ /* on return: */ /* LSCOM = list of 2 * tensor orders for succesive couplings */ /* followed tensor orders for spins */ /* return true if delta QN not consistent with L */ int *lsscom; int i, llj, lld, llmax, jdif, jsum, maxspin, ierr, nm2; /* check tensor sum rules for N */ ierr = -1; llmax = 0; lld = ld << 1; llj = ln << 1; jdif = iscom[nspin] - jscom[nspin]; if (jdif < 0) jdif = -jdif; if (lld >= llj) { if (llj < jdif) return ierr; llmax = kd << 1; if (llmax < lld) llmax = lld; } else { if (lld < jdif) return ierr; llmax = kd << 1; if (llmax < llj) llmax = llj; } jsum = iscom[nspin] + jscom[nspin]; if (llmax > jsum) return ierr; if (nspin == 0) { lscom[0] = llj; return 0; } maxspin = si1; if (si1 > itsym) maxspin = nspin; lsscom = &lscom[nspin]; lsscom[0] = llj; for (i = 1; i <= nspin; ++i) { lsscom[i] = 0; } if (si1 > 0) { lsscom[si1] = 2; if (si2 > 0) lsscom[si2] += 2; } /* generate list of tensor orders, checking quantum numbers along the way */ nm2 = nspin - 2; for (i = 0; i <= nm2; ++i) { if (i == itsym) { /* special for Itot */ if (si1 <= itsym) { llj = 0; } else if (si2 <= itsym || si1 == si2) { llj = lsscom[si1]; } i = nm2; } else if (i >= maxspin) { /* operator complete */ llj = llf; } else { llj -= lsscom[i + 1]; /* reduce order by spin operator order */ if (llj < 0) llj = -llj; } lscom[i] = llj; jdif = iscom[i] - jscom[i]; if (jdif < 0) jdif = -jdif; if (llj < jdif) return ierr; jsum = iscom[i] + jscom[i]; if (llj > jsum) return ierr; } lscom[nspin - 1] = llf; if (llf != 0) maxspin = nspin; return maxspin; } /* getll */ int getmask(const int *xbra, const int *xket, const int kd, const int ldel, const int loff, const int alpha) { /* CREATE MASK FOR DIRCOS */ /* KD IS TOTAL K QUANTIUM CHANGE */ int ldif, lbra, lket, mask, kbra, kket, nitot, kk, ksbra, ksket; mask = 7; nitot = glob.nitot; if (nitot >= 3) { /* check Itot symmetry */ ksbra = xbra[XISYM]; ksket = xket[XISYM]; kk = (xbra[XIQN] ^ xket[XIQN]) >> 1; if ((loff & MIDEN) != 0) { /* Identity operator */ if (kk != 0) return 0; } else if (EVEN(kk)) { /* A1|A1, A1|E, E|E, A2|A2, ..*/ if ((loff & MSYM2) != 0) return 0; } else { /* kk is ODD */ if ((loff & MSYM2) == 0) return 0; /* */ if (is_esym[ksbra] == 0 && is_esym[ksket] == 0) return 0; /* == 0 */ } mask = 0; kk = ksket - ksbra; if (alpha == 0) { if (kk == 0) mask = 5; } else if (alpha < nitot) { if (MOD(kk + alpha, nitot) == 0) mask = 1; if (MOD(kk - alpha, nitot) == 0) mask += 4; } else { /* quasi-diagonal */ if (MOD(ksbra + ksket - alpha, nitot) == 0) mask = 2; } if (mask == 0) return mask; } /* mask bit 0 (val = 1): kbra = kket + kd, lbra = lket + ldel */ /* mask bit 2 (val = 4): kbra = kket - kd, lbra = lket - ldel */ /* mask bit 1 (val = 2): kbra = kd - kket, lbra = ldel - lket */ lbra = xbra[XLVAL]; lket = xket[XLVAL]; ldif = lbra - lket; if (ldel == 0) { /* operator diagonal in l */ if (ldif != 0) mask &= 2; /* clear bits 0 and 2 */ } else { /* operator off-diagonal in l */ if (ldif == 0) { mask &= 2; /* clear bits 0 and 2 */ } else if (ldif != ldel) { mask &= 6; /* clear bit 0 */ } else if (-ldif != ldel) { mask &= 3; /* clear bit 2 */ } } if ((mask & 2) == 0) return mask; /* bit 1 not set */ mask &= 5; /* clear bit 1 temporarily */ /* quasi-diagonal, KD = KBRA + KKET */ if (lbra + lket != ldel) return mask; kbra = xbra[XKBGN]; kket = xket[XKBGN]; kk = kd - kbra - kket; if (kk < 0) return mask; if (EVEN(xket[XIQN])) kk -= 2; if (EVEN(xbra[XIQN])) kk -= 2; if (kk < 0) return mask; return (mask + 2); } /* getmask */ double rmatrx(ld, lv, ixcom, jxcom) const int ld, lv; const int *ixcom, *jxcom; { /* find N corrections to reduced matrix elements */ double dlsq, dnsq, dn, del, tmp; int nsum, ndif, ni, nj, lmin, lt; tmp = 1.; if (ld == lv) return tmp; ni = ixcom[XNVAL]; nj = jxcom[XNVAL]; nsum = ni + nj + 1; dn = del = (double) nsum; dnsq = dn * dn; if (ld > lv) { /* contracting tensor order */ lmin = lv; lt = ld; if (lmin == 0) { lmin = 1; tmp = ni * (double) (ni + 1) / dn; } } else { /* expanding tensor order */ lmin = ld; lt = lv; if (lmin == 0) { lmin = 1; tmp = ni * (double) (ni + 1) * dn; } } if (lt >= nsum) return 0.; ndif = ni - nj; if (ndif != 0) del = (double) (ndif * ndif); while (lt > lmin) { dlsq = (double) (lt * lt); tmp *= 0.25 * (dnsq - dlsq); if (ndif != 0) tmp -= tmp * del / dlsq; --lt; } return sqrt(tmp); } /* rmatrx */ int symnsq(inq, ins, iscom, jscom, zval) const int inq, ins, *iscom, *jscom; double *zval; { /* subroutine to find powers of N*(N+1) and N.S */ /* on entry: */ /* INQ= power of N*(N+1) */ /* INS= power of N.S */ /* ISCOM,JSCOM= list of angular momenta *2 */ /* ISPN = S *2 */ /* Z= matrix multipling factor */ /* on return: */ /* Z= matrix factor */ double x, dotns, dotnsp, zr, zl; int modf, ispn, nni, nnj, j; unsigned int ipwr; nni = iscom[nspin]; nnj = jscom[nspin]; modf = nni - nnj; if (inq <= 0) modf = 0; if (modf == 0 && ins == 0) return 0; zl = zr = 0.5 * (*zval); if (modf != 0) { /* CORRECT ELE FOR COMMUTATION */ if (nnj == 0) { zr = 0.; } else { x = nnj * (double) (nnj + 2) / (nni * (double) (nni + 2)); ipwr = (unsigned int) inq; for (;;) { if (ODD(ipwr)) zr *= x; if ((ipwr >>= 1) == 0) break; x *= x; } } } if (ins > 0) { j = iscom[0]; ispn = iscom[nspin + 1]; x = (double) (ispn * (ispn + 2)); dotns = 0.125 * ((j - nni) * (double) (j + nni + 2) - x); j = jscom[0]; ispn = jscom[nspin + 1]; x = (double) (ispn * (ispn + 2)); dotnsp = 0.125 * ((j - nnj) * (double) (j + nnj + 2) - x); ipwr = (unsigned int) ins; for (;;) { if (ODD(ipwr)) { zl *= dotns; zr *= dotnsp; } if ((ipwr >>= 1) == 0) break; dotns *= dotns; dotnsp *= dotnsp; } } *zval = zl + zr; return 0; } /* symnsq */ int symksq(ikq, ksi, ksj, n, wk, ix, jx) const int ikq, ksi, ksj, n; double *wk; short *ix, *jx; { /* subroutine to multiply right and left by a power of K*K */ /* on entry: */ /* IKQ= power of K*K */ /* KSI,KSJ= starting values of K */ /* N= length of vector WK */ /* WK= vector of a sub-diagonal */ /* IX,JX= relative indices of elements in WK */ /* note: if IKQ < 0 then operator is only define for K=-1-IKQ */ /* on return: */ /* WK= modified vector */ double akisq, akjsq, xi, xj; int kcmp, i, ki, kj, zflg; unsigned int nkq; kcmp = -1 - ikq; zflg = 0; for (i = 0; i < n; ++i) { xi = wk[i]; ki = ksi + ((int) ix[i] << 1); kj = ksj + ((int) jx[i] << 1); if (ikq > 0) { nkq = (unsigned int) ikq; if (ki != kj) { zflg = n; xj = xi; akisq = (double) (ki * ki); akjsq = (double) (kj * kj); for (;;) { if (ODD(nkq)) { xi *= akisq; xj *= akjsq; } if ((nkq >>= 1) == 0) break; akisq *= akisq; akjsq *= akjsq; } wk[i] = 0.5 * (xi + xj); } else if (ki > 0) { zflg = n; akisq = (double) (ki * ki); for (;;) { if (ODD(nkq)) xi *= akisq; if ((nkq >>= 1) == 0) break; akisq *= akisq; } wk[i] = xi; } else { wk[i] = 0.; } } else if (ki == kcmp && ki == kj) { ix[0] = ix[i]; jx[0] = jx[i]; wk[0] = xi; return 1; } } return zflg; } /* symksq */ int dpmake(nsize, dp, t, n, wk, idx, jdx, isunit) const int nsize, n; const BOOL isunit; double *dp; const double *t, *wk; const short *idx, *jdx; { /* find derivative contribution from sub-diagonal */ const double *pt; double ele, tz; int i, k, iz, jz; if (isunit) { iz = idx[0]; jz = jdx[0]; if (iz != jz) { /* off-diagonal unit matrix */ for (k = 0; k < nsize; ++k) { dp[k] = 2. * ddot(n, &t[iz], 1, &t[jz], 1); iz += nsize; jz += nsize; } } else { /* diagonal unit matrix */ for (k = 0; k < nsize; ++k) { pt = &t[iz]; dp[k] = ddot(n, pt, 1, pt, 1); iz += nsize; } } } else { dcopy(nsize, &zero, 0, dp, 1); for (i = 0; i < n; ++i) { iz = idx[i]; jz = jdx[i]; ele = wk[i]; if (iz != jz) { ele *= 2; for (k = 0; k < nsize; ++k) { dp[k] += t[iz] * t[jz] * ele; iz += nsize; jz += nsize; } } else { for (k = 0; k < nsize; ++k) { tz = t[iz]; dp[k] += tz * tz * ele; iz += nsize; } } } } return 0; } /* dpmake */ int tensor(zval, iscom, jscom, lscom, imap, npair, alpha) double *zval; const int iscom[], jscom[], lscom[], imap[]; int npair, alpha; { /* subroutine to find spherical tensor spin coupling coefficients */ /* on entry: */ /* Z= initial multiplier */ /* NS= number of spins */ /* ISCOM,JSCOM= list of angular momenta *2 */ /* LLCOM= list of tensor orders *2 */ /* LSCOM = list of tensor orders for spins *2 */ /* on return: */ /* Z= modified multiplier */ double zsq, z; int isgn, i, jjf, jji, nnf, llj, nni, lln, lls, ssf, ssi, ix, jj[9]; if (npair <= 0) return 0; z = (*zval); isgn = 0; zsq = 1.; for (i = 0; i < npair; ++i) { ix = imap[i + i]; lln = lscom[ix]; nni = iscom[ix]; nnf = jscom[ix]; ix = imap[i + i + 1]; lls = lscom[ix]; ssi = iscom[ix]; ssf = jscom[ix]; if (alpha >= 0 && ix == itptr) { ix = itsym + nspin + 1; getzitot(&z, lls, iscom[ix], &lscom[ix], &iscom[itsym], &jscom[itsym], alpha, glob.nitot); i = nspin - 1; } llj = lscom[i]; jji = iscom[i]; jjf = jscom[i]; if (lln != 0) { if (llj == 0) { /* T * U */ isgn += jji + ssi + nnf; z *= c6jj(nni, lln, nnf, ssf, jji, ssi); } else if (lls == 0) { /* U = 1 */ isgn += jjf + nni + ssi + llj; zsq *= jji + 1; zsq *= jjf + 1; z *= c6jj(nni, llj, nnf, jjf, ssi, jji); } else { /* T X U */ isgn += 2; zsq *= 2.5; ix = (jjf + 1) * (llj + 1); zsq *= ix; zsq *= jji + 1; jj[0] = nnf; jj[1] = lln; jj[2] = nni; jj[3] = ssf; jj[4] = lls; jj[5] = ssi; jj[6] = jjf; jj[7] = llj; jj[8] = jji; z *= c9jj(jj); } } else if (lls != 0) { /* T = 1 */ isgn += jji + nni + ssf + llj; zsq *= jji + 1; zsq *= jjf + 1; z *= c6jj(ssi, llj, ssf, jjf, nni, jji); } } if (ODD2(isgn)) z = -z; *zval = z * sqrt(zsq); return 0; } /* tensor */ int setint(lu, ifdiag, nsav, ndip, idip, isimag) FILE *lu; BOOL *ifdiag; const int ndip; int *nsav, *isimag; bcd_t *idip; { /* subroutine to initialize intensity data */ SVIB *pvib1, *pvib2; short *iiv1, *iiv2; double zfac; unsigned int ijv; int i, icase, isym, jsym, ivdif, ii, iv1, iv2, lv1, lv2, ifac, lv, nsym; int si1, ibcd, ndecv, nbcd, ipty, iflg, kd, ld, ifc, imag, ldel, jj; bcd_t itmp; BOOL bad_dip; *ifdiag = ((glob.idiag + 1) != 0); *nsav = 1; free(dipinfo); dipinfo = (SDIP *) mallocq((size_t) ndip * sizeof(SDIP)); nddip = ndip; dipinfo[0].fac = 1.; dipinfo[0].flg = (short) 0; dipinfo[0].kd = C0; dipinfo[0].ld = C0; dipinfo[0].ldel = C0; dipinfo[0].fc = C0; ifac = glob.vibfac + 1; ndecv = glob.vibdec; nbcd = (int) idip[0] & 0x7f; for (i = 0, ibcd = 0; i < ndip; ++i, ibcd += nbcd) { ii = (int)idip[ibcd + ndecv + 1]; si1 = (ii >> 4) & 0x0f; ii &= 0x0f; if (ndecv > 2) ii = 100 * ii + bcd2i(idip[ibcd + 3]); if (ndecv > 1) ii = 100 * ii + bcd2i(idip[ibcd + 2]); itmp = idip[ibcd + 1]; isym = (int) itmp & 0x0f; ii = 10 * ii + ((int)(itmp >> 4) & 0x0f); iv2 = ii / ifac; iv1 = ii - iv2 * ifac; icase = bcd2i(idip[ibcd + ndecv + 2]); ivdif = iv1 - iv2; if (ivdif < 0) { ii = iv1; iv1 = iv2; iv2 = ii; } pvib1 = &vinfo[iv1]; pvib2 = &vinfo[iv2]; lv1 = pvib1->lvupper; lv2 = pvib2->lvupper; if (ODD(lv1)) --iv1; if (ODD(lv2)) --iv2; ipty = ((lv1 ^ lv2) & 2) >> 1; lv1 = pvib1->lvqn; lv2 = pvib2->lvqn; /* MAKE SURE VIBRATIONS ARE ORDERED */ bad_dip = TRUE; ifc = 0; zfac = 1.; do { /* one pass */ /* CONVERT SYMMETRY */ if (isym > 3 || icase > MAXINT) break; if (iv1 >= glob.nvib) break; if (iv2 >= glob.nvib) break; if ((lv1 + lv2) < 0) break; if (pvib1->spt != pvib2->spt) break; if (pvib1->gsym != pvib2->gsym) break; if (pvib1->ewt[0] != pvib2->ewt[0]) break; if (pvib1->ewt[1] != pvib2->ewt[1]) break; if (glob.oblate == FALSE) isym = revsym[isym]; iflg = 1; if (isym == 0) { ld = 0; zfac = 0.009274; } else { ld = 1; iflg |= MELEC; } imag = 0; kd = isoddk[isym]; /* ld = 0,1; kd = 0, 1, 1, 0 */ if (icase > 0) { switch (icase) { case 1: /* ld = 2,2; kd = 0, 1, 1, 2 */ ld = 2; if (isym == 3) kd = 2; break; case 2: /* ld = 2,2; kd = 2, 1, 1, 2 */ case 11: case 12: ld = 2; if (kd == 0) kd = 2; if (icase == 11) { iflg |= MINOQ; ++imag; /* commutator */ } else if (icase == 12) { ifc = -1; ++imag; isym = 3 - isym; } break; case 5: if (isym != 0) { iflg |= MINOQ; ++imag; /* commutator */ } break; case 6: /* ld = 0,3; kd = 2, 3, 3, 2 */ if (isym != 0) ld = 3; kd += 2; break; case 7: case 8: ifc = icase * 10 + si1 - 70; if (ifc > 10) { ifc = 9 - ifc; ++imag; isym = 3 - isym; } si1 = 0; break; case 10: iflg |= MINOQ; /* double commutator */ break; } if ((iflg & MINOQ) != 0 && isym == 0) break; } if (icase == 11 && isym == 0) break; for (ii = 0; ii < 4; ++ii) { jsym = ii ^ isym; if (pvib1->wt[ii] != pvib2->wt[jsym]) { jj = pvib1->wt[4]; if (jj == 0) break; if (pvib1->wt[ii] != pvib2->wt[jsym ^ jj]) break; } } if (ii < 4) break; if (pvib1->wt[4] != pvib2->wt[4]) break; imag += ld; ldel = lv1 - lv2; nsym = setgsym((int)pvib1->gsym); if (ldel != 0) { if (icase == 5 && isym == 0) break; if (ldel < 0) ldel = -ldel; if (ivdif < 0) /* DELTA K * DELTA L < 0 */ ldel = -ldel; } else if (icase == 5 && isym == 0) { if (lv1 == 0 || si1 > 0) break; iflg |= MLZ; ld = 1; } if (ODD(imag)) { iflg |= MODD; } else if (isym > 0) { if (iv1 == iv2 && ldel == 0) break; } if (glob.nitot >= 3) { if (MOD(kd - ldel, glob.nitot) != 0) break; /* require alpha == 0 */ if (si1 <= itsym) iflg |= MIDEN; /* identity operator */ } lv = 1; if (si1 > 0) { /* check spin properties */ if (ld == 0 && isym != 0) break; iiv1 = pvib1->spt; iiv2 = pvib2->spt; if (checksp(TRUE, si1, 0, iiv1, iiv2, &zfac) != 0) break; lv = iflg & MELEC; /* simple electric dipole, lv = 2 */ } if (ld == kd && ld < lv) ld -= (ld - lv) & 0x7fe; /* subtract even number */ isimag[i] = (isym > 0)? 1: 0; dipinfo[i].fac = zfac; dipinfo[i].flg = (short) iflg; dipinfo[i].kd = (signed char) kd; dipinfo[i].ld = (signed char) ld; dipinfo[i].ldel = (signed char) ldel; dipinfo[i].fc = (signed char) ifc; bad_dip = FALSE; } while (FALSE); if (bad_dip) { putbcd(sbcd, NSBCD, &idip[ibcd]); fprintf(lu, " WARNING: dipole %3d %s has no matrix elements\n", (i + 1), sbcd); dipinfo[i].fac = 0.; dipinfo[i].flg = (short)0; dipinfo[i].kd = C0; dipinfo[i].ld = C0; dipinfo[i].ldel = C0; dipinfo[i].fc = C0; iv2 = iv1 = glob.vibfac; isym = 0; } ii = iv2 * ifac + iv1; ijv = ((unsigned int) ii << 2) + (unsigned int) isym; if (ndecv < 3) { idip[ibcd + 5] = idip[ibcd + ndecv + 2]; idip[ibcd + 4] = idip[ibcd + ndecv + 1]; idip[ibcd + 3] = (bcd_t)0; } else { idip[ibcd + 4] &= 0xf0; idip[ibcd + 3] = (bcd_t) ((ijv >> 16) & 0xff); } idip[ibcd + 2] = (bcd_t) (ijv >> 8); idip[ibcd + 1] = (bcd_t) ijv; if (si1 == 0) idip[ibcd + 4] &= 0x0f; } return 0; } /* setint */ int intens(iblk, isiz, jblk, jsiz, ndip, idip, dip, s) const int iblk, isiz, jblk, jsiz, ndip; const bcd_t *idip; const double *dip; double *s; { /* function to calculate intensities */ /* on entry: */ /* IBLK,JBLK= block numbers */ /* ISIZ,JSIZ= block sizes */ /* IDIP= dipole id codes , IDIP(0)= number of dipoles */ /* DIP= dipole values */ /* on return: */ /* S= transition dipole matrix in original basis */ /* INTENS returns true is S is not 0 */ static int idipoff; SDIP *pdip; double dd; long ndms; unsigned int ijv; int ncos, ifup, isym, i, k, lv, n, icase, ibase, mkd, ibcd, nbcd; int jbase, kbgni, kbgnj, nblki, nblkj, ixtmp, iret, kd, ld, npair, ldel; int ix, jx, nx, si1, si2, iff, jff, nni, nnj, iwt[5], jwt[5], ndecv; int ivv, jvv, ixx, jxx, ifc, ksym, kl, ioff, joff, alpha, nitot, dipoff; BOOL isunit; bcd_t bijv1, bijv2, bijv3; if (ndmx <= 0) { puts("working vectors not allocated"); exit(EXIT_FAILURE); } dipoff = idipoff; if (dipoff >= nddip) dipoff = 0; idipoff = dipoff + ndip; iret = 0; cgetv[0].cblk = 0; cgetv[1].cblk = 0; /* get quantum information */ nblki = getqq(iblk, &iff, iwt, ibkptr, ikmin, ivs); ioff = glob.maxblk; nblkj = getqq(jblk, &jff, jwt, &ibkptr[ioff], &ikmin[ioff], &ivs[ioff]); ixx = iff - jff; if (ixx < -2 || ixx > 2) return 0; ixx = iff + jff; if (ixx < 2 || ODD(ixx)) return 0; if (iwt[0] != jwt[0]) return 0; if (iwt[3] != jwt[3]) return 0; nitot = glob.nitot; alpha = 0; bijv3 = (bcd_t) 0; /* clear dipole matrix */ dclr(isiz, jsiz, s, 1); /* loop over sub-blocks */ ndms = isiz; ndecv = glob.vibdec; nbcd = (int) idip[0] & 0x7f; for (ixx = 0; ixx < nblki; ++ixx) { ibase = ibkptr[ixx]; n = ibkptr[ixx + 1] - ibase; kbgni = ikmin[ixx]; ixtmp = ivs[ixx]; getqs(ixtmp, iff, n, kbgni, ixcom, iscom, &ivv); nni = iscom[0]; for (jxx = 0; jxx < nblkj; ++jxx) { joff = jxx + ioff; jbase = ibkptr[joff]; n = ibkptr[joff + 1] - jbase; kbgnj = ikmin[joff]; ixtmp = ivs[joff]; getqs(ixtmp, jff, n, kbgnj, jxcom, jscom, &jvv); nnj = jscom[0]; nx = ((nni > nnj) ? nni : nnj) >> 1; i = (ivv < jvv) ? ivv : jvv; ijv = (unsigned int) (ivv + jvv + i * glob.vibfac); ijv = (ijv << 2) + blksym(ixcom, jxcom); bijv1 = (bcd_t) ijv; bijv2 = (bcd_t)(ijv >> 8); if (ndecv == 3) bijv3 = (bcd_t)(ijv >> 16); for (i = 0, ibcd = 0; i < ndip; ++i, ibcd += nbcd) { /* set dipole for transition */ if (bijv1 != idip[ibcd + 1]) continue; if (bijv2 != idip[ibcd + 2]) continue; if (bijv3 != idip[ibcd + 3]) continue; pdip = &dipinfo[dipoff + i]; kl = pdip->flg; if ((kl & MINOQ) != 0 && nni == nnj) continue; isym = (int)(bijv1 & 3); si1 = (int) (idip[ibcd + 4] >> 4) & 0x0f; si2 = 0; icase = bcd2i(idip[ibcd + 5]); kd = (int) pdip->kd; ifc = (int) pdip->fc; ldel = (int) pdip->ldel; lv = 1; ld = (int) pdip->ld; if (si1 > 0) lv = kl & MELEC; ksym = 0; if (nitot >= 3 && kd > 0) { alpha = nitot; ksym = 1; } do { if (ksym == 0) { alpha = 0; } /* find matrix elements */ mkd = getmask(ixcom, jxcom, kd, ldel, kl, alpha); if (mkd == 0) continue; npair = getll(2, ld, lv, kd, si1, si2, lscom, iscom, jscom); if (npair < 0) continue; ifup = 1 - (kl & MELEC); ncos = dircos(ixcom, jxcom, ld, kd, ndmx, wk, idx, jdx, ifup, kl, mkd, &isunit); if (ncos == 0) continue; if (ncos < 0) { puts("DIRCOS WORKING VECTOR TOO SHORT IN INTENS"); exit(EXIT_FAILURE); } dd = dip[i] * pdip->fac; /* correct for special cases */ switch (icase) { case 3: if (nni == nnj) { k = nx * (nx + 1); } else { k = nx * nx; } dd *= k; break; case 4: symksq(1, kbgni, kbgnj, ncos, wk, idx, jdx); break; case 5: case 11: if (ifup > 0) break; if (nni < nnj) dd = -dd; dd *= nx; break; case 9: if (nni == nnj) { k = nx * (nx + 1); dd *= k; } else { k = nx * nx; dd *= k + 1; } dd *= k; break; case 10: dd *= nx * nx; break; case 7: case 8: case 12: specfc(ifc, ivv, jvv, kd, kbgni, kbgnj, ncos, wk, idx, jdx); break; } /* correct for reduced matrix of N */ if (ld != lv) dd *= rmatrx(ld, lv, ixcom, jxcom); /* couple dipoles through the spins */ tensor(&dd, iscom, jscom, lscom, ismap, npair, alpha); for (k = 0; k < ncos; ++k) { ix = idx[k] + ibase; jx = jdx[k] + jbase; s[ix + jx * ndms] += dd * wk[k]; } iret = iff + 1; } while (--ksym >= 0); } } } return iret; } /* intens */ int setopt(luin, nfmt, itd, nbcd, namfil) FILE *luin; int *nfmt, *itd, *nbcd; char *namfil; { /* this subroutine reads in quantum quantities */ /* on entry: */ /* LUIN is the fortran unit number to read from */ /* NFMT is the maximum number of quanta */ /* on return: */ /* NFMT= number of quantum number formats for catalog */ /* ITD= 2 for linear, 3 for non-linear */ /* NAMFIL contains a file name on which to find parameter names */ /* SETOPT = number of cards read, -1 for end-of-file */ #define NVEC 10 #define NDECSPIN 11 #define MAXIAX 11 static int isymv[] ={2,2,2,2,2,2,3,4,4,5,6,6}; static int iaxv[] ={2,1,2,3,4,5,2,4,5,2,4,5}; SVIB *pvinfo, *pvinfom; double rvec[NVEC], rvec0[NVEC], vsym; size_t nl; int ivib, lopt, i, ii, j, k, nvibm, isym, iend, ivsym; int iwtpl, iwtmn, nt, ncards, knnmin, knnmax, si2, iax; int spinsgn, ewt0, ewt1, nspinqn; int ivwt[3]; bcd_t bcdspin[NDECSPIN]; char card[82], ctyp; cgetv[0].cblk = cgetv[1].cblk = 0; /* reset store for getqn */ cjjini(); glob.lsym = TRUE; glob.esym = TRUE; glob.esymdec = 100; rvec0[0] = 1; rvec0[1] = 0; rvec0[2] = (double)MAXN_DIRCOS; rvec0[3] = 0; rvec0[4] = 1; rvec0[5] = 1; rvec0[6] = 1; rvec0[7] = 0; rvec0[8] = 99; rvec0[9] = 0; ncards = nvibm = nt = si2 = nspinqn = 0; *itd = 2; bcdspin[0] = (bcd_t)NDECSPIN; /* read option cards */ do { if (fgetstr(card, 82, luin) <= 0) return -1; ++ncards; ctyp = C0; if (isalpha((int) card[0])) ctyp = card[0]; else if (isalpha((int) card[1])) ctyp = card[1]; iend = getbcd(card, bcdspin, NDECSPIN); spinsgn = NEGBCD(bcdspin[0]); dcopy(NVEC, rvec0, 1, rvec, 1); rvec[7] = 0.; pcard(&card[iend], rvec, NVEC, NULL); lopt = (int) rvec[0]; knnmin = (int) rvec[1]; knnmax = (int) rvec[2]; iax = (int) rvec[4]; iwtpl = (int) rvec[5]; iwtmn = (int) rvec[6]; vsym = rvec[7]; ewt0 = (int) rvec[8]; ewt1 = 0; if (vsym < -0.5) { ivsym = -1; vsym = 0.; } else if (vsym > 0.5 && vsym < 2.e+15) { ivsym = 1; } else { ivsym = 0; vsym = 0.; } ivib = lopt; if (ivib < 0) ivib = -ivib; if (knnmax < 0 || knnmax > MAXN_DIRCOS) knnmax = MAXN_DIRCOS; if (knnmin > knnmax) knnmax = knnmin; if (knnmin != knnmax) *itd = 3; if (ncards == 1) { dcopy(NVEC, rvec, 1, rvec0, 1); strcpy(namfil, "sping.nam"); if (ctyp != C0) namfil[4] = ctyp; glob.ixz = (int) rvec[3]; glob.idiag = (int) rvec[9]; if (glob.idiag > 3) { glob.idiag = 0; } else if (glob.idiag < 0) { glob.idiag = -1; } if (ivib > 80) { if (ivib > MAXVIB) ivib = MAXVIB; if (ivib > 80 && sizeof(int) == 2) ivib = 80; } else if (ivib <= 0) { ivib = 1; } if (glob.nvib > 1) { free(vinfo); vinfo = NULL; } vinfo1.spt = sptzero; if (ivib > 1) { nvibm = ivib - 1; nl = (size_t) ivib * sizeof(SVIB); vinfo = (SVIB *) mallocq(nl); memcpy(vinfo, &vinfo1, sizeof(SVIB)); } else { vinfo = &vinfo1; } for (k = 0; k <= 4; ++k) vinfo->wt[k] = 0; vinfo->ewt[0] = vinfo->ewt[1] = 0; glob.nvib = ivib; if (ivib <= 9) { glob.vibfac = 9; glob.vibdec = 1; k = 4; /* 16 */ } else if (ivib <= 99) { glob.vibfac = 99; glob.vibdec = 2; k = 7; /* 128 */ } else { glob.vibfac = 999; glob.vibdec = 3; k = 10; /* 1024 */ } glob.msshft = (unsigned int)(k + 2); /* msshft = shifts needed for vib and symmetry */ glob.msmask = (int)(1 << (unsigned int) k) - 1; ivib = 0; glob.oblate = (lopt < 0); glob.nqnn = 3; if (spinsgn != 0) glob.nqnn = 1; if (ewt0 < 0) glob.esymdec = 1000; } else if (ivib >= glob.nvib) { /* if IVIB out of range read another card */ continue; } pvinfo = &vinfo[ivib]; k = iax; if (iax < 0) iax = -iax; if (iax > MAXIAX) iax = 0; if (iax <= 3 && glob.oblate) iax = revsym[iax]; isym = isymv[iax]; iax = iaxv[iax]; pvinfo->gsym = (short) (isym << 1); /* set up minimum K values */ pvinfo->knmax = (short) knnmax; pvinfo->knmin[0] = 0; pvinfo->knmin[1] = 1; pvinfo->knmin[2] = 1; pvinfo->knmin[3] = 2; if (k < 0) { /* iax was negative */ if (isym >= 3) { pvinfo->knmin[0] = -2; pvinfo->knmin[3] = -2; } pvinfo->gsym |= 1; } if (knnmax != 0 && glob.nqnn == 1) glob.nqnn = 2; /* set ewt */ if (ewt0 < 0) ewt0 = -ewt0; if (ewt0 >= glob.esymdec) { ewt1 = ewt0 / glob.esymdec; ewt0 -= glob.esymdec * ewt1; } if (isym < 3) { ewt0 = glob.esymdec - 1; glob.esym = FALSE; } else if (isym == 4 && iax == 5) {/* 4-fold B entry has no E */ ewt0 = -1; } if (isym > 3 && pvinfo->wt[4] < 4) { pvinfo->ewt[0] = pvinfo->ewt[1] = (short) 0; pvinfo->wt[4] |= 4; } if (isym == 6) { /* six-fold, iax = 4,5 */ pvinfo->ewt[iax-4] = (short) ewt0; } else if (ewt0 >= 0) { pvinfo->ewt[0] = pvinfo->ewt[1] = (short) ewt0; } pvinfo->lvupper = 0; if (ewt1 != 0) { if (ewt1 > 10) { pvinfo->lvupper = 2; ewt1 -= 10; } if (isym >= 3 && ewt1 > isym) ewt1 = MOD(ewt1 - 1, isym) + 1; pvinfo->lvupper += (short) (ewt1 << 2); pvinfo->knmin[0] = 0; pvinfo->knmin[3] = 0; /* 0,1,1,0 */ } pvinfo->lvqn = (short) ewt1; for (i = 0; i <= 3; ++i) { k = pvinfo->knmin[i]; if (k < 0 && knnmin == 0) continue; while (k < knnmin) k += 2; pvinfo->knmin[i] = (short) k; } k = getsp(bcdspin, pvinfo); if (nspinqn < k) nspinqn = k; pvinfo->nqn = (short) k; if (ncards == 1) { /* fill in defaults with wt = 0 */ pvinfom = vinfo; for (i = 1; i <= nvibm; ++i) { ++pvinfom; memcpy(pvinfom, vinfo, sizeof(SVIB)); } ivib = -glob.nvib; } setwt(pvinfo, ivib, iax, iwtpl, iwtmn, vsym); /* set up weights */ } while (ivsym < 0); /* end of reading */ setsp(); /* set up format */ if (glob.nqnn == 3) *itd = 3; glob.nqn = glob.nqnn; glob.vibfmt = (nvibm != 0); glob.iqfmt0 = 0; if (glob.vibfmt) { glob.nqn += 1; glob.iqfmt0 = 1000; } nt = glob.nqn; k = nt + nspinqn; glob.nqn = k; if (k > (*nfmt)) k = nt + 2; glob.maxqn = k; glob.nqn0 = nt; glob.iqfmt0 += 100 * nt + k; /* make sure +/- l values come in pairs */ pvinfo = pvinfom = vinfo; k = 0; for (i = 0; i <= nvibm; ++i) { pvinfo->nqn += (short) nt; pvinfo->wt[4] &= 3; j = pvinfo->wt[4]; if (j != 0) { if (EVEN(pvinfo->gsym)) { pvinfo->wt[4] = 0; /* no Itot */ } else if (pvinfo->wt[0] == pvinfo->wt[j] && pvinfo->wt[1] == pvinfo->wt[j ^ 1]) { pvinfo->wt[4] = 0; } } j = pvinfo->lvqn; if (k != 0) { /* previous l was non-zero and first of pair */ isym = (int)((unsigned int)pvinfo->gsym >> 1); if (isym >= 3 && MOD(j + k, isym) != 0) break; if (k > j) { if (j == 0) break; pvinfom->lvqn = (short) (-j); if (pvinfo->knmin[0] == 0) pvinfo->knmin[0] = 2; if (pvinfo->knmin[3] == 0) pvinfo->knmin[3] = 2; } else { pvinfo->lvqn = (short) (-k); if (pvinfom->knmin[0] == 0) pvinfom->knmin[0] = 2; if (pvinfom->knmin[3] == 0) pvinfom->knmin[3] = 2; } for (j = 0; j <= 4; ++j) { if (pvinfo->wt[j] != pvinfom->wt[j]) break; } if (j <= 4) break; if (pvinfo->gsym != pvinfom->gsym) break; if (pvinfo->ewt[0] != pvinfom->ewt[0])break; if (pvinfo->ewt[1] != pvinfom->ewt[1])break; pvinfo->lvupper += 1; j = 0; } k = j; pvinfom = pvinfo; ++pvinfo; } if (k != 0) { puts("L DOUBLETS SHOULD BE IN PAIRS"); exit(EXIT_FAILURE); } k = 0; pvinfo = vinfo; for (ivib = 0; ivib <= nvibm; ++ivib) { nt = pvinfo->nspstat; for (isym = 0; isym < 4; ++isym) { if (getwt(pvinfo, isym, 0, ivwt) <= 0) continue; for (ii = 1; ii <= nt; ++ii) { if (getwt(pvinfo, isym, ii, ivwt) > 0) ++k; } } ++pvinfo; } if (glob.nbkpj > 0) { free(moldv); moldv = NULL; free(blkptr); blkptr = NULL; } glob.nbkpj = k; nl = (size_t) k * sizeof(int); moldv = (int *) mallocq(nl); moldv[0] = 0; nl += sizeof(int); blkptr = (int *) mallocq(nl); blkptr[0] = 0; k = 0; pvinfo = vinfo; for (ivib = 0; ivib <= nvibm; ++ivib) { nt = pvinfo->nspstat; for (isym = 0; isym < 4; ++isym) { if (getwt(pvinfo, isym, 0, ivwt) <= 0) continue; for (ii = 1; ii <= nt; ++ii) { if (getwt(pvinfo, isym, ii, ivwt) > 0) { blkptr[k] = k; moldv[k] = isym + (ivib << 2) + (ii << glob.msshft); ++k; } } } ++pvinfo; } blkptr[k] = k; pvinfo = NULL; *nbcd = (NDECPAR + 1) + glob.vibdec; if (sizeof(int) == 2 && glob.vibdec > 2) { puts("too many vibrations for 16-bit integer computer"); ncards = 0; } *nfmt = 1; if (nspinqn != 0 && nvibm != 0) *nfmt = setfmt(&k, -1); return ncards; } /* setopt */ int setfmt(iqnfmt, nfmt) int *iqnfmt, nfmt; { SVIB *pvinfo; short *jjs; int iqfmtq, iqfmt0, i, j, k, iv, si2, ff, nset, nspinv, nvib; iqfmt0 = 0; nvib = glob.nvib; if (nfmt > 0 && nfmt < nvib) nvib = nfmt; pvinfo = vinfo; for (iv = 0; iv < nvib; ++iv) { iqfmtq = glob.iqfmt0; i = setgsym((int) pvinfo->gsym); jjs = pvinfo->spt; nset = jjs[0]; nspinv = nset - 1; if (nspinv > nspin) nspinv = nspin; jjs += nset; ff = 200 - jjs[0]; if ((int) pvinfo->nqn > glob.maxqn) { iqfmtq += 10 * (ff & 1) + 4000; } else { si2 = 0; k = glob.maxqn - glob.nqn0 - 1; for (i = 0; i < k; ++i) { if (i == itsym && itsym < itptr) { i = itptr; j = 0; si2 = si2 << 1; } if (i < nspinv) { j = jjs[i + 1]; if (i < itsym) j += ff; } else { j = ff; } si2 = (si2 + (j & 1)) << 1; } si2 += (ff & 1); if (si2 > 9) si2 &= 7; if (itptr < nspin) iqfmtq += 2000; if (itsym < itptr) iqfmtq += 4000; iqfmtq += 10 * si2; } if (iv == 0) iqfmt0 = iqfmtq; if (nfmt < 0) { if (iqfmtq != iqfmt0) return nvib; } else { iqnfmt[iv] = iqfmtq; } ++pvinfo; } if (nfmt <= 0) return 1; return glob.maxqn; } int setblk(lu, npar, idpar, par, nbkpf, negy) FILE *lu; const int npar; int *nbkpf, *negy; bcd_t *idpar; const double *par; { /* this subroutine initializes quantum quantities */ /* on entry: */ /* LU= logical unit for listing */ /* IDPAR= parameter id */ /* NBKPF= maximum value of F */ /* NEGY = max dimensioned energy in main program */ /* on return: */ /* NBKPF= blocks per F */ /* NEGY= possibly downsized NEGY to account for local storage */ static char csym[] = { 'x', 'c', 'b', 'a' }; static int nvwt[] = {1,1,1,2,3,2,3}; SVIB *pvinfo; SPAR *spar_now; short *jjs; double ai; size_t nl; int nodd, isym, jsym, ixzt, i, k, kl, lt, n, kdiag, nvibm, neven, nsize, alpha; int kd, ld, ii, jj, im, kk, ni, nj, iv, jv, ivmin, ix, lv, lx, npair, ldel; int knnmin, ntstat, si1, si2, iff, ikq, jjt, ins, neuler, sznz, iff0; int iwt, nsym, gsym, maxblk, maxsblk, itmp, nd, ifc, nf; unsigned int ivsym; int ivwt[3], jvwt[3]; char ctmp; cgetv[0].cblk = cgetv[1].cblk = 0; /* reset store for getqn */ kdiag = glob.idiag; if (pasort(lu, npar, idpar, par) == 0) glob.idiag = -1; if (npar <= 0) { if (glob.nvib > 1) free(vinfo); vinfo = &vinfo1; glob.nvib = 0; if (glob.nbkpj > 0) { free(moldv); moldv = &zmoldv; free(blkptr); blkptr = &zblkptr; glob.nbkpj = 0; } if (glob.maxblk > 0) { free(ivs); ivs = &zivs; free(ikmin); ikmin = &zikmin; free(ibkptr); ibkptr = &zibkptr; glob.maxblk = 0; } if (ndmx > 0) { free(idx); idx = &zidx; free(jdx); jdx = &zjdx; free(wk); wk = &zwk; ndmx = 0; } return 0; } nvibm = glob.nvib - 1; n = (*nbkpf); if (n > 0) { iff0 = n << 1; pvinfo = vinfo; for (i = 0; i <= nvibm; ++i) { jjs = pvinfo->spt; k = jjs[0]; /* select largest N relative to F */ ni = (iff0 - jjs[k]) >> 1; if ((int) pvinfo->knmax > ni) pvinfo->knmax = (short) ni; ++pvinfo; } pvinfo = NULL; } pvinfo = NULL; iff0 = 200; /* find interactions */ ntstat = glob.nbkpj; sznz = ndmax = 0; if (ntstat <= 0) return 0; nf = nspin - 1; for (ii = 0; ii < ntstat; ++ii) { im = moldv[ii]; iff = iff0; kk = getqs(im, iff, -1, -1, ixcom, iscom, &iv); if (ODD(iscom[nspin])) { /* special for odd spin multiplicity (make N integer) */ kk = getqs(im, --iff, 0, -1, ixcom, iscom, &iv); } isym = ixcom[XSYM]; ni = ixcom[XNVAL]; pvinfo = &vinfo[ixcom[XVIB]]; gsym = pvinfo->gsym; nsym = setgsym(gsym); getwt(pvinfo, isym, kk, ivwt); for (jj = 0; jj < ii; ++jj) { jjt = moldv[jj]; kk = getqs(jjt, iff, 0, -1, jxcom, jscom, &jv); jsym = jxcom[XSYM]; nj = jxcom[XNVAL]; nd = ni - nj; pvinfo = &vinfo[jxcom[XVIB]]; if (pvinfo->gsym != (short) gsym) continue; getwt(pvinfo, jsym, kk, jvwt); iwt = nvwt[nsym]; for (k = 0; k < iwt; ++k) { if (ivwt[k] != jvwt[k]) break; } if (k < iwt) continue; /* check for connection restrictions on N */ ixzt = glob.ixz; if (ODD(ixzt) && iscom[nspin] != jscom[nspin]) continue; for (k = 0; k < nf; ++k) { if (k == itsym) k = itptr; ixzt = ixzt >> 1; /* check for matching spin multiplicity*/ if (ODD(iscom[k] + jscom[k])) break; /* check for connection restrictions on spin */ if (ODD(ixzt) && iscom[k] != jscom[k]) break; } if (k < nf) continue; ivmin = (iv < jv) ? iv : jv; itmp = iv + jv + ivmin * glob.vibfac; ivsym = blksym(ixcom, jxcom) + ((unsigned int)itmp << 2); for (spar_now = spar_head[ivmin]; spar_now != NULL; spar_now = spar_now->next) { if (spar_now->ipsym < ivsym) break; if (spar_now->ipsym > ivsym) continue; if (((int) spar_now->flags & MNSQ) != 0) continue; if (sznz < 0) sznzfix(sznz, ni, nj, ixcom, jxcom, iscom, jscom); kl = idpars(spar_now, &ikq, &neuler, <, &ld, &kd, &ins, &si1, &si2, &sznz, &ifc, &alpha, &ldel); if (ODD(neuler)) continue; if (sznz > 0) { sznz = sznzfix(sznz, ni, nj, ixcom, jxcom, iscom, jscom); } if (getmask(ixcom, jxcom, kd, ldel, kl, alpha) == 0) continue; npair = getll(0, ld, lt, kd, si1, si2, lscom, iscom, jscom); if (npair < 0) continue; if (kd > 0) glob.idiag = kdiag; if (nd > ndmax) ndmax = nd; kd = blkptr[ii] - blkptr[jj]; if (kd != 0) { /* II and JJ are coupled */ if (kd > 0) { /* rename ptr II to ptr JJ */ kk = blkptr[ii]; ix = jj; } else { /* rename ptr JJ to ptr II */ kk = blkptr[jj]; ix = ii; } glob.idiag = kdiag; for (k = 0; k <= ii; ++k) { if (blkptr[k] == kk) blkptr[k] = blkptr[ix]; } } break; } /* loop over parameters */ } /* jj loop */ } /* ii loop */ /* block connections now found */ switch (glob.idiag) { case -1: fputs("NO DIAGONALIZATION NEEDED\n", lu); break; case 0: fputs("ENERGY SORT OF WANG SUB-BLOCKS\n", lu); break; case 1: fputs("EIGENVECTOR SORT OF STATES\n", lu); break; case 2: fputs("ENERGY FOLLOWS FIRST ORDER WITHIN WANG SUB-BLOCK\n", lu); break; case 3: fputs("ENERGY FOLLOWS FIRST ORDER WITHIN V AND SPIN SUB-BLOCKS\n", lu); break; } if (glob.oblate) { fputs("OBLATE ROTOR\n", lu); ctmp = csym[1]; csym[1] = csym[3]; csym[3] = ctmp; } else { fputs("PROLATE ROTOR\n", lu); } if (glob.nqnn == 1) { fputs("LINEAR MOLECULE QUANTA, K SUPPRESSED\n", lu); } else if (glob.nqnn == 2) { fputs("SYMMETRIC TOP QUANTA\n", lu); } fputs(" V KMIN KMAX WTPL WTMN ESYMWT NSYM SPINS\n", lu); pvinfo = vinfo; jj = 0; for (iv = 0; iv <= nvibm; ++iv) { knnmin = pvinfo->knmin[0]; ii = pvinfo->knmin[1]; if (ii < knnmin) knnmin = ii; ii = pvinfo->knmin[2]; if (ii < knnmin) knnmin = ii; ii = pvinfo->knmin[3]; if (ii < knnmin) knnmin = ii; if (knnmin < 0) knnmin = 0; k = pvinfo->gsym; isym = k >> 1; ii = (int)(pvinfo->lvupper >> 2); if ((pvinfo->lvupper & 2) != 0) ii += 10; if (ODD(k)) { isym = -isym; jj = 1; } ii = ii * glob.esymdec + pvinfo->ewt[0]; fprintf(lu, " %4d %4d %4d %4d %4d %6d %4d", iv, knnmin, pvinfo->knmax, pvinfo->wt[0], pvinfo->wt[3], ii, isym); jjs = pvinfo->spt; ni = jjs[0]; for (i = 1; i < ni; ++i) { ai = jjs[i] * 0.5; fprintf(lu, "%5.1f", ai); } fputc('\n', lu); if (EVEN(isym)) { /* even */ i = ii; ii -= pvinfo->ewt[0]; if (isym != 4) ii += pvinfo->ewt[1]; if (i != ii || pvinfo->wt[0] != pvinfo->wt[2] || pvinfo->wt[3] != pvinfo->wt[1]) { fprintf(lu, " %4d %4d %4d %4d %4d %6d B\n", iv, knnmin, pvinfo->knmax, pvinfo->wt[2], pvinfo->wt[1], ii); } } ++pvinfo; } if (jj != 0) { fprintf(lu, "I(TOT) IS LAST QUANTUM BEFORE F AND %s\n", "IS INDICATED BY NEGATIVE VALUE OF NSYM ABOVE"); } /* bubble sort block number */ lx = ntstat; while (lx > 1) { lv = 0; im = 0; for (i = 1; i < lx; ++i) { if (blkptr[im] > blkptr[i] || (blkptr[im] == blkptr[i] && moldv[im] > moldv[i])) { lv = i; jj = blkptr[i]; blkptr[i] = blkptr[im]; blkptr[im] = jj; jj = moldv[i]; moldv[i] = moldv[im]; moldv[im] = jj; } im = i; } lx = lv; } fputs("BLOCK - WT - SYM - V - TSP - N - other quanta (rel. to F=0 )\n", lu); /* convert block label to pointer */ i = blkptr[0]; blkptr[0] = 0; n = nsize = neven = nodd = maxblk = maxsblk = 0; for (jj = 0; jj <= ntstat; ++jj) { if (blkptr[jj] != i) { ii = jj - blkptr[n++]; if (maxblk < ii) maxblk = ii; blkptr[n] = jj; if (neven > nsize) nsize = neven; if (nodd > nsize) nsize = nodd; if (jj == ntstat) break; i = blkptr[jj]; neven = 0; nodd = 0; } jjt = moldv[jj]; jjt = getqs(jjt, 0, -1, 0, ixcom, iscom, &ii); isym = ixcom[XSYM]; iv = ixcom[XVIB]; pvinfo = &vinfo[iv]; iwt = getwt(pvinfo, isym, jjt, ivwt); kk = pvinfo->knmin[isym]; if (kk < 0) { ii = ixcom[XISYM]; /* get spin symmetry */ kk = isym; if (is_esym[ii] < 0) kk += 2; kk &= 2; } ii = pvinfo->knmax - kk + 2; if (ii < 0) ii = 0; ii >>= 1; /* size of even block */ neven += ii; if (ii > maxsblk) maxsblk = ii; kk = pvinfo->knmin[3 - isym]; if (kk < 0) { ii = ixcom[XISYM]; /* get spin symmetry */ kk = isym; if (is_esym[ii] >= 0) kk += 2; kk &= 2; } ii = pvinfo->knmax - kk + 2; if (ii < 0) ii = 0; ii >>= 1; /* size of odd block */ nodd += ii; if (ii > maxsblk) maxsblk = ii; --jjt; fprintf(lu, "% 5d %4d %c %4d %4d", n + 1, iwt, csym[isym], iv, jjt); ni = nspin; for (ii = 0; ii < nspin; ++ii) { ai = iscom[ni] * 0.5; fprintf(lu, " %5.1f", ai); if (ni == itsym && itsym < itptr) ii = itptr; ni = ii; } fputc('\n', lu); } ++maxblk; glob.maxblk = maxblk; fprintf(lu, "Maximum Dimension for Hamiltonian = %d\n", nsize); fflush(lu); *nbkpf = n; glob.nbkpj = n; if (nsize < (*negy)) *negy = nsize; if (ndmx > 0) { free(ivs); ivs = NULL; free(ikmin); ikmin = NULL; free(ibkptr); ibkptr = NULL; free(idx); idx = NULL; free(jdx); jdx = NULL; free(wk); wk = NULL; ndmx = 0; } ii = maxsblk + maxsblk - 1; ndmx = nsize; if (ndmx < ii) ndmx = ii; nl = (size_t) (nsize + ndmx) * sizeof(double); wk = (double *) mallocq(nl); wk[0] = zero; /* allocate space for sub-block information */ nl = (size_t) (ndmx + 1) * sizeof(short); idx = (short *) mallocq(nl); idx[0] = 0; jdx = (short *) mallocq(nl); jdx[0] = 0; maxblk = maxblk + maxblk; nl = maxblk * sizeof(short); ibkptr = (short *) mallocq(nl); ibkptr[0] = 0; ikmin = (short *) mallocq(nl); ikmin[0] = 0; nl = (size_t) maxblk * sizeof(int); ivs = (int *) mallocq(nl); ivs[0] = 0; ii = glob.vibfac + 1; return (ii * ii); } /* setblk */ int pasort(lu, npar, idpar, par) FILE *lu; const int npar; bcd_t *idpar; const double *par; { /* initialize spin factors and arrange operators */ /* on entry: */ /* NSIZE= block size and dimension */ /* IDPAR= list of parameter identifiers ( element 0 is length ) */ static int idmy[NDXCOM], initl; static bcd_t idunit[NDECPAR]; SVIB *pvib1, *pvib2; SPAR *spar_now, *spar_last, *spar_match, *spar_free; short *iiv1, *iiv2; double dtmp; size_t nl; int ifac, iv1d, ipar, iv12q, ikqp, insp, njqt, ityi, isym, neulerp, i, lt; int kk, ltp, si1, si2, iv1, iv2, lv1, lv2, kdp, ldp, iv12, ikq, alpha, gsym; int ins, neuler, sznz, sznzp, si1p, si2p, nt, kd, ld, ivdif, ibcd, ldel; int iret, idflags, ifc, ifcp, ilim, kl, klp, nitot, ndecv, nbcd, imag; int alphap, ldelp; unsigned int ivsym; bcd_t *idval; BOOL isused, first; if (glob.parinit > 0) { for (i = 0; i < glob.parinit; ++i) { spar_free = spar_head[i]; while (spar_free != NULL) { spar_now = spar_free->next; free(spar_free); spar_free = spar_now; } spar_head[i] = spar_free; } spar_free = spar_now = NULL; glob.parinit = 0; } if (npar <= 0) { free(ipder); ipder = &zipder; initl = 0; return 0; } /* find reduced spin matrix elements */ spfac[0] = 0.; spfac2[0] = 0.; spfac[1] = sqrt(1.5); spfac2[0] = 0.; nt = glob.mxspin; for (i = 2; i <= nt; ++i) { dtmp = 0.5 * i; /* dtmp = I */ spfac[i] = sqrt(dtmp * (dtmp + 1.) * (i + 1)); spfac2[i] = 0.25 * sqrt((dtmp + 1.5) / (dtmp - 0.5)) / dtmp; } /* initialize SPECFC and DIRCOS */ specfc(0, 0, 0, 0, 0, 0, 0, &zero, &szero, &szero); dircos(idmy, idmy, 0, 0, 0, &zero, &szero, &szero, 0, MODD, 0, &first); ifac = glob.vibfac + 1; ndecv = glob.vibdec; ilim = ifac * ifac - 1; if (initl > 0) { free(ipder); ipder = NULL; } nl = (size_t) npar *sizeof(int); ipder = (int *) mallocq(nl); initl = npar; kk = 1; ityi = 1; ipder[0] = 0; ibcd = 0; nbcd = (int)idpar[0] & 0x7f; for (i = 1; i < npar; ++i) { ibcd += nbcd; if (NEGBCD(idpar[ibcd]) == 0) { ipder[i] = kk++; ityi = i + 1; } else { ipder[i] = -ityi; } } glob.nfit = kk; /* .. code parameter id for power of N*(N+1) and equivalence */ /* .. NJQ = POWER + 1 (first occurrance of power) */ /* .. IPTYP= ITY + 1 (first occurrance of cosine) */ /* .. ITY = sequence number for hybrid operators */ /* .. negate NJQ and IPTYP for successive occurrance */ /* .. IP points to parameter */ glob.parinit = glob.nvib; for (i = 0; i < glob.nvib; ++i) { spar_head[i] = NULL; } pvib1 = pvib2 = vinfo; nsqmax = ityi = iv1d = iret = iv1 = iv2 = nitot = gsym = 0; ipar = npar - 1; ibcd = ipar * nbcd; spar_free = NULL; isused = FALSE; while (ipar >= -1) { ityi = 0; first = TRUE; do { /* repeat until parameter subtypes exhausted */ ivdif = 0; if (ipar < 0) { /* make sure there is a unit operator */ idval = idunit; iv12q = ilim; } else { idval = &idpar[ibcd + 1]; iv12q = bcd2i(idval[0]); if (ndecv > 1) iv12q += bcd2i(idval[1]) * 100; if (ndecv > 2) iv12q += bcd2i(idval[2]) * 10000; idval += ndecv; iv2 = iv12q / ifac; iv1 = iv12q - iv2 * ifac; ivdif = iv1 - iv2; if (ivdif < 0) { i = iv1; iv1 = iv2; iv2 = i; iv12q = iv2 * ifac + iv1; } } if (iv12q == ilim) { iv1 = iv2 = iv1d; ivdif = 0; iv12 = iv1d * (ifac + 1); } else { if (iv1 >= glob.nvib) break; iv12 = iv12q; } if (idval[1] == (bcd_t)0x91 && idval[0] == (bcd_t)0) { /* special for rho */ if (iv1 == iv2) { isused = TRUE; dtmp = par[ipar]; specfc(0, iv1, iv1, 0, 0, 0, 1, &dtmp, &szero, &szero); } break; } pvib1 = &vinfo[iv1]; pvib2 = &vinfo[iv2]; gsym = pvib1->gsym; if ((short) gsym != pvib2->gsym) break; gsym = setgsym(gsym); nitot = glob.nitot; if (spar_free == NULL) { /* allocate more structures */ spar_free = (SPAR *) mallocq(sizeof(SPAR)); } spar_now = spar_free; spar_now->flags = 0; spar_now->alpha = C0; spar_now->mldel = C0; ityi = idpari(idval, ityi, spar_now); if (ityi == 0) break; njqt = (int) spar_now->njq; if (njqt > nsqmax) nsqmax = njqt; isym = (int) spar_now->ipsym; idpars(spar_now, &ikq, &neuler, <, &ld, &kd, &ins, &si1, &si2, &sznz, &ifc, &alpha, &ldel); iiv1 = pvib1->spt; iiv2 = pvib2->spt; if (si1 > 0 && checksp(first, si1, si2, iiv1, iiv2, &spar_now->zfac) != 0) break; if (sznz > 0) { /* setup for SzNz */ if (checksp(first, 1, 0, iiv1, iiv2, &spar_now->zfac) != 0) break; spar_now->zfac *= 0.5; } first = FALSE; /* factor of 2 for B-C, etc */ if (kd != 0 && isym == 0) spar_now->zfac *= 2.; /* setup flags */ idflags = 1; imag = (kd > ld)? kd: ld; if (ifc < 0) ++imag; if (si2 < 0) ++imag; /* commutator is odd */ if (pvib1->wt[4] != pvib2->wt[4]) break; lv1 = pvib1->lvqn; lv2 = pvib2->lvqn; if (lv1 + lv2 < 0) break; if (ODD(pvib1->lvupper)) { --iv12; --iv1; } if (ODD(pvib2->lvupper)) { iv12 -= ifac; --iv2; } /* setup for l-doubling */ ldel = lv1 - lv2; if (ldel != 0) { /* add flag for KL symmetry when off diagonal in l */ if (ldel < 0) ldel = -ldel; if (ivdif < 0) /* DELTA K * DELTA L < 0 */ ldel = -ldel; glob.lsym = FALSE; } else if (lv1 != 0) { /* change symmetry for Lz operator */ i = ((pvib1->lvupper ^ pvib2->lvupper) & 2) >> 1; if (ODD(imag + i)){ isym = 3 - isym; ++imag; idflags |= MLZ; } } spar_now->mldel = (signed char)ldel; if (ODD(imag)) { if (iv1 == iv2 && ldel == 0) break; idflags |= MODD; } if (nitot < 3) { /* calculate alpha */ if (gsym >= 3 && MOD(kd - ldel, gsym) != 0) continue; /* matrix element breaks symmetry */ alpha = 0; } else { alpha = MOD(ldel - kd, nitot); if (si1 <= itsym) { if (alpha != 0) continue; idflags |= MIDEN; /* identity spin operator */ } if (alpha < 0) alpha += nitot; if (alpha == 0 || (alpha + alpha) == nitot) { if (ityi > 160) continue; if (ityi > 80 && kd == 0 && lv1 == 0 && lv2 == 0) continue; } else { glob.esym = FALSE; spar_now->zfac *= 0.5; } if (ityi >= 80) { i = ityi / 80; /* i = 3,2,1,0 */ if (ODD2(i)){ if (ODD(isym + imag)) spar_now->zfac = -spar_now->zfac; idflags |= MSYM2; isym = 3 - isym; } if (ODD(i)){ alpha += nitot; /* special for quasi-diagonal */ } } } if (ikq != 0 || neuler != 0 || sznz != 0 || ifc != 0) idflags |= MNOUNIT; /* no unit matrix */ spar_now->alpha = (unsigned char)alpha; if (isym == 0 && iv1 == iv2 && kd != 0) iret = 1; ivsym = (unsigned int) isym + ((unsigned int)iv12 << 2); spar_now->ipsym = ivsym; spar_match = spar_last = NULL; for (spar_now = spar_head[iv2]; spar_now != NULL; spar_now = spar_now->next) { if (spar_now->ipsym <= ivsym) break; spar_last = spar_now; } if (EVEN(neuler)) { /* not Euler denominator parameter */ kl = idflags & MMASK; for ( /* no init */ ; spar_now != NULL; spar_now = spar_now->next) { if (spar_now->ipsym < ivsym) break; spar_last = spar_now; klp = idpars(spar_now, &ikqp, &neulerp, <p, &ldp, &kdp, &insp, &si1p, &si2p, &sznzp, &ifcp, &alphap, &ldelp); klp &= MMASK; /* check for same K dependence */ if (ikq == ikqp && ld == ldp && kd == kdp && kl == klp && sznz == sznzp && ifc == ifcp && alpha == alphap && ldel == ldelp && neuler == 0 && neulerp == 0) { /* check for same operator except power of N*(N+1) */ if (lt == ltp && ins == insp && si1 == si1p && si2 == si2p) { idflags |= MNSQ; } else { if ((idflags & MNSQ) != 0) break; } spar_match = spar_now; } else { /* operator with different K dependence */ if ((idflags & MNSQ) != 0) break; } } /* end loop over previous parameters */ } isused = TRUE; if (spar_match != NULL) { if (ipar < 0) break; /* already have a unit operator */ idflags |= MCOS_OK; spar_last = spar_match; } spar_now = spar_free; spar_free = NULL; if (spar_last != NULL) { spar_now->next = spar_last->next; spar_last->next = spar_now; } else { spar_now->next = spar_head[iv2]; spar_head[iv2] = spar_now; } spar_last = NULL; spar_now->ip = ipar; spar_now->flags = (short) idflags; } while (ityi > 1); i = ipar; if (iv12q == ilim) { ++iv1d; if (iv1d >= glob.nvib) --ipar; } else { --ipar; } if (i > ipar) { if (!isused) { putbcd(sbcd, NSBCD, &idpar[ibcd]); fprintf(lu, " WARNING: parameter %6d %s has no matrix elements\n", (i + 1), sbcd); } iv1d = 0; ibcd -= nbcd; isused = FALSE; } } /* end ipar loop */ if (spar_free != NULL) free(spar_free); if (glob.esym) glob.lsym = FALSE; return iret; } /* pasort */ int idpari(idval, itp, pspar) bcd_t *idval; int itp; SPAR *pspar; { /* subroutine to parse parameter identifier for complex parameter type */ /* on entry: */ /* IDVAL= parameter ID / vibrational field */ /* ITP= value of sub-parameter type for hybrid operator */ /* on return: */ /* pspar modified */ /*0 1 2 3 4 5 6 7 8 9 */ static int itpnxt[10] = { 0, 0, 0, 2, 6, 6, 2, 0, 9, 0 }; static int itpop[10] = { 0, 0, 0, 10, 40, 40, 10, 1, 0, 1 }; static int zfacv[10] = { 1, 2, 6, 3, -8, 8, -6, 2, 4, -4 }; int ins, ity, nsx, isy, itysav; int itmp, i, iphaz, ldt, lnt, kdt, ksq, si1, si2, exval; bcd_t ibtmp; iphaz = 0; if (itp > 0) { /* get current sub-parameter type from previous */ if (itp < 10) { itp = itpnxt[itp]; if (itp == 0) return 0; } else { iphaz = itp >> 4; itp = itpnxt[itp & 15]; } } pspar->zfac = 1.; pspar->euler = C0; pspar->msznz = C0; pspar->fc = C0; exval = 0; ibtmp = idval[0]; nsx = (int)(ibtmp & 0x0f); ksq = (int)(ibtmp >> 4); ity = bcd2i(idval[1]); itysav = ity * 10 + ksq; ibtmp = idval[2]; ins = (int)(ibtmp & 0x0f); si1 = (int)(ibtmp >> 4); if (ins > 0 && nspin == 0) return 0; if (ins >= 5 && itsym == 0) return 0; ibtmp = idval[3]; si2 = (int)(ibtmp & 0x0f); itmp = (int)(ibtmp >> 4); if (si2 > si1) { /* make sure SI1 > SI2 */ i = si1; si1 = si2; si2 = i; if (si2 == 0) si2 = -1; } if (si1 > nspin) return 0; if (ity > 90) { if (si1 != 0 || ins >= 5) return 0; } ibtmp = idval[4]; if (ibtmp != (bcd_t)0) { if (ibtmp >= (bcd_t) 0x60) return 0; itmp += 10 * (int)(ibtmp & 0x0f); exval = (int) ibtmp & 0xf0; } if (itmp != 0) { i = (itmp - 1) / 10; itmp -= (i >> 1) * 10; if (ODD(i)) itmp = 9 - itmp; pspar->fc = (char) itmp; } pspar->msi1 = (char) si1; pspar->msi2 = (char) si2; pspar->njq = (unsigned char) nsx; /* exclude S * S and I * S from aa, bb, cc operators */ if (itp == 2) { if (si2 > 0) { itp = itpnxt[itp]; if (itp == 0 && iphaz == 0) return 0; } } if (itp != 0) { if (itp == 6 && si1 == 0) { itp = 8; pspar->njq = (unsigned char) (nsx + 1); } } else { if (iphaz != 0) { if (ins >= 5) iphaz -= 1; else if (glob.nitot >= 3) iphaz -= 5; } else { if (ins >= 5) iphaz = 4; /* iphaz = 4,3,2,1,0 */ if (glob.nitot >= 3) iphaz += 15; } if (ity == 0) { itp = 1; } else if (ity <= 3) { itp = ity; if (glob.oblate) itp = revsym[itp]; itp += 2; /* ITP=3,4,5 for A,B,C */ if (itp == 3 && si1 == 0) itp = 7; } else { itp = 1; itmp = ity; if (ity >= 92) itmp = (ity - 92) >> 1; /* negate odd powers of P+**2 */ if (itmp <= 11 && EVEN(itmp)) pspar->zfac = -(pspar->zfac); } } isy = ity / 20; if (isy > 3) { isy = 0; } else if (isy > 0) { if (glob.oblate) { itysav += (revsym[isy] - isy) * 200; } else { isy = revsym[isy]; } } if ((int) pspar->fc < 0) isy = 3 - isy; /* fix for sin coeff. */ pspar->ipsym = (unsigned int) isy; if (itp > 1) { itysav += itpop[itp] - ity * 10; pspar->zfac *= 2. / (double) zfacv[itp]; } itp += iphaz << 4; if (ins >= 5) { ins -= 5; if(iphaz >= 5) iphaz = MOD(iphaz, 5); pspar->msznz = (unsigned char) (iphaz + 1); } pspar->mins = (unsigned char) ins; ity = itysav / 10; ksq = itysav - ity * 10; ldt = kdt = 0; if (ity <= 3) { if (ity != 0) ldt = 2; } else if (ity < 12) { ldt = kdt = ity + ity - 6; } else if (ity < 20) { kdt = ity + ity - 22; ldt = kdt + 1; } else if (ity < 80) { isy = ity / 20; ldt = kdt = ity - isy * 20 + 1; if (ODD((isy >> 1) + ldt)) --kdt; } else if (ity < 90) { /* K energies coded as 8n with K=10*n+KSQ */ if (exval != 0) return 0; ksq = -1 - 10 * (ity - 80) - ksq; } else { if (exval != 0) return 0; i = ity - 90; if (ODD(i)) { if (ksq > 1) { /* augmented Euler series */ --i; ksq += 8; kdt = i; } else { pspar->ipsym = 3; /* symmetry = 3 to put first in list */ } } else { kdt = i; } pspar->euler = (unsigned char) (i + 2); } pspar->ksq = (char) ksq; if (si1 == 0) { lnt = 0; } else if (si2 <= 0) { if (isy == 0 || ldt == 0) si2 = 0; else if (si2 < 0) --ldt; lnt = 1; } else if (ldt < 2) { lnt = 0; } else { lnt = 2; } pspar->mln = (unsigned char) lnt; pspar->mkdel = (unsigned char) kdt; /* try to use raising operator */ if (ldt == kdt && ldt > lnt) ldt -= (ldt - lnt) & 0x7e; /* subtract even number */ pspar->mld = (unsigned char) ldt; if (exval != 0) { /* check for Euler restrictions */ ibtmp = (bcd_t) exval; exval = exval >> 3; if (idval[1] == (bcd_t) 0 && idval[2] == (bcd_t) 0 && idval[3] == (bcd_t) 0 && idval[4] == ibtmp) { /* energy Euler term */ if (idval[0] == (bcd_t) 0) return 0; /* bad energy-like term */ } else if (exval == 2) { exval = 12; /* operator Euler term */ } pspar->euler = (unsigned char) exval; } return itp; } /* idpari */ int checksp(const BOOL first, int si1, int si2, const short *iiv1, const short *iiv2, double *zfac) { int ii, iip; ii = iip = 0; if (si2 > 0) { if ((short) si2 < iiv1[0]) ii = iiv1[si2]; if ((short) si2 < iiv2[0]) iip = iiv2[si2]; if (ii == iip) { if (ii == 0) return 2; if (si2 == si1) { if (ii < 2) return 3; *zfac *= spfac2[ii]; } *zfac *= spfac[ii]; } else { /* ii != iip */ if (ODD(ii + iip)) return 4; /* check multiplicity */ if (si2 > itsym) return 4; } } if (si2 != si1){ ii = iip = 0; if ((short) si1 < iiv1[0]) ii = iiv1[si1]; if ((short) si1 < iiv2[0]) iip = iiv2[si1]; if (ii == iip) { if (ii == 0) return 5; *zfac *= spfac[ii]; } else { /* ii != iip */ if (ODD(ii + iip)) return 6; /* check multiplicity */ if (si1 > itsym) return 6; } } if (first && si1 > itsym){ iip = 0; if (glob.nitot >= 3) iip = si1 - itsym; if (si1 == si2) { if (iip > 1) return 1; setzitot(2, 0, 2, ii, glob.nitot); /* quadrupole */ } else if (si2 > itsym) { if (iip > 2) return 1; setzitot(1, 1, 0, ii, glob.nitot); /* 2-spin product */ setzitot(1, 1, 2, ii, glob.nitot); } else { if (iip > 1) return 1; setzitot(1, 0, 1, ii, glob.nitot); /* 1-spin vector */ } } return 0; } /* checksp */ int setwt(pvinfov, ivib, iax, iwtpl, iwtmn, vsym) SVIB *pvinfov; const int ivib, iax, iwtpl, iwtmn; double vsym; { /* set weights */ /* on entry: */ /* IVIB= number of vibrational states * number of wang states */ /* IAX= axis of symmetry */ /* IWTPL= weight for even rotational states */ /* IWTMN= weight for odd rotational states */ /* VSYM= symmetry code for vibrational states */ /* on return: */ /* IWT= weights for vibration rotation states */ SVIB *pvinfo; int n, nvsym, i, k; short ii, jj, itmp; n = 1; nvsym = 0; if (ivib < 0) { if (vsym > 0.5) nvsym = 1; n = -ivib; } pvinfo = pvinfov; for (i = 0; i < n; ++i) { ii = (short) iwtpl; jj = (short) iwtmn; if (nvsym != 0) { k = (int) (fmod(vsym, 10.) + 0.5); vsym = (vsym - k) / 10.; if (vsym < 0.) nvsym = 0; if (ODD(k)) { itmp = ii; ii = jj; jj = itmp; } } switch (iax) { case 1: /* WEIGHTS FOR Z */ pvinfo->wt[0] = ii; pvinfo->wt[1] = jj; pvinfo->wt[2] = jj; pvinfo->wt[3] = ii; pvinfo->wt[4] = 2; ii = jj; break; case 2: /* WEIGHTS FOR B */ pvinfo->wt[0] = ii; pvinfo->wt[1] = jj; pvinfo->wt[2] = ii; pvinfo->wt[3] = jj; pvinfo->wt[4] = 3; break; case 3: /* WEIGHTS FOR X */ pvinfo->wt[0] = ii; pvinfo->wt[1] = ii; pvinfo->wt[2] = jj; pvinfo->wt[3] = jj; pvinfo->wt[4] = 2; break; case 4: /* A1, A2 only */ pvinfo->wt[0] = ii; pvinfo->wt[3] = jj; pvinfo->wt[4] = 3; break; case 5: /* B1, B2 only */ pvinfo->wt[1] = jj; pvinfo->wt[2] = ii; pvinfo->wt[4] = 3; break; } /* end switch */ ++pvinfo; } return 0; } /* setwt */ int getwt(pvinfo, isym, iispin, ivwt) SVIB *pvinfo; const int isym, iispin; int *ivwt; { int iwt, jsym, ndm, nsym, kodd; short *jjs; if (iispin == 0){ iwt = 0; if (pvinfo->knmax >= pvinfo->knmin[isym] || pvinfo->knmax >= pvinfo->knmin[3 - isym]) iwt = 1; return iwt; } jsym = isym; iwt = pvinfo->wt[4]; nsym = pvinfo->gsym >> 1; if (iwt != 0) { /* resolve Itot symmetries */ jjs = pvinfo->spt; ndm = jjs[0]; if (ODD2(jjs[iispin * ndm + nspin - nsym + 1])) jsym ^= iwt; } ivwt[0] = pvinfo->wt[jsym]; iwt = ivwt[0]; ivwt[1] = 0; ivwt[2] = 0; if (nsym >= 3) { kodd = isoddk[jsym]; if (nsym == 4){ if (kodd != 0) { ivwt[1] = pvinfo->ewt[0]; /* E state for odd K */ iwt = ivwt[1]; } else { ivwt[2] = pvinfo->wt[jsym ^ 2]; /* check B weights */ if (iwt == 0) iwt = ivwt[2]; } } else { ivwt[1] = pvinfo->ewt[kodd]; if (nsym == 6 && kodd != 0) { ivwt[0] = 0; ivwt[2] = pvinfo->wt[jsym ^ 2]; /* check B weights */ iwt = ivwt[2]; } if (iwt == 0) iwt = ivwt[1]; } } return iwt; } /* getwt */ int setgsym(const int gsym) { static int oldgsym = -1; static int nsym; int k, kk; if (gsym == oldgsym) return nsym; oldgsym = gsym; nsym = gsym >> 1; is_esym[0] = 0; for (k = 1, kk = nsym - 1; k < kk; ++k, --kk) { is_esym[k] = 1; is_esym[kk] = -1; } if (k == kk) is_esym[k] = 0; if (ODD(gsym)) { itptr = nspin - 2; itsym = nspin - nsym; glob.nitot = nsym; } else { itptr = itsym = nspin + nspin + 1; glob.nitot = 0; } /* set up nominal spin coupling map */ for (k = 0; k < nspin; ++k) { ismap[k + k] = k - 1; /* last N,J,F1 .. */ if (k < itsym) { ismap[k + k + 1] = k + nspin + 1; /* S,I1,I2 .. */ } else { ismap[k + k + 1] = itptr; /* Itot */ break; } } ismap[0] = nspin; /* fix up position of N */ return nsym; } /* setgsym */ int getsp(ispnx, pvinfo) const bcd_t *ispnx; SVIB *pvinfo; { /* subroutine to set up spin structure */ /* on entry: */ /* ISPNX = spin code */ /* IVIB = vibrational quantum number */ /* on return: */ /* JJS array of spin combinations */ /* (this array is initialized mostly in setsp) */ /* */ SSP *ssp_now; short *jjs; size_t nl; int i, ii, jj, nspinv, idec, nspnx, nsstat, icmp, nset, isbig; int nitot, iret; short iis[MAXSPIN+1]; bcd_t itmp; jj = -1; nspnx = (int)(ispnx[0] & 0x7f); isbig = 0; itmp = ispnx[1]; if ((int)(itmp & 0x0f) == 0) { jj = (int)(itmp >> 4); isbig = 1; } /* decode spins */ nspinv = 0; nl = 1; idec = isbig; for (i = 1; i <= MAXSPIN; ++i) { if (isbig == 0 && jj >= 0) { ii = jj; jj = -1; } else { ++idec; if (idec >= nspnx) break; itmp = ispnx[idec]; ii = (int)(itmp & 0x0f); if (isbig != 0) ii = ii * 10 + jj; jj = (int)(itmp >> 4); } if (ii == 0) break; if (ii >= MAXII) { ii = 1; } else if (ii > 1) { nspinv = i; } nl *= (size_t) ii; --ii; iis[i] = (short) ii; if (ii > glob.mxspin) glob.mxspin = ii; } nsstat = (int) nl; i = (int)(nsstat << glob.msshft); if (i < 0 || (size_t)(i >> glob.msshft) != nl) { puts("spin problem too big"); exit(EXIT_FAILURE); } nitot = 0; iret = nspinv; i = pvinfo->gsym; if (ODD(i)) { nitot = i >> 1; if (nitot > nspinv) { iret += nitot; if (iret > MAXSPIN) iret = MAXSPIN; for (i = nspinv + 1; i <= iret; ++i) iis[i] = 0; nspinv = iret; } if (nitot > 3) { iret -= (nitot - 3); } } pvinfo->nspstat = nsstat; nset = nspinv + 1; iis[0] = (short) nset; ssp_now = &ssp_head; do { /* check previous combinations for match */ icmp = 0; jjs = ssp_now->sspt; for (i = 0; i < nset; ++i) { icmp = iis[i] - jjs[i]; if (icmp != 0) break; } if (icmp == 0) icmp = nitot - ssp_now->nitot; /* match obtained, so return */ if (icmp == 0) { pvinfo->spt = jjs; return iret; } ssp_now = ssp_now->next; } while (ssp_now != NULL); nl = (size_t) (nsstat + 1) * (size_t) nset * sizeof(short); jjs = (short *) mallocq(nl); nl = sizeof(SSP); ssp_now = (SSP *) mallocq(nl); ssp_now->next = ssp_head.next; ssp_head.next = ssp_now; ssp_now->sspt = jjs; ssp_now->ssize = nsstat; ssp_now->nitot = nitot; jjs[0] = (short) nset; for (i = 1; i < nset; ++i) { jjs[i] = iis[i]; } pvinfo->spt = jjs; return iret; } /* getsp */ void setsp(void) { /* subroutine to set up spin structure */ /* let: jjs = array of spin combinations in struct ssp.sspt */ /* jjs[0] = NSET = NSPIN + 1 = size of the set */ /* {N,J,F1..F} or {N,J..Itot,F} */ /* jjs[1]..jjs[NSPIN] = S, I1, I2.. */ /* jjs[k*NSET]..jjs[NSPIN-1+k*NSET] = 2*(N-F),2*(J-F)..*/ /* jjs[NSPIN+k*NSET] = minimum value of 2*F +(1) */ /* */ SSP *ssp_now; ITMIX *pitmix0; EITMIX *pitmix; int nset, nspinv, ispsum, nn, iv, itot, ibgn, iend, ns, ns0, ifhalf; int minf0, minft, jj, ival, ibase, i, itmp, knt, nlim, nitot, ii; int itsym0, itptr0; short *iis, *jjs, *jjs0; nspin = 0; ssp_now = &ssp_head; /* head has no spins */ while ((ssp_now = ssp_now->next) != NULL) { /* loop over spin sets */ iis = ssp_now->sspt; nspinv = iis[0] - 1; for (iv = nspinv; iv > nspin; --iv) { if (iis[0] > 0) { nspin = iv; break; } } } ssp_now = &ssp_head; /* head has no spins */ while ((ssp_now = ssp_now->next) != NULL) { /* loop over spin sets */ iis = ssp_now->sspt; nset = iis[0]; nspinv = nset - 1; if (nspinv > nspin) nspinv = nspin; ispsum = 0; for (iv = 1; iv <= nspinv; ++iv) { ispsum += iis[iv]; } ifhalf = (ispsum & 1); iend = iis[nspinv]; ival = iend + 1; ns = ival; pitmix0 = NULL; pitmix = NULL; nitot = ssp_now->nitot; if (nitot == 0) { nlim = nspinv; ibgn = iend; itsym0 = itptr0 = nspin + nspin + 1; } else { /* Itot coupling */ itptr0 = nspin - 2; itsym0 = nspin - nitot; nlim = itsym0 + 1; ii = iend; for (i = nlim; i < nspinv; ++i) { if ((int) iis[i] != ii) break; ns *= ival; } if (i < nspin) { puts("spins under Itot should be the same"); exit(EXIT_FAILURE); } pitmix0 = get_itmix(ii, nitot); if (nitot < 3) pitmix0 = NULL; iend = nitot * ii; ibgn = (iend & 1); } knt = ssp_now->ssize; ns0 = knt / ns; jjs = iis + nset; for (nn = -ispsum; nn <= ispsum; nn += 2) { /* select N quanta */ minf0 = -ifhalf; if (nn < minf0) minf0 = nn; /* nn = 2 * (N - F) */ if (pitmix0 != NULL) pitmix = pitmix0->mix; for (itot = ibgn; itot <= iend; itot += 2) { ns = ns0 * (itot + 1); if (pitmix != NULL && itot > ibgn) pitmix += 1; for (iv = 0; iv < ns; ++iv) { /* try all possible spins */ jj = nn; jjs[0] = (short) nn; ival = iv; minft = minf0; for (i = 1; i < nlim; ++i) { ibase = iis[i] + 1; itmp = ival; ival /= ibase; itmp -= ival * ibase; jj += itmp - iis[i]; /* jj = J + N - S - 2F */ if (jj < minft) minft = jj; jj += itmp; /* jj = 2 * (J - F) */ jjs[i] = (short) jj; } jj += ival - itot; if (jj < minft) minft = jj; jj += ival; if (jj == 0) { /* this possibility gives F=0 */ if (itptr0 < nspinv) jjs[itptr0 + 1] = (short) itot; jj = ifhalf - minft; minft = jj + (jj & 1); jjs[nspinv] = (short) minft; if (pitmix != NULL){ jjs[itsym0 + 1] = pitmix->qsym[0]; for (i = itsym0 + 2; i <= itptr0; ++i) jjs[i] = 0; nitot = pitmix->n; jjs0 = jjs; for (i = 1; i < nitot; ++i) { jjs += nset; --knt; for (jj = 0; jj <= nspinv; ++jj) jjs[jj] = jjs0[jj]; jjs[itsym0 + 1] = pitmix->qsym[i]; } } if (--knt <= 0) { ns0 = 0; break; } jjs += nset; } } /* end loop over all spins */ } /* end loop over Itot */ } /* end loop over N */ } /* end loop over spin types */ } /* setsp */ int getqq(iblk, f, iwtb, sbkptr, kmin, vs) const int iblk; int *f, *iwtb, *vs; short *sbkptr, *kmin; { /* subroutine to get sub-block structure and major quantum numbers */ /* on entry: */ /* MXBLK=biggest block size */ /* IBLK= block number */ /* on return: */ /* GETQQ= number of sub-blocks */ /* F= F quantum number *2 */ /* IWTB= statistical weight */ /* SBKPTR= pointer to beginning index for sub-block */ /* KMIN= minimum K for each sub-block */ /* VS= vibration and symmetry for each sub-block */ SVIB *pvinfo; short *jjs; int ibgn, iend, ksym, i, n, nsize, isblk, nsblk, k, kk, nn, iv; int mss, nset, ff; unsigned int mvs; isblk = iblk - 1; if (glob.nbkpj <= 0) return 0; ff = isblk / glob.nbkpj; isblk -= glob.nbkpj * ff; /* form remainder */ ff = ff + ff; /* get information for sub-blocks */ nsize = nsblk = 0; ibgn = blkptr[isblk]; iend = blkptr[isblk + 1]; for (i = ibgn; i < iend; ++i) { mvs = (unsigned int) moldv[i]; ksym = (int) mvs & 3; iv = (int) (mvs >> 2) & glob.msmask; mss = (int) (mvs >> glob.msshft); pvinfo = &vinfo[iv]; if (i == ibgn) setgsym((int)pvinfo->gsym); jjs = pvinfo->spt; nset = jjs[0]; jjs += nset * mss; if (ff < jjs[nset - 1]) continue; nn = jjs[0] + ff; n = nn >> 1; if (ODD(n)) ksym = 3 - ksym; k = pvinfo->knmin[ksym]; if (k < 0) { kk = MOD(jjs[itsym + 1] >> 2, glob.nitot); /* get spin symmetry */ k = ksym; if (is_esym[kk] < 0) k += 2; k &= 2; } kk = pvinfo->knmax; if (kk > n) kk = n; kk -= k; if (kk >= 0) { sbkptr[nsblk] = (short) nsize; kmin[nsblk] = (short) k; vs[nsblk] = (int) mvs; nsize += 1 + (kk >> 1); if (nsblk == 0) { *f = ff - (nn & 1); k = pvinfo->wt[4]; if (k != 0 && ODD2(jjs[itsym + 1])) mvs ^= k; mvs &= 3; iwtb[0] = pvinfo->wt[mvs]; iwtb[1] = pvinfo->ewt[isoddk[mvs]]; iwtb[2] = pvinfo->wt[mvs ^ 2]; iwtb[3] = pvinfo->gsym; iwtb[4] = pvinfo->nqn; } ++nsblk; } } if (nsblk > 0) sbkptr[nsblk] = (short) nsize; return nsblk; } /* getqq */ int getqn(iblk, indx, maxqn, iqn, idgn) const int iblk, indx, maxqn; short *iqn; int *idgn; { /* subroutine to get quantum numbers */ /* on entry: */ /* IBLK= block number */ /* INDX= index within block */ /* on return: */ /* IQN= list of quantum numbers */ /* IDGN= degeneracy */ /* NCOD= seraching aid to find state within a Wang block */ static int last; int ibgn, iend, iqsp, isym, i, k, nsblk, ioff, lower, iv, iqf; int n, ncod, ldgn, ix, kq, ivbase, lv; short mgsym, ngsym; /* get sub-block info from store or GETQQ */ /* look at last used then oldest */ if (iblk == cgetv[0].cblk) { cgetq = cgetv; last = 0; ioff = 0; nsblk = cgetq->cnblk; } else if (iblk == cgetv[1].cblk) { cgetq = &cgetv[1]; last = 1; ioff = glob.maxblk; nsblk = cgetq->cnblk; } else if (last != 0) { cgetq = cgetv; nsblk = getqq(iblk, &cgetq->cff, cgetq->cwt, ibkptr, ikmin, ivs); if (nsblk == 0) { *idgn = 0; return glob.nqn; } last = 0; ioff = 0; cgetq->cnblk = nsblk; cgetq->cblk = iblk; cgetq->csblk = 0; } else { cgetq = &cgetv[1]; ioff = glob.maxblk; nsblk = getqq(iblk, &cgetq->cff, cgetq->cwt, &ibkptr[ioff], &ikmin[ioff], &ivs[ioff]); if (nsblk == 0) { *idgn = 0; return glob.nqn; } last = 1; cgetq->cnblk = nsblk; cgetq->cblk = iblk; cgetq->csblk = ioff; } /* check for request of size */ if (indx <= 0) { *idgn = ibkptr[nsblk + ioff]; return glob.nqn; } /* search for sub-block */ ix = indx - 1; ibgn = cgetq->csblk; iend = nsblk + ioff; if (ix < ibkptr[ibgn]) ibgn = ioff; for (i = ibgn + 1; i < iend; ++i) { if (ix < ibkptr[i]) break; } cgetq->csblk = ibgn = i - 1; /* assemble quanta */ ncod = ibkptr[i] - ibkptr[ibgn]; ldgn = cgetq->cwt[0]; ngsym = (short) setgsym(cgetq->cwt[3]); iqf = cgetq->cff; iqsp = getqs(ivs[ibgn], iqf, 0, 0, ixcom, iscom, &ivbase) - 1; iv = ixcom[XVIB]; isym = ixcom[XSYM]; n = ixcom[XNVAL]; k = ix - ibkptr[ibgn]; kq = ikmin[ibgn] + (k << 1); if (ngsym >= 3) { mgsym = (short)MOD(kq - ixcom[XLVAL] + ixcom[XISYM] + ngsym, ngsym); if (mgsym != 0) { if ((mgsym + mgsym) == ngsym) { /* B symmetry */ ldgn = cgetq->cwt[2]; } else { /* E symmetry */ if (glob.esym) isym = 4; ldgn = cgetq->cwt[1]; } } } lv = ixcom[XLVAL]; if (lv != 0) { /* l-doubled state */ if (glob.lsym) isym = 4; ncod = 1; if (kq == 0 && isym == 3) { if (iv == ivbase) ++iv; else --iv; } } if (glob.nqnn == 2) { if (ODD(isym)) kq = -kq; ncod = 1; } else { if (isym <= 3) lower = (n + isym) & 1; else /* degenerate state */ lower = (lv > 0)? 1: 0; if (glob.oblate) { iqn[2] = (short) kq; kq = n - kq + lower; ncod = -ncod; } else { iqn[2] = (short) (n - kq + lower); } } iqn[0] = (short) n; iqn[1] = (short) kq; k = glob.nqnn; if (glob.vibfmt) iqn[k++] = (short) iv; if (cgetq->cwt[4] > maxqn) { iqn[k++] = (short) iqsp; iqn[k] = (short) ((iqf + 1) >> 1); } else { for (i = 0; k < maxqn; ++i) { iqn[k++] = (short) ((iscom[i] + 1) >> 1); if (i == itsym && ngsym > 3) i += ngsym - 3; } } ldgn *= iqf + 1; if (ldgn < 0) ldgn = 0; *idgn = ldgn; return ncod; } /* getqn */ int dircos(xbra, xket, ld, kd, ncmax, direl, ibra, iket, ifup, loff, mask, isunit) const int *xbra, *xket; const int ld, kd, ncmax, ifup, loff; int mask; double *direl; short *ibra, *iket; BOOL *isunit; { /* SUBROUTINE TO CALCULATE DIRECTION COSINE ELEMENTS */ /* XBRA,XKET INTEGER VECTORS CONTAINING: */ /* 0: DIMENSION OF SUB-BLOCK */ /* 1: SYMMETRY CODE FOR BLOCK (0=A,1=Bx,2=By,3=Bz) */ /* 2: 'N' QUANTUM NUMBER */ /* 3: BEGINNING K */ /* 4: L VALUE */ /* 5: VIBRATIONAL QUANTA */ /* 6: SPIN SYMMETRY */ /* L IS TENSOR ORDER OF COSINE */ /* KD IS TOTAL K QUANTIUM CHANGE */ /* DIREL IS A VECTOR OF DIRECTION COSINE ELEMENTS */ /* NCMAX IS MAXIMUM LENGTH OF DIREL */ /* IBRA,IKET INDICES OF MATRIX ELEMENTS IN DIREL (START WITH ZERO) */ /* IFUP =0 IF NO UPPER TRIANGLE TO BE CALC. */ /* (LOFF AND MSYM2) IF OPERATOR HAS ALTERNATE SYMMETRY UNDER ITOT */ /* (LOFF AND MLZ) IF OPERATOR INCLUDES LZ */ /* (LOFF AND MODD) IF OPERATOR IS ODD-ORDER */ /* L=0 KD=0 DIREL=1. */ /* IFUP <0 IF ELECTRIC DIPOLE */ /* ODD ORDER PARAMETERS ARE IMAGINARY AND */ /* EVEN ARE REAL IF NOT ELECTRIC DIPOLE */ /* FOR ELECTRIC DIPOLE ODD ORDER IS REAL AND EVEN ARE IMAGINARY */ /* L=1 KD=0 DIREL=PHI(Z) */ /* L=1 KD=1 DIREL=PHI(X) OR PHI(Y) */ /* L=2 KD=0 DIREL=0.5*(3* PHI(Z)**2 -1.) */ /* L=2 KD=1 DIREL=PHI(Z)*PHI(X)+PHI(X)*PHI(Z) OR */ /* =PHI(Z)*PHI(Y)+PHI(Y)*PHI(Z) */ /* L=2 KD=2 DIREL=2*(PHI(X)**2 - PHI(Y)**2) OR */ /* =PHI(X)*PHI(Y)+PHI(Y)*PHI(X) */ /* IF KD > L THEN (P-)**(KD-L) OPERATOR MULTIPLIED BY */ /* L,L OPERATOR */ /* WANG FUNCTIONS MULTIPLIED BY i**(SYMMETRY) TO MAKE ALL */ /* ODD ANGULAR MOMENTA IMAG. AND ALL EVEN ONES REAL */ /* IFC = FOURIER SERIES, > 0 IS COS, < 0 IS SIN */ /* MASK = 3-BIT MASK TO ENABLE MATRIX ELEMENT SUB-TYPE */ /* BEST RESULTS ON SPEED IF N FOR BRA CHANGES SLOWEST */ /***********************************************************************/ #define NFAC_DIRCOS 20 static double ff[MAXN_DIRCOS], fac[NFAC_DIRCOS]; double ele, elex, eletmp; int kdel, kbra, kbgn, nbgn, kket, isgn; int ncase, isum, kbra0, kdel2, kket0, kket2, i, k; int n, ikdel, kkdel, isbra, nnbra, nsbra, isket, idif; int ii, il, kk, lld, ir, ks, nnket, nsket, isrev; int iphaz, isym, nbra, nket, kbra0x, kket0x; BOOL matsym; kbra0 = xbra[XKBGN]; kbra0x = xbra[XIQN] & 1; kket0 = xket[XKBGN]; kket0x = xket[XIQN] & 1; nsbra = xbra[XDIM]; isbra = xbra[XSYM]; nsket = xket[XDIM]; isket = xket[XSYM]; isym = (isbra ^ isket) & 3; if ((loff & MODD) != 0) { if (ncmax == 0) { /* initialize */ fac[0] = 0.5; for (lld = 1; lld < NFAC_DIRCOS; ++lld) { ele = lld - 0.5; fac[lld] = fac[lld - 1] * sqrt(lld / ele); } fac[0] = 1. / fac[1]; return 0; }; isym = 3 - isym; } /****** set up state phases, isum is power of i */ isrev = isym >> 1; /* isym is odd if operator is imaginary */ isum = 4 + isket - isbra + (isym & 1); if (ifup < 0) /* correct to make dipoles real */ ++isum; if (ODD(isum) && xbra[XVIB] < xket[XVIB]) isum += 2; /* imaginary vibration phase assumed */ isgn = isum >> 1; if ((loff & MLZ) != 0 && xket[XLVAL] < 0) ++isgn; /* Lz operator found */ ncase = 0; /*********** TRY SCALAR ******************************/ if (kd == 0 && ld == 0 && (mask & 5) != 0) { mask &= 2; ii = 0; ncase = nsket; n = nsbra; kdel = kbra0 - kket0; kk = kbra0; if (kdel != 0) { if (kdel > 0) { kdel = kdel >> 1; ncase -= kdel; if (ncase < 0) ncase = 0; } else { ii = (-kdel) >> 1; kdel = -ii; n -= ii; if (n < 0) n = 0; kk = kket0; } } if (ncase > n) ncase = n; if (ncase > 0) { idif = ncmax - ncase; if (idif < 0) return idif; if (kk == 0) { /* kk = max beginning k */ if (kbra0x == kket0x) { if (kbra0x == 0) mask = 0; *isunit = TRUE; kk = -1; } else { *isunit = FALSE; mask = 0; /* kk = 0 */ } } else { /* kk > 0 */ *isunit = TRUE; mask = 0; } ele = 1.; if (ODD(isgn+isrev)) ele = -ele; for (k = 0; k < ncase; ++k) { direl[k] = ele; ibra[k] = (short) ii; iket[k] = (short) (ii + kdel); ++ii; } if (kk == 0) direl[0] *= fac[0]; } if (mask == 0) return ncase; /* no quasi-diagonal delta K = 0 */ } /************** end Scalar *********/ *isunit = FALSE; nbra = xbra[XNVAL]; nket = xket[XNVAL]; kdel = kd; ikdel = kdel - ld; if (ikdel <= 0) { kkdel = kdel; } else { k = kbra0 + ikdel + ((nsbra - 1) << 1); ffcal(nbra, k, ff); kkdel = ld; } kdel2 = kkdel + kkdel; lld = ld + ld; nnbra = nbra + nbra; nnket = nket + nket; if (ld > 0) { ele = (double) (nnbra + 1); if (nnket != nnbra) ele = sqrt(ele * (nnket + 1)); if (kkdel != 0) ele *= fac[ld]; isgn += nbra + kket0 + kkdel; } else { ele = 1.; if (kdel != 0) ele = 0.5; } matsym = TRUE; k = mask & 5; if (kdel == 0){ if (k != 0) mask = (mask & 2) + 1; matsym = FALSE; } else { if (k == 1 || k == 4) matsym = FALSE; } if (matsym) { for (k = 1; k <= 5; ++k) { if (xbra[k] != xket[k]) { matsym = FALSE; break; } } } if (ODD(isgn + isrev)) ele = -ele; for (iphaz = 1; iphaz <= 4; iphaz = iphaz << 1) { /* iphaz = 1,2,4 */ if (mask < iphaz) break; if (iphaz != 2) { /* iphaz = 1: +KDEL, iphaz = 4: -KDEL */ /* START CALCULATION OF */ /* KBRA=KBRA0+ (IL+ i)*KINC IL >= 0 */ /* KKET=KKET0+ (IR+ i)*KINC IR >= 0, i=0,... */ /* KBRA=KKET+KDEL */ if ((mask & iphaz) == 0) continue; if (iphaz == 4) { if (matsym) break; kdel = -kdel; kkdel = -kkdel; } isum = kdel + kket0 - kbra0; /* ISUM = (IL-IR)*KINC */ if (isum >= 0) { il = isum >> 1; ir = 0; kbra = kdel + kket0; kket = kket0; } else { il = 0; ir = (-isum) >> 1; kbra = kbra0; kket = kbra0 - kdel; } n = nsket - ir; i = nsbra - il; if (n > i) n = i; if (n <= 0) continue; nbgn = ncase; ncase += n; idif = ncmax - ncase; if (idif < 0) return idif; kket2 = kket + kket; if (ld == 0) { for (n = nbgn; n < ncase; ++n) { direl[n] = ele; ibra[n] = (short) il; iket[n] = (short) ir; ++il; ++ir; } } else { for (n = nbgn; n < ncase; ++n) { kk = -kdel2 - kket2; direl[n] = ele * c3jj(nnbra, lld, nnket, kk, kdel2, kket2); kket2 += 4; ibra[n] = (short) il; iket[n] = (short) ir; ++il; ++ir; } } il = ir = 1; if (kbra == 0) il = kbra0x; if (kket == 0) ir = kket0x; if (ODD(il + ir)) /* CORRECT FOR K=0 */ direl[nbgn] *= fac[0]; if (ikdel > 0) { /* ADD ON P- OR P+ OPERATOR TO LEFT */ kbgn = kbra; if (kdel > 0) kbgn -= ikdel; for (n = nbgn; n < ncase; ++n) { eletmp = direl[n]; kk = kbgn + ikdel; for (k = kbgn; k < kk; ++k) { eletmp *= ff[k]; } direl[n] = eletmp; kbgn += 2; } } } else { /* iphaz == 2 */ /* DO QUASI-DIAGONAL PART */ /* K=1,...,KDEL-1 */ /* KBRA=KBRA0+ (IL - i)*KINC IL >= n-1 */ /* KKET=KKET0+ (IR + i)*KINC IR >= 0 , i=0 .. n-1 */ /* KDEL=KKET+KBRA */ kdel2 = -kdel2; if (ODD(isrev)) ele = -ele; if ((mask & 2) == 0) continue; isum = kdel - kbra0 - kket0; if (isum < 0) continue; isum = isum >> 1; /* isum = (IL + IR) */ ir = 0; /* IR = 0 */ if (isum >= nsbra) ir = isum - nsbra + 1; else if (kket0x == 0) ++ir; n = isum + 1; if (n > nsket) n = nsket; else if (kbra0x == 0) --n; n -= ir; if (n <= 0) continue; il = isum - ir; if (matsym) n = (il - ir + 2) >> 1; kket = kket0 + (ir << 1); kbra = kdel - kket; nbgn = ncase; ncase += n; idif = ncmax - ncase; if (idif < 0) return idif; if (ODD2(isbra + nnbra)) { elex = -ele; } else { elex = ele; } kket2 = kket + kket; for (n = nbgn; n < ncase; ++n) { ibra[n] = (short) il; iket[n] = (short) ir; eletmp = elex; if (ld != 0) { ks = -kdel2 - kket2; eletmp *= c3jj(nnbra, lld, nnket, ks, kdel2, kket2); kket2 += 4; } for (kk = 0; kk < ikdel; ++kk) { k = kbra - kk; if (k <= 0) { eletmp *= ff[-k]; } else { eletmp *= ff[k - 1]; } } direl[n] = eletmp; --il; ++ir; kbra -= 2; } /* loop over elements */ } } /* loop over 3 cases */ if (matsym && ncase > 0 && ifup != 0) { n = ncase - 1; if (ibra[n] != iket[n]) ++n; k = ncase; ncase += n; idif = ncmax - ncase; if (idif < 0) return idif; for (i = 0; i < n; ++i, ++k) { direl[k] = direl[i]; ibra[k] = iket[i]; iket[k] = ibra[i]; } } return ncase; } /* dircos */ int ffcal(nff, kff, ff) const int nff, kff; double *ff; { static int nlast = -1; static int klast = 0; static double sq; int k; if (nff != nlast) { nlast = nff; sq = (double) nff; sq *= nff + 1; k = 0; } else { if (klast >= kff) return 0; k = klast; } klast = kff; if (klast > nff) klast = nff; while (k < klast) { ff[k] = sqrt(sq); ++k; sq -= k + k; } return 0; } /* ffcal */