1      logical function sym_atom(geom, iat, q1)
2      implicit none
3      integer geom              ! Geometry handle [input]
4      integer iat               ! Atom index [input]
5      double precision q1       ! Constituency number [output]
6c
7c     Return true if (iatom) is the lexically highest atom
8c     symmetry equivalent atoms. If true also return the
9c     constituency factor q1 (= no. of symmetry equivalent atoms)
10c
11c     As an optimization inline sym_center_map ... the only reason
12c     geomP.fh and mafdecls.fh are included here ... the interface
13c     is compatible with that of the actual routine.
14c
15#include "nwc_const.fh"
16#include "geomP.fh"
17#include "mafdecls.fh"
18      integer iat_new
19      integer nops, op
20      integer n                 ! Counts no. of equivalent pairs
21      integer sym_number_ops, sym_center_map_inline
22      external sym_number_ops
23#include "itri.fh"
24      sym_center_map_inline(geom,iat,op) = int_mb(op - 1 +
25     $     sym_center_map_index(geom) + (iat-1)*sym_num_ops(geom))
26c
27      q1 = 0.0d0
28      sym_atom = .false.
29c
30c     Loop thru operations in the group and map to new pairs
31c
32      nops = sym_number_ops(geom)
33      n = 1                     ! Identity always counts
34      do op = 1, nops
35c
36c     Map centers
37c
38         iat_new = sym_center_map_inline(geom, iat, op)
39c
40c     Compare index
41c
42         if (iat .lt. iat_new) then
43            return
44         else if (iat .eq. iat_new) then
45            n = n + 1
46         endif
47      end do
48c
49      q1 = dble(nops + 1) / dble(n)
50      sym_atom = .true.
51c
52      end
53      logical function sym_shell(basis, ishell, q1)
54C$Id$
55      implicit none
56#include "errquit.fh"
57#include "bas.fh"
58#include "geom.fh"
59#include "global.fh"
60      integer basis             ! Basis set handle [input]
61      integer ishell            ! Shell index [input]
62      double precision q1       ! Constituency number [output]
63c
64c     Return true if (ishell) is the lexically highest shell
65c     of symmetry equivalent shells. If true, also return the
66c     constituency factor q1 (= no. of symmetry equivalent shells).
67c
68      integer ice
69      integer geom
70      logical sym_atom
71      external sym_atom
72c
73      sym_shell = .false.
74c
75c     Get geometry handle, centers where shells are located and
76c     number of group operations
77c
78      if (.not. bas_geom(basis, geom)) call errquit
79     $     ('sym_shell: bas_geom?', 0, BASIS_ERR)
80      if (.not. bas_cn2ce(basis, ishell, ice)) call errquit
81     $     ('sym_shell_pair: bas_cn2ce', ishell, BASIS_ERR)
82c
83      sym_shell = sym_atom(geom, ice, q1)
84c
85      end
86      logical function sym_shell_pair(basis, ishell, jshell, q2)
87C$Id$
88      implicit none
89#include "errquit.fh"
90#include "bas.fh"
91#include "geom.fh"
92#include "global.fh"
93      integer basis             ! Basis set handle [input]
94      integer ishell, jshell    ! Shell indices [input]
95      double precision q2       ! Constituency number [output]
96c
97c     Return true if (ishell,jshell) is the lexically highest
98c     pair of symmetry equivalent shells. If true, also return the
99c     constituency factor q2 (= no. of symmetry equivalent pairs).
100c
101c     This routine uses the exchange symmetry ishell <-> jshell
102c     and incorporates a factor of two into q2 to account for
103c     this.
104c
105      integer ice, jce
106      integer geom
107      logical sym_atom_pair
108      external sym_atom_pair
109c
110      sym_shell_pair = .false.
111      if (ishell.lt.jshell) return
112c
113c     Get geometry handle, centers where shells are located and
114c     number of group operations
115c
116
117*      write(6,*) ' sym_shell_pair ', ga_nodeid(), ishell, jshell
118*      call util_flush(6)
119*      call ga_sync()
120
121      if (.not. bas_geom(basis, geom)) call errquit
122     $     ('sym_shell_pair: bas_geom?', 0, BASIS_ERR)
123      if (.not. bas_cn2ce(basis, ishell, ice)) call errquit
124     $     ('sym_shell_pair: bas_cn2ce', ishell, BASIS_ERR)
125      if (.not. bas_cn2ce(basis, jshell, jce)) call errquit
126     $     ('sym_shell_pair: bas_cn2ce', jshell, BASIS_ERR)
127c
128*      write(6,*) ' sym_shell_pair centers', ga_nodeid(), ice, jce
129*      call util_flush(6)
130*      call ga_sync()
131c
132      sym_shell_pair = sym_atom_pair(geom, ice, jce, q2)
133      if (ishell .ne. jshell) q2 = q2 + q2
134c
135*      write(6,*) ' sym_shell_pair return ', ga_nodeid(),
136*     $     sym_shell_pair, q2
137*      call util_flush(6)
138*      call ga_sync()
139c
140      end
141      logical function sym_atom_pair(geom, iat, jat, q2)
142      implicit none
143      integer geom              ! Geometry handle [input]
144      integer iat, jat          ! Atom indices [input]
145      double precision q2       ! Constituency number [output]
146c
147c     Return true if (iatom,jatom) is the lexically highest
148c     pair of symmetry equivalent atom. If true also return the
149c     constituency factor q2 (= no. of symmetry equivalent pairs)
150c
151c     This routine uses the exchange symmetry iat <-> jat
152c     but does not incorporate any factors into q2 to account for
153c     this (q2 is point group symmetry only).
154c
155c     As an optimization inline sym_center_map ... the only reason
156c     geomP.fh and mafdecls.fh are included here ... the interface
157c     is compatible with that of the actual routine.
158c
159#include "nwc_const.fh"
160#include "geomP.fh"
161#include "mafdecls.fh"
162      integer iat_new, jat_new, ijat, ijat_new
163      integer nops, op
164      integer n                 ! Counts no. of equivalent pairs
165      integer sym_number_ops, sym_center_map_inline
166      external sym_number_ops
167#include "itri.fh"
168      sym_center_map_inline(geom,iat,op) = int_mb(op - 1 +
169     $     sym_center_map_index(geom) + (iat-1)*sym_num_ops(geom))
170c
171      q2 = 0.0d0
172      sym_atom_pair = .false.
173      if (iat .lt. jat) return
174      ijat = itri(iat,jat)
175c
176c     Loop thru operations in the group and map to new pairs
177c
178      nops = sym_number_ops(geom)
179      n = 1                     ! Identity always counts
180      do op = 1, nops
181c
182c     Map centers
183c
184         iat_new = sym_center_map_inline(geom, iat, op)
185         jat_new = sym_center_map_inline(geom, jat, op)
186c
187c     Compare canonical indices
188c
189         ijat_new = itri(iat_new, jat_new)
190         if (ijat .lt. ijat_new) then
191            return
192         else if (ijat .eq. ijat_new) then
193            n = n + 1
194         endif
195      end do
196c
197      q2 = dble(nops + 1) / dble(n)
198      sym_atom_pair = .true.
199c
200      end
201      logical function sym_atom_quartet(geom, iat, jat, kat, lat, q4)
202      implicit none
203      integer geom              ! Geometry handle [input]
204      integer iat, jat, kat, lat ! Atom indices [input]
205      double precision q4       ! Constituency number [output]
206c
207c     Return true if (iatom,jatom,katom,latom) is the lexically highest
208c     quartet of symmetry equivalent atoms. If true also return the
209c     constituency factor q4 (= no. of symmetry equivalent quartets)
210c
211c     This routine uses the standard three exchange symmetries
212c
213c     (iat<->jat) <-> (kat<->lat)
214c
215c     ... it is easy to extend it to not use these
216c
217#include "nwc_const.fh"
218#include "geomP.fh"
219#include "mafdecls.fh"
220      integer iat_new, jat_new, kat_new, lat_new
221      integer nops, op, ij, kl, ijkl, ijkl_new, ij_new, kl_new
222      integer n                 ! Counts no. of equivalent pairs
223      integer sym_number_ops, sym_center_map_inline
224      external sym_number_ops
225#include "itri.fh"
226      sym_center_map_inline(geom,iat,op) = int_mb(op - 1 +
227     $     sym_center_map_index(geom) + (iat-1)*sym_num_ops(geom))
228c
229      sym_atom_quartet = .false.
230c
231c     Assume that the code is looping thru a unique set of
232c     indices but maybe not with conventional canonical order
233c     ... so don't do the canonical check here
234*      if (iat .lt. jat) return  ! Labels must be in canonical order
235*      if (iat .lt. kat) return
236*      if (kat .lt. lat) return
237*      if (iat .eq. kat .and. jat .lt. lat) return
238c
239      q4 = 0.0d0
240      ij   = itri(iat,jat)
241      kl   = itri(kat,lat)
242      ijkl = itri(ij, kl)
243c
244c     Loop thru operations in the group and map to new pairs
245c
246      nops = sym_number_ops(geom)
247      n = 1                     ! Identity always counts
248      do op = 1, nops
249c
250c     Map centers
251c
252         iat_new = sym_center_map_inline(geom, iat, op)
253         jat_new = sym_center_map_inline(geom, jat, op)
254         kat_new = sym_center_map_inline(geom, kat, op)
255         lat_new = sym_center_map_inline(geom, lat, op)
256c
257c     Compare canonical indices
258c
259         ij_new   = itri(iat_new,jat_new)
260         kl_new   = itri(kat_new,lat_new)
261         ijkl_new = itri(ij_new, kl_new)
262c
263         if (ijkl .lt. ijkl_new) then
264            return
265         else if (ijkl .eq. ijkl_new) then
266            n = n + 1
267         endif
268      end do
269c
270      q4 = dble(nops + 1) / dble(n)
271      sym_atom_quartet = .true.
272c
273      end
274      logical function sym_shell_quartet(basis,
275     $     ishell, jshell, kshell, lshell, q4)
276      implicit none
277#include "errquit.fh"
278#include "bas.fh"
279#include "geom.fh"
280#include "global.fh"
281      integer basis             ! Basis set handle [input]
282      integer ishell, jshell    ! Shell indices [input]
283      integer kshell, lshell    ! Shell indices [input]
284      double precision q4       ! Constituency number [output]
285c
286      integer ice, jce, kce, lce
287      integer geom
288      logical sym_atom_gen_quartet
289      external sym_atom_gen_quartet
290c
291      q4 = 0.0d0
292      sym_shell_quartet = .false.
293      if (ishell.lt.jshell) return
294      if (kshell.lt.lshell) return
295c
296c     Get geometry handle, centers where shells are located and
297c     number of group operations
298c
299      if (.not. bas_geom(basis, geom)) call errquit
300     $     ('sym_shell_pair: bas_geom?', 0, BASIS_ERR)
301      if (.not. bas_cn2ce(basis, ishell, ice)) call errquit
302     $     ('sym_shell_pair: bas_cn2ce', ishell, BASIS_ERR)
303      if (.not. bas_cn2ce(basis, jshell, jce)) call errquit
304     $     ('sym_shell_pair: bas_cn2ce', jshell, BASIS_ERR)
305      if (.not. bas_cn2ce(basis, kshell, kce)) call errquit
306     $     ('sym_shell_pair: bas_cn2ce', kshell, BASIS_ERR)
307      if (.not. bas_cn2ce(basis, lshell, lce)) call errquit
308     $     ('sym_shell_pair: bas_cn2ce', lshell, BASIS_ERR)
309c
310      sym_shell_quartet = sym_atom_gen_quartet(
311     $     geom, ice, jce, kce, lce, q4)
312c
313      end
314      logical function sym_atom_gen_quartet(geom,
315     $     iat, jat, kat, lat, q4)
316      implicit none
317#include "errquit.fh"
318      integer geom              ! Geometry handle [input]
319      integer iat, jat, kat, lat ! Atom indices [input]
320      double precision q4       ! Constituency number [output]
321c
322#include "nwc_const.fh"
323#include "geomP.fh"
324#include "mafdecls.fh"
325      integer iat_new, jat_new, kat_new, lat_new
326      integer nops, op, ij, kl, ij_new, kl_new
327      integer n                 ! Counts no. of equivalent pairs
328      integer sym_number_ops, sym_center_map_inline
329      external sym_number_ops
330#include "itri.fh"
331      sym_center_map_inline(geom,iat,op) = int_mb(op - 1 +
332     $     sym_center_map_index(geom) + (iat-1)*sym_num_ops(geom))
333c
334      q4 = 0.0d0
335      sym_atom_gen_quartet = .false.
336c
337      if (iat .lt. jat) return  ! Labels must be in canonical order
338      if (kat .lt. lat) return
339c
340      ij   = itri(iat,jat)
341      kl   = itri(kat,lat)
342c
343c     Loop thru operations in the group and map to new pairs
344c
345      nops = sym_number_ops(geom)
346      n = 1                     ! Identity always counts
347      do op = 1, nops
348c
349c     Map centers
350c
351         iat_new = sym_center_map_inline(geom, iat, op)
352         jat_new = sym_center_map_inline(geom, jat, op)
353         kat_new = sym_center_map_inline(geom, kat, op)
354         lat_new = sym_center_map_inline(geom, lat, op)
355c
356c     Compare canonical indices
357c
358         ij_new   = itri(iat_new,jat_new)
359         kl_new   = itri(kat_new,lat_new)
360c
361         if (ij .lt. ij_new) return
362         if (ij.eq.ij_new .and. kl.lt.kl_new) return
363         if (ij.eq.ij_new .and. kl.eq.kl_new) then
364            n = n + 1
365         endif
366      end do
367c
368      q4 = dble(nops+1) / dble(n)
369      if (abs(q4-nint(q4)).gt.1e-12) call errquit
370     $     ('sym_atom_gen_quartet: not divisible', 0, BASIS_ERR)
371      sym_atom_gen_quartet = .true.
372c
373      end
374