program mind
c-----------------------------------------------------------------------
c     Interactively displays molecules on the screen.
c
c   Written by :
c 
c           Julian Tirado-Rives	(julian@doctor.chem.yale.edu)
c           James Blake 	(jim@doctor.chem.yale.edu)
c           Department of Chemistry
c           Yale University
c
c Permission to use, copy, and distribute this software and its
c documentation for any purpose without fee is hereby granted, 
c provided that the above notice appear in all copies.
c
c Revision History
c
c     3.0     added h-bond determination and plotting.
c     3.1     added "set display {normal|reversed|clear}" toggle for 
c             normal display, clear atoms/filled bonds, and clear atoms/
c             clear bonds.
c-----------------------------------------------------------------------
      parameter (natom = 10000)
      parameter (nbond = 50000)
c
      integer start,end
      logical wireframe,ball_and_stick,need_sort,label_num,
     &        label_cha,first_draw,need_lim,perspective,
     &        autoscl,showhb,havehb
      logical more_data,round_bonds
      logical color_ws,mono_ws,PostScript,save_color,save_mono
      character*2 command
      character*80 line,filename,cmd_line,plotfile
c
      common/j_f_b/more_data
      common/const/pi,twopi,degtorad
      common/draw_log/wireframe,ball_and_stick,need_sort,label_num,
     &                label_cha,need_lim,perspective,round_bonds,
     &                autoscl,showhb,havehb
      common/mol_num/numat,num(natom),co(3,natom),atrad(natom)
      common/font_num/ifont
      common /device_type/ color_ws,mono_ws,PostScript
      common/commands/num_cmd,cmd_line(100)
c
      autoscl = .true.
      pi = acos(-1.0e0)
      twopi=2.0e0*pi
      degtorad=pi/180.0e0
c
c     Initializing the screen of the terminal
c
      call plots(0,1)
c
c     Beginning of loop :
c
      do while (.true.)

         if (first_draw) then
            command='DR'
            first_draw=.false.
         else
            read (*,'(A)',err=999,end=999) line
            call parse_command(line,command,angx,angy,angz)
         endif
c
         if (command.eq.'QU') then
            goto 999
         else if (command.eq.'NE') then
            close(9)
            read(*,'(A)') filename
            num_cmd = 0
            call init_mind(.true.,filename,command)
            if (num_cmd.gt.0) call do_macro(command,num_cmd,cmd_line)
         else if (command.eq.'MO') then
            if (more_data) then
               num_cmd = 0
               call init_mind(.false.,filename,command)
            if (num_cmd.gt.0) call do_macro(command,num_cmd,cmd_line)
               command='DR'
            else
               command='  '
               call notify('No additional data in this file. ')
            endif
         else if (command.eq.'TU') then
            call rot_mol(co,numat,angx,angy,angz,.true.)
            command='DR'
C
C        Write out a plot file and send it to the LaserWriter.
C
         else if (command.eq.'PL') then
            save_color = color_ws
            save_mono  = mono_ws
            color_ws   = .false. 
            mono_ws    = .false. 
            PostScript = .true.
            if (autoscl) call auto_scale
            call plots(1,2)
            call newpen(1)
            call fontstyle(ifont)
            call master_draw
            call notify('Generating plot ... ')
            call plot(0.0,0.0,999)
            color_ws = save_color
            mono_ws = save_mono 
            PostScript = .false.
            call notify('Plotting on the LaserWriter... ')
            call print_it
            call notify('Plot finished. ')
            command='  '
C
C        Write out a plot file and MoVe it to the named file.
C
         else if (command.eq.'FI') then
            read(*,'(A)') plotfile
            call limits(plotfile,start,end)
            ipl_len = end - start + 1
            save_color = color_ws
            save_mono  = mono_ws
            color_ws   = .false.
            mono_ws    = .false.
            PostScript = .true.
            if (autoscl) call auto_scale
            call plots(1,2)
            call newpen(1)
            call fontstyle(ifont)
            call master_draw
            call notify('Generating plot ... ')
            call plot(0.0,0.0,999)
            color_ws = save_color
            mono_ws = save_mono
            PostScript = .false.
            call notify('Moving file ... ')
            call move_it(plotfile,ipl_len)
            call notify('Finished. ')
            command='  '
c
c        Make a stereo plot by writing out two different plot files
c        "left" and "right" with rotations about the y-axis.
c
         else if (command.eq.'SV'.or.command.eq.'SP') then
            plotfile = 'left'
            ipl_len = 4
            save_color = color_ws
            save_mono  = mono_ws
            color_ws   = .false.
            mono_ws    = .false.
            PostScript = .true.
            call notify('Generating left plot ... ')
            call rot_mol(co,numat,0.0,-4.0,0.0,.true.)
            call plots(1,2)
            call newpen(1)
            call fontstyle(ifont)
            call master_draw
            call plot(0.0,0.0,999)
            call notify('Moving file ... ')
            call move_it(plotfile,ipl_len)
c
            call notify('Generating right plot ... ')
            plotfile = 'right'
            ipl_len = 5
            call rot_mol(co,numat,0.0,8.0,0.0,.true.)
            call plots(1,2)
            call newpen(1)
            call fontstyle(ifont)
            call master_draw
            call plot(0.0,0.0,999)
            call notify('Moving file ... ')
            call move_it(plotfile,ipl_len)
c
            color_ws = save_color
            mono_ws = save_mono
            PostScript = .false.
            call rot_mol(co,numat,0.0,-4.0,0.0,.true.)
            call notify('Sending files to mugs ... ')
            if(command .eq.'SV') then
               call stereo_view
            else if(command .eq.'SP') then
               call stereo_pub
            endif
            call notify('Stereo plot generated. ')
            command='  '
         else if (command.eq.'CR') then
            ang_inc = 2.0
            angle = ang_inc
            do while (angle.le.360.0)
               call rot_mol(co,numat,0.0,ang_inc,0.0,.true.)
               call plots(1,1)
               call master_draw
               angle = angle + ang_inc
            enddo
            angle = 360.0 - (angle - ang_inc)
            call rot_mol(co,numat,0.0,angle,0.0,.true.)
            call plots(1,1)
            call master_draw
            command='  '
         endif
         if (command.eq.'DR') then
c
c           Set the scale :
c
            if (autoscl) call auto_scale
c
c           Clearing the screen :
c  
            call plots(1,1)
c
c           selecting pen 1 (The color black)
c
            call newpen(1)
c
c           Branch to master drawing subroutine :
c
            call master_draw
c
c           Loop back to dialog :
c
            command='  '
         endif
      enddo
c
c     Termination :
c
  999 continue 
      call plot(0.0,0.0,666)
      call exit(0)
      end
c----------------------------------------------------------------------
      subroutine readin
c----------------------------------------------------------------------
c
c          This subroutine does all the initial reading from the input
c     file.   It makes extensive use of the free-format reading library
c     previously developed.
c
      parameter (natom = 10000)
      parameter (nbond = 50000)
c
      real rea,min_z,max_z
      integer nint,nrea,ncha,inva,start,end
      character*80 line,form,cmd_line
      character*30 cha
      character*4 atlabel, at_sym
      logical continue
      logical show_atom
      logical wireframe,ball_and_stick,need_sort,label_num,label_cha,
     &        need_lim,perspective,more_data,round_bonds,
     &        autoscl,showhb,havehb
c
      dimension inva(40),rea(40),cha(40)
c
      common/j_f_b/more_data
      common/draw_spec/rad_fact,bond_fact,top_z,bottom_z
      common/draw_log/wireframe,ball_and_stick,need_sort,label_num,
     &                label_cha,need_lim,perspective,round_bonds,
     &                autoscl,showhb,havehb
      common/frame/scale,xori,yori,max_z,min_z
      common/mol_num/numat,num(natom),co(3,natom),atrad(natom)
      common/multi_mol/mol_numat(100),num_mol
      common/mol_char/atlabel(natom)
      common/per_tbl_num/radii(105)
      common/per_tbl_char/at_sym(105)
      common/show/show_atom(natom)
      common/commands/num_cmd,cmd_line(100)
c
c     Get from the first line the general specifications.
c
      more_data = .false.
      read(9,'(A)',err=999,end=999) line
      call freeread(line,nint,nrea,ncha,inva,rea,cha,form)
      rad_fact=rea(1)
      bond_fact=rea(2)
c
c     Read the atom specifications.
c
      i=0
      in_at = 0
      num_mol = 1
      num_cmd = 0
      top_z=-9999.9
      bottom_z=-top_z
      continue=.true.
      do while (continue)
           read(9,'(A)',end=1000) line
           if (line(1:4).eq.'----') then
                continue = .false.
                more_data = .true.
                goto 1000
           endif
         if (line(1:1).ne.'#') then
           cha(1)(1:1) = ' '
           call freeread(line,nint,nrea,ncha,inva,rea,cha,form)
           if (cha(1)(1:1).eq.'#') then
c
c             the line was marked as a comment, ignore
c
           else if (cha(1)(1:1).eq.'%') then
c
c             this line is used as a separator between multiple
c             molecules to be shown at the same time in the screen.
c
              mol_numat(num_mol) = i-in_at
              in_at = i
              num_mol = num_mol+1
           else if (cha(1)(1:1).eq.'$') then
c
c             the line was marked as a command, save in an array :
c
              call limits(line,start,end)
              num_cmd = num_cmd+1
              cmd_line(num_cmd)=line((start+1):end)
           else
c
c               the line may be an atom specification.
c
                if (nrea.ge.3) then
                     i=i+1
                     num(i)=inva(1)
                     do j=1,3
                          co(j,i)=rea(j)
                     enddo
                     if (co(3,i).gt.top_z) top_z = co(3,i)
                     if (co(3,i).lt.bottom_z) bottom_z = co(3,i)
                     if (ncha .ge. 1) then
                        atlabel(i)=cha(1)(1:4)
                     else
                        atlabel(i)=at_sym(num(i))
                     endif
                     atrad(i)=rad_fact*radii(num(i))
                     show_atom(i)=.true.
                 endif
           endif
         endif
      enddo
 1000 numat=i
      mol_numat(num_mol) = i-in_at
      min_z = bottom_z
      max_z = top_z
      return
  999 more_data = .false.
      numat = 0
      return
      end
c-----------------------------------------------------------------------
      subroutine center
c-----------------------------------------------------------------------
c
      parameter (natom = 10000)
      parameter (nbond = 50000)
c
      common/mol_num/numat,num(natom),co(3,natom),atrad(natom)
      common/draw_spec/rad_fact,bond_fact,top_z,bottom_z
c
      xmax=co(1,1)
      ymax=co(2,1)
      zmax=co(3,1)
      xmin=xmax
      ymin=ymax
      zmin=zmax
c
      do i=2,numat
           xmax=max(xmax,co(1,i))
           xmin=min(xmin,co(1,i))
           ymax=max(ymax,co(2,i))
           ymin=min(ymin,co(2,i))
           zmax=max(zmax,co(3,i))
           zmin=min(zmin,co(3,i))
      enddo
c
      xcenter=(xmax+xmin)*0.5
      ycenter=(ymax+ymin)*0.5
      zcenter=(zmax+zmin)*0.5
c
      do i=1,numat
           co(1,i)=co(1,i)-xcenter
           co(2,i)=co(2,i)-ycenter
           co(3,i)=co(3,i)-zcenter
      enddo
      top_z = top_z - zcenter
      bottom_z = bottom_z - zcenter
      return
      end
c-----------------------------------------------------------------------
      subroutine init_mind(FULL,filename,command)
c-----------------------------------------------------------------------
      parameter (natom = 10000)
      parameter (nbond = 50000)
c
      logical wireframe,ball_and_stick,need_sort,label_num,
     &        label_cha,first_draw,need_lim,perspective,
     &        autoscl,showhb,havehb
      logical exists,more_data,round_bonds
      logical FULL
      character*2 command
      character*80 filename
c
      common/j_f_b/more_data
      common/draw_log/wireframe,ball_and_stick,need_sort,label_num,
     &                label_cha,need_lim,perspective,round_bonds,
     &                autoscl,showhb,havehb
      common /hyd_bonds/ hbndlength,hbangle,hbnd_d(natom/3),
     &			 hbnd_h(natom/3),hbnd_a(natom/3),
     &			 nhbdon,nhbhyd,nhbacc
      common/other_bonds/initial,mod_bonds(nbond)
      common/mol_num/numat,num(natom),co(3,natom),atrad(natom)
      common/perspective/eye_to_scr,scr_to_top
      common/font_num/ifont
      common/useratt/user_scale
c
c
      initial = 1
      if (FULL) then
         inquire(file=filename,exist=exists)
         round_bonds = .false.
         user_scale = 1.0
         if(.not. exists) return 
         open(9,file=filename,form='formatted',status='old')
         first_draw=.false.
         perspective=.false.
         wireframe=.true.
         ifont = 0
         call fontstyle(ifont)
         ball_and_stick=.false.
         label_num=.false.
         label_cha=.false.
      endif
      call init_labels
      call init_arcs
      call init_lines
c
      need_sort  = .true.
      need_lim   = .true.
      showhb     = .false.
      havehb     = .false.
      hbndlength = 2.5e0
      hbangle    = 120.0e0
c
c     Read the input data :
c
      call notify('Reading data ... ')
      call readin
      if (numat .eq. 0) then
         close(9)
         command = '  '
         call notify('No additional data in this file. ')
         return
      endif          
      if (.not.more_data) close(9)
      call notify('Generating connectivity lists ... ')
      call center
      call get_limits(3,zin,zax)
      scr_to_top = 0.0e0
      eye_to_scr = 2.0e0*(zax - zin)
c
c     Determine the connectivity :
c
      call get_bonds
      call bond_at_list
      call notify('Generating HB lists ...')
      call hbpotlist
      call notify('Full view. ')
      return
      end
c----------------------------------------------------------------------
      block data general
c----------------------------------------------------------------------
c
      common /grid/grid_on
      logical grid_on
      integer reversed
      data grid_on, reversed / .false., 0 /
c
c     Initialization of the "PERiodic TaBLe" commons.
c
      character*4 at_sym
      common/per_tbl_num/radii(105)
      common/per_tbl_char/at_sym(105)
      common/per_color/iatcolor(105), reversed
c
c Covalent Radii + 10%
c
      data radii /
     &    0.370,0.350,1.353,0.990,0.902,0.820,0.790,0.750,0.792,0.781,
     &    1.694,1.496,1.298,1.221,1.166,1.090,1.040,1.078,2.233,1.914,
     &    1.584,1.452,1.342,1.298,1.287,1.287,1.276,1.265,1.287,1.375,
     &    1.386,1.342,1.320,1.276,1.254,1.232,2.376,2.101,1.782,1.595,
     &    1.474,1.430,1.397,1.375,1.375,1.408,1.474,1.628,1.584,1.551,
     &    1.540,1.496,1.463,1.441,2.585,2.178,1.859,1.815,1.815,1.804,
     &    1.793,1.782,2.035,1.771,1.749,1.749,1.738,1.727,1.716,1.716,
     &    1.716,1.584,1.474,1.430,1.408,1.386,1.397,1.430,1.474,1.639,
     &    1.628,1.617,1.606,1.606,1.595,1.595,1.650,1.650,1.760,1.815,
     &    15*0.0 /
      data at_sym/'H   ','He  ','Li  ','Be  ','B   ','C   ','N   ',
     &            'O   ','F   ','Ne  ','Na  ','Mg  ','Al  ','Si  ',
     &            'P   ','S   ','Cl  ','Ar  ','K   ','Ca  ','Sc  ',
     &            'Ti  ','V   ','Cr  ','Mn  ','Fe  ','Co  ','Ni  ',
     &            'Cu  ','Zn  ','Ga  ','Ge  ','As  ','Se  ','Br  ',
     &            'Kr  ','Rb  ','Sr  ','Y   ','Zr  ','Nb  ','Mo  ',
     &            'Tc  ','Ru  ','Rh  ','Pd  ','Ag  ','Cd  ','In  ',
     &            'Sn  ','Sb  ','Te  ','I   ','Xe  ','Cs  ','Ba  ',
     &            'La  ','Ce  ','Pr  ','Nd  ','Pm  ','Sm  ','Eu  ',
     &            'Gd  ','Tb  ','Dy  ','Ho  ','Er  ','Tm  ','Yb  ',
     &            'Lu  ','Hf  ','Ta  ','W   ','Re  ','Os  ','Ir  ',
     &            'Pt  ','Au  ','Hg  ','Tl  ','Pb  ','Bi  ','Po  ',
     &            'At  ','Rn  ','Fr  ','Ra  ','Ac  ','Th  ',
     &         15*'    '/
c
c The colors are:
c  1       black
c  2       red
c  3       yellow
c  4       green
c  5       cyan
c  6       blue
c  7       magenta
c  8       white
      data iatcolor / 4, 7, 6, 6, 5, 1, 6, 2, 7, 7, 
     &                3, 6, 6, 7, 7, 3, 5, 2, 3, 4,
     &                85*6 /
      end