From 7784c0c0dff8d4232cec86d1a77eecf4312d879f Mon Sep 17 00:00:00 2001 From: Robert Harrison Date: Thu, 4 May 2000 21:18:40 +0000 Subject: [PATCH] symmetry orbital infor for abelian groups --- src/symmetry/sym_abelian.F | 430 ++++++++++++++++++++++++++++++++++++- 1 file changed, 429 insertions(+), 1 deletion(-) diff --git a/src/symmetry/sym_abelian.F b/src/symmetry/sym_abelian.F index 8fcf0f3967..94c93ce421 100644 --- a/src/symmetry/sym_abelian.F +++ b/src/symmetry/sym_abelian.F @@ -1,6 +1,6 @@ logical function sym_abelian_group(geom) * -* $Id: sym_abelian.F,v 1.2 1999-07-27 21:06:58 d3e129 Exp $ +* $Id: sym_abelian.F,v 1.3 2000-05-04 21:18:40 d3g681 Exp $ * implicit none integer geom ! [input] @@ -19,3 +19,431 @@ c sym_abelian_group = inp_match(nab,.false.,group,ab,ind) c end + subroutine sym_abelian_bas_info(basis, + $ sh_op_map, bf_op_map, bf_phase, + $ sh_n_uq_op, bf_n_uq_op, + $ sh_uq_op, bf_uq_op, + $ sh_n_uq, bf_n_uq, + $ sh_uq, bf_uq, + $ char_tab, + $ bf_per_ir, + $ bf_per_ir_cum, + $ bf_so_ir, + $ zir, + $ so_uq_to_lab, so_lab1, so_lab2) + implicit none +#include "bas.fh" +#include "mafdecls.fh" +#include "sym.fh" + integer basis ! [input] + integer sh_op_map(8,*), bf_op_map(8,*) ! All other arguments are [output] + double precision bf_phase(8,*) + integer sh_n_uq_op(*), bf_n_uq_op(*) + integer sh_uq_op(8,*), bf_uq_op(8,*) + integer sh_n_uq, bf_n_uq + integer sh_uq(*), bf_uq(*) + double precision char_tab(8,8) + integer bf_per_ir(8) + integer bf_per_ir_cum(8) + integer bf_so_ir(8,*) + character*8 zir(8) + integer so_uq_to_lab(2,8,*) + integer so_lab1(3,*) + integer so_lab2(8,*) +c +c Returns: +c +c sh_op_map(R,ish) the result of applying operator R to shell ish +c . R|ish> -> |jsh> +c . sh_op_map(R,ish) = jsh +c +c bf_op_map(R,i) and bf_phase(R,i) are the result of applying operator R +c to basis function i +c . R|i> -> (+/-)|j> +c . bf_op_map(R,i) = j +c . bf_phase(R,i) = (+/-)1 +c +c sh_n_uq_op(ish) is the number of operations that produce +c unique results when applied to shell ish. +c +c bf_n_uq_op(i) is the number of operations that produce unqiue +c results when applied to basis function i (same as that of ish). +c This is also the number of irreps that can be generated by +c applying projectors to this function (1 per irrep). +c +c sh_uq_op(Ruq,ish) is the list of operators that produce unique +c results for shell ish +c . sh_uq_op(Ruq=1..sh_n_uq_op(ish),ish) = R +c +c bf_uq_op(Ruq,i) is the list of operators that produce unique +c results for basis function i (same as that of ish) +c . bf_uq_op(Ruq=1..bf_n_uq_op(i),i) = R +c +c sh_n_uq is the number of symmetry unique shells +c +c bf_n_uq is the number of symmetry unique basis functions +c +c sh_uq(ish_uq=1..sh_n_uq) is the ish_uq'th symmetry unique shell +c (lexically highest symmetry related shell) +c +c bf_uq(iuq=1..bf_n_uq) is the iuq'th symmetry unique basis function +c (lexically highest symmetry related basis function) +c +c char_tab(op,ir) is the character of the operation op for irreducible +c representation ir (+/- 1). +c +c bf_per_ir(ir) is the no. of basis functions in this irrep +c +c bf_per_ir_cum(ir) is the sum of all basis function in irreps. +c lower than ir +c . bf_per_ir_cum(ir) = sum(k=1..ir-1) bf_per_ir(k) +c +c bf_so_ir(iruq=1..bf_n_uq_op(i),i) is the list of irreps that +c can be generated from function i (with 1 function per irrep.) +c +c zir(ir) returns the name of the irrep +c +c (see below for decsription of the labelling of SOs) +c so_uq_to_lab(1,ir,iuq) = lab1 = SO label for projector for ir +c . acting on the iuq'th unique basis function +c so_uq_to_lab(2,ir,iuq) = lab2 = second type of SO label +c +c so_lab1(1,lab1) = irrep of given SO +c so_lab1(2,lab1) = lab2 = mapping between labelling schemes +c so_lab1(3,lab1) = iuq = unique basis function that is generator +c +c so_lab2(ir,lab2) = lab1 = mapping between labelling schemes +c +c Unlike the rest of the symmetry code this routine INCLUDES +c the identity as the first operator, whereas the rest of +c the symmetry code does not include the identity. +c + integer geom, nop, l_a, k_a, l_b, k_b, nbf, nsh + integer i, j, ish, jsh, ir, op, opuq, ilo, ihi, iuq + integer lab1, lab2(8) + integer nbfir + double precision c(8), val + double precision ddot + external ddot +c + if (.not. bas_geom(basis,geom)) + $ call errquit('sym_abelian_bas_info: basis?', basis) + if (.not. bas_numbf(basis,nbf)) + $ call errquit('sym_abelian_bas_info: basis?', basis) + if (.not. bas_numcont(basis,nsh)) + $ call errquit('sym_abelian_bas_info: basis?', basis) +c + if (.not. sym_abelian_group(geom)) + $ call errquit('the point group is not abelian',0) +c + nop = sym_number_ops(geom) + 1 + if (.not. ma_push_get(mt_dbl, nbf, 'abel a', l_a, k_a)) + $ call errquit('sym_abel_bas_info: ma?', nbf) + if (.not. ma_push_get(mt_dbl, nbf, 'abel b', l_b, k_b)) + $ call errquit('sym_abel_bas_info: ma?', nbf) +c +c Construct the effect of the group operations on shells and +c basis functions using the existing routine to apply an +c operator to an MO vector ... brute force but correct. +c + do ish = 1, nsh + sh_op_map(1,ish) = ish ! identity + if (.not. bas_cn2bfr(basis, ish, ilo, ihi)) + $ call errquit('sym_abelian_bas_info: basis?', basis) + + do i = ilo, ihi + bf_op_map(1,i) = i ! identity + bf_phase(1,i) = 1.0d0 + call dfill(nbf, 0.0d0, dbl_mb(k_a), 1) + dbl_mb(k_a+i-1) = 1.0d0 + do op = 2,nop + call dfill(nbf, 0.0d0, dbl_mb(k_b), 1) +* write(6,*) ' BF ', i, ' OP ', op +* call output(dbl_mb(k_a), 1, nbf, 1, 1, nbf, 1, 1) + call sym_movecs_apply_op(basis, op-1, + $ dbl_mb(k_a), dbl_mb(k_b)) +* call output(dbl_mb(k_b), 1, nbf, 1, 1, nbf, 1, 1) + do j = 1, nbf + if (abs(abs(dbl_mb(k_b+j-1))-1.0d0) .lt. 1d-12) then +* write(6,*) ' J ', j + bf_op_map(op,i) = j + bf_phase(op,i) = dbl_mb(k_b+j-1) + goto 10 + end if + end do + call errquit('sym_abelian_bas_info: _op_map failed',i) + 10 continue + if (.not. bas_bf2cn(basis, j, jsh)) + $ call errquit('sym_abelian_bas_info: basis?', basis) +* write(6,*) ' JSH ', jsh + sh_op_map(op,ish) = jsh + end do + end do + end do +c + if (.not. ma_chop_stack(l_a)) call errquit + $ ('sym_abelian: ma corrupted', 0) +c + write(6,*) ' Abelian group shell symmetry _op_map' + do ish = 1, nsh + write(6,11) ish,(sh_op_map(op,ish),op=1,nop) + 11 format(i5,3x,8i5) + end do + write(6,*) ' Abelian group basis function symmetry _op_map' + do i = 1, nbf + write(6,11) i,(nint(bf_op_map(op,i)*bf_phase(op,i)),op=1,nop) + end do +c +c Determine for each shell/bf the operations that generate +c all symmetry related shells. +c + call ifill(nsh, 0, sh_n_uq_op, 1) + call ifill(nbf, 0, bf_n_uq_op, 1) + do ish = 1, nsh + do op = 1, nop + do opuq = 1, sh_n_uq_op(ish) + if (sh_op_map(sh_uq_op(opuq,ish),ish) .eq. + $ sh_op_map(op,ish)) goto 20 + end do + sh_n_uq_op(ish) = sh_n_uq_op(ish) + 1 + sh_uq_op(sh_n_uq_op(ish),ish) = op + if (.not. bas_cn2bfr(basis, ish, ilo, ihi)) + $ call errquit('sym_abelian_bas_info: basis?', basis) + do i = ilo, ihi + bf_n_uq_op(i) = bf_n_uq_op(i) + 1 + bf_uq_op(bf_n_uq_op(i),i) = op + end do + 20 continue + end do + end do +c + write(6,*) ' Abelian symmetry unique operations by shell' + do ish = 1, nsh + write(6,22) ish, (sh_uq_op(opuq,ish), opuq=1,sh_n_uq_op(ish)) + 22 format(i5,2x,8i3) + end do + write(6,*) ' Abelian symmetry unique operations by function' + do i = 1, nbf + write(6,22) i, (bf_uq_op(opuq,i), opuq=1,bf_n_uq_op(i)) + end do +c +c Determine the lexically highest symmetry related shells/bf +c + sh_n_uq = 0 + bf_n_uq = 0 + do ish = 1, nsh + do op = 1, nop + if (ish .lt. sh_op_map(op,ish)) goto 30 + end do + sh_n_uq = sh_n_uq + 1 + sh_uq(sh_n_uq) = ish + if (.not. bas_cn2bfr(basis, ish, ilo, ihi)) + $ call errquit('sym_abelian_bas_info: basis?', basis) + do i = ilo, ihi + bf_n_uq = bf_n_uq + 1 + bf_uq(bf_n_uq) = i + end do + 30 continue + end do +c + write(6,*) ' Abelian symmetry unique shells', sh_n_uq + write(6,33) (sh_uq(ish),ish=1,sh_n_uq) + 33 format(15i5) + write(6,*) ' Abelian symmetry unique functions', bf_n_uq + write(6,33) (bf_uq(i),i=1,bf_n_uq) +c +c Symmetry orbitals (SO) are formed by applying a projector +c to the generating functions (the symmetry unique or lexically +c highest symmetry related functions). +c +c The number of symmetry related functions (bf_n_uq_op(i)) +c is also the number of irreducible representations spanned +c by the symmetry related functions (Abelian groups only). +c + call sym_abelian_char_tab(geom, char_tab, zir) + call ifill(8, 0, bf_per_ir, 1) +c +c Loop thru generators, compute the character +c and hence the irreps spanned. +c + do iuq = 1, bf_n_uq ! Loop over generators + i = bf_uq(iuq) +c + do op = 1, nop + if (bf_op_map(op,i) .eq. i) then + c(op) = bf_n_uq_op(i)*bf_phase(op,i) + else + c(op) = 0.0d0 + endif + end do +c + nbfir = 0 + do ir = 1, nop + val = ddot(nop, c, 1, char_tab(1,ir), 1) / dble(nop) + if (abs(val-1.0d0) .lt. 1d-12) then + nbfir = nbfir + 1 + bf_so_ir(nbfir,iuq) = ir + bf_per_ir(ir) = bf_per_ir(ir) + 1 + else if (val .lt. -1d-12) then + call errquit('sym_abelian_bas_info: neg. proj?',0) + end if + end do + end do +c +c Check no. of bf per irrep against standard routine +c + call sym_bas_irreps(basis, .false., bf_per_ir_cum) + do ir = 1, nop + if (bf_per_ir_cum(ir) .ne. bf_per_ir(ir)) call errquit + $ ('sym_abelian_bas_info: miscount of functions', ir) + end do +c + nbfir = 0 + do ir = 1, nop + bf_per_ir_cum(ir) = nbfir + nbfir = nbfir + bf_per_ir(ir) + end do + write(6,123) ' bf_per_ir_cum ', (bf_per_ir_cum(ir),ir=1,nop) + write(6,123) ' bf_per_ir ', (bf_per_ir(ir),ir=1,nop) + 123 format(1x,a,8i5) +c + write(6,*) ' Symmetry of SOs in natural order ' + nbfir = 1 + do iuq = 1, bf_n_uq + i = bf_uq(iuq) + write(6,124) nbfir, (zir(bf_so_ir(ir,iuq)),ir=1,bf_n_uq_op(i)) + 124 format(i5,2x,8(a4,1x)) + nbfir = nbfir + bf_n_uq_op(i) + end do +c +c The SO coefficients of the functions arsing from basis funciton i +c are given by the projector +c +c . |ir> = (1/h(i)) sum(unique ops for i) chartab(op,ir) |op * i> +c +c with +c +c . |op * i> = bf_phase(op,i) * |bf_op_map(op,i)> +c . h(i) = number of unique ops for i. +c +c Finally, we must compute the various ways of numbering the SOs +c +c The natural order (lab1) is as above and is used when generating +c integrals over the SOs, but a more useful order (lab2) orders +c the functions within each irrep. +c +c lab2(1:8) = 0 +c loop over unique shells +c . loop over functions in shell +c . loop over irreps +c . if (non-zero from projector) then +c . lab1 = lab1 + 1 +c . lab2(ir) = lab2(ir) + 1 +c +c i.e, lab1 = lab1(irrep, shell component, unique shell). +c +c lab2 will be the index within the irrep. of the SO ... to get the +c index within the whole basis you need to add bf_per_ir_cum(ir). +c +c We need maps between the two labeling schemes. +c +c Given the AO basis function number of the generator and the irrep +c . what is the natural SO number (lab1) +c . what is the sorted SO number (lab2) +c +c Given lab1 +c . what is the ir +c . what is the lab2 +c . what is the generator AO basis function +c +c Given ir and lab2 +c . what is the lab1 +c + call ifill(8,0,lab2,1) + call ifill(2*8*nbf,-99999999,so_uq_to_lab,1) + call ifill(3*nbf,-99999999,so_lab1,1) + call ifill(8*nbf, -99999999,so_lab2,1) + lab1 = 0 + do iuq = 1, bf_n_uq + i = bf_uq(iuq) + do opuq = 1, bf_n_uq_op(i) + ir = bf_so_ir(opuq,iuq) +c + lab1 = lab1 + 1 + lab2(ir) = lab2(ir) + 1 +c + so_uq_to_lab(1,ir,iuq) = lab1 + so_uq_to_lab(2,ir,iuq) = lab2(ir) + so_lab1(1,lab1) = ir + so_lab1(2,lab1) = lab2(ir) + so_lab1(3,lab1) = iuq + so_lab2(ir,lab2(ir)) = lab1 + end do + end do +c + write(6,*) ' Map from unique bf to lab1 and lab2' + do iuq = 1, bf_n_uq + i = bf_uq(iuq) + write(6,77) iuq, + $ (so_uq_to_lab(1,bf_so_ir(ir,iuq),iuq), + $ so_uq_to_lab(2,bf_so_ir(ir,iuq),iuq), + $ ir=1,bf_n_uq_op(i)) + 77 format(i5,8(2x,2i4)) + end do + write(6,*) ' Map from lab1 to ir, lab2 and unique' + do lab1 = 1, nbf + write(6,78) lab1, (so_lab1(i,lab1), i=1,3) + 78 format(i5,2x,3i4) + end do + write(6,*) ' Map from ir+lab2 to lab1 ' + do ir = 1, 8 + do i = 1, bf_per_ir(ir) + write(6,79) ir, i, so_lab2(ir,i) + 79 format(2i5,2x,i5) + end do + end do +c + end + subroutine sym_abelian_char_tab(geom, char_tab, zir) + implicit none +c + integer geom ! [input] + double precision char_tab(8,8) ! [output] char_tab(op,irrep) + character*8 zir(8) +c +c Return the character table conveniently dimensioned +c and with the classes/operations reordered so that +c they are in the correct order. +c + integer maxireps, maxops + parameter (maxireps = 20, maxops=120) + integer op_class_index(maxops) + integer nop, iclass_dim(maxireps) + character*8 zclass(maxireps), group_name + double precision chars(maxireps*maxireps) +c + integer ir, op, class, nir + logical sym_char_table + external sym_char_table +c + call sym_op_classify(geom, op_class_index) + call sym_group_name(geom,group_name) + if (.not. sym_char_table(group_name, nop, nir, iclass_dim, + $ zir, zclass, chars)) + $ call errquit('sym_abelian_char_table: no char table',0) +c + call dfill(8*8, 0.0d0, char_tab, 1) +c + do ir = 1, nir + char_tab(1,ir) = chars((ir-1)*nir+1) + do op = 2, nop + class = op_class_index(op-1) ! Note -1 + char_tab(op,ir) = chars((ir-1)*nir+class) + end do + end do +c + write(6,*) ' Character table in natural order ' + call output(char_tab, 1, nop, 1, nir, 8, 8, 1) +c + end +