PROGRAM A22cgGS2r
C CCC for at_22 potential N=66  CCCCCCCCCCCCCCCCCCCCCCCC
C W.Quapp  30.09.2004  for RGF=Newton Trajectory
C GS nach: Baron Peters et al. JCP 120 (2004) 7877-7886
c script version - modular structure 
c Projector
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      Integer NDIM,L
      PARAMETER(NDIM=66,L=25)
      REAL*8 RotX(NDIM),RotY(NDIM),RotZ(NDIM),X(NDIM)
      REAL*8 PP(NDIM,NDIM),RR(NDIM,NDIM)
      REAL*8 Y(L+1,NDIM),RsN,ProMat(NDIM,NDIM),Rs(NDIM)
      Integer N,I,J,JJ,LA,Lha,Lnum
       OPEN(9,FILE='at22chain.dat') 
       OPEN(10,FILE='at22Rs.dat')
       OPEN(11,FILE='at22param.dat',status='old') 
       OPEN(12,FILE='at22proj.dat') 
       OPEN(13,FILE='at22point.dat') 
c  Rs: search direction, first direction between minima
C   2.version: new at every step, however: 
C   only after 3 steps, moderate adaption ! 
       N=NDIM 
       read(11,*) JJ, LA
       JJ=L+1
       rewind 9
       Do 31, J=1,L+1
31     read(9,50) (Y(J,I),I=1,N)
c  projector direction: 
       Lha= 3
       Lnum= 1
       IF(LA .gt. Lha) Lnum=LA-Lha
       Do 28 J=1,N
28     Rs(J)=Y(Lnum,J) - Y(L+1,J) 
c   Normierung
      RsN=0.0d0
      Do 32 J=1,N
32    RsN= RsN+ Rs(J)*Rs(J)
      RsN=DSQRT(RsN)
      If(RsN.eq.0.0d0) then
       write(6,*) ' '
       write(6,*) ' '
       write(6,*) 'STOP ,  search direction error '
       write(6,*) ' '
       write(6,*) ' '
       stop 
      endif 
      Do 33 J=1,N    
33    Rs(J)=Rs(J)/RsN
        rewind 10
        write(10,50)(Rs(I),I=1,N)
        rewind 13
        READ(13,50)(X(I),I=1,N)

c rotation zero directions
      do 55,I=1,N/3
           JJ=3*(I-1)
        RotX(JJ+1)=0.0d0
        RotX(JJ+2)=-X(JJ+3)
        RotX(JJ+3)= X(JJ+2)
c
        RotY(JJ+1)= X(JJ+3)
        RotY(JJ+2)=0.0d0 
        RotY(JJ+3)=-X(JJ+1)
c
        RotZ(JJ+1)=-X(JJ+2)
        RotZ(JJ+2)= X(JJ+1)
        RotZ(JJ+3)=0.0d0                  
55     continue
C   Projektionsoperator  !! P_rot is left,  P_Rs is right ??

      call PROJECTOR(N,RotX,RR)
      call PROJECTOR(N,RotY,PP)
      call matmult(PP,N,RR,N,ProMat,N)
      call PROJECTOR(N,RotZ,RR) 
      call matmult(RR,N,ProMat,N,PP,N)
ccc
ccc      call PROJECTOR(N,Rs,ProMat)
ccc
       call PROJECTOR(N,Rs,RR)
       call matmult(RR,N,PP,N,ProMat,N)
cc
      rewind 12
      Do 37, J=1,N
37    WRITE(12,50)(ProMat(J,I),I=1,N)
50    FORMAT(1X,3F20.10)
      Stop
      End
c
c ccccccccccccccccccccccccccccccccccccccccccccc
       SUBROUTINE PROJECTOR(N,Rs,ProMat) 
c  projector to vector Rs of dim N in the
c  Projection Matrix ProMat of dim (N,N) 
      INTEGER N
      REAL*8 Rs(N),ProMat(N,N),TN
      INTEGER J,JJ
c
      TN=0.0d0
      Do 32 J=1,N
32    TN= TN+ Rs(J)*Rs(J)
      TN=DSQRT(TN)
c  normalization  
      Do 33 J=1,N 
33    Rs(J)=Rs(J)/TN
C  Projektionsoperator
      Do 34 J =1,N  
      Do 34 JJ=1,N  
      ProMat(J,JJ)=0.0d0
34    continue      
      Do 35 J=1,N   
      ProMat(J,J)=1.0d0
      Do 36 JJ=1,N     
      ProMat(J,JJ)= ProMat(J,JJ) - Rs(J)*Rs(JJ)
36    continue
35    continue
      RETURN  
      END     
ccc
      subroutine matmult(m1,d1,m2,d2,mout,dim)
c  matrix multiplication:
c    m1[d1,dim]*m2[dim,d2]=mout[d1,d2]
      real*8 m1,m2,mout 
      integer d1,d2,dim,i,j,k
      dimension m1(d1,dim),m2(dim,d2),mout(d1,d2)
c
      do 10,i=1,d1
          do 20,j=1,d2
            mout(i,j)=0.0d0
            do 30,k=1,dim  
             mout(i,j)=mout(i,j)+m1(i,k)*m2(k,j)
   30       continue
  20     continue   
 10   continue      
      return        
      end