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. 3.I.2010                             ----- 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  1.01.10:  Port to CVF
c  3.01.09:  Port to IVF
C
C----------------------------------------------------------------------------
C
C   NOTES:  This is an IBM-PC bound program and it is tied to the graphics
C           package incorporated with Compaq Visual Fortran
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     I N S T A L L A T I O N:
C-----------------------------------------------------------------------------
C
C     1/ place VECTOR.EXE and PMIFST.CFG in directory C:\ROT
C     2/ using Windows Explorer send a shortcut to VECTOR to the desktop
C     3/ add C:\ROT to the PATH if it is planned to launch from the command
C        prompt
C
C-----------------------------------------------------------------------------
C     R U N N I N G:
C-----------------------------------------------------------------------------
C
C   - If the C:\ROT directory has been added to the PATH then VECTOR can
C     be called from the Command prompt window set to any directory.  
C     This directory is treated as the default directory for any output.
C
C   - Drag and drop operation is possible if VECTOR icon has been placed
C     on the desktop.  Drag the appropriate VIBCA output file to VECTOR 
C     and then just press the ENTER key in response to the file name
C     question.  The output will go to the c:\ROT directory.
C
C-----------------------------------------------------------------------------
C    C O M P I L A T I O N:
C-----------------------------------------------------------------------------
C
C    This version will only compile satisfactorily with
C    Compaq Visual Fortran 6.50 (and possibly with not too distant earlier
C    versions of Microsoft Powerstation Fortran)
C    Note that there is compatibility with older versions in that
C
C    Compilation is now to be for QWIN graphics - this necessitates explicit
C    programming out of several unnecessary frills, but results in smoother
C    launch of the program than is possible with the STANDARD graphics as used 
C    previously.
C
C-------------------------------- 
C    Command line compilation:
C-------------------------------- 
C
C    Simplest compilation for the local machine:
C
C        df -static -libs=qwin -fpscomp:filesfromcmd pmifst.for
C
C    Optimised compilation for any PENTIUM:
C
C        df -check:all -nodebug -traceback -arch=pn1 -tune=pn1
C           -fast -static -libs=qwin -fpscomp:filesfromcmd pmifst.for
C
C    Other processor options are pn2,pn3,pn4,k6_2,k7
C
C-------------------------------- 
C    Visual Studio compilation:
C-------------------------------- 
C
C    FORTRAN:  /compile_only /fpscomp:filesfromcmd
C              /libs:qwin /nologo /nopdbfile /optimize:3 /traceback /tune:pn1 
C              /architecture:pn1 /static
C
C    LINK:     kernel32.lib /nologo /subsystem:windows /pdb:none 
C              /machine:IX86 /out:"Debug/pmifst.exe"
C
C
C-------------------------------------------------------- 
C    Visual Studio Net compilation for Intel Visual Fortran:
C--------------------------------------------------------
C
C    FORTRAN:  /nologo /G5 /fpscomp:filesfromcmd /warn:unused 
C              /Qsave /module:"Debug/" /object:"Debug/" /traceback 
C              /libs:qwin /c
C
C    LINK:     /OUT:"Debug/vector.exe" /NOLOGO /SUBSYSTEM:WINDOWS 
C              /MACHINE:IX86 kernel32.lib
C
C_____________________________________________________________________________
C
C   Standard colours:
C
C     0 - black         4 - dark red        8 - dark grey    12 - red
C     1 - dark blue     5 - dark purple     9 - blue         13 - purple
C     2 - dirty green   6 - dirty yellow   10 - green        14 - yellow
C     3 - blue/green    7 - grey           11 - light blue   15 - white
C
C===========================================================================
c
C...Initialization commands for graphics.  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
      USE DFLIB
c
      PARAMETER (ntextc=0, ntextb=7, nplotc=15, nplotb=1)
      RECORD /rccoord/curpos
      RECORD /xycoord/ixy
      RECORD /wxycoord/wxy
      logical*2 true
C
      PARAMETER (maxat=25, maxvec=3*maxat, maxpts=maxat+maxvec,
     *           maxbon=6*maxpts,true=.true.)
c
      COMMON /BLK0/ATOMS,VBA,TRANS,PRC,RJA,VCA,INITIA,NATOMS
     1 /BLK1/CORD,CORPLT,AXEND,axendr,ASPECT,ROTMAT,NATMS
     1 /BLK2/W,IBOND,NBONDS
      COMMON /PARBLK/RANGEA,RANGEB,AMIN,BMIN
      COMMON /limits/wxy,maxx,maxy,LINOFS,curpos,ixy,
     * 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,AMAX,BMAX,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
      CON=0.174532925199D-1
      ITSAVE=1
      linofs=14
C
C...start graphics
C
      call startg(iconf)
      dummy4=passdirkeysqq(.true.)
c
C...HEADER
C
      numfonts = INITIALIZEFONTS ( )
      fontnum = SETFONT ('t''Arial''h75w25ei')
      dummy4=setbkcolor(ntextb)
      call clearscreen($gclearscreen)
c
      NBOTL=110
      dummy=setcolor(int2(15))
      CALL MOVETO (INT2( 0), INT2(NBOTL), ixy)
      dummy=lineto(INT2( maxx), INT2(NBOTL))
      CALL MOVETO (INT2( 0), INT2(0), ixy)
      dummy=lineto(INT2( maxx), INT2(0))
      dummy=setcolor(int2(8))
      dummy=floodfill(int2(1),int2(1),int2(15))
c
      nvert=(NBOTL-75)/2
      nhor=(maxx-540)/2
      dummy=setcolor(int2(11))
      CALL MOVETO (INT2(nhor), INT2(nvert), ixy)
      CALL OUTGTEXT('VECTOR')
      dummy=setcolor(int2(9))
      CALL MOVETO (INT2(nhor+1), INT2(nvert+1), ixy)
      CALL OUTGTEXT('VECTOR')
      dummy=setcolor(int2(1))
      CALL MOVETO (INT2(nhor+2), INT2(nvert+2), ixy)
      CALL OUTGTEXT('VECTOR')
c
      nvert=nbotl-40
      nhor=nhor+220

      fontnum = SETFONT ('t''Arial''h24w8e')
      dummy=setcolor(int2(7))
      CALL MOVETO (INT2(nhor),  INT2(nvert-1), ixy)
      CALL OUTGTEXT('Normal Coordinate Displacement'//
     *              ' Vectors')
      dummy=setcolor(int2(0))
      CALL MOVETO (INT2(nhor+1),  INT2(nvert), ixy)
      CALL OUTGTEXT('Normal Coordinate Displacement'//
     *              ' Vectors')
C
C
      dummy=setcolor(int2(15))
      CALL MOVETO (INT2(    0), INT2(NBOTL+ 32), ixy)
      dummy=lineto(INT2( maxx), INT2(NBOTL+ 32))
      dummy=setcolor(int2(8))
      dummy=floodfill(int2(1),INT2(NBOTL+30),int2(15))
c
      dummy=setcolor(int2(0))
      CALL MOVETO (INT2(    0), INT2(NBOTL+ 32), ixy)
      dummy=lineto(INT2( maxx), INT2(NBOTL+ 32))
      dummy=setcolor(int2(7))
      CALL MOVETO (INT2(    0), INT2(NBOTL+  1), ixy)
      dummy=lineto(INT2( maxx), INT2(NBOTL+  1))    
c
      fontnum = SETFONT ('t''Arial''h20w10')
      dummy=setcolor(int2( 0))
      CALL MOVETO (INT2( 11), INT2(NBOTL+ 7), ixy)
      CALL OUTGTEXT('version 3.I.2010')
      CALL MOVETO (INT2( maxx-169), INT2(NBOTL+ 7), ixy)
      CALL OUTGTEXT('Zbigniew KISIEL')
      dummy=setcolor(int2(15))
      CALL MOVETO (INT2( 10), INT2(NBOTL+ 6), ixy)
      CALL OUTGTEXT('version 3.I.2010')
      CALL MOVETO (INT2( maxx-170), INT2(NBOTL+ 6), ixy)
      CALL OUTGTEXT('Zbigniew KISIEL')
c
C...Warn of missing config file
C
70    nrlin=nint(real(NBOTL)/(real(maxy)/real(myrows)))+5
      call settextposition(nrlin,1,curpos)
      if(iconf.eq.0)then
        dummy=setcolor(int2(12))
        fontnum = SETFONT ('t''Arial''h18w9e')
        CALL MOVETO (INT2(100), INT2(200), ixy)
        CALL OUTGTEXT(
     *    'Configuration file C:\ROT\PMIFST.CFG was not found:')
        CALL MOVETO (INT2(250), INT2(220), ixy)
        fontnum = SETFONT ('t''Arial''h18w9i')
        CALL OUTGTEXT(
     *    'default sized window of 800x540 pixels will be used')
        call settextposition(15,1,curpos)
      endif
c
      WRITE(emplin,'(79(1H ))')
1     format(i5)
c
      dummy=settextcolor(int2(ntextc))
      dummy4=setbkcolor(ntextb)
C
C   Input of data
C
3003  FORMAT(I5)
C
6001  WRITE(*,3040)
3040  FORMAT(1X/' Name of VIBCA results file:  ',$)
      dummy=displaycursor($gcursoron)
      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
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
C
      dummy4=setbkcolor(ntextb)
      dummy=displaycursor($gcursoroff)
      call clearscreen($gclearscreen)
      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
      if(myrows.eq.25.or.myrows.eq.30)then
        mult=1
      else
        mult=2
      endif
      ASPECT=1.d0
c     ASPECT=maxy/(0.75*maxx)
      CALL setlinestyle(int2(#FFFF))
      IH=maxx/2
      IV=maxy/2
c
c...bounding box
      dummy4=setbkcolor(nplotb)
      CALL setviewport(int2(1),int2(linofs+1),
     *                 int2(maxx-1),int2(maxy-mult*linofs-1))
      dummy=setwindow(TRUE,-1.d0,-1.d0,1.d0,1.d0)
      dummy=setcolor(int2(15))
      CALL  moveto_w(-1.d0,-1.d0,wxy)
      dummy=lineto_w(-1.d0, 1.d0)
      dummy=lineto_w( 1.d0, 1.d0)
      dummy=setcolor(int2(0))
      dummy=lineto_w( 1.d0,-1.d0)
      dummy=lineto_w(-1.d0,-1.d0)
c
      CALL setviewport(int2(2),int2(linofs+2),
     *                 int2(maxx-2),int2(maxy-mult*linofs-2))
      dummy=setwritemode($GPSET)
      CALL clearscreen($GVIEWPORT)
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(int2(0),int2(linofs+1),
     *                 int2(maxx-IH),int2(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(int2(maxx-IH),int2(linofs+1),
     *                 maxx,int2(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(int2(IHH),int2(maxy-linofs-1-IV),
     *                 int2(IHH+IH),int2(maxy-linofs-1))
      dummy=setwindow(true,amin,bmin,amax,bmax)
      CALL view(1,3,1,0)
c
      dummy=settextcolor(int2(11))
      CALL settextposition(3,2,curpos)
      CALL outtext('A-B')
      CALL settextposition(3,int2(mycols/2),curpos)
      CALL outtext('B-C')
      CALL settextposition(int2(myrows-2),2,curpos)
      CALL outtext('A-C')
C
500   IK=INKEY(N)
      IF(IK.NE.13)GOTO 500
      dummy=displaycursor($gcursoron)
c
      dummy4=setbkcolor(ntextb)
      dummy=settextcolor(int2(ntextc))
      call clearscreen($GCLEARSCREEN)
C
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
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)then
        dummy=setexitqq(qwin$exitnopersist)
        STOP
      endif
      dummy=displaycursor($gcursoroff)
      CALL ROTATE
      dummy=setexitqq(qwin$exitnopersist)
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/CORD,CORPLT,AXEND,axendr,ASPECT,ROTMAT(3,3),NATMS
     *       /BLK2/WT,IBOND,NBONDS
      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/CORD,CORPLT,AXEND,axendr,ASPECT,ROTMAT(3,3),NATMS
      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 startg(iconf)
c
C   This routine uses QWIN graphics and techniques from the CLEANWIN programming
C   example for CVF6 to avoid the full-screen startup of standard graphics,
C   while preserving a simple frame.  
C   Note the use of the WIN32 routines MoveWindow, UpdateWindow, GetWindowLong,
C   SetWindowLong, GetHWndQQ - their operation and parameter values are not
C   really understood!
c
      USE DFLIB
      USE DFWIN
c
      RECORD /rccoord/curpos
      RECORD /xycoord/ixy
      RECORD /wxycoord/wxy
      INTEGER*2 maxx,maxy,linofs,mymode,myrows,mycols,ifc,ifm,ifhb,ifb,
     *          idelta,ivista
      integer*4 dummy4,i
      character fntnam*30,line*80
      COMMON /limits/wxy,maxx,maxy,linofs,curpos,ixy,
     *               mymode,myrows,mycols
      type (windowconfig)wc
      type (qwinfo)win
      logical status
C
C
C...Determine the type of Windows (using a WIN32 function).  Official Windows 
C   version numbers are (os.dwMajorVersion,os.dwMinorversion,...):  
C   
C   Server 2008     6.0.6001            ME      4.90
C   Vista           6.0.6000            NT 4    4.0.1381
C   Server 2003     5.2                 98      4.10
C   XP              5.1.2600            95      4.00
C   2000            5.0.2195   
c      
      type (t_OSVERSIONINFO)os    
      os.dwOSVersionInfoSize=sizeof(os)
      status=GetVersionEx( os )
c     write(*,*)os.dwPlatformId,os.dwMajorVersion,os.dwMinorversion
c     ntsys=os.dwPlatformId-1
C
C...An alternative way to determine the type of Windows, but this has the 
C   the disadvantage of a flashing COMMAND PROMPT screen
C
c     status=SYSTEMQQ('ver>opsys.ver')
c     open(3,file='opsys.ver',status='unknown')
c     read(3,'(A)')line
c     read(3,'(A)')line
c     if(line(1:7).eq.'Windows')then
c       ntsys=0
c     else
c       ntsys=1
c       if(line(1:20).eq.'Microsoft Windows XP')ntsys=2
c     endif
c     close(3)
c     status=SYSTEMQQ('del opsys.ver')
c
c...set the principal window parameters, as hardcoded below and
c   specified in the PMIFST.CFG file
c
      wc.numtextcols=80
      wc.numtextrows=30
c
      open(3,file='c:\rot\pmifst.cfg',status='old',err=12)
7       read(3,'(a)')line
        if(line(1:1).eq.'!')goto 7
        read(line,5)wc.numxpixels
        read(3,'(a)')line
        read(line,5)wc.numypixels
        read(3,'(a)')line
        fntnam=line(36:65)
5       format(35x,i4)
c
        wc.fontsize=QWIN$EXTENDFONT
        wc.numcolors=-1
        wc.extendfontname=trim(fntnam)//char(0)
        wc.extendfontsize=-1
c
        wc.extendfontattributes=0
8       read(3,'(a)')line
        if(line(1:1).eq.'!')goto 9
        read(line,5)iattr
        if(iattr.lt.1.or.iattr.gt.15)then
           write(*,10)iattr
10         format(1x//' Extended font attribute from PMIFST.CFG is',i5,
     *                ', which is illegal (1-15 allowed)'//
     *                ' **** TRY AGAIN! *****'//)
           pause
           stop
        endif
      if(iattr.eq. 1)wc.extendfontattributes=wc.extendfontattributes
     *               +QWIN$EXTENDFONT_NORMAL
      if(iattr.eq. 2)wc.extendfontattributes=wc.extendfontattributes
     *               +QWIN$EXTENDFONT_UNDERLINE
      if(iattr.eq. 3)wc.extendfontattributes=wc.extendfontattributes
     *               +QWIN$EXTENDFONT_BOLD
      if(iattr.eq. 4)wc.extendfontattributes=wc.extendfontattributes
     *               +QWIN$EXTENDFONT_ITALIC
c
      if(iattr.eq. 5)wc.extendfontattributes=wc.extendfontattributes
     *               +QWIN$EXTENDFONT_FIXED_PITCH
      if(iattr.eq. 6)wc.extendfontattributes=wc.extendfontattributes
     *               +QWIN$EXTENDFONT_VARIABLE_PITCH
c
      if(iattr.eq. 7)wc.extendfontattributes=wc.extendfontattributes
     *               +QWIN$EXTENDFONT_FF_ROMAN
      if(iattr.eq. 8)wc.extendfontattributes=wc.extendfontattributes
     *               +QWIN$EXTENDFONT_FF_SWISS
      if(iattr.eq. 9)wc.extendfontattributes=wc.extendfontattributes
     *               +QWIN$EXTENDFONT_FF_MODERN
      if(iattr.eq.10)wc.extendfontattributes=wc.extendfontattributes
     *               +QWIN$EXTENDFONT_FF_SCRIPT
      if(iattr.eq.11)wc.extendfontattributes=wc.extendfontattributes
     *               +QWIN$EXTENDFONT_FF_DECORATIVE
c
      if(iattr.eq.12)wc.extendfontattributes=wc.extendfontattributes
     *               +QWIN$EXTENDFONT_ANSI_CHARSET
      if(iattr.eq.13)wc.extendfontattributes=wc.extendfontattributes
     *               +QWIN$EXTENDFONT_DEFAULT_CHARSET
      if(iattr.eq.14)wc.extendfontattributes=wc.extendfontattributes
     *               +QWIN$EXTENDFONT_SYMBOL_CHARSET
      if(iattr.eq.15)wc.extendfontattributes=wc.extendfontattributes
     *               +QWIN$EXTENDFONT_OEM_CHARSET
      goto 8
c
9     close(3)
      iconf=1
      goto 11
C
C...open default sized window for a 800x600 screen if the configuration file
C   cannot be opened
C
12    wc.numxpixels=800
      wc.numypixels=540
      wc.extendfontname='Courier New'
      iconf=0
C
C...Kill menu
C
11    DO I = 7,1,-1
      STATUS= DELETEMENUQQ(I, 0)
      END DO
C
C...Kill status bar
C
      i = clickqq( QWIN$STATUS )
C
C...Kill scroll bars and unwanted features (note that the title seems possible
C   only on the (killed) daughter window and not on the framewindow
C
      i = GetWindowLong( GetHWndQQ(QWIN$FRAMEWINDOW), GWL_STYLE )
      i = ior( iand( i, not(WS_THICKFRAME) ), WS_BORDER )
      i = iand( i, not(WS_MAXIMIZEBOX) )
      k = SetWindowLong( GetHWndQQ(QWIN$FRAMEWINDOW), GWL_STYLE, i )    
C
      i = GetWindowLong( GetHWndQQ(0), GWL_STYLE )
      i = ior(iand( i, not(WS_CAPTION.or.WS_SYSMENU.or.WS_THICKFRAME)), &
     &    WS_BORDER)
      k = SetWindowLong( GetHWndQQ(0), GWL_STYLE, i )     
c
c...Position window - for compatibility with small pixel size screens make the
c   top and left edge of the bounding frame disappear
c   
      ifxed= GetSystemMetrics(sm_cxfixedframe)                          x thickness of frame around nonzizable window
      ifyed= GetSystemMetrics(sm_cyfixedframe)                          y thickness of frame around nonzizable window
c     win.x = -ifxed
c     win.y = -ifyed
      win.x =  ifxed
      win.y =  ifyed
c
c...Correct sizing parameters to take account of removal of the menu bar and of 
c   the caption area for the child window - these vary with system font size and
c   are augmented with emprirically established IDELTA fudge parameter as well
c   as IVISTA parameter for VISTA
C   
C
C  W98,ME small fonts:        ifc=19 ifm=19 idelta=5 -> sum=43
C  W2000, small=100% fonts:   ifc=19 ifm=19 idelta=5 -> sum=43
C  W2000, large=125% fonts:   ifc=24 ifm=24 idelta=4 -> sum=52
C  W2000, custom=132% fonts:  ifc=25 ifm=25 idelta=3 -> sum=53
C  WindowsXP,small fonts:     ifc=26 ifm=20 idelta=2 -> sum=48
C  WindowsXP, 125% fonts:     ifc=32 ifm=25 idelta=2                    32,25,1,1,3,3
C  Vista, 125% fonts:         ifc=25 ifm=25                             25,25,1,1,3,3
C 
C
      ifc  = GetSystemMetrics(sm_cycaption)                             ht of caption area
      ifm  = GetSystemMetrics(sm_cymenu)                                ht of single line menu bar
      ifhb = GetSystemMetrics(sm_cxborder)                              width of border
      ifb  = GetSystemMetrics(sm_cyborder)                              height of border
c     write(*,*)ifc,ifm,ifb,ifhb,ifxed,ifyed
c     pause
c      
      if(os.dwMajorVersion.eq.6)then
        ivista=10
      else
        ivista=0
      endif
c
      idelta=4
      if(ifc.eq.19.and.ifc.eq.19)idelta=5                                             
      if(ifc.eq.25.and.ifc.eq.25)idelta=3                                             
      if(ifc.ge.26)idelta=2                                             
c           
      win.w = wc.numxpixels-(2*ifhb+ivista)
      win.h = wc.numypixels-(ifc+ifm+idelta+ivista)
c
      win.type=qwin$set
      dummy4 = setwsizeqq(qwin$framewindow,win)
      status = getwsizeqq(QWIN$FRAMEWINDOW,QWIN$SIZECURR, win)
c
      wc.numtextcols=80
      wc.numtextrows=30
      wc.title=' 'C
c 
      status=setwindowconfig(wc)
      if(.not.status)status=setwindowconfig(wc)
C
C...Magical Windows incantations to make style set above real (without
C   these commands the active window does not expand to the size of the 
C   program framewindow)
C     
      i = MoveWindow( GetHWndQQ(0), -1, -1, 0, 0, int4(.TRUE.) ) 
      call clearscreen($GCLEARSCREEN)
      status = UpdateWindow(GETHANDLEFRAMEQQ())    
c
C  pixel limits on x and y axes (0,maxx), (0,maxy)
c
      maxx=wc.numxpixels-1
      maxy=wc.numypixels-1
      myrows=wc.numtextrows
      mycols=wc.numtextcols
      linofs=nint(real(maxy)/real(myrows))+1
c
c     write(*,*)maxx,maxy
c     write(*,*)myrows,mycols
c     pause
c
      return
      end
C
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
      USE DFLIB
C
      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/wxy,maxx,maxy,LINOFS,curpos,ixy,
     * mymode,myrows,mycols,blue,red
      COMMON /BLK1/CORD,CORPLT,AXEND,axendr,ASPECT,ROTMAT,NATMS
     *       /BLK2/WT,IBOND,NBONDS
      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(int2(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(int2(11))
      DO 4 N=1,5,2
        CALL setlinestyle(int2(#8888))
        call moveto_w(AXEND(N,IX),AXEND(N,IY),wxy)
        dummy=lineto_w(AXEND(N+1,IX),AXEND(N+1,IY))
4     CONTINUE
      CALL setlinestyle(int2(#FFFF))
      dummy=setcolor(int2(15))
      IF(IAXIS.EQ.2)THEN
        dummy=setcolor(int2(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(int2(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(int2(12))
        ELSE
          DUMMY=SETCOLOR(int2(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(int2(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
      CALL setlinestyle(int2(#FFFF))
      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(int2(0))
        ELSE
          dummy=setcolor(int2(11))
        ENDIF
        if(maxy.lt.700)then
          CALL PLOTCH(1,1,ISTART,JSTART,1)
        else
          istart=istart-5
          jstart=jstart+3
          CALL PLOTCH(2,2,ISTART,JSTART,1)
        endif
51    CONTINUE
      dummy=setcolor(int2(15))
C
50    RETURN
      END
C
C_____________________________________________________________________________
C
      SUBROUTINE ROTATE
C
      USE DFLIB
C      
      RECORD /rccoord/curpos
      RECORD /xycoord/ixy
      RECORD /wxycoord/wxy
      logical*2 true
c
      PARAMETER (maxat=25, maxvec=3*maxat, maxpts=maxat+maxvec,
     *           maxbon=6*maxpts,true=.true.)
      PARAMETER (ntextc=0, ntextb=7, nplotc=15, nplotb=1)
      PARAMETER (nhelpc=0, nhelpb=7)
C
      COMMON /limits/wxy,maxx,maxy,LINOFS,curpos,ixy,
     * mymode,myrows,mycols,blue,red
      COMMON /BLK1/CORD,CORPLT,AXEND,axendr,ASPECT,ROTMAT,NATMS
     *       /BLK2/WT,IBOND,NBONDS
     *     /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
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
C
C...Complete plot refresh from here:
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   dummy=SETTEXTCOLOR(int2(ntextc))
      dummy4= setbkcolor(ntextb)
      CALL clearscreen($GCLEARSCREEN)
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=1.d0
c     ASPECT=maxy/(0.75*maxx)
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
      CALL setlinestyle(int2(#FFFF))
C
C...Write top and bottom of display
C
      if(myrows.eq.25.or.myrows.eq.30)then
        mult=1
      else
        mult=2
      endif
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(18),
     * CHAR(27),CHAR(26)
600   FORMAT('Mode',I3,':',F7.1,'cm-1, ',A10,13X,A1,'=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
C...bounding box
C
      dummy4=setbkcolor(nplotb)
      CALL setviewport(int2(1),int2(linofs+1),
     *                 int2(maxx-1),int2(maxy-mult*linofs-1))
      dummy=setwindow(TRUE,-1.d0,-1.d0,1.d0,1.d0)
      dummy=setcolor(int2(15))
      CALL  moveto_w(-1.d0,-1.d0,wxy)
      dummy=lineto_w(-1.d0, 1.d0)
      dummy=lineto_w( 1.d0, 1.d0)
      dummy=setcolor(int2(0))
      dummy=lineto_w( 1.d0,-1.d0)
      dummy=lineto_w(-1.d0,-1.d0)
      dummy=setcolor(int2(nplotc))
c
      CALL setviewport(int2(2),int2(linofs+2),
     *                 int2(maxx-2),int2(maxy-mult*linofs-2))
      dummy=setwindow(TRUE,amin,bmin,amax,bmax)
      dummy=setwritemode($GPSET)
      CALL clearscreen($GVIEWPORT)
C
C...Plot the molecule
C
      CALL VIEW(IA,IB,ICONTR,1)                                         <-----
C
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
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)
c
      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.'U')GOTO 800
      IF(KK.EQ.'H')GOTO 500
c
c...toggle atom numbering (with N)
c
      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
c
c...toggle inertial axes (with A)
c
      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
c
c...zero previous rotations (with Z)
c
      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 (with O)
c
      IF(KK.EQ.'O')THEN
        CALL clearscreen($GVIEWPORT)
        DUMMY= SETTEXTCOLOR(int2( 7))
        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)
c
        write(*,'(1x/)')
        WRITE(*,502)
        CALL setviewport(int2(0),int2(linofs+1),
     *                   int2(maxx-1),int2(maxy-linofs-1))
        CALL moveto(0,0,ixy)
        dummy=lineto(int2(0),int2(maxy-2*linofs))
        CALL setviewport(int2(1),int2(linofs+1),
     *                   int2(maxx-1),int2(maxy-linofs-1))
453     IK=INKEY(N)
        IF(IK.EQ.0)GOTO 453
        GOTO 3
      ENDIF
C
C...Parameter check (with P)
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 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
C
C...Exit
C
20    DUMMY= SETTEXTCOLOR(int2(ntextc))
      dummy4 = setbkcolor(ntextb)
      CALL settextposition(1,1,curpos)
      CALL outtext(emplin)
      DUMMY=SETTEXTCOLOR(int2(15))
      dummy4=setbkcolor (12)
      WRITE(OUTSTR,'(A)')'ARE YOU SURE ?'
      CALL settextposition(1,1,curpos)
      CALL outtext(outstr)
916   IK=INKEY(N)
      KK=CHAR(IK)
c
      IF(KK.EQ.'Y'.OR.KK.EQ.'y')then
        dummy=settextcolor(int2(ntextc))
        dummy4=setbkcolor(ntextb)
        call clearscreen($gclearscreen)
        GOTO 4
      endif
C
      IF(KK.EQ.'N'.OR.KK.EQ.'n')then
        DUMMY= SETTEXTCOLOR(int2(ntextc))
        dummy4 = setbkcolor(ntextb)
        CALL settextposition(1,1,curpos)
        CALL outtext(toplin)
        GOTO 2
      endif
      GOTO 916
C
C
C...Help screen (with H)
C
500   DUMMY=SETTEXTCOLOR(int2(nhelpc))
      dummy4= setbkcolor(nhelpb)
      CALL clearscreen($GCLEARSCREEN)
      CALL settextposition(1,1,curpos)
      WRITE(*,232)
232   FORMAT(1x/'___________A V A I L A B L E   O P T I O N S',
     * 36(1h_))
      write(*,234)
234   FORMAT(1x/
     *'            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
501   IK=INKEY(N)
      IF(IK.EQ.0)GOTO 501
      GOTO 333
C
C...Rotation about the horizontal axis (with E)(with X)
C
5     IF(IK.EQ.69)ANGLE=-ANGLE
      CALL TRANSF(ANGLE,IB,IC,0)                                        <----- 
      GOTO 3
C
C...Rotation about the vertical axis (with E)(with X)
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 (with W)(with R)
C
7     IF(IK.EQ.87)ANGLE=-ANGLE
      CALL TRANSF(ANGLE,IA,IB,0)                                        <----- 
      GOTO 3
C
C...Toggle weighted/unweighted eigenvectors (with U)
C
800   IF(IWTUNW.EQ.1)THEN
        IWTUNW=2
      ELSE
        IWTUNW=1
      ENDIF
      GOTO 901
C
C...Toggle through eigenvectors (with up arrow) (with down arrow)
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)                                        <-----
      dummy=SETTEXTCOLOR(int2(ntextc))
      dummy4= setbkcolor(ntextb)
      GOTO 33
C
C...Reverse phase (with -)
C
700   IPHASE=-IPHASE
      GOTO 901
C
C...Change vector scaling (with left arrow) (with right arrow)
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/CORD,CORPLT,AXEND,axendr,ASPECT,ROTMAT(3,3),NATMS
      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
      USE DFLIB
C
      PARAMETER (ntextc=0, ntextb=7, nplotc=15, nplotb=1)
      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/wxy,maxx,maxy,LINOFS,curpos,ixy,
     * mymode,myrows,mycols,blue,red
      COMMON /BLK1/CORD,CORPLT,AXEND,axendr,ASPECT,ROTMAT,NATMS
      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 ))')
      dummy4=setbkcolor(ntextb)
      dummy=settextcolor(int2(ntextc))
      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
      USE DFLIB
C
      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/wxy,maxx,maxy,LINOFS,curpos,ixy,
     * mymode,myrows,mycols,blue,red
      COMMON /BLK1/COR,CORPLT,AXEND,axendr,ASPECT,ROTMAT,NATMS
      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
      USE DFLIB
C
      RECORD /rccoord/curpos
      RECORD /xycoord/ixy
      RECORD /wxycoord/wxy
C
      COMMON /limits/wxy,maxx,maxy,LINOFS,curpos,ixy,
     * 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),NOCHAR,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(NOCHAR,'(A1)')CHAR(0)
      SPACE=' '
      DO 100 NCHARS=80,1,-1
100     IF(TEXT(NCHARS).NE.NOCHAR.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 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_____________________________________________________________________________
c
      integer*2 function INKEY(N2)
c
c   By L Pszczolkowski:
c
c...This emulates    for MSF PS1.0 the INKEY function of Z.Czumaj which
c   in turn emulated for MSF5.0    the INKEY function from IIUWGRAF graphics
c   library for the Hercules card
c
c   The function GETCHARQQ returns the ASCII character if the corresponding
c   key was pressed.  If function or direction key was pressed then 0
c   or hex E0 is returned and another call to GETCHARQQ is required to
c   get the extended code of the character
c
      USE DFLIB
c
      INTEGER*2 IK,n2
      CHARACTER*1 KK
c
      KK=GETCHARQQ()
      IK=ICHAR(KK)
      IF(IK.EQ.0 .OR. IK.EQ.224 ) THEN
          KK=GETCHARQQ()
          IK=-ICHAR(KK)
      ENDIF
      n2=ik
      INKEY=IK
      END
C_____________________________________________________________________________
C_____________________________________________________________________________
C
C    Routines used by VECTOR:
C
C    BONLOC        - locate bonds
C    SCALIN        - axis scaling
C    graphicsmode  - set graphics
C    VIEW          - projection on one of the principal planes
C    ROTATE        - main graphics screen
C    TRANSF        - carry out the rotation transformation
C    PARCHE        - parameter checks
C    BONDL         - evaluate bond length (used by PARCHE)
C    ANGLE         - evaluate bond angle (used by PARCHE)
C    DANGLE        - evaluate dihedral angle (used by PARCHE)
C    INNUM         - input and checks of validity of two digit number
C    PLOTCH        - plotter type character generator
C    MOMINA        - calculate principal moments of inertia
C    EIGEN         - matrix diagonalisation
C    INKEY         - keyboard input
C_____________________________________________________________________________
C_____________________________________________________________________________