PROGRAM DODC
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (MX=600)
      DIMENSION R(3,MX),rc(3),ic(MX),wc(MX),
     1a(3),b(3),c(3),bb(3),u(3)
      logical lconfig
c
      inquire(file='DODC.CON',exist=lconfig)
      open(30,file='dodc.scr')
      WRITE(6,*)' Do DC.SC for DCO.'
      OPEN(2,FILE='FILE.X',STATUS='OLD',err=999)
      read(2,*)
      read(2,*)N
      if(N.gt.MX)call report('Too many atoms.')
      do 1 i=1,N
1     READ(2,*)idumm,(R(J,i),J=1,3)
      close(2)
      write(6,*)' FILE.X read,',N,' atoms'
      open(3,file='DC.SC')
      if(lconfig)then
       open(4,file='DODC.CON')
       iin=4
      else
       iin=5
      endif
      write(6,*)'How many groups:'
      read(iin,*)NG
      write(30,*)NG,' groups'
      write(3,*)NG,' groups'
      do 2 IG=1,NG
      write(6,*)IG,' . group:'
      write(3,*)IG,' . group:'
      write(6,*)'Number of atoms determining center and their numbers:'
      read(iin,*)nc,(ic(i),i=1,nc)
      write(30,*)nc,(ic(i),i=1,nc)
      write(6,*)'Input ',nc,' center weights:'
      read(iin,*)(wc(i),i=1,nc)
      write(30,*)(wc(i),i=1,nc)
      do 3 i=1,3
      rc(i)=0.0d0
      do 3 ii=1,nc
3     rc(i)=rc(i)+R(i,ic(ii))*wc(ii)
      sum=0.0d0
      do 32 ii=1,nc
32    sum=sum+wc(ii)
      do 31 i=1,nc
31    rc(i)=rc(i)/sum
      write(6,1000)(rc(i)/10.0d0,i=1,3)
      write(3,1000)(rc(i)/10.0d0,i=1,3)
1000  format(3f15.6,' nm')
      write(6,*)'Number of transitions:'
      read(iin,*)ns
      write(30,*)ns
      ns=ns+1
      write(6,*)ns,' states'
      write(3,*)ns,' states'
      write(3,*)
      write(3,*)'0 0.0'
      do 4 is=2,ns
      write(6,*)'Frequency of transition ',is-1,' (cm-1) :'
      read(iin,*)w
      write(30,*)w
      write(6,1001)is,w/8065.54093D0,w
1001  format(i3,f15.6,' eV = ',f15.6,' cm-1')
4     write(3,1001)is,w/8065.54093D0,w
      write(3,*)
      do 5 is=2,ns
      write(6,*)'Dipole strength of transition ',is-1,' (Debye^2) :'
      read(iin,*)d
      write(30,*)d
      ud=dsqrt(d)
      write(6,*)'Four atom numbers defining plane:'
      write(6,*)'    (Vectors a=1->2, b=3->4)'
      read(iin,*)i1,i2,i3,i4
      write(30,*)i1,i2,i3,i4
      do 6 ix=1,3
      a(ix)=R(ix,i2)-R(ix,i1)
6     bb(ix)=R(ix,i4)-R(ix,i3)
      call norm(a)
      call vp(bb,a,c)
      call norm(c)
      call vp(a,c,b)
      write(6,*)'Deviation from a = ',i1,' -> ',i2,' (deg):'
      read(iin,*)alpha
      write(30,*)alpha
      alpha=alpha*3.141592653589d0/180.0d0
      ua=ud*cos(alpha)
      ub=ud*sin(alpha)
      do 7 ix=1,3
7     u(ix)=ua*a(ix)+ub*b(ix)
      write(3,1002)is,(u(ix),ix=1,3)
5     write(6,1002)is,(u(ix),ix=1,3)
1002  format(i3,3f15.6,' 0.0 0.0 0.0')
2     continue
      close(3)
      write(6,*)'DC.SC written'
      if(lconfig)close(iin)
      close(30)
      stop
999   call report('FILE.X not found')
      stop
      end
c
      SUBROUTINE REPORT(RS)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      CHARACTER*(*) RS
      WRITE(6,3000)
3000  FORMAT(/,80(1H*))
      WRITE(6,*)RS
      WRITE(6,3001)
3001  FORMAT(80(1H*),/,/,'PROGRAM STOPPED')
      STOP
      END

      SUBROUTINE norm(a)
      REAL*8 a(3),an
      an=dsqrt(a(1)**2+a(2)**2+a(3)**2)
      a(1)=a(1)/an
      a(2)=a(2)/an
      a(3)=a(3)/an
      return
      end

      SUBROUTINE vp(a,b,c)
      REAL*8 a(3),b(3),c(3)
      c(1)=a(2)*b(3)-a(3)*b(2)
      c(2)=a(3)*b(1)-a(1)*b(3)
      c(3)=a(1)*b(2)-a(2)*b(1)
      return
      end