PROGRAM a162cg2r
C CCC for poly ALA15  potential N=486  CCCCCCCCCCCCCCCCCCCCCCCC
C W.Quapp 16.03.2005/07.07.2006  for RGF=Newton Trajectory
C GS nach: Baron Peters et al. JCP 120 (2004) 7877-7886
C    und   WQ, JCP 2005 
c script version - modular structure
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
c                   Projector                          c
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      Integer NDIM,L,N3
      PARAMETER(NDIM=486,L=12,N3=162)
cc
      REAL*8 RotX(NDIM),RotY(NDIM),RotZ(NDIM),X(NDIM)
      REAL*8 PP(NDIM,NDIM),RR(NDIM,NDIM)
      REAL*8 Y(L,NDIM),RsN,ProMat(NDIM,NDIM),Rs(NDIM)
      REAL*8 xm,ym,zm,xyzmass(N3),gesmass
      Integer N,I,J,JJ,LA,Lha,Lnum,Natoms(N3) 
      Character*2  CharAtoms(N3)
c
      Data Natoms/6,6,7,6,6,8,6,7,6,6,8,6,7,6,6,8,
     * 6,7,6,6,8,6,7,6,6,8,6,7,6,6,8,6,
     * 7,6,6,8,6,7,6,6,8,6,7,6,6,8,6,7,
     * 6,6,8,6,7,6,6,8,6,7,6,6,8,6,7,6,
     * 6,8,6,7,6,6,8,6,7,6,6,8,6,7,6,8,
     * 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
     * 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
     * 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
     * 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
     * 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
     * 1,1/ 
c
       OPEN(9,FILE='a162chain.dat') 
       OPEN(10,FILE='a162Rs.dat')
       OPEN(11,FILE='a162param.dat',status='old') 
       OPEN(12,FILE='a162proj.dat') 
       OPEN(13,FILE='a162point.dat') 
       OPEN(33,FILE='a162poig.dat')
      OPEN(44,FILE='a162proto.txt',status='unknown',access='append')
c
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
       rewind 13
       READ(13,*)(X(I),I=1,N)
       rewind 9
       Do 31, J=1,L
       read(9,*) J1
       read(9,*)
       Do 31, JJ=1,N3
31     read(9,555) CharAtoms(JJ),(Y(J,3*(JJ-1)+I),I=1,3)
555    Format(1X,A3, 3F16.8)
C
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  overwrite the new calculated X(I) in the chain: 
       Do 331, I=1,N
331     Y(LA,I)= X(I)

       WRITE(44,*) 'guess point of 1.predictor of START a162cgGS1 '
       rewind 33 
       Do 7, J=1,N3
       write(44,*) Natoms(J), ( Y(LA,3*(J-1)+I ),I=1,3)
 7     WRITE(33,*) Natoms(J), ( Y(LA,3*(J-1)+I ),I=1,3)
       WRITE(33,*) '  '
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
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,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,*) ' '
       write(6,*) ' '
       write(44,*) ' '
       write(44,*) 'STOP ,  search direction error '
       write(44,*) ' '
       write(44,*) ' '
       stop 
      endif 
      Do 33 J=1,N    
33     Rs(J)=Rs(J)/RsN
       rewind 10
       write(10,50)(Rs(I),I=1,N)
       WRITE(44,*) ' Projector direction Rs '
      Do 47, J=1,N3
47     write(44,555) CharAtoms(J),(Rs(3*(J-1)+I),I=1,3)
c   reduce point to center of mass
        xm=0.0d0
        ym=0.0d0
        zm=0.0d0
        gesmass=N3
C mass weighting 
C       gesmass=0.0
      do 43,I=1,N3
c        IF(Natoms(I).eq.1) then
         xyzmass(I)=1.0d0
c         else
c         xyzmass(I)=Natoms(I)*2.0d0
c         ENDIF
c         gesmass=gesmass+ xyzmass(I)
43     continue
      do 44,I=1,N3        
           JJ=3*(I-1)      
         xm= xm + X(JJ+1)/ gesmass
         ym= ym + X(JJ+2)/ gesmass
         zm= zm + X(JJ+3)/ gesmass
44      continue
       write(6,*) 'old center of mass:', xm, ym, zm 
       do 45,I=1,N3         
         JJ=3*(I-1)      
         X(JJ+1)=-xm + X(JJ+1) 
         X(JJ+2)=-ym + X(JJ+2) 
         X(JJ+3)=-zm + X(JJ+3) 
45      continue
 
c rotation zero directions
      do 55,I=1,N3
           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  Projektion operator !!  P_rot is left,  P_Rs is right 
C                        umgekehrt funktioniert es !! 

      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
      call PROJECTOR(N,Rs,ProMat)
ccc
      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 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