1      subroutine getmdp(natoms,heat,izo,ido,istatz,cc,ianc,
2     &              bl,alph,bet,ibl,ialph,ibet,imap,ianz,iz,
3     &              c,cz,alpha,beta,ian)
4
5c this is really getmop
6
7      implicit double precision (a-h,o-z), integer ( i-n)
8      integer getlin
9      character*137 string
10      character*2 catom,catomt,tolowf,iel
11      logical first,ottest, oerror,zmstar,getnxt,gnreal
12
13      common /zmfrst/ ihaszm, nz, mxzat
14
15      dimension ian(*),c(3,*),cz(3,*),alpha(*),beta(*)
16      dimension ianc(*),cc(3,*)
17      dimension bl(*),alph(*),bet(*),ibl(*),ialph(*),ibet(*),
18     &          imap(*),ianz(*),iz(4,*)
19      dimension izt(4),iztt(3,3),zzz(3),idm(3),tmp(3)
20      dimension iel(100)
21      character*137 coplin
22      character*137 line
23      common /curlin/ line
24      common /rdwr/   iun1,iun2,iun3,iun4,iun5
25      data iel/'bq',
26     &         'h ', 'he',
27     &         'li', 'be', 'b ', 'c ', 'n ', 'o ', 'f ', 'ne',
28     &         'na', 'mg', 'al', 'si', 'p ', 's ', 'cl', 'ar',
29     &         'k ', 'ca',
30     &                     'sc', 'ti', 'v ', 'cr', 'mn',
31     &                     'fe', 'co', 'ni', 'cu', 'zn',
32     &                     'ga', 'ge', 'as', 'se', 'br', 'kr',
33     & 'rb','sr','y ','zr','nb','mo','tc','ru','rh','pd','ag','cd',
34     & 'in','sn','sb','te','i ','xe','cs','ba','la','ce','pr','nd',
35     & 'pm','sm','eu','gd','tb','dy','ho','er','tm','yb','lu','hf',
36     & 'ta','w ','re','os','ir','pt','au','hg','tl','pb','bi','po',
37     & 'at','rn','fr','ra','ac','th','pa','u ','np','pu','am','cm',
38     & 'bk','cf','x '/
39
40      istatz = 0
41      ifnd = 0
42      irwd = 0
43      toang = 0.52917706d0
44
455     continue
46
47      first = .true.
48      if (izo.eq.1) then
49          nzz = nz
50      else
51          nz = 0
52      endif
53      maxnz = mxzat
54      ottest=.true.
55      heat = 0.0d0
56
5710    continue
58      zmstar = .false.
59      nztt = 0
60      do while (getlin(0).eq.1)
61          if (icdex(line,'[AMB').ne.0) return
62          i1 = index(line,'(')
63          i2 = index(line,')')
64          if (i1.ne.0.and.i2.ne.0.and.i2.gt.i1) then
65               line = line(1:i1-1)//line(i2+1:)
66          endif
67          if (index(line,'VARIABLE        FUNCTION').ne.0) then
68               nword = 3
69               if (index(line,'SECOND').ne.0) nword = 4
70               if (getlin(0).eq.0) goto 10
71               do ii=1,nword
72                   ktype = nxtwrd(string,nstr,itype,rtype)
73               end do
74               if (ktype.eq.3) then
75                   heat = rtype
76               endif
77               goto 10
78          endif
79          if (index(line,'HEAT OF FORMATION').ne.0) then
80             do while(.true.)
81                ktype = nxtwrd(string,nstr,itype,rtype)
82                if (ktype.eq.3) heat = rtype
83                if (ktype.eq.0.or.ktype.eq.3) goto 10
84             end do
85          endif
86          if (index(line,'POINT   POTENTIAL').ne.0.or.
87     &        index(line,'POINT  POTENTIAL').ne.0) then
88               in = 2
89               if (index(line,'FEMTOSECONDS').ne.0) in = 3
90               if (getlin(0).eq.0) goto 10
91c               do ii=1,4
92c Armando Juan Navarro Vazquez says total energy not interresting
93c Use potential energy
94c
95               do ii=1,in
96                   ktype = nxtwrd(string,nstr,itype,rtype)
97               end do
98               if (ktype.eq.3) heat = rtype
99               goto 10
100          endif
101          if (index(line,'FINAL GEOMETRY OBTAINED').ne.0) then
102             call redel(line,4)
103             coplin = line
104             iws = 0
105             ktype = -1
106             do while (ktype.ne.0)
107                ktype = nxtwrd(string,nstr,itype,rtype)
108                if (ktype.ne.0) iws = iws + 1
109             end do
110             line = coplin
111             if (iws.eq.7) then
112               natoms = 0
113               do while(.true.)
114                 ktype = nxtwrd(string,nstr,itype,rtype)
115                 if (ktype.eq.1.and.(nstr.eq.1.or.nstr.eq.2)) then
116                      if (nstr.eq.1) then
117                         catomt(1:1) = string(1:1)
118                         catomt(2:2) = ' '
119                      else
120                         catomt = string(1:2)
121                      endif
122                 else
123                      if (natoms.ne.0) then
124                         istatz = 2
125                         return
126                      endif
127                      goto 100
128                 endif
129
130                 iatmp = 0
131                 catom = tolowf(catomt)
132                 do i=1,100
133                    if ( catom .eq. iel(i) ) iatmp = i - 1
134                 end do
135                 if (catom.eq.'xx') iatmp = 99
136                 if (catom.eq.'tv') iatmp = 99
137                 if (iatmp.eq.0.or.iatmp.gt.100) goto 100
138
139                 natoms = natoms + 1
140                 ianc(natoms) = iatmp
141
142                 do ii=1,3
143                    ktype = nxtwrd(string,nstr,itype,rtype)
144                    if (ktype.eq.3) then
145                       cc(ii,natoms) = rtype / toang
146                    else
147                       goto 100
148                    endif
149                    ktype = nxtwrd(string,nstr,itype,rtype)
150                    if (ktype.ne.2)  goto 100
151                 end do
152                 if (getlin(0).ne.1) return
153               end do
154             endif
155          endif
156          if (index(line,'scf done: ').ne.0) then
157             do while(.true.)
158                ktype = nxtwrd(string,nstr,itype,rtype)
159                if (ktype.eq.3) heat = rtype
160                if (ktype.eq.0.or.ktype.eq.3) goto 10
161             end do
162          endif
163          do i=1,3
164             idm(i) = 0
165             zzz(i) = 0.0d0
166             izt(i) = 0
167          end do
168          ktype = nxtwrd(string,nstr,itype,rtype)
169          if (ktype.eq.2) then
170             if (itype.eq.1.and.first) zmstar = .true.
171             if (zmstar) then
172                 ktype = nxtwrd(string,nstr,itype,rtype)
173             else
174                 goto 100
175             endif
176          endif
177c element label
178          if (ktype.eq.1.and.(nstr.eq.1.or.nstr.eq.2)) then
179               if (nstr.eq.1) then
180                  catomt(1:1) = string(1:1)
181                  catomt(2:2) = ' '
182               else
183                  catomt = string(1:2)
184               endif
185          else
186               goto 100
187          endif
188
189c check which zmstar type
190
191          if (zmstar.and.nztt.eq.0) then
192               izmtyp = 0
193               if (gnreal(tmp,3,.false.)) izmtyp = 1
194          endif
195
196          num = nztt
197          if (num.gt.3) num = 3
198          if (.not.zmstar.or.(zmstar.and.izmtyp.eq.1)) num = 3
199
200          getnxt = .true.
201          do i=1,num
202             if (getnxt) ktype = nxtwrd(string,nstr,itype,rtype)
203
204c bl,aplh,bet
205
206             getnxt = .true.
207             if (ktype.eq.2.or.ktype.eq.3) then
208                if (ktype.eq.2) then
209                   zzz(i) = 1.0d0*itype
210                else
211                   zzz(i) = rtype
212                endif
213             else
214                if (zmstar) then
215                   if (ktype.eq.0) then
216c                     read in connectivity
217                      goto 200
218                   else
219                      goto 100
220                   endif
221                else
222                   goto 100
223                endif
224             endif
225
226c optimize flag 1-3
227
228             ktype = nxtwrd(string,nstr,itype,rtype)
229             if (zmstar) then
230                if (ktype.eq.1) then
231                   if (nstr.eq.1) then
232                      if (string(1:1).eq.'*'.or.
233     &                    string(1:1).eq.'+') then
234                          idm(i) = 1
235                      else
236                          goto 100
237                      endif
238                   else
239                      goto 100
240                   endif
241                else
242                   getnxt = .false.
243                   idm(i) = 0
244                endif
245             else
246                if (ktype.eq.2) then
247                   idm(i) = abs(itype)
248                else
249                   goto 100
250                endif
251             endif
252          end do
253200       continue
254c
255c connectivity
256c
257          num = nztt
258          if (num.gt.3) num = 3
259          if (.not.zmstar) num = 3
260
261          do i=1,num
262             if (getnxt) ktype = nxtwrd(string,nstr,itype,rtype)
263             getnxt = .true.
264             if (ktype.eq.2) then
265                izt(i) = itype
266             endif
267          end do
268
269c check first three atoms
270
271          if (izo.ne.1) then
272             ioke = 1
273             if (nztt.le.2) then
274                do i=1,3
275                   iztt(i,nztt+1) = izt(i)
276                end do
277             endif
278             if (nztt.le.2) then
279                if (iztt(1,1).ne.0) ioke = 0
280                if (iztt(2,1).ne.0) ioke = 0
281                if (iztt(3,1).ne.0) ioke = 0
282             endif
283             if (nztt.ge.1.and.nztt.le.2) then
284                if (iztt(1,2).le.0) ioke = 0
285                if (iztt(2,2).ne.0) ioke = 0
286                if (iztt(3,2).ne.0) ioke = 0
287             endif
288             if (nztt.eq.2) then
289                if (iztt(1,3).le.0) ioke = 0
290                if (iztt(2,3).le.0) ioke = 0
291                if (iztt(3,3).ne.0) ioke = 0
292             endif
293             if (ioke.ne.1) then
294                istatz = 0
295                first = .true.
296                goto 5
297             endif
298          endif
299
300c check atom4 and upwards
301
302          iztot = 0
303          do i=1,3
304             iztot = iztot + izt(i)
305          end do
306          if (nztt.gt.4.and.iztot.eq.0) then
307             istatz = 0
308             first = .true.
309             goto 5
310          endif
311
312300       continue
313          if (zmstar.and.heat.eq.0.0d0.and.ifnd.eq.0) then
314             ifnd = 1
315             goto 100
316          endif
317c
318c copy temporary variables into z-mat
319c
320          iatmp = 0
321          catom = tolowf(catomt)
322          do i=1,100
323             if ( catom .eq. iel(i) ) iatmp = i - 1
324          end do
325          if (catom.eq.'xx') iatmp = 99
326          if (catom.eq.'tv') iatmp = 99
327          if (iatmp.eq.0.or.iatmp.gt.100) goto 100
328          first = .false.
329          if (nztt.ge.1) istatz = 1
330          nz = nz + 1
331          nztt = nztt + 1
332          ianz(nz) = iatmp
333          do i=1,3
334             iz(i,nz) = izt(i)
335          end do
336
337          iz(4,nz) = 0
338          bl(nz) = zzz(1)
339          ibl(nz) = idm(1)
340          alph(nz) = zzz(2)
341          ialph(nz) = idm(2)
342          bet(nz) = zzz(3)
343          ibet(nz) = idm(3)
344      end do
345      first = .false.
346
347100   if (first) goto 10
348      if (izo.eq.1) then
349          if (istatz.eq.0) nz = nzz
350          return
351      endif
352      if (istatz.eq.1) then
353c          call prtzm
354          ihaszm = 1
355          call stocc(maxnz,nz,iz,bl,alph,bet,alpha,beta,ierror)
356          if (ierror.lt.5) then
357             call stoc(maxnz,nz,0,0,0,ianz,iz,bl,alph,bet,
358     &       ottest,natoms,ian,c,cz,imap,alpha,beta,
359     &       oerror,.true.,.true.)
360             if (oerror) then
361                istatz = 0
362                return
363             endif
364          else
365             istatz = 0
366             return
367          endif
368          do i=1,natoms
369             do j=1,3
370                cc(j,i) = c(j,i)
371             end do
372             ianc(i) = ian(i)
373          end do
374      else
375          if (ifnd.eq.1.and.irwd.eq.0.and.ido.eq.1) then
376             irwd = 1
377             call rewfil
378             goto 5
379          endif
380      endif
381
382      return
383      end
384