PROGRAM gfac
      implicit none
      integer*4 N0,NOAT,I,J,K,ia,jb,l,NOSC,ix,id,ig
      PARAMETER (N0=4500)
      REAL*8 AAT(N0,3,3),aii,BM,sum,R(3,N0),chatot,gn(3,3),
     4BOHR,am(N0),char(N0),cm(3),P(N0,3,3),eps
c     
      BOHR=0.52917715d0
      WRITE(*,*)' Rotational g-factor from AAT'
      WRITE(*,6666)
6666  FORMAT(/,
     1'  INPUT FILES: FILE.X   ... geometry',/,
     1'               FTRY.INP ... masses',/,
     1'               FILE.TEN ... TENSORS',/)

      open(15,file='FILE.X',status='old')
      read(15,*)
      read(15,*)NOAT
      do 1 i=1,NOAT
      read(15,*)j,(R(ix,i),ix=1,3)
      do 1 ix=1,3
1     R(ix,i)=R(ix,i)/BOHR
      close(15)
      write(6,*)NOAT,' atoms'

      open(15,file='FILE.TEN',status='old')
      read(15,*)
C
C  READ DIPOLE DERIVATIVES = TENSOR P(A,B)
C
      DO 10 L=1,NOAT
      DO 10 J=1,3
10    READ (15,*)  (P(L,J,I),I=1,3)
C
C  READ TENSORS I(A,B) AND J(A,B) FOR THE 0 0 0        ORIGIN:
C    (=ELECTRONIC & NUCLEAR TERMS)
C
      DO 220 L=1,NOAT
      DO 220 J=1,3
 220  READ (15,*) (AAT(L,J,I),I=1,3)
      DO 230 L=1,NOAT
      DO 230 J=1,3
230   READ (15,*)
C
C  READ VELOCITY REPRESENTATION OF DIPOLE DERIVATIVES
C
      DO 100 L=1,NOAT
      DO 100 J=1,3
100   READ (15,*)
      CLOSE(15)

      WRITE(*,*)' AAT read from FILE.TEN'
C
      OPEN(15,FILE='FTRY.INP')
      READ(15,*)
      READ(15,*)NOSC
      DO 222 I=1,NOSC
222   READ(15,*) 
C
C     Read-in atomic charges for FPC:
      READ(15,*) (CHAR(I),I=1,NOAT)
      CHATOT=0.0d0
      DO 77 I=1,NOAT
77    CHATOT=CHATOT+CHAR(I)
      WRITE(*,777)CHATOT
777   FORMAT(' TOTAL CHARGE',F20.6)
C
C     Read-in atomic masses:
      READ(15,*)
c     READ(15,1040)(am(I),I=1,NOAT)
      READ(15,*)(am(I),I=1,NOAT)
      BM=0.0d0
      DO 102 IA=1,NOAT
102   BM=BM+am(IA)
      WRITE(*,778)BM
778   FORMAT(' TOTAL MASS  ',F20.6)
      CLOSE(15)
      WRITE(*,*)' FTRY.INP READ IN'

      do 3 ix=1,3
      cm(ix)=0.0d0
      do 31 ia=1,NOAT
31    cm(ix)=cm(ix)+R(ix,ia)*am(ia)
3     cm(ix)=cm(ix)/BM
      write(6,606)cm
606   format(/' Mass center:',3f12.6,/)

      do 4 i=1,NOAT
      do 4 ix=1,3
4     R(ix,i)=R(ix,i)-cm(ix)

      DO 2201 L=1,NOAT
      DO 2201 J=1,3
      DO 2201 I=1,3
      SUM=0.0d0
      DO 2202 ID=1,3
      DO 2202 IG=1,3
2202  SUM=SUM-0.25d0*EPS(I,IG,ID)*CM(IG)*P(L,J,ID)
2201  AAT(L,J,I)=AAT(L,J,I)+SUM
      WRITE(6,*)' AAT shifted to mass center'

      do 2 ia=1,3
      do 2 jb=1,3
      j=jb+1
      if(j.gt.3)j=1
      k=j+1
      if(k.gt.3)k=1
      aii=0.0d0
      sum=0.0d0
      do 21 i=1,NOAT
      aii=aii+r(j,i)**2*am(i)+r(k,i)**2*am(i)
21    sum=sum+(AAT(i,j,ia)*R(k,i)-AAT(i,k,ia)*R(j,i))
2     gn(ia,jb)=sum/aii
      write(6,600)((gn(ia,jb),ia=1,3),jb=1,3),
     1(gn(1,1)+gn(2,2)+gn(3,3))/3.0d0
600   format('Rotational g-tensor:',/,/,
     13(3F12.6,/),/,'Isotropic:',/,/,f12.6,/)

      STOP
      END

      FUNCTION EPS(I,J,K)
      REAL*8 EPS
      INTEGER*4 I,J,K
      EPS=0.0d0
      IF (I.EQ.1.AND.J.EQ.2.AND.K.EQ.3)EPS= 1.0d0
      IF (I.EQ.1.AND.J.EQ.3.AND.K.EQ.2)EPS=-1.0d0
      IF (I.EQ.2.AND.J.EQ.3.AND.K.EQ.1)EPS= 1.0d0
      IF (I.EQ.2.AND.J.EQ.1.AND.K.EQ.3)EPS=-1.0d0
      IF (I.EQ.3.AND.J.EQ.1.AND.K.EQ.2)EPS= 1.0d0
      IF (I.EQ.3.AND.J.EQ.2.AND.K.EQ.1)EPS=-1.0d0
      RETURN                                                          
      END