PROGRAM A22ntGS2
C CCC for Ala  potential N=60  CCCCCCCCCCCCCCCCCCCCCCCC
C W.Quapp June 2007  for Newton Trajectory
C GS nach: Baron Peters et al. JCP 120 (2004) 7877-7886
C    und   WQ, JCP 122, iss17, (2005) 174106 (11 pages)
c script version - modular structure 
c Projector 
c
c Spezialfall einer schraegen Richtung fuer RS 
c  in 83,FILE='ala2Seite.dat'
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      Integer NDIM,L,N3
      PARAMETER(NDIM=58,L=30,N3=22)
      REAL*8 Y(L+1,NDIM),RsN,ProMat(NDIM,NDIM),Rs(NDIM)
      REAL*8 torad,btoa
      Integer N,I,J,JJ,LA,Lha,Lnum,Natoms(N3) 
      Character*2 CharAtoms(N3)
      Character*5 Zcoor(NDIM+2)
c CHASS:
      Data Zcoor/'ch2  ','cc3  ','cch3 ','oc4  ','occ4 ','dih4 ',
     1   'hc5  ','hcc5 ','dih5 ','hc6  ','hcc6 ','dih6 ','nc7  ',
     2   'ncc7 ','dih7 ','cn8  ','cnc8 ','dih8 ','cc9  ','ccn9 ',
     3   'dih9 ','oc10 ','occ10','dih10','hc11 ','hcn11','dih11',
     7   'hn12 ','hnc12','dih12','cc13 ','ccn13','dih13','hc14 ',
     8   'hcc14','dih14','hc15 ','hcc15','dih15','hc16 ','hcc16',
     7   'dih16','nc17 ','ncc17','dih17','cn18 ','cnc18','dih18',
     8   'hc19 ','hcn19','dih19','hn20 ','hnc20','dih20','hc21 ',
     7   'hcn21','dih21','hc22 ','hcn22','dih22'/
      Data CharAtoms/'h','c','c','o','h','h','n','c','c','o','h','h',
     1               'c','h','h','h','n','c','h','h','h','h'/
      Data Natoms/1,6,6,8,1,1,7,6,6,8,1,1,6,1,1,1,7,6,1,1,1,1/
      OPEN(9,FILE='ala2chain.dat') 
      OPEN(10,FILE='ala2Rs.dat')
      OPEN(11,FILE='ala2param.dat',status='old') 
      OPEN(12,FILE='ala2proj.dat') 
      OPEN(25,FILE='ala2told.dat')
      OPEN(83,FILE='ala2Seite.dat')
      OPEN(33,FILE='ala2poig.dat')
      OPEN(44,FILE='ala2protocol.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 ! 
C
       N=NDIM 
       read(11,*) JJ, LA
       JJ=L
cccccccccccccccccccccccccccccccccccccc       
C for GE only !!   not for usual NTs c 
ccc !!      IF(LA.gt.3) GOTO 3             
cccccccccccccccccccccccccccccccccccccc
       rewind 9
       Do 311, J=1,L+1
        read(9,50) (Y(J,I),I=1,N)
311     read(9,50) 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        rewind 33 
        Do 7, J=1,5
 7      WRITE(33,*) Zcoor(J), Y(LA,J)
        Do J=6,N
        WRITE(33,*) Zcoor(J+1), Y(LA,J)
        enddo
        WRITE(33,*) '  '
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  projector direction:
c       Lha=  3
        Lnum= 1
c       IF(LA .gt. 6) Lnum=LA-Lha
ccc 2006 Seitenrichtung:cccccccccc
c  c     if(LA.lt.7) then
        READ(83,*)(Rs(I),I=1,6)
        READ(83,*)(Rs(I),I=6,N)
        goto 300
ccc       else
c        Do 28 J=1,N
c28      Rs(J)= Y(L+1,J) - Y(Lnum,J)
c  c      endif
cccccccccccccccccccccccccccccccccccccccc
c  Normierung
c  (bohr,radian) <<--(angstroem,degree) of internal coordinates
        torad=dacos(-1.0d0)/180.0d0
        btoa = 0.52917706d0
c  stretch, bend, torsion alternately ( for GamessUS )
c  note: fix dih6 !!!  
         Rs(1)=Rs(1)/btoa
         Rs(2)=Rs(2)/btoa
         Rs(3)=Rs(3)*torad
         Rs(4)=Rs(4)/btoa
         Rs(5)=Rs(5)*torad
         If(N.gt.5) then
          Do 31 I=2,N3-4
          Rs(3*I)=Rs(3*I)/btoa
          Rs(3*I+1)=Rs(3*I+1)*torad
31        Rs(3*I+2)=Rs(3*I+2)*torad
         Endif
          Rs(57)=Rs(57)/btoa
          Rs(58)=Rs(58)*torad
300    continue
          If(LA.eq.2) write(44,*) ' first Rs is '
          If(LA.eq.2) write(44,50) (Rs(I),I=1,5)
          If(LA.eq.2) write(44,50) (Rs(I),I=6,N)
c   Normierung: in Chartesian Metric (old part ... )
      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
      goto 10
cc approximate GE by RS=Tangent of former step
cc for GE to first eigenvalue only !!
3     IF(LA.gt.3) then
       Write(44,*) ' RS=TANGENT for GE search '
       READ(25,*)(Rs(I),I=1,N)
      endif
10    rewind 10
      write(10,50)(Rs(I),I=1,N)
cccccccccccccccccccccccccccccccccccccccccccccc
      call PROJECTOR(N,Rs,ProMat)
cccccccccccccccccccccccccccccccccccccccccccccc
      rewind 12
      Do 37, J=1,N
      WRITE(12,50)(ProMat(J,I),I=1,N)
37    WRITE(12,50)
50    FORMAT(1X,3F20.10)
52    FORMAT(1X,2A,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
      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     
ccccccccccccccccccccccccccccccccccccccccccccccccccc
       SUBROUTINE PROJECTION(N,VEC,ProMat)
c  projection of vector VEC of dim N
c  by Projection Matrix ProMat of dim(N,N)
      REAL*8 VEC(N),ProMat(N,N),Proj(N)
      INTEGER N,J,JJ
      Do 1 J=1,N
      Proj(J)=0.0d0
      Do 2 JJ=1,N
2     Proj(J)=Proj(J) + Vec(JJ)*ProMat(J,JJ)
1     Continue
      Do 3 J=1,N
3     Vec(J) = Proj(J)
      RETURN
      END
cccccccccccccccccccccccccccccccccccccccccccccccccccc
      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