program mdip
      implicit none
      integer*4 NQ,NOAT,K,I,J,NC
      real*8 BOHR,CM,CLIGHT,rc(3),u0(3),m0(3),ur(3),um(3),
     1c,a1,a2,mp(3),me(3),d(3),mm(3),sp
      real*8,allocatable:: R(:,:),E(:),u(:,:),m(:,:)
      integer*4,allocatable:: la(:)
      character*80 s80
c     investigates origin of magnetic dipoles
      BOHR=0.52917715d0
      CM=219474.0d0
      CLIGHT=137.035999139d0

      OPEN(17,FILE='F.INP',STATUS='OLD')
      READ(17,*)NQ,NOAT,NOAT
      allocate(R(3,NOAT),E(NQ),u(NQ,3),m(NQ,3))
      DO 3 K=1,NOAT
      READ(17,*)R(1,K),(R(I,K),I=1,3)
      DO 3 I=1,3
3     R(I,K)=R(I,K)/BOHR
      READ(17,*)
      DO 4 K=1,NOAT
      DO 4 I=1,NQ
4     READ (17,*)
      READ(17,*)
      READ(17,*)(E(I),I=1,NQ)
      do 1 I=1,NQ
      E(I)=E(I)/CM
1     read(17,*)(u(I,J),J=1,3),(m(I,J),J=1,3)
      close(17)
      WRITE(*,*)' S-matrix read in'

      open(17,file='MDIP.OPT')
2     read(17,1700,end=99,err=99)s80
1700  format(a80)
      if(s80(1:4).eq.'NCEN')then
       read(17,*)NC
       allocate(la(NC))
       read(17,*)(la(I),I=1,NC)
       do 5 J=1,3
       rc(J)=0.0d0
       do 6 I=1,NC
6      rc(J)=rc(J)+R(J,la(I))
5      rc(J)=rc(J)/dble(NC)
       write(6,601)(rc(J)*BOHR,J=1,3)
601    format(' Center: ',3F12.4,' A',/,
     1 ' Dipoles at common origin and at origin at the center:')
       do 7 I=1,NQ
       c=0.5d0*E(I)/CLIGHT
       do 8 J=1,3
       u0(J)=u(I,J)
8      m0(J)=m(I,J)
       call vp(u0,rc,ur)
c      m' = m + (1/2) w r x u
       do 9 J=1,3
9      m0(J)=m0(J)+c*ur(J)
7      write(6,602)I,E(I)*CM,(m(I,J),J=1,3),(m0(J),J=1,3)
602    format(i4,F12.2,2(3G12.4,'  '))
      endif
      goto 2
99    close(17)

      write(6,603)
603   format(' Minimal positions /A   new dip  Aold Anew:')
      do 11 I=1,NQ
      c=0.5d0*E(I)/CLIGHT
      do 10 J=1,3
      d(J)=u(I,J)
      mm(J)=m(I,J)
      u0(J)=d(J)
10    m0(J)=mm(J)
      a1=dsqrt(sp(mm,mm))
      call vnorm(u0)
      call vnorm(m0)
c     parallel (mp) and perpedicular (me) parts w/r to u:
      do 12 J=1,3
      mp(J)=u0(J)*sp(u0,mm)
12    me(J)=mm(J)-mp(J)
c     um: shift direction
      call vp(d,me,um)
      call vnorm(um)
c     shiftL
      do 13 J=1,3
13    rc(J)=um(J)*dsqrt(sp(me,me)/sp(d,d))/c
c     new moment:
      call vp(d,rc,ur)
      do 14 J=1,3
14    mm(J)=mm(J)+c*ur(J)
      a2=dsqrt(sp(mm,mm))
11    write(6,604)I,E(I)*CM,(rc(J),J=1,3),mm,a1,a2
604   format(i3,F10.2,3f10.3,5G11.3)

      end

      subroutine vp(a,b,c)
      implicit none
      real*8 a(*),b(*),c(*)
      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


      function sp(a,b)
      implicit none
      real*8 a(*),b(*),sp
      sp=a(1)*b(1)+a(2)*b(2)+a(3)*b(3)
      return
      end

      subroutine vnorm(k)
      implicit none
      real*8 k(*),k2,sp
      k2=dsqrt(sp(k,k))
      k(1)=k(1)/k2
      k(2)=k(2)/k2
      k(3)=k(3)/k2
      return
      end