PROGRAM A22cgGS2
C CCC for ALA  potential N=66  CCCCCCCCCCCCCCCCCCCCCCCC
C W.Quapp 30.11.2004/09.06.2006 for RGF=Newton Trajectory
C GS nach: Baron Peters et al. JCP 120 (2004) 7877-7886
c script version - modular structure  
c Projector (for additional use of z_mat coordinates)
c use of  grad in any point which is not stationary
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      Integer NDIM,L,N3
      PARAMETER(NDIM=66,L=32,N3=22)
      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,gesmass,xyzmass(N3)
      Character*2 CharAtoms(N3)
      Integer N,I,J,JJ,J1,LA,Lha,Lnum,Natoms(N3),I1,I2
      Data  Natoms/1,6,1,6,1,8,7,6,1,6,1,6,8,7,6,1,1,1,1,1,1,1/  
c
       OPEN(9,FILE='ala22chain.dat') 
       OPEN(10,FILE='ala22Rs.dat')
       OPEN(11,FILE='ala22param.dat',status='old') 
       OPEN(12,FILE='ala22proj.dat') 
       OPEN(13,FILE='ala22point.dat')
       OPEN(33,FILE='ala22poig.dat')
       OPEN(44,FILE='ala22proto.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)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  overwrite the new calculated X(I) in the chain: 
       Do 331, I=1,N
331     Y(LA,I)= X(I)
cc       WRITE(44,*) ' guess point of predictor '
        rewind 33 
        Do 7, J=1,N3
cc      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,*) '  '
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccc use skew special direction:
c       rewind 83
c       READ(83,*)(Y(L,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,J)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cc read gradient vector to be the search direction 
c        rewind 15
c        Do 17, J=1,N3
c17      READ(15,*)I1,I2,(Rs(3*(J-1)+I),I=1,3)
c and pull back the mass weighting of Gaussian03 gradient 
c      do 30,I=1,N3
c       JJ=3*(I-1)
c       xyzmass(I)=Dsqrt( Natoms(I)*2.0d0 ) 
c       IF(Natoms(I).ne.1) then
c         Rs(JJ+1) =  Rs(JJ+1)/ xyzmass(I)
c         Rs(JJ+2) =  Rs(JJ+2)/ xyzmass(I)
c         Rs(JJ+3) =  Rs(JJ+3)/ xyzmass(I)
c       ENDIF
c 30     continue
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c   Norm
      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(44,*) 'STOP ,  search direction error '
       write(6,*) ' '
       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=0.0d0
      do 43,I=1,N/3
c      IF(Natoms(I).eq.1) then
             xyzmass(I)=1.0d0
ccc  Natoms(I)
c         else
c             xyzmass(I)=1.0d0
ccc  Natoms(I)*2.0d0
c         ENDIF
ccc       gesmass=gesmass+ xyzmass(I)
43     continue
       gesmass=N3
      do 44,I=1,N3        
           JJ=3*(I-1)      
         xm= xm + X(JJ+1) *xyzmass(I)/ gesmass
         ym= ym + X(JJ+2) *xyzmass(I)/ gesmass
         zm= zm + X(JJ+3) *xyzmass(I)/ gesmass
44      continue
        write(6,*) 'old center of mass: ( ', xm, ym, zm, ')'
        write(44,*)'old c. o. m.: ( ', 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,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, umgedreht
c                                    funktioniert es leidlich 
      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)
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 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