c$DEBUG
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C    VECTOR - Visualisation of normal coordinate displacement vectors
C             calculated with VIBCA (option IFVCT=1)
C
C
C             The eigenvectors can be toggled through, scaled, rotated,
C             reversed, and plotted as weighted/unweighted.
C
C             Output can be written to file VECTOR.COR which gives a route
C             to plotting the previewed mode by means of some other program.
C             VECTOR.COR is also viewable directly with PMIFST.
C             The moments of inertia/rotational constants displayed while
C             generating this file may be of some diagnostic use, but the user
C             is warned not to put too much trust in them.
C
C
C       ver. 20.III.01                            ----- Zbigniew KISIEL -----
C                          __________________________________________________
C                         | Institute of Physics, Polish Academy of Sciences |
C                         | Al.Lotnikow 32/46, Warszawa, POLAND              |
C                         |                             kisiel@ifpan.edu.pl  |
C                         |     http://info.ifpan.edu.pl/~kisiel/prospe.htm  |
C_________________________/--------------------------------------------------
C
C  Modification history:
C
c     2.93:  Creation
c  2.09.98:  Addition of output to file from display mode
C 13.12.98:  Modification of vdW radius selection in BONLOC
c 27.03.99:  Addition of some colour
c 20.03.01:  Transfer of various mods from PMIFST
C
C----------------------------------------------------------------------------
C
C   NOTES:  This is an IBM-PC bound program based on the graphics package
C           first incorporated with MSF 5.0.  See the routine GRAPHICSMODE
c           for more details.
C
C           Dimensioning is based on the limit MAXAT for the number of
c           atoms in the molecule, change this in all PARAMETER
c           statements if necessary.
C_____________________________________________________________________________
C
C   When molecule is displayed under the rotation option, various actions
C   can be selected by the following keys:
C
C     U - toggle mass weighted/unweighted eigenvectors
C  up/down arrows - toggle through eigenvectors
C  left/right arrows - vary vector scaling
C     - - reverse vector phase
C   E,X - rotation about the horizontal axis
C   S,D - rotation about the vertical axis
C   W,R - rotation about the axis perpendicular to the plane of the screen
C         (the six keys above are case sensitive - lower case gives slow and
C          upper case fast rotation)
C     N - display numbers of atoms
C     A - display principal axes
C     O - generate output file VECTOR.COR with effects of motion along the
C         current mode on rotational parameters
C     P - parameter checks: display is switched to a mode allowing checks of
C         values of bond lengths and bond angles
C         In this mode:
C         - B selects calculation of bond length to be followed by two atom
C           numbers
C         - A selects calculation of bond angle to be followed by three atom
C           numbers, a,b or c can be input instead of the third number for
C           calculation of angle between bond and a principal axis
C         - D selects calculation of dihedral angle to be followed by
C           four atom numbers
C         - E exits back to previous screen
C  <cr> - termination of graphics display
C_____________________________________________________________________________
C
C
C...Initialization commands for graphics (FGRAPH.FI should only be included
C   into the .MAIN block, FGRAPH.FD should be declared in every program
C   block calling graphics routines).  The three structured
C   variables contain coordinates:
C     curpos.row and curpos.col - cursor coordinates (INTEGER*2)
C     ixy.xcoord and ixy.ycoord - pixel coordinates (INTEGER*2)
C     wxy.wx and wxy.wy - window coordinates (REAL*8)
C
      INCLUDE 'FGRAPH.FI'                                               MSF5+PS1
c     INCLUDE 'FLIB.FI'                                                 PS1
      INCLUDE 'FGRAPH.FD'                                               MSF5+PS1
c     USE MSFLIB                                                        PS4
      RECORD /rccoord/curpos
      RECORD /xycoord/ixy
      RECORD /wxycoord/wxy
C
C...The code in the two INCLUDE files clashes with IMPLICIT statements - so
C   when using MSF5.0 graphics the recomendation is to type everything
C   explicitly (check with -4Yg compilation switch)
C
      PARAMETER (maxat=25, maxvec=3*maxat, maxpts=maxat+maxvec,
     *           maxbon=6*maxpts)
c
      COMMON /BLK0/INITIA,NATOMS,ATOMS,VBA,TRANS,PRC,RJA,VCA
     1 /BLK1/NATMS,CORD,CORPLT,AXEND,axendr,ASPECT,ROTMAT
     1 /BLK2/IBOND,NBONDS,W
      COMMON /PARBLK/RANGEA,RANGEB,AMIN,BMIN
      COMMON /limits/maxx,maxy,LINOFS,curpos,ixy,wxy,
     * mymode,myrows,mycols,blue,red
      COMMON /lines/toplin,botlin,emplin,COMENT
      COMMON /EVECT/EW(maxvec,maxat,3),EU(maxvec,maxat,3),
     * FMODE(maxvec),NVEC
      REAL*8 CORD(maxpts,3),W(maxpts),VBA(3),VCA(3),TRANS(3,3),RJA(3),
     * PRC(3),CORDX(maxpts),CORDY(maxpts),CORDZ(maxpts),ROTMAT(3,3),
     * CORPLT(maxpts,3),ATOMS(maxpts,3),AXEND(6,3),axendr(6,3)
      REAL*8 CON,RANGEA,RANGEB,RANGEC,
     1 SMALLA,SMALLB,SMALLC,ASPECT,
     2 AMIN,BMIN,EW,EU,FMODE
      INTEGER*2 dummy,maxx,maxy,mymode,myrows,mycols,LINOFS
      INTEGER*4 dummy4,blue,red
      INTEGER*2 IBOND(maxbon,2),NATOMS(maxpts,4),INKEY,N
      CHARACTER*79 toplin,botlin,emplin,COMENT,LINE
      CHARACTER*25 FILNAM
C
      WRITE(emplin,'(79(1H ))')
      CON=0.174532925199D-1
      ITSAVE=1
      linofs=14
C
C   Input of data
C
      CALL clearscreen($GCLEARSCREEN)
      WRITE(*, 3000)
3000  FORMAT(1X/'  ',76(1H_)/' |',T79,'|'/
     * ' |     V E C T O R  -  Visualisation of normal',
     * ' coordinate displacement',T79,'|'/
     * ' |                     vectors from VIBCA',T79,'|'/
     * ' |',76(1H_),'|'/'  version 20.III.2001',T64,'Zbigniew KISIEL'/)
3003  FORMAT(I5)
C
6001  WRITE(*,3040)
3040  FORMAT(1X/' Name of VIBCA results file:  ',$)
      READ(*,'(A)',ERR=6001)FILNAM
      OPEN(2,FILE=FILNAM,STATUS='OLD',ERR=6001)
C
C   INPUT OF VIBCA RESULTS FILE
C
C...coordinates (ignore dummy atoms)
C
6600  read(2,'(a)',end=6601)line
      if(line(1:7).ne.' ======')goto 6600
      read(2,'(a)')line
      coment=line
6500  READ(2,'(A)')LINE
      IF(LINE(2:22).EQ.'PRINCIPAL COORDINATES')THEN
        READ(2,'(A)')LINE
        READ(2,'(A)')LINE
        READ(2,'(A)')LINE
        GOTO 6020
      ENDIF
      GOTO 6500
C
6020  N=0
6501  READ(2,'(A)')LINE
      IF(LINE(4:4).EQ.' ')THEN
        N=N+1
        READ(LINE,6506)NNNN,CORDX(N),CORDY(N),CORDZ(N),W(N)
        IF(W(N).EQ.0.D0)N=N-1
        GOTO 6501
      ENDIF
6506  FORMAT(1X,I2,7X,3(F12.6,2X),F12.6)
      NATMS=N
      write(*,'(1x/i6,a)')natms,' atoms'
      DO 6511 N=1,NATMS
        WRITE(*,'(10X,4F10.4)')CORDX(N),CORDY(N),CORDZ(N),W(N)
6511  CONTINUE
C
C...eigenvectors, skip the rotation+translation modes marked by VIBCA
C   as 'Mode 0'
C
6502  READ(2,'(A)')LINE
      IF(LINE(2:9).NE.'Mode  1:')GOTO 6502
C
      BACKSPACE(2)
      NVEC=0
C
6503  READ(2,'(A)')LINE
      IF(LINE(2:5).EQ.'Mode')THEN
        NVEC=NVEC+1
        READ(2,*)FMODE(NVEC),J,(EW(NVEC,1,I),I=1,3),(EU(NVEC,1,I),I=1,3)
        DO 6504 N=2,NATMS
          READ(2,*)J,(EW(NVEC,N,I),I=1,3),(EU(NVEC,N,I),I=1,3)
6504    CONTINUE
        GOTO 6503
      ENDIF
      write(*,'(1x,i5,a)')nvec,' eigenvectors'
      CLOSE(2)
C
C
C...Location of bonds, determination of scaling parameters and of points
C   defining principal axes
C
      DO 6505 N=1,NATMS
        CORD(N,1)=CORDX(N)
        CORD(N,2)=CORDY(N)
        CORD(N,3)=CORDZ(N)
        CORD(N+NATMS,1)=CORD(N,1)+EW(1,N,1)
        CORD(N+NATMS,2)=CORD(N,2)+EW(1,N,2)
        CORD(N+NATMS,3)=CORD(N,3)+EW(1,N,3)
        W(N+NATMS)=0.000001D0
6505  CONTINUE
C
      CALL BONLOC
      write(*,'(1x,i5,a)')nbonds,' bonds'
      NATMS=NATMS+NATMS
C
      DO 3500 N=1,NATMS
        CORPLT(N,1)=CORD(N,1)
        CORPLT(N,2)=CORD(N,2)
        CORPLT(N,3)=CORD(N,3)
3500  CONTINUE
C
      DO 510 N=1,6
        DO 510 NN=1,3
          AXEND(N,NN)=0.D0
510   CONTINUE
      CALL SCALIN(1,RANGEA,SMALLA)
      AXEND(1,1)=SMALLA
      AXEND(2,1)=SMALLA+RANGEA
      CALL SCALIN(2,RANGEB,SMALLB)
      AXEND(3,2)=SMALLB
      AXEND(4,2)=SMALLB+RANGEB
      CALL SCALIN(3,RANGEC,SMALLC)
      AXEND(5,3)=SMALLC
      AXEND(6,3)=SMALLC+RANGEC
      IF(RANGEC.LT.2.D0)THEN
        AXEND(5,3)=-1.D0
        AXEND(6,3)=1.D0
      ENDIF
      do 511 n=1,6
        do 511 nn=1,3
          axendr(n,nn)=axend(n,nn)
511   continue
C
C
C...Plots of projections on the three principal planes
C
6014  WRITE(*, 4001)
4001  FORMAT(1X/' Projections on the three principal planes ?  ',$)
      READ(*,3003,ERR=6014)IFLAG
      IF(IFLAG.NE.1)GO TO 4023
C
C...initialise graphics - see comments in routine ROTATE for details
      call graphicsmode()
      ASPECT=maxy/(0.75*maxx)
      dummy=setactivepage(1)
      dummy=setvisualpage(1)
      pixcol=getcolor()
      dummy4 = setbkcolor( BLUE )
$NODEBUG
      CALL setlinestyle(#FFFF)
$DEBUG
      CALL moveto(0,linofs,ixy)
      dummy=lineto(maxx,linofs)
      CALL moveto(0,maxy-linofs,ixy)
      dummy=lineto(maxx,maxy-linofs)
      WRITE(toplin,'(a)')coment
      WRITE(botlin,'(2A)')'Projections on the three principal planes',
     1 '                  Press ENTER to exit'
      CALL settextposition(1,1,curpos)
      CALL outtext(toplin)
      CALL settextposition(myrows,1,curpos)
      CALL outtext(botlin)
C
      IH=maxx/2
      IV=maxy/2
C
C...A-B plane
      CALL SCALIN(1,RANGEA,SMALLA)
      CALL SCALIN(2,RANGEB,SMALLB)
      AMEAN=SMALLA+0.5D0*RANGEA
      BMEAN=SMALLB+0.5D0*RANGEB
      RANGEA=1.05D0*RANGEA
      RANGEB=1.05D0*RANGEB
      PIXELR=REAL(IV)/REAL(IH)
      IF(RANGEB*ASPECT/RANGEA.LE.PIXELR)THEN
        RANGEB=RANGEA*PIXELR/ASPECT
      ELSE
        RANGEA=RANGEB*ASPECT/PIXELR
      ENDIF
      AMIN=AMEAN-RANGEA*0.5D0
      AMAX=AMEAN+RANGEA*0.5D0
      BMIN=BMEAN-0.5D0*RANGEB
      BMAX=BMEAN+0.5D0*RANGEB
      CALL setviewport(0,linofs+1,maxx-IH,linofs+1+IV)
      dummy=setwindow(.TRUE.,amin,bmin,amax,bmax)
      CALL view(1,2,1,0)
C
C...B-C plane
      CALL SCALIN(2,RANGEA,SMALLA)
      CALL SCALIN(3,RANGEB,SMALLB)
      AMEAN=SMALLA+0.5D0*RANGEA
      BMEAN=SMALLB+0.5D0*RANGEB
      RANGEA=1.05D0*RANGEA
      RANGEB=1.05D0*RANGEB
      IF(RANGEB*ASPECT/RANGEA.LE.PIXELR)THEN
        RANGEB=RANGEA*PIXELR/ASPECT
      ELSE
        RANGEA=RANGEB*ASPECT/PIXELR
      ENDIF
      AMIN=AMEAN-RANGEA*0.5D0
      AMAX=AMEAN+RANGEA*0.5D0
      BMIN=BMEAN-0.5D0*RANGEB
      BMAX=BMEAN+0.5D0*RANGEB
      CALL setviewport(maxx-IH,linofs+1,maxx,linofs+1+IV)
      dummy=setwindow(.TRUE.,amin,bmin,amax,bmax)
      CALL view(2,3,1,0)
C
C...A-C plane
      CALL SCALIN(1,RANGEA,SMALLA)
      CALL SCALIN(3,RANGEB,SMALLB)
      AMEAN=SMALLA+0.5D0*RANGEA
      BMEAN=SMALLB+0.5D0*RANGEB
      RANGEA=1.05D0*RANGEA
      RANGEB=1.05D0*RANGEB
      IF(RANGEB*ASPECT/RANGEA.LE.PIXELR)THEN
        RANGEB=RANGEA*PIXELR/ASPECT
      ELSE
        RANGEA=RANGEB*ASPECT/PIXELR
      ENDIF
      AMIN=AMEAN-RANGEA*0.5D0
      AMAX=AMEAN+RANGEA*0.5D0
      BMIN=BMEAN-0.5D0*RANGEB
      BMAX=BMEAN+0.5D0*RANGEB
      IHH=maxx/2-IH/2
      CALL setviewport(IHH,maxy-linofs-1-IV,IHH+IH,maxy-linofs-1)
      dummy=setwindow(.TRUE.,amin,bmin,amax,bmax)
      CALL view(1,3,1,0)
C
500   IK=INKEY(N)
      IF(IK.NE.13)GOTO 500
      dummy=setvideomode($DEFAULTMODE)
C
C...Rotations of the molecule and structure checks
C
4023  WRITE(*,'(/'' Rotations and structure checks ?    '',$)')
      READ(*,3003,ERR=4023)IFLAG
      IF(IFLAG.NE.1)STOP
      CALL ROTATE
C
6601  STOP
      END
C
C_____________________________________________________________________________
C
      SUBROUTINE BONLOC
C
C   BONLOC locates bonds in the molecule ie. atoms that are closer together
C   than the sum of their covalent radii + 0.1 A
C
C      - numbers of atoms involved in an identified bond are placed in IBOND,
C        NBONDS is the number of bonds identified
C      - BONLOC will only recognise atoms of H,C,N,P,O,S,F,CL,BR,I in their
C        stable isotopic species

      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (NRADII=46)
      PARAMETER (maxat=25, maxvec=3*maxat, maxpts=maxat+maxvec,
     *           maxbon=6*maxpts)
c
      COMMON /BLK1/NATMS,CORD,CORPLT,AXEND,axendr,ASPECT,ROTMAT(3,3)
     *       /BLK2/IBOND,NBONDS,WT
      COMMON /EVECT/EW(maxvec,maxat,3),EU(maxvec,maxat,3),
     *              FMODE(maxvec),NVEC
      REAL*8 CORD(maxpts,3),WT(maxpts),CORPLT(maxpts,3),
     * AXEND(6,3),axendr(6,3),ASPECT
      INTEGER*2 IBOND(maxbon,2),MASS(NRADII)
      REAL*4 CORAD(NRADII)
C
      DATA MASS/  1,    2,  10,  11,  12,  13,  14,  15,  16,  17,
     *  18,  19,  28,  29,  30,  31,  32,  33,  34,  35,  37,
     *  70,  72,  73,  74,  76,  75,  77,  78,  80,  82,  79,  81,
     * 116, 117, 118, 119, 120, 121, 123, 125, 126, 128, 130, 127,
     *  40/
      DATA CORAD/0.37,0.37,0.81,0.81,0.77,0.77,0.74,0.74,0.71,0.71,
     * 0.71,0.72,1.17,1.17,1.17,1.10,1.04,1.04,1.04,0.99,0.99,
     * 1.22,1.22,1.22,1.22,1.22,1.21,1.17,1.17,1.17,1.17,1.14,1.14,
     * 1.40,1.40,1.40,1.40,1.40,1.41,1.41,1.37,1.37,1.37,1.37,1.33,
     * 0.1/
C
10    FORMAT(1X,'MASS OF ATOM.',I2,' = ',F8.4,'  IS NOT INCLUDED',
     1 ' IN BONLOC RECOGNITION TABLE')
C
      NBONDS=0
      DO 9 NSTART=2,NATMS
        NNN=NSTART-1
        MMASS=IDINT(WT(NNN)/0.5D0)*0.5+0.5
C
C...treatment of dummy atoms
        IF(MMASS.EQ.0)THEN
          NBONDS=NBONDS+1
          IBOND(NBONDS,1)=NNN
          IBOND(NBONDS,2)=NNN
          GOTO 9
        ENDIF
C
C...treatment of normal atom pairs, first identify the current atom, then
C   go through remaining atoms and determine whether any are closer than sum
C   of covalent radii
        DO 3 NN=1,NRADII
          IF(MMASS.EQ.MASS(NN))GO TO 4
3       CONTINUE
        WRITE(*, 10)NNN,WT(NNN)
4       RAD=CORAD(NN)
        DO 1 N=NSTART,NATMS
          MMASS=IDINT(WT(N)/0.5)*0.5+0.5
          IF(MMASS.EQ.0)GOTO 1
          DO 2 NN=1,NRADII
            IF(MMASS.EQ.MASS(NN))GO TO 5
2         CONTINUE
          WRITE(*, 10)N,WT(N)
5         DIST=RAD+CORAD(NN)
          SEP=SQRT((CORD(NNN,1)-CORD(N,1))**2+(CORD(NNN,2)-CORD(N,2))
     1     **2 +(CORD(NNN,3)-CORD(N,3))**2)
          IF(SEP.LT.DIST+0.1)GO TO 6
          GO TO 1
6         NBONDS=NBONDS+1
          IBOND(NBONDS,1)=NNN
          IBOND(NBONDS,2)=N
1       CONTINUE
9     CONTINUE
C
C...bonds representing normal coordinate vectors
C
      DO 11 NN=1,NATMS
        NBONDS=NBONDS+1
        IBOND(NBONDS,1)=NN
        IBOND(NBONDS,2)=NN+NATMS
11    CONTINUE
C
      RETURN
      END
C
C_____________________________________________________________________________
C
      SUBROUTINE SCALIN(NAXIS,RANGEX,SMALLX)
C
C   Determination of minimum cooordinate value and the range of coordinate
C   values along a given principal axis
C
C
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (maxat=25, maxvec=3*maxat, maxpts=maxat+maxvec,
     *           maxbon=6*maxpts)
      COMMON /BLK1/NATMS,CORD,CORPLT,AXEND,axendr,ASPECT,ROTMAT(3,3)
      REAL*8 CORD(maxpts,3),CORPLT(maxpts,3),AXEND(6,3),axendr(6,3),
     * ASPECT
C
      BIGX=-10.0
      SMALLX=10.0
C
      DO 1 N=1,NATMS
        IF(CORPLT(N,NAXIS).GT.BIGX)BIGX=CORPLT(N,NAXIS)
        IF(CORPLT(N,NAXIS).LT.SMALLX)SMALLX=CORPLT(N,NAXIS)
1     CONTINUE
C
      RANGEX=1.1d0*(BIGX-SMALLX)
      smallx=1.05d0*smallx
C
      RETURN
      END
C_____________________________________________________________________________
C
      SUBROUTINE graphicsmode()
C
C   This routine sets the graphics mode on the PC and
C   determines the pixel limits on x and y axes (0,maxx), (0,maxy)
C
      INCLUDE  'FGRAPH.FD'                                              MSF5+PS1
c     USE MSFLIB                                                        PS4
C
      RECORD /rccoord/curpos
      RECORD /xycoord/ixy
      RECORD /wxycoord/wxy
      INTEGER*2 dummy,maxx,maxy,linofs,mymode,myrows,mycols
      INTEGER*4 blue,red
      RECORD /videoconfig/myscreen
      COMMON /limits/maxx,maxy,linofs,curpos,ixy,wxy,
     * mymode,myrows,mycols,blue,red
C
C...Set graphics mode.
C
      CALL getvideoconfig(myscreen)
C
C...this setting was an old default (EGA and 25 rows of text)
c         mymode=$ERESCOLOR
c         myrows=25
C
C...this setting should have widest compatibility (640*480, 30 line VGA)
c         mymode=$VRES256COLOR
          mymode=$VRES16COLOR
          myrows=30
C
C...this is for 800*600 SVGA, usually 16 colors is better than 256
c         mymode=$SRES256COLOR
c         mymode=$SRES16COLOR
c         myrows=37
c
          dummy = setvideomoderows ( mymode, myrows )
          myrows=dummy
          if(myrows.ne.25.and.myrows.ne.30.and.myrows.ne.37)then
          write(*,'(1x//i5,a,a//)')myrows,
     *     ' <--- this is the number of lines in this mode',
     *           ' and it is not yet supported'
        pause
      endif
C
      IF( dummy .EQ. 0 ) STOP 'Error:  cannot set graphics mode'
C
C...Determine parameters of the mode
C
      CALL getvideoconfig( myscreen )
      maxx = myscreen.numxpixels - 1
      maxy = myscreen.numypixels - 1
      mycols=myscreen.numtextcols
      linofs=nint(real(maxy+1)/real(myrows))
c
c...patch for different use of SETBKCOLOR by PS4 from PS1 (and MSF5)
c
      blue=$BLUE
      red=$RED
      if($BLUE.ne.#2a0000)then
        blue=1
        red=4
      endif
c
      dummy4 = setbkcolor( BLUE )
      call clearscreen($GCLEARSCREEN)
C
      RETURN
      END
C_____________________________________________________________________________
C
      SUBROUTINE VIEW(IX,IY,ICONTR,ICLEAR)
C
C   Graphics plot of rectangular projection of molecular structure onto
C   one of the principal planes
C
C   NOTE: the pixel origin is in the top lh. corner with Y increasing
C      downwards, but in the floating point window declared with the
C      option .TRUE. the Y-axis increases in the usual up direction.
C
C         IX, IY - principal axes to be used for X and Y coordinates of
C                  the plot (=1,2,3 for a,b,c resp.)
C         ICONTR = 0 draw the molecule
C                  +1 draw the molecule and principal axes (do not label)
C                  +2 draw the molecule and principal axes (label axes)
C                  +10 number atoms
C         ICLEAR = 0 do not clear display prior to drawing
C                  1 clear display
C
      INCLUDE  'FGRAPH.FD'                                              MSF5+PS1
c     USE MSFLIB                                                        PS4
      RECORD /rccoord/curpos
      RECORD /xycoord/ixy
      RECORD /wxycoord/wxy
C
      PARAMETER (maxat=25, maxvec=3*maxat, maxpts=maxat+maxvec,
     *           maxbon=6*maxpts)
c
      COMMON /TXTBUF/LINE
      COMMON /limits/maxx,maxy,LINOFS,curpos,ixy,wxy,
     * mymode,myrows,mycols,blue,red
      COMMON /BLK1/NATMS,CORD,CORPLT,AXEND,axendr,ASPECT,ROTMAT
     *       /BLK2/IBOND,NBONDS,WT
      COMMON /PARBLK/RANGEX,RANGEY,XMIN,YMIN
      REAL*8 CORD(maxpts,3),WT(maxpts),CORPLT(maxpts,3),
     * AXEND(6,3),axendr(6,3),XMIN,YMIN
      REAL*8 ASPECT,RANGEX,RANGEY,x1,y1,x2,y2,DELX,DELY,ROTMAT(3,3)
      INTEGER*2 dummy,maxx,maxy,mymode,myrows,mycols,LINOFS
      INTEGER*4 blue,red
      INTEGER*2 IBOND(maxbon,2)
      INTEGER*4 ISTART,JSTART
      CHARACTER LINE*80,AXNAM(3)
      DATA AXNAM/'A','B','C'/
      arlen=1.0
      arangl=30.
      TANTH=tan(arangl/57.29578D0)
C
      iaxis=icontr-(icontr/10)*10
      inumbr=icontr/10
      IYSPAN=maxy-2*linofs
      IF(ICLEAR.EQ.1)CALL clearscreen($GVIEWPORT)
C
C...Plot cirles at atom positions (small circles normally, large circles
C   for atom numbers)
C
      DELX=RANGEX/50.D0
      DELY=RANGEY/50.D0
      DUMMY=SETCOLOR(15)
      DO 3 N=1,NATMS
        IF(WT(N).LT.0.1D0)GOTO 3
        IF(inumbr.LT.1)THEN
          X1=CORPLT(N,IX)-DELX*0.3D0
          Y1=CORPLT(N,IY)+DELY*0.3D0
          X2=CORPLT(N,IX)+DELX*0.3D0
          Y2=CORPLT(N,IY)-DELY*0.3D0
        ELSE
          X1=CORPLT(N,IX)-DELX
          Y1=CORPLT(N,IY)+DELY
          X2=CORPLT(N,IX)+DELX
          Y2=CORPLT(N,IY)-DELY
        ENDIF
        dummy=ellipse_w($GFILLINTERIOR,X1,Y1,X2,Y2)
3     CONTINUE
C
C...Plot and label axes
C
      IF(IAXIS.LT.1)GOTO 110
      dummy=setcolor(11)
      DO 4 N=1,5,2
$NODEBUG
        CALL setlinestyle(#8888)
$DEBUG
        call moveto_w(AXEND(N,IX),AXEND(N,IY),wxy)
        dummy=lineto_w(AXEND(N+1,IX),AXEND(N+1,IY))
4     CONTINUE
$NODEBUG
          CALL setlinestyle(#FFFF)
$DEBUG
      dummy=setcolor(15)
      IF(IAXIS.EQ.2)THEN
        dummy=setcolor(11)
        DO 5 N=2,6,2
          ISTART=maxx*(AXEND(N,IX)-XMIN)/RANGEX-4
          JSTART=IYSPAN-IYSPAN*(AXEND(N,IY)-YMIN)/RANGEY+4
          WRITE(LINE,'(A)')AXNAM(N/2)
          CALL PLOTCH(2,2,ISTART,JSTART,-1)
          CALL PLOTCH(2,2,ISTART+1,JSTART,-1)
5       CONTINUE
        dummy=setcolor(15)
      ENDIF
C
C...Plot bonds, atom bonds with continuous line, bonds to vectors
C   with red line, dummy atoms not bonded
C
110   DO 8 N=1,NBONDS
        IF(IBOND(N,1).EQ.IBOND(N,2))GOTO 8
C
c...bond colour
        IF(WT(IBOND(N,2)).le.0.1D0)THEN
          DUMMY=SETCOLOR(12)
        ELSE
          DUMMY=SETCOLOR(15)
        ENDIF
c
c...plot bond
        CALL moveto_w(CORPLT(IBOND(N,1),IX),
     1                 CORPLT(IBOND(N,1),IY),wxy)
        dummy=lineto_w(CORPLT(IBOND(N,2),IX),CORPLT(IBOND(N,2),IY))
c
c...plot arrowhead
        IF(WT(IBOND(N,2)).le.0.1D0)THEN
          AA= corplt(ibond(n,2),ix)-corplt(ibond(n,1),ix)
          IF(aa.ne.0.d0)THEN
            RMB= (corplt(ibond(n,2),iy)-corplt(ibond(n,1),iy))/
     *           (corplt(ibond(n,2),ix)-corplt(ibond(n,1),ix))
          ELSE
            RMB=1.D+10
          ENDIF
          if(1.d0-RMB*TANTH.ne.0.d0)then
            rm=(TANTH+RMB)/(1.d0-RMB*TANTH)
          else
            rm=1.D+10
          endif
          if(1.d0+RMB*TANTH.ne.0.d0)then
            rm1=(-TANTH+RMB)/(1.d0+RMB*TANTH)
          else
            rm1=1.D+10
          endif
c
          DY=dsqrt(rm**2/(1.d0+rm**2))*arlen
          IF(RM.ne.0.0D0)then
            DX=DY/RM
          else
            DX=arlen
          endif
          DY=DY*dely
          DX=DX*delx
          X1=corplt(ibond(n,2),ix)+DX
          Y1=corplt(ibond(n,2),iy)+DY
          X2=corplt(ibond(n,2),ix)-DX
          Y2=corplt(ibond(n,2),iy)-DY
          AA=(X1-corplt(ibond(n,1),ix))**2+
     *       (Y1-corplt(ibond(n,1),iy))**2
          BB=(X2-corplt(ibond(n,1),ix))**2+
     *       (Y2-corplt(ibond(n,1),iy))**2
          if(BB.lt.AA)then
            X1=X2
            Y1=Y2
          endif
          CALL moveto_w(CORPLT(IBOND(N,2),IX),CORPLT(IBOND(N,2),IY),wxy)
          dummy=lineto_w(X1,Y1)
c
          DY=dsqrt(rm1**2/(1.d0+rm1**2))*arlen
          IF(RM1.ne.0.0D0)then
            DX=DY/RM1
          else
            DX=arlen
          endif
          DY=DY*dely
          DX=DX*delx
          X1=corplt(ibond(n,2),ix)+DX
          Y1=corplt(ibond(n,2),iy)+DY
          X2=corplt(ibond(n,2),ix)-DX
          Y2=corplt(ibond(n,2),iy)-DY
          AA=(X1-corplt(ibond(n,1),ix))**2+
     *       (Y1-corplt(ibond(n,1),iy))**2
          BB=(X2-corplt(ibond(n,1),ix))**2+
     *       (Y2-corplt(ibond(n,1),iy))**2
          if(BB.lt.AA)then
            X1=X2
            Y1=Y2
          endif
          CALL moveto_w(CORPLT(IBOND(N,2),IX),CORPLT(IBOND(N,2),IY),wxy)
          dummy=lineto_w(X1,Y1)
        endif
c
8     CONTINUE
      DUMMY=SETCOLOR(15)
C
C...Label atoms, first floating point coordinates of atoms have to be
C   converted to pixel coordinates, then the label is written using the
C   plotter character routine
C
$NODEBUG
          CALL setlinestyle(#FFFF)
$DEBUG
      IF(INUMBR.LT.1)GOTO 50
      DO 51 K=1,NATMS
        ISTART=maxx*(CORPLT(K,IX)-XMIN)/RANGEX
        JSTART=IYSPAN-IYSPAN*(CORPLT(K,IY)-YMIN)/RANGEY+2
        IF(K.LT.10)THEN
          WRITE(LINE,'(I1)')K
          ISTART=ISTART-1
        ELSE
          WRITE(LINE,'(I2)')K
          ISTART=ISTART-4
        ENDIF
        IF(WT(K).GT.0.5D0)THEN
          dummy=setcolor(0)
        ELSE
          dummy=setcolor(11)
        ENDIF
        CALL PLOTCH(1,1,ISTART,JSTART,1)
51    CONTINUE
      dummy=setcolor(15)
C
50    RETURN
      END
C
C_____________________________________________________________________________
C
      SUBROUTINE ROTATE
C
      INCLUDE  'FGRAPH.FD'                                              MSF5+PS1
c     USE MSFLIB                                                        PS4
      RECORD /rccoord/curpos
      RECORD /xycoord/ixy
      RECORD /wxycoord/wxy
c
      PARAMETER (maxat=25, maxvec=3*maxat, maxpts=maxat+maxvec,
     *           maxbon=6*maxpts)
C
      COMMON /limits/maxx,maxy,LINOFS,curpos,ixy,wxy,
     * mymode,myrows,mycols,blue,red
      COMMON /BLK1/NATMS,CORD,CORPLT,AXEND,axendr,ASPECT,ROTMAT
     *       /BLK2/IBOND,NBONDS,WT
     *     /PARBLK/RANGEA,RANGEB,AMIN,BMIN
      COMMON /lines/toplin,botlin,emplin,COMENT
      COMMON /EVECT/EW(maxvec,maxat,3),EU(maxvec,maxat,3),
     *              FMODE(maxvec),NVEC
      COMMON /MOMIA/cartdf,pmi
      REAL*8 CORD(maxpts,3),CORPLT(maxpts,3),WT(maxpts),ROTMAT(3,3),
     * AXEND(6,3),axendr(6,3),ASPECT,CARTDF(maxvec,4),PMI(3,3,2)
      REAL*8 ANGLE,RANGEA,RANGEB,SMALLA,SMALLB,AMIN,AMAX,BMIN,BMAX,
     1 AMEAN,BMEAN,EW,EU,FMODE
      INTEGER*2 dummy,maxx,maxy,mymode,myrows,mycols,LINOFS
      INTEGER*4 dummy4,blue,red
      INTEGER*2 IBOND(maxbon,2),INKEY,N
      CHARACTER*79 toplin,botlin,emplin,COMENT
      CHARACTER KK,OUTSTR*21,IWU(2)*10,nam(3)
C
      DATA IWU/'weighted  ','unweighted'/
      data nam/'a','b','c'/
      icontr=2
C
1     WRITE(*,'(/'' Specify the projection plane IA,IB'',t37,
     *        ''- use 1,2,3 for a,b,c:''//
     *        t31,''..... '',$)')
      READ(*,*,ERR=1)IA,IB
      IF(IA.EQ.IB)GOTO 1
      IF(IA.LT.1.OR.IA.GT.3)GOTO 1
      IF(IB.LT.1.OR.IB.GT.3)GOTO 1
      IC=6-(IA+IB)
      NVECTR=1
      IWTUNW=1
      NAT=NATMS/2
      SCALE=1.d0
      IPHASE=1
      call transf(0.d0,ia,ib,-1)
C
C---P-L-O-T---, first set up graphics conditions
C
C...Determine the graphics card and set the highest resolution graphics mode.
C   Use the second graphics page (page 1), this is convenient with Hercules
C   but may have to be changed to page 0 with some boards.  Also determine
C   pixel colour for use in erasing points from previous pass
C
333   call graphicsmode()
C
C...It is assumed that the ratio of physical screen dimensions Y/X is 0.75.
C   The pixel scaling ratio along the Y-coordinate for 1:1 X:Y dispersion
C   is maxy/(0.75*maxx) and is a unity for EGA and VGA, and 0.645 for
C   Hercules.
C
      ASPECT=maxy/(0.75*maxx)
C
c...Mode of operation
      dummy=setactivepage(1)
      dummy=setvisualpage(1)
      pixcol=getcolor()
      dummy4 = setbkcolor( BLUE )
C
C...The setlinestyle parameter is INTEGER*2 and for continuous line has to
C   be 65535 - this is identified as an error by $DEBUG since I*2 range is
C   +-32767, hence have to $NODEBUG locally
C
$NODEBUG
      CALL setlinestyle(#FFFF)
$DEBUG
C
C...Write top and bottom of display
C
      CALL moveto(0,linofs,ixy)
      dummy=lineto(maxx,linofs)
      dummy=lineto(maxx,maxy-linofs)
      dummy=lineto(0,maxy-linofs)
      dummy=lineto(0,linofs)
C
      WRITE(botlin,'(2A)')'Rotations: S-D E-X W-R Z       |O|utput ',
     1 '|A|xes |N|umbers |P|arameters   |H|elp'
C
33    CALL settextposition(1,1,curpos)
      WRITE(toplin,600)NVECTR,FMODE(NVECTR),IWU(IWTUNW),CHAR(24),
     * CHAR(25),CHAR(27),CHAR(26)
600   FORMAT('Mode',I2,':',F7.1,'cm-1, ',A10,13X,2A1,'=mode ',
     * 2A1,'=scale wtd|U|nwtd |-|phase')
      CALL outtext(toplin)
      CALL settextposition(myrows,1,curpos)
      CALL outtext(botlin)
C
C...Plot selected projection.  The first step is scaling and since this code
C   is executed on each rotation the scaling varies continuously for maximum
C   screen filling
C
3     CALL SCALIN(IA,RANGEA,SMALLA)
      CALL SCALIN(IB,RANGEB,SMALLB)
      AMEAN=SMALLA+0.5D0*RANGEA
      BMEAN=SMALLB+0.5D0*RANGEB
      RANGEA=1.05D0*RANGEA
      RANGEB=1.05D0*RANGEB
      PIXELR=REAL(maxy-2*linofs)/REAL(maxx)
      IF(RANGEB*ASPECT/RANGEA.LE.PIXELR)THEN
        RANGEB=RANGEA*PIXELR/ASPECT
      ELSE
        RANGEA=RANGEB*ASPECT/PIXELR
      ENDIF
      AMIN=AMEAN-RANGEA*0.5D0
      AMAX=AMEAN+RANGEA*0.5D0
      BMIN=BMEAN-0.5D0*RANGEB
      BMAX=BMEAN+0.5D0*RANGEB
C
C...Define graphics viewport and floating point window for automatic
C   coordinate to pixel conversion.  Note that setviewport
C   puts the origin of the pixel axes at (0,linofs+1).
C
      CALL setviewport(1,linofs+1,maxx-1,maxy-linofs-1)
      dummy=setwindow(.TRUE.,amin,bmin,amax,bmax)
C
C...Plot molecule
      CALL VIEW(IA,IB,ICONTR,1)
C
C---O-P-T-I-O-N-S---
C
2     IK=INKEY(N)
      IF(IK.EQ.13)GOTO 20
      IF(IK.EQ.-72.OR.IK.EQ.-80)GOTO 900
      IF(IK.EQ.-75.OR.IK.EQ.-77)GOTO 710
      IF(CHAR(IK).EQ.'-')GOTO 700
      IF(IK.LT.65)GOTO 2
      ANGLE=5.
      IF(IK.GT.90)THEN
        ANGLE=1.
        IK=IK-32
      ENDIF
      KK=CHAR(IK)
      IF(KK.EQ.'E'.OR.KK.EQ.'X')GOTO 5
      IF(KK.EQ.'S'.OR.KK.EQ.'D')GOTO 6
      IF(KK.EQ.'W'.OR.KK.EQ.'R')GOTO 7
      IF(KK.EQ.'N')then
        if(icontr/10.lt.1)then
          icontr=icontr+10
        else
          icontr=icontr-10
        endif
        CALL VIEW(IA,IB,icontr,1)
      endif
      IF(KK.EQ.'A')then
        if(icontr-(icontr/10)*10.eq.0)then
          icontr=icontr+2
        else
          icontr=icontr-2
        endif
        CALL VIEW(IA,IB,icontr,1)
      endif
      IF(KK.EQ.'U')GOTO 800
      IF(KK.EQ.'H')GOTO 500
      IF(KK.eq.'Z')then
        call transf(0.d0,ia,ib,-1)
        goto 901
      endif
c
c...output of cartesian positions of tips of displacement vectors at current
c   weighting, phase and scaling.  These are followed by cartesians
c   for atoms
c
      IF(KK.EQ.'O')THEN
        CALL clearscreen($GVIEWPORT)
        CALL settextposition(3,1,curpos)
        open(4,file='vector.cor',status='unknown')
        write(4,'(a)')coment
        write(4,'(i5)')nat
c
c...calc for molecule with atoms moved by the current displacement vector
        do 444 n=nat+1,natms
          write(4,445)n-nat,corplt(n,1),corplt(n,2),corplt(n,3),
     *          wt(n-nat)
          cartdf(n-nat,4)=wt(n-nat)
          cartdf(n-nat,1)=corplt(n,1)
          cartdf(n-nat,2)=corplt(n,2)
          cartdf(n-nat,3)=corplt(n,3)
444     continue
445     format(i4,4f16.8)
        call momina(nat,1)
        write(4,447)nvectr,fmode(nvectr),iwu(iwtunw)
447     format('|'/'|__Arrowhead coordinates for Mode',I2,':',
     *     F7.1,'cm-1, ',a//'Coordinates for atoms:')
c
c...calc for unmoved molecule
        do 446 n=1,nat
          write(4,445)n,corplt(n,1),corplt(n,2),corplt(n,3),
     *          wt(n)
          cartdf(n,4)=wt(n)
          cartdf(n,1)=corplt(n,1)
          cartdf(n,2)=corplt(n,2)
          cartdf(n,3)=corplt(n,3)
446     continue
        write(4,454)
        write(*,454)
454     format(1x/20x,
     *     '+displacement       unmoved           difference')
        call momina(nat,2)
        write(4,'(a)')'  Moments of inertia:'
        write(*,'(a)')'  Moments of inertia:'
        do 449 n=1,3
          write(4,448)nam(n),pmi(n,1,1),pmi(n,1,2),pmi(n,1,1)-pmi(n,1,2)
          write(*,448)nam(n),pmi(n,1,1),pmi(n,1,2),pmi(n,1,1)-pmi(n,1,2)
449     continue
448     format(12x,a,3f18.8)
        write(4,'(a)')'  Rotational constants:'
        write(*,'(a)')'  Rotational constants:'
        do 451 n=1,3
          write(4,448)nam(n),pmi(n,2,1),pmi(n,2,2),pmi(n,2,1)-pmi(n,2,2)
          write(*,448)nam(n),pmi(n,2,1),pmi(n,2,2),pmi(n,2,1)-pmi(n,2,2)
451     continue
        write(4,'(a)')'  Planar moments:'
        write(*,'(a)')'  Planar moments:'
        do 452 n=1,3
          write(4,448)nam(n),pmi(n,3,1),pmi(n,3,2),pmi(n,3,1)-pmi(n,3,2)
          write(*,448)nam(n),pmi(n,3,1),pmi(n,3,2),pmi(n,3,1)-pmi(n,3,2)
452     continue
        close(4)
        write(*,'(a,$)')char(7)
         call delay(int2(10))
        write(*,'(a,$)')char(7)
c
        write(*,'(1x)')
        WRITE(*,502)
        CALL setviewport(0,linofs+1,maxx-1,maxy-linofs-1)
        CALL moveto(0,0,ixy)
        dummy=lineto(0,maxy-2*linofs)
        CALL setviewport(1,linofs+1,maxx-1,maxy-linofs-1)
453     IK=INKEY(N)
        IF(IK.EQ.0)GOTO 453
        GOTO 3
      ENDIF
C
C...Parameter check
C
      IF(KK.EQ.'P')THEN
        if(icontr/10.lt.1)icontr=icontr+10
        CALL VIEW(IA,IB,icontr,1)
        CALL PARCHE
        GOTO 33
      ENDIF
      GOTO 2
C
C...Exit
C
20    dummy4 = setbkcolor( RED )
      CALL settextposition(1,1,curpos)
      CALL outtext(emplin)
      WRITE(OUTSTR,'(A)')'ARE YOU SURE ?'
      CALL settextposition(1,1,curpos)
      CALL outtext(outstr(1:14))
      CALL settextposition(1,20,curpos)
      WRITE(*,'(1X,A1,$)')CHAR(7)
916   IK=INKEY(N)
      KK=CHAR(IK)
      IF(KK.EQ.'Y'.OR.KK.EQ.'y')GOTO 915
      IF(KK.NE.'N'.AND.KK.NE.'n')GOTO 916
      dummy4 = setbkcolor( BLUE )
      CALL settextposition(1,1,curpos)
      CALL outtext(toplin)
      GOTO 2
915   dummy=setvideomode($DEFAULTMODE)
      GOTO 4
C
C...Help screen
C
500   dummy4 = setbkcolor( BLUE )
      CALL clearscreen($GCLEARSCREEN)
      CALL settextposition(2,1,curpos)
      WRITE(*,232)
232   FORMAT(   '  SUMMARY OF COMMANDS ACTIVE IN GRAPHICS MODE:'/
     *          '  ============================================'//
     *'            U - toggle mass weighted/unweighted eigenvectors'/
     *'       arrows - change mode (up/down) or vector scaling ',
     * '(left/right)'/
     *'  E/N S/D W/R - rotation of molecule about principal axes'/
     *'            Z - zero previous rotations'/
     *'          A/N - toggle principal axes/atom numbering'/
     *'            H - display this help screen'//
     *'            P - check structural parameters. When in this mode'/
     *'                B selects bond length, A bond angle and D dihe-'/
     *'                dral angle. Once parameter type has been')
      write(*,233)
233   format(
     *'                selected numbers of atoms are to be typed in, ',
     * 'one digit'/
     *'                numbers are to be terminated with <ENTER>, two'/
     *'                digit numbers need no termination'/12x,
     *'O - output vector cartesians for current mode, phase,'/
     *'                weighting and scaling to VECTOR.COR'//
     *'  Exit from parameter check mode is with E, exit from graphics '/
     *'  mode is with <ENTER> followed by Y'/)
      WRITE(*,502)
502   FORMAT(1X,40X,' Press ENTER to return to diagram     ',$)
c      CALL setviewport(0,linofs+1,maxx-1,maxy-linofs-1)
c      CALL moveto(0,0,ixy)
c      dummy=lineto(0,maxy-2*linofs)
c      CALL setviewport(1,linofs+1,maxx-1,maxy-linofs-1)
501   IK=INKEY(N)
      IF(IK.EQ.0)GOTO 501
      GOTO 333
C
C...Rotation about the horizontal axis
C
5     IF(IK.EQ.69)ANGLE=-ANGLE
      CALL TRANSF(ANGLE,IB,IC,0)
      GOTO 3
C
C...Rotation about the vertical axis
C
6     IF(IK.EQ.83)ANGLE=-ANGLE
      CALL TRANSF(ANGLE,IA,IC,0)
      GOTO 3
C
C...Rotation about the axis perpendicular to the projection plane
C
7     IF(IK.EQ.87)ANGLE=-ANGLE
      CALL TRANSF(ANGLE,IA,IB,0)
      GOTO 3
C
C...Toggle weighted/unweighted eigenvectors
C
800   IF(IWTUNW.EQ.1)THEN
        IWTUNW=2
      ELSE
        IWTUNW=1
      ENDIF
      GOTO 901
C
C...Toggle through eigenvectors
C
900   IF(IK.EQ.-72)THEN
        NVECTR=NVECTR+1
        IF(NVECTR.GT.NVEC)NVECTR=1
      ELSE
        NVECTR=NVECTR-1
        IF(NVECTR.LT.1)NVECTR=NVEC
      ENDIF
C
901   RR=IPHASE*SCALE
      DO 902 N=1,NAT
        IF(IWTUNW.EQ.1)THEN
          CORD(N+NAT,1)=CORD(N,1)+RR*EW(NVECTR,N,1)
          CORD(N+NAT,2)=CORD(N,2)+RR*EW(NVECTR,N,2)
          CORD(N+NAT,3)=CORD(N,3)+RR*EW(NVECTR,N,3)
        ELSE
          CORD(N+NAT,1)=CORD(N,1)+RR*EU(NVECTR,N,1)
          CORD(N+NAT,2)=CORD(N,2)+RR*EU(NVECTR,N,2)
          CORD(N+NAT,3)=CORD(N,3)+RR*EU(NVECTR,N,3)
        ENDIF
902   CONTINUE
904   DO 903 N=1,NATMS
        DO 903 NN=1,3
          CORPLT(N,NN)=CORD(N,NN)
903   CONTINUE
      do 905 n=1,6
        do 905 nn=1,3
          axend(n,nn)=axendr(n,nn)
905   continue
      call transf(0.0d0,IA,IB,1)
      GOTO 33
C
C...Reverse phase
C
700   IPHASE=-IPHASE
      GOTO 901
C
C...Change vector scaling
C
710   IF(IK.EQ.-75)THEN
        SCALE=SCALE*0.75D0
      ELSE
        SCALE=SCALE*1.333333333333333D0
      ENDIF
      GOTO 901
C
4     RETURN
      END
C
C_____________________________________________________________________________
C
      SUBROUTINE TRANSF(ANGLED,IX,IY,ITYPR)
C
c   ITYPR = -1 zero cumulative rotation matrix ROTMAT
c         =  0 rotate by ANGLED in IX,IY plane and add this rotation
c              to ROTMAT
c         =  1 do complete 3D transformation according to current ROTMAT
C
C
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (maxat=25, maxvec=3*maxat, maxpts=maxat+maxvec,
     *           maxbon=6*maxpts)
      COMMON /BLK1/NATMS,CORD,CORPLT,AXEND,axendr,ASPECT,ROTMAT(3,3)
      REAL*8 CORD(maxpts,3),CORPLT(maxpts,3),AXEND(6,3),axendr(6,3),
     * ASPECT,R(3,3),T(3,3),TT(3)
C
      CON=0.174532925199D-1
      ANGLE=ANGLED*CON
      S=DSIN(ANGLE)
      C=DCOS(ANGLE)
C
c...simple in-plane rotation and update cumulative rotation matrix
c
      if(itypr.eq.0)then
        DO 2 N=1,3
          DO 2 NN=1,3
            IF(N.EQ.NN)THEN
              R(N,NN)=1.D0
            ELSE
              R(N,NN)=0.D0
            ENDIF
2       CONTINUE
        R(IX,IX)= C
        R(IX,IY)=-S
        R(IY,IX)= S
        R(IY,IY)= C
C
        DO 3 N=1,3
          DO 3 NN=1,3
            T(N,NN)=0.D0
            DO 3 NNN=1,3
              T(N,NN)=T(N,NN)+R(N,NNN)*ROTMAT(NNN,NN)
3       CONTINUE
c
        DO 4 N=1,3
          DO 4 NN=1,3
            ROTMAT(N,NN)=T(N,NN)
4       CONTINUE
c
        DO 1 N=1,NATMS
          A=C*CORPLT(N,IX)-S*CORPLT(N,IY)
          B=S*CORPLT(N,IX)+C*CORPLT(N,IY)
          CORPLT(N,IX)=A
          CORPLT(N,IY)=B
1       CONTINUE
c
        DO 12 N=1,6
          A=C*AXEND(N,IX)-S*AXEND(N,IY)
          B=S*AXEND(N,IX)+C*AXEND(N,IY)
          AXEND(N,IX)=A
          AXEND(N,IY)=B
12      CONTINUE
      endif
c
C...zero cumulative rotation matrix
C
      if(itypr.eq.-1)then
        DO 10 N=1,3
          DO 10 NN=1,3
            IF(N.EQ.NN)THEN
              ROTMAT(N,NN)=1.D0
            ELSE
              ROTMAT(N,NN)=0.D0
            ENDIF
10      CONTINUE
      ENDIF
C
C...3D rotation according to cumulative rotation matrix
C
      if(itypr.eq.1)then
        DO 51 N=1,NATMS
          DO 52 nn=1,3
            tt(nn)=rotmat(nn,1)*corplt(n,1)+rotmat(nn,2)*corplt(n,2)+
     *             rotmat(nn,3)*corplt(n,3)
52        CONTINUE
          DO 53 nn=1,3
            corplt(n,nn)=tt(nn)
53        CONTINUE
51      CONTINUE
c
        DO 54 N=1,6
          DO 55 nn=1,3
            tt(nn)=rotmat(nn,1)*axend(n,1)+rotmat(nn,2)*axend(n,2)+
     *             rotmat(nn,3)*axend(n,3)
55        CONTINUE
          DO 56 nn=1,3
            axend(n,nn)=tt(nn)
56        CONTINUE
54      CONTINUE
      ENDIF
C
      RETURN
      END
C
C_____________________________________________________________________________
C
      SUBROUTINE PARCHE
C
C   Check of structural parameters (bond lengths and bond angles) using
C   a scrollable window for listing the results
C
      INCLUDE  'FGRAPH.FD'                                              MSF5+PS1
c     USE MSFLIB                                                        PS4
      RECORD /rccoord/curpos
      RECORD /xycoord/ixy
      RECORD /wxycoord/wxy
C
      PARAMETER (maxat=25, maxvec=3*maxat, maxpts=maxat+maxvec,
     *           maxbon=6*maxpts)
c
      COMMON /limits/maxx,maxy,LINOFS,curpos,ixy,wxy,
     * mymode,myrows,mycols,blue,red
      COMMON /BLK1/NATMS,CORD,CORPLT,AXEND,axendr,ASPECT,ROTMAT
      COMMON /PARBLK/RANGEA,RANGEB,AMIN,BMIN
      COMMON /lines/toplin,botlin,emplin,COMENT
      CHARACTER KK,outlin*20,TEXT*20,RTEXT*21
      REAL*8 BONDL,ANGLE,DANGLE,CORD(maxpts,3),CORPLT(maxpts,3),
     * AXEND(6,3),axendr(6,3),ROTMAT(3,3)
      REAL*8 ANG,CON,COSANG,RA,RANGEA,RANGEB,
     1 ASPECT,AMIN,BMIN
      INTEGER*2 dummy,maxx,maxy,mymode,myrows,mycols,LINOFS
      INTEGER*4 blue,red
      INTEGER*2 NPOS,NBUF,inkey
      CHARACTER*79 toplin,botlin,emplin,COMENT
      CHARACTER*36 RES(50),CURRES,BLANK
C
      CON=0.174532925199D-1
      NBUF=0
      rtext=CHAR(24)//CHAR(25)//'=scrollback  |E|xit'
      WRITE(BLANK,'(36(1H ))')
      CALL settextposition(myrows,1,curpos)
      CALL outtext(emplin)
C
C
C   INPUT LOOP
C
C...selection of B (bond), A (angle), D (dihedral angle),
C   scrollback or E (exit)
C
3     CALL settextposition(1,1,curpos)
      CALL outtext(emplin)
      CALL settextposition(1,1,curpos)
      CALL outtext('|B|ond,|A|ngle,|D|ihangle?')
      CALL settextposition(1,58,curpos)
      CALL outtext(rtext)
33    IK=INKEY(dummy)
      IF(IK.EQ.0)GOTO 33
      IF(IK.GT.90)IK=IK-32
      IF(IK.EQ.-72.OR.IK.EQ.-80)GOTO 8
      KK=CHAR(IK)
      IF(KK.NE.'A'.AND.KK.NE.'B'.AND.KK.NE.'D'.AND.KK.NE.'E')GOTO 4
      IF(KK.EQ.'E')GOTO 5
C
      CALL settextposition(1,1,curpos)
      CALL outtext(emplin)
      CALL settextposition(1,1,curpos)
      CALL outtext('---- Type in numbers of atoms or a,b,c ----')
      CALL settextposition(myrows,1,curpos)
      CALL outtext(BLANK)
      CALL settextposition(myrows,1,curpos)
      IF(KK.EQ.'B')THEN
        TEXT='----------Bond('
      ELSE IF(KK.EQ.'A')THEN
        TEXT='---------Angle('
      ELSE IF(KK.EQ.'D')THEN
        TEXT='Dihedral Angle('
      ENDIF
      CALL outtext(TEXT)
      CURRES=TEXT
      NPOS=16
C
C...read first two atom numbers, and calculate the relevant bond length
      CALL INNUM(N,NPOS,CURRES)
      CALL settextposition(myrows,NPOS,curpos)
      CALL outtext(',')
      CURRES=CURRES(1:(NPOS-1))//','
      NPOS=NPOS+1
      CALL INNUM(NN,NPOS,CURRES)
      IF(N.EQ.NN)GOTO 3
      IF(N.LT.0.OR.NN.LT.0)GOTO 3
      IF(KK.EQ.'B')THEN
        RA=BONDL(CORD(N,1),CORD(N,2),CORD(N,3),
     *           CORD(NN,1),CORD(NN,2),CORD(NN,3))
        WRITE(outlin,3104)RA
3104    FORMAT(') =',F8.5,$)
        CALL settextposition(myrows,NPOS,curpos)
        CALL outtext(outlin)
        GOTO 6
      ENDIF
C
C...third atom number for angle and dihedral angle
C   (if angle selected then a,b or c are also allowed on input for
C    determining angles between bonds and principal axes)
C
      CALL settextposition(myrows,NPOS,curpos)
      CALL outtext(',')
      CURRES=CURRES(1:(NPOS-1))//','
      NPOS=NPOS+1
      CALL INNUM(NNN,NPOS,CURRES)
      IF(NNN.EQ.N.OR.NNN.EQ.NN)GOTO 3
      IF(NNN.LT.0.AND.KK.NE.'A')GOTO 3
      IF(KK.EQ.'A')THEN
C
C...angle between bond and axis
        IF(NNN.LT.0)THEN
          NNN=-NNN
          RA=BONDL(CORD(N,1),CORD(N,2),CORD(N,3),
     *             CORD(NN,1),CORD(NN,2),CORD(NN,3))
          COSANG=(CORD(NN,NNN)-CORD(N,NNN))/RA
          ANG=DACOS(COSANG)/CON
C
C...angle between two bonds
        ELSE
         ANG=ANGLE(CORD(N,1),CORD(N,2),CORD(N,3),
     *           CORD(NN,1),CORD(NN,2),CORD(NN,3),
     *           CORD(NNN,1),CORD(NNN,2),CORD(NNN,3))
        ENDIF
C
        WRITE(outlin,3103)ANG
3103    FORMAT(') =',F8.3,$)
        CALL settextposition(myrows,NPOS,curpos)
        CALL outtext(outlin)
        GOTO 6
      ENDIF
C
C...fourth atom number for dihedral angle
      CALL settextposition(myrows,NPOS,curpos)
      CALL outtext(',')
      CURRES=CURRES(1:(NPOS-1))//','
      NPOS=NPOS+1
      CALL INNUM(NNNN,NPOS,CURRES)
      IF(NNNN.EQ.N.OR.NNNN.EQ.NN.OR.NNNN.EQ.NNN)GOTO 3
      IF(NNNN.LT.0)GOTO 3
      ANG=DANGLE(CORD(N,1),CORD(N,2),CORD(N,3),
     *           CORD(NN,1),CORD(NN,2),CORD(NN,3),
     *           CORD(NNN,1),CORD(NNN,2),CORD(NNN,3),
     *           CORD(NNNN,1),CORD(NNNN,2),CORD(NNNN,3))
      WRITE(outlin,3103)ANG
      CALL settextposition(myrows,NPOS,curpos)
      CALL outtext(outlin)
      GOTO 6
C
C...Beep for input error and return to option selection
C
4     CALL settextposition(1,1,curpos)
      WRITE(*,'(1X,A1,$)')CHAR(7)
      GOTO 3
C
C...Update the scrollback buffer
C
6     CURRES=CURRES(1:(NPOS-1))//outlin
      DO 7 N=19,1,-1
        RES(N+1)=RES(N)
7     CONTINUE
      RES(1)=CURRES
      NBUF=NBUF+1
      GOTO 3
C
C...Scroll through the results buffer
C
8     IF(IK.EQ.-72)THEN
        NBUF=NBUF+1
        IF(NBUF.GT.50)NBUF=1
      ELSE
        NBUF=NBUF-1
        IF(NBUF.LT.1)NBUF=50
      ENDIF
      CALL settextposition(myrows,43,curpos)
      CALL outtext(BLANK)
      CALL settextposition(myrows,43,curpos)
      CALL outtext(RES(NBUF))
      GOTO 3
C
5     RETURN
      END
C
C_____________________________________________________________________________
C
      FUNCTION BONDL(A1,A2,A3,B1,B2,B3)
C
C   Bond (A.B)
C
      REAL*8 A1,A2,A3,B1,B2,B3,BONDL
C
      BONDL=DSQRT((A1-B1)**2+(A2-B2)**2+(A3-B3)**2)
C
      RETURN
      END
C_____________________________________________________________________________
C
      FUNCTION ANGLE(A1,A2,A3,B1,B2,B3,C1,C2,C3)
C
C   Angle (A.B.C) using cosine rule, returned in degrees
C
      REAL*8 A1,A2,A3,B1,B2,B3,C1,C2,C3,CONV,ANGLE,BONDL
      REAL*8 BL1,BL2,BL3
      CONV=1.745329252D-2
C
      BL1=BONDL(A1,A2,A3,B1,B2,B3)
      BL2=BONDL(B1,B2,B3,C1,C2,C3)
      BL3=BONDL(A1,A2,A3,C1,C2,C3)
      ANGLE=(BL1**2+BL2**2-BL3**2)/(2.D0*BL1*BL2)
      IF(DABS(ANGLE).GT.1.D0)ANGLE=DSIGN(1.D0,ANGLE)
      ANGLE=DACOS(ANGLE)/CONV
C
      RETURN
      END
C_____________________________________________________________________________
C
      FUNCTION DANGLE(A1,A2,A3,B1,B2,B3,C1,C2,C3,D1,D2,D3)
C
C   Dihedral angle (A.B.C.D), returned in degrees
C
      REAL*8 A1,A2,A3,B1,B2,B3,C1,C2,C3,D1,D2,D3,
     *      CONV,DANGLE,ANGLE,BONDL
      REAL*8 ANGL1,ANGL2,P1,P2,P3,Q1,Q2,Q3,R12S,R23,
     * R34S,BL1,BL2,DOTP
      CONV=1.745329252D-2
C
      ANGL1=ANGLE(A1,A2,A3,B1,B2,B3,C1,C2,C3)-90.D0
      ANGL2=ANGLE(B1,B2,B3,C1,C2,C3,D1,D2,D3)-90.D0
      R12S=BONDL(A1,A2,A3,B1,B2,B3)*DSIN(ANGL1*CONV)
      R23=BONDL(B1,B2,B3,C1,C2,C3)
      R34S=BONDL(C1,C2,C3,D1,D2,D3)*DSIN(ANGL2*CONV)
C
      P1=A1-B1+R12S*(C1-B1)/R23
      P2=A2-B2+R12S*(C2-B2)/R23
      P3=A3-B3+R12S*(C3-B3)/R23
      Q1=D1-C1-R34S*(C1-B1)/R23
      Q2=D2-C2-R34S*(C2-B2)/R23
      Q3=D3-C3-R34S*(C3-B3)/R23
C
      BL1=DSQRT(P1**2+P2**2+P3**2)
      BL2=DSQRT(Q1**2+Q2**2+Q3**2)
      DOTP=P1*Q1+P2*Q2+P3*Q3
C
      DANGLE=DOTP/(BL1*BL2)
      IF(DABS(DANGLE).GT.1.D0)DANGLE=DSIGN(1.D0,DANGLE)
      DANGLE=DACOS(DANGLE)/CONV
C
      RETURN
      END
C_____________________________________________________________________________
C
      SUBROUTINE INNUM(N,NPOS,CURRES)
C
C   Input and check of validity of a 1 to 2 digit integer number;
C      N - takes on the value of the input atom number (or -1,-2,-3 for
C          a,b,c axis in the case of angle between a bond and a principal
C          axis)
C   NPOS - starting column for writing of next number to bottom line of
C          screen, updated by INNUM after the number is written
C CURRES - buffer containing the current result line
C
      INCLUDE  'FGRAPH.FD'                                              MSF5+PS1
c     USE MSFLIB                                                        PS4
      RECORD /rccoord/curpos
      RECORD /xycoord/ixy
      RECORD /wxycoord/wxy
C
      PARAMETER (maxat=25, maxvec=3*maxat, maxpts=maxat+maxvec,
     *           maxbon=6*maxpts)
c
      COMMON /limits/maxx,maxy,LINOFS,curpos,ixy,wxy,
     * mymode,myrows,mycols,blue,red
      COMMON /BLK1/NATMS,COR,CORPLT,AXEND,axendr,ASPECT,ROTMAT
      CHARACTER TEXT*2,KK,CURRES*36
      INTEGER*2 dummy,maxx,maxy,mymode,myrows,mycols,LINOFS
      INTEGER*4 blue,red
      INTEGER*2 inkey,NDIG,IK,NPOS
      REAL*8 COR(maxpts,3),CORPLT(maxpts,3),AXEND(6,3),axendr(6,3),
     * ASPECT,ROTMAT(3,3)
      NDIG=0
C
C...Accept input
C
      CALL settextposition(myrows,NPOS,curpos)
5     IK=inkey(dummy)
      IF(IK.EQ.0)GOTO 5
      KK=CHAR(IK)
C
C...identify axis
      IF(KK.EQ.'A'.OR.KK.EQ.'a')THEN
        CALL outtext('a')
        CURRES=CURRES(1:(NPOS-1))//'a'
        NPOS=NPOS+1
        N=-1
        GOTO 4
      ENDIF
      IF(KK.EQ.'B'.OR.KK.EQ.'b')THEN
        CALL outtext('b')
        CURRES=CURRES(1:(NPOS-1))//'b'
        NPOS=NPOS+1
        N=-2
        GOTO 4
      ENDIF
      IF(KK.EQ.'C'.OR.KK.EQ.'c')THEN
        CALL outtext('c')
        CURRES=CURRES(1:(NPOS-1))//'c'
        NPOS=NPOS+1
        N=-3
        GOTO 4
      ENDIF
C
C...identify number (one or two digit)
      IF(IK.EQ.13.AND.NDIG.EQ.1)THEN
        N=N1
        KK=CHAR(N+48)
        IF(N.GT.NATMS)THEN
          WRITE(*,'(1X,A1,$)')CHAR(7)
          NDIG=0
          GOTO 5
        ENDIF
        CALL outtext(KK)
        CURRES=CURRES(1:(NPOS-1))//KK
        NPOS=NPOS+1
        GOTO 4
      ENDIF
      IF(KK.LT.'0'.OR.KK.GT.'9')GOTO 5
      IF(NDIG.EQ.0)THEN
        N1=IK-48
        NDIG=1
        GOTO 5
      ENDIF
      IF(NDIG.EQ.1)THEN
        NDIG=2
        N=N1*10+(IK-48)
        IF(N.GT.NATMS)THEN
          WRITE(*,'(1X,A1,$)')CHAR(7)
          NDIG=0
          GOTO 5
        ENDIF
        WRITE(TEXT,'(I2)')N
        CALL outtext(TEXT)
        CURRES=CURRES(1:(NPOS-1))//TEXT
        NPOS=NPOS+2
        GOTO 4
      ENDIF
C
4     RETURN
      END
C
C_____________________________________________________________________________
C
      SUBROUTINE PLOTCH(XSCALE,YSCALE,ISTART,JSTART,IDIR)
C
C    This is a simple plotter type character generator which plots the
C    specified character string LINE starting from absolute pixel coordinates
C    (ISTART,JSTART)  The font is separately scalable in horizontal and
C    vertical directions and choice of upright type or italics is possible.
C
C    LINE - text string to be written
C    XSCALE - integer scaling factor for the horizontal character dimension
C    YSCALE - integer scaling factor for the vertical character dimension
C    ISTART, JSTART - starting coordinates for placing the text in units of
C               the output device
C    IDIR   - direction of writing the character string
C             = +-1 left to right, = +-2 from bottom to top
C             (negative values specify italics)
C
C
      INCLUDE  'FGRAPH.FD'                                              MSF5+PS1
c     USE MSFLIB                                                        PS4
      RECORD /rccoord/curpos
      RECORD /xycoord/ixy
      RECORD /wxycoord/wxy
C
      COMMON /limits/maxx,maxy,LINOFS,curpos,ixy,wxy,
     * mymode,myrows,mycols,blue,red
      COMMON /TXTBUF/LINE
      EQUIVALENCE (LINE,TEXT(1))
      EQUIVALENCE(CHARPL(1,1),ACHAR(1,1)),(CHARPL(1,11),BCHAR(1,1)),
     1         (CHARPL(1,21),CCHAR(1,1)),(CHARPL(1,31),DCHAR(1,1)),
     1         (CHARPL(1,41),ECHAR(1,1)),(CHARPL(1,51),FCHAR(1,1)),
     1         (CHARPL(1,61),GCHAR(1,1)),(CHARPL(1,71),HCHAR(1,1)),
     1         (CHARPL(1,81),JCHAR(1,1)),(CHARPL(1,91),KCHAR(1,1))
      CHARACTER*80 LINE
      CHARACTER TEXT(80),NUL,SPACE
      INTEGER*2 dummy,maxx,maxy,mymode,myrows,mycols,LINOFS
      INTEGER*4 blue,red
      INTEGER*4 XSCALE,YSCALE,ISTART,JSTART,IDIR
      INTEGER*2 NPLOT(95),IX,IY,IXX,IYY
      INTEGER*2 ACHAR(22,10),BCHAR(22,10),CCHAR(22,10),DCHAR(22,10),
     1 ECHAR(22,10),FCHAR(22,10),GCHAR(22,10),HCHAR(22,10),
     1 JCHAR(22,10),KCHAR(22,5),CHARPL(22,95)
C
      SLANT=0.
      IDIREC=IABS(IDIR)
      IXSTRT=ISTART
      IYSTRT=JSTART
      WRITE(NUL,'(A1)')CHAR(0)
      SPACE=' '
      DO 100 NCHARS=80,1,-1
100     IF(TEXT(NCHARS).NE.NUL.AND.TEXT(NCHARS).NE.SPACE)GOTO 101
C
101   GO TO (31,32),IDIREC
31    SCALEX=XSCALE
      SCALEY=YSCALE
      ISPACE=3*XSCALE
      IF(IDIR.LT.0)SLANT=0.5*SCALEY
      GO TO 30
32    SCALEX=YSCALE*1.5
      SCALEY=XSCALE
      IF(SCALEY.GT.1.)SCALEY=SCALEY/1.5
      ISPACE=2*XSCALE
      IF(IDIR.LT.0)SLANT=0.334*SCALEY
30    CONTINUE
C
C   CHARACTER DEFINITIONS - NPLOT contains the number of points in character
C      definition, CHARPL contains the definitions
C
C   ASCII codes from 32 to 126 are supported with their usual meanings
C   with the following exceptions:
C           ^ ( 94) - shiift up for superscripts
C           ~ (126) - shift down for subscripts
C           ` ( 96) - one half (1 over 2) character
C           @ ( 64) - arrow to the left ie. <-
C
C   _!"#$%&'()   *+,-./0123   456789:;<=   >?@ABCDEFG   HIJKLMNOPQ
C
C   RSTUVWXYZ[   \]^_`abcde   fghijklmno   pqrstuvwxy   z{|}~
C
      DATA NPLOT/6,4,6,8,10,1,11,6,4,4,
     *           10,4,6,2,5,2,9,5,7,9,
     4           5,8,11,4,11,11,10,11,3,4,
     >           3,10,10,5,11,8,7,6,5,10,
     H           6,6,5,7,3,5,4,9,7,11,
     R           9,10,4,6,3,5,4,5,4,4,
     \           2,4,1,2,6,9,10,8,11,9,
     f           6,10,7,7,6,7,5,10,7,9,
     p           10,10,6,8,6,8,3,11,4,6,
     z           4,9,2,9,1/
C
      DATA ACHAR/0,0,1,6,3,6,4,5,4,4,3,3,0,0,0,0,0,0,0,0,0,0,            space
     ! 0,6,0,2,100,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
     " 1,6,1,5,0,4,103,6,3,5,2,4,0,0,0,0,0,0,0,0,0,0,
     # 1,6,1,0,103,6,3,0,100,4,4,4,100,2,4,2,0,0,0,0,0,0,
     $ 2,6,2,0,100,1,3,1,4,2,3,3,1,3,0,4,1,5,4,5,0,0,
     % 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
     & 4,0,0,4,0,5,1,6,2,5,2,4,0,2,0,1,1,0,2,0,4,2,
     ' 1,5,0,5,0,6,1,6,1,4,0,3,0,0,0,0,0,0,0,0,0,0,
     ( 2,7,0,5,0,1,2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
     ) 0,7,2,5,2,1,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0/
      DATA BCHAR/1,5,2,4,2,2,1,1,103,5,2,4,102,2,3,1,100,3,4,3,0,0,      *
     + 0,3,4,3,102,5,2,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
     , 1,0,0,0,0,1,1,1,1,-1,0,-2,0,0,0,0,0,0,0,0,0,0,
     - 0,3,4,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
     . 0,0,0,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
     / 0,0,4,6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
     O 1,0,3,0,4,1,4,5,3,6,1,6,0,5,0,1,1,0,0,0,0,0,
     1 0,0,2,0,101,0,1,6,0,5,0,0,0,0,0,0,0,0,0,0,0,0,
     2 4,0,0,0,4,4,4,5,3,6,1,6,0,5,0,0,0,0,0,0,0,0,
     3 0,6,4,6,2,4,3,4,4,3,4,1,3,0,1,0,0,1,0,0,0,0/
      DATA CCHAR/3,0,3,6,0,3,0,2,4,2,0,0,0,0,0,0,0,0,0,0,0,0,            4
     5 4,6,0,6,0,3,3,3,4,2,4,1,3,0,0,0,0,0,0,0,0,0,
     6 3,6,2,6,0,4,0,1,1,0,3,0,4,1,4,2,3,3,1,3,0,2,
     7 0,6,4,6,1,2,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
     8 1,3,3,3,4,2,4,1,3,0,1,0,0,1,0,2,1,3,0,4,0,5,
     9 1,0,2,0,4,2,4,5,3,6,1,6,0,5,0,4,1,3,3,3,4,4,
     : 0,4,0,5,1,5,1,4,0,4,100,0,1,0,1,1,0,1,0,0,0,0,
     ; 0,4,0,5,1,5,1,4,0,4,101,0,0,0,0,1,1,1,1,-1,0,-2,
     < 3,6,0,3,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
     = 0,4,4,4,100,2,4,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0/
      DATA DCHAR/0,6,3,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,            >
     ? 0,5,1,6,3,6,4,5,4,4,3,3,2,3,2,2,102,0,2,0,0,0,
     @ 0,7,1,8,1,5,100,3,3,3,100,1,1,1,2,0,0,-2,2,-2,0,0,
     A 0,0,2,6,4,0,101,2,3,2,0,0,0,0,0,0,0,0,0,0,0,0,
     B 0,3,3,3,4,4,4,5,3,6,0,6,0,0,3,0,4,1,4,2,3,3,
     C 4,5,3,6,1,6,0,5,0,1,1,0,3,0,4,1,0,0,0,0,0,0,
     D 0,0,3,0,4,1,4,5,3,6,0,6,0,0,0,0,0,0,0,0,0,0,
     E 4,0,0,0,0,6,4,6,100,3,2,3,0,0,0,0,0,0,0,0,0,0,
     F 0,0,0,6,4,6,100,3,2,3,0,0,0,0,0,0,0,0,0,0,0,0,
     G 4,5,3,6,1,6,0,5,0,1,1,0,3,0,4,1,4,3,2,3,0,0/
      DATA ECHAR/0,0,0,6,100,3,4,3,104,6,4,0,0,0,0,0,0,0,0,0,0,0,       H
     I 0,0,2,0,101,0,1,6,100,6,2,6,0,0,0,0,0,0,0,0,0,0,
     J 0,0,2,0,3,1,3,6,0,6,0,0,0,0,0,0,0,0,0,0,0,0,
     K 0,0,0,6,100,2,2,3,4,6,102,3,4,0,0,0,0,0,0,0,0,0,
     L 0,6,0,0,4,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
     M 0,0,0,6,3,3,6,6,6,0,0,0,0,0,0,0,0,0,0,0,0,0,
     N 0,0,0,6,4,0,4,6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
     O 1,0,3,0,4,1,4,5,3,6,1,6,0,5,0,1,1,0,0,0,0,0,
     P 0,0,0,6,3,6,4,5,4,3,3,2,0,2,0,0,0,0,0,0,0,0,
     Q 1,0,3,0,4,1,4,5,3,6,1,6,0,5,0,1,1,0,103,1,4,-1/
      DATA FCHAR/0,0,0,6,3,6,4,5,4,3,3,2,0,2,102,2,4,0,0,0,0,0,         R
     S 4,6,1,6,0,5,0,4,1,3,3,3,4,2,4,1,3,0,0,0,0,0,
     T 2,0,2,6,100,6,4,6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
     U 0,6,0,1,1,0,3,0,4,1,4,6,0,0,0,0,0,0,0,0,0,0,
     V 0,6,2,0,4,6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
     W 0,6,0,0,3,3,6,0,6,6,0,0,0,0,0,0,0,0,0,0,0,0,
     X 0,6,4,0,100,0,4,6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
     Y 2,0,2,3,0,6,104,6,2,3,0,0,0,0,0,0,0,0,0,0,0,0,
     Z 4,0,0,0,4,6,0,6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
     [ 2,7,0,7,0,-1,2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0/
      DATA GCHAR/0,6,4,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,            \
     ] 0,-1,2,-1,2,7,0,7,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
     ^ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
     _ 0,-2,4,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
     ` 0,5,1,5,1,6,0,6,0,4,1,3,0,0,0,0,0,0,0,0,0,0,
     a 4,1,3,0,1,0,0,1,0,3,1,4,4,4,4,1,5,0,0,0,0,0,
     b 0,6,0,0,100,1,1,0,3,0,4,1,4,3,3,4,1,4,0,3,0,0,
     c 4,1,3,0,1,0,0,1,0,3,1,4,3,4,4,3,0,0,0,0,0,0,
     d 4,6,4,1,5,0,104,1,3,0,1,0,0,1,0,3,1,4,3,4,4,3,
     e 3,0,1,0,0,1,0,3,1,4,3,4,4,3,4,2,0,2,0,0,0,0/
      DATA HCHAR/1,0,1,5,2,6,3,5,100,3,2,3,0,0,0,0,0,0,0,0,0,0,         f
     g 0,-2,3,-2,4,-1,4,4,1,4,0,3,0,1,1,0,3,0,4,1,0,0,
     h 0,0,0,6,100,3,1,4,3,4,4,3,4,0,0,0,0,0,0,0,0,0,
     i 0,4,1,4,1,0,100,0,2,0,101,6,1,6,0,0,0,0,0,0,0,0,
     j 0,-2,1,-2,2,-1,2,4,102,6,2,6,0,0,0,0,0,0,0,0,0,0,
     k 0,0,0,6,100,1,2,2,4,4,102,2,4,0,0,0,0,0,0,0,0,0,
     l 0,6,1,6,1,0,100,0,2,0,0,0,0,0,0,0,0,0,0,0,0,0,
     m 0,0,0,4,100,3,2,4,3,4,3,0,103,3,5,4,6,4,6,0,0,0,
     n 0,0,0,4,100,3,2,4,3,4,4,3,4,0,0,0,0,0,0,0,0,0,
     o 1,0,3,0,4,1,4,3,3,4,1,4,0,3,0,1,1,0,0,0,0,0/
      DATA JCHAR/0,-2,0,4,100,1,1,0,3,0,4,1,4,3,3,4,1,4,0,3,0,0,        p
     q 4,-2,4,4,104,3,3,4,1,4,0,3,0,1,1,0,3,0,4,1,0,0,
     r 0,0,0,4,100,2,2,4,3,4,4,3,0,0,0,0,0,0,0,0,0,0,
     s 4,4,1,4,0,3,1,2,3,2,4,1,3,0,0,0,0,0,0,0,0,0,
     t 1,6,1,1,2,0,3,0,100,4,2,4,0,0,0,0,0,0,0,0,0,0,
     u 0,4,0,1,1,0,3,0,4,1,4,4,104,1,5,0,0,0,0,0,0,0,
     v 0,4,2,0,4,4,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
     w 0,4,0,1,1,0,2,0,3,1,3,4,103,1,4,0,5,0,6,1,6,4,
     x 0,0,4,4,100,4,4,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
     y 0,4,2,0,104,4,2,0,1,-2,0,-2,0,0,0,0,0,0,0,0,0,0/
      DATA KCHAR/0,4,4,4,0,0,4,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,           z
     { 3,7,2,7,1,6,1,4,0,3,1,2,1,0,2,-1,3,-1,0,0,0,0,
     | 0,6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
     } 0,-1,1,-1,2,0,2,2,3,3,2,4,2,6,1,7,0,7,0,0,0,0,
     ~ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/
C
C
C   character identification and writing loop
C
      ILIFT=0
      DO 3 N=1,NCHARS
        IK=ICHAR(TEXT(N))
        IF(IK.GT.32.AND.IK.LE.126)GOTO 6
        IF(IK.NE.32)GO TO 3
C
C...... blank character
        GO TO (33,34),IDIREC
33      IXSTRT=IXSTRT+5*SCALEX
        GO TO 3
34      IYSTRT=IYSTRT-5*SCALEY
        GO TO 3
C
6       NN=IK-31
        IF(IK.NE.126.AND.IK.NE.94)GOTO 5
C
C...... subscripts and superscripts
        ISHIFT=4
        IF(IK.EQ.126)ISHIFT=-4
        GOTO (115,116),IDIREC
115     IYSTRT=IYSTRT-ISHIFT*SCALEY
        IXSTRT=IXSTRT+ISPACE*SCALEX*ISHIFT/16
        GOTO 3
116     IXSTRT=IXSTRT-ISHIFT*SCALEX
        IYSTRT=IYSTRT-ISPACE*SCALEY*ISHIFT/16
        GOTO 3
C
5       NPTS=NPLOT(NN)
C
C..... move to beginning of character
C
        GO TO (35,36),IDIREC
35      IX=IXSTRT+CHARPL(1,NN)*SCALEX+CHARPL(2,NN)*SLANT
        IY=IYSTRT-CHARPL(2,NN)*SCALEY
        GO TO 37
36      IX=IXSTRT-CHARPL(2,NN)*SCALEX
        IY=IYSTRT-CHARPL(1,NN)*SCALEY-CHARPL(2,NN)*SLANT
37      call moveto(ix,iy,ixy)
C
C...... go through points defining the character
C
114     DO 7 MN=3,2*NPTS,2
          I=CHARPL(MN,NN)
          J=CHARPL(MN+1,NN)
          IF(I.LT.90)GOTO 110
          ILIFT=1
          I=I-100
110       GO TO (38,39),IDIREC
38        IXX=IXSTRT+I*SCALEX+J*SLANT
          IYY=IYSTRT-J*SCALEY
          GO TO 40
39        IXX=IXSTRT-J*SCALEX
          IYY=IYSTRT-I*SCALEY-J*SLANT
C
C...... move pen to next point of character definition
C
40        IF(ILIFT.NE.1)GOTO 111
          call moveto(IXX,IYY,ixy)
          ILIFT=0
          GOTO 7
111       dummy=lineto(IXX,IYY)
C
7       CONTINUE
C
C   Backtrack for remainder of character 8 (16 defining points)
C
      IF(IK.NE.56)GOTO 112
      IK=32
      NN=1
      NPTS=6
      GOTO 114
C
C   Find starting position for next character
C
112     LARGE=0
        DO 9 MMNN=1,2*NPTS,2
          J=CHARPL(MMNN,NN)
          IF(J.GT.90)J=J-100
9         IF(J.GT.LARGE)LARGE=J
        GO TO(41,42),IDIREC
41      IXSTRT=IXSTRT+ISPACE+LARGE*SCALEX
        GO TO 3
42      IYSTRT=IYSTRT-ISPACE-LARGE*SCALEY
C
3     CONTINUE
C
C
      RETURN
      END
C_____________________________________________________________________________
C
      SUBROUTINE delay( wait )
C
C   Pause for the number of 1/100 seconds specified by WAIT.
C
      INTEGER*2 wait,  kount, dummy
      INTEGER*2 tick1, tick2, tick
c
      kount = 0
      CALL GETTIM( dummy, dummy, dummy, tick1)
      DO WHILE( kount .LT. wait )
         CALL GETTIM( dummy, dummy, dummy,tick2)
         tick = tick2 - tick1
         IF( tick .LT. 0 ) tick = tick + 100
         tick1 = tick2
         kount = kount + tick
      END DO
c
      RETURN
      END
C
C_____________________________________________________________________________
C
      SUBROUTINE MOMINA(natms,icalc)
C
C     Principal moments of inertia of a molecule from its cartesians and
c     masses supplied via CARTDF in the COMMON block /MOMIA/
C
C     natms = number of atoms
c     icalc = 1 or 2 defines the output data set (last index of PMI)
c       PMI = first index defines a,b,c
c           = second index defines moment of inertia, rotational constant,
c             or planar moment
c           = third index defines number of data set
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      PARAMETER (maxat=25, maxvec=3*maxat, maxpts=maxat+maxvec,
     *           maxbon=6*maxpts)
c
      COMMON /MOMIA/CARTDF(maxvec,4),pmi(3,3,2)
      REAL*8 CORD(maxvec,3),A(6),EIGVEC(3,3)
      CON=0.1745329252199D-1
C
C...Transformation to centre of mass coordinates
C
6000  WW=0.0D0
      DO 991 I=1,NATMS
        WW=WW+CARTDF(I,4)
        cord(i,1)=cartdf(i,1)
        cord(i,2)=cartdf(i,2)
        cord(i,3)=cartdf(i,3)
991   CONTINUE
      FMX=0.0D0
      FMY=0.0D0
      FMZ=0.0D0
      DO 992 N=1,NATMS
        FMX=FMX+CARTDF(N,4)*CORD(N,1)
        FMY=FMY+CARTDF(N,4)*CORD(N,2)
        FMZ=FMZ+CARTDF(N,4)*CORD(N,3)
992   CONTINUE
      DO 993 N=1,NATMS
        CORD(N,1)=CORD(N,1)-FMX/WW
        CORD(N,2)=CORD(N,2)-FMY/WW
        CORD(N,3)=CORD(N,3)-FMZ/WW
993   CONTINUE
C
C   Principal moments of inertia, rotational constants and planar moments
C
      DO 994 N=1,6
        A(N)=0.0D0
994   CONTINUE
      DO 995 N=1,NATMS
        A(1)=A(1)+CARTDF(N,4)*(CORD(N,2)**2+CORD(N,3)**2)
        A(2)=A(2)-CARTDF(N,4)*CORD(N,1)*CORD(N,2)
        A(3)=A(3)+CARTDF(N,4)*(CORD(N,1)**2+CORD(N,3)**2)
        A(4)=A(4)-CARTDF(N,4)*CORD(N,1)*CORD(N,3)
        A(5)=A(5)-CARTDF(N,4)*CORD(N,2)*CORD(N,3)
        A(6)=A(6)+CARTDF(N,4)*(CORD(N,1)**2+CORD(N,2)**2)
995   CONTINUE
C
      CALL EIGEN(A,EIGVEC,3,0)
C
      PMI(1,1,icalc)=A(6)
      PMI(2,1,icalc)=A(3)
      PMI(3,1,icalc)=A(1)
      PMI(1,2,icalc)=505379.07D0/PMI(1,1,icalc)
      PMI(2,2,icalc)=505379.07D0/PMI(2,1,icalc)
      PMI(3,2,icalc)=505379.07D0/PMI(3,1,icalc)
      PMI(1,3,icalc)=0.5d0*
     1               (pmi(2,1,icalc)+pmi(3,1,icalc)-pmi(1,1,icalc))
      PMI(2,3,icalc)=0.5d0*
     1               (pmi(1,1,icalc)+pmi(3,1,icalc)-pmi(2,1,icalc))
      PMI(3,3,icalc)=0.5d0*
     1               (pmi(1,1,icalc)+pmi(2,1,icalc)-pmi(3,1,icalc))
C
      RETURN
      END
C
C_____________________________________________________________________________
C
        SUBROUTINE EIGEN(A,R,N,MV)
C
C   MATRIX DIAGONALIZATION ROUTINE FOR SYMMETRIC MATRICES FROM THE
C   'SSP' LIBRARY
C
C   - Jacobi type algorithm
C   - matrix to be diagonalized has to be in symmetric storage mode
C     (A and R are presently dimensioned for 10x10 matrices)
C
        DIMENSION A(55),R(100)
      DOUBLE PRECISION A,R,ANORM,ANRMX,THR,X,Y,SINX,SINX2,COSX,
     1                 COSX2,SINCS,RANGE,DSQRT,DABS
C
C
C       GENERATE IDENTITY MATRIX
C
5       RANGE=1.0D-12
      IF(MV-1) 10,25,10
   10 IQ=-N
      DO 20 J=1,N
      IQ=IQ+N
      DO 20 I=1,N
      IJ=IQ+I
      R(IJ)=0.0d0
      IF(I-J) 20,15,20
   15 R(IJ)=1.0d0
   20 CONTINUE
C
C        COMPUTE INITIAL AND FINAL NORMS (ANORM AND ANORMX)
C
   25 ANORM=0.0d0
      DO 35 I=1,N
      DO 35 J=I,N
      IF(I-J) 30,35,30
   30 IA=I+(J*J-J)/2
      ANORM=ANORM+A(IA)*A(IA)
   35 CONTINUE
      IF(ANORM) 165,165,40
40      ANORM=1.414213562d0*DSQRT(ANORM)
      ANRMX=ANORM*RANGE/FLOAT(N)
C
C        INITIALIZE INDICATORS AND COMPUTE THRESHOLD, THR
C
      IND=0
      THR=ANORM
   45 THR=THR/FLOAT(N)
   50 L=1
   55 M=L+1
C
C        COMPUTE SIN AND COS
C
   60 MQ=(M*M-M)/2
      LQ=(L*L-L)/2
      LM=L+MQ
   62 IF(DABS(A(LM))-THR) 130,65,65
   65 IND=1
      LL=L+LQ
      MM=M+MQ
      X=0.5d0*(A(LL)-A(MM))
   68 Y=-A(LM)/DSQRT(A(LM)*A(LM)+X*X)
      IF(X) 70,75,75
   70 Y=-Y
   75 SINX=Y/DSQRT(2.0d0*(1.0d0+(DSQRT(1.0d0-Y*Y))))
      SINX2=SINX*SINX
   78 COSX=DSQRT(1.0d0-SINX2)
      COSX2=COSX*COSX
      SINCS =SINX*COSX
C
C        ROTATE L AND M COLUMNS
C
      ILQ=N*(L-1)
      IMQ=N*(M-1)
      DO 125 I=1,N
      IQ=(I*I-I)/2
      IF(I-L) 80,115,80
   80 IF(I-M) 85,115,90
   85 IM=I+MQ
      GO TO 95
   90 IM=M+IQ
   95 IF(I-L) 100,105,105
  100 IL=I+LQ
      GO TO 110
  105 IL=L+IQ
  110 X=A(IL)*COSX-A(IM)*SINX
      A(IM)=A(IL)*SINX+A(IM)*COSX
      A(IL)=X
  115 IF(MV-1) 120,125,120
  120 ILR=ILQ+I
      IMR=IMQ+I
      X=R(ILR)*COSX-R(IMR)*SINX
      R(IMR)=R(ILR)*SINX+R(IMR)*COSX
      R(ILR)=X
  125 CONTINUE
      X=2.0d0*A(LM)*SINCS
      Y=A(LL)*COSX2+A(MM)*SINX2-X
      X=A(LL)*SINX2+A(MM)*COSX2+X
      A(LM)=(A(LL)-A(MM))*SINCS+A(LM)*(COSX2-SINX2)
      A(LL)=Y
      A(MM)=X
C
C        TESTS FOR COMPLETION
C
C        TEST FOR M = LAST COLUMN
C
  130 IF(M-N) 135,140,135
  135 M=M+1
      GO TO 60
C
C        TEST FOR L = SECOND FROM LAST COLUMN
C
  140 IF(L-(N-1)) 145,150,145
  145 L=L+1
      GO TO 55
  150 IF(IND-1) 160,155,160
  155 IND=0
      GO TO 50
C
C        COMPARE THRESHOLD WITH FINAL NORM
C
  160 IF(THR-ANRMX) 165,165,45
C
C        SORT EIGENVALUES AND EIGENVECTORS
C
  165 IQ=-N
      DO 185 I=1,N
      IQ=IQ+N
      LL=I+(I*I-I)/2
      JQ=N*(I-2)
      DO 185 J=I,N
      JQ=JQ+N
      MM=J+(J*J-J)/2
      IF(A(LL)-A(MM)) 170,185,185
  170 X=A(LL)
      A(LL)=A(MM)
      A(MM)=X
      IF(MV-1) 175,185,175
  175 DO 180 K=1,N
      ILR=IQ+K
      IMR=JQ+K
      X=R(ILR)
      R(ILR)=R(IMR)
  180 R(IMR)=X
  185 CONTINUE
C
      RETURN
      END
C_____________________________________________________________________________
C_____________________________________________________________________________