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