       !****************************************************************************************************
       !                          A Subroutine of Main, evaluating the system energy                       *
       !****************************************************************************************************
       
       Subroutine EnergyEvaluating ( PDBid_A, PDBid_B, SITE, Nstp, Nitr, rn, Engy )
       
       Implicit None
       
! ---- Variables : --------------------------------------------------------------------------------- !

       CHARACTER   PDBid_A*4
       CHARACTER   PDBid_B*4
              
       Integer     SITE
       
       Integer     Nstp, Nitr
       
       Real        rn(1:Nstp, 1:SITE, 1:3)    ! The coordinates of the strucutre in the transition pathway
       
       Real        Engy(1:Nstp)
       
! ---- Variables :  <o> <o> <o> <o> <o> <o> <o> <o> <o> <o> <o> <o> <o>   Read in [Energy_Function.inp]       
       
       Character*120  Chartemp
       
       Real        gama  ! The Force constant
       Real        rcut  ! The cut-off distance
       
       Real        beta  !  The parameters in the energy function, barrier reduction 
       Real        alpa  !  The parameters in the energy function, equilibrium difference  
       
! ---- Variables :   <E>.<E>.<E>.<E>.<E>.<E>.<E>.<E>.<E>.<E>.<E>.<E>.<E>.  Calculate Eng1 & Eng2
       
       Real        r0A(1:SITE,1:3)  ! The coordinates at equilibrium position A
       Real        r0B(1:SITE,1:3)  ! The coordinates at equilibrium position B
       
       Real        Eng1(1:Nstp)
       Real        Eng2(1:Nstp)
       
       Real        rij, rijA, rijB
       
! ---- Variables :  *| *| *| *| *| *| *| *| *| *| *| *| *| *| *| *| *| *|  Output the energy        
       
       Real        df(1:SITE, 1:3)         ! Deformation vectors 
       Real        vc(1:SITE, 1:3)         ! [V]
       
       Real        DotPduct3D
       
       Real        dotDV         ! dot product of df and vc
       
       Real        absdf         ! the magnitude of deformation vector df
       
       Real        rate(1:Nstp)  ! the transition rate, V.D/|D|
       
       Real        divs(1:Nstp)  ! diverse from div
       
       Real        rpv(1:SITE, 1:3)    ! the projection of V on D
       
       Real        rti(1:SITE, 1:3)    ! The temp coordinates of transition states
       
       Real        rmsd                ! Function to calculate the rmsd
       
! ---- Variables : --------------------------------------------------------------------------------- !
       
	   Character ACE*2 ! the string will be stored here 
       
       Integer     n, k, ni, nj, ct1, ct2
       
       
       
       
! ==== Body of Subroutine EnergyEvaluating ========================================================= !     
       
! <o> <o> <o> <o> <o> <o> <o> <o> <o> <o> <o> <o> <o> <o> <o> <o> <o>  Read in [Energy_Function.inp]       
       
        Open(130, File='Energy_Function.inp')

        Read(130,*) CharTemp
        Read(130,*) CharTemp
        Read(130,*) gama
        
        Read(130,*) CharTemp
        Read(130,*) CharTemp
        Read(130,*) rcut        

        Read(130,*) CharTemp
        Read(130,*) CharTemp
        Read(130,*) beta        
        
        Read(130,*) CharTemp
        Read(130,*) CharTemp
        Read(130,*) alpa        
                
        Close (130) ! 'Energy_Function.inp' 
        
! <E>.<E>.<E>.<E>.<E>.<E>.<E>.<E>.<E>.<E>.<E>.<E>.<E>.<E>.<E>.<E>.<E>. Calculate Eng1 & Eng2
        
        Do n=1, SITE
           
           r0A(n,1) = rn(1,n,1)
           r0A(n,2) = rn(1,n,2)
           r0A(n,3) = rn(1,n,3)

           r0B(n,1) = rn(Nstp,n,1)
           r0B(n,2) = rn(Nstp,n,2)
           r0B(n,3) = rn(Nstp,n,3)
        
        Enddo ! n 
        
        ! ------------------------- - - - -
        
        ct1=0
        ct2=0
        
        Do k=1, Nstp
        
            Eng1(k)=0.0
            Eng2(k)=0.0
            
            Do ni=1, SITE-1
            Do nj=ni+1, SITE
                
                rij  = Sqrt( (rn(k,ni,1)-rn(k,nj,1))**2 + (rn(k,ni,2)-rn(k,nj,2))**2 + (rn(k,ni,3)-rn(k,nj,3))**2  )
                rijA = Sqrt( ( r0A(ni,1)-r0A(nj,1) )**2 + ( r0A(ni,2)-r0A(nj,2) )**2 + ( r0A(ni,3)-r0A(nj,3) )**2  )
                rijB = Sqrt( ( r0B(ni,1)-r0B(nj,1) )**2 + ( r0B(ni,2)-r0B(nj,2) )**2 + ( r0B(ni,3)-r0B(nj,3) )**2  )
                
                If( rijA<= rcut ) Then
                    
                    Eng1(k) = Eng1(k)   +   0.5 * gama * ( rij - rijA )**2
                    ct1=ct1+1
                EndIf

                If( rijB<= rcut ) Then
                    
                    Eng2(k) = Eng2(k)   +   0.5 * gama * ( rij - rijB )**2
                    ct2=ct2+1
                EndIf
                                
            Enddo ! nj
            Enddo ! ni
        
        Enddo ! k
        
!        Write(*,*) ct1, ct1/(1.0*Nstp)
!        Write(*,*) ct2, ct2/(1.0*Nstp)
        
        ! ======================== = = = = = =
        
        Do k=1, Nstp
            
            Engy(k) =  0.5 * (     Eng1(k) + (Eng2(k)-alpa)   -  Sqrt(   ( Eng1(k) - (Eng2(k)-alpa) )**2  + 4.0 * beta * beta )      ) 
        
        Enddo ! k

! *| *| *| *| *| *| *| *| *| *| *| *| *| *| *| *| *| *| *| *| *| *| *| Output the energy        
        
        Open(131, File='['//PDBid_A//'-'//PDBid_B//']_Energy Evaluation['//ACE(Nitr,2)//'].dat' )
        
        ! # - # - # - # - # - # - # - # - # - # - # - # -  -   -    - rate(k)
        
        Do n=1, SITE
            df(n,1)  =  r0B(n,1) - r0A(n,1)
            df(n,2)  =  r0B(n,2) - r0A(n,2)
            df(n,3)  =  r0B(n,3) - r0A(n,3)
        Enddo ! n
        
        absdf =  DotPduct3D( df, df, SITE ) 
        
        Do k=1, Nstp
        
            Do n=1, SITE
                vc(n,1) =  rn(k,n,1) - r0A(n,1)
                vc(n,2) =  rn(k,n,2) - r0A(n,2)
                vc(n,3) =  rn(k,n,3) - r0A(n,3)
            Enddo
            
            dotDV = DotPduct3D( vc, df, SITE )
            
            rate(k) = dotDV/absdf
        
        Enddo ! k
        
        ! + * + * + * + * + * + * + * + * + * + * + * + *  *   *    * divs(k)
        
        Do k=1, Nstp
            
            Do n=1, SITE
                rpv(n,1) = r0A(n,1) + rate(k) * df(n,1)
                rpv(n,2) = r0A(n,2) + rate(k) * df(n,2)
                rpv(n,3) = r0A(n,3) + rate(k) * df(n,3)
            Enddo
            
            Do n=1, SITE
                rti(n,1) = rn(k,n,1)
                rti(n,2) = rn(k,n,2)
                rti(n,3) = rn(k,n,3)
            Enddo
            
            divs(k) = rmsd(SITE, rpv, rti)
            
        Enddo ! 
        
        ! # * # * # * # * # * # * # * # * # * # * # * # *  *   *    * Output
        
        Write(131, '(A)') '        Step    Rate      diversity          Energy         Eng1           Eng2'
        Do k=1, Nstp
!             Write(131,*) k, rate(k), divs(k), Engy(k), Eng1(k), Eng2(k)
             Write(131,'(I5, 5(F12.6,4x) )') k, rate(k), divs(k), Engy(k), Eng1(k), Eng2(k)
        Enddo ! k
        
        Close(131) ! '['//PDBid_A//'-'//PDBid_B//']_Energy Evaluation'//'.dat'
        
! =====E==N===D====================================================================================== !
       Return
       End            
       
       
       
       
       
       
       
       
       
       
       
        