comp-chem/NWChem
1
0
Fork 0
mirror of https://github.com/nwchemgit/nwchem.git synced 2026-07-28 22:25:48 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 1

symmetry orbital infor for abelian groups

Browse source
This commit is contained in:
Robert Harrison 2000-05-04 21:18:40 +00:00
parent 8bd0f18fdf
commit 7784c0c0df
1 changed files with 429 additions and 1 deletions
Download patch file Download diff file Expand all files Collapse all files

430
src/symmetry/sym_abelian.F
View file

@ -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