1c***********************************************************************
2c  Subroutine sym_nwc: a 'friend' of geom
3c
4c Driver routine:Attributes
5c
6c   calls gensym chain:  (creates matrix representaions of the
7c                         symmetry operators)
8c   calls dosymops:      (applies operatorst to symmetry distinct atoms
9c                         prints and stores expanded lists)
10c   This routine has direct access to geom private include files and
11c   thus common blocks.
12c
13c   calls sym_map: creates atom pair mapping relations based upon the
14c                  symmetry operations and deposits symmetry operations
15c                  and mapping tables to rtdb for future retrieval.
16c
17c
18c   call sym_init_inv_op: to make list of inverse operations
19c
20c                                          A.C.H. 5/12/94
21c                 Interfaces to NWCHEM     A.C.H. 10/10/94 + RJH
22c
23c***********************************************************************
24      subroutine sym_nwc(geom,rtdb,nata,oprint,scale,threquiv,nops)
25C$Id$
26      implicit real*8 (a-h,o-z)
27#include "errquit.fh"
28      parameter(maxops=192)
29      Parameter (EPS=1.D-14)
30
31#include "mafdecls.fh"
32#include "geom.fh"
33#include "nwc_const.fh"
34#include "geomP.fh"
35#include "rtdb.fh"
36
37      logical LResult, oprint
38      integer geom,rtdb
39      dimension cdist(3),cang(3)
40      double precision mass_new(max_cent), charge_new(max_cent)
41      double precision invnucexp_new(max_cent)
42C      external pertbl
43c
44c--> parameters needed by sym chain
45c
46      itype  = isystype(geom)
47      numgrp = group_number(geom)
48      numset = setting_number(geom)
49c
50c   lattice vectors & angles (cell constants)
51c   values should have been read in and initialized properly for
52c   each system dimension
53c
54      do 1000 i=1,3
55        cdist(i)=lattice_vectors(i,geom)
56        cang(i) =lattice_angles(i,geom)
57 1000 continue
58c
59c    generate symmetry operators for whatever group (crystals, surfaces
60c    polymers or molecules) have been requested.
61c
62      if ((itype.eq.3).and.(numset.gt.2)) then
63          !write(*,*) "generating unsual groups"
64          call gensym_extra(itype,group_number(geom),numset,
65     $                      sym_ops(1,1,geom),
66     $                      nops,oprint, group_name(geom))
67      else
68         call gensym(itype,numgrp,numset,sym_ops(1,1,geom),
69     $               nops,oprint, group_name(geom),geom,rtdb)
70      end if
71
72c
73c     scratch space needed by dosymops: tags, coordinates
74c     space on stack: pointer list, points to the assymetric (distinct) atoms
75c     apply symmops to coordinate list
76c
77      LResult = MA_Push_Get(MT_Byte,nata*(nops+1)*16,'scratch labels',
78     &     latags, iatags)
79      LResult = MA_Push_Get(MT_Dbl,nata*(nops+1)*3,'scratch atoms',
80     &     lcoord, icoord) .and. lresult
81      Lresult = MA_Push_Get(MT_Int,nata,'distinct atoms',ldatom,idatom)
82     $     .and. lresult
83      if (.not. lresult) call errquit('sym_nwc: dosymops ma failed', 0,
84     &       MA_ERR)
85c
86      call dosymops(sym_ops(1,1,geom),nops,coords(1,1,geom),
87     $     tags(1,geom),nata,
88     &     itype,Byte_MB(iatags),Dbl_MB(icoord),nareal,nata*(nops+1),
89     $     cdist,cang,Int_MB(idatom),
90     $     charge(1,geom), charge_new, geom_mass(1,geom),
91     $     mass_new, geom_invnucexp(1,geom),invnucexp_new,threquiv)
92c
93      if(nareal.gt.max_cent) then
94        write(6,*)' Too many atoms, increase nw_max_atom in ',
95     &       'util/nwc_const.fh error in sym_nwc : ', nareal
96        call errquit('sym: too many atoms', 1, INPUT_ERR)
97      endif
98c
99c     Allocate memory for the center map on the heap since it persists
100c
101      if (nops .gt. 0) then
102      if (.not. MA_Alloc_Get(MT_Int,nops*nareal,'atom sym map',
103     &     lmscr,imscr)) call errquit
104     $     ('sym_nwc: ma_alloc_get of sym map failed', nops*nareal,
105     &       MA_ERR)
106      else
107       lmscr = -1
108       imscr = 99999999
109      endif
110c
111c--> derive mapping table for atom pairs permuted by sym operations
112c--> NOTE: at this point coords will contain cartesian coordinates
113c          and dbl_mb(icoord) has the fractional coordinates
114c          (unless system is molecular in which case is also cartesian)
115c
116      if (nops .gt. 0) then
117        if (.not. MA_Push_Get(MT_Dbl,nops*3,'xnew scratch',
118     &       lxnew,ixnew)) call errquit
119     $       ('sym_nwc: ma for xnew scratch', nops*3, MA_ERR)
120      else
121        ixnew = 0
122      endif
123      if (.not. MA_Push_Get(MT_INT,nareal,'ilbl',
124     &     llbl,ilbl)) call errquit
125     $     ('sym_nwc: ma for ilbl', nareal, MA_ERR)
126c
127      call sym_map(itype, Byte_MB(iatags), sym_ops(1,1,geom),
128     $     nops, nareal, Int_Mb(imscr),Dbl_Mb(ixnew),Dbl_Mb(icoord),
129     $     oprint,INT_MB(ilbl), threquiv)
130c
131      LResult = MA_Pop_Stack(llbl)
132      if (nops .gt. 0) then
133        LResult = MA_Pop_Stack(lxnew)
134      endif
135c
136c     Insert new information into the geometry data structures
137c     Geom will hold the cartesian coords ... fractional ones
138c     remain in icoord.
139c     Note that the memory for the map (lmscr) is now freed in geom_destroy
140c
141      call sym_put_in_geom(geom, nata, nareal,
142     $     Byte_MB(iatags), int_mb(idatom), dbl_mb(icoord),
143     $     lmscr, imscr, nops, charge_new, mass_new, invnucexp_new)
144c
145c ... sloppy way for the moment to fix units for 2-d case
146c
147      if(itype.eq.2) then
148         do i=1,nareal
149            coords(3,i,geom)=coords(3,i,geom)*scale
150         enddo
151      endif
152c
153c---> print gmat,amat,vol evaluated by geom_nd routines
154c
155c
156c---> print gmat,amat,vol evaluated by geom_nd routines
157c
158c
159c--> Print Amatrix, G matrix, volume for 3-d and 2-d cells.
160c
161c      if(itype.eq.3) then
162c         write(*,9)
163c         do 513 i=1,3
164c            write(*,15) (amatrix(i,j,geom), j=1,3)
165c 513     continue
166c
167c         write(*,21)
168c         do 517 i=1,3
169c            write(*,15) (amatrix_inv(i,j,geom), j=1,3)
170c 517     continue
171c
172c         write(*,10)
173c         do 130 i=1,3
174c            write(*,15) (metric_matrix(i,j,geom), j=1,3)
175c 130     continue
176c
177c         write(*,17) volume_direct(geom)
178c
179c      elseif(itype.eq.2) then
180c         write(*,9)
181c         do 630 i=1,3
182c            write(*,15) (amatrix(i,j,geom), j=1,3)
183c 630     continue
184c         write(*,21)
185c         do 631 i=1,3
186c            write(*,15) (amatrix_inv(i,j,geom), j=1,3)
187c 631     continue
188c
189c         write(*,10)
190c         do 640 i=1,3
191c            write(*,15) (metric_matrix(i,j,geom), j=1,3)
192c 640     continue
193c
194c         write(*,18) volume_direct(geom)
195c      endif
196c
197c--> Print cartesian coordinates of unit cell.
198c
199c      if(itype.eq.3) then
200c        write(*,30)
201c        l=0
202c        do 510 i=1,nareal
203c           write(*,20) i,tags(i,geom),(coords(j,i,geom), j=1,3)
204c           l=l+3
205c 510    continue
206c         write(*,14)
207c   convert z-coords of slab to bohr from angstroms
208c      elseif(itype.eq.2) then
209c         do i=1,nareal
210c            coords(3,i,geom)=coords(3,i,geom)*angstrom_to_au
211c         enddo
212c         write(*,30)
213c         l=0
214c         do 512 i=1,nareal
215c            write(*,20) i,tags(i,geom),(coords(j,i,geom), j=1,3)
216c            l=l+3
217c 512     continue
218c         write(*,14)
219c      endif
220c
221c--> final clean up of memory: gensym and dosymops activities
222c
223      LResult = MA_Pop_Stack(ldatom)
224      LResult = MA_Pop_Stack(lcoord)
225      LResult = MA_Pop_Stack(latags)
226c
227      call sym_init_inv_op(geom)
228c
229c--> formats
230c
231 9    format(//20x,
232     &     'CRYSTALLOGRAPHIC TO CARTESIAN MATRIX (A.U.)',/)
233 10   format(//,23x,'METRIC MATRIX:DIRECT LATTICE (A.U.)',/)
234 14   format(/,80('-'),//)
235 15   format(20x,3(f12.6,2x))
236 17   format(/,19x,'UNIT CELL VOLUME = ',f12.6,3x,'BOHR**3')
237 18   format(/,19x,'UNIT CELL AREA = ',f12.6,3x,'BOHR**3')
238 20   format(10x,i5,2x,a16,3(f10.6,2x))
239 21   format(//20x,
240     & 'CARTESIAN TO CRYSTALLOGRAPHIC MATRIX (A.U.)',/)
241 30   format(/,15x,'--------- CARTESIAN COORDINATES (A.U.) ---------',
242     &/)
243      return
244      end
245
246
247
248