1*
2* $Id$
3*
4
5      SUBROUTINE vib_RDINP(
6     &    HESS,HESST,COORD,ATMASS,ZAN,OLDLBLAT,IZMAT,NELS,
7     &    PROJEC,ZEROPE,HESOUT,INTERN,
8     &    rtdb,first_pass)
9C
10C This routine reads the user input and the hessian matrix from
11C      the appropriate place.  The default is to use the "default"
12C      masses and to read the analytic hessian from tape 10.
13C
14      IMPLICIT NONE ! REAL*8 (A-H,O-Z)
15#include "errquit.fh"
16#include "mafdecls.fh"
17#include "stdio.fh"
18#include "rtdb.fh"
19#include "geom.fh"
20#include "bas.fh"
21      INTEGER NELS
22      LOGICAL SETMAS,PROJEC,ZEROPE,HESOUT,INTERN
23      logical dum_log
24      integer natom, nat3, nhess, nhesst
25      logical first_pass
26      logical o_new_masses
27      DOUBLE PRECISION HESS(NAT3,NAT3)  ! Hessian matrix
28      DOUBLE PRECISION HESST(NHESST)    ! Lower triangle of hessian matrix
29      DOUBLE PRECISION COORD(3,NATOM)   ! Atom x,y,z coordinates
30      DOUBLE PRECISION ATMASS(NATOM)    ! Atoms mass
31      DOUBLE PRECISION ZAN(NATOM)       ! Nuclear charge (i.e. atomic number)
32      INTEGER OLDLBLAT(2,NATOM)     ! Atom character labels
33      INTEGER IZMAT(NELS)
34      integer num_new_masses
35      integer geom_index
36      double precision new_mass
37      character*32 mass_rtdb_id
38      character*16 mass_tag
39      character*16 tag
40      character*22 element
41      character*2 symbol
42      double precision q
43      double precision xyz(3)
44      COMMON /cvib_HESS/ NATOM,NAT3,NHESS,NHESST   ! Hessian information
45      integer numans
46      double precision ams, wave
47      COMMON /cvib_SETCON/ AMS(36),WAVE,NUMANS     ! setup parameters
48c
49      integer h_lblb, i_lblb, geom
50      integer il_cnt, iii, ijunk, jjj
51c
52      integer rtdb
53      HESOUT = .FALSE.            ! do not write Hessian to tape 10
54      o_new_masses = .false.
55*
56      if (first_pass) WRITE(luout,*)' Vib: Default input used '
57      PROJEC = .FALSE.
58      ZEROPE = .FALSE.
59      HESOUT = .FALSE.
60      INTERN = .FALSE.
61      if (rtdb_get(rtdb,'vib:project',MT_LOG,1,dum_log)) then
62        projec = dum_log
63        write(luout,*)' vib:project option set to',projec
64      endif
65      if (rtdb_get(rtdb,'vib:zero point energy',MT_LOG,1,dum_log)) then
66        zerope = dum_log
67        write(luout,*)' vib:zero point energy option set to',zerope
68      endif
69      if (zerope) projec = .true.
70*
71* over ride all input to do normal and projected analysis
72*
73      if (first_pass) then
74        projec = .false.
75        zerope = .false.
76      else
77        projec = .true.
78        zerope = .true.
79      endif
80*.. hessian read in vib_vib hess and hesst set there check it here
81      call vib_chkhess(hess,nat3,first_pass)
82* set labels, zan and coords information and mass information
83c
84      if (.not.ma_push_get(MT_BYTE,(2*natom),' labels for vib',
85     &    h_lblb,i_lblb))
86     &    call errquit('vib_rdinp ma_get for labels failed',911, MA_ERR)
87      if (.not.geom_create(geom,'geometry')) call errquit
88     &    ('vib_rdinp: error creating geometry',911, GEOM_ERR)
89      if (.not.geom_rtdb_load(rtdb,geom,'geometry')) call errquit
90     &    ('vib_rdinp: error loading geometry',911, RTDB_ERR)
91c
92c see if ANY mass is overwritten by freq/vib input block
93c
94      if (rtdb_get
95     &    (rtdb,'vib:remass:count',mt_int,1,num_new_masses)) then
96        o_new_masses = .true.
97      else
98        num_new_masses = 0
99        o_new_masses = .false.
100      endif
101      il_cnt = 0
102      do iii = 1,natom
103        if(.not.geom_cent_get(geom,iii,tag,xyz,q))
104     &      call errquit('vib_rdinp: geom_cent_get failed',911,
105     &       GEOM_ERR)
106        zan(iii) = q
107        coord(1,iii) = xyz(1)
108        coord(2,iii) = xyz(2)
109        coord(3,iii) = xyz(3)
110        if(.not.geom_tag_to_element(tag,symbol,element,ijunk)) then
111          if (symbol.ne.'bq')
112     &        call errquit('vib_rdinp: tag2elem fail',911, GEOM_ERR)
113        endif
114        byte_mb(il_cnt+i_lblb) = symbol(1:1)
115        il_cnt = il_cnt+1
116        byte_mb(il_cnt+i_lblb) = symbol(2:2)
117        il_cnt = il_cnt+1
118      enddo
119
120
121*
122* this loop structure separated out to do one pass over data in rtdb.
123* set the incore geometry with any mass information from the freq/vib
124* input.
125*
126      if (o_new_masses) then
127        do jjj = 1,num_new_masses
128          write(mass_rtdb_id,11111)jjj          ! lexical index
129          if (rtdb_get(rtdb,mass_rtdb_id,mt_int,1,geom_index)) then
130            write(mass_rtdb_id,11113)jjj
131            if (.not.rtdb_get(rtdb,mass_rtdb_id,mt_dbl,1,new_mass))
132     &          call errquit
133     &          ('vib_rdinp: rtdb get failed for mass',911, RTDB_ERR)
134            if (.not.geom_mass_set(geom,geom_index,new_mass))
135     &          call errquit
136     &          ('vib_rdinp: geom_mass_set failed for mass',911,
137     &       GEOM_ERR)
138          else  ! check for tag input
139            write(mass_rtdb_id,11112)jjj
140            if (rtdb_cget(rtdb,mass_rtdb_id,1,mass_tag)) then
141              write(mass_rtdb_id,11113)jjj
142              if (.not.rtdb_get(rtdb,mass_rtdb_id,mt_dbl,1,new_mass))
143     &            call errquit
144     &            ('vib_rdinp: rtdb get failed for mass',911, RTDB_ERR)
145* tag input found loop over all geometry atoms for a possible match
146              do iii = 1,natom
147                if(.not.geom_cent_get(geom,iii,tag,xyz,q))
148     &              call errquit('vib_rdinp: geom_cent_get failed',911,
149     &       GEOM_ERR)
150                if (bas_do_tags_match(mass_tag,tag)) then
151                  if (.not.geom_mass_set(geom,iii,new_mass))
152     &                call errquit
153     &                ('vib_rdinp: geom_mass_set failed for mass',911,
154     &       GEOM_ERR)
155                endif
156              enddo
157            else
158              if(.not.rtdb_print(rtdb,.true.)) call errquit('911',911,
159     &       RTDB_ERR)
160              write(luout,*)' expected mass input for index ',jjj
161              write(luout,*)' It was neither index or tag based '
162              call errquit
163     &            ('vib_rdinp: fatal error for mass informaition',911,
164     &       RTDB_ERR)
165            endif
166          endif
167        enddo
168      endif
169*
170* best place to set mass is in geometry object
171* modified in-core geometry object with new masses from
172* vibrational/freq input now reset's masses.
173*
174      setmas = geom_masses_get(geom,natom,atmass)
175      if (.not.setmas)
176     &    call errquit
177     &    ('vib_rdinp: geom_masses_get failed',911, GEOM_ERR)
178*
179      if (.not.geom_destroy(geom)) call errquit
180     &    ('vib_rdinp: geom_destroy failed',911, GEOM_ERR)
181C
182C Write out atom information read here
183C
184      if (first_pass) then
185        WRITE(6,10000)
186C
187        il_cnt = 0
188        DO 00300 III = 1,NATOM
189          WRITE(6,10001)
190     &        byte_mb(il_cnt+i_lblb),
191     &        byte_mb(il_cnt+i_lblb+1),
192     &        III,
193     &        COORD(1,III),COORD(2,III),COORD(3,III),
194     &        ATMASS(III)
195          il_cnt = il_cnt + 2
19600300   CONTINUE
197        WRITE(6,10002)
198      endif
199      if (.not.ma_pop_stack(h_lblb)) call errquit
200     &    ('vib_rdinp ma_pop failed',911, MA_ERR)
201      RETURN ! leave routine
20210000 FORMAT(//,1X,28('-'),' Atom information ',28('-'),/,
203     +5X,'atom',4X,'#',8X,'X',14X,'Y',14X,'Z',12X,'mass',/,
204     +1X,74('-'))
20510001 FORMAT(4x,a1,a1,3x,I5,4(1PD15.7))
20610002 FORMAT(1X,74('-')///)
20711111 format('vib:remass:',i4,':lexi:')
20811112 format('vib:remass:',i4,':tags:')
20911113 format('vib:remass:',i4,':mass:')
210C
211      END
212
213
214