1      subroutine getpod(icom,ifd,iff,doscl,ioatms,ioadd,
2     &                  coo,ianz,iaton,iconn,
3     &                  scal,scali,smag,natc,ichx,icrtp)
4
5      implicit double precision (a-h,p-x),integer (i-n),logical (o)
6      parameter (numat1=20000)
7      parameter (mxcon=10)
8      parameter (mxel=100)
9      parameter (mxonh=100)
10      parameter (numatm=2000)
11      parameter (small=1.0d-4)
12      common /pnthlp/ ipoints,ipnt
13      common /athlp/ iatoms, mxnat
14      common /gamori/imap(numat1),cstand(3,numat1),istand(numat1),
15     &               msucc
16      common /oniomh/ ioni,nion,ionih(mxonh),natonh(mxonh),
17     &                iomap(numatm),icntat(mxonh),fct(mxonh),
18     &                xyzi(3,mxonh)
19      common /gauori/ nzm,nso,nio,nzo,ioropt,ifor,
20     &                ixyz98,iopr,isymm,irc,imp2,icntp,itd
21      common /align/  vecs(3,3),nscnd,iscst,ialtyp,iocnt
22      common /zmpart/ ipart,imn,imx,idcur
23      common /rdwr/   iun1,iun2,iun3,iun4,iun5
24      common /geocnv/ fdum3(10),idum(7),ieav,ifrav,mxpnt
25      common /mopver/ mopopt
26      common /animo/  movie
27
28      character lstr*137
29      character*6 astr
30      character*8 refcod
31      integer first,prev,next,last,getmf,doscl
32      logical dummy,dozme
33      logical ctoz,molpot,elpot,chpot
34      common /choic/  iftyp,isbin,ctoz,molpot,elpot,chpot
35      common /getpnt/ irtype,ipdbon,ipdbgro,ifav,ioxyz,
36     &                iconv,ircus,dozme
37      common /slagau/ ihasd,isgau,ido5d,ido7f,ido9g,ihasg
38      common /cllab/  iclon,iclpnt(4)
39      common /ecce/iecce
40      character*23 stmp
41      dimension dum(3),coo(3,*),ianz(*),iaton(*),iconn(mxcon+1,*)
42
43c iftyp      file type
44c
45c 1          mopac
46c 2          gamess
47c 3          gamus
48c 4          gauss
49c 5          adfin
50c 6          chemx
51c 7          cpmd
52c 8          qchem
53c 9          orca
54c 10         xyz
55c 11         tinker + small
56c 12         tinker protien
57c 13         mol file
58c 15         nwchem
59
60cd     write(iun3,'(a)')'enter subroutine getpoi'
61
62      first = -1
63      next  = -2
64      last  = -3
65      prev  = -4
66      tol = small
67      idoprt = 0
68
69      icomm = icom
70
71c     parse variables to C used in getpoi
72
73      call parpoi(nzm,nso,nio,nzo,ioropt,ifor,ixyz98,iopr,isymm,irc,
74     &            imp2,icntp,msucc,ioni,mopopt,isbin,irtype,
75     &            ipdbgro,ifav,ioxyz,iconv,ircus,nscnd,iscst,ialtyp)
76
77c for gaussian
78
79      stmp = 'Input orientation:'
80      ns = 18
81      if (ixyz98.gt.0.and.(ioropt.eq.2)) then
82          stmp = 'Standard orientation:'
83          ns = 22
84      endif
85
86      if (dozme) then
87c since dumzm resets the ipart variable, we need the help variable dopart
88c this is because otherwise conpdb downstream goes wrong
89          if (ipart.eq.1) idoprt = ipart
90          ifav = 0
91          call dumzm(coo,ianz,iatoms)
92          goto 100
93      endif
94
95c     mopac binary .gpt or GRAPH file (msf excluded by
96c                                      iftyp.ne.6: .not.chemx)
97
98      if (isbin.eq.1.and.iftyp.ne.6) iatoms = 0
99
100c     convert mopac coordinates to a.u.
101c     vdwr are in angstrom, vrad in a.u.
102
103      if (iftyp.eq.1.and.iconv.eq.1) call xyzcoo(1,1,0)
104
105c     CPMD
106
107      if (iftyp.eq.7.and.ipoints.eq.0.and.
108     &    (ioxyz.eq.1.or.irtype.eq.4)) then
109          call xyzcoo(1,0,0)
110          call haszm(.false.)
111      endif
112
113c     GAMESS xyz optimise or FREQ job
114
115      if (iftyp.eq.2.and.ipoints.eq.0.and.
116     &    (ioxyz.eq.1.or.irtype.eq.4)) then
117          call xyzcoo(1,0,0)
118          call haszm(.false.)
119      endif
120
121c     GAUSSIAN single point or FREQ job, the first point
122
123      if (iftyp.eq.4.and.
124     &    (ipoints.le.1.or.irtype.eq.4).and.icomm.eq.-1) then
125          call rewfil
126          call gaupoi(istat)
127          if (istat.eq.-1) then
128             ifav = 0
129          else
130             ifav = 1
131          endif
132          if ((istat.eq.0.or.irtype.eq.4).and.ioni.ne.1)
133     &          call xyzcoo(1,0,0)
134          if (istat.eq.0) call haszm(.false.)
135      endif
136
137      if (icomm.eq.prev) then
138          if ( ipnt.le.1 ) return
139          ipnt = ipnt -1
140          icomm = ipnt
141      endif
142
143      ifav = 0
144
145c     NOT {GAMESS-UK,GAMESS-US,GAUSSIAN or CPMD}
146c     AND iatoms=0: EXIT
147
148      if (iatoms.eq.0.and..not.
149     &    ((iftyp.ge.2.and.iftyp.le.4).or.iftyp.eq.7)) then
150          call inferr('No coordinates found !',0)
151cd         write(iun3,'(a)')'leave subroutine getpoi'
152          return
153      endif
154
155c Mopac2007 AUX file
156
157      if (iftyp.eq.1.and.mopopt.eq.4.and.ipoints.ge.1) then
158
159          if (icomm .eq. first) then
160             call rewfil
161             call search(lstr,'HEAT_OF_FORM_UPDATED',istat)
162             call redel(lstr,2)
163             if (istat.eq.1) call gtapnt(iatoms,coo,istat)
164             if (istat.eq.1) ipnt = 1
165          elseif (icomm .eq. next) then
166             call search(lstr,'HEAT_OF_FORM_UPDATED',istat)
167             call redel(lstr,2)
168             if (istat.eq.1) call gtapnt(iatoms,coo,istat)
169             if (istat.eq.1) then
170                ipnt = ipnt +1
171             else
172                call rewfil
173                call search(lstr,'HEAT_OF_FORM_UPDATED',istat)
174                call redel(lstr,2)
175                if (istat.eq.1) call gtapnt(iatoms,coo,istat)
176                if (istat.eq.1) ipnt = 1
177             endif
178          elseif (icomm.gt.0) then
179             call rewfil
180             do i=1,icomm
181                call search(lstr,'HEAT_OF_FORM_UPDATED',istat)
182                call redel(lstr,2)
183                if (istat.eq.1) call gtapnt(iatoms,coo,istat)
184                if (istat.eq.1) then
185                    dummy = .true.
186                else
187                    dummy = .false.
188                endif
189             end do
190             ipnt = icomm
191          endif
192c         convert coo from angs to au
193          call xyzcoo(0,1,0)
194      endif
195
196      if (iftyp.eq.1.and.iconv.eq.0.and.ipoints.ge.1) then
197          if (icomm .eq. first) then
198             call rewfil
199             call getmop(iatoms,heat,0,1,istat)
200             if (istat.eq.1) ipnt = 1
201          elseif (icomm .eq. next) then
202             call getmop(iatoms,heat,0,0,istat)
203             if (istat.eq.1) then
204                ipnt = ipnt +1
205             else
206                call rewfil
207                call getmop(iatoms,heat,0,0,istat)
208                if (istat.eq.1) ipnt = 1
209             endif
210          elseif (icomm.gt.0) then
211             call rewfil
212             idomopf = 1
213             do i=1,icomm
214                call getmop(iatoms,heat,0,idomopf,istat)
215                if (istat.eq.1) then
216                    dummy = .true.
217                else
218                    dummy = .false.
219                endif
220                idomopf = 0
221             end do
222             ipnt = icomm
223          endif
224      endif
225
226      if ((iftyp.ge.10.or.(iftyp.eq.5.and.ihasg.ge.1)).and.
227     &    ipoints.gt.1) then
228
229          if (icomm .eq. first) then
230
231             if (iftyp.eq.5) then
232                call rewmf
233                call srchmf(lstr,'[GEOMETRIES]',istat)
234             else
235                call rewfil
236             endif
237
238             if (iftyp.eq.5.and.ihasg.eq.2) then
239                call getzm(iatoms,0,0,istat)
240                if (istat.eq.1) ipnt = 1
241             else
242                if (iftyp.eq.11.or.iftyp.eq.12) then
243                   call gettnk(igttnk,0,ipdbon,iff,iheat,heat)
244                   if (igttnk.eq.1) ipnt = 1
245                elseif (iftyp.eq.13) then
246                   call getmol(igetmo,0)
247                   if (igetmo.eq.1) ipnt = 1
248                else
249                   call getxyz(igetxy,heat,0)
250                   if (igetxy.eq.1) ipnt = 1
251                endif
252             endif
253
254             if (iftyp.eq.5) then
255                iolin = getmf()
256                call rdfc(ipnt,istat)
257                if (istat.eq.1) then
258                   ifav = 1
259                else
260                   ifav = 0
261                endif
262                call putmf(iolin)
263             endif
264
265          elseif (icomm .eq. next.and..not.
266     &            (iftyp.eq.5.and.ifrav.eq.1)) then
267             if (iftyp.eq.5.and.ihasg.eq.2) then
268                call getzm(iatoms,0,0,istat)
269                if (istat.eq.1) then
270                   dummy = .true.
271                else
272                   dummy = .false.
273                endif
274             else
275                if (iftyp.eq.11.or.iftyp.eq.12) then
276                   call gettnk(igttnk,0,ipdbon,iff,iheat,heat)
277                   if (igttnk.eq.1) then
278                      dummy = .true.
279                   else
280                      dummy = .false.
281                   endif
282                elseif (iftyp.eq.13) then
283                   call getmol(igetmo,0)
284                   if (igetmo.eq.1) then
285                      dummy = .true.
286                   else
287                      dummy = .false.
288                   endif
289                else
290                   call getxyz(igetxy,heat,0)
291                   if (igetxy.eq.1) then
292                       dummy = .true.
293                   else
294                       dummy = .false.
295                   endif
296                endif
297             endif
298             if (dummy) then
299                ipnt = ipnt +1
300             else
301                if (iftyp.eq.5) then
302                   call rewmf
303                   call srchmf(lstr,'[GEOMETRIES]',istat)
304                else
305                   call rewfil
306                endif
307                if (iftyp.eq.5.and.ihasg.eq.2) then
308                   call getzm(iatoms,0,0,istat)
309                   if (istat.eq.1) then
310                      dummy = .true.
311                   else
312                      dummy = .false.
313                   endif
314                else
315                   if (iftyp.eq.11.or.iftyp.eq.12) then
316                      call gettnk(igttnk,0,ipdbon,iff,iheat,heat)
317                      if (igttnk.eq.1) then
318                         dummy = .true.
319                      else
320                         dummy = .false.
321                      endif
322                   elseif (iftyp.eq.13) then
323                      call getmol(igetmo,0)
324                      if (igetmo.eq.1) then
325                         dummy = .true.
326                      else
327                         dummy = .false.
328                      endif
329                   else
330                      call getxyz(igetxy,heat,0)
331                      if (igetxy.eq.1) then
332                          dummy = .true.
333                      else
334                          dummy = .false.
335                      endif
336                   endif
337                endif
338                if (dummy) ipnt = 1
339             endif
340          elseif (icomm.gt.0.or.
341     &           (icomm.eq.next.and.iftyp.eq.5.and.ifrav.eq.1)) then
342             if (icomm.eq.next) then
343                ipntt = ipnt + 1
344                if (ipntt.gt.ipoints) ipntt = ipoints
345             endif
346             if (icomm.gt.0) ipntt = icomm
347             if (iftyp.eq.5) then
348                call rewmf
349                call srchmf(lstr,'[GEOMETRIES]',istat)
350             else
351                call rewfil
352             endif
353             do i=1,ipntt
354                if (iftyp.eq.5.and.ihasg.eq.2) then
355                   call getzm(iatoms,0,0,istat)
356                   if (istat.eq.1) then
357                      dummy = .true.
358                   else
359                      dummy = .false.
360                   endif
361                else
362                   if (iftyp.eq.11.or.iftyp.eq.12) then
363                      if (iff.eq.7) then
364
365c                        speed things up ambfor MD trajectories
366
367                         if (i.ne.ipntt) then
368                            call gtheat(igttnk,iheat,heat)
369                         else
370                            call search(lstr,'[AMBFOR]',istat)
371                            if (istat.eq.1) then
372                                call bckfil
373                                call gettnk(
374     &                               igttnk,0,ipdbon,iff,iheat,heat)
375                            else
376                                igttnk = 0
377                            endif
378                         endif
379                      else
380                         call gettnk(igttnk,0,ipdbon,iff,iheat,heat)
381                      endif
382                      if (igttnk.eq.1) then
383                         dummy = .true.
384                      else
385                         dummy = .false.
386                      endif
387                   elseif (iftyp.eq.13) then
388                      call getmol(igetmo,0)
389                      if (igetmo.eq.1) then
390                         dummy = .true.
391                      else
392                         dummy = .false.
393                      endif
394                   else
395                      call getxyz(igetxy,heat,0)
396                      if (igetxy.eq.1) then
397                          dummy = .true.
398                      else
399                          dummy = .false.
400                      endif
401                   endif
402                endif
403             end do
404             if (iftyp.eq.5) then
405                call rdfc(ipntt,istat)
406                if (istat.eq.1) then
407                   ifav = 1
408                else
409                   ifav = 0
410                endif
411             endif
412             ipnt = ipntt
413          endif
414      endif
415
416      if (iftyp.eq.8) then
417          if ( icomm .eq. first ) then
418             call rewmf
419             ipnt = 1
420          elseif ( icomm .eq. next ) then
421             if ( ipnt.ge.ipoints ) return
422             ipnt = ipnt +1
423          elseif ( icomm .eq. last ) then
424             call rewmf
425             ipnt = ipoints
426          elseif (icomm.gt.0) then
427             call rewmf
428             ipnt = icomm
429          endif
430          call qcxyz(0,.false.,ipnt,istat)
431          if (istat.eq.-1) then
432               ifav = 0
433          else
434               ifav = 1
435          endif
436      endif
437
438      if (iftyp.eq.9) then
439          if ( icomm .eq. first ) then
440             call rewmf
441             ipnt = 1
442          elseif ( icomm .eq. next ) then
443             if ( ipnt.ge.ipoints ) return
444             ipnt = ipnt +1
445          elseif ( icomm .eq. last ) then
446             call rewmf
447             ipnt = ipoints
448          elseif (icomm.gt.0) then
449             call rewmf
450             ipnt = icomm
451          endif
452          call orcxyz(0,ipnt,istat)
453          if (istat.eq.-1) then
454               ifav = 0
455          else
456               ifav = 1
457          endif
458      endif
459
460      if (iftyp.eq.15) then
461          if ( icomm .eq. first ) then
462             if (iecce.eq.1) then
463                call rewfil
464             else
465                call rewmf
466             endif
467             ipnt = 1
468          elseif ( icomm .eq. next ) then
469             if ( ipnt.ge.ipoints ) return
470             ipnt = ipnt +1
471          elseif ( icomm .eq. last ) then
472             if (iecce.eq.1) then
473                call rewfil
474             else
475                call rewmf
476             endif
477             ipnt = ipoints
478          elseif (icomm.gt.0) then
479             if (iecce.eq.1) then
480                call rewfil
481             else
482                call rewmf
483             endif
484             ipnt = icomm
485          endif
486          if (iecce.eq.1) then
487             call enwxyz(0,ipnt,istat)
488          else
489             call nwxyz(0,ipnt,istat)
490          endif
491      endif
492
493      if (iftyp.eq.2.and.ipoints.gt.0) then
494          if ( icomm .eq. first ) then
495             call rewfil
496             ipnt = 1
497          elseif ( icomm .eq. next ) then
498             if ( ipnt.ge.ipoints ) return
499             ipnt = ipnt +1
500          elseif ( icomm .eq. last ) then
501             call rewfil
502             ipnt = ipoints
503          elseif (icomm.gt.0) then
504             call rewfil
505             ipnt = icomm
506          endif
507          if (ioxyz.eq.1) then
508             itmp = ipnt - 1
509             call gampoi(itmp,istat,ioxyz)
510          else
511             call gampoi(ipnt,istat,ioxyz)
512          endif
513          if ( istat .eq. -1 ) then
514               ifav = 0
515          else
516               ifav = 1
517          endif
518c          if (msucc.eq.0.and.ioxyz.eq.0) ifav = 0
519      endif
520
521      if (iftyp.eq.4.and.ipoints.gt.1) then
522          if ( icomm .eq. first ) then
523             call rewfil
524             ipnt = 1
525          elseif ( icomm .eq. next ) then
526             if ( ipnt.ge.ipoints ) return
527             ipnt = ipnt +1
528          elseif ( icomm .eq. last ) then
529             call rewfil
530             do i=1,ipoints
531                 if ((nzm.eq.nzo.or.2*nzm.eq.nzo).and.nzm.gt.0) then
532                    call search(lstr,'Z-MATRIX (ANGSTROMS',istat)
533                 elseif (nzo.gt.0) then
534                    call search(lstr,'Z-Matrix orientation:',istat)
535                 else
536                    call search(lstr,stmp(1:ns),istat)
537                 endif
538             end do
539             call bckfil
540             ipnt = ipoints
541          elseif (icomm.gt.0) then
542             call rewfil
543             do i=1,icomm
544                 if ((nzm.eq.nzo.or.2*nzm.eq.nzo).and.nzm.gt.0) then
545                    call search(lstr,'Z-MATRIX (ANGSTROMS',istat)
546                 elseif (nzo.gt.0) then
547                    call search(lstr,'Z-Matrix orientation:',istat)
548                 else
549                    call search(lstr,stmp(1:ns),istat)
550                 endif
551             end do
552             call bckfil
553             ipnt = icomm
554          endif
555          call gaupoi(istat)
556          if ( istat .eq. -1 ) then
557               ifav = 0
558          else
559               ifav = 1
560          endif
561      endif
562
563      if (iftyp.eq.3.and..not.(ipoints.le.1.and.iatoms.gt.0) ) then
564          if ( icomm .eq. first ) then
565             call rewfil
566             ipnt = 1
567             call gamupt(ipnt,istat)
568          elseif ( icomm .eq. next ) then
569             if ( ipnt.ge.ipoints ) return
570             ipnt = ipnt +1
571             call gamupt(ipnt,istat)
572          elseif ( icomm .eq. last ) then
573             call rewfil
574             do i=1,ipoints
575                call gamupt(ipoints,istat)
576             end do
577             call bckfil
578             ipnt = ipoints
579          elseif (icomm.gt.0) then
580             call rewfil
581             do i=1,icomm
582                call gamupt(icomm,istat)
583             end do
584             call bckfil
585             ipnt = icomm
586          endif
587          ifav = 1
588      endif
589
590      if (iftyp.eq.7.and..not.(ipoints.le.1.and.iatoms.gt.0)
591     &     .and.irtype.ne.3 ) then
592c          iunt = iun2
593c          iun2 = iun5
594          if ( icomm .eq. first ) then
595             call rewfil
596             ipnt = 1
597             call cpmdpt(ipnt,istat)
598          elseif ( icomm .eq. next ) then
599             if ( ipnt.ge.ipoints ) return
600             ipnt = ipnt +1
601             call cpmdpt(ipnt,istat)
602          elseif ( icomm .eq. last ) then
603             call rewfil
604             do i=1,ipoints
605                call cpmdpt(ipoints,istat)
606             end do
607             call bckfil
608             ipnt = ipoints
609          elseif (icomm.gt.0) then
610             call rewfil
611             do i=1,icomm
612                call cpmdpt(icomm,istat)
613             end do
614             call bckfil
615             ipnt = icomm
616          endif
617          ifav = 1
618c          iun2 = iunt
619      endif
620
621      if (iftyp.eq.7.and..not.(ipoints.le.1.and.iatoms.gt.0)
622     &     .and.irtype.eq.3 ) then
623          if (icomm.ne.first) then
624             iunt = iun2
625             iun2 = iun5
626          endif
627          if ( icomm .eq. first ) then
628             call rewfil
629             ipnt = 1
630             call cpmdpt(ipnt,istat)
631          elseif ( icomm .eq. next ) then
632             if ( ipnt.ge.ipoints ) return
633             ipnt = ipnt +1
634             call cpmdptdyn(ipnt,istat)
635          elseif ( icomm .eq. last ) then
636             call rewfil
637             do i=1,ipoints
638                call cpmdptdyn(ipoints,istat)
639             end do
640             call bckfil
641             ipnt = ipoints
642          elseif (icomm.gt.0) then
643             call rewfil
644             do i=1,icomm
645                call cpmdptdyn(icomm,istat)
646             end do
647             call bckfil
648             ipnt = icomm
649          endif
650          ifav = 1
651          if (icomm.ne.first) iun2 = iunt
652      endif
653
654      if (((icrtp.le.3.and.icrtp.ge.1).or.icrtp.eq.5)
655     &     .and.ipoints.gt.1) then
656          if (icomm.eq.first) then
657             call rewfil
658             if (icrtp.eq.3) then
659                call rfbio(0,1,istat)
660             elseif (icrtp.eq.1) then
661                call rdchx(0,3,0,0,0,istat,1)
662             elseif (icrtp.eq.5) then
663                call getxdt(1,natc,istat)
664             else
665                call rfdat(0,istat,refcod)
666             endif
667             ipnt = 1
668          elseif ( icomm .eq. next ) then
669             if ( ipnt.ge.ipoints ) return
670             ipnt = ipnt +1
671             if (icrtp.eq.3) then
672                call rfbio(0,0,istat)
673             elseif (icrtp.eq.1) then
674                call rdchx(0,3,0,0,0,istat,1)
675             elseif (icrtp.eq.5) then
676                call getxdt(ipnt,natc,istat)
677             else
678                call rfdat(0,istat,refcod)
679             endif
680          elseif (icomm.gt.0) then
681             call rewfil
682             if (icrtp.eq.3) then
683                call rfbio(0,1,istat)
684             elseif (icrtp.eq.1) then
685                call rdchx(0,3,0,0,0,istat,1)
686             elseif (icrtp.eq.5) then
687                call getxdt(icomm,natc,istat)
688             else
689                call rfdat(0,istat,refcod)
690             endif
691             if (icrtp.ne.5) then
692                do i=1,icomm-1
693                   if (icrtp.eq.3) then
694                      call rfbio(0,0,istat)
695                   elseif (icrtp.eq.1) then
696                      call rdchx(0,3,0,0,0,istat,1)
697                   else
698                      call rfdat(0,istat,refcod)
699                   endif
700                end do
701             endif
702             ipnt = icomm
703          endif
704          if (.not.(icrtp.eq.3.and.istat.eq.1))
705     &       call fdat(ifd,0,0,0,0,0)
706      endif
707
708      if (ipdbon.eq.1.and.ipdbgro.eq.1) then
709          if (icomm.eq.first) then
710             ipnt = 1
711             call gtfrm(ipnt)
712          elseif ( icomm .eq. next ) then
713             if ( ipnt.ge.ipoints ) return
714             ipnt = ipnt + 1
715             call gtfrm(ipnt)
716          elseif (icomm.gt.0) then
717             call rewfil
718             call gtfrm(icomm)
719             ipnt = icomm
720          endif
721      endif
722
723100   continue
724
725      if (ioadd.eq.0) then
726        if (.not.(
727     &    ( (iftyp.eq.2.or.iftyp.eq.3).and.ipoints.eq.0
728     &      .and.(ioxyz.eq.1.or.irtype.eq.4) ).or.
729     &    (iftyp.eq.4.and.ipoints.le.1).or.
730     &    (icrtp.eq.3.and.ipoints.gt.1).or.
731     &    (ipdbon.eq.1.and.dozme))) call docent
732      endif
733
734      if (iftyp.eq.2.and.ipoints.gt.0.and.ioxyz.eq.0.and.msucc.eq.1
735     &    .and.(.not.dozme))
736     &    call rotmom(ipnt,ifav)
737
738
739      ibeg = 1
740      if (ioadd.eq.1) ibeg = ioatms + 1
741
742      do i=ibeg,iatoms
743         if (ipdbgro.eq.1) then
744            if (ianz(i).eq.100.and.ipnt.gt.1) then
745               if (i.gt.iclpnt(1)+7) iaton(i) = 0
746            else
747               iaton(i) = 1
748            endif
749         else
750            iaton(i) = 1
751         endif
752      end do
753
754      if (ipdbon.eq.1) then
755         if (ialtyp.eq.1) then
756            if (dozme) then
757               call convar(coo,iconn,ianz,nscnd,iscst)
758               call pmfass(1,0)
759               call totpmf(dum(1))
760               call upsco()
761            endif
762         else
763            if (dozme) then
764               if (idoprt.eq.0) call conpdb
765               call chkbck(0)
766            else
767               if (((iftyp.eq.11.or.iftyp.eq.12).and.iff.eq.7)
768     &             .or.ipdbgro.eq.1) call chkbck(1)
769            endif
770         endif
771      else
772         if ((.not.(iftyp.eq.6.or.ichx.ne.0).and.
773     &        .not.(iftyp.eq.5.and.ihasg.eq.1).and.
774     &        .not.(iftyp.ge.10.and.iftyp.lt.15)
775     &       ).or.(dozme.and..not.ichx.eq.1)) call doconn
776      endif
777
778      if (ipdbon.eq.0) then
779         if (.not.ichx.eq.1.and.iff.ne.7) call dohcon(0)
780      else
781         call domcon(1,1)
782      endif
783c
784c  get scale
785c
786
787      if (ioadd.eq.0) then
788         if (doscl.eq.1.and..not.(icrtp.eq.3.and.ipoints.gt.1).and.
789     &       .not.(ipdbon.eq.1.and.dozme)) then
790            call doscal
791            if (dozme.and.iatoms.lt.7) scali = 3.5d0
792            scal = scali * 2.4d0 * smag
793         endif
794      endif
795
796      if (ifav.eq.1) call parfc
797
798      call drwgeo
799
800      if (movie.eq.0.or.(movie.ne.0.and.ipnt.eq.ipoints)) then
801         if (ctoz) call allzmt(ipdbon)
802         if (.not.dozme) call upzme
803      endif
804
805      dozme = .false.
806
807
808      if (icrtp.eq.2.and.ipoints.gt.1) then
809          lstr = refcod
810          call inferr(lstr,0)
811cd         write(iun3,'(a)')'leave subroutine getpoi'
812          return
813      elseif (icrtp.eq.5.and.ipoints.gt.1) then
814          call zerstr(ipnt,astr,6,1)
815          lstr = 'point '//astr
816          call inferr(lstr,0)
817          return
818      endif
819
820      if (.not.( iftyp.eq.7.or.
821     &          (iftyp.ge.2.and.iftyp.le.4).or.
822     &          (iftyp.eq.1.and.iconv.eq.0).or.
823     &          (iftyp.ge.10.and.ipoints.gt.1)
824     &         ).or.dozme) then
825cd        write(iun3,'(a)')'leave subroutine getpoi'
826         return
827      endif
828
829      if (ipoints.eq.0) then
830         lstr = ' '
831         if (iftyp.ne.1) then
832             lstr = 'Single point'
833         endif
834      elseif (ipnt.eq.1) then
835         lstr = 'First point'
836      elseif (ipnt.eq.ipoints) then
837         lstr = 'Last point'
838      else
839         if (ioxyz.eq.1) then
840            call zerstr(ipnt-1,astr,6,1)
841         else
842            call zerstr(ipnt,astr,6,1)
843         endif
844         lstr = 'point '//astr
845      endif
846
847      if (iftyp.eq.7) lstr=lstr(1:22)//'CPMD output file'
848
849      call inferr(lstr,0)
850
851cd     write(iun3,'(a)')'leave subroutine getpoi'
852      return
853      end
854
855      subroutine docond(coo,iconn,ianz)
856      implicit double precision (a-h,p-x),integer (i-n),logical (o)
857      parameter (mxcon=10)
858      common /athlp/ iatoms, mxnat
859      dimension coo(3,*),ianz(*),iconn(mxcon+1,*)
860
861      call convar(coo,iconn,ianz,iatoms,0)
862
863      return
864      end
865
866      subroutine convar(coo,iconn,ianz,iatoms,ioff)
867      implicit double precision (a-h,p-x),integer (i-n),logical (o)
868      parameter (mxel=100)
869      parameter (mxcon=10)
870      common /elmcom/ vdwr(mxel),vrad(mxel),icol(mxel)
871      common /rdwr/  iun1,iun2,iun3,iun4,iun5
872      dimension coo(3,*),iconn(mxcon+1,*),ianz(*),ctemp(3)
873
874      adjus=0.52917706d0
875
876      do i=ioff+1,ioff+iatoms
877
878         k = 0
879         iconn(1,i) = 0
880         iani = ianz(i)
881
882         do j=ioff+1,ioff+iatoms
883
884            ianj = ianz(j)
885
886            if (iani.gt.0.and.ianj.gt.0.and.
887     &          iani.le.100.and.ianj.le.100) then
888
889               dmaxsq = (vdwr(iani) + vdwr(ianj))/adjus
890               dmaxsq = dmaxsq*dmaxsq
891
892               dijsq = 0
893               do l=1,3
894                  ctemp(l) = coo(l,i) - coo(l,j)
895                  dijsq = dijsq + ctemp(l)*ctemp(l)
896               end do
897
898               if (i.ne.j.and.iani.ne.99.and.ianj.ne.99) then
899                 if (dijsq.lt.dmaxsq) then
900                     k = k + 1
901                     if (k.le.mxcon) then
902                        iconn(1,i) = k
903                        iconn(k+1,i) = j
904                     else
905                        write(iun3,*)
906     &                     'more than mxconn connections found'
907                     endif
908                 endif
909               endif
910
911            endif
912
913         end do
914      end do
915
916      return
917      end
918
919      subroutine dohcod(ihflag,coo,ianz,iconn,isurf)
920      implicit double precision (a-h,p-x),integer (i-n),logical (o)
921      parameter (mxcon=10)
922      parameter (mxel=100)
923      parameter (nacc=6)
924      parameter (ndon=4)
925      common /athlp/ iatoms, mxnat
926      logical hon
927      common /hcon/ iacc(nacc),idon(ndon),hdmin,hdmax,hamin,hamax,hon
928      common /rdwr/  iun1,iun2,iun3,iun4,iun5
929      logical isacc,isdon
930      real angle,hamin,hamax
931      dimension ctemp(3),isel(4)
932      dimension coo(3,*),ianz(*),iconn(mxcon+1,*),isurf(*)
933
934      call domcon(1,1)
935
936      if (.not.hon) return
937
938      hdmin2 = hdmin*hdmin
939      hdmax2 = hdmax*hdmax
940
941c
942c Hydrogen bonds
943c
944      do i=1,iatoms
945         do j=i+1,iatoms
946           if (ihflag.eq.0.or.
947     &     (ihflag.gt.0.and.(isurf(i).eq.1.or.isurf(j).eq.1))) then
948            iai = ianz(i)
949            iaj = ianz(j)
950            dijsq = 0
951            do l=1,3
952               ctemp(l) = coo(l,i) - coo(l,j)
953               dijsq = dijsq + ctemp(l)*ctemp(l)
954            end do
955            if ((iai.eq.1.and.isacc(iaj)).or.
956     &          (iaj.eq.1.and.isacc(iai))) then
957               if (dijsq.gt.hdmin2.and.dijsq.lt.hdmax2) then
958                  isel(1) = 0
959                  if (iai.eq.1) then
960                      if (iconn(1,i).gt.0) then
961                         iai = ianz(iconn(2,i))
962                         if (isdon(iai)) then
963                            isel(1) = iconn(2,i)
964                            isel(2) = i
965                            isel(3) = j
966                         endif
967                      endif
968                  else
969                      if (iconn(1,j).gt.0) then
970                         iaj = ianz(iconn(2,j))
971                         if (isdon(iaj)) then
972                            isel(1) = iconn(2,j)
973                            isel(2) = j
974                            isel(3) = i
975                         endif
976                      endif
977                  endif
978                  if (isel(1).ne.0) then
979                    call intcor(intc,angle,isel,3)
980                    if (intc.eq.1) then
981                     if (abs(angle).gt.hamin.and.
982     &                    abs(angle).lt.hamax) then
983                         k = iconn(1,i)
984                         if (k.lt.mxcon) then
985                            iconn(1,i) = iconn(1,i) + 1
986                            iconn(iconn(1,i)+1,i) = -j
987                         else
988                            write(iun3,*)
989     &                        'more than mxconn connections found'
990                         endif
991                         k = iconn(1,j)
992                         if (k.lt.mxcon) then
993                            iconn(1,j) = iconn(1,j) + 1
994                            iconn(iconn(1,j)+1,j) = -i
995                         else
996                            write(iun3,*)
997     &                        'more than mxconn connections found'
998                         endif
999                     endif
1000                    endif
1001                  endif
1002               endif
1003            endif
1004           endif
1005         end do
1006      end do
1007
1008      if (ihflag.eq.1) call clrsrf
1009
1010      return
1011      end
1012
1013      subroutine nohcod(iconn)
1014      implicit double precision (a-h,p-x),integer (i-n),logical (o)
1015      parameter (mxcon=10)
1016      common /athlp/ iatoms, mxnat
1017      dimension icn(mxcon+1),iconn(mxcon+1,*)
1018
1019c
1020c Deactivate Hydrogen bonds
1021c
1022      do i=1,iatoms
1023          inn = 0
1024          do j=1,iconn(1,i)
1025             if (iconn(1+j,i).gt.0) then
1026                inn = inn + 1
1027                icn(inn) = iconn(1+j,i)
1028             endif
1029          end do
1030          iconn(1,i) = inn
1031          do j=1,inn
1032             iconn(1+j,i) = icn(j)
1033          end do
1034      end do
1035
1036      return
1037      end
1038
1039      logical function isacc(iatnr)
1040      implicit double precision (a-h,p-x),integer (i-n),logical (o)
1041      parameter (nacc=6)
1042      parameter (ndon=4)
1043      logical hon
1044      common /hcon/ iacc(nacc),idon(ndon),hdmin,hdmax,hamin,hamax,hon
1045      real hamin,hamax
1046
1047      isacc = .false.
1048
1049      do i=1,nacc
1050          if (iatnr.eq.iacc(i)) isacc = .true.
1051      end do
1052
1053      return
1054      end
1055
1056      logical function isdon(iatnr)
1057      implicit double precision (a-h,p-x),integer (i-n),logical (o)
1058      parameter (nacc=6)
1059      parameter (ndon=4)
1060      logical hon
1061      common /hcon/ iacc(nacc),idon(ndon),hdmin,hdmax,hamin,hamax,hon
1062      real hamin,hamax
1063
1064      isdon = .false.
1065
1066      do i=1,ndon
1067          if (iatnr.eq.idon(i)) isdon = .true.
1068      end do
1069
1070      return
1071      end
1072
1073      subroutine doscad(t,coo,zv,pincr,scal,scali,smag,iscupd)
1074      implicit double precision (a-h,p-z),integer (i-n),logical (o)
1075      common /athlp/ iatoms, mxnat
1076      common /gracom/  uscl,colscd,colscpd,ivdwpl
1077      dimension t(3),coo(3,*)
1078
1079      scali = 0.0d0
1080      do i=1,iatoms
1081         dist = 0.0d0
1082         do j=1,3
1083             dt = coo(j,i)-t(j)
1084             dist = dist + dt*dt
1085         end do
1086         if (dist.gt.scali) scali = dist
1087      end do
1088      if (scali.eq.0.0d0) then
1089         scali = 1.0d0
1090      else
1091         scali = dsqrt(scali)
1092      endif
1093      scal = scali * 2.4d0 * smag
1094
1095      if (scali.gt.15.0d0) then
1096         colscpd = 1.3d0
1097      else
1098         colscpd = .4d0
1099      endif
1100
1101      colscd = 2.0d0
1102
1103      if (iscupd.eq.1) then
1104         zv = scali
1105         iscupd = 0
1106      endif
1107
1108      pincr = 0.02d0*scali
1109
1110      return
1111      end
1112
1113      subroutine redcod(iconn)
1114      implicit double precision (a-h,p-x),integer (i-n),logical (o)
1115      parameter (mxcon=10)
1116      common /athlp/ iatoms, mxnat
1117      dimension icnn(mxcon),iconn(mxcon+1,*)
1118
1119      do j=1,iatoms
1120         ibnds = 0
1121         do i=1,iconn(1,j)
1122            ioke = 1
1123            do k=1,i-1
1124               if (iconn(k+1,j).eq.iconn(i+1,j)) ioke = 0
1125            end do
1126            if (ioke.eq.1.and.ibnds.lt.mxcon) then
1127               ibnds = ibnds + 1
1128               icnn(ibnds) = iconn(i+1,j)
1129            endif
1130         end do
1131         iconn(1,j) = ibnds
1132         do i=1,ibnds
1133            iconn(i+1,j) = icnn(i)
1134         end do
1135      end do
1136
1137      return
1138      end
1139
1140      subroutine domcod(ipt,idoall,monmod,iconn)
1141      implicit double precision (a-h,p-x),integer (i-n),logical (o)
1142      parameter (mxcon=10)
1143      parameter (maxdm=20)
1144      common /distmn/ rdm(maxdm),idmon(2,maxdm),ndm
1145      integer srfmap, srfloc
1146      common /hlpsrf/ npts1,npts2,npts3,srfmap,srfloc,ifogl,itsrf
1147      real rout
1148      dimension iconn(mxcon+1,*)
1149
1150
1151      if (idoall.eq.1) then
1152         ibeg = 1
1153         iend = ndm
1154      else
1155         ibeg = ipt
1156         iend = ipt
1157      endif
1158
1159      do i=ibeg,iend
1160         if (monmod.eq.2) then
1161            call intcor(intc,rout,idmon(1,i),2)
1162         else
1163            intc = 1
1164         endif
1165         if (intc.eq.1) then
1166            if (monmod.eq.2) rdm(i) = dble(rout)*0.52917706d0
1167            if (ifogl.ne.1) then
1168               i1 = idmon(1,i)
1169               i2 = idmon(2,i)
1170               k = iconn(1,i1)
1171               if (k.lt.mxcon) then
1172                  iconn(1,i1) = k + 1
1173                  iconn(k+2,i1) = -i2
1174               endif
1175               k = iconn(1,i2)
1176               if (k.lt.mxcon) then
1177                  iconn(1,i2) = k + 1
1178                  iconn(k+2,i2) = -i1
1179               endif
1180            endif
1181         endif
1182      end do
1183
1184      call ogmon
1185
1186      return
1187      end
1188
1189      subroutine chkbcd(ncalf,ihet,reson,iaton,iatclr,ision)
1190      implicit double precision (a-h,p-x),integer (i-n),logical (o)
1191      integer fancy,shade,atcol,dolabs,persp,fyesno,backb
1192      common /displ/ fancy,shade,atcol,dolabs,persp,irtcel,
1193     &               ifd,fyesno,backb,logo
1194      common /cllab/  iclon,iclpnt(4)
1195      common /boxion/ ionon,ionclr
1196      parameter (mxheta=150)
1197      integer reson
1198      dimension ihet(*),reson(*),iaton(*),iatclr(*)
1199
1200      if (backb.eq.1) then
1201         call actcal(1)
1202         call ribbs
1203         if (ihet(1).eq.1) call acthel(1,0,0,0)
1204         if (ihet(2).eq.1) call acthel(1,1,0,0)
1205         if (ihet(3).eq.1) call acthel(1,2,0,0)
1206         if (ihet(4).eq.1) call acthel(1,3,0,0)
1207         do i=5,mxheta
1208            j = -i+1
1209            if (ihet(i).eq.1) then
1210               call actami(j,0,1,0)
1211            endif
1212         end do
1213         if (ionon.eq.0) then
1214            call actami(ision,ionclr,1,0)
1215         endif
1216         do i=1,ncalf
1217            if (reson(i).eq.1) then
1218               call actami(i,0,1,0)
1219            endif
1220         end do
1221         if (iclon.eq.1) then
1222            icllow = iclpnt(1)
1223            do i=0,7
1224               iaton(icllow+i) = 1
1225               iatclr(icllow+i) = 11
1226            end do
1227         endif
1228      endif
1229
1230      return
1231      end
1232