1C $Id$
2      Subroutine nwc_gen (x,v,dx,fmat,hess,xs,natom,n3,n3tm,ndbar,op)
3C
4C  Run NWChem task to get potential and first and second derivatives
5C   for geometry x.
6C
7      Implicit None
8#include "errquit.fh"
9#include "nwc_drdyP.fh"
10#include "geom.fh"
11#include "mafdecls.fh"
12#include "printlevels.fh"
13#include "rtdb.fh"
14#include "stdio.fh"
15      Integer natom, n3, n3tm, ndbar, op
16      Double Precision v
17      Double Precision x(n3),dx(n3),fmat(n3tm,n3),
18     &      hess(ndbar),xs(n3), velocities(3,n3)
19*
20      Integer ndima
21      Parameter (ndima=3)
22      Double Precision amat(ndima,ndima),
23     &      rr(ndima,ndima),rrs(ndima,ndima),
24     *   rri(ndima,ndima)
25*
26      Integer geom
27      Integer i,j,jj
28      double precision threquiv
29      character*16 groupname
30      character*255 dummy_file
31      Integer ncenter_B4_autosym, ncenter, nata
32      Integer current_print_level,nops
33      logical oautosym
34      logical status_rtdb, status_ignore
35*     logical copy_sym, copy_c1, oautosym
36*
37
38      call util_print_get_level(current_print_level)
39      call util_print_set_level(print_none)
40*
41      if (op.eq.DRDY_CODE_HESSIAN) then
42*
43* ... delete finite difference files that may exist.
44        call util_file_name('hess',  .false., .false.,dummy_file)
45        call util_file_unlink(dummy_file)
46        call util_file_name('fd_ddipole',.false.,.false.,dummy_file)
47        call util_file_unlink(dummy_file)
48      endif
49*
50      if (.not.geom_create(geom,'geometry'))
51     &      call errquit('nwc_gen: geom_create failed',911,
52     &       GEOM_ERR)
53      if (.not.geom_set_user_units(geom,'a.u.'))
54     &      call errquit('nwc_gen: geom_set_user_units failed',911,
55     &       GEOM_ERR)
56      if (.not.geom_cart_set(geom,natom,atomic_labels,x,atomic_charge))
57     &      call errquit('nwc_gen: geom_cart_set failed',911,
58     &       GEOM_ERR)
59      if (.not.geom_masses_set(geom,natom,nwcmass))
60     &      call errquit('nwc_gen:geom_masses_set failed',911,
61     &       GEOM_ERR)
62      if (.not.geom_atomct_set(geom,natom,atomic_atomct))
63     &      call errquit('nwc_gen:geom_atomct_set failed',911,
64     &       GEOM_ERR)
65*      if (.not.geom_rtdb_store(my_rtdb,geom,'geometry'))
66*     &      call errquit('nwc_gen: geom_rtdb_store failed',911,
67*     &       RTDB_ERR)
68c
69c Now figure out if we are going to use autosym and if so, get the
70c value for the threshold and take care of the symmetry
71c
72      status_rtdb   = rtdb_parallel(.false.)
73
74      if (.not. rtdb_get(my_rtdb,'drdy:autosym',mt_log,1,oautosym))
75     &  call errquit('nwc_gen: problem getting autosym from rtdb',555,
76     &       RTDB_ERR)
77      if (oautosym) then
78        if (.not.rtdb_get(my_rtdb,'drdy:threquiv',mt_dbl,1,threquiv))
79     &   call errquit('nwc_gen:problem getting threquiv from rtdb',555,
80     &       RTDB_ERR)
81        ncenter_B4_autosym = natom
82        ncenter = natom
83        call dcopy(n3,x,1,dx,1)   ! temporary use of dx
84        call dcopy(ncenter,atomic_charge,1,copy_charge,1)
85        call dcopy(ncenter,atomic_atomct,1,copy_atomct,1)
86        do i = 1,ncenter
87          copy_labels(i) = atomic_labels(i)
88        enddo
89        call dcopy(3*n3,0.d0,0,velocities,1)
90        call geom_auto_sym(my_rtdb,geom,dx,
91     &      copy_charge,copy_labels,copy_atomct,ncenter,
92     &      threquiv,groupname,velocities)
93        if (op.ne.DRDY_CODE_SPENERGY)
94     &    write(luout,10) groupname
95        if (groupname(1:2).ne."C1") then
96        if (geom_group_set(geom,groupname)) then
97          if (.not.geom_cart_set(geom,ncenter,copy_labels,dx,
98     &         copy_charge))
99     &      call errquit('nwc_gen: geom_cart_set failed', 0,
100     &       GEOM_ERR)
101          if (ncenter_B4_autosym .ne. ncenter) call errquit
102     &       ('nwc_gen: autosym bug : number of atoms wrong',
103     &         ncenter, INPUT_ERR)
104        else
105          write(luout,*) ' autosym detected unknown group ',
106     &                      groupname
107          call errquit('nwc_gen: autosym: invalid group',0, INPUT_ERR)
108        endif
109c
110c     Apply system and symmetry info to the list of
111c     unique centers build mapping tables set up coord lists
112c
113        nata = ncenter
114        call sym_nwc(geom,my_rtdb,nata,.false.,1.0d00,threquiv,nops,
115     %       .false.)
116c
117c     Check that if we used autosym that we ended up with the
118c     same no. of atoms ... if we don't then autosym and nwchemsym
119c     don't agree on the orientation of point group elements
120c
121        if (.not. geom_ncent(geom,ncenter)) call errquit
122     $       ('nwc_gen: geom_cent?',0, GEOM_ERR)
123        if (ncenter_B4_autosym .ne. ncenter) call errquit
124     $       ('nwc_gen: autosym bug : too many atoms',ncenter,
125     &       INPUT_ERR)
126c
127c     Force exact symetry on the coordinates
128c
129        call sym_geom_project(geom, threquiv)
130        endif
131      endif
132c
133c store as default geometry
134c
135      if (.not.geom_rtdb_store(my_rtdb,geom,'geometry'))
136     &      call errquit('nwc_gen: geom_rtdb_store failed',911,
137     &       RTDB_ERR)
138      if (.not.geom_destroy(geom))
139     &      call errquit('nwc_gen: geom_destroy failed',911, GEOM_ERR)
140      call util_print_set_level(current_print_level)
141*      write(luout,*)' x after  ',x
142c
143c Actually do the deed
144c
145      status_ignore   = rtdb_parallel(status_rtdb)
146      status_rtdb   = rtdb_parallel(.false.)
147
148      if ((op.ge.DRDY_CODE_SPENERGY).and.
149     &    (op.le.DRDY_CODE_HESSIAN)) then
150        call drdy_synch(op,'nwc_gen')
151      else
152        write(luout,*) 'Trying to calculate an unknown op code'
153        call errquit('nwc_gen: Unknown op code',op, INPUT_ERR)
154      endif
155      status_ignore   = rtdb_parallel(status_rtdb)
156      status_rtdb   = rtdb_parallel(.false.)
157c
158c Get the appropriate values to hand back
159c
160      call drdy_nwc_get_energy(my_rtdb,v)
161      if (op.ge.DRDY_CODE_GRADIENT) then
162         call drdy_nwc_get_coords(my_rtdb,xs)
163         call drdy_nwc_get_gradient(my_rtdb,dx)
164      endif
165      if (op.eq.DRDY_CODE_HESSIAN)
166     &   call drdy_nwc_get_hessian(my_rtdb,hess)
167*
168C
169C  Get transformation matrix (rotation matrix amat) from the gaussian
170C     standard orientation (xs) to the original orientation (x)
171      if (op.ge.DRDY_CODE_GRADIENT) then
172        call drdy_rotmat(x,xs,amat,rr,rrs,rri,natom,ndima)
173        call drdy_rotg(dx,amat,rr,natom,n3)
174      endif
175C
176      if (op.eq.DRDY_CODE_HESSIAN) then
177        jj = 0
178        do i = 1,n3
179          do j = 1,i
180            jj = jj + 1
181            fmat(i,j) = hess(jj)
182            fmat(j,i) = hess(jj)
183          enddo
184        enddo
185C  Transform hessian matrix from gaussian standard orientation to the
186C     original orientation
187        call drdy_rotf(fmat,amat,rr,natom,n3,n3tm)
188      endif
189*
190      status_ignore = rtdb_parallel(status_rtdb)
191      return
19210    Format (1x,'Symmetry group is ',a3)
1931000  Format (1x,79a1)
1941001  Format (' ')
1951002  Format (1x,2i2)
1961003  Format (1x,5a1,1p3e20.10)
197      End
198