program cpropag3
      implicit none
      integer*4 nat,n,i,ix,iy,n27,nats,ia,ja,iccub,ice,ie,iofe,nx,ny,nz,
     1ex,ey,ez,cxp,cyp,czp,jp,iargc,iia,natb,nb,nbx,nby,nbz,io,ioa,
     1iofio,iofip,iofjo,iofjp,ip,ipa,jo,joa,jpa
      real*8 a(3),b(3),c(3),ax,ay,az
      real*8 ,allocatable::r(:),rb(:),f(:,:),P(:,:,:),AA(:,:,:),fb(:,:)
      integer*4 ,allocatable::pp(:,:),ind(:),q(:)
      character*80 s80
      LOGICAL LAPT

      LAPT=.false.
      write(6,600)
600   format(' From a crystal piece distributes force field',/,
     1       '           to a new structure',/,/,
     1       '       Input:  CELL.TXT   cell vectors',/,
     1       '               CELL.X     cell geometry',/,
     1       '               FILE.X     geometry',/,
     1       '               FILE.FC    force field',/,/,
     1       '               FILE.TEN   tensors',/,/,
     1       '      Output:  BIG.FC     new force field',/,
     1       '               BIG.X      new geometry',/,/)

      if(iargc().ne.3)call report('Usage: cpropag3 Nx Ny Nz')

      call getarg(1,s80)
      read(s80,*)nbx
      call getarg(2,s80)
      read(s80,*)nby
      call getarg(3,s80)
      read(s80,*)nbz
      write(6,601)nbx,nby,nbz
601   format(i2,' x ',i2,' x ',i2)

      open(9,file='CELL.TXT')
      read(9,*)a
      read(9,*)b
      read(9,*)c
      close(9)

      open(9,file='CELL.X')
      read(9,*)
      read(9,*)nats
      close(9)

      open(9,file='FILE.X')
      read(9,*)
      read(9,*)nat
      n=3*nat
      allocate(r(n),f(n,n),P(nat,3,3),AA(nat,3,3),q(nat))
      do 1 i=1,nat
1     read(9,*)q(i),(r(ix+3*(i-1)),ix=1,3)
      close(9)

      CALL READFF(n,f)

      n27=nat/nats
      allocate(pp(n27,3))
      pp(1,1)=0

c     determine central cube:
      ice=iccub(nat,r,n27,pp,a,b,c)
      iofe=nats*(ice-1)

c     new structure:
      natb=nats*nbx*nby*nbz
      write(6,*)natb,' atoms in the BIG file'
      nb=3*natb
      allocate(fb(nb,nb),ind(natb),rb(nb))
      fb=0.0d0

      OPEN(20,FILE='BIG.X')
      write(20,601)nbx,nby,nbz
      write(20,*)natb
      iia=0
      do 9 ex=0,nbx-1
      do 9 ey=0,nby-1
      do 9 ez=0,nbz-1
c     position of cell (ex ey ez):
      ax=dble(ex)*a(1)+dble(ey)*b(1)+dble(ez)*c(1)
      ay=dble(ex)*a(2)+dble(ey)*b(2)+dble(ez)*c(2)
      az=dble(ex)*a(3)+dble(ey)*b(3)+dble(ez)*c(3)
      do 9 i=1,nats
      iia=iia+1
      ie =i+iofe
      ind(iia)=ie
      rb(3*iia-2)=r(3*ie-2)+ax
      rb(3*iia-1)=r(3*ie-1)+ay
      rb(3*iia  )=r(3*ie  )+az
9     write(20,20)q(ie),(rb(3*iia+ix),ix=-2,0)
20    format(i2,3f12.6,' 0 0 0 0 0 0 0 0.0')
      close(20)
      write(6,*)' BIG.X was written.'
   

c     old cells                    io - jo:
      do 2 io=1,n27
      iofio=nats*(io-1)
      do 2 jo=1,n27
      iofjo=nats*(jo-1)

      write(6,602)io,jo
 602  format(2i3,1h:)
c     relative position of old cells:
      nx=pp(jo,1)-pp(io,1)
      ny=pp(jo,2)-pp(io,2)
      nz=pp(jo,3)-pp(io,3)

c     ice...ic to all other pairs ip     jp:
c                              iofip  iofjp
c     new cells:
      ip=0
      do 2 ex=0,nbx-1
      do 2 ey=0,nby-1
      do 2 ez=0,nbz-1
      ip=ip+1
      iofip=nats*(ip-1)
      jp=0
      do 2 cxp=0,nbx-1
      do 2 cyp=0,nby-1
      do 2 czp=0,nbz-1
      jp=jp+1
      iofjp=nats*(jp-1)
      if(nx.eq.cxp-ex.and.ny.eq.cyp-ey.and.nz.eq.czp-ez)then
       write(6,603)ip,jp
 603   format(' to ',2i3)
       do 8 ia=1,nats
       ioa=ia+iofio
       ipa=ia+iofip
       do 8 ja=1,nats
       joa=ja+iofjo
       jpa=ja+iofjp
       do 8 ix=1,3
       do 8 iy=1,3
 8     fb(ix+3*(ipa-1),iy+3*(jpa-1))=f(ix+3*(ioa-1),iy+3*(joa-1))
      endif
 2    continue

      OPEN(20,FILE='BIG.FC')
      CALL WRITEFF(nb,fb)
      CLOSE(20)
      write(6,*)'BIG.FC written.'

      CALL READTEN(NAT,P,AA,LAPT,r)

      call WRITETEN(P,AA,rb,nat,natb,IND)

      end

      SUBROUTINE READTEN(NAT,P,A,LAPT,X)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION P(NAT,3,3),A(NAT,3,3),X(*)
      LOGICAL LAPT
      real*8,allocatable::AI(:,:,:),AJ(:,:,:),PV(:,:,:),R(:,:)
      allocate(AI(NAT,3,3),AJ(NAT,3,3),PV(NAT,3,3),R(3,NAT))

      BOHR=0.52917705993d0
      DO 3 L=1,NAT
      DO 3 I=1,3
3     R(I,L)=X(I+3*(L-1))/BOHR

      OPEN(15,FILE='FILE.TEN',STATUS='OLD')
      READ (15,*)
C
C     Read dipole derivatives:
      DO 10 L=1,NAT
      DO 10 J=1,3                                                     
10    READ (15,*) (P(L,J,I),I=1,3)
      WRITE(6,*)' Dipole derivatives read - in '
C
C     If no axial tensor, then construct its polar part only:
      IF(LAPT)THEN
       DO 2203 L=1,NAT
       DO 2203 J=1,3
       DO 2203 I=1,3
       SUM=0.0d0
       DO 2204 ID=1,3
       DO 2204 IG=1,3
2204   SUM=SUM+0.25d0*EPS(I,IG,ID)*R(IG,L)*P(L,J,ID)
       AJ(L,J,I)=0.0d0
       PV(L,J,I)=P(L,J,I)
2203   AI(L,J,I)=SUM
       WRITE(6,*)' APT part of AAT constructed'
       GOTO 1
      ENDIF

      DO 220 L=1,NAT                                                 
      DO 220 J=1,3                                                    
220   READ (15,*) (AI(L,J,I),I=1,3)
      DO 230 L=1,NAT                                                 
      DO 230 J=1,3                                                    
230   READ (15,*) (AJ(L,J,I),I=1,3)                                   
      DO 100 L=1,NAT                                                 
      DO 100 J=1,3                                                    
100   READ (15,*) (PV(L,J,I),I=1,3)                                   
      WRITE(6,*)' VCD parameters read in ...'
C     Distributed-origin gauge:
1     CALL VCDD0(AI,PV,P,R,NAT)   
      WRITE(6,*)'DOG used'
C     
C     Total axial tensor = electronic + nuclear:
      DO 4 L=1,NAT
      DO 4 I=1,3
      DO 4 J=1,3
4     A(L,I,J)=AI(L,I,J)+AJ(L,I,J)
C
C     Put axial tensor in the origin on the atom L:
      DO 2201 L=1,NAT
      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)*R(IG,L)*P(L,J,ID)
2201  A(L,J,I)=A(L,J,I)+SUM
      WRITE(6,*)' AAT shifted to atomic origin'

      CLOSE(15)
      WRITE(6,*)' FILE.TEN read ...'
      RETURN
      END
C
      SUBROUTINE VCDD0(VCD,VEL,DIP,C,NAT)
C     STOLEN FROM CADPAC, INDICES ARE DIFFERENT !!, VEL WITH NUCLEI
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      INTEGER*4 ALPHA,BETA,GAMMA,DELTA,LAMBDA                           
      DIMENSION VCD(NAT,3,3),VEL(NAT,3,3),DIP(NAT,3,3),                  
     1C(3,NAT)
      DO 90 BETA =1,3                                                 
      DO 90 GAMMA=1,3                                                 
      DO 90 DELTA=1,3                                                 
      SKEW=0.25D0*EPS(BETA,GAMMA,DELTA)                               
      IF (DABS(SKEW).GT.1.0D-10)THEN                                          
      DO 80 ALPHA=1,3                                                 
      DO 80 LAMBDA=1,NAT                                              
80    VCD(LAMBDA,ALPHA,BETA)=VCD(LAMBDA,ALPHA,BETA)                   
     1 -SKEW*C(GAMMA,LAMBDA)                                          
     2 *(VEL(LAMBDA,ALPHA,DELTA)-DIP(LAMBDA,ALPHA,DELTA))             
      ENDIF                                                           
90    CONTINUE
      RETURN
      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                                                             
C 

      SUBROUTINE WRITETEN(P,A,X,NAT,NATN,IND)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION P(NAT,3,3),A(NAT,3,3),X(*),IND(*)
      real*8 ,allocatable::R(:,:),AB(:,:,:)

      allocate(R(3,NATN))
      BOHR=0.52917705993d0
      DO 3 L=1,NATN
      DO 3 I=1,3
3     R(I,L)=X(I+3*(L-1))/BOHR

      allocate(AB(NATN,3,3))
      DO 2201 Ln=1,NATN
      L=IND(Ln)
      DO 2201 J=1,3
      DO 2201 I=1,3
      AB(Ln,J,I)=A(L,J,I)
      DO 2201 ID=1,3
      DO 2201 IG=1,3
2201  AB(Ln,J,I)=AB(Ln,J,I)+0.25d0*EPS(I,IG,ID)*R(IG,Ln)*P(L,J,ID)
      WRITE(6,*)' AAT shifted back to laboratory system'

      OPEN(15,FILE='BIG.TEN')
      WRITE(15,1500) NATN,NATN-6,0
1500  FORMAT(3I5)

      DO 10 Ln=1,NATN
      L=IND(Ln)
      DO 10 J=1,3
10    WRITE(15,1501) (P(L,J,I),I=1,3),L
1501  FORMAT(3F14.8,I5)

      DO 220 L=1,NATN
      DO 220 J=1,3                                                    
220   WRITE(15,1501) (AB(L,J,I),I=1,3),L

      Z0=0.0d0
      DO 230 L=1,NATN                                                 
      DO 230 J=1,3                                                    
230   WRITE(15,1501) (Z0,I=1,3),L

      DO 100 Ln=1,NATN
      L=IND(Ln)
      DO 100 J=1,3
100   WRITE(15,1501) (P(L,J,I),I=1,3),L

      CLOSE(15)
      WRITE(6,*)' BIG.TEN written ...'
      RETURN
      END
      
      SUBROUTINE READFF(N,FCAR)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION FCAR(N,N)
      OPEN(20,FILE='FILE.FC',STATUS='OLD')
      N1=1
1     N3=N1+4
      IF(N3.GT.N)N3=N
      DO 130 LN=N1,N
130   READ(20,17)(FCAR(LN,J),J=N1,MIN(LN,N3))
      N1=N1+5
      IF(N3.LT.N)GOTO 1
17    FORMAT(4X,5D14.6)
      DO 31 I=1,N
      DO 31 J=I+1,N
   31 FCAR(I,J)=FCAR(J,I)
      WRITE(*,*)' Cartesian FF read in ... '
      CLOSE(20)
      RETURN
      END

      SUBROUTINE WRITEFF(N,FCAR)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION FCAR(N,N)
      N1=1
1     N3=N1+4
      IF(N3.GT.N)N3=N
      DO 130 LN=N1,N
130   WRITE(20,17)LN,(FCAR(LN,J),J=N1,MIN(LN,N3))
      N1=N1+5
      IF(N3.LT.N)GOTO 1
17    FORMAT(I4,5D14.6)
      RETURN
      END
 

      function iccub(nat,r,n27,ip,wi,wj,wk)
c     determine central cube
      implicit none
      integer*4 nat,nat1,ix,iu,ii,iccub,ice,i,n27,ip(n27,3)
      real*8 cet(3),r(*),xi,yi,zi,d,wi(3),wj(3),wk(3),down,
     1vi,vj,vk,vii,vij,vik,vjj,vjk,vkk,sp,vu(3),a,b,c
      real*8,allocatable::cei(:,:)
      allocate(cei(n27,3))

      nat1=nat/n27
c     center of all and of each one:
      do 40 ix=1,3
      cet(ix)=0.0d0
      do 42 iu=1,n27
      cei(iu,ix)=0.0d0
      ii=nat1*(iu-1)
      do 41 i=ii+1,ii+nat1
      cet(   ix)=cet(   ix)+r(3*(i-1)+ix)
41    cei(iu,ix)=cei(iu,ix)+r(3*(i-1)+ix)
42    cei(iu,ix)=cei(iu,ix)/dble(nat1)
40    cet(ix)=cet(ix)/dble(nat)

      write(6,600)cet
600   format(/,' CRYSTAL.X center (A):',7x,3f7.2,/,' Cube centers (A):')
      write(6,601)(iu,(cei(iu,ix),ix=1,3),iu=1,n27)
601   format(3(i3,':',3f7.2))
      d=1.0d10
      ice=0
      do 43 iu=1,n27
      xi=cei(iu,1)-cet(1)
      yi=cei(iu,2)-cet(2)
      zi=cei(iu,3)-cet(3)
      if(dsqrt(xi**2+yi**2+zi**2).lt.d)then
       d=dsqrt(xi**2+yi**2+zi**2)
       ice=iu
      endif
43    continue

      if(ice.eq.0)then
       call report('Central cube not found')
      else
       write(6,602)ice,d
602    format(' Central cube:',i3,f10.2,' A from the center')
       iccub=ice
      endif

c     relative positions to the central cell
      do 1 iu=1,n27
      do 2 ix=1,3
2     vu(ix)=cei(iu,ix)-cei(ice,ix)
      vi=sp(vu,wi)
      vj=sp(vu,wj)
      vk=sp(vu,wk)
      vii=sp(wi,wi)
      vij=sp(wi,wj)
      vik=sp(wi,wk)
      vjj=sp(wj,wj)
      vjk=sp(wj,wk)
      vkk=sp(wk,wk)
      down=(vkk*vii-vik**2)*(vjj*vii-vij**2)
     1-(vjk*vii-vij*vik)*(vjk*vii-vik*vij)
c     write(6,*)vi,vj,vk,vii,vij,vik,vjj,vjk,vkk,down
      c=(vi*(vjk*vij-vjj*vik)+vj*(vij*vik-vjk*vii)+vk*(vjj*vii-vij**2))*
     1vii/down
      b=(vj*vii-vi*vij-c*(vjk*vii-vik*vij))/(vjj*vii-vij**2)
      a=(vi-b*vij-c*vik)/vii
      write(6,603)iu,a,b,c
 603  format(i3,3f5.1,$)
      if(mod(iu,4).eq.0)write(6,*)
      ip(iu,1)=nint(a)
      ip(iu,2)=nint(b)
1     ip(iu,3)=nint(c)
      if(mod(n27,4).ne.0)write(6,*)

      return
      end

      SUBROUTINE report(s)
      character*(*) s
      write(6,*)s
      stop
      end

      function sp(v1,v2)
      real*8 v1(*),v2(*),sp        
      sp=v1(1)*v2(1)+v1(2)*v2(2)+v1(3)*v2(3)
      return
      end